Alpine plants are on the move: Quantifying distribution shifts of Australian alpine plants through time

Alpine plant species’ distributions are thought to have been shifting to higher elevations in response to climate change. By moving upslope, species can occupy cooler and more suitable environments as climate change warms their current ranges. Despite evidence of upslope migration in the northern hemisphere, there is limited evidence for elevational shifts in southern hemisphere plants. Our study aimed to determine if alpine plants in Australia have migrated upslope in the last 2 to 6 decades.


| INTRODUC TI ON
Many species are predicted to undergo shifts in distribution in response to changes in climate (Parmesan & Yohe, 2003;Peterson et al., 2011;Thuiller et al., 2005). In line with these predictions, several northern hemisphere empirical studies have found evidence for both plant and animal species shifting their ranges poleward or upslope (Grabherr et al., 1994;Jump et al., 2012;Koide et al., 2017;Konvicka et al., 2003;Liang et al., 2018;Moritz et al., 2008;Parolo & Rossi, 2008). However, much less is known about the response of species in the southern hemisphere. We should not assume that southern species are showing similar responses to those in the north, as the southern hemisphere has a lower rate of climate change (Friedman et al., 2013) and very different species composition from that of the northern hemisphere (Box, 2002). Our study will address this knowledge gap by quantifying how Australian alpine plants have shifted their distributions over the last 20-60 years.
Alpine areas are habitats for many endangered and endemic species that are highly susceptible to climatic changes (Dirnböck et al., 2011). Along with their biological significance, alpine environments are socially, culturally and economically important (Costin et al., 2000;Hughes, 2011). Within these highly important ecoregions, rapid warming has been found to be amplified at higher elevations (Beniston et al., 1997;McGowan et al., 2018;Pepin et al., 2015).
There are already noticeable impacts stemming from this increased warming including habitat loss and population decline in plants and animals (Grabherr et al., 2010;Halloy & Mark, 2003;Rehnus et al., 2018). It has been predicted that distribution shifts in alpine flora and fauna will occur in the Australian alpine region, including poleward range shifts and upward range shifts in elevation (Pickering et al., 2004). Several studies have quantified the impacts of climate change in Australian alpine areas including increases in seedling establishment of Eucalyptus pauciflora subsp. niphophila at lower elevations since 1970, which may lead to downhill range shifts in this species (Green & Venn, 2012). Studies have also found changes in species composition including increased shrub cover and decreased graminoid cover with experimental warming in this region (Wahren et al., 2013); as well as observed shifts in species' phenology with a few species flowering earlier in the year than in the past (Gallagher et al., 2009) and increases in the invasion of non-native species (Scherrer & Pickering, 2001). However, our study is the first to test for elevational range shifts across a wide range of Australian alpine plant species.
Since average temperatures decrease by 6.5°C per 1000 m of elevation gained (Wallace & Hobbs, 2006), it is predicted that most alpine species will move upslope in response to climate change to remain within suitable climatic conditions (Parmesan & Yohe, 2003;Pickering et al., 2004). However, northern hemisphere studies have revealed substantial variation between species in both the direction and the speed of range shifts (Lenoir & Svenning, 2015;Rumpf et al., 2018). Our first goal was to test the hypothesis that species' mean elevations are shifting uphill and to quantify the proportion of Australian alpine plants undergoing uphill range shifts, the proportion of species whose ranges are staying static and the proportion whose ranges are shifting downhill.
In addition to quantifying species' mean elevational shifts, we quantified changes in the leading and trailing edges of species ranges. As the climate warms, it is expected that more favourable conditions will arise at the upper edge of species' ranges (Breshears et al., 2008). These conditions might allow for increased plant establishment opportunities at higher elevations, resulting in an upslope shift of species' upper range margin. Likewise, the trailing edge of each species will encounter less favourable temperatures, potentially resulting in reduced recruitment and/or increased mortality of established plants (Breshears et al., 2008). Thus, our second hypothesis was that species' upper and lower range edges have shifted upslope over time. However, it is not clear whether changes in range edges will be more or less dramatic than changes in the centre of species' ranges. Some studies suggest that plants and animals will be most sensitive to climatic changes at the fringes of their range (Lesica & McCune, 2004) and that species are typically more abundant in the centre of their range, declining gradually towards range margins (Lenoir & Svenning, 2015). However, other studies have shown that species are more abundant in suboptimal conditions (Dallas et al., 2017) or have non-linear relationships between their abundance and environmental suitability (Baer & Maron, 2020).
Understanding how climate change affects the extent of a species' elevational range (i.e. the difference between the maximum and minimum elevation at which the species occurs) is also important because species with contracting elevational ranges will be at a higher risk of becoming threatened or possibly going extinct (La Sorte & Jetz, 2010). In line with projections from the northern hemisphere (Koide et al., 2017;La Sorte & Jetz, 2010), our third hypothesis was that alpine species' elevational ranges would contract through time.
Finally, we asked whether the rate of species' distribution changes was related to growth form. Evidence from Europe shows that species with shorter life cycles (e.g. herbs and grasses) have undergone larger changes in their distribution in response to climate change over the last few decades than have species with slower life cycles (trees and shrubs; Lenoir et al., 2008). Our fourth and final hypothesis was that herbaceous monocotyledons have shifted their distributions more over time than herbaceous dicotyledons, which in turn have shifted more than woody dicotyledons.

| Location
Our study focused on alpine communities in Kosciuszko National Park, in southern New South Wales, Australia (Figure 1a-c). The tree line in this region is between 1800 and 1900 m above sea level (Pickering & Venn, 2013). Over the past 60 years, the mean temperature in this region has increased by 0.02°C per year (Hennessy et al., 2003;Sritharan et al., 2021; Figure S1a) and snow regimes are F I G U R E 1 Site location of our study within (a) Australia in (b) Kosciuszko National Park, which is also used as the boundary for selecting historic species records from the Atlas of Living Australia. (c) Shows the transects taken during fieldwork within Kosciuszko to determine modern distribution of alpine species. (d) Flowchart of the methods used to clean species occurrence data collected from the Atlas of Living Australia. Records were first scanned and removed based on three criteria, followed by the scanning and removal of whole species from the data based on four criteria. Data at the end of screening contained 36 species and 1860 records which were then used in further analyses and hypothesis testing in the study. already showing significant changes, with snow depth ( Figure S1b) and snow cover duration both declining (Bhend et al., 2012).

| Data collection
To obtain historic species distribution data, we downloaded plant occurrence records (including herbarium data and field observations) containing geographic data from the Spatial Portal on the with occurrence records in the region (Figure 1d). We then excluded any species occurrence records that did not have accurate taxonomic data (i.e. individuals had to be identified to species level or be historically differentiated if there were subspecies, Figure 1d).
Records were also excluded if they had an inaccurate or absent date and/or location data (i.e. >50 m of "Coordinate Uncertainty" for ALA records), or if there were ≥30 years between records (to avoid isolated points having undue influence on analyses). These exclusions resulted in the total number of species with records decreasing from 2393 to 1442 ( Figure 1d). We then restricted our list further to 389 species as we excluded a large number of species that did not have at least 25 records of occurrence (to ensure there were sufficient replicates for analysis) or if their data did not span until at least the year 2000 (to allow enough time to quantify a response).
Species were excluded if they had a minimum elevation recorded below 1350 m (i.e. non-alpine). We also excluded species that were considered difficult to identify taxonomically, to ensure that all records were of the correct species. After applying these final data quality filters our study included 36 species, spanning 15 families ( Figure 1d and see Supporting Information, Table S1, for full species list). These species are a subset of representative species that are common in the area (pers. obs. F. Hemmings; J. Auld), which is also clear in the fact that they were observed commonly in historic records. These 36 species also cover the most species-rich families of the local flora above the tree line (including Asteraceae, Ericaceae, Ranunculaceae, Poaceae, Orchidaceae, Cyperaceae and Juncaceae) as well as some less common families (including Lycopodaceae and Droseraceae) (Doherty et al., 2015).
We collected modern distribution data during the 2019-2020 spring-summer season to determine the current range of a subset of the species obtained from historic data (see Supporting Information, Table S1, for species that were sampled in the modern field survey). Modern data were collected along 11 transects following walking trails and roads within Kosciuszko National Park

