Climate change threatens biodiversity directly by influencing biophysical variables that drive species' geographic distributions and indirectly through socio-economic changes that influence land use patterns, driven by global consumption, production and climate. To date, no detailed analyses have been produced that assess the relative importance of, or interaction between, these direct and indirect climate change impacts on biodiversity at large scales. Here, we apply a new integrated modelling framework to quantify the relative influence of biophysical and socio-economically mediated impacts on avian species in Vietnam and Australia. We find that socio-economically mediated impacts on suitable ranges are largely outweighed by biophysical impacts, but global shifts of production are likely to result in adverse impacts on habitats worldwide. By translating economic futures and shocks into spatially explicit predictions of biodiversity change, we now have the power to analyse in a consistent way outcomes for nature and people of any change to policy, regulation, trading conditions or consumption trend at any scale from sub-national to global. We present a novel framework for integrated macro-economic, land use, and biodiversity change modelling that permits quantitative analysis of questions critical to land use and biodiversity outcomes under broader socio-economic narratives, but also very specific policy 3 scenarios. We are now in a position to analyse the impacts of diverse domestic and international policy settings on land use and biodiversity, including changes to trade agreements and other economic shocks. Applying this new framework, we provide a first assessment of the relative magnitude of socio-economically and biophysically mediated climate change impacts on biodiversity in Vietnam and Australia. Climate change affects biodiversity through a multitude of pathways. There is pervasive evidence that climate change directly affects environmental conditions that are related to the climatic niches of many taxa, with the potential of significant shifts in their distributional ranges or even the total extinction of species (1,2). However, climate change also affects biodiversity through indirect human-mediated impacts: it drives the loss of livelihoods and displacement (3) and affects food and commodity production systems through its impacts on land productivity and human health (4,5) and environmental suitability for different land uses (6,7). Resulting global transitions of land use patterns are set to drive habitat conversion and may have dramatic impacts on biodiversity (8)(9)(10). While there are some examples of studies examining synergistic effects of land use and climate change on species (11,12), large-scale assessments of biodiversity change in response to climate change have tended to look only at direct impacts of climate change on biophysical conditions or habitat loss and fragmentation alone (8). Analyses that couple direct biophysical impacts on species with indirect socio-economic impacts via consumption, commodity, and land use change are sorely needed to fill important gaps in our knowledge of interactions between land use and climate change (10), to foster a more holistic understanding of the impacts of climate change, and to support the design of cross-sectoral adaptation and mitigation strategies (13). No single model of drivers of change in biodiversity and ecosystem services can capture all relevant dynamics at a high level of detail and there is an increasing awareness of the urgency to consider interactions between direct and indirect drivers of change under future scenarios to characterise prospects and management options for biodiversity and ecosystem services (13). Coupling demographic, economic and biophysical models has the potential to advance understanding and improve representation of synergies between direct and indirect drivers in biodiversity modelling, and to discover non-linear system behaviours that may not be apparent when considering drivers in isolation (13). Here, we contribute to the recent advances in integrated assessment modelling (14)(15)(16)(17) by applying an integrated modelling framework to compare the relative influence of direct biophysical and indirect socio-economic climate change impacts on the distribution and extent of suitable ranges for avian species in Vietnam and Australia . Recent advances in computable general equilibrium (CGE) modelling (18,19) bring unprecedented power to parametrise the impacts of climate change on commodity consumption and production patterns at very high commodity and temporal resolution across the global economy. We combine this economic modelling power with state-of-the-art land use change modelling to spatially downscale commodity demand changes caused by climate change(4) into changes of land use patterns. The spatial realisation of changing land use patterns varies with changes in the suitability of land for particular uses and is thereby also driven by climate change (7,20). Commodity demand changes are projected annually and land use predicted in 10-year time steps, producing decadal time-series maps of land use. Maps are integrated with climate change predictions into a biodiversity impact assessment using species distribution models (SDMs) (21)(22)(23)(24) . SDMs, fitted to current climate, land use, and other environmental variables (Supplementary Table 1) are extrapolated to conditions in 2070 under a range of climate and land use scenarios. Predictions of relative likelihood of occurrence are thresholded to examine changes in the ecologically suitable ranges for 1282 bird species in Vietnam and Australia (21)(22)(23)(24). Study area. We focussed our analysis on Vietnam and Australia because the countries provide unique socio-economic contexts, while hosting a similar number of bird species that are vulnerable, endangered or critically endangered (25,26). Due to the country's small size and limited number of occurrences, SDMs for Vietnam were built using data from a 30 x 30degree tile that comprises large parts of Southeast Asia. This enabled us to capture the occurrence of bird species present in Vietnam, across a much broader range of environmental variables, enabling better prediction of likely occurrence under future climates. Climate Change. We chose two alternative representative concentration pathways (RCP) that represent two extremes of the expected radiative forcing levels, RCP2.6 and RCP8.5 (27). For each pathway, we acquired the 2070 predictions of the first 19 bioclimatic variables downscaled as part of WorldClim Version 1.4(28) from 15 GCM of 5 (CMIP5) (29).
In order to account for variation of GCM predictions under different models, we determined the cell-wise first and third quartiles, as well as the medians of each of the 19 variables across the 15 GCM (Supplementary Table 2 for a complete list of the GCM used in this analysis). Main results were derived by predicting land use and species distributions under the medians of these variables. We also predicted both land use and species distributions under the first and third quartiles to approximate the range of outcomes for species across all 15 GCM (Supplementary Fig. 1). In CGE models, we included the parametrization of both climate change pathways proposed by Roson & Satori(4). We developed an inter-temporal Global Trade and Analysis Project (GTAP) model (30) to simulate changes in production under different climate change scenarios. CGE models use input-output-tables derived from national economic census data. These tables represent the inputs required in each economic sector to produce outputs and meet household and government demands (both nationally and internationally), which in turn is affected by prices and thus supply. Sectors are linked within each national economy, but also between economies. Producers in each country can produce various commodities to sell domestically (to household and government) and internationally to foreign households and governments. The households and governments generate their income from selling (to producers) productive input factors (land, capital, labour, etc.) and through taxes. In our version of GTAP, the total land area (land endowment) from which allocations are made to crop sectors (land requirements) can be changed in the baseline. Therefore, land supply is not necessarily fixed, as is the case in most other GTAP models. Estimations within GTAP are carried out relative to this baseline supply and we convert relative changes in agricultural land requirements to absolute changes in cropland by using their respective shares in the total harvested area for a by-sector-weighting of the average relative change of all classes. This weighted average change is applied on the current area under cropland to derive a future value. This means there is a direct proportional link between changes in land requirements and changes in the total area of agricultural land and the total area under agricultural land can change at the expense or to the benefit of other classes. Our inter-temporal GTAP model uses the GTAP 9 database (31), which is subdivided into 139 regions and 57 commodity sectors (31) and extends the GTAP model by replacing the recursive dynamic module of the current GTAP model with a forward-looking dynamic (or inter-temporal) module, where the producer can optimise profits overtime (32,33). More than just a trade model, the inter-temporal GTAP model allows optimal investment behaviours, in which producers in each country can adjust their decisions based on the impacts from both past and foreseeable future events. Agents in the model can react to future threats long before their full realisation (33). This makes the model a perfect tool for the simulation of future phenomena like climate change. Climate change impacts are modelled in our GTAP model following the work by Roson & Satori(4), in which impacts are realised as shocks to land supply and agricultural and labour productivity. The reduction in endowments of productive land and productivity negatively affect the production of commodities. Agricultural commodities are expected to be the most affected. With production shrinking more in some commodities than others, the price will adjust to balance the demand and supply of commodities. As a result, we will see a substitution effect between domestically produced products and their competitive imports along with a substitution effect in factors of production (such as land), balancing demands between sectors. Unlike the Kompas et al. (33) approach, which relied on a one-step simulation approach, here we apply a multi-step simulation approach allowing the shocks to be applied into smaller successive intervals combined with extrapolation techniques to further enhance the simulation accuracy (see Horridge et al. (34) and Pearson(35) for the details on multi-steps CGE solution methods). The solution of the inter-temporal GTAP model in this paper has been carried out within a parallel computing platform (19,36) with the use of PETSC (37)(38)(39) and HSL (40) libraries. sector at 2016 levels(41) (Supplementary Table 3) and calculated a weighted average change in land requirements for all cropland classes in each time step. Changes in urban land were estimated using estimates of urban population changes (42) and adjusting the amount of land under this class, assuming that urban population density remains steady through time. Future applications of this work will establish links between land-use classes related to forestry and livestock-raising, as has been demonstrated recently (17). We predicted land use maps under both pathways in 10-year time steps, using an R implementation (R package 'lulcc'( 43)) of the Conversion of Land Use and its Effects at Small regional extents (CLUE-S) model proposed by Verburg et al. (44). First, we determined the local suitability for different land uses through logistic regression of land use against the linear combination of a range of biophysical and socio-economic drivers in Generalised Linear Models (GLMs), from 15,000 randomly sampled pixels in each region (Supplementary Table 1 for a detailed list of data, Supplementary Figure 3 for effect sizes of predictors in each GLM). The selection of variables for land use suitability models was based on work by Verburg et al (45). ). Correlation analysis eliminated highly correlated predictor pairs (Spearmen's rank correlation coefficient  0.7), always keeping the predictor whose highest correlation with any other remaining predictor was smaller, to maximise independent information retained in the final set. The final predictor sets were checked against literature (46,47). We discarded a small number of predictors using cross-validated Lasso penalisation in the 'glmnet' R-package (48) and used the reduced predictor sets to build GLM and predict to future timesteps by interpolating GCM-predicted WorldClim variables in 10-year time steps under RCPs 2.6 and 8.5 (Supplementary Table 2 for used GCM models).
GLM predictions produced probability maps that represent the landscape's potential suitability for each land use class, under consideration of changes in the included variables. Transitions between classes were restricted according to a transition matrix that specified which transitions were possible (Supplementary Table 4). We specified conversion elasticities of each class (the amount by which land-use can shift without changing the total area it occupies) based on literature (43,44). To satisfy the CGE-projected changes in land endowments, changes were allocated in areas with the highest suitability for each land use, until estimated land area demands were met (49). Competition between land uses is handled in CLUE-S by allocating the land-use with the highest predicted suitability in a given cell, accounting for conversion elasticity and allowed transitions. We masked category I and II protected areas (7), precluding these areas from land-use changes. Since there was no CGE-modelled future demand for herbaceous ground vegetation and shrubland, as well as the forest classes, the overall amount of area allocated to those land uses was simply what was left over from the uses prescribed by agricultural and other demands. The proportional allocation between each of these residual categories was determined based on their mean predicted suitability in the landscape. All land use simulations were made using GCM-predicted first and third quartiles, as well as the medians of bioclimatic variables. Correlative species distribution models (SDM) can predict responses of species to changing environmental conditions by extrapolating from the covariate space in which they were observed (21)(22)(23)(24). The MaxEnt software package (50) (ver. 3.3.3k) was used to fit SDMs for 656 bird species in Australia and 739 bird species in Vietnam, using presence-only data from the (GBIF) (51). We filtered records to retain observations from 1950-2018 and records with more than or equal to 20 occurrence points across the model-fitting area. We included a range of biophysical, topographic and socio-economic predictors as well as land use (Supplementary Table 1). Correlation analysis eliminated highly correlated predictor pairs (Spearmen's rank correlation coefficient  0.7), using the same method as for land use predictors. We ensured through literature review that the final predictor sets were ecologically meaningful to avian species across taxa (52-55) at our aspired scale. We kept 9 predictors for Australia and 10 predictors for Vietnam, including 5 and 6 climate predictors respectively. Sampling bias is a pervasive issue particularly affecting presence-only data that is often sampled opportunistically. We accounted for sampling bias by estimating the intensity of sampling effort in response to demographic and topographic variables (56), and using this map of sampling effort to probabilistically select background points. By selecting background points proportional to sampling bias the effect of sampling effort on the location of presence records is largely eliminated as a form of bias (57). Variables used in the bias models were protected area status, distance to roads, distance to built-up areas and roughness. Predictions were made using the estimated quartiles and medians of bioclimatic variables and the according land use maps that were also predicted under quartiles and medians. We controlled overfitting by determining the permutation importance of predictors and dropping predictors with a value < 1%. Test AUC were estimated via 5-fold cross-validation of each model and final models built on all available records. Species for which only uninformative models were fitted (AUC < 0.7) were excluded (58). We recorded the log ratio of the respective number of cells with relative likelihoods predicted above MaxEnt's MaxSSS threshold (59) (where the sum of model sensitivity and specificity is maximised) between present (2018) and future time step (2070) as a measure of change . In Australia, we constrained this change estimation for each species to bioregions containing records of the species, and adjacent bioregions (60). All data preparation and modelling for land use and SDMs was conducted in R(61), using packages 'lulcc' (43) for land use simulations and 'dismo' (62) for MaxEnt models. All analyses and spatial predictions of the land use model and SDM were performed at 0.5 arcminute resolution; approximately 1 km at the equator. SDM building and predictions were computationally expensive and required up to 50 GB of working memory on 12 parallel cores. Direct biophysical impacts dominate changing range sizes. For birds in both regions, we forecast major declines in ecologically suitable ranges, with severity of loss scaling with the severity of climate change . Under RCP 8.5, a much higher number of species would be expected to experience decreases of more than half of Characteristics Curve values (AUCs) (63) of models in the two regions as well as the respective number of models (n) retained (AUC > 0.7) (58). AUC provides a measure of a model's discriminatory performance in terms of how well test predictions discriminate between occupied and unoccupied locations (58,63) . d, e -Fractions of models in which a predictor was used. Full names and definitions of all predictors can be found in Supplementary The total output of most agricultural crop sectors in both regions was predicted to decrease more with increasing climate change. In particular, in Vietnam, sectors such as oil seeds and plant-based fibres shrink by up to 20% under RCP 8.5 . The land requirements for each sector generally increase in proportion to the overall output of each sector. This is due to climate change impacts on crop yields as parametrised in the CGE-model: reductions in land productivity mean that more land is required to maintain sector outputs. Accordingly, in both countries, even while total outputs tend to decrease, land requirements of agricultural sectors remain approximately the same, or increase slightly . The changes in land requirements for crop lead to an increase in cropland of <0.5% of the total land area in both regions under RCP 8.5 and a very slight decrease in Australia under RCP 2.6 . Increases in urban land in both countries were modelled on FAOSTAT estimates of urban population growth (41).
In Australia, land use changes occur locally and are concentrated in coastal areas along the north-east, south and west of the continent, although some changes also occur further inland . In Vietnam, land use change is higher overall, with a particular concentration of change in the central-southern and northern coastal areas of the country, that also approximately coincide with the country's major river deltas . Given that the distributions of most species are constrained, aggregated, and not random, small percentage changes in land use at the national scale still have significant impacts on some species locally . For example, species losing more than 10% of their currently suitable range under indirect impacts of RCP 8.5 in Vietnam include the vulnerable chestnut-necklaced partridge (Arborophila charltonii) and the near-threatened yellow-billed nuthatch (Sitta solangiae). These declines are highly localised and predominantly occur in the centre-south of the country . Direct climate change impacts are more severe: 324 and 362 species lose at least 10% of their suitable ranges under direct impacts of RCP 2.6 and RCP 8.5 respectively, with areas particularly affected across taxa under RCP 8.5 in the northern highlands and the central eastern parts of the country . Among the species losing more than 95% of their current suitable range under the direct impacts of RCP 8.5 are the Chinese thrush (Turdus mupinensis) and the critically endangered white-rumped vulture (Gyps bengalensis) loses more than 5%. A higher number of species are affected by the direct impacts of climate change, with areas predicted to suffer particularly high suitable range declines along the southern and eastern coasts, the southwest and the southeast of the continent. In Australia, 188 and 230 species are expected to lose more than 10% of their suitable range under RCP 2.5 and RCP 8.5 respectively. Amongst the Australian species losing more than 95% of their suitable range under the direct impacts of RCP 8.5 are the kalkadoon grasswren (Amytornis ballarae) and the Australasian pitpit (Anthus Australis) and a number of other species now categorised as of least concern . This highlights the potential dangers of climate change to species that we do not yet consider under threat, but for which extinction debts are accruing (64). We predict economic change under two climate change scenarios, keeping all other aspects of the global economy at current baseline values. This way, we could capture and isolate the effect of climate change on the economy. However, this approach omits other socioeconomic processes that could affect supply and demand, such as population growth, changes to economic growth, energy efficiency, and shifts in social demands. These other factors may impact habitat and biodiversity through agricultural expansion, deforestation or urbanisation. While this study was designed to assess the net effects of direct and indirect climate change impacts on species as a first case study introducing our integrated assessment framework, these factors will be incorporated in future iterations that include an even more comprehensive CGE parametrization (i.e. full CGE baseline scenarios with socio-economic pathway narratives (69) and integration of climate models in CGE analysis) and through improvements to current CGE methods by including, for example, stochastic effects of natural disasters in the CGE modelling. Biodiversity model assumptions may underplay potential indirect impacts. We chose not to produce SDMs for species with less than 20 occurrence records, to avoid the inflation of AUC for range-restricted species and species with very low prevalence (70) and to assure sufficient discrimination between presences and background points (71,72). Rare or spatially restricted species can be more vulnerable to localised impacts such as habitat loss (73), but these effects are difficult to capture when biodiversity data are poor. We assumed unlimited dispersal ability in Vietnam and dispersal ability constrained to bioregions adjacent to those containing observation records in Australia. Disconnected patches of potential habitat outside of observed ranges (but within adjacent bioregions in Australia) were counted in future predictions, regardless of whether those areas were functionally linked (by suitable or traversable habitats) to the observed range and thus within the dispersal range of species. This may lead to an over-estimation of habitat utilisation and a commensurate underestimation of both direct and indirect impacts of climate change, particularly for taxa with a low dispersal ability that rely on small pockets of habitat within their range and are unable to reach disconnected patches of potential habitat. The importance of connectivity as a key component of habitat structure is well known (74) and crucial for population viability in many species with low dispersal ability (75,76). The parametrization of species' dispersal ability and explicit modelling of landscape structure in response to land use change would allow for the inclusion of these fragmentation effects. This may be particularly important when our framework is extended to non-avian species . Exported biodiversity impacts. While we found that total agricultural sector outputs decrease in both Vietnam and Australia, decreases in land productivity mean that land use in production for some agricultural commodities were predicted to increase slightly. We assumed a global economic equilibrium in which commodities can be substituted through trade between regions, thus implying that global demand for land-based commodities is serviced by regions that benefit from a comparative advantage under climate change. Where comparative advantage is due to increases in land productivity (land use efficiency), additional land may not be required to increase outputs. However, where this advantage is due to other economic mechanisms and not driven by the cost of converting additional land for production, more land may be allocated to agricultural or other commodity production, increasing habitat loss. For example, in Canada, our CGE model predicted an increase of wheat sector output by over 37% under RCP 8.5, while land endowments increase by only 14% due to increases in land use efficiency. In India, wheat output is estimated to increase by 8% under RCP 8.5, while land endowments to the wheat sector increase by 6 %, suggesting much lower land use efficiency in India than in Canada (see Supplementary Fig.
2 for a global, country-wise mapping of projected changes in sector outputs and land endowments of the wheat sector). Despite lower land use efficiency, wheat production in India still grows, because growth is economically feasible as long as it is not limited by factors arising from the sector's context in domestic and international markets. In both countries, increases in land use lead to agricultural expansion, but in Canada more wheat can be produced per unit land and areas lost to wheat farming are likely to be much smaller per produced unit than in India. Nonetheless, if wheat production occurs in parts of Canada that were previously in, for example, natural prairie, then significant biodiversity losses may occur. Our framework provides in-depth insight into the links between sectors and regions and allows for a better understanding of global shifts in land requirements, enabling the finescale identification of hotspots for production, agricultural expansion and ultimately habitat destruction under consideration of the global economic processes. In this first implementation of our framework we could capture and quantify principal relationships between climate change, the global economy, land use and avian habitat. Future uses of our approach could include regional and global biodiversity assessments following individual policy shocks, such as the introduction or abolishment of taxes or international trade deals, or could seek to capitalize on existing narratives of socio-economic futures and climate change pathways (so-called Shared Socio-economic Pathways) (69) to parametrise climate adaptation policies, sustainable development goals and other aspects of socio-political transitions within the CGE modelling. Expanding consideration of biodiversity to include non-avian taxa and explicitly dealing with the role of connectivity and dispersal will enable a more comprehensive assessment of biodiversity impacts under socio-economic change. A key feature of our approach is that it provides opportunity to downscale country-level commodity demands to spatial explicit land use changes and biodiversity impacts, enabling a more meaningful analysis of the habitat and biodiversity implications of economic shocks or the implications of trade than have previously been possible. Better integration of models and scenarios of biodiversity is required to guide evidencebased climate adaptation strategies and to chart progress toward sustainable development goals (77) This work 04795. The authors acknowledge contributions made throughout the research phase by B. Bryan, P. Lentini, H. Kujala and J. Elith. This work 04795. The authors acknowledge contributions made throughout the research phase by B. Bryan, P. Lentini, H. Kujala and J. Elith. SK led the modelling with contributions from HPV, TK, NG, NC, BW. BW, TK and SK conceptualised the analysis framework. SK led the writing of the manuscript with contributions from all authors. The authors declare no competing interests. Sources for data used in land use and species distribution modelling are listed in Supplementary Information. Direct download links to these data sets are available in the code repository accompanying this study (see below). We provide outputs of the CGE, landuse and species distribution models in a data repository (DOI: 10.5281/zenodo.3517724). The GTAP database that underpins GGE modelling is available from GTAP under license. However, we provide input files and commented CGE modelling code that contains details of parameters settings for global economic models and detailed commodity demand output tables for each of the scenarios modelled upon publication. Components of the GTAP database can be provided on request. All R-code for land use and species distribution modelling is available online (DOI: 10.5281/zenodo.2703618). The tablo-code for the CGE modelling will be provided in the code repository accompanying this study upon publication. Authors: Natasha C. R.