Climate velocities and lagged species elevational 3 shifts in mountain ranges 4

Mountain ranges support concentrations of climate-endangered endemic species, and are potential refugia for species retreating from the lowlands under 25 anthropogenic climate change. Predicting the outcome for biodiversity requires 26 knowledge of whether species are shifting uphill at the same rate as temperature 27 isotherms (i.e. whether they are successfully tracking the velocity of climatic 28 changes) 1 . Here, we provide a global assessment of the velocity of climate change 29 in mountain ranges: applying thermal dynamic theory, deriving moist adiabatic 30 lapse rates (MALR) using local surface temperature and water vapor. MALR 31 varied substantially around the world, from 3 to 9°C cooling per km elevation 32 increase. Consider the rate of terrestrial surface warming from 1971 to 2015, 24 33 regions can be identified as exhibiting high velocities where the isotherms have 34 shifted more than one standard deviation of the global mean value (> 8.45 m yr - 35 1 ). High velocities are typically found in relatively dry parts of the world, but also 36 occur in wet regions with low lapse rates, such as in Northern Sumatra, Western 37 Guiana Shield, Northern Andes, Costa Rica, Nepal, and Madagascar. Analysis of 38 biodiversity data in relation to mountain-specific velocities revealed more cases 39 of tracking between species and isotherms than previously suggested 2 and more likely occurred at lower climate velocity. Nevertheless, upslope migrations of 41 montane species have generally been lagging behind climate velocity. Such lags 42 could continue to effect change even if the climate were to stabilize immediately. Reducing emissions would be expected to minimize lags, as well as slow the 44 velocities of warming and required responses everywhere.

was involved, and the value often used is 9.8°C km −1 ). However, if the air condenses 72 moisture as it cools, it gains some heat of condensation, which slows the rate of 73 cooling. Thus, realized or "moist adiabatic" lapse rates (MALR) are often lower (e.g. 74 6.5°C km −1 ) than the adiabatic lapse rate. As a consequence, temperature lapse rates 75 are determined by water vapor changes and latent heat release, which are linked to 76 surface temperature, elevation and moisture content in the region of interest 16 . 77 Importantly, the limited data available indicates that temperature lapse rates 78 differ between regions and seasons 16,17 . For example, lapse rates range from 4.1°C 79 km −1 to 6.8°C km −1 for different times and seasons along the slopes of the central 80 Himalayas 18 . In the Washington Cascades mountains in North America, the values of 81 lapse rates vary from 2.5 km −1 to 7.5°C km −1 in late summer and spring 19 . 82 Nonetheless, the mechanisms underlying temperature lapse rate are well established 83 in thermodynamic theory, with MALR being a function of local temperature and 84 vapor pressure 20 . Here, we use this knowledge to assess the velocity of mountain 85 climate changes on a global scale, and then evaluate whether species are keeping pace 86 with shifting isotherms. Mountain areas were represented by grids with 0.5° 87 resolution of mean elevation above 1,000 m 21 (Figs. 1,2). The velocity is represented 88 by the vertical isotherm shifts (i.e. how far an isothermal line rises from the surface), 89 calculated by dividing the rate of warming in each mountain grid by its MALR, 90 derived from terrain-surface conditions. We estimate MALR (usually noted as Г ) as: 91 where γ is the mixing ratio of the mass of water vapor per unit mass of dry air, which 93 is also influenced by air pressure; T represents the air temperature (other parameters pressure at terrain surface from 2011 to 2015 from CRU TS 3.24 to calculate MALR 96 and consider the rate of warming from 1971 to 2015 for each mountain grid. The 97 MALR formula predicts that higher surface temperature or water vapor leads to 98 curvilinear decreases in MALR (Fig. 1a). Given the widely varied temperature and 99 water vapor differences between mountain grids (Fig. 1 b,c), our estimates show that 100 temperature lapse rates of global mountains vary from 3 to 9°C km -1 (Fig 1d). 101 Given these MALR estimates and the rates of surface temperature change ( Fig. 2  102 a,b), the climate velocities in global mountains vary considerably from -16.67 m yr -1 103 (Saltillo, Mexico) to +16.80 m yr -1 (Mashhad, Iran), with global average elevational 104 isotherm increases of +5.42 ± 0.03 m yr -1 in mountain ranges (Fig. 2c; Extended Data 105 Fig. 1). We identify high climate-velocity mountains where the isotherms have shifted 106 more than one standard deviation of the global mean value (higher than 8.45 m yr -1 ) 107 ( Fig. 2d). Defined in this way, about 16% of mountainous areas are exposed to high 108 climate velocity and many are well-recognized biodiversity hotspots, such as 109 Northern Sumatra, Hengduan, Nepal, Southern Ghats, Madagascar,Mediterranean 110 basin, Northern Sahara, Brazilian Highlands, Northern Andes, Costa Rica, and 111 Western America and Mexico (Supplementary Table 1). Other high velocity 112 mountains may also be of concern, including most of the Middle East mountains 113 (which includes endemics) and central Asia, Siberia, and the Alaska-Yukon region 114 Table 1). The approach delineates that high climate velocity can be 115 due to high rates of temperature warming or to low lapse rates. It is intuitive that 116 higher warming rate results in high climate velocity (Extended Data Fig. 2a) and this 117 is often observed in drier regions with low mean annual precipitation (MAP), possibly 118 due to the limited heat capacity of arid regions [22][23][24] (Supplementary Table 1). Yet, high 119 climate velocities are also observed in mountain ranges that have relatively modest 120 levels of warming but low lapse rates, which can occur at areas with high water vapour pressure and/or high surface temperatures (Extended Data Fig. 2b,and 122 Supplementary pressure, interacting with mountain topographies and levels of regional warming. 127 We also applied our framework to mountains on islands that harbor a high 128 proportion of endemic species and provide refuges for shifting species otherwise 129 constrained by the ocean 25,26 . We define islands as landmasses that are smaller than 130 Australia and surrounded by water 21 . We examined the velocity of climate change in 131 island mountains globally (only 14 islands have mountain grid cells higher than 1,000 132 m, based on our dataset) and found a mean shift of 3.35 m yr -1 , which was lower than 133 the global average (5.42 m yr -1 ). Two island regions-Madagascar and Japan -were 134 found to have high climate velocity (>8.45 m yr -1 ) ( Fig. 3a-c). In Japan, high climate 135 velocity is mainly caused by the surface temperature warming (Fig. 2b  Our study provides estimates of the rates of elevational shifts required by species 175 and whether they tracked (kept pace with) temperature changes precisely or lagged 176 behind. We compared these rate estimates with observed elevational shifts in montane 177 species, taken from peer-reviewed articles reporting multi-species elevational 178 redistributions (Supplementary Table 2). The global variation in predicted elevation 179 shifts in different mountain ranges shows that rates of expected species uplifts are 180 much lower (e.g. in Italy and France) than predicted in the original paper using widely 181 applied empirical values (range from 5.5°C km -1 to 6.5°C km -1 ) or a globally-182 averaged lapse rate of 5.5°C km -1 (Extended Data Fig. 6), and this may partially 183 explain previous findings that mountain species appeared to be generally lagging 184 behind climate changes 2 . To assess the probability of species tracking climate 185 velocity, considering each taxon at each region as one data-point, we conducted 186 bootstrapping to control the sample size effect and compare the probability of tracking 187 in relation to climate velocity (Fig. 4a). Each dataset was subsampled and the mean 188 shift was compared with its corresponding climate velocity, using Wilcoxon signed-189 rank test at significance level of 0.05. If the shifts of species are not significantly 190 different from the climate velocity, we consider the species successfully track the 191 climate velocity. We repeated the procedure for 1000 times (i = 1000) and used the Our estimation of global temperature lapse rates, based on the MALR formula 219 considering the latent heat release and water vapor changes, provides a heuristic 220 understanding of climate velocity in global mountains. However, many other 221 mechanisms, such snow albedo, radiative flux changes, surface heat loss and aerosols, 222 also influence the energy balance regimes, making a direct estimate of climate 223 conditions and climate change in extremely difficult in mountain regions 10 . endangered species, our study helps identify vulnerable regions with high climate 226 velocity, which we suggest are priority regions for conservation. Given this 227 vulnerability, extensive monitoring networks for both mountain climate and biological 228 impacts are urgently needed. 229

Methods 230
The climatic data sources and the calculation of MALR 231 The climatic data, including mean annual temperature and water vapor pressure, and 232 the corresponding global digital elevation model were derived from the gridded CRU 233 TS3.24 database (0.5° resolution), which we averaged over every 5 years. Both mean 234 annual temperature and water vapor pressure were derived from local weather stations 235 and subsequently averaged across coarse spatial extent to obtain the final values 24 . 236 MALR at each grid was generated by the MALR formula: 237 where Г is the moist adiabatic lapse rate (K/m), denotes the Earth's gravitational 239 acceleration (9.8076 m/s 2 ), denotes the heat of water vaporization (2,501,000 240 J/kg), denotes the specific gas constant of dry air (287 J kg −1 K −1 ), denotes 241 the dimensionless ratio of the specific gas constant of dry air to the specific gas 242 constant for water vapour (0.622), denotes the specific heat of dry air at constant 243 pressure (1,005 J kg −1 K −1 ), and denotes the air temperature (K). is the mixing 244 ratio of the mass of water vapor to the mass of dry air: 245 where represents the water vapor pressure of the air and represents the pressure 247 of the air. Here, was derived from the Barometric formula: 248 where denotes the static pressure (101,325.00 pascals), denotes the molar 250 mass of Earth's air (0.0289644 kg/mol), denotes the universal gas constant for air 251 (8.31432 N m mol −1 K −1 ), ℎ denotes the height above sea level (meters), and The preprocessing of climatic variables (from monthly data to annual data) and the 254 calculation of basic climate velocity and climate velocity was computed by pySpark. 255 Islands are defined as landmasses that are smaller than Australia and surrounded by 256 water 21 . Here, the input dataset is not the digital elevation model (DEM) from CRU 257 but from SRTM, which is more conservative on defining terrestrial area (smaller 258 terrestrial regions), so it is better for island detection (the islands near shores are not 259 connected to continents). Greenland is not included as it is not surrounded by the 260 ocean in the dataset. These analyses were run in Wolfram Mathematica 9. Though we 261 only present the results when the anthropogenic warming accelerated (Extended Data 262 The surface aspect of each pixel was calculated by using the Surface Aspect Tool in 267 ArcGIS Pro (the license of 3D analyst is required), and the input digital elevation 268 model (DEM) is from CRU in order to match the climatic data. Two characters 269 (elongation and orientation) were calculated for each mountain ranges previously 270 defined in the literature 36 in Wolfram Mathematica 9. Elongation is defined as 1-271 (smallest axis of the best-fit ellipse / largest axis of the best-fit ellipse), and 272 orientation is computed as the angle between the largest axis and the horizontal axis. 273 The orientation of only those mountains having elongation value greater than 0.5 were 274 further analyzed. To provide a better alignment with the definition of global 275 mountains described in the main text, we then analyzed climatic data from elevations 276 higher than 1,000 m (a.s.l.) within these 'expert-identified' mountains. For the clarity, 277 we summarized the results based on categories of surface aspects and orientations of to statistically quantify the differences across different aspects and orientations of 280 mountains. A Wilcoxon signed-rank test was then applied in paired-dataset 281 comparisons as a post Hoc analysis. 282 283

Biological datasets 284
We adapted published studies providing range-or boundary-shift information based 285 on an exhaustive literature review (Supplementary Table 2 Nevada] 38 , the number of samples in two investigating periods were 1,168 V.S. 292 29,174); (2) the information provided in a large-scale (exceed 5 x 5 degree on the 293 map) research was insufficient, so we could not divide the records into geographical 294 regions (e.g., 39 ). Along with the raw datasets we used in our analysis, statistically 295 summarized information for all literature is also provided. The standard error was not 296 provided for studies with (1) only one record in a region or (2) insufficient 297 information from the original paper. 298 299

Corresponding climate velocity to the biological data 300
To pair biological records with the climate velocity derived from MALR, climate 301 velocities were assigned or statistically summarized based on the following criteria: 302 (1) if the spatial scale of a research was less than 1 x 1 degree (i.e. distributed within a 303 grid), the climate velocity at the same corresponding grid was used (eg. 40 ); (2) if the 304 scale of a study exceeded 1 x 1 degree (i.e., distributed on multiple grids), the regional climate velocity was statistically calculated (provided as mean, variance, and sample 306 size; e.g. 41 ); (3) if the study region exceeds 5 x 5 degrees, the biological records were 307 grouped into different regions (e.g. 8 ), and the climate velocity of each region was 308 derived based on criteria 1 and 2; (4) for studies encompassing multiple regions and 309 periods, the corresponding spatiotemporal information was used in order to derive the 310 most accurate climate velocity (eg., 42 ); (5) if a study investigated the same region 311 multiple times in different temporal periods, the same geographical information was 312 used (e.g. 43 ). 313

314
The probability of species tracking climate velocity-comparing biological data 315 and climate velocity 316 In order to statistically compare the biological data with corresponding climate 317 velocity, we calculated the probability of species tracking climate velocity after 318 bootstrapping the data to meet a relatively constant sample size across regions and 319 taxa. Some studies did provide detailed raw data (Supplementary Table 2), but for 320 those only reporting statistical results-such as mean and variance/standard 321 deviation/error (Extended Data Fig. 9)-we applied different probability distribution 322 functions, normal distribution as well as non-normal distributions (log-normal and 323 Student-t distribution), in our statistical analysis according to the statistical 324 characteristics (mean and standard deviation/error) provided in the original reports. 325 The nature of non-normally-distributed frequency in our dataset ( Supplementary Fig.  326 3) was taken cared here, and the final probability reported was averaged from all three 327 approaches with different probabilities. 328 The probability of species tracking climate velocity was then calculated as follows: 329 First of all, we used the bootstrap technique to subsample the dataset to control the 330 inconsistencies induced by having different sample sizes across studies. For each taxon in each region, we set the sample size to n and drew n records (n in Fig. 4a). If 332 the total number of records for that taxon in that region is smaller than n, all records 333 were used. For those only reporting statistical results (21 out of 47), we applied 334 different probability functions to generate the drawn value as described in the 335 previous paragraph. The drawn biological data were then compared to the 336 corresponding climatic velocity using the non-parametric method-Wilcoxon signed-337 rank test-because many datasets did not satisfy the normal distributed assumption 338 ( Supplementary Fig. 3). This procedure (draw and comparison) was then iterated 339 1,000 times (i in Fig. 4a), and we calculated the number of iterations at which the 340 biological data showed no significant difference to the corresponding climate velocity 341 (i.e., did not meet the significant level, 0.05; p in Fig. 4 which implies the ability of 342 tracking climate change; the lower p, the higher ability) and divided it by the total 343 number of iterations (1,000; i in Fig. 4a); the result was the probability of species 344 tracking climate velocity. A logistic-type (probit) function was then applied to 345 estimate the probability curve. We also ran a sensitivity analysis by setting different 346 values for n (30, 40, 50, 60, 70, 80, 90, and 100), and the results indicated that n does 347 not influence the probability of species tracking climate velocity (Supplementary 348 Table 3), so we decided to set n = 30 to fairly address the small-sample-size research. 349 The data processing and statistical analysis in this section were done in R. 350 351

Probability of species tracking using different lapse rate calculations 352
The straightforward thought to test how better climate velocity derived from MALR 353 (with the consideration of water vapor) is tracked by species should be: applying 354 regression for water vapor and the residual that derived from the regression of (MALR formula and Fig. 1a), the analysis based on the concept of regression should 358 not work. In order to bridge the method, which is more familiar to readers, to our 359 formal analysis (which is going to be described later), we firstly provide a set of 360 intuitive but potentially biased histograms to show the different explanation powers 361 among different lapse rate calculations (Extended Data Fig. 10a). The residual 362 between mean observed shifting rate and the mean climate velocity derived from 363 different lapse rate calculations for a taxon are shown. We can still find that the 364 residual histogram of MALR is more normally distributed than that of others (i.e. it 365 explains biological dataset better than others), yet please note that this method 366 neglects the variance of different records of a taxon. Consequently, in order to 367 formally compare different lapse rate calculations, the method with probability scope 368 and subsampling is applied. Based on the results generated from the previous section,  Plant Ecol 195, 179-196 (2008). Science 320, 1768-1771 (2008). detected on a regional scale. Global Ecol Biogeogr 12, 403-410 (2003).