| Statistical analysis
All analyses were performed in R Studio, version 3.6.1 (R Core Team, To determine if mean elevation was shifting upslope across all species in our study, we extracted the estimated model coefficients of the slopes from the individual species' regressions. To control for differences in species' sample sizes and standard errors, we used these in a univariate meta-analysis weighted by the sampling variance (v i ) for each species (i). This analysis was performed using the function rma.uni in the metafor package (Viechtbauer, 2010). Species that had a maximum elevation at or very near the peak of the highest mountain in Kosciuszko (2228 m) could have had a confounding effect on the results as they are unable to shift higher in elevation.
To assess the magnitude of this problem, we performed a metaanalysis with a binary moderator (predictor) variable for whether the species occurred above or below 2200 m. Since there was no significant effect for the location of the species above or below 2200 m (R 2 = 0.06, p = 0.08), we were able to present the mean elevational trend for all species in the results and include all species in the meta-analysis.

| Hypothesis 2: that species' upper and lower range edges are shifting upslope over time
Only 18 of the total 36 species with historic maximum elevations below 2200 m were included in this analysis (because species at the top of the available space cannot move any further upslope). The remaining 18 species with historic maximum elevations above 2200 m were included only in the analyses of whether lower limits were shifting upward over time as we were unable to measure their upper range edge shifts for these species as there is no higher land available for them to colonize. We performed quantile regressions using the function rq in the quantreg package (Koenker, 2020)

| Hypothesis 3: that species ranges are contracting over time
We used only species with historic occurrence below 2200 m (18 of 36 species). We used the anova function in base R (R Core Team, 2019) to determine whether the upper (95th) and lower (5th) quantiles were significantly different, which indicates range contraction or expansion. We conducted a p-adjustment due to multiple comparison testing, using the Holm-Bonferroni method in the p.adjust function in base R (R Core Team, 2019).

| Hypothesis 4: that herbaceous monocotyledons have shifted their distributions more over time than herbaceous dicotyledons, which in turn have shifted more than woody dicotyledons
We categorized each species into one of three growth form types (see Table S1 for classifications). For our species, there were 14 herbaceous monocotyledons, 16 herbaceous dicotyledons and 5 woody dicotyledons. One species (Lycopodium fastigiatum) is a club moss that differs greatly from the other species in growth form, so we excluded this outlier from our analysis of growth form. We performed a meta-analysis on the mean linear regression of all species with a moderator (predictor) variable for growth form.

| RE SULTS
Almost a third (11 of 36) of the species in our study have undergone a significant upslope shift through time (p < .05, Figures 2 and 3 year downward in mean elevation ( Figure 3). However, none of these downslope shifts remained significant after adjusting for multiple hypothesis testing. The remaining 21 species (58%) showed no significant shift over time (p > .05). When combining all species mean elevational shifts in a meta-analysis, we found an average range shift of approximately 1.12 m increase in elevation per year, but this upslope trend was not significantly different from zero (Z = 1.23, p = .22; Figure S2).
Ten of thirty-six species showed a significant upslope shift in the lower range edge (5th quantile) of their distribution through time (p < .05, Figures 2 and 3 Figure 3), but this trend did not remain significant after p-value adjustment. The total average for all significant lower range edge shifts was an upslope shift of 4.66 m per year (Table S2). The remaining 25 species showed no significant movement over time (p > .05). Counter to our prediction, we found no significant overall shift in lower range edge through time (meta-analysis, Z = 0.25, p = .80; Figure S3).
We were only able to analyse the change in the upper range edges of the 18 species whose distributions did not historically reach the maximum elevation of the available habitat. Four of these eighteen species showed a significant upslope shift in their upper range edge (95th quantile, p < .05, Table S1), with an average of 5.96 m per year (Table S2). However, only one of these results  (Table S2), but the trend did not remain significant after p-adjustment. The overall average for significant shifts in the lower range edge was 3.85 m per year (Table S2). The remaining 13 species showed no significant shift in their upper range edge (p > .05, Figure 2). The  Figure S4).
Two of the eighteen species for which we could analyse both the upper and lower range edges showed a significant difference in the slope between their leading (95th quantile) and trailing (5th quantile) edges (p < .05, Table S1), indicating a change in elevational range over time.
Orites lancifolius showed an increase in elevational range over time and this increase remained significant after P-adjustment (Figure 2). The range of Brachyscome stolonifera had contracted over time, but this change was not significant after accounting for multiple hypothesis testing.
Four of the eighteen species whose ranges were already at the top of the available topography showed significant range contractions (i.e. their lower quantile had a significant positive slope, p < .05, Figure 3). The results for Luzula acutifolia and Stackhousia pulvinaris remained significant after adjusting for multiple hypothesis testing. Unexpectedly, one species (Colobanthus affinis) showed a significant range expansion (lower quantile negatively significant, p < .05, Figure 3), but this result did not remain significant after P-adjustment. The remaining 13 of 18 species' ranges have neither contracted nor expanded over time (p > .05).
Although we excluded species that had ≥30 years between data points to mitigate the problem of isolated points having undue influence on the analyses, we found that some species still showed large breaks in occurrence data through time, particularly before the year 2000 (Figures 2 and 3). In a post hoc analysis, we selected species that had three or fewer points of data before 1985 and no data between 1985 and 2000 (i.e. a large period of time between historic occurrence and modern occurrence data). We then removed these few data before 1985 and reanalysed each elevation regression (mean, 95th quantile and 5th quantile) and found that for four species (Brachyscome stolonifera, Carex hypandra, Luzula acutifolia subsp. nana and Pappochroma setosum) the shifts in mean elevation were no longer significant after removing these data points (Table S3) and for one species, Plantago glacialis, the shift in the upper quantile of elevation became significant after removing the outlying data points (Table S3) have shifted their ranges uphill over the last 2 to 6 decades ( Figure 5). That is, some southern hemisphere species seem to be mitigating the effects of climate change by shifting to cooler habitats in a similar manner as their northern hemisphere counterparts (Parolo & Rossi, 2008) and species in the tropics (Jiménez-García et al., 2021;Morueta-Holme et al., 2015). Although this is a welcome news, upslope migration is not a long-term solution for Australian species. Australia has much lower mountain peaks than Europe or North America, so species have a limited amount of higher ground to colonize. For example, if current trends continue, the elevation of the lower range edge of Astelia psychrocharis will reach the peak elevation in NSW in under 50 years. This species is already classified as endangered in Victoria (VicFlora, 2016), where the peaks are all below 2000 m.
Four of our study species showed a significant downslope shift through time (Figures 2, 3 and 5). One potential explanation for downhill shifts is changing biotic pressures. First, warming may allow more competitive species from lower elevations to colonize higher elevations. Increased competitive pressure at high elevations could reduce the realized niches of some higher-elevation species, resulting in a range contraction and a corresponding downhill shift of the centre of the range ( Figure S5). However, all four of the species that showed significant downslope shifts in our study historically occurred at the mountain peak, so we were not able to quantify changes in their upper range limits and thus cannot test this idea. Second, as the dominant species at lower elevations become physiologically less well suited to the changed climatic conditions, their competitive ability may decrease, potentially allowing higher-elevation species that were previously competitively excluded from a region to establish (although the high elevation species may have to compete for these habitats with species moving in from lower elevations). To test these ideas, future empirical analyses could look at changes in community composition and structure, rather than focussing solely on species' ranges. Another possible explanation for downslope shifts is that cold air drainage, cold sinks or temperature inversions may occur at intervals along elevational gradients and provide species with favourable habitat at lower elevations (Wallace & Hobbs, 2006). Other topographic factors such as the slope and aspect of the mountain may also have an influence on species' elevation range shifts and future studies could repeat field sampling in the transects where we collected modern occurrence data. This would enable researchers to gain long-term data of individual species' distributions on particular slopes/aspects of the mountain region in order to determine whether topography plays a role in elevational shifts. Finally, along with the obvious increase in temperature through time due to global climate change, there have been changes in other aspects of climate such as precipitation, snow cover duration and water availabilitychanges which could also be driving these downslope range shifts (Lenoir et al., 2010). Improving understanding of which species are moving downhill and the mechanisms underpinning these responses is clearly a worthwhile direction for future work.
Many species showed no shift in elevation over time ( Figure 5).
This may result from barriers to migration such as competition, resource availability and changing snowmelt regimes (Cannone et al., 2007;Pickering et al., 2004). Species in some ecosystems may cope F I G U R E 3 Individual species' relationships between year and elevation for species that had a historic maximum elevation above 2200 m. Solid blue lines represent significant linear regressions on the mean (p < .05). Green lines represent lower (5th) quantile regressions, with dashed lines showing non-significant quantile regressions (p > .05) and solid lines showing significant quantile regressions (p < .05). As only species with historic maximum elevations below 2200 m were included in this analysis (because species at the top of the available space cannot move any further upslope), there are no upper (95th) quantile regressions in this analysis or figure. The following symbols indicate outcomes that remained significant after a Bonferroni-Holm P-value adjustment: linear regressions (*) and lower quantile regressions (+) with increasing temperatures through plasticity or rapid evolution (Bonamour et al., 2019;Merilä, 2012). However, a recent study found evidence for morphological change through time in only 2 of 21 Australian native alpine species (Sritharan et al., 2021), suggesting that rapid evolution is not a common response to climate change in our study area. Microtopography and shifts in short horizontal distances could also be an alternative to elevational shifts as they can result in temperature differences similar to that of large elevational gradients (Scherrer & Körner, 2009).
We speculated that some species might persist by shifting poleward (as seen in other systems, particularly in the northern hemisphere; Parmesan, 2006;Virkkala & Lehikoinen, 2014). There are small areas of alpine habitat in Victoria and Tasmania, to the south of our study site (Pickering et al., 2004 could re-sample our data with the historic data found and this would provide greater evidence for species distribution shifts due to climate change in this important ecoregion. We suspect that the coming decades will reveal even more Australian species on the move. Thirteen of the species in our study have also been monitored for long-term trends in distribution shifts (albeit over a shorter time frame than our current study using historic records and on five points of elevation rather than a continuous gradient) on Mount Clarke (

| Hypothesis 3: that species ranges are contracting over time
Our study identified six species that have undergone range contractions. This is important, as one of the strongest predictors for the threat of extinction in a species is the size of its geographical range (Staude et al., 2020). Species that are restricted to a smaller elevational range may also be less resilient to environmental pressures such as invasive species and human activity such as recreational skiing, hiking and mountain biking (Costin et al., 2000). Species with smaller ranges also tend to have a lower local abundance (Gaston et al., 2000), which leads to an increased sensitivity to environmental changes. A potential direction for future research would be pairing local abundance of species within their elevational ranges and their elevational range shifts to identify possible areas of species loss.
Despite the benefits of migration away from unfavourable climate conditions, it is vital to consider the many implications on the species and their pre-existing alpine communities. Shifts in elevations of species may result in new interactions between species that were historically spatially separate (Alexander et al., 2018;Descombes et al., 2020) and the new dynamics may have effects on existing communities (Dormann & Woodin, 2002). These novel interactions have unknown impacts on species competition and survival (Walther, 2010 (Lenoir et al., 2008). However, we found no significant difference in elevation shift among growth forms ( Figure 5). This result is consistent with modelling by Lustenhouwer et al. (2017), in which the longer time until the occurrence of reproduction of larger species was counterbalanced by their greater dispersal distances and fecundity, resulting in no difference in migration rate between annual herbs, perennial herbs and shrubs.
A meta-analysis of northern hemisphere studies revealed that species had an average rate of elevation movement of 1.11 m year −1 upslope (Chen et al., 2011), which is very similar to the average rate of change in our study (1.95 m year −1 , Table S2). However, there were substantial differences between the proportion of species in our study that had shifted in elevation over time and the proportion of species showing elevational shifts in Europe and North America. In the northern hemisphere, the highest proportion of species are moving upward in elevational ranges; with 60-65% of species shifting upslope, 20-25% species shifting downslope and 10-15% of species not shifting in any direction (Chen et al., 2011;Gibson-Reinemer & Rahel, 2015;Lenoir et al., 2010), whereas our study revealed that in the Australian alpine region ~30% of species are shifting upslope, ~10% are shifting downslope and ~60% are not shifting at all. The lower proportion of species changing in our study may be due to the shorter time frame used (median ~33 years in our study cf. ~50-65 years in northern studies; Gibson-Reinemer & Rahel, 2015;Parmesan & Yohe, 2003). Another possibility is that northern hemisphere species are more likely to be responding to climate change because of the faster warming they have experienced (Friedman et al., 2013).

| CON CLUS IONS
Alpine ecosystems are considered to be particularly at risk under future climate change due to their high climate sensitivity and large numbers of endemic and threatened species (Dirnböck et al., 2011).
Identifying which species are shifting their distributions can help us to predict which species would benefit most from higher intervention conservation methods such as assisted colonization to help species migrate across isolated and fractured landscapes (Parmesan & Hanley, 2015). Our findings also suggest that future climatic warming could result in more species elevational shifts, which in turn could lead to large changes in ecological community dynamics, for instance, by changing the suite of interacting species present in different regions. We still have much to learn about the likely impact of this world changing process on our precious alpine plant communities.

CO N FLI C T O F I NTE R E S T
All authors have no conflicts of interest to declare.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13494.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data are publicly accessible through Dryad (https://doi. org/10.5061/dryad.r4xgx d2d4). Analysis code is freely available in the supplementary information and on GitHub (https://github.com/ SEver ingha m/Alpin e-Plant -Migra tion).