Potential for invasion of traded birds under climate and land-cover change.

Humans have moved species away from their native ranges since the Neolithic, but globalization accelerated the rate at which species are being moved. We fitted more than half million distribution models for 610 traded bird species on the CITES list to examine the separate and joint effects of global climate and land-cover change on their potential end-of-century distributions. We found that climate-induced suitability for modelled invasive species increases with latitude, because traded birds are mainly of tropical origin and much of the temperate region is 'tropicalizing.' Conversely, the tropics are becoming more arid, thus limiting the potential from cross-continental invasion by tropical species. This trend is compounded by forest loss around the tropics since most traded birds are forest dwellers. In contrast, net gains in forest area across the temperate region could compound climate change effects and increase the potential for colonization of low-latitude birds. Climate change has always led to regional redistributions of species, but the combination of human transportation, climate, and land-cover changes will likely accelerate the redistribution of species globally, increasing chances of alien species successfully invading non-native lands. Such process of biodiversity homogenization can lead to emergence of non-analogue communities with unknown environmental and socioeconomic consequences.

INTRODUCTION

Biological invasions are one of the top‐five greatest threatening factors affecting native biodiversity and ecosystem functioning worldwide (Pyšek et al., 2020; Vilà & Hulme, 2017). Risks associated with invasive species have also been documented in the primary economic sectors, such as agriculture, forestry, fisheries (Paini et al., 2016), and also in human health (Stoett et al., 2019). Invasive risk assessment protocols, aimed at listing and ranking species by their likelihood to invade and cause harm, have become a main tool for supporting policy‐makers and managers alike in their decisions regarding how to handle these species whenever they cause damage (Vilà et al., 2019). Invasive risk protocols may be implemented such that they provide the full analytical package of risk. However, they can also cover subsets of the protocol, for example, by focusing on invader impacts (Strubbe et al., 2019). While the exact method and approach for risk assessments will differ with the specific objectives and scope of the exercise (Vanderhoeven et al., 2017), quantifying species' overall invasion risk will generally include spatially explicit projections of changes in environmental suitability for species establishment and spread. Devising adequate strategies for countering the negative effects of invasive species, therefore, requires the use of predictive models for estimating the spatial dynamics of area‐suitability to species invasion (Araújo & Peterson, 2012).

While trade has been historically constrained by geography, a globalized world, preferential trade flows rather than distance correlate with invasion risk (Menchetti & Mori, 2014; Ribeiro, Bingre, et al., 2022). Besides propagule pressure—or the frequency with which areas are exposed to potential colonization events by non‐native species—studies have documented that suitability to invasion by non‐native birds increases with human disturbances (Cardador & Blackburn, 2020), climate similarity between native and invaded ranges (Stuart et al., 2004), trade fluxes (Reino et al., 2017), or several combinations of more than one of these as well as other factors (Cardador & Blackburn, 2020; Vicente et al., 2019).

Although the effects of changes in land cover (Reino et al., 2018), climate (Thuiller et al., 2011), and their interactions (Fordham et al., 2018; Jetz et al., 2007; Taheri et al., 2021) on native species distributions have been extensively investigated, the associated question of how the two drivers combined and in isolation will alter the future suitability of species invasions has been largely understudied. Climate change, for example, is a well‐known driver of species redistributions (Pecl et al., 2017; Pereira et al., 2010) and its impacts on invasive species have been discussed in general terms (Walther et al., 2009). Yet, not many studies have examined how different climate change scenarios might affect changes in invasion suitability. Those that did tend to find that impacts are non‐uniform across species (e.g., Baquero et al., 2021; Bezeng et al., 2017; Peterson et al., 2008) with, perhaps, the exception of climate change extremes that are generally thought to enhance invasions by promoting the transport of propagules into new regions, and by decreasing the resistance of disturbed native communities to establishment (Diez et al., 2012). Likewise, many studies have shown that intensification of land‐use and associated changes in vegetation cover are critical factors governing distributional dynamics of invasive species (Cardador & Blackburn, 2019).

Few studies, however, have prospectively examined how future changes in land cover might affect biological invasions. One reason is that modelling the effect of land‐cover changes on species distributions poses several technical problems that are more difficult to address than the familiar models examining climate effects (Thuiller, Araújo, & Lavorel, 2004). A recent review on climatic niches of 434 invasive species from 86 studies found that only about 1% of the studies accounted for species' land‐cover preferences (Liu et al., 2020). Another review of invasive bird species distribution modelling, showed that nearly half of the studies combined both climate and land‐cover predictors (Liu et al., 2020) (12 out of 26 studies). However, they were mostly investigating suitability of range expansion in the present without investigation of future projected land‐cover and climate change effects. The few studies that estimated changes in invasion suitability based on the joint analysis of land‐cover change and climate change, have generally found that the combined changes increase invasion suitability when compared with climate change alone (Di Febbraro et al., 2019). Given the magnitude of land‐cover changes in the past decades (Di Febbraro et al., 2019) and the increased changes projected for the 21st century (Hurtt et al., 2020), enquiring how land‐cover change might interact with climate change in altering the regional suitability to invasion by introduced species is a critical first step towards shaping invasive species policy and management.

We provide a comprehensive global analysis of the joint effects of future climate change and land‐cover change on the global distribution patterns of suitability to invasion by several commonly traded birds. Specifically, using a large ensemble of 12 species distributions models and 17 global climate models, we fitted over half‐million models for 610 internationally traded birds listed within Appendix II of the Convention on International Trade in Endangered Species (CITES). The purpose of this list is to identify species of conservation concern that have been commonly traded and require restrictions to reduce pressure on their capture for trade. Nevertheless, the high market demand for many of these species has led to several of them having been transported and established as non‐native to the host regions; a first and necessary step before developing invasive behavior. The models were initially projected considering the SSP2‐45, SSP3‐70, and SSP5‐85 emissions scenarios. While scenario SSP2‐45 closely matches the Paris Agreement goal of “limiting global temperature increase to well below 2°C, while pursuing efforts to limit the increase to 1.5°C", few researchers consider it as plausible outcome of current political trajectories with extreme scenario SSP5‐85 being more likely (Schwalm et al., 2020). Hence, climate change analyses were focused on middle ground scenario SSP3‐70 and the extreme scenario SSP5‐85, with results for scenario SSP2‐45 being provided as supplementary information. Based on literature review, the selected bird species were associated with the predominant known land‐cover types where they inhabit within their native and invaded ranges. Dynamic models of trends in land‐cover change under SSP3‐70 an SSP5‐85 (Hurtt et al., 2020) were then used to infer the extent of the projected changes in habitat loss and gain for each species' preferential land‐cover type towards the end of the century. Finally, the joint examination of trends in climate and land‐cover changes were undertaken to identify coldspots and hotspots of changes in CITES‐listed traded birds' invasion suitability around the world.

MATERIALS AND METHODS

Data

Species data

We obtained the list of bird species known to be traded from Appendix II of CITES (Convention on International Trade in Endangered Species) database (available online at https://trade.cites.org) (Figure S1). The CITES Trade Database is managed by UNEP‐WCMC (United Nations Environmental Program, World Conservation Monitoring Centre), and records trade between 183 signatory countries, holding over 18 million records of trade in CITES‐listed wildlife. For birds, which are among the most heavily traded taxa worldwide (Reino et al., 2017), CITES Appendices cover about 1700 out of ca. 2600 bird species known to be traded internationally (~65%; Reino et al., 2017). CITES trade data are often used as broadly representative snapshots of the global legal trade in wildlife (Can et al., 2019), and covers some of the most invasive bird genera, such as nearly all of Psittaciformes, while capturing most of the trade volumes of Passeriformes as well (Cardador et al., 2016; Reino et al., 2017). Although created to list species that “may be affected by trade and become endangered if trade is left unregulated”, CITES Appendix II includes 130 bird species that were previously recorded as introduced in the wild, 50 of which have become established in non‐native regions (Dyer, Cassey, et al., 2017). However, the CITES trade database relies on information communicated by governments, which are not free of errors or biases. For example, governments might occasionally fail to correctly report transactions, species might be misidentified, or traded amounts poorly estimated (Reino et al., 2017). Despite caveats, CITES represents the only global legally binding convention addressing international wildlife trade in a structured and verifiable manner, constituting a valuable source of information to assess the relationships between conservation risks and international trade (Hierink et al., 2020; Phelps et al., 2010). CITES Appendix II thus provides a good starting point for examining potential non‐native range dynamics of commonly traded species.

We cross‐checked the list with birds in the comprehensive Copenhagen Global Avian Distributional Database mapped at a 1 × 1 latitude‐longitude grid (Rahbek et al., 2012). Maps represent a conservative extent‐of‐occurrence of the native breeding ranges based on museum specimens, published sight records, and spatial distribution of habitats between documented records, which have subsequently been validated by ornithological experts. We retained species for modelling with at least 15 records. Our final list was composed of 610 bird species.

Climate variables

We obtained 19 bioclimatic variables from WorldClim (version 2.1) related to temperature and precipitation at 10 min spatial resolution for both the baseline (1970–2000) and future period (2060–2080) (Hijmans et al., 2005). For the future, we considered the shared socio‐economic pathway scenarios (SSP2‐45 SSP3‐70 and SSP5‐85), encompassing a lower, middle, and upper boundary of the range of scenarios; because the lower boundary is now considered unlikely, most analysis focus on the middle and upper boundary with results for models assuming SSP2.45 being provided in the supplementary material. Future climate projections were derived from an ensemble of 22 different global circulation models (GCMs): ACCESS‐CM2; ACCESS‐ESM1‐5; BCC‐CSM2‐MR; CanESM5; CanESM5‐CanOE; CMCC‐ESM2; CNRM‐CM6‐1; CNRM‐CM6‐1‐HR; CNRM‐ESM2‐1; EC‐Earth3‐Veg; EC‐Earth3‐Veg‐LR; GISS‐E2‐1‐G; GISS‐E2‐1‐H; INM‐CM4‐8; INM‐CM5‐0; IPSL‐CM6A‐LR; MIROC‐ES2L; MIROC6; MPI‐ESM1‐2‐HR; MPI‐ESM1‐2‐LR; MRI‐ESM2‐0; UKESM1‐0‐LL. To make climate projections consistent with the species data, climate data were aggregated into 1‐degree grid cells by averaging climate scores originally at a finer resolution to the coarser resolution grid used for species distributions data.

Depending on the number of presence records, we selected a maximum of five uncorrelated variables for each species using a hybrid approach consisting, firstly, in examining pairwise variable correlations and discarding clearly collinear variables and, secondly, using of the variance‐inflation factor (VIF) test (Naimi et al., 2014) to examine if variables were strongly correlated thus incurring in familiar problems of multicollinearity (Dormann et al., 2013; Segurado et al., 2006). In addition, we quantified the importance of all variables to explain each species' distribution, measuring the performance of a set of species distribution models fitted using the individual variable as predictor. Using a stepwise procedure, we then excluded the less important variable from each pair showing the highest pairwise correlation (>0.7). In addition, variables with VIF >10 were excluded from the stepwise procedure. In the end, we selected the five variables with the highest individual importance for each species.

Climatic niche modelling

We analyzed the climatic niches of the 610 selected traded bird species for which we fitted species distributions models using the sdm R package (Naimi & Araújo, 2016). Using 12 machine learning algorithms, a bootstrapping resampling procedure with five replications, and a fully factorial combination of the selected predictor variables, we fitted a maximum of 960 models per species that resulted in ~585,000 models for all species. The fully factorial exploration of the variable combinations allows understanding the importance of variables to explain species distributions. Since some of the chosen modelling algorithms required a minimum of two variables, all the combinations with at least two variables were considered in the fully factorial exploration. Furthermore, the maximum number of variables in the combinations were defined according to the sample size (i.e., number of presence records) for a species. A major problem with modelling species with a high number of predictor variables in relation to low sample size is that it can lead to a model overfitting (Fielding & Bell, 1997). To avoid an overfitting, no more than n/5 predictors were included in the final model (where n is the total number of species presence records), thus for a species with at least 25 presences we examined the factorial combination of up to 5 predictor variables.

Given that the available species distributions data cannot be strictly interpreted as presence and absence, we generated pseudo‐absence (background) records for each species separately. The extent of the study area can also affect model results (Thuiller, Brotons, et al., 2004), thus the region for drawing background samples should include, but not exceed, all areas that are accessible to the bird species (Araújo et al., 2019). The zoogeographic regions where the occurrence records of species are located were considered to represent their native territory in a broad sense. These regions are defined not just based on species co‐occurrences but also on knowledge about their shared phylogenetic history. That is, biogeographical regions identify faunas that have not been in contact during previous historical climate changes. Hence, the assumption is made that a species whose native range overlaps with the biogeographical region would have either had access to currently unoccupied areas of that region, or at least coexisted with faunas that were distributed more broadly across the region. We used the updated map of the vertebrate zoogeographic regions of the world to specify the native area for each species (Holt et al., 2013). The background records were then generated within the study area for each species, using a sample size equal to 70% of the 1‐degree grid cells within the area. Sample size was thus similar for all species.

To develop climate niche models, with fitted 12 commonly used methods with the sdm R package (Naimi & Araújo, 2016), including the generalized linear model (GLM), generalized additive model (GAM), classification and regression trees (CART), boosted regression trees (BRT), random forests (RF), multiple discriminant analysis (MDA), support vector machine (SVM), multi‐layer perceptron neural networks (MLP), maximum entropy (Maxent), multivariate adaptive regression splines (MARS), maximum likelihood estimation for presence‐only (Maxlike), and bioclimatic envelope (bioclim). Resampling by bootstrapping with five replications was used to generate the training and test datasets. Bootstrapping uses sampling with replacement to draw a dataset with the same sample size as the original dataset used for training the models. The records that are not selected in the training dataset are then identified and used for evaluating the models (i.e., test dataset). For each variable combination, we fitted the 12 ecological niche models for each replication and evaluated them for their performance using the test dataset. We used the area under curve (AUC) of receiver operating characteristic (ROC) plot and the true skill statistic (TSS) to measure the predictive performance of models (Fielding & Bell, 1997). A ROC curve plots sensitivity values (true positive fraction) on the y‐axis against ‘1 – specificity’ values (false positive fraction) for all thresholds on the x‐axis. AUC is a threshold‐independent metric that varies from 0 to 1 and provides a single measure of model performance. AUC values under 0.5 indicate discrimination worse than chance; a score of 0.5 implies random predictive discrimination; and a score of 1 indicates perfect discrimination. TSS is calculated as “sensitivity + specificity ‐1” and ranges from −1 to +1, where +1 indicates perfect agreement, a value of 0 implies agreement expected by chance, and a value of less than 0 indicates agreement worse than chance.

For each species, we excluded any model that performed worse than a set predictive performance threshold (i.e., AUC < 0.7). The remaining models were used to characterize the potential distribution of the species in both current and future times. For this purpose, each individual model was used to predict the potential distribution of a species into the climatic variables for the current time and also to project the distribution into the future time.

Ensemble forecasting of climate suitability of traded CITES‐listed birds

For each species, we used an ensemble of models (Araújo & New, 2007) to calculate and predict a consensus potential climatic distribution of each species in the baseline period across the globe. Consensus among ensembles of models has been shown to reduce uncertainty and increase predictive accuracy of bird species distributions models under climate change (Araújo et al., 2005) and to approach predicted ability expected even when considering more complex mechanistic models (Fordham et al., 2018). Consensus was achieved through AUC‐weighted mean across all models (Garcia et al., 2012). The same approach was used to project the distribution into the future, representing a consensus climatic niche model across GCMs. To infer individual species suitability scores across the non‐native ranges, two additional steps were undertaken. First, for both current and future times, the estimated suitability scores over pixels located in native geographical areas (the biogeographical regions used as the study area) were set to 0. Second, pixels with climatic conditions beyond the range of current climatic conditions in native geographical areas were also set to 0. For the latter, we found the range (minimum and maximum) of the climatic variables that were selected to train ecological niche models for each species over the study area. We then identified the pixels with values out of the range for at least one of the selected climatic variables. Resulting maps in the current and future times were then inferred as the change in climate suitability for the occurrence of each species. We stacked the maps for all the 610 species and calculated the mean over all the values at each pixel to characterize the overall likelihood that modelled traded avian species would find suitable conditions for spread in both current and future times. The subtraction of the potential modelled distributions at the future time from the current time at each pixel characterized the changes (delta) in likelihood of expansion due to climate change.

The workflow for fitting and projecting ensemble models is described in Figure S2, and the R script for reproducing the analysis is provided in the supplementary materials 1 and 2.

Land‐cover change analysis

We used the harmonized land‐use dataset to examine how changes in land cover might affect future avian redistribution (Hurtt et al., 2020). The dataset uses several models to characterize fractional land‐cover patterns in a 0.25° × 0.25° grid resolution between 850 and 2100 (Hurtt et al., 2020). To make the land‐cover change analysis consistent with climate‐change analysis, we used scenarios SSP3‐70 and SSP5‐85, consistently with the climate change analysis, built with the Asia‐Pacific Integrated assessment Model (AIM), and REMIND (REgional Model of Investment and Development) and MAGPIE (Model of Agricultural Production and its Impacts on the Environment), respectively, that simulate a pattern of land‐cover changes compatible with scenarios chosen (Kriegler et al., 2017). We generated outputs for 2000 and 2080 to match the end‐period of the baseline and the future periods respectively. The land classes considered were primary forested land (primf), primary non‐forested land (primn), potentially forested secondary land (secdf), potentially non‐forested secondary land (secdn), pasture resulting from the combination of classes of managed pasture (pastr) and rangelands (range), urban land (urban), annual crops resulting from the combination of the classes of C3 and C4 annual crops (c3ann and c4 ann) with C3 nitrogen‐fixing crops (c3nfx), perennial crops resulting from the combination of C3 and C4 perennial crops (c3per and c4per). For visual inspection of the distribution of land cover classes see Figure S3. Aggregation of some land‐cover classes was necessary because the disaggregated classes were not discriminative of the bird preferences.

We identified the land‐cover classes that are suitable for each one of the 610 traded species considered based on expert knowledge of LR, JR, and DS supported by analysis of the “Birds of the World” (Billerman et al., 2020) (see supplementary material 1). Land‐cover associations were mainly based on described habitat preferences, though information on species conservation status was used as well, especially to assess tolerance to human‐modified habitats such as urban areas, pastures, annual crops, and plantations, in tandem with visual inspection of “eBird bird” occurrence maps and the global distribution of land cover as given by Figure S3. For each species, we then examined changes in the availability of suitable land‐cover classes across the grid cells of the world. Then, by quantifying the proportion of suitable and unsuitable land cover for each species, in every grid cell, in the two time periods, we were able to estimate if suitability of land cover for individual species is expected to remain stable, decrease, or increase towards the end‐of‐the‐century.

RESULTS

Most traded CITES‐listed birds are of tropical origin (see Figure S1). Unsurprisingly, stacked modelled climate suitability for these species shows high scores around tropical and adjacent subtropical areas followed by Mediterranean and semi‐arid areas, then decreasing towards high‐latitude temperate and boreal climates. This general pattern of decreasing climate suitability with latitude was recorded for both the baseline and end‐of‐century periods (Figure 1a,b). However, when temporal trends are examined in a detail, a worsening of climate conditions for traded birds emerges in the tropics, semi‐arid and warmer areas of the temperate zone (e.g., Mediterranean region, southern half of USA), whereas an improvement of climate conditions for these species emerges in colder temperate and boreal zones, including the southern tip of South America, southeastern coast of Australia, New Zealand, northern USA, Canada, and much of central and northern Eurasia (Figure 2a). Such trends are true for all climate change scenarios considered (for scenario SSP2‐45 see Figure S4) but are more pronounced for scenarios involving greater emissions.

FIGURE 1

Stacked mean modeled climate suitability for selected 610 traded bird species across invadable ranges in the baseline period (a); and in the future, for the socioeconomic pathways (SSPs) scenarios of SSP3.70 (b) and SSP5.85 (c). Map lines delineate study areas and do not necessarily depict accepted national boundaries.

FIGURE 2

Changes in climate suitability of traded bird species towards the end of the century. (a) Projected change in climate suitability of CITES species across invadable ranges between the baseline period and the future. Warm colours indicate increases in suitability to invasion, cold colours indicate decreases in suitability. Map lines delineate study areas and do not necessarily depict accepted national boundaries. (b) Magnitude of potential flows in bird invasion between biogeographical regions of the World, sensu Holt et al., 2013. Left nodes are native ranges, right nodes are invaded ranges. The size of left nodes indicates the number of traded bird species available in the region, whereas the size of right nodes indicates the number of species moving into that region. The thickness of the arrows represents the number of species projected to find suitable climate across the invadable non‐native region of the world by 2060–2080 (the left graph corresponds to SSP370, and the right graph corresponds to SSP585). Red arrows represent increases in numbers of species projected to have increased climate suitability in the invaded range with regards to the baseline period, whereas blue arrows represent projected decreases in suitability.

When modelled climate suitabilities are aggregated at the level of biogeographic regions, disregarding the actual fluxes of birds driven by trade, most potential climatically driven fluxes of traded birds should be among regions with similar climates (Figure 2b). Since most birds are of tropical provenance, all other things being equal, most climate‐similarity‐induced fluxes should be among regions with tropical climates. However, climate change is projected to lead to homogenization of the differences in climate suitabilities between temperate and tropical regions. A trend of increased climate suitability for the spread of bird species native from the Afrotropic, Australasia, the Indo‐Malay, and the Neotropics is thus expected towards the Nearctic and the Palearctic regions (Figure 2b). This trend is accompanied by decreased climate suitability for across tropical bird colonization apart from Australasia which shows increased climate suitability towards the Neotropics. Again, trends that are qualitatively consistent across emissions scenarios.

Suitability of land‐cover and land‐cover change can also affect avian range expansions, but unlike climate suitability and climate change, the pattern lacks a clear latitudinal gradient. Nearly 70% and more than 60% of the CITES‐listed birds have a strong association with primary forests and secondary forests, respectively, followed by non‐forested primary lands (Figure 3). Because of this association with forests in their native regions, the boreal region stands out as having highly favorable land cover for many traded bird species as this region is largely covered by forests (Figure 4). For the same reason, equatorial areas of Africa, such as the Congo basin countries, southeast Asia, such as Papua New Guinea and Malaysia, and the Amazon basin are also highly suitable for the traded birds. Prior adaptation to human‐modified environments is also known to positively influence avian invasion success (Cardador & Blackburn, 2020; Strubbe et al., 2015) causing strongly urbanized regions, such as parts of the eastern United States of America and southern China to be prone to avian invasions. Psittaciformes, a heavily traded taxon (Reino et al., 2017) almost exclusively nesting in tree hollows, exemplify how traded birds can benefit from such land‐cover changes, as some parrots and parakeets released in urban environments readily find suitable nesting trees in parks and large gardens, allowing them to first establish and then spread into more natural environments (Hernández‐Brito et al., 2020).

FIGURE 3

Traded‐bird‐species‐land‐cover associations. The bars represent the number of birds associated with the different land‐cover categories. The numbers are percentages of species associated with the different types of land cover. Because species can be associated with more than one type of land cover the sum of percentage scores is greater than 100.

FIGURE 4

Estimated land‐cover induced suitability to invasion by selected CITES‐traded bird species in 2020 and 2080. (a) Land cover suitability for invasive species in 1990 and 2080, with two emissions scenarios; (b) Delta land‐cover suitability obtained by subtracting suitability scores in 2080 to scores in 1990. Map lines delineate study areas and do not necessarily depict accepted national boundaries.

When the effects of future projected land‐cover changes are examined (Figures S5 and S6), the estimated suitability for CITES birds across the boreal region is maintained (Figure 4). The highest increase is expected along the Brahmaputra River in the Chinese Himalayas, due to projected expansion of secondary forest replacing existing mountainous pastures (Figures S5 and S6). The Congo basin countries are projected to lose suitability for the bird species examined across much of their central area, owing to primary forest conversion into pasture and to a lesser extent to secondary forest. Eastern Africa is also losing favorability for birds, given the conversion of non‐forested primary land (i.e., savanna) into non‐forested secondary land and annual crops. Likewise, reductions in suitability of land cover are projected across southern and eastern areas of Brazil, broadly coinciding with the Atlantic Forest and the Cerrado biomes, both undergoing intensive forest‐to‐cropland conversion. As for climate, trends are consistent across scenarios but with intensification of changes with the scenario SSP5‐85.

The effects of interactions between climate change and land‐cover change on the traded invasive species distribution dynamics are difficult to quantify. However, by comparing projected climate change and land‐cover change effects together, we uncover potential hotspots and coldspots of changes in suitability for invasion across two dimensions of environmental change (Figure 5). Much of the temperate region across the northern hemisphere witnessed improvements in the projected climate and land‐cover suitability for traded bird species. Near the tropics, the pattern tends to be reversed because deforestation and climate change reduces the amount of suitable habitat for the species. Our modeling has thus identified areas where both climate and land cover (as a surrogate for habitat) would indicate the suitability of establishment for the traded bird species. Yet, the extent to which introduced birds will be able to colonize and persist in these regions is also dependent on their dispersal capacity and biotic interactions with other species (Araújo & Peterson, 2012; Soberón, 2010); an issue that we do not address herein.

FIGURE 5

Covariation in projected trajectories of climate change and land‐cover change effects on selected CITES‐traded bird species invasion suitabilities by 2080. Positive values in both axes represent increased suitability for invasion with regards to the baseline period, while negative values represent decreased suitabilities. Map lines delineate study areas and do not necessarily depict accepted national boundaries.

DISCUSSION

The past five decades witnessed a large increase in the world's human population only matched by the intensification of global trade (Kentor, 2001). The globalization of trade boosted biological invasions, either as a direct consequence of animals being intentionally moved out of their native ranges for commercial purposes (Reino et al., 2017), or as an indirect consequence of species being unintentionally moved and introduced in new territories (Costello et al., 2007; Vall‐llosera & Cassey, 2017). While future dynamics of biological invasions are uncertain and contingent on sociopolitical (Ribeiro, Bingre, et al., 2022) and human behavioral (Ribeiro, Araújo, et al., 2022) aspects that are difficult to anticipate, there are clear trends that can inform prospective studies. For example, most recent fluxes of traded birds follow a latitudinal gradient with most invasive bird species being of tropical origin. The temperate region has been a major receptor of traded birds due to historical and socioeconomic reasons, as demand for tropical birds has always come mainly from affluent consumers in the Nearctic, Western Palearctic and Eastern Palearctic (Dyer, Cassey, et al., 2017). More recently, and following the wild bird‐trade ban imposed by the European Union, there has been an increase in trade between tropical countries across biogeographical realms (Reino et al., 2017). How socioeconomic dynamics will drive future tendencies in the supply and demand of traded birds is unknown (Ribeiro et al., 2019; Ribeiro, Bingre, et al., 2022), but the impact of such dynamics in the successful establishment and invasion of alien birds is, as shown here, likely to become at least partially constrained by climate change and land‐cover change.

Climate has acted as a limiting filter for widespread invasion of tropical species (Abellán et al., 2017; Strubbe & Matthysen, 2009). Cold winters in the higher latitudes with temperate and boreal climates and reduced precipitation in the lower latitudes with temperate climates, have constrained the ability of alien species of tropical origin to establish in larger numbers. But the projected gradual “tropicalization” of higher‐latitude regions, along with expanding forested areas across the region could lead to increased invasion potential, as traded birds are overwhelmingly forest dwellers. Simultaneously, a decrease in the suitability of cross‐tropical‐regions invasions is also expected. Such patterns might counter the negative impacts of the recent increase in the trade of birds across the tropics as few tropical bird species have invaded non‐native biogeographical realms, with the exception of a number of parakeet species that have established in temperate climate cities (Dyer, Redding, & Blackburn, 2017; Redding et al., 2019).

Species distribution models have well‐known limitations and uncertainties (Botkin et al., 2007; Urban et al., 2016), but the biogeographical‐level projections provided herein—indicating a tropicalization of much of higher latitudes—are consistent with general trends identified using other methodologies (e.g., Garcia et al., 2014). While analysis of individual species model results would likely expose occasional inaccuracies (e.g., Broennimann et al., 2007), the rigorous analytical pipeline developed herein, involving a large ensemble of well‐tested species distributions and climate models, is bound to limit inaccuracies to its current maximum (e.g., Araújo et al., 2005). The links of dependence of birds with land‐cover were extracted from the literature, but the resolution of species distributions and land cover do not match. Additionally, land‐cover is a limited surrogate for detailed wildlife‐habitat relationships (Barton et al., 2014). For example, when the forest class is associated with a species preference, it neglects that forests can be very different and that species might be dependent on specific types of forests, not just any forest. Inferred increased suitability of invasion owing to forest expansion in high latitudes should thus be interpreted cautiously, while reductions of suitability to invasion of forest‐dwellers across the tropics, owing to deforestation, are probably more generalizable.

More severe are uncertainties related with the examination of the joint effects of land‐cover and climate change, which are unlikely to be linear. That is, for any given unit of increase in climate change we assume, for lack of better assumption, the same unit of impact in land‐cover change. Not just the units between land‐cover and climate change are different (Thuiller, Araújo, & Lavorel, 2004), but the rate with which species respond to the changes in each one of these environmental changes is also likely to differ (Triviño et al., 2013). For example, responses to climate change are often delayed, generating lagged responses (Devictor et al., 2008; Forero‐Medina et al., 2011), whereas land‐cover transformations can have immediate impacts. At present, we do not have more sophisticated approaches to jointly examine the interactions between land‐cover and climate change impacts, but the simple analysis provided is a first approximation to the problem.

Additionally, changes in socioeconomic needs, cultural preferences and ease of trade may result in different species being selected for transport over time. For example, game‐bird trade dominated between the 14th and 18th century being gradually replaced by Passeriformes and Psittaciformes, mainly as pets (Dyer, Cassey, et al., 2017). While absolute trade volumes are still dominated by passerine birds and parrots and allies, in the 1995–2015 CITES database, Psittaciformes are still the most frequently traded followed by raptors such as Accipitriformes and Falconiformes, and owl species (Strigiformes). These changes are at least partly driven by preferred species traits related to physical attractiveness, such as coloration and song (Su et al., 2014), price, as in the case of the social status associated with falconry in the Middle East (Soorae et al., 2008), or cultural phenomena such as the popularity in Japan of ‘bird cafes’ showcasing exotic raptors and owls (Vall‐Llosera & Su, 2019) or the possible possible effect of Harry Potter movies on demand for owls (Vesper, 2017). Predicting the identity of species that will be traded in the future is fraught with uncertainties. The sensitivity of traded species to environmental change will depend very much on their identity (i.e., trait) and origin (i.e., climate niche). Our results uncover the importance of explicitly considering the joint effects of future changes in climate and land cover for delivering robust forecasts of invasion risk. However, forecasts are based on projections and assume that trends in the identity and origin of traded birds remains constant. Should these trends change, so would forecasts.

Biodiversity‐risk protocols increasingly use models to provide estimates of species probabilities of occurrence against a number of threatening factors (Williams & Araújo, 2000, 2002). More recently, the framework has been extended to deal with projected changes in the quantities of interest, such as estimated changes in species' climatic suitability (Araújo et al., 2011) or land cover (Araújo et al., 2008). For example, in the context of risks to invasion, the European Union's main Regulation (1143/2014) for managing biological invasions mandates assessing invasion risks under ‘foreseeable climate change conditions. However, invasion risk is difficult to predict, given its dependence on multiple factors acting in synergy. Nonetheless, factors beyond climate change, such as concurrent changes in land cover are typically not considered by existing studies. This can lead to overestimating risk when climate suitability is not matched by availability of suitable habitat (Hof et al., 2011), or underestimating risk when interactions of climate change and land‐cover change interact positively in the suitability for invasion. Our study delivers spatially explicit projections of estimated changes in climate and land‐cover suitability for 610 traded birds, thus providing a crucial basis to anticipate regions with varying propensity to invasion. Effective biosecurity measures will benefit from increased uptake of modeling scenarios that consider both climate and land‐cover change. Invasive species risk analyses are typically performed for a specific (set of) species identified as potential invaders by, e.g., horizon‐scanning exercises or early‐warning networks in given geographical areas (Turbé et al., 2017). Our modeling results can be directly used for a first assessment of the establishment and geographical spread risk components of any risk analysis procedure aimed at any of the traded bird species we consider here. Alternatively, bespoke finer‐grained forecasts of where invasive species are most likely to colonize can easily be obtained by leveraging our modeled suitable areas as input for dispersal‐explicit models simulating population dynamics across landscapes (Fordham et al., 2012, 2013), while using invasion histories of previously established invaders as proxy for the strength and influence of biotic interactions on invasive spread (Lovell et al., 2021).

AUTHOR CONTRIBUTIONS

Miguel B. Araújo and Luís Reino conceptualized and planned the research. B.N. conducted the analysis. Joana Ribeiro, Luís Reino, César Capinha, and Diederik Strubbe provided the habitat data for species. Carsten Rahbek provided the bird range data. Babak Naimi and Miguel B. Araújo wrote the manuscript with contributions from all authors.

CONFLICT OF INTEREST

No competing interests.

Data S1

Click here for additional data file.