Predicted alteration of vertebrate communities in response to climate‐induced elevational shifts

Climate change is driving species to migrate to novel areas as current environments become unsuitable. As a result, species distributions have shifted uphill in montane ecosystems globally. Heterogeneous dispersal rates among shifting species could result in complex changes to community assemblages. For example, interspecific differences in dispersal ability could lead to the disruption, or creation, of species interactions and processes within communities, likely amplifying the impact of climate change on ecosystems. Here, we studied the dispersal success of vertebrate species in a tropical montane ecosystem under a climate‐induced uphill shift and assessed the derived impacts on community structures.

Forecasting and mitigating the impacts of climate change presents a major challenge for ecologists and conservationists, but we have generally been unsuccessful at generating comprehensive predictions or unambiguous recommendations for mitigation (Gilman et al., 2010). Due to the intrinsic challenges of studying species' responses to climate change, predictions have systematically focused on the direct effect of a changing environment on individual species.
Studies on the response of single species to climate change are useful when studying endangered or invasive species, for which targeted conservation approaches need to be implemented (Kulhanek et al., 2011;Sousa-Silva et al., 2014). However, those studies usually do not provide insights into the impacts of climate change at the ecosystem level, where community"re-shuffling" could play a critical role in the response of ecosystems to the effects of climate change (Gibson-Reinemer et al., 2015;Gilman et al., 2010;Tylianakis et al., 2008).
In open communities, under climatic-induced geographical shift, dispersal differences among species could dissolve existing species interactions, eliminating some processes and creating new ones (Ackerly, 2003;Schweiger et al., 2008). Disrupted and novel interactions will likely alter abundances, distributions and extinction probabilities of species (Carpenter et al., 2008;Deutsch et al., 2008;Williams & de la Fuente, 2021). In this regard, metacommunity dynamics will be dependent on dispersal restrictions under climate change (Gilman et al., 2010), causing extinctions of low-dispersal species and favouring the range expansion of high-dispersal species.
The present study examined the effects of climate-induced uphill elevation shifts in a tropical montane ecosystem by simulating the shift of 7613 community assemblages derived from empirical observations (Williams et al., 2010). We hypothesised that (1) species dispersal success under climate-induced uphill shifts would be determined by the interaction among species' dispersal ability, landscape composition and global warming. We predicted that dispersal success would decrease with elevation due to an unavoidable increase in the isolation of patches as elevation increases, creating the so-called "montane islands" (Brown, 1971(Brown, , 1978. Furthermore, we expected that the declining pattern in dispersal success with elevation would be exacerbated by global warming, as climate change is expected to restrict species shifts (Gilman et al., 2010) by reducing thermo-suitable habitat for movement across the landscape (Meade et al., 2018). To model the likelihood of species shifting uphill, we adopted a spatially explicit approach that accounted for the interaction between species dispersal ability and landscape composition.
The effect of global warming on dispersal success was analysed by implementing a sequential increase in landscape resistance based on six different climate change scenarios. Given that the likelihood of species shifting depends on species dispersal ability, we also hypothesised that (2) heterogeneous dispersal rates among shifting species should result in changes in community assemblages, potentially altering species co-existence. We measured potential changes in community assemblages and species co-existence along the elevational gradient as the pairwise temporal change in β-diversity and probabilistic change in species co-occurrence. The predicted increasing tendency of temporal β-dissimilarity between assemblages and the loss of existing patterns of species co-occurrence are expected to escalate with elevation, derived from an increasing local population extinction rate towards mountaintops, that is, a classic pattern of an "escalator to extinction" (Marris, 2007;Urban, 2018).

| Study area
The Australian Wet Tropics World Heritage Area is composed of mixed tropical rainforest in an area of approximately 7000 km 2 ( Figure 1a). The elevation varies from sea level to highlands at species interactions, population dynamics and species potential to adapt to a changing environment.

K E Y W O R D S
community reshuffling, escalator to extinction, global warming, landscape ecology, range shift, species co-occurrence, species extinction, tropical mountain 1000 m, with isolated peaks reaching up to 1620 m (Nix & Switzer, 1991). The total rainforest area increases with elevation towards the midlands (600-700 m; Figure 1b), which are mainly dominated by tablelands. As elevation increases from the midlands, the total rainforest area decreases rapidly towards the isolated mountaintops, with <2% of the total rainforest area extending above 1200 m in elevation ( Figure 1b). Compared to other habitats, rainforest coverage increases with elevation, being the predominant habitat type at high altitudes ( Figure S1).

| Rainforest patches and elevational bands
The elevational gradient of the Australian Wet Tropics was segregated into 15 bands of 100 m width using a 3-arc-second (~90 m 2 ) resolution digital elevation model (DEM; Figure 1a). The reclassified elevation raster was then cropped using rainforest and adjacent wet sclerophyll (rainforest buffer areas) forest types selected from the broad vegetation groups of Queensland (Neldner et al., 2019). The resultant raster was thinned using a pixel clumping function, applying a "rook's case" direction to define what was considered adjacent. This approach was followed to reduce the noise produced by small isolated patches in the agricultural areas (e.g. crop edges).
We defined the patches used to identify community assemblages as clumps of 100 pixels or more (~0.009 km 2 ). This filtering algorithm reduced the total number of patches at low elevations but did not affect the total rainforest area ( Figure S2), optimizing the number of biologically relevant patch sizes for vertebrate communities without significantly affecting the total area used in the simulation ( Figures   S2 and S3). Using the same algorithm, we defined a second raster with a more relaxed thinning, removing clumps of less than 10 pixels (~0.0009 km 2 ). This second raster was used to define the transition layer for movement analyses. Thus, patches that might be too small to contain.

| Transition layer for movement analysis
We developed a resistance layer for every elevational band and climate change scenario combination (n = 90). The transition layer developed for movement analysis was defined using a thermal resistance surface ( Figure S4). The thermal resistance surface was created using a regionally adapted thermal raster of annual mean temperature (Storlie et al., 2013). Values in the temperature raster were reclassified into resistance units, with a resistance unit defined as a change of 0.25°C in annual mean temperature (i.e. the maximum thermal difference between adjacent rainforest patches). Thermal layers were elevation specific, and the resistance gradient was calculated using the increasing deviance from the mean temperature of the elevational bands (e.g. a pixel in the temperature raster with one-degree deviation from the mean band's temperature was transformed to a pixel of four resistance units). We assigned a fixed resistance value to the matrix (i.e. cleared land) of 50% higher than the highest thermal resistance value associated with a rainforest patch, favouring the avoidance of the matrix when possible. Transition values within the transition layer were calculated using the mean thermal resistance of adjacent pixels in the thermal resistance surface, using a "rook's case" direction to define what pixels were considered adjacent. In addition to the transition layer based on current temperature, we developed five extra transition layers for each elevational band, defined by five global warming scenarios, from +1°C to +5°C with 1°C steps.

| Community assemblages
Community assemblages were defined at the established rainforest patches (see section 2.2.1) by extracting the predicted potential habitat from species distribution models (SDMs) for 202 rainforest vertebrates in the region (Tables S1 and S2; Figure S5). These SDMs were developed and critically assessed by Williams et al. (2010) to refine the ecological realism of each species' model (Reside et al., 2019). The refinement of species potential distribution was based on expert opinion, species records, literature-based evidence and empirical model validation (i.e. targeted field expeditions with approximately 11,650 independent standardized surveys).

| Species dispersal capacity
Species dispersal capacity was categorized into six groups from low to very high dispersal based on the potential dispersal index developed by Williams et al. (2010) (Table S1; Figure S6). The index was calculated using characteristics known to influence the capacity for dispersal: flight, known movement sourced from literature and body size. To calculate the potential dispersal index for frogs, a further category was included to identify those species which have a freshwater tadpole dispersal stage (Williams et al., 2010). Low dispersal capacity was defined by a maximum travel distance of 0.5 units in the transition layer, and very high dispersal was defined by a maximum of 16 units, with each category able to disperse double than the previous category (i.e. 0.5, 1, 2, 4, 8 and 16 units).

Least-cost path simulation
Using each combination of elevational band and transition layer, we calculated the least-cost path between all the patches' centroids from an elevational band and all the patches' centroids in the next higher elevational band using Dijkstra's algorithm (Dijkstra, 1959).
This process simulated the cost of shifting uphill for species in each patch of forest in each elevational band, resulting in a large matrix of least-cost paths containing all possible paths between rainforest patches ( Figure S4).

Modelling dispersal success
We defined dispersal success as the percentage of species within an elevational band able to colonize a higher elevational band. Dispersal was considered successful when the minimum least-cost path to a higher elevational band was smaller than a threshold, defined as species dispersal capacity (see Section 2.2.4; Figure S4). For example, if the minimum least-cost path between patch A and patch B in the transition layer was four units, only species with a dispersal capacity of four or higher would be able to colonize patch B ( Figure S4). Thus, species extinction probability at a given elevational band would result from the combination of local extinction probability of populations in all the patches at that elevation.

Change in total area
In addition to the dispersal success, we also measured the total change in area for each combination of elevation, global warming scenario and dispersal capacity. Total area change was calculated as the difference between the total area colonized and the total area lost due to the elevational shift. In this regard, species with a high dispersal capacity would benefit from the elevational shift when the area colonized was larger than the source area. In contrast, low-dispersal species would reduce their total habitat when their inability to shift to higher elevations resulted in local population extinction at the source patch.

Temporal β-diversity
The variation in community composition was measured as the pairwise temporal change in β-dissimilarity between assemblages in the baseline and elevational shift scenarios along the elevational gradient. Rainforest vertebrate species were segregated for analyses into the main four taxa (i.e. birds, mammals, frogs and reptiles) to account for potential intertaxon variation associated with inherent differences in dispersal capacity ( Figure S6). The overall change in dissimilarity was measured as Jaccard similarity coefficient (β jac ).
We further explored the partition of the measured dissimilarity into its turnover (β jtu ) and nestedness (β jne ) components (Baselga, 2010(Baselga, , 2012. The partition in the temporal dissimilarity allowed us to distinguish whether the change in dissimilarity was driven by substitution of species (β jtu ) or loss/gain of species in a nested pattern (β jne ). We used the ratio between β jne and β jac to describe the relative contribution of each component in the change in β-dissimilarity (Dobrovolski et al., 2012). Ratio values above 0.5 indicated that the change in dissimilarity was driven mainly by nestedness, while a ratio below 0.5 suggested that the pattern in dissimilarity was mainly driven by spatial turnover. Additionally, we followed Baselga et al. (2015)'s approach to determine whether the observed pairwise temporal βdissimilarity and partition values were significantly different from random expectations using a FE (Fixed patches total; Equiprobable species total) null model (Ulrich & Gotelli, 2007). This null model was simulated by reshuffling the species presence in patches while maintaining the species' frequencies constant, creating vertebrate metacommunities where species randomly go extinct or colonize higher elevations from a common regional species pool. This null model was used to randomize each elevational shift scenario at each elevational band, which were then used to recalculate the pairwise temporal change in β-dissimilarity. Null models were then compared to the observed change in β-dissimilarity using Kolmogorov-Smirnov tests.
This process was iterated 1000 times for each combination of taxa, elevational band and climate scenario using 1000 null models to assess the robustness of our simulation.

Spatial β-diversity
Spatial β-diversity studies have shown that pairwise β-diversity and multiple-site β-diversity metrics can reveal different patterns (Baselga, 2013;Baselga et al., 2015). Thus, in addition to pairwise temporal β-dissimilarity analyses, we calculated multiple-site spatial dissimilarities for each scenario and baseline by computing 1000 multiple-site dissimilarity values using 25 randomly sampled patches at each elevational band. We used the dissimilarity distributions derived from this analysis and compared the resultant β-dissimilarity values with their temporal pairwise counterparts.
Given that the difference in results using the temporal and spatial approaches showed almost an identical simulation outcome, we have only presented the pairwise temporal dissimilarity findings in the results section, as our primary focus was to study changes at the population level (i.e. pairwise comparison among assemblages).
Spatial β-diversity results are included in Figures S7 and S8.

Species co-occurrence
We measured the impact of the simulated elevational shifts on species co-occurrence using a probabilistic model, comparing expected versus observed species co-occurrence at each combination of taxa, elevational band and climate change scenario (Veech, 2013).
A positive association between two species was determined when the observed frequency of co-occurrence was significantly large and greater than expected, while a negative association was defined by a frequency of co-occurrence significantly small and less than expected (Veech, 2013).

Sensitivity analyses
We tested the sensitivity of our simulation to different parametrization of species dispersal capacity. We developed six scenarios where we altered species-specific dispersal capacity index by ±5, 10 and 50% to explore whether our inferences about assemblage changes were markedly driven by the initial species' dispersal ability rather than by the interspecific differences in dispersal capacity.

| Hypothesis 1: Dispersal success under climate-induced uphill shifts
The simulated climate-induced uphill shift in the Australian Wet Tropics suggested a high rate of local population extinction (i.e. the loss of a species in all the assemblages within an elevational band) and a strong reduction in the total area occupied by species Total change in suitable area segregated by species dispersal capacity. Colours depict the sequential increase in global warming from the current climate (+0°C; dark blue) to +5°C scenario (yellow), with a +1°C step. The horizontal black line at 0% represents the inflection point between "winner" (increase in total area) and "loser" species. Subpanels show species dispersal capacity increased with dispersal capacity in all the scenarios, with some species benefitting from the elevational shift when the source area was smaller than the colonised novel habitat (total area change >0%).

| Change in pairwise temporal β-diversity
The effects of a heterogeneous dispersal success across species, elevation and global warming scenarios were observed in the community assemblages of the region (Figure 3; Figure 4). The pairwise temporal variation in assemblages between baseline and elevational shift scenarios resulted in a substantial increase in β-dissimilarity along the elevational gradient ( Figure 3). As expected, the temporal variation in β-dissimilarity was dominated by the nestedness component (β jne /β jac , mean ± SD = 0.86 ± 0.02; Figure S9), implying a strong pattern of species loss between the baseline and the elevational shift scenarios (Figure 4). This pattern was particularly marked on mountaintops, with a reduced number of species able to reach areas at the mountain peaks across the region, indicated by a mass local extinction of species attempting to shift from 1200-to 1300-m and from 1300-to 1400-m elevational bands (Figure 4). The pattern in community changes showed high consistency across taxa (Figures 3 and 4).
However, the increasing change in assemblages and local extinction rate tended slower for species with high dispersal capacity (i.e. birds and mammals; Figure S6) than for low-dispersal species (i.e. frogs and reptiles; Figure S6), highlighting the expected influence of the interaction between landscape resistance and species' dispersal capacity.
The null model approach adopted to validate the robustness of our results showed high consistency for all taxa. The observed overall dissimilarity (β jac ) and nestedness component (β jne ) were significantly different from random expectations in all 1000 trials across scenarios, elevation and taxa (Kolmogorov-Smirnov D > 0.15, p < .03; Table S3). However, results did not significantly differ from null models at mountaintops for reptiles, mammals and frogs (D < 0.26, p > .32 for β jac ; D < 0.50, p > .11 for β jne ; Table S3). The lack of significance in the overall change in β-dissimilarity and the nestedness component at 1300 m compared to the null model for these taxa could be an artefact of the combination of the low number and high isolation of patches at mountaintops. The high level of isolation combined with a low number of rainforest patches at mountaintops resulted in most assemblages (89%) being unable to shift (minimum least-cost path >highest dispersal capacity). Thus, the variation in temporal β-dissimilarity at the 1300-m elevational band for these taxa was derived from a very low number of assemblages (n = 3), increasing the chances of getting a similar result when assemblage composition was assigned randomly.

| Change in species co-occurrence
The increasing rate of species local extinction across elevation was reflected in the predicted changes in species co-existence, showing high consistency across taxa. The trend in positive co-occurrence change was particularly accentuated at high elevations (>900 m), suggesting that the lower local extinction rate in the lowlands and midlands had a marginal influence on positive co-existence ( Figure 5). However, the pattern in negative co-existence (i.e. the F I G U R E 3 Pairwise change in temporal β-dissimilarity across taxa in the Australian Wet Tropics. The figure depicts the overall change in β-dissimilarity (β jac ) across the elevational gradient. Colours depict the sequential increase in global warming from the current climate (+0°C; dark blue) to +5°C scenario (yellow), with a +1°C step F I G U R E 4 Local extinction rate derived from an uphill range shift in the Australian Wet Tropics. Values represent the percentage of species loss across taxa, elevational gradient and climate change scenarios. Colours depict the sequential increase in global warming from the current climate (+0°C; dark blue) to +5°C scenario (yellow), with a +1°C step number of pair of species with a high probability of not occurring in the same community) showed higher sensitivity to low rates of local extinction than the positive co-occurrence trend ( Figures 5   and 6), highlighting the potential disruption of species negative coexistence even at lower elevations.

| Sensitivity analyses
Results from sensitivity analyses reinforced the robustness of our findings. Changes in dispersal capacity parametrization did not substantially influence the outcome of the simulation (Figures S10 and S11).

| The impact of climate-induced uphill shift on communities and ecosystems
Globally, empirical studies suggest that patterns in climate-induced uphill shifts could be general across taxa (Chen et al., 2011).
However, the ecological implications at the ecosystem level resulting from community distributional shifts have not been broadly considered. We found that an uphill shift in the Australian Wet Tropics would significantly change community assemblages and species co-existence. The observed changes in community structures evidenced the effect of heterogeneous dispersal potential among shifting species. The influence of dispersal restriction across elevation resulted in a clear pattern of increasing temporal β-dissimilarity due to a marked local extinction rate, leading to a decline in species co-occurrence across elevations. The severity of the impact of species shifting in elevation in tropical montane forests was especially noticeable in the highlands, with many species being unable to shift to the isolated mountaintops, which could lead to mass local extinctions. The ecosystem resilience to environmental changes could be further compromised due to the potential loss of species co-existence and ecosystem processes, denoted by a marked decline in species co-occurrence along the elevational gradient. However, some species could benefit from climate-induced elevational shifts, as shown by high-dispersal species shifting to larger areas.

| Escalator to extinction: response of montane communities to climate change across the elevational gradient
Our findings align with those predicted from uphill shifts, showing a classic example of the "escalator to extinction" (Marris, 2007;Urban, 2018). As evidenced in this study, the impact of climate change escalated exponentially towards high elevations, leading to the local extinction of 55% on average of the species pushed towards mountaintops. The increasing local extinction rate with elevation suggested a marked homogenization of assemblages, denoted by a change in β-dissimilarity in a nested pattern close to 1 towards mountain peaks, and a marked reduction in species co-existence, defined by an exponential decrease in species positive co-occurrence towards mountaintops. These findings might reflect the potential differences in community stability across elevations, with upland F I G U R E 5 Change in positive species co-occurrence compared to the current number of probable co-existing species (baseline) across the elevational gradient. Values below the baseline threshold (horizontal dashed line) represent a decrease in co-occurrence. Colours depict the sequential increase in global warming from the current climate (+0°C; dark blue) to +5°C scenario (yellow), with a +1°C step F I G U R E 6 Change in negative species co-occurrence (the number of pair of species with a high probability of not occurring in the same community) compared to the current number of probable co-existing species (baseline) across the elevational gradient. Values below the baseline threshold (horizontal dashed line) represent a decrease in co-occurrence. Colours depict the sequential increase in global warming from the current climate (+0°C; dark blue) to +5°C scenario (yellow), with a +1°C step assemblages showing lower plasticity to a changing environment than their counterparts at lower elevations (Freeman, Scholer, et al., 2018;Williams & de la Fuente, 2021). Our results are in accordance with projected impacts of climate change in mountain ecosystems (Mamantov et al., 2021;Manes et al., 2021;Williams et al., 2003), warning of severe deterioration of montane communities, which host unique ecosystem services (Korner & Spehn, 2019), and a disproportionate amount of global biodiversity (Jetz et al., 2004;Quintero & Jetz, 2018). For example, many upland species in the Australian Wet Tropics are highly specialized to high elevations and rarely occur in the lowlands (Nix & Switzer, 1991;Williams et al., 1995Williams et al., , 2010. The marked local extinction rate at high elevations predicted across the Australian Wet Tropics could mean the extinction of upland specialists at the regional scale, even with the most conservative elevational shift scenario. Furthermore, the significant change in community assemblages derived from a strong pattern of local extinctions could lead to significant changes in ecological processes and interactions (Gilman et al., 2010;Tylianakis et al., 2008), potentially altering key biotic interactions and further diminishing ecosystem resilience.

| Simulation assumption: a forced uphill shift
Species in montane ecosystems can show a variety of responses to climate change (Mccain & King, 2014;Tingley et al., 2012). The potential for species to adapt or even move downhill in elevation has been previously reported based on demographic responses and changes in occupancy (Chen et al., 2011;Moritz et al., 2008;Neate-Clegg et al., 2021;Tingley et al., 2012;Williams & de la Fuente, 2021). However, the underlying mechanisms for species heterogeneous range shifts and demographic fluctuations could also be related to lags in response to changes in the environmental conditions (Alexander et al., 2018;Bertrand et al., 2011;Fordham et al., 2016;La Sorte & Jetz, 2012;Lu et al., 2021), rather than species adaptability to climate change. Additionally, even with the potential for rapid adaptation reported for some taxa ( extinct. This implied that the overall change in β-dissimilarity was unavoidably driven by nestedness (i.e. β jne /β jac > 0.5; Figure S9), as the simulation entailed a net loss of species due to local extinction or no change at all (i.e. β jac = 0).

| Limitations and perspectives: towards a mechanistic approach
Our results provided compelling evidence for significant species loss due to the direct impacts of climate change, affecting ecosystems at a regional scale. The findings derived from this study highlight the importance of considering community-level elevational shift processes when studying the impact of climate change. However, some limitations in our simulation are worth noting. This study would have benefited from information on species-specific adaptation capacity, species interactions and explicit demographic dynamics (Williams et al., 2008). By incorporating potential species adaptation to climate change, we could have accounted for idiosyncrasies in species shifts (Tingley et al., 2012), allowing for a deeper understanding of species-specific impacts of climate-induced distributional shifts.
Additionally, a broader knowledge of species interactions in rainforest ecosystems would have enabled the explicit estimation of the impacts of community "re-shuffling" on disruptions of key ecosystem processes and the creation of novel interactions. Moreover, a deeper understanding of explicit migration, death, and reproductive success is needed to accurately forecast the likelihood of survival in the colonized patches. In this regard, an explicit simulation of individuals' movement within patches could account for the potential under or overestimation of migration rate induced by patch size (e.g. differences between small and big species moving through a large patch), potentially providing a more robust measure of shift rates. Combined, explicit demography dynamics, adaptation capacity and species interactions would provide a more mechanistic approach required to inform conservation decisions at the population level. Finally, we note the potential influence of synchronic habitat shifts that could affect the dispersal likelihood and establishment of migrating vertebrates. However, synchronous vegetation movement is likely to have a negligible influence on tropical vertebrate migrations, as vegetation tend to show longer delays in response to environmental changes (Alexander et al., 2018;Lu et al., 2021) compared to the fast shift shown by vertebrates (Freeman & Freeman, 2014;Williams & de la Fuente, 2021). Thus, even though evidence worldwide (Chen et al., 2011), and specifically in the Australian Wet Tropics (Williams & de la Fuente, 2021), shows that climate-induced distributional shifts in montane ecosystems tend uphill, our results could have been slightly different if the incorporation of species resilience, population dynamics, adaptation potential, and synchronization in habitat shift were possible. Nonetheless, we are confident that the overall pattern in assemblage changes would have followed a similar trajectory, highlighting the marked vulnerability of tropical highland communities.

| General conclusion
This study found marked and consistent impacts of climate-induced uphill shifts on community assemblages in the region, indicated by an escalating increase in temporal β-dissimilarity and the associated escalating decline in species co-occurrence. We also discovered that an uphill elevational shift in montane rainforests of the Australian Wet Tropics could severely impact upland ecosystems.
The predicted potential mass extinction of species in the highlands of the Wet Tropics would result in the collapse of upland montane ecosystems, threatening many endemic species solely restricted to high elevations. In conclusion, our results highlight the importance of considering the community "reshuffling" process resulting from multispecies synchronic distributional shifts. Furthermore, understanding the accumulation of individual changes in community assemblages is crucial to provide sensitive predictions at the ecosystem level, providing a more informative research input to conservation planning and management.

ACK N OWLED G EM ENTS
The authors thank James Cook University, Skyrail Foundation, Abriculture, Queensland Department of Environment and Science, Wet Tropics Management Authority, Earthwatch Institute, Rainforest-CRC, Marine and Tropical Research Facility and the National Environmental Research Program for funding support. The authors also thank the many volunteers, students and assistants who have contributed to the fieldwork over the years.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available in