Identifying island safe havens to prevent the extinction of the World’s largest lizard from global warming

Abstract The Komodo dragon (Varanus komodoensis) is an endangered, island‐endemic species with a naturally restricted distribution. Despite this, no previous studies have attempted to predict the effects of climate change on this iconic species. We used extensive Komodo dragon monitoring data, climate, and sea‐level change projections to build spatially explicit demographic models for the Komodo dragon. These models project the species’ future range and abundance under multiple climate change scenarios. We ran over one million model simulations with varying model parameters, enabling us to incorporate uncertainty introduced from three main sources: (a) structure of global climate models, (b) choice of greenhouse gas emission trajectories, and (c) estimates of Komodo dragon demographic parameters. Our models predict a reduction in range‐wide Komodo dragon habitat of 8%–87% by 2050, leading to a decrease in habitat patch occupancy of 25%–97% and declines of 27%–99% in abundance across the species' range. We show that the risk of extirpation on the two largest protected islands in Komodo National Park (Rinca and Komodo) was lower than other island populations, providing important safe havens for Komodo dragons under global warming. Given the severity and rate of the predicted changes to Komodo dragon habitat patch occupancy (a proxy for area of occupancy) and abundance, urgent conservation actions are required to avoid risk of extinction. These should, as a priority, be focused on managing habitat on the islands of Komodo and Rinca, reflecting these islands’ status as important refuges for the species in a warming world. Variability in our model projections highlights the importance of accounting for uncertainties in demographic and environmental parameters, structural assumptions of global climate models, and greenhouse gas emission scenarios when simulating species metapopulation dynamics under climate change.


| INTRODUC TI ON
Tropical island endemics are highly vulnerable to global change (Sodhi, Koh, Brook, & Ng, 2004). Island species often possess attributes, such as low dispersal rates, small population sizes, low heterozygosity, and strong local adaptation, which constrain the ecological and evolutionary responses needed to counter rapid environmental change (Fordham & Brook, 2010;Frankham, 1998). The impact of climate change on tropical island species may be exacerbated because tropical species are particularly susceptible to increases in temperature associated with global warming (Tewksbury, Huey, & Deutsch, 2008) and sea-level rise is projected to flood low-lying island habitats, leading to a marked deterioration in quality and availability of habitats (Menon, Soberón, Li, & Peterson, 2010;Wetzel, Kissling, Beissmann, & Penn, 2012). Consequently, the impacts of anthropogenic climate change are expected to be disproportionately large for tropical archipelagos, such as Indonesia, which harbor extremely high levels of endemism (Sodhi et al., 2004).
The Komodo dragon Varanus komodoensis is the world's largest lizard and a unique and iconic island-endemic species. They occupy a highly restricted range distribution (Ciofi & De Boer, 2004) and are considered "Vulnerable" by the International Union for Conservation of Nature (IUCN, 2017). However, the last IUCN Red List assessment for the species was done over 20 years ago (World Conservation Monitoring Centre, 1996), and since this time, there has been considerable research on trends in Komodo dragon population abundance and range area (Imansyah, Jessop, Ciofi, & Akbar, 2008;Jessop et al., 2007Jessop et al., , 2018Purwandana et al., 2014Purwandana et al., , 2015. Therefore, a reassessment of the present-day threat status of Komodo dragons is overdue. Populations of Komodo dragons currently persist on five islands in southeastern Indonesia (Figure 1). Fewer than 3,000 individuals live in an area of ~600 km 2 , split across the islands of Komodo, Rinca, Nusa Kode, and Gili Motang within the World Heritage-listed Komodo National Park (KNP; Purwandana et al., 2014). On the largest island of Flores (13,540 km 2 ), outside of KNP's boundaries, range contraction and population declines have been documented, mainly due to expansion of human settlements, illegal hunting of prey species, and forest clearance for agriculture (Ariefiandy, Purwandana, Natali, et al., 2015;Ciofi & De Boer, 2004). Komodo dragons on Flores are currently only found in a small number of habitat fragments along the north and west coasts of the island. Only ~80 km 2 of potential Komodo dragon habitat is protected for conservation reflecting these islands' status as important refuges for the species in a warming world. Variability in our model projections highlights the importance of accounting for uncertainties in demographic and environmental parameters, structural assumptions of global climate models, and greenhouse gas emission scenarios when simulating species metapopulation dynamics under climate change.

K E Y W O R D S
climate change, conservation management, demographic model uncertainty, extinction risk, population viability, sea-level rise F I G U R E 1 The current-day geographic range of Komodo dragons is spread across five islands in East Nusa Tenggara, Indonesia. Islands within the Komodo National Park are indicated by orange circles. Orange shading shows the location of reserves on Flores purposes on Flores (Ariefiandy, Purwandana, Nasu, Surahman, & Jessop, 2015;Ciofi & De Boer, 2004).
Komodo dragons have persisted for millennia in their current range despite interactions with early hominids and substantial climatic-and sea-level-related habitat changes during the last glacial maxima (36-17 ka;van den Bergh et al., 2009) and subsequent Holocene warming. However, climate change models project that over the next century, Indonesia will experience unprecedented rates of both temperature rise and reduced rainfall (Boer & Faqih, 2004;Diffenbaugh & Giorgi, 2012), leading to a prolonged dry season with increased fire frequency and decreased soil moisture (Fernandes et al., 2017). This is projected to cause a contraction of mesic forest cover and an expansion of drier vegetation communities, such as savannah woodland (Bickford, Howard, Ng, & Sheridan, 2010). This vegetation transformation is likely to negatively impact Komodo dragons by altering habitat and prey availability with impacts on both survival and reproduction (Jessop et al., 2006(Jessop et al., , 2007. In addition, rising sea levels are likely to inundate the low-lying valleys (Heaney, 1991) that currently support the highest densities of Komodo dragons, leading to a permanent loss of their preferred lowland habitat (Purwandana et al., 2014). When combined with the existing issues of human-induced habitat loss, this could be disastrous for the species.
Given the Komodo dragon's potential susceptibility to climate change and its economic importance to Indonesia through tourism (Walpole & Goodwin, 2000;Walpole & Leader-Williams, 2002), onground resource managers urgently need near-term projections of Komodo dragon metapopulation dynamics to inform climate change conservation policies. Our specific aims were to: (a) estimate the effects of projected regional temperature increase and sea-level rise on the habitat suitability, metapopulation structure, and total population size of Komodo dragons between 2010 and 2050, (b) determine whether the risk of extirpation differed across island populations, and (c) investigate the effect of structural uncertainties in global climate models on projections of extinction risk for Komodo dragons.
The novelty of our approach to range dynamics modeling under climate change is the inclusion of multiple sources of uncertainty in a process-based model, including important structural uncertainties in coupled atmosphere-ocean general circulation models (AOGCMs), which are often ignored in ecological studies. Accounting for AOGCM structural uncertainties in climate change projections is important, because assumptions regarding climate sensitivity (the equilibrium global mean warming for a doubling of CO 2 concentration) and aerosol forcing can have an equally strong influence on projections of regional climate change as the choice of greenhouse gas emission scenario (Fordham, Wigley, & Brook, 2011;IPCC, 2013). We show that although choice of AOGCM structure causes large uncertainty in ecological model projections of future range and abundance of F I G U R E 2 Schematic diagram illustrating our method for coupling ecological niche (ENM) and stochastic population models to create a nichepopulation model (NPM), which was used to predict the future range dynamics of Komodo dragons. The presence and environmental data (step 1) were used in the ENM (step 2) to project range-wide, annual, habitat suitability under current and future climate conditions (step 3). Habitat suitability projections from the ENM (2010-2050) were used to define the spatial habitat patch structure (step 4) for stochastic metapopulation models (step 5) that simulated spatiotemporal abundance and occupancy patterns across all five islands of the species' global distribution (the NPM, step 6). Simulations accounted for environmental and demographic parameter uncertainty in the NPM (steps 7 and 8), with sensitivity analysis used to identify parameters exerting the greatest influence on the model projections (step 9) Komodo dragons, global warming is forecast to cause the extirpation of entire island populations in the near future.

| ME THODS
To predict the potential effects of future climate change on Komodo dragons, we used a recently developed approach that couples an ecological niche model (ENM, also known as a species distribution model) and a metapopulation simulation model (spatially explicit population viability analysis), to create a coupled niche-population model (NPM). The NPM was used to simulate landscape-scale population processes, including dispersal and source-sink dynamics, across the species' geographic range, allowing projections of metapopulation abundance, island population extirpation rates, and habitat patch occupancy under various climate change scenarios (Figure 2; Fordham, Akçakaya, Araújo, Keith, & Brook, 2013;Fordham, Mellin, et al., 2013).
The NPM approach has several advantages over alternative methods, such as pattern-based ENMs or nonspatial population viability analyses. First, it permits integration of dispersal and metapopulation dynamics into forecasts of a species' geographic range, providing more ecologically realistic predictions of a species' response to climate change (Anderson et al., 2009). Second, it yields direct estimates of extinction risk, in addition to vulnerability measures based on projected changes in range area, habitat patch occupancy, and total habitat suitability .
Third, this approach inherently incorporates demographic responses to multiple, and often synergistic, processes of global change (Fordham, Resit Akçakaya, et al., 2012). Together, these advantages mean that, where data permit, coupled NPMs provide better predictions of range and population dynamics in response to climate change (Fordham et al., 2018).
Komodo dragons have been the focus of intensive long-term and large-scale population monitoring efforts, generating temporal and spatial occupancy and demographic data (Ariefiandy, Purwandana, Natali, et al., 2015;Ciofi et al., 2007;Purwandana et al., 2014Purwandana et al., , 2015Sastrawan & Ciofi, 2002). We used these data to parameterize our NPM and make predictions of climate change impacts on the future abundance and distribution of Komodo dragons. We allowed uncertainty in key demographic parameters (such as survival rates, dispersal, and density-dependant population growth rates), and in ecological niche model projections of habitat suitability, to propagate through to model projections of future Komodo dragon range and abundance (Fordham, Haythorne, & Brook, 2016). We accounted for important structural assumptions in AOGCM climate change projections, as well as greenhouse gas emission scenario, on predictions of extinction risk for Komodo dragons. To do this, we considered six different plausible future climate scenarios from an ensemble of AOGCMs (Fordham et al., 2011). Each scenario was based on different combinations of global greenhouse gas emissions and climate model structural parameters for climate sensitivity to CO 2 concentration and aerosol forcing. AOGCM sensitivity is classified by the amount of temperature change resulting from a doubling of atmospheric CO 2 concentration since industrialization (Wigley et al., 2009)

| Ecological niche model
We modeled the current distribution of Komodo dragons using geographically referenced occurrence data (N = 4,028; see Appendix S1, Table S1.1), intersected with climate and landscape variables at a 1 × 1-km grid-cell resolution (see step 1 in Figure 2). We built a consensus ensemble ecological niche model (ENM; step 2 in Figure 2) using the R package biomod2 (Marmion, Parviainen, Luoto, Heikkinen, & Thuiller, 2009;Thuiller, Lafourcade, Engler, & Araújo, 2009). We used four different modeling algorithms (generalized linear models, generalized additive models, generalized boosted-regression models, and MAXENT) with interactions allowed, but depth limited to 1 in the GLM and GAM models, and 3 in the GBM models. We used an ensemble ENM approach because it avoids potentially important biases inherent with using single ENM models (Araújo et al., 2019;Araújo & New, 2007). We weighted the contribution of each of the individual ENMs in the ensemble using the True Skill Statistic (TSS) model evaluation scores. In doing so, individual models that had better evaluation scores contributed more to the ensemble model output (Marmion et al., 2009).
Our correlative ENM approach is based on the fundamental assumption that the distribution of a species, in approximate equilibrium with its environment, is a good indicator of its ecological requirements (Elith & Leathwick, 2009). Use of proximal predictors of ecological suitability (i.e., those that directly influence geographic distributions or demographic mechanisms) is the gold standard for reducing uncertainty in projections from species distribution models (Araújo et al., 2019). However, distal predictors, including landscape and some climate variables, can be used as proxies for proximal constraints such as habitat and prey availability (Elith & Leathwick, 2009). Because robust spatiotemporal data were not available to generate some proximal predictors of Komodo dragon geographic range (e.g., vegetation community and prey abundance), we used available distal surrogates (Appendix S1). In this context, temperature can be considered both a proximal and distal predictor variable, in that it has both a direct (proximal) influence on Komodo dragon physiology (and therefore fitness) and behavior, and an indirect (distal) influence on vegetation (habitat) and prey distributions.
We evaluated the ENM using both independent evaluation (by means of a separate dataset for Komodo dragon occupancy collected with camera traps) and cross-validation. The independent data provide a more stringent test of the predictive ability of the ENM for the current day (Elith & Leathwick, 2009). The calibration and evaluation of the ENM are described in detail in Appendix S1. The ensemble ENM was used to project future Komodo dragon habitat suitability at annual time steps (up to 2100) under AOGCM-projected changes in temperature and sea level for the six "plausible future" scenarios (step 3 in Figure 2). The resulting ENM projections of habitat suitability under the six "plausible futures," when taken together, accounted for uncertainties in greenhouse gas emission policy and AOGCM structural assumptions (see 2.1.2).
The six plausible future time series of Komodo dragon habitat suitability generated by the ensemble ENM were then used as a proxy for upper abundance (carrying capacity, K) in the coupled NPM (VanDerWal, Shoo, Johnson, & Williams, 2009; see step 4 in Figure 2). To relate ENM model projections of habitat suitability to K, we used an iterative approach to optimize the relationship between these two variables based on target values of current abundance (after Anderson et al., 2009) for each island taken from previously published studies (see Appendix S1). We only projected up to 2050 so that we avoided extrapolating to novel climates (Elith, Kearney, & Phillips, 2010;Elith & Leathwick, 2009) where mean temperatures after 2050 exceeded those of the ENM training dataset. This approach, while widely used as a method for avoiding extrapolation, assumes that the covariance structure of predictors remains constant through time. However, covariance structures can, in some instances, change through time (Jackson, Betancourt, Booth, & Gray, 2009).

| Current climate and environmental spatial data
Through a literature review and consulting with species experts (including authors on this paper), we identified temperature, topography, and distance from the coast as the landscape and climate variables likely to have the greatest impact on Komodo dragon distribution either directly or indirectly (by, e.g., influencing the distribution of their prey species).
Temperature directly affects the physiology of Komodo dragons, limiting occurrence in cooler, higher (above ~500 m) areas, as well as the hottest lowland areas (Auffenberg, 1981), and is also likely to indirectly affect distribution via its impact on vegetation communities and prey species distributions. Distance from the coast and slope are good proxies for land/vegetation cover (and therefore prey availability) because higher and steeper areas drive orographic rainfall patterns, which in turn affect vegetation type. Topographic slope can also act as a physical barrier to Komodo dragon movement.
Higher altitudes (>600 m elevation) and moist forest habitats constrain the distribution of the Komodo dragon, causing them to rarely be found more than 5-6 km from the shoreline (Auffenberg, 1981).
Conversely, flatter and drier areas are commonly characterized by lowland savannah woodland and lowland deciduous forests (<300 m elevation), the latter being the most preferred habitat for Komodo dragons (Auffenberg, 1981) and their prey .
We modeled temperature as a function of elevation, based on the moist adiabatic lapse rate. We first calculated a long-term sea-level average temperature baseline of 26.38°C (n = 56 years of observational data). We then applied an elevation lapse rate of 6.1°C decline per 1,000 m increase in altitude (Harris et al., 2014;Raxworthy et al., 2008) using a 250-m horizontal resolution digital elevation model (Shuttle Radar Topography Mission [SRTM]; Jarvis, Reuter, Nelson, & Guevara, 2008). We then upscaled the temperature raster to 1-km grid-cell resolution by multicell averaging. Further methodological details are available in Appendix S1. We used mean rather than minimum or maximum temperature in our models because there is low diurnal and seasonal variance in temperature in our study region (due to its equatorial location) and because longterm mean temperature was the metric most reliably recorded by the few weather stations located in the region.
Topographic slope was calculated as the maximum rate of change in elevation from each cell to all eight neighboring cells using the SRTM elevation raster and the slope tool in ArcGIS (v 10.3.1). We upscaled slope from the original SRTM cell size (250 m) to the 1-km grid-cell resolution by summing the number of smaller constituent cells (n = 16) that had a slope value greater than the 75th percentile for slope across the entire study area. Higher values in the final slope raster, therefore, indicate that a greater proportion of the cell area consisted of steep slopes.
We calculated distance from the nearest coastline to each grid cell in the study region using the 1 × 1-km resolution elevation raster. Cells with a mean elevation of <0 m were classified as the sea, while cells with mean elevation >0 m were assigned to land, enabling us to define the island's coastlines at the boundary between land and sea cells. We then calculated the distance from the centroid of each cell to the nearest point on the coastline using the R package "sp" (Bivand, Pebesma, & Gomez-Rubio, 2013). The distance to coast value for each grid cell was recalculated at each time step to account for sea-level rise (see details in section 2.1.2).
Our choice of spatial predictors was influenced, to some extent, by the availability of additional robust predictors of climate and environmental conditions. For example, commonly used estimates of rainfall and temperature, those based on the interpolation of meteorological data (i.e., downloadable from WorldClim2), were investigated and found to be highly uncertain for our study region (Fick & Hijmans, 2017; specifically their Figure 4). This is likely due to the complex topography and the sparseness of weather stations in the region. We explored the relationship between WorldClim 2 projections for rainfall and temperature and found them to be strongly correlated. The correlations for these datasets in our study area were between −0.73 and −0.94 (depending on the specific area of interest), indicating that as temperature increases, precipitation decreases. This relationship is likely to be due to the approach used to spatially interpolate the climate data, whereby a splining algorithm that uses elevation (which is highly correlated with temperature) as a covariate is used to spatially project precipitation (Fick & Hijmans, 2017). These exploratory analyses are shown in Appendix S1 (Table S1.2) and provide strong justification for not including precipitation as a predictor variable in an SDM that already contains temperature (Dormann et al., 2013). Furthermore, some of the ENM algorithms we applied in biomod2 perform poorly with colinear predictors (Dormann et al., 2012).
We explored the utility of several land use and land cover layers, but were unable to find one that reliably separated suitable and unsuitable land cover classes for Komodo dragons (these exploratory investigations are discussed in detail in Appendix S1, Figures S1.1 and S1.2). In addition, there were no datasets available that included future predicted land use/land-use change for our study region. These are required for generating projections of habitat suitability from the ENM, unless the current spatial pattern in land use is assumed to remain static over time, which is unrealistic (Hof et al., 2018). The implications of not including land-use data in this study are addressed in the Discussion.

| Projected changes in temperature and sea level
We considered six different plausible future climate scenarios using  Table 1), which account for uncertainty in the structural assumptions associated with AOGCMs.
We used an ensemble of seven atmosphere-ocean general circulation models (AOGCMs) to generate robust projections of regional changes in temperature, for each of the six plausible futures, on a 2.5-by-2.5-degree grid. The ensemble included a mixture of globally (GFDLCM20, MIROCMED, MRI-232A, UKHADCM3) and regionally (BCCRBCM2, CCSM-30, CSIRO-30) skillful AOGCMs (Harris et al., 2014). Skill was assessed based on the ability of the AOGCM to replicate observed climatic conditions (Fordham et al., 2011). Future temperature and sea-level change time series were generated for three (2.5° × 2.5°) grid cells delineating the study region.
For each of the six plausible future climate scenarios (Table 1) (1-km grid cells) using the "change factor" empirical method, whereby the low-resolution change from a AOGCM is added directly to a high-resolution baseline of observed climatology (Hulme, Raper, & Wigley, 1995). Similarly, we applied the annual sea-level forecasts to the land elevation raster by subtracting the projected amount of sea-level rise from the elevation value at each 1-km grid cell. Thus, the future spatial layers of distance from coast explicitly take account of forecast sea-level rise by removing

| Capture-mark-recapture analyses
We developed capture-mark-recapture (CMR) models in "Rmark" (Laake, 2013) to estimate annual survival rates for each age class of Komodo dragons (hatchlings, juveniles, and adults), across the different island populations (step 5 in Figure 2). These analyses used data collected from more than 1,000 marked individuals over six capture occasions between 2003 and 2013 (Purwandana et al., 2014.  Laver et al., 2012). Age class thresholds were adjusted to account for this interisland variation in age and size of maturity (Appendix S1, Table S1.3). The CMR survival analysis is described in detail in Appendix S1.

| Coupled niche-population model (NPM)
We simulated the metapopulation dynamics of Komodo dragons across the species' entire current-day range using stage-based stochastic matrix models in RAMAS GIS v5 (Akcakaya & Root, 2005) and a coupled niche-population modeling (NPM) framework . We used these models to project changes in range-wide and island-specific abundance, patch occupancy, and risk of extinction under the six plausible future climate scenarios (step 6 in Figure 2). Our NPM framework assumes that estimates of apparent survival rates (from CMR models of different island populations) capture responses to varying environmental conditions, such as prey availability, which is known to vary between islands Ariefiandy, Purwandana, Coulson, Forsyth, & Jessop, 2013;Jessop et al., 2007;Laver et al., 2012;Purwandana et al., 2014Purwandana et al., , 2015Purwandana et al., , 2016.
The spatial structure of the modeled metapopulation (the size and location of subpopulations within and across island populations) was based on annual predictions of habitat suitability from the ENM (step 4 in Figure 2) transformed to carrying capacity (K).
The transformation used a formula that was optimized iteratively to ensure that the baseline abundance on each island approximated recent island-based upper abundance estimates (see Appendix S1, Table S1.4; Jessop et al., 2007;Purwandana et al., 2014;Purwandana et al., 2015). This optimization approach is widely used when configuring NPMs (Anderson et al., 2009) and addresses the problem that the relationship between ENM suitability and upper abundance is not necessarily linear (Thuiller et al., 2014 Jessop et al., 2018;Sastrawan & Ciofi, 2002; Appendix S1, Figure S1.3). A detailed description of the approach used to delineate the spatial structure of the metapopulation is in Appendix S1. The entire Komodo dragon metapopulation for the present day, thus consisted of 43 panmictic populations, spread across five islands.
We assumed that dispersal between islands does not occur based on recent genetic, CMR, and radiotelemetry evidence (Jessop et al., 2018  . Hatchlings were excluded from the density-dependent calculation because they are arboreal and are therefore not competing for the same resources as the older life stages. A detailed description of the population model parameterization is given in Appendix S1. Demographic processes were simulated using a stage-structured, postbreeding, female-only stochastic population model (step 6 in Figure 2). We modeled only females because there is no strong evidence that vital rates differ between sexes and the species is polygamous (Purwandana et al., 2014. The stochastic models were parameterized using stage class-specific survival rates (see CMR analysis), fecundity rates, and dispersal rates. The maximum rate of intrinsic population growth (Rmax) for all populations was estimated to be 1.18 based on our own Pradel models of population growth and previous studies by Purwandana et al. (2014). We accounted for interisland differences in survival rates based on our CMR estimates (see Appendix S1 for details).

| Simulations
We ran six baseline spatiotemporally explicit NPMs, one for each of the six plausible future climate scenarios (described in Table 1 Figure 2 and Appendix S1, Table S1.5). We varied these parameter values to account for documented and suspected uncertainties in demographic (e.g., vital rates, density dependence, Rmax) and environmental parameters (e.g., uncertainty in the association between projected habitat suitability and upper abundance). Each NPM was run at annual time steps from 2010 to 2050, with a 20-year burn-in period, using 2010 environmental conditions to achieve an equilibrium metapopulation structure (step 8 in Figure 2). Each run of the NPM included 1,000 simulations that accounted for environmental and demographic stochasticity (which resulted in a total of 1,200,000 NPM iterations: 200 Latin hypercube iterations of the six plausible future NPMs, each run 1,000 times with stochasticity).

| Niche-Population Model (NPM) output
Under each climate change scenario, we report the median (50th percentile of all model simulations) number of female Komodo dragons projected to be alive each year (2010 to 2050), and the annual relative change in habitat patch occupancy to 2050 (bounded by the 5th and 95th percentiles from all simulations) for the entire metapopulation and for each island separately. We also report the median expected minimum abundance (EMA), which is a continuous integrated metric reflecting risks of both a decline in abundance and extinction risk (McCarthy & Thompson, 2001). In addition to these outputs, we also report the results of a sensitivity analysis using SARDM, which was done for two greenhouse gas emission scenarios (RCP 8.5 and RCP 2.6), both with midrange climate sensitivity and aerosol forcing (step 9 in Figure 2).

| Ecological niche model (ENM)
Evaluation of the ensemble ENM showed that the model had good

| Capture-mark-recapture analyses
Annual survival rates differed across the four island populations

| Coupled niche-population model (NPM)
Our NPM predicted a decline in (female) Komodo dragon metapopulation size by 2050 under all six plausible future climate change scenarios (Figure 4). Global greenhouse gas emission scenario (i.e.,   Results from a sensitivity analysis that tested the influence of not applying the habitat quality multiplier to small island populations showed that even if we had treated all reserves as equal, the fate of the smaller island populations on Kode and Motang would have remained bleak. There was no difference in predicted extirpations under the reference scenarios and a difference in EMA of just 2-3 individuals under the policy scenarios when the multiplier was applied to these populations (Appendix S4, Table S4.2).

| D ISCUSS I ON
The aims of our study were to estimate the effects of projected  the probability of occurrence decreasing rapidly beyond 2 km from the coastline, owing to a nonlinear interaction between temperature and distance from the coast. This relationship between temperature and distance from coast is likely to strongly influence the distribution and composition of vegetation communities in the study region (Auffenberg, 1981), influencing habitat suitability for Komodo drag- these two habitats, which generally contain lowland deciduous forests that grow in a narrow band, extending a maximum of ~5 km inland from the coast (Auffenberg, 1980(Auffenberg, , 1981Forth, 2010 Purwandana et al., 2016).
A reduction in the quality of vegetation is being observed across much of the current distribution of Komodo dragons as a possible result of recent climate change, causing lower prey availability, contributing to reduced population densities (Ariefiandy, Purwandana, Natali, et al., 2015). Projected increases in temperature under the  (Auffenberg, 1981;Heaney, 1991).  Heller and Zavaleta (2009). These include expanding upland reserve networks and translocating dispersal-limited species.
Our work identified populations on the two largest protected islands in KNP (Komodo and Rinca) to fare best under climate change.
Therefore, it is imperative that habitats on these islands be managed to reflect their status as important refuges for Komodo dragons (Kearney, Adams, Fuller, Possingham, & Watson, 2018). We suggest that park managers should act now to limit any processes that could exacerbate future habitat loss or degradation in these areas, such as the expansion of ecotourism (Boakes, Fuller, & McGowan, 2018;Hole et al., 2009 Several other key processes were not considered in our study, which, if accounted for, would most likely synergize with climate change to worsen the future extinction risk of Komodo dragons (Fordham & Brook, 2010). The approach we have used here assumes that biotic interactions remain constant through time in a shifting climate, which is unlikely to be true (Fordham, Akcakaya, et al., 2013). Accounting for biotic interactions directly in the model, and the potential for these interactions to decouple, would likely amplify extinction risk (Lurgi, López, & Montoya, 2012).
For example, projected decreases in annual rainfall for eastern Indonesia in the 21st century are likely to increase aridity in the region, leading to a reduction in primary productivity and major structural changes to vegetation communities (Naylor, Battisti, Vimont, Falcon, & Burke, 2007). Such changes would have complex effects on Komodo dragon distribution and abundance, ranging from direct (environmental thermal heterogeneity) to indirect (vegetation composition change and altered fire regimes) effects.
Unfortunately, the precipitation variables commonly used in ENMs (e.g., those downloadable from WorldClim2; Fick & Hijmans, 2017) are highly uncertain for our study region and therefore could not be used (see Methods and Appendix S1).
Komodo dragons are sensitive to land-use change that has occurred in nonprotected habitats on Flores over recent decades (Ciofi & De Boer, 2004). Therefore, ongoing habitat loss and conversion could exacerbate forecast reductions in Komodo dragon abundance and distribution on Flores, which is the only island in the distribution of Komodo dragons where habitat loss continues. Moreover, spatial data were not available for the distribution and abundance of the multiple prey species Komodo dragons feed on across their range, preventing us from directly modeling species interactions and assuming that prey availability is captured implicitly in ENM estimates of habitat suitability (Dormann et al., 2012) and CMR estimates of survival. While we were not able to model the impact of changes in prey availability or current and future land use due to a lack of suitably resolved spatial datasets, this should be prioritized in future research if the spatial layers and data become available.

| CON CLUS IONS AND IMPLIC ATIONS
Indonesia's biodiversity, with its high levels of endemism and numerous insular populations, is expected to incur severe impacts from global change (Fordham & Brook, 2010;Sodhi et al., 2004). Our projections for an iconic, and highly range-restricted large-bodied, predatory lizard support this general assertion and predict that Komodo dragons will be subject to a dramatic population decline and possible extinction due to global warming by 2050 (unless the most optimistic future climate scenarios are realized). The prospect of countering the effects of climate change on Komodo dragons, alongside other agents of rapid global change, is clearly a daunting task facing Indonesia. We strongly advocate that national and provincial conservation agencies act now to address impending climate change impacts on the islands within the Komodo dragon's range, bearing in mind that protected areas on the islands of Komodo and Rinca are likely to provide the best "future-proofed" habitats for Komodo dragons. Without quick action to mitigate climate change impacts in KNP and on the island of Flores, we risk committing Komodo dragons-a globally iconic species-to extinction.

ACK N OWLED G M ENTS
We are grateful to the Directorate General of Forest Protection NCAR is supported in part by the US National Science Foundation.
The Coobowie writing retreat, funded by the University of Adelaide's Environment Institute, supported the revision of this manuscript.

CO N FLI C T O F I NTE R E S T
None declared.   Tsoar, A., & Kadmon, R. (2006). Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill