# Unexpected westward range shifts in European forest plants link to nitrogen deposition

Pieter Sanczuk, Kris Verheyen, Jonathan Lenoir, Florian Zellweger, Jonas Lembrechts, Francisco Rodríguez-Sánchez, Lander Baeten, Markus Bernhardt-Römermann, Karen de Pauw, Pieter Vangansbeke, et al.

► To cite this version:

Pieter Sanczuk, Kris Verheyen, Jonathan Lenoir, Florian Zellweger, Jonas Lembrechts, et al. Unex-pected westward range shifts in European forest plants link to nitrogen deposition. Science, 2024, 386 (6718), pp.193-198. 10.1126/science ado0878. hal-04734750

HAL Id: hal-04734750

https://u-picardie.hal.science/hal-04734750v1

Submitted on 10 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Unexpected westward range shifts in European forest plants links to nitrogen deposition

2 Pieter Sanczuk¹*, Kris Verheyen ¹, Jonathan Lenoir², Florian Zellweger³, Jonas J. Lembrechts⁴,⁵, Francisco

3 Rodriguez-Sanchez⁵, Lander Baeten¹, Markus Bernhardt-Römermann⁷,⁸, Karen De Pauw¹, Pieter

4 Vangansbeke¹, Michael P. Perring⁹,¹⁰, Imre Berki¹¹, Anne Bjorkman¹²,¹³, Jörg Brunet¹⁴, Markéta

5 Chudomelová¹⁵, Emiel De Lombaerde¹, Guillaume Decocq², Thomas Dirnböck¹⁶, Tomasz Durak¹⁷,

6 Caroline Greiser¹⁸,¹⁹, Radim Hédl¹⁵,²⁰, Thilo Heinken²⁰, Ute Jandt²²,²³, Bogdan Jaroszewicz²⁴, Martin

7 Kopecký²⁵,²⁶, Dries Landuyt¹, Martin Macek²⁵, František Máliš²⁷,²⁸, Tobias Naaf²⁹, Thomas A. Nagel³⁰, Petr

8 Petřík³¹,³², Kamila Reczyńska³³, Wolfgang Schmidt³⁴, Tibor Standovár³⁵, Ingmar Staude²³,³⁶, Krzysztof

9 Świerkosz³⁷, Balázs Teleki³⁸, Thomas Vanneste¹, Ondrej Vild¹⁵, Donald Waller³⁹, Pieter De Frenne¹

10 *corresponding author: Pieter Sanczuk; Pieter.Sanczuk@UGent.be

11 Affiliations:

12 ¹Forest & Nature Lab, Department of Environment, Ghent University, Melle-Gontrode, Belgium

13 ²UMR CNRS 7058 "Ecologie et dynamique des systèmes anthropisés" (EDYSAN), Université de Picardie

14 Jules Verne, Amiens Cedex 1, France

15 ³Forest Resources and Management, Swiss Federal Research Institute WSL, Birmensdorf, Switzerland

16 ⁴Research Center on Plants and Ecosystems (PLECO), University of Antwerp, Wilrijk, Belgium

17 ⁵Ecology & Biodiversity (E&B), Utrecht University, Utrecht, the Netherlands

18 ⁶Universidad de Sevilla, Departamento de Biología Vegetal y Ecología, Sevilla, Spain

19 ⁷Friedrich Schiller University Jena, Institute of Ecology and Evolution, Jena, Germany

20 ⁸iDiv, German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Leipzig, Germany

21 ⁹Environment Centre Wales, UKCEH (UK Centre for Ecology and Hydrology), Bangor, UK

22 ¹⁰The UWA Institute of Agriculture, The University of Western Australia, Perth, Australia

23 ¹¹Faculty of Forestry, Institute of Environmental and Earth Sciences, University of Sopron, Sopron, Hungary

25 ¹²Department of Biological and Environmental Sciences, University of Gothenburg, Gothenburg, Sweden

26 ¹³Gothenburg Global Biodiversity Centre, Gothenburg, Sweden

27 ¹⁴Southern Swedish Forest Research Centre, Swedish University of Agricultural Sciences, Lomma, Sweden

28 ¹⁵Department of Vegetation Ecology, Institute of Botany of the Czech Academy of Sciences, Brno, Czech Republic

30 ¹⁶Ecosystem Research and Environmental Information Management, Environment Agency Austria, Vienna, Austria

32 ¹⁷Institute of Biology, University of Rzeszów, Rzeszów, Poland

33 ¹⁸Department of Physical Geography, Stockholm University, Stockholm, Sweden
<table>
  <tr>
    <td>34</td>
    <td><sup>19</sup>Department of Forest Ecology and Management, Swedish University of Agricultural Sciences, Umeå, Sweden</td>
  </tr>
  <tr>
    <td>35</td>
    <td><sup>20</sup>Department of Botany, Palacký University in Olomouc, Olomouc, Czech Republic</td>
  </tr>
  <tr>
    <td>36</td>
    <td><sup>21</sup>Institute of Biochemistry and Biology, University of Potsdam, Potsdam, Germany</td>
  </tr>
  <tr>
    <td>37</td>
    <td><sup>22</sup>Institute of Biology/Geobotany and Botanical Garden, Martin-Luther-University Halle-Wittenberg, Halle/Saale, Germany</td>
  </tr>
  <tr>
    <td>38</td>
    <td><sup>23</sup>German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena Leipzig, Leipzig, Germany</td>
  </tr>
  <tr>
    <td>39</td>
    <td><sup>24</sup>Faculty of Biology, Bialowieza Geobotanical Station, University of Warsaw, Bialowieza, Poland</td>
  </tr>
  <tr>
    <td>40</td>
    <td><sup>25</sup>Department of Geoecology, Institute of Botany of the Czech Academy of Sciences, Průhonice, Czech Republic</td>
  </tr>
  <tr>
    <td>41</td>
    <td><sup>26</sup>Faculty of Forestry and Wood Sciences, Czech University of Life Sciences Prague, Praha 6 - Suchdol, Czech Republic</td>
  </tr>
  <tr>
    <td>42</td>
    <td><sup>27</sup>Department of Phytology, Technical University in Zvolen, Zvolen, Slovakia</td>
  </tr>
  <tr>
    <td>43</td>
    <td><sup>28</sup>National Forest Centre, Zvolen, Slovakia</td>
  </tr>
  <tr>
    <td>44</td>
    <td><sup>29</sup>Leibniz Centre for Agricultural Landscape Research (ZALF), Muencheberg, Germany</td>
  </tr>
  <tr>
    <td>45</td>
    <td><sup>30</sup>Department of forestry and renewable forest resources, Biotechnical Faculty, University of Ljubljana, Ljubljana, Slovenia</td>
  </tr>
  <tr>
    <td>46</td>
    <td><sup>31</sup>Department of vegetation ecology, Czech Academy of Sciences, Institute of Botany, Brno, Czech Republic</td>
  </tr>
  <tr>
    <td>47</td>
    <td><sup>32</sup>Department of Ecology, Faculty of Environmental Sciences, Czech University of Life Sciences Prague, Czech Republic</td>
  </tr>
  <tr>
    <td>48</td>
    <td><sup>33</sup>Independent Researcher, Wrocław, Poland</td>
  </tr>
  <tr>
    <td>49</td>
    <td><sup>34</sup>Department of Silviculture and Forest Ecology of the Temperate Zones, University of Goettingen, Goettingen, Germany</td>
  </tr>
  <tr>
    <td>50</td>
    <td><sup>35</sup>Department of Plant Systematics, Ecology and Theoretical Biology, Institute of Biology, ELTE Eötvös Loránd University, Budapest, Hungary</td>
  </tr>
  <tr>
    <td>51</td>
    <td><sup>36</sup>Institute of Biology, Leipzig University, Leipzig, Germany</td>
  </tr>
  <tr>
    <td>52</td>
    <td><sup>37</sup>Museum of Natural History, University of Wrocław, Wrocław, Poland</td>
  </tr>
  <tr>
    <td>53</td>
    <td><sup>38</sup>HUN-REN-UD Biodiversity and Ecosystem Services Research Group, University of Debrecen, Debrecen, Hungary</td>
  </tr>
  <tr>
    <td>54</td>
    <td><sup>39</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>55</td>
    <td><sup>40</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>56</td>
    <td><sup>41</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>57</td>
    <td><sup>42</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>58</td>
    <td><sup>43</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>59</td>
    <td><sup>44</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>60</td>
    <td><sup>45</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>61</td>
    <td><sup>46</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>62</td>
    <td><sup>47</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>63</td>
    <td><sup>48</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
  <tr>
    <td>64</td>
    <td><sup>49</sup>Botany, University of Wisconsin - Madison, Madison, USA</td>
  </tr>
</table>
**Abstract:** Climate change is commonly assumed to induce species' range shifts towards the poles. Yet, other environmental changes may affect the geographical distribution of species in unexpected ways. Here we quantify multi-decadal shifts in the distribution of European forest plants, and link these shifts to key drivers of forest biodiversity change: climate change, atmospheric deposition (nitrogen and sulphur) and forest canopy dynamics. Surprisingly, westward distribution shifts were 2.6 times more likely than northward ones. Not climate change, but nitrogen-mediated colonization events, possibly facilitated by the recovery from past acidifying deposition, best explain westward movements. Biodiversity redistribution patterns appear complex and are more likely driven by the interplay among several environmental changes than due to the exclusive effects of climate change alone.

**Key words:** Acidification, atmospheric pollution, eutrophication, nitrogen deposition, climate change, conservation policy, isotherms, forestREplot, forest ecosystems, species range shifts, sulphur deposition, understorey plants

**One Sentence Summary:** Nitrogen deposition rather than climate change explains unexpected westward range centroid shifts of European forest plants
79

**Main text**

One of the most prominent biogeographical changes of the 21st century is the large-scale redistribution of plants and animals in response to changes in the climate system (1). Warming temperatures are causing many terrestrial species to move towards higher latitudes and elevations, resulting in a reordering of species' distributions (1–3) and the emergence of novel communities (4). Empirical evidence has been reported for a wide range of ecosystems and taxa - from poleward and upslope range shifts in temperate regions (5, 6) and high-latitude boreal biomes (7) to upslope shifts in mountain vegetation (8) – suggesting an emerging link with anthropogenic climate warming (2, 9).

According to the most recent global synthesis (1), terrestrial species are shifting towards higher latitudes at an average rate of 1.11 km year⁻¹. This trend, however, lacks statistical significance, possibly because estimates are often blurred by variation in methodological attributes (1, 10). Alternatively, species' redistributions in geographical directions that are orthogonal (i.e., west-east oriented) or even inverse (e.g., equatorward) to the moving isotherms are less likely to be detected from commonly studied range boundary shifts along thermal transects of latitude and elevation alone (11, 12). Indeed, other prominent environmental changes such as atmospheric (nitrogen and sulphur) deposition and forest disturbances show spatial patterns that are weakly correlated to the geographic direction of climate change (13–15), and can also influence demographic processes of colonization and local extinction (1, 12, 16). To what extent these other environmental changes are contributing to species range shifts remains largely unquantified (17–19).

Here we quantify the rate and geographic direction of range shifts in 266 European forest understory plant species using multi-decadal community data collected in mature forest stands across 2,954 resurveyed semi-permanent vegetation plots (20) (Fig. 1A). Plant community data were derived from baseline surveys recorded between 1933 and 1994 and paired resurveys carried out after the baseline surveys between 1987 and 2017 (median [min – max] inter-survey interval: 39 [13 – 67] years). We quantified shifts of species' distributions within the spatial extent of the study area based on range centroids, i.e., the abundance-weighted geometric center of a species' distribution (Fig. S1). In contrast to the more frequent quantification of range boundary shifts at the trailing or leading edges, analyzing centroid shifts allows to obtain more robust estimates of the magnitude and geographic direction of complex distribution shifts (6, 11, 21). This is important, because range shift estimates from leading and trailing edges alone are more prone to bias from stochastic processes or low sample sizes that may blur overall biogeographical trends (11).

109

**The rate and geographic direction of centroid shifts**
To calculate centroid shift of each species, we first located the position of the abundance-weighted range centroid at the time of the baseline survey and the resurvey, and assessed the magnitude (i.e., the distance) and geographical direction (i.e., the bearing) of the centroid shift over time. Centroid shifts were expressed as the absolute shift rate (km year⁻¹) as well as the projected south-north (km north year⁻¹) and west-east (km east year⁻¹) rate (schematically explained in Fig. S1). Centroid shifts were calculated for the 266 species that were recorded in ≥1% of the plots to increase robustness of the estimates. The directionality (i.e., angular dispersion of the directions of centroids shifts) across all species was tested using the Rayleigh's *r* coefficient, a circular regression coefficient which quantifies how uniform and isotropic the directions of centroid shifts are (Rayleigh's *r* = 1 if all species are moving in the same direction while Rayleigh's *r* = 0 with random directional movements, i.e., anisotropic, *meaning that* directions of shift can be drawn from a uniform circular distribution).

Species' centroid shifts were first compared to the velocity and direction of climate change realized over the course of the study period. Spatially explicit climate change velocities were calculated by climate analogue mapping (22, 23), an approach that is theoretically equivalent to the mapping of species centroid shifts (Fig. S1). In contrast to the frequent calculation of climate-change velocities based on gradients of isotherms alone (3, 24), climate analogue mapping allows to consider consolidated changes of multiple bioclimatic variables at the same time. For example, we here simultaneously consider changes in maximum growing-season temperatures, minimum winter temperatures and growing-season precipitation as one measure of the climate change velocity between the baseline survey and resurvey periods (Fig. S2). This is highly relevant because plants respond not only to warming temperatures but also to alterations in precipitation regimes. In this method, for all resurveyed vegetation plots, a grid search (at 4 km × 4 km resolution) was performed to map all raster cells within the study area in which the climatic conditions in the resurvey period are similar (i.e., show no statistical difference) to a given plot's climate during the baseline period (i.e., 'analogue climate conditions'). For each plot, we then located the position of the nearest raster cell with analogue climatic conditions to calculate the velocity and geographic direction of climate change over time. Identical to the centroid shifts, the climate change velocity for each plot was expressed as the absolute shift velocity (km year⁻¹), as well as the projected south-north (km north year⁻¹) and west-east (km east year⁻¹) velocity. The directionality of climate analogue shifts was tested using Rayleigh's *r* coefficient as described above.

Centroid shifts across the 266 understory plant species varied between 0.006 km year⁻¹ (*Symphytum cordatum*) to 18.27 km year⁻¹ (*Abies alba* seedlings), and occurred at a mean rate of 3.56 (5% – 95% quantile: 0.39 – 9.80) km year⁻¹ (Fig. 2A, Data S1). Surprisingly, two-thirds of the *studied* plant species showed directional shifts along the west-east axis (Rayleigh's *r* = 0.23; *df* = 265; *p* < 0.05). Most of these
shifts were westward (39% of species), but we also noted many eastward shifts (23%). Southward shifts (23%) were more frequent than northward shifts (15%). Westward range **centroid shifts** were thus **2.6 times** more likely than the northward range shifts expected in response to climate change. The average south-north rate of centroid shifts was slow but significantly equatorward (-0.63 [-4.30 – 2.89] km north year⁻¹; one-sample t-test: t = -4.36, *df* = 265, *p* < .001), while the rate of west-east shifts was **1.8 times** faster and significantly westward (-1.17 [-6.95 – 4.17] km east year⁻¹; one-sample t-test: t = -4.90, *df* = 265, *p* < .001). The observed rates of centroid shifts **towards each cardinal direction** were **minimum 62% (southward shifts) and maximum 70% (eastward shifts)** faster than expected by chance as confirmed by a null-model approach (**Fig. S5**).

The climate significantly changed over the course of the study period in 2,949 of the 2,954 resurveyed vegetation plots (99.8%) based on climate analogue mapping. Maximum growing-season temperatures increased by an average of 1.59 [1.15 – 2.21] °C. Climate change took place at an average absolute velocity of 0.66 (0.07 – 1.67) km year⁻¹ (**Fig. 2B**). Unsurprisingly, the dominant geographic direction of climate change was north (40% of the plots; Rayleigh’s *r* = 0.29; *df* = 2,948; *p* < 0.05). The south-north velocity of climate change was also significantly poleward (0.24 [-0.72 – 1.94] km north year⁻¹; one-sample t-test: t = 18.15, *df* = 2,948, *p* < .001). The west-east climate change velocity was marginal but significantly eastward (0.06 [-0.92 – 0.72] km east year⁻¹; one-sample t-test: t = 5.55, *df* = 2,948, *p* < .001) **and thus opposite to the most common cardinal direction of centroid shifts of European forest plants**.

Neither the geographic direction nor the velocity of climate change was reflected in the species' centroid shifts, **hinting towards the importance of other environmental changes**. Two prominent alternative drivers of forest plant community changes are elevated atmospheric inputs of nitrogen and forest canopy cover changes (16, 25–28), **with eutrophying effects of nitrogen complicated by recovery dynamics from past acidification caused by the combined deposition of nitrogen and sulphur compounds** (29, 30). These drivers show spatial patterns not confounded with the velocity of climate change (across all plots, pairwise Spearman correlations between the south-north velocity of climate change *versus* the rates of nitrogen deposition and forest canopy cover changes were only 0.04 and 0.01, respectively).

To better understand the potential drivers associated with **the reported changes in the geographical distribution of species**, we related the rates of centroid shifts to the average rate of nitrogen deposition between the baseline survey and resurvey (kg N ha⁻¹ year⁻¹) extracted from **atmospheric deposition maps** at 0.1° resolution (c. 8 km × 8 km within the study area) (**Fig. 1A**), and to the observed rate of forest canopy cover change at each site (% canopy cover increase year⁻¹), while also accounting for the velocity of climate change (the absolute, south-north and west-east velocity) in a linear mixed-effects modeling framework. **We furthermore tested for the potential confounding effect of past acidifying deposition, considering the**
known adverse effects on European forest plant communities (31). The acidifying deposition rate, however, was highly correlated to the rate of nitrogen deposition over the course of the study due to partially shared emission sources (Spearman correlation: 0.87; n = 2,954 plots), and their individual effects are therefore difficult to tease apart in an observational study. We calculated species-specific experienced rates of atmospheric (nitrogen and acidifying) deposition and forest canopy cover changes as the average rate across all plots where the species was observed, weighted by the species' original abundance in the baseline time period survey (Fig. S6 for a data flow chart).

Model outputs show that the absolute rate of centroid shifts was weakly but exclusively linked to the rate of forest canopy change, with greater opening of the canopy enhancing centroid shifts (Fig. 3). The velocity of climate change was not associated to the rates of centroid shifts. In contrast, the rate of nitrogen deposition was significantly linked to the west-east rate of centroid shifts, with species that initially experienced a lower nitrogen deposition rate across their distributions showing faster westward shifts (Fig. 3A). Variation partitioning revealed that the nitrogen deposition rate rather than the climate change velocity explained most of the variation in the species' centroid shifts, albeit the proportion of variation explained was small (Figs. 3B, S7). The estimated effects of past acidifying deposition on species centroid shifts were nearly identical to the effects of the spatially correlated nitrogen deposition rate (Fig. S8). We are therefore unable to distinguish with certainty whether centroid shifts were brought about by eutrophying rather than changes in acidifying deposition, or a combination of both. In either case, however, atmospheric deposition rates – and not the climate change velocities – were the superior predictors of westward species movements.

## Colonization and extinction centroids

To shed light on the mechanisms driving centroid shifts, we decomposed centroid shifts into shifts attributed to the individual contribution of colonization and local extinction. For this analysis, we introduce the concept of colonization centroids (the centroid of plots newly colonized by a species, abundance-weighted by the percentage cover in the resurvey) and extinction centroids (the centroid of plots in which a species became extinct, abundance-weighted by the percentage cover in the baseline survey). Colonization and extinction centroids were expressed as the projected distance from the species' baseline centroid position in each geographic direction (*km north* and *km east*). From a biogeographical point of view, longer distances reflect that colonization or local extinction events took place farther away from the baseline centroid, suggesting that these processes occurred mostly in one preferred direction (schematically explained in Fig. S1C).

The average absolute distance of colonization centroids (202.20 [28.30 – 478.46] km; *n* species = 202) was larger than the average absolute distance of extinction centroids (82.22 [4.22 – 249.48] km; *n* species =
246; Figs. 4, S9, S10). Colonization centroids were also more isotropic across species (Rayleigh's *r* = 0.32; *df* = 201; *p* < 0.05) than extinction centroids (Rayleigh's *r* = 0.11; *df* = 245; *p* < 0.05). This suggests that colonization events were happening more in one preferred direction and occurred further away from baseline range centroids (i.e., closer to range boundaries) compared to local extinctions.

Colonization along the west-east axis was most closely related to the nitrogen deposition rate (Figs. S11, S12). Because the observed dominant direction of species' colonization was westward (39% of the species colonized west, 21% east, 28% south and 12% north), colonization occurred more frequently for species that initially experienced a lower rate of nitrogen depositions across their distributions. This pattern is possibly associated to the westward colonization of nitrogen generalist species that can take advantage of eutrophic conditions such as observed in large parts of western Europe (Figs. 1A, S13). Indeed, linking the colonization centroids to each species' ecological indicator value for nitrogen niche width (with larger values indicating generalist species with a broader niche) (32) revealed that for nitrogen generalists in particular, those that initially occurred in areas with a lower rate of nitrogen deposition have taken advantage to move more westward (Figs. 5, S14). Nitrogen generalists that already occurred in areas with higher nitrogen deposition (western Europe, Figs. 1A, S13) tended to remain in place without necessarily moving westward. More specialist species, i.e. those with narrow niche widths for nitrogen and that often also have smaller range sizes (33), have shown lower colonization rates across temperate Europe, allowing generalist species to replace specialists (14). Also the decreasing levels of acidifying deposition (since the peak in the 1980s (34)) may have facilitated the recovery of species' ranges in formerly polluted regions (30). Using our observational data, we cannot fully disentangle these recovery effects following past acidification caused by both nitrogen and sulphur pollution from dynamics of eutrophication chiefly involving nitrogen deposition. However, eutrophication may be the most likely driver, because (i) we show that west-east colonization distances were statistically better linked to nitrogen deposition and the species' nitrogen niche width than to acidifying deposition and the species acidity niche width (Fig. S15); and (ii) the vast majority of the vegetation plots appear to be relatively well buffered against soil acidification (Fig. S16). Regardless of whether the driver of westward colonization chiefly involves eutrophying deposition or facilitated by the recovery from past acidification, forest plant species native to regions with lower deposition rates are more vulnerable to unanticipated range shifts in response to atmospheric pollution – a key finding for forest biodiversity conservation policy.

Local extinction events along the south-north axis were preferentially located southward relative to the species' baseline range centroid position, and thus closer to species' warm range limits (25% south versus 13% north) (Fig. 4). This trend was significantly associated with climate change and a higher rate of nitrogen deposition (Figs. S11, S12). Local extinction events along the west-east axis, however, occurred
more often (28% east and 31% west). Eastward local extinctions occurred more commonly in species that experienced a higher rate of nitrogen deposition across their distribution. Such nitrogen-mediated local extinctions were amplified when forest canopies became more open. Velocities of climate change also interacted with the rate of forest canopy cover change, in that local extinctions due to climate change occurred more often in forests where the canopy cover decreased. This confirms the importance of tree canopies buffering the impacts of environmental changes (35).

Our findings suggest that **atmospheric deposition** and forest canopy cover dynamics interact to determine how forest plant species are shifting their ranges, and that these environmental changes induce shifts that can be independent from isotherm shifts. This contradicts the idea that species have shifted ranges mainly in response to warming air temperatures. Rather, other environmental changes, especially rates of **atmospheric deposition** and forest canopy cover dynamics, have likely induced unexpected westward range shifts in European forest plants. Although it remains unclear whether the effect of atmospheric deposition chiefly involves eutrophication or a recovery effect from past acidification due to both nitrogen and sulphur pollution, our findings point at nitrogen deposition as the most likely driver explaining the westward range shifts in European forest plants. Since the continued success of the United Nations air convention (CLRTAP) and the European Union Emissions Ceiling Directive (NECD) in reducing nitrogen and sulphur emission levels, prospective trends in climate change and atmospheric deposition are unlikely parallel, with climate change outpacing the effects of atmospheric deposition on **future** species' range shifts. Accurate and recent species range shift data will be key to adequately anticipate the respective impacts of climate change and atmospheric deposition on biodiversity and ecosystem functioning. It is already clear, however, that biodiversity redistribution patterns appear complex and are more likely driven by the interplay among several environmental changes than due to the exclusive effect of climate change alone.
263 **References:**

264 1. J. Lenoir, R. Bertrand, L. Comte, L. Bourgeaud, T. Hattab, J. Murienne, G. Grenouillet, Species better track climate warming in the oceans than on land. *Nat. Ecol. Evol.* **4**, 1044–1059 (2020).

266 2. I. Chen, J. K. Hill, R. Ohlemüller, D. B. Roy, C. D. Thomas, Rapid Range Shifts of Species of Climate Warming. *Science* (80-.). **333**, 1024–1027 (2011).

268 3. M. T. Burrows, D. S. Schoeman, L. B. Buckley, P. Moore, E. S. Poloczanska, K. M. Brander, C. Brown, J. F. Bruno, C. M. Duarte, B. S. Halpern, J. Holding, C. V Kappel, W. Kiessling, M. I. O’Connor, J. M. Pandolfi, C. Parmesan, F. B. Schwing, W. J. Sydeman, A. J. Richardson, The Pace of Shifting Climate in Marine and Terrestrial Ecosystems. *Science* (80-.). **334**, 652–656 (2011).

272 4. J. W. Williams, S. T. Jackson, Novel climates, no-analog communities, and ecological surprises. *Front. Ecol. Environ.* **5**, 475–782 (2007).

274 5. C. Parmesan, N. Ryrholm, C. Stefanescu, J. K. Hill, C. D. Thomas, H. Descimon, B. Huntley, L. Kaila, J. Kullberg, T. Tammaru, W. J. Tennent, J. A. Thomas, M. Warren, Poleward shifts in geographical ranges of butterfly species associated with regional warming. *Nature* **399**, 579–583 (1999).

277 6. J. Lenoir, P. A. Marquet, P. De Ruffray, H. Brisse, A Significant Upward Shift in Plant Species Optimum Elevation During the 20th Century. *Science* (80-.). **320**, 1768–1772 (2008).

279 7. R. Virkkala, A. Lehikoinen, Patterns of climate-induced density shifts of species: poleward shifts faster in northern boreal birds than in southern birds. *Glob. Chang. Biol.* **20**, 2995–3003 (2014).

281 8. S. Dullinger, A. Gattringer, W. Thuiller, D. Moser, N. E. Zimmermann, A. Guisan, W. Willner, C. Plutzar, M. Leitner, T. Mang, M. Caccianiga, T. Dirnböck, S. Ertl, A. Fischer, J. Lenoir, J. Svenning, A. Psomas, D. R. Schmatz, Extinction debt of high-mountain plants under twenty-first-century climate change. *Nat. Clim. Chang.* **2**, 6–9 (2012).

285 9. C. Parmesan, G. Yohe, A globally coherent fingerprint of climate change impacts across natural systems. *Nature* **421**, 37–42 (2003).

287 10. M. A. Rubenstein, S. R. Weiskopf, R. Bertrand, S. L. Carter, L. Comte, M. J. Eaton, C. G. Johnson, J. Lenoir, A. J. Lynch, B. W. Miller, T. L. Morelli, M. A. Rodriguez, A. Terando, L. M. Thompson, Climate change and the global redistribution of biodiversity: substantial variation in empirical support for expected range shifts. *Environ. Evid.* **12**, 1–21 (2023).

291 11. L. P. Shoo, S. E. Williams, J. Hero, Detecting climate change induced range shifts: Where and how should we be looking? *Austral Ecol.* **31**, 22–29 (2006).

293 12. J. D. Ash, T. J. Givnish, D. M. Waller, Tracking lags in historical plant species' shifts in relation to regional climate change. *Glob. Chang. Biol.* **23**, 1305–1315 (2017).

295 13. T. L. Greaver, C. M. Clark, J. E. Compton, D. Vallano, A. F. Talhelm, C. P. Weaver, L. E. Band, J. S. Baron, E. A. Davidson, C. L. Tague, E. Felker-Quinn, J. A. Lynch, J. D. Herrick, L. Liu, C. L. Goodale, K. J. Novak, R. A. Haeuber, Key ecological responses to nitrogen are altered by climate change. *Nat. Clim. Chang.* **6**, 836–843 (2016).

299 14. I. R. Staude, D. M. Waller, M. Bernhardt-Römermann, A. D. Bjorkman, J. Brunet, P. De Frenne, R. Hédl, U. Jandt, J. Lenoir, F. Máliš, K. Verheyen, M. Wulf, H. M., Pereira, P. Vangansbeke, A. Ortmann-Ajkai, R. Pielech, I. Berki, M. Chudomelová, G. Decocq, T. Dirnböck, T. Durak, T. Heinken, B. Jaroszewicz, M. Kopecký, M. Macek, M. Malicki, T. Naaf, T. A. Nagel, P. Petřík, K. Reczyńska, F. H. Schei, W. Schmidt, T. Standovár, K. Świerkosz, B. Teleki, H. Van Calster, O. Vild, L. Baeten, Replacements of small- by large-ranged species scale up to diversity loss in Europe's temperate forest biome. *Nat. Ecol. Evol.* **4**, 802–808 (2020).

306 15. C. Senf, R. Seidl, Mapping the forest disturbance regimes of Europe. *Nat. Sustain.* **4**, 63–70 (2021).
307 16. P. Sanczuk, K. De Pauw, E. De Lombaerde, M. Luoto, C. Meeussen, S. Govaert, T. Vanneste, L. Depauw, J. Brunet, S. A. O. Cousins, C. Gasperini, P. Hedwall, G. Iacopetti, J. Lenoir, J. Plue, F. Selvi, F. Spicher, J. Uria-diez, K. Verheyen, P. Vangansbeke, P. De Frenne, Microclimate and forest density drive plant population dynamics under climate change. *Nat. Clim. Chang.* **13**, 840–847 (2023).

311 17. S. M. Crimmins, S. Z. Dobrowski, J. A. Greenberg, J. T. Abatzoglou, A. R. Mynsberge, Changes in Climatic Water Balance Drive Downhill Shifts in Plant Species’ Optimum Elevations. *Science* (80–). **331**, 324–327 (2011).

314 18. D. Jacobsen, The dilemma of altitudinal shifts: caught between high temperature and low oxygen. *Front. Ecol. Environ.* **18**, 211–218 (2020).

316 19. J. Lenoir, J. Ge, A. Guisan, P. Vittoz, T. Wohlgemuth, N. E. Zimmermann, S. Dullinger, H. Pauli, W. Willner, J. Svenning, Going against the flow: potential mechanisms for unexpected downslope range shifts in a warming climate. *Ecography (Cop).* **33**, 295–303 (2010).

319 20. D. M. Olson, E. Dinerstein, E. D. Wikramanayake, N. D. Burgess, G. V. N. Powell, E. C. Underwood, J. A. D. Amico, I. Itoua, H. E. Strand, J. C. Morrison, J. Loucks, T. F. Allnutt, T. H. Ricketts, Y. Kura, J. F. Lamoreux, W. Wesley, P. Hedao, K. R. Kassem, Terrestrial Ecoregions of the World: A New Map of Life on Earth. *Bioscience* **51**, 933–938 (2001).

323 21. J. Vanderwal, H. T. Murphy, A. S. Kutt, G. C. Perkins, B. L. Bateman, J. J. Perry, A. E. Reside, Focus on poleward shifts in species’ distribution underestimates the fingerprint of climate change. *Nat. Clim. Chang.* **2**, 1–5 (2012).

326 22. M. C. Fitzpatrick, R. R. Dunn, Contemporary climatic analogs for 540 North American urban areas in the late 21st century. *Nat. Commun.* **614**, 1–7 (2019).

328 23. J. W. Williams, S. T. Jackson, J. E. Kutzbach, Projected distributions of novel and disappearing climates by 2100 AD. *PNAS* **104**, 5738–5742 (2007).

330 24. S. R. Loarie, P. B. Duffy, H. Hamilton, G. P. Asner, C. B. Field, D. D. Ackerly, The velocity of climate change. *Nature* **462**, 1052–1055 (2009).

332 25. P. De Frenne, B. J. Graae, F. Rodríguez-Sánchez, A. Kolb, O. Chabrerie, G. Decocq, H. De Kort, A. De Schrijver, M. Diekmann, O. Eriksson, R. Gruwez, M. Hermy, J. Lenoir, J. Plue, D. A. Coomes, K. Verheyen, Latitudinal gradients as natural laboratories to infer species’ responses to temperature. *J. Ecol.* **101**, 784–795 (2013).

336 26. J. Segar, H. M. Pereira, L. Baeten, M. Bernhardt-römermann, P. De Frenne, N. Fernández, F. S. Gilliam, J. Lenoir, A. Ortmann-ajkai, K. Verheyen, D. Waller, B. Teleki, J. Brunet, M. Chudomelová, G. Decocq, T. Dirnböck, R. Hédl, T. Heinken, B. Jaroszewicz, M. Kopecký, Martin Macek, F. Máliš, T. Naaf, A. Orczewska, K. Reczynska, W. Schmidt, J. Šebesta, A. Stachurska-Swakoń, T. Standovár, K. Swierkosz, O. Vild, M. Wulf, I. R. Staude, Divergent roles of herbivory in eutrophying forests. *Nat. Commun.* **13**, 1–10 (2022).

341 27. K. Verheyen, L. Baeten, P. De Frenne, M. Bernhardt-ro, J. Cornelis, G. Decocq, H. Dierschke, O. Eriksson, T. Heinken, M. Hermy, P. Hommel, K. Kirby, G. Walther, M. Wulf, G. Verstraeten, Driving factors behind the eutrophication signal in understorey plant communities of deciduous temperate forests. *J. Ecol.* **100**, 352–365 (2012).

345 28. P.-O. Hedwall, J. U.-Diez, J. Brunet, L. Gustafsson, A.-L. Axelsson, J. Strengbom, Interactions between local and global drivers determine long-term trends in boreal forest understorey vegetation. *Glob. Ecol. Biogeogr.* **30**, 1765–1780 (2021).

348 29. N. Duarte, L. H. Pardo, M. J. Robin-Abbott, Susceptibility of Forests in the Northeastern USA to Nitrogen and Sulfur Deposition: Critical Load Exceedance and Forest Health. *Water, Air Soil Pollut.*, 1–21 (2013).

350 30. E. Tipping, J. A. C. Davies, P. A. Henrys, S. G. Jarvis, S. M. Smart, Long-term effects of atmospheric deposition on British plant species. *Environ. Pollut.* **281**, 117017 (2021).
352 353 G. Riofrio-Dillon, R. Bertrand, J.-C. Gégout, Toward a recovery time: forest herbs insight related to anthropogenic acidification. *Glob. Chang. Biol.* **18**, 3383–3394 (2012).

354 355 356 357 J. Dengler, F. Jansen, O. Chusova, E. Hüllbusch, M. P. Nobis, K. Van Meerbeek, I. Axmanová, H. H. Bruun, M. Chytrý, R. Guarino, G. Karrer, K. Moeys, T. Raus, M. J. Steinbauer, Y. Didukh, M. Diekmann, T. Englisch, E. Fernández-pascual, P. Julve, G. Nakhutsrishvili, W. A. Ozinga, E. Ruprecht, U. Šilc, J.-P. Theurillat, F. Gillet, Ecological Indicator Values for Europe (EIVE) 1.0. *Veg. Classif. Surv.* **4**, 7–29 (2023).

358 359 360 361 S. Kambach, J. Lenoir, G. Decocq, E. Welk, G. Seidler, S. Dullinger, J. Gégout, A. Guisan, H. Pauli, J. Svenning, P. Vittoz, T. Wohlgemuth, N. E. Zimmermann, H. Bruelheide, Of niches and distributions: range size increases with niche breadth both globally and regionally but regional estimates poorly relate to global estimates. *Ecography (Cop.).* **42**, 467–477 (2019).

362 363 P. Grennfelt, A. Engleryd, M. Forsius, Ø. Hov, H. Rodhe, E. Cowling, Acid rain and air pollution: 50 years of progress in environmental science and policy. *Ambio* **49**, 849–864 (2020).

364 35. P. De Frenne, Novel light regimes in European forests. *Nat. Ecol. Evol.* **8**, 196–202 (2024).

365 366 J. T. Abatzoglou, S. Z. Dobrowski, S. A. Parks, K. C. Hegewisch, TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958 – 2015. *Sci. Data* **5**, 1–12 (2018).

367 368 EMEP, EMEP Status Report 1/2023 “Transboundary Particulate Matter, Photo-Oxidants, Acidifying and Eutrophying Components” (2023).

369 38. European Environment Agency. https://www.eea.europa.eu/en/datahub/datahubitem-view/.

370 371 372 373 374 375 376 P. Sanczuk, K. Verheyen, J. Lenoir, F. Zellweger, J. Lembrechts, R.-S. Francisco, L. Baeten, M. Bernhardt-Römermann, K. De Pauw, P. Vangansbeke, M. Perring, I. Berki, A. Bjorkmann, J. Brunet, M. Chudomelová, E. De Lombaerde, G. Decocq, T. Dirnböck, T. Durak, C. Greiser, R. Hédl, T. Heinken, U. Jandt, B. Jaroszewicz, M. Kopecky, D. Landuyt, M. Macek, F. Máliš, T. Naaf, T. Nagel, P. Petřík, K. Reczynska, W. Schmidt, T. Standovar, I. Staude, K. Świerkosz, B. Teleki, T. Vanneste, O. Vild, D. Waller, P. De Frenne, Data and Code associated with original research article: Nitrogen deposition mediates westward range shifts in European forest plants [Dataset]. *Dryad*, doi: doi.org/10.5061/dryad.4b8gthtmt (2024).

377 40. P. De Frenne, F. Rodriguez-Sanchez, D. A. Coomes, L. Baeten, G. Verstraeten, M. Vellend, M. Bernhardt-Romermann, C. D. Brown, J. Brunet, J. Cornelis, G. M. Decocq, H. Dierschke, O. Eriksson, F. S. Gilliam, R. Hedl, T. Heinken, M. Hermy, P. Hommel, M. A. Jenkins, D. L. Kelly, K. J. Kirby, F. J. G. Mitchell, T. Naaf, M. Newman, G. Peterken, P. Petrik, J. Schultz, G. Sonnier, H. Van Calster, D. M. Waller, G.-R. Walther, P. S. White, K. D. Woods, M. Wulf, B. J. Graae, K. Verheyen, Microclimate moderates plant responses to macroclimate warming. *PNAS* **110**, 18561–18565 (2013).

383 384 385 386 387 F. Zellweger, P. De Frenne, J. Lenoir, P. Vangansbeke, K. Verheyen, L. Baeten, R. Hédl, I. Berki, J. Brunet, H. Van Calster, M. Chudomelová, G. Decocq, T. Dirnböck, T. Durak, T. Heinken, B. Jaroszewicz, M. Kepecky, F. Mális, M. Macek, M. Malicki, T. Naaf, T. A. Nagel, A. Ortmann-ajkai, P. Petrik, R. Pielech, K. Reczyn, W. Schmidt, R. Pielech, K. Reczynska, O. Vild, M. Wulf, D. Coomes, Forest microclimate dynamics drive plant responses to warming. *Science (80–).* **368**, 772–775 (2020).

388 42. S. A. Chamberlain, E. Szöcs, taxize: taxonomic search and retrieval in R. *F1000Research* **2**, 1–26 (2013).

389 390 J. Lenoir, J. Svenning, Climate-related range shifts – a global multidimensional synthesis and new research directions. *Ecography (Cop.).* **38**, 15–28 (2015).

391 392 D. D. Breshears, T. E. Huxman, H. D. Adams, C. B. Zou, J. E. Davison, Vegetation synchronously leans upslope as climate warms. *PNAS* **105**, 11591–11592 (2008).

393 394 45. Edzer Pebesma, Simple Features for R: Standardized Support for Spatial Vector Data. *R J.* **10**, 439–446 (2018).

395 396 R. J. Hijmans, geosphere: Spherical trigonometry. Retrieved from https://cran.r-project.org/package=geosphere. (2017).
397 47. R. J. Hijmans, J. Van Etten, raster: Geographic analysis and modeling with raster data. R package version 2.0-12. doi: http://cran.r-project.org/package=raster (2012).

399 48. Microsoft, S. Weston, foreach: Provides Foreach Looping Construct. R package version 1.5.1. (2020). https://cran.r-project.org/package=foreach.

401 49. Microsoft, S. Weston, doParallel: Foreach Parallel Adaptor for the "parallel" Package. R package version 1.0.16. (2020). https://cran.r-project.org/package=doParallel.

403 50. M. Tennekes, J. Nowosad, J. Gombin, S. Jeworutzki, K. Russell, R. Zijdeman, J. Clouse, R. Lovelace, J. Muenchow, tmap: Thematic Maps. *J. Stat. Softw.* **84**, 1–39 (2018).

405 51. C. Duprè, C. Stevens, T. Ranke, A. Bleekers, C. Peppler-Lisbach, D. Gowing, N. Dise, E. Dorland, R. Bobbink, M. Diekmann, Changes in species richness and composition in European acidic grasslands over the past 70 years: the contribution of cumulative atmospheric nitrogen deposition. *Glob. Chang. Biol.* **16**, 344–357 (2010).

409 52. D. Nieto-lugilde, J. Lenoir, S. Abdulhak, D. Aeschimann, S. Dullinger, J. Gégout, A. Guisan, H. Pauli, J. Renaud, J. Theurillat, W. Thuiller, J. Van Es, P. Vittoz, W. Willner, T. Wohlgemuth, N. E. Zimmermann, J. Svenning, Tree cover at fine and coarse spatial grains interacts with shade tolerance to shape plant species distributions across the Alps. *Ecography (Cop.)*. **38**, 578–589 (2015).

413 53. D. Bates, M. Mächler, B. M. Bolker, S. C. Walker, Fitting Linear Mixed-Effects Models Using lme4. *J. Stat. Softw.* **67**, 1–29 (2015).

415 54. C. F. Dormann, J. Elith, S. Bacher, C. Buchmann, G. Carl, G. Carr, J. R. Garc, B. Gruber, B. Lafourcade, P. J. Leit, M. Tamara, C. McClean, P. E. Osborne, B. S. Der, A. K. Skidmore, D. Zurell, S. Lautenbach, Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. *Ecography (Cop.)*. **27**, 27–46 (2013).

419 55. K. Barton, MuMln: Multi-Model Inference. (2017).

420 56. S. Nakagawa, H. Schielzeth, A general and simple method for obtaining R2 from generalized linear mixed-effects models. *Methods Ecol. Evol.* **4**, 133–142 (2013).

422 57. R Core Team. *A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing (Vienna, Austria, 2021).

424 58. B. Wen, H. Blondeel, L. Baeten, M. P. Perring, L. Depauw, S. L. Maes, L. De Keersmaeker, H. Van Calster, M. Wulf, T. Naaf, K. Kirby, M. Bernhardt-Römermann, T. Dirnböck, F. Máliš, M. Kopecký, O. Vild, M. Macek, R. Hédl, M. Chudomelová, J. Lenoir, J. Brunet, T. A. Nagel, K. Verheyen, D. Landuyt, Predicting trajectories of temperate forest understorey vegetation responses to global change. *For. Ecol. Manage.* **566** (2024).

429 59. B. Ulrich, "An ecosystem approach to soil acidification" in *Ulrich B., Sumner M.E. (Eds). Soil Acidity*. (Springer US, New York, 1991), pp. 28–79.

431
432 Acknowledgements

This paper is an outcome of the sREplot working group supported by sDiv, the Synthesis Centre of the
German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig (DFG FZT 118). **Funding:**
PS, PDF, and PV received funding from the European Research Council (ERC) under the European Union's
Horizon 2020 research and innovation programme (ERC Starting Grant FORMICA 757833) and Ghent
University grant BOF23/GOA/019. FZ received funding from the Swiss National Science Foundation
(project 193645). JJL received funding from the Flemish Research Organization (FWO, OZ7828, OZ7916
and OZ7792). KDP received funding from the Flemish Research Organization (FWO, ASP035-19 De
Pauw). IB was supported via grant EFOP-3.6.1.-16-2016-00018. MC was supported by the postdoctoral
fellowship L200052302 of the Czech Academy of Sciences, grant 21-11487S by the Czech Science
Foundation.MK, RH, MC, OV and MM were supported by the Czech Academy of Sciences (RVO
67985939). TD received funding from the eLTER PLUS project of the European Union Horizon 2020
(INFRAIA-01-2018-2019). CG received funding from FORMAS [project 2021-01993]. RH was supported
by the Czech Science Foundation project 21-11487S.FM was supported by the project APVV-19-0319. PP
was supported by the long-term development project of the Czech Academy of Sciences (RVO 67985939).
FRS was supported by VI PPIT – US. **Author contributions:** PS and PDF conceived the ideas and
designed the methodology, and with contributions of JL, JLL and FR; all authors collected data; PS
analysed data in collaboration with PDF; PS led the writing of the manuscript in collaboration with PDF,
and with significant contributions of KV, JL, FZ, JJL, FR, LB, MBR, KDP, PV and MPP. All authors
contributed critically to the draft and gave final approval for publication. **Competing interests:** None of
the authors have a conflict of interest. **Data and materials availability:** Raw macroclimate data can be
accessed via the TerraClimate database (36). Historical nitrogen and sulphur deposition data can be
accessed via the EMEP/CEIP 2023 Present state of emission data-base (37). Shapefile of biogeographical
regions are provided by the European Environment Agency (38). Species' ecological indicator values of
nitrogen niche width are provided by Dengler et al. (2023) (32). Raw data on species centroid shifts,
environmental change data, and R scripts to reproduce the methods, analyses and source code of all figures
are available on Dryad (39).

459
460 Supplementary Materials

461 Methods

462 Figs. S1 to S16

463 Data S1

464
465 **Figure Legends**

466 **Figure 1 | Spatial and environmental gradients of the study.** (A) Map of modelled nitrogen (N) 467 deposition rate (sum of oxidized and reduced wet and dry deposition expressed in kg N ha⁻¹ year⁻¹; dry 468 deposition accounted for deciduous forest surface) at 0.1° resolution for the reference year 2000, and 469 distribution of the 2,954 resurveyed vegetation plots (grey crosses, spatially jittered for clarity) across the 470 European temperate forest biome (shaded green background) (20). (B) Observed environmental variation 471 across the 2,954 vegetation plots (grey dots) of three key drivers of forest biodiversity over the course of 472 the study period investigated here: nitrogen deposition rate (total of oxidized and reduced wet and dry N), 473 climate-change velocity (realized changes of both temperature and precipitation, expressed in km north 474 year⁻¹) and the rate of canopy cover change (average annual rate; expressed in % canopy cover increase 475 year⁻¹). In all boxplots, we present the median (horizontal line), 1st and 3rd quantile (lower and upper hinges), 476 and 1.5 times the inter-quartile-range (whiskers). Half violin plots represent the density distributions of the 477 environmental change values. The grey dashed lines represent no changes (not shown for nitrogen 478 deposition rates). Negative values in the case of climate and canopy cover change indicate southward 479 velocities and canopy opening, respectively.

480 **Figure 2 | Rate and geographic direction of species centroid shifts and climate analogue shifts.** (A) 481 Rate and geographic direction of species range centroid shifts (n species = 266). (B) Velocity and 482 geographic direction of climate analogue shifts (n plots = 2,954). In all graphs, the Rayleigh's *r* statistic 483 represents a test of uniformity that compares the bearings of shifts to a uniform circular distribution (null 484 hypothesis). Larger values indicate more directional shifts. Asterisk (*) indicates significant deviations 485 from the null hypothesis (p < 0.05). See Figs. S3, S4 for results of the analyses including rare species and 486 per biogeographic region.

487 **Figure 3 | Effects of environmental changes on centroid shifts based on the most parsimonious model.** 488 (A) Results of the mixed-effects models (n species = 266) indicating coefficient estimates and 95% 489 confidence intervals (CI) of the effects of the velocity of climate change (km year⁻¹, km north year⁻¹, km 490 east year⁻¹), average nitrogen deposition rate between the baseline survey and resurvey (kg N ha⁻¹ year⁻¹), 491 and rate of canopy change (% cover increase year⁻¹), as well as all pairwise interactions (indicated by ‘×’) 492 on the absolute rate of centroid shifts (km year⁻¹) and the south-north rate (km north year⁻¹; negative values 493 indicate southward shifts) and west-east rate (km east year⁻¹; negative values indicate westward shifts) 494 extracted from the most parsimonious model structure (empty rows were not included in the final selected 495 model). All the predictor variables were z-transformed to increase comparability. Rates of canopy opening 496 (negative values of canopy change) are associated with greater absolute rates of centroid shifts. The west- 497 east rate was exclusively linked to nitrogen deposition, with faster westward shifts in species with lower 498 rates of nitrogen deposition across their distributions. Model fit is presented as the proportion of variation 499 explained by the fixed effect (marginal R², R²<sub>m</sub>) and the proportion of variation explained by the fixed and 500 random effects (conditional R², R²<sub>c</sub>). Models accounted for plant growth form as random intercept (five 501 levels: forbs, graminoids, pteridophytes, shrubs and trees). (B) Results of the variation partitioning analyses 502 representing the individual contribution of each environmental predictor. Bar plots are proportional to the 503 variation explained by the unique contribution of each fixed effect (expressed as R²<sub>m</sub>). In all graphs, 504 estimates and error bars represent the median value and 2.5 – 97.5 percentiles across 1,000 bootstrap 505 samples. See Fig. S7 for results on the analyses that also included rare species and Fig. S8 for a direct 506 comparison with the estimated effects of acidifying deposition.

507 **Figure 4 | Distance and geographic direction of colonization and extinction centroids.** (A) Distance 508 and geographic direction of colonization centroids (n species = 202). (B) Distance and geographic direction 509 of extinction centroids (n = 246). Longer distances reflect that colonization or local extinction are happening 510 in one preferred direction relative to the baseline centroid. Rayleigh's *r* statistic represents a test of 511 uniformity that compares the bearings of shifts to a uniform circular distribution (null hypothesis). Larger
values indicate more directed shifts. Asterisk (*) indicates significant deviations from the null hypothesis ($p < 0.05$). Results of the analyses that also included rare species are presented in Fig. S9. Results for each biogeographic region are presented in Fig. S10.

**Figure 5** | Effects of nitrogen deposition rates and species nitrogen niche width on colonization centroids in the most parsimonious models. Results of the mixed-effects model testing for the interaction effect between the average nitrogen deposition rate between the baseline survey and resurvey across each species' distribution (kg N ha⁻¹ year⁻¹) and species' ecological indicator values for nitrogen niche width (an index integrating the intra- and inter-regional variability in the nitrogen niche, with higher values indicating more generalist species and smaller values specialist species) on the south-north and west-east colonization centroids. Negative distances indicate southward (A) or westward colonization (B). Interaction effects between species nitrogen niche width and the average nitrogen deposition rate between the baseline survey and resurvey across each species' distribution are plotted along the color gradient. Model predictions were plotted to a maximum of 35 kg N ha⁻¹ year⁻¹ to avoid extrapolation uncertainty for deposition values where observations were scarce. Modelling results without outlier data ($n = 3$ data points) are provided in Fig. S13. The most parsimonious model of south-north colonization was an intercept-only model (blue solid line). Westward colonization depended on species nitrogen niche width. Nitrogen generalists that initially occurred in areas with lower rates of nitrogen deposition moved more westward. Colonization in the more specialist species was generally suppressed, irrespective of the nitrogen deposition rate. In all plots, the red dashed line represents the zero-line.
531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

Supplementary Materials for

Unexpected westward range shifts in European forest plants links to nitrogen deposition

Pieter Sanczuk¹*, Kris Verheyen¹, Jonathan Lenoir², Florian Zellweger³, Jonas J. Lembrechts⁴, Francisco

Rodriguez-Sanchez⁵, Lander Baeten¹, Markus Bernhardt-Römermann⁶,⁷, Karen De Pauw¹, Pieter

Vangansbeke¹, Michael P. Perring⁸,⁹, Imre Berki¹⁰, Anne Bjorkman¹¹,¹², Jörg Brunet¹³, Markéta

Chudomelová¹⁴, Emiel De Lombaerde¹, Guillaume Decocq², Thomas Dirnböck¹⁵, Tomasz Durak¹⁶,

Caroline Greiser¹⁷,¹⁸, Radim Hédl¹⁴,¹⁹, Thilo Heinken²⁰, Ute Jandt²¹,²², Bogdan Jaroszewicz²³, Martin

Kopecký²⁴,²⁵, Dries Landuyt¹, Martin Macek²⁴, František Máliš²⁶,²⁷, Tobias Naaf²⁸, Thomas A. Nagel²⁹,

Petr Petřík³⁰,³¹, Kamila Reczyńska³², Wolfgang Schmidt³³, Tibor Standovár³⁴, Ingmar Staude²²,³⁵,

Krzysztof Świerkosz³⁶, Balázs Teleki³⁷, Thomas Vanneste¹, Ondrej Vild¹⁴, Donald Waller³⁸, Pieter De

Frenne¹

*corresponding author : Pieter Sanczuk ; Pieter.Sanczuk@UGent.be

The PDF file includes:

Methods

Figs. S1 to S16

Data S1
553

**Methods**

554

1. Study area and vegetation data

555

556

557

558

559

560

561

562

563

564

565

566

567

We compiled a database including 2,954 permanent and semi-permanent vegetation plots distributed in mature forest stands across five biogeographic regions (Atlantic [n plots = 322], Alpine [266], Boreal [21], Continental [2049] and Pannonian [296]) (Fig. 1; www.forestREplot.ugent.be) (40, 41). Community composition data were derived from baseline surveys recorded between 1933 and 1994 and paired resurveys carried out after the baseline surveys between 1987 and 2017 (median [min – max] inter-survey interval: 39 [13 – 67] years). All plots were located in ancient (such that they have not been cleared for any other land use since at least the 18th century) mature forest stands and excluded heavily managed plantations. Vegetation plots did not experience major disturbances between the baseline surveys and the resurveys (e.g., no replanting or clearcutting). All plots were permanently marked (i.e., physical mark permanently visible in the field) or quasi-permanent. Plot sizes for botanical surveys ranged between 25 and 1300 m². The presence and abundance (as the percentage ground cover, visually estimated) of all vascular plant species in the understorey (defined as all plants < 1.3 meter tall; thus including seedlings of tree and shrub species), was recorded in all plots. Plant taxonomy was standardized with the R package *taxize* (42).

568

569

2. Centroid shifts

570

2.1. Rate (speed and direction) of range centroid shifts

571

572

573

574

575

576

577

578

579

580

Species centroid shifts are brought about by the combined effects of local colonization (a new arrival in the community of a given plot), changes in abundance, and local extinction, and these processes occur across the entire distribution of species (14, 43, 44). We quantified the rate of range shifts based on the actual change, over time, in the geographical coordinates of the range centroids, i.e., the abundance-weighted geometric center of a species' distribution across the vegetation plot network. In contrast to the more frequent quantification of range shifts from either the trailing or leading edges, analyzing centroid shifts allows us to obtain more robust estimates of the magnitude and geographic direction of complex distributional shifts (6, 11, 12, 21). This is important, since range shift estimates based on either leading or trailing edges can be affected by stochastic processes and low sample sizes that may blur the overall biogeographical trends (11).

581

582

Species' centroid shifts were quantified across all vegetation plots across the study area. To avoid bias attributed to rare species observed in only a few of the plots, we quantified range centroid shifts for all
species that occurred in more than 1% of the plots. Across the entire plot network, this resulted in a final set of 266 species. However, as a sensitivity analysis, range centroid shifts (and all downstream analyses) were also quantified for all species that occurred in more than 0.1% of the plots (N species = 596, for the entire plot network) (results in Figs. S3, S7). Centroid shifts for each separately biogeographic region (defined by the European Environment Agency; www.eea.europa.eu) were also quantified by running the analyses on the subset of vegetation plots within each region (results in Figs. S4, S9). It is important to note here that the overall trends in centroid shifts across the entire study area are not directly comparable to the analyses for each biogeographic region separately, since likely disjunct shifts in different regions are averaged in the overall trends across the entire plot network.

We quantified the geographic displacement of centroids over time (schematically presented in Fig. S1) by means of the geographic distance and bearing (a circular variable expressed in decimal degrees ranging from 0° to 360°) between species' centroid positions at the time of the baseline survey ($t_1$) and at the time of the resurvey ($t_2$) (12, 21). First, for each species, the mean latitude and longitude was calculated at $t_1$ and $t_2$ weighted by the percentage cover (abundance-weighted) to locate the range centroids at the time of both surveys. Next, the geographic distance (minimal geodesic distance) between the centroids at $t_1$ and $t_2$ was quantified with the R package *sf* (45). The bearing was quantified as the geographical direction of travel along a rhumb line (so-called 'loxodrome', true course) between the centroids at $t_1$ and $t_2$ using the R package *geosphere* (46).

The distance and bearing of biotic centroid shifts were used to derive three response variables:

The absolute rate of centroid shifts ($km \ year^{-1}$) was calculated as the geographic distance between the centroid at $t_1$ and $t_2$ (i.e. the length of the vector) divided by the average inter-survey interval in all plots where species $sp_i$ was recorded.

The projected south-north and west-east rate ($km \ north \ year^{-1}$ and $km \ east \ year^{-1}$) of centroid shifts was calculated as the rate of centroid shift multiplied by the cosine (south-north rate) or the sine (west-east rate) of the geographic direction of shift.

The directionality (i.e., the angular dispersion of the bearings) of centroid shifts was quantified based on a Rayleigh's test of uniformity. In the Rayleigh's test, bearings of shifts are compared to a uniform circular distribution (i.e., the null hypotheses reflecting random shifts) to calculate Rayleigh's *r* statistic. This statistic represents the directionality of a set of vectors, ranging from zero (anisotropic distribution, meaning
random directional movements) to one (isotropic distribution; all vectors are centered to one single direction).

## 2.2. Distance of colonization and extinction centroids

To mechanistically unravel the processes at play in driving range centroid shifts, we decomposed centroid shifts of each species into shifts attributed to the effects of colonization (here, when a species is observed in a plot at $t_2$ but not in $t_1$), and shifts attributed to the effects of local extinction (when a species is observed in a plot at $t_1$ but not in $t_2$). For this analysis, we introduce the concept of colonization centroid- that is the centroid of newly colonized plots by a species (abundance-weighted by the percentage cover at $t_2$) – and extinction centroid – that is the abundance-weighted centroid of plots where a species became locally extinct (abundance-weighted by the percentage cover at $t_1$) (Fig. S1). Species-specific colonization and extinction centroids are expressed as the projected distance in each geographic direction (*km north* and *km east*) relative to the species’ original centroid at $t_1$. The northward and eastward distances of the colonization and extinction centroids were calculated analogous to the rate of centroid shifts, but not accounting for the inter-survey interval. From a biogeographical point of view, longer distances reflect that colonization or local extinction events took place farther away from the baseline centroid, suggesting that these processes are happening in a preferred direction relative to the baseline centroid position (schematically explained in Fig. S1C). The directionality of the colonization and extinction centroids was quantified based on the Rayleigh’s *r* statistic.

To avoid bias attributed to rare events of species colonization or local extinction in only a few of the plots, we quantified both centroid types for all species that colonized or went extinct in more than 1% of the plots. Across the entire plot network, this resulted in a final set of 202 species (colonization centroids) and 246 species (extinction centroids). However, as a sensitivity analysis, both centroids (and all downstream analyses) were also quantified for all species that colonized or went locally extinct in more than 0.1% of the plots (n species = 527 and 542 for colonization and local extinction, respectively) (results in Figs. S3, S7).

## 2.3. Null model of random movements

To assess to what extent the centroid shifts are independent from noise attributed to, for instance, the spatial distribution of the vegetation plots, we applied a conservative null model approach as an additional sensitivity analysis. For each species, we calculated baseline and resurvey centroids from resampled (with replacement) plots in the vegetation data set. The baseline centroid was calculated based on a resample of
all plots where the species was observed in the baseline survey, with a sample size equal to the original data set in the baseline period. The resurvey centroid was calculated based on a sample of all plots where the species was observed in both the resurvey and baseline survey (thus, including also newly colonized plots and plots where the species went locally extinct), with a sample size equal to the original data set in the resurvey period. We thus assumed that the species could have moved randomly within its distribution (here, conservatively defined as all plots where the species was observed). The null model was iterated 1,000 times for each species that was observed in >1% of the plots. Based on this sensitivity analysis, we conclude that the observed rate of centroid shifts was higher and the main directions less dispersed than would be expected at random owing to, for instance, the spatial distribution of the vegetation plots (i.e., shifts towards each cardinal direction were 62-70% faster than expected based on randomized shifts within each species' distribution; Fig. S5).

### 3. Environmental change variables

The absolute and projected velocities of species' centroid shifts and colonization and extinction centroids were linked to (i) the velocity of climate change (ii) nitrogen deposition rate and (iii) the rate of forest canopy cover change.

#### 3.1. The velocity of climate change

Spatially explicit climate change velocities were estimated by climate analogue mapping (4, 22), referring to the spatial change of climatic conditions (integrating both temperature and precipitation) between the baseline surveys and the resurveys. Climate analogue mapping is a statistical approach that maps all areas that have similar climatic conditions relative to a certain location's climate. Climate analogue mapping allowed us to quantify changes in the geographical distribution of climate over time, while also considering the geographic variation in climatic conditions at the landscape scale due to e.g., topoclimate. In contrast to the frequent calculation of climate change velocities based on thermal gradients and isotherms alone (3, 24), climate analogue mapping allowed us to consider consolidated shifts of multiple bioclimatic variables, i.e. the maximum growing-season temperatures, mean minimum winter temperatures and total growing-season precipitation and represented the climatic norms over two ten-year periods (the baseline period [1958-1967] and resurvey period [2007-2016]; time interval of 49 years) (Fig. S2). Gridded monthly climate data were retrieved from the TerraClimate database (36) and analysed at 1/24th degree native resolution (approximately 4 km in the study area). We identified climate analogues for all 2,954 resurveyed
vegetation plots. For each plot, we performed a grid search at 4 km resolution to map all raster cells within the study area (N cells = 518,205) that had similar climatic conditions in the resurvey period compared to a plot's climate in the baseline period. Climate 'analogy' was assessed by testing, for each bioclimatic variable separately, whether the climatic variation of any location in the resurvey period could be drawn from the same distribution as the climate variation for a given plot in the baseline period using non-parametric two-sample Wilcoxon rank sum tests. We considered significant differences (*p*<.05) as novel climates and non-significant differences (*p*≥.05) as analogue climates. Climate analogues were finally assessed by overlaying the climate analogue areas for each bioclimatic variable separately (see example in Fig. S2).

The velocity of climate change was subsequently calculated as the geographic distance (in *km*) between the plot's coordinates and the cell coordinates of the nearest climate analogue, divided by the time interval. Quantifying the velocity of climate change based on the nearest climate analogue methods is particularly useful in the context of this study, as climate analogue shifts are conceptually identical to biotic centroid shifts: estimates of the absolute velocity of climate change (*km year<sup>-1</sup>), as well as the projected velocities along the south-north (*km north year<sup>-1</sup>*) and west-east (*km east year<sup>-1</sup>*) axes were obtained by using the exact formulas as the ones used to estimate the rate of centroid shifts. The directionality of climate analogue shifts was quantified based on the Rayleigh's *r* statistic.

The calculation of climate change velocities was performed in R making use of the packages *raster* (47) and *geosphere* (Hijmans 2019). Parallel computation was implemented using the R packages *foreach* (48) and *doParallel* (49). Maps were produced using the R package *tmap* (50).

Note that we here map climate analogues based on statistical testing as an alternative to the classification of climate analogue mapping based on continuous dissimilarity matrices because (1) it can provide information on the contribution of each bioclimatic variable separately; and (2) it allows us to consider the interannual climatic variability (ICV) in both the baseline and resurvey periods, which is biologically very relevant. For example, locations with stable climatic conditions (with low ICV) are mapped as novel climatic conditions even under small changes in the climate system, while areas with a highly variable climates (with large ICV) are more resilient and will not experience biologically relevant changes even under relatively large oscillations in the climate system.

## 3.2. Nitrogen and sulphur deposition rate

We quantified the nitrogen (N) and sulphur (S) deposition rates for each plot using the EMEP gridded database (https://emep.int/mscw/mscw_moddata.html), providing modelled nitrogen and sulphur
deposition data at 0.1° native resolution (~8 km × 8 km within the study area). For each plot in the vegetation network, the rate of nitrogen (oxidized and reduced) and sulphur (oxidized) deposition was quantified as the total wet and dry deposition (dry deposition accounted for deciduous temperate forest surface) between the baseline survey and the resurvey, divided by the inter-survey time interval (expressed in total *kg ha⁻¹ year⁻¹*). Values of annual nitrogen deposition before the year 1990 were obtained based on the reference year 2000 and corrected for the time period using the period-specific correction factors published in (51). Yearly values of nitrogen depositions from 1990 onwards were directly retrieved from the EMEP data bases. Sulphur and nitrogen deposition contribute to the acidifying deposition rate. This rate (calculated as: *kg N ha⁻¹ year⁻¹/14 + (kg S ha⁻¹ year⁻¹/32.06)×2*, Keq *ha⁻¹ year⁻¹*) was very strongly correlated (Spearman correlation: 0.87) to the nitrogen deposition rate across all 2,954 plots **due to partially shared emission sources**.

### 3.3. Rate of forest canopy cover change

For each plot in the study area, forest canopy structure was quantified *in situ* as the total cover of the tree canopies (>7 m) in the baseline survey and resurvey based on visual estimations. Species-specific tree cover estimates were summed in each plot. The rate of forest canopy cover change was quantified by subtracting the total canopy cover in the resurvey period from the total canopy cover in the baseline survey divided by the inter-survey time interval (*percentage canopy cover increase year⁻¹*).

### 3.4. Nitrogen **and** acidity (**reaction**) niche width

Species' nitrogen **and** acidity (**reaction**) niche width values were extracted from the EIVE-database (32), presenting ecological indicator values in European plant species (a numerical index between zero and ten; with higher values indicating more generalist species). The niche width metrics are continuous measures bounded between 0 and 10 that integrates the intra-regional and inter-regional variability in each species' nitrogen and **acidity** niche. In contrast to classical Ellenberg indicator values, which describe niche optima on ordinal scales, niche width data extracted from the EIVE-data base (32) have several important advantages for this study: (i) species with a more narrow nitrogen niche tend to be more specialized (33); (ii) the biogeographical implication is that species with a larger niche width often also have larger ranges (33); (iii) the continuous nature of the index facilitates its use within conventional linear mixed-effects modelling.

### 4. Statistical modeling
733 4.1. Linear mixed-effects models

734 We ran a set of linear mixed-effects models (LMM) to relate (i) the rate of centroid shifts (absolute, south-
735 north and west-east rate) and (ii) the northward and eastward distance of colonization and extinction
736 centroids to three global environmental change drivers that are expected to potentially alter range dynamics
737 in understorey species (1, 14, 16, 52): the velocity of climate change (the absolute velocity [km year⁻¹]), and
738 the projected south-north and west-east velocity [km north year⁻¹ and km east year⁻¹]), the nitrogen
739 deposition rate (kg N ha⁻¹ year⁻¹ of reduced and oxidized wet and dry deposition of nitrogen) and the rate of
740 forest canopy cover change (% canopy cover increase year⁻¹). The LMMs included all possible two-way
741 interactions and accounted for plant growth form (random intercept, including five levels: forbs,
742 graminoids, pteridophytes, shrubs and trees) as a random intercept and were conducted with the R package
743 lme4 (53). Because centroid shifts are quantified in a species-specific manner while predictor variables
744 describe the environment at the plot-level, we calculated species-specific abundance-weighted mean values
745 for each predictor to match the observational unit of the model (i.e., the average environmental change
746 values across all plots where species sp occurred, weighted for its original abundance (12); see Fig. S6 for
747 a data flow chart). The absolute velocity of biotic centroid shifts and the absolute velocity of climate
748 analogue shifts were square root transformed to obtain normality. Pairwise Spearman correlations (r) among
749 the environmental change predictors were acceptable (median |r| = 0.09; maximum |r| = 0.46 for nitrogen
750 deposition rate and the absolute velocity of climate change) (54). For each model, we selected the most
751 parsimonious (based on Akaike information criterion with small-sample correction [AICc]) using the R
752 package MuMln (55), with the restricted maximum likelihood argument to ‘FALSE’. Once the best model
753 structure was selected, we set the restricted maximum likelihood argument to ‘TRUE’ for exact coefficient
754 estimation (53). Model coefficient estimates (mean values) and 95% confidence intervals (2.5 and 97.5
755 percentiles) were generated by iterating the single best model structure on 1,000 bootstrap samples.
756 Significance was considered when the 95% confidence interval did not include zero. Model fit was assessed
757 as the percentage of variance explained by the fixed effects (marginal R²; R²ₘ) and the percentage of
758 variance explained by both fixed and random effects (conditional R²; R²ₑ) following (56). Prior to
759 modelling, all predictor variables were z-transformed to allow a better comparison of the predictors’ effect
760 sizes.

761 To test whether the colonization centroids were associated to the species nitrogen and acidity niche width
762 and the rates of atmospheric (nitrogen and acidifying, respectively) deposition (and their respective pairwise
763 interaction effects), we ran a LMM with plant growth form as a random intercept term. The most
764 parsimonious model was selected as described above.
765 4.2. Variation partitioning

766 In a final step, the unique contribution of each of the three environmental change variables to the full model
767 was analysed in a variation partitioning analysis. To obtain the variation explained uniquely by each focal
768 environmental change variable, we subtracted the variation explained by the fixed effects (marginal R²;
769 R²_m) of the partial model, i.e., the full model minus one (out of the three) focal environmental change
770 predictor variables, from the variation explained by the fixed effects in the full model. The procedure was
771 repeated on 1,000 bootstrap replicates and the average proportion variation explained (average R²_m) and
772 95% confidence intervals were calculated. All analyses were performed in R version 4.2.2 (57)
<figure><img src="image_2.png" /><figcaption>**Figure S1** | Schematic illustration of centroid shifts and climate analogue shifts. (A) Schematic overview and calculation of the rate of centroid shifts. (B) Schematic overview and calculation of the velocity of climate change based on climate analogue shifts. (C) Conceptual framework and definition of centroid shifts, colonization and extinction centroids for a certain species (sp.i). From a biogeographical point of view, longer distances reflect that colonization or local extinction events took place farther away</figcaption></figure>

774

775

776

777

778

779
780 from the baseline centroid, suggesting that these processes are happening in a preferred direction relative

781 to the baseline centroid position.
<figure><img src="image_3.png" /><figcaption>**Figure S2** | Climate analogue mapping integrating three bioclimate variables. Climate analogue mapping at 4 km × 4 km resolution for the climatic conditions in the baseline (average climate values in the period 1958-1967) and resurvey period (average climate values in the period 2007-2016) in two vegetation plots. Climate analogues for a given plot (black square) in the baseline period (left) and shifts in climate analogues in the resurvey relative to the climatic conditions in the baseline period (right) for three bioclimate variables separately: maximum growing-season temperatures (orange), minimum winter temperatures (blue) and total growing-season precipitation (grey), and the overlayed climate analogue areas of all three bioclimatic variable together (red). The velocity of climate change was calculated as the distance to the nearest climate analogue raster cell (thus, within the red area) over time.</figcaption></figure>

782

783

784

785

786

787

788

789

790

791
A
Species centroid shifts

<figure><img src="image_4.png" /></figure>

B
Climate analogue shifts

<figure><img src="image_5.png" /></figure>

**Figure S3** | Rate and geographic direction of centroid shifts and climate analogue shifts. (A) Rate and geographic direction of species centroid (including rare species; *n* species = 596). (B) Rate and geographic direction climate analogue shifts (*n* plots = 2,954). In all graphs, the Rayleigh's *r* statistic represents a test of uniformity that compares the bearings of shifts to a uniform circular distribution (null hypothesis). Larger values indicate more directional shifts. Asterisk (*) indicates significant deviations from the null hypothesis (*p* < 0.05).

792

793

794

795

796

797

798
<figure><img src="image_6.png" /><figcaption>**Figure S4** | Rate and geographic direction of centroid shifts and climate analogue shifts per biogeographic region. (A) Rate and bearing of species centroid shifts (n species Alpine = 173; Atlantic = 144; Boreal = 41; Continental = 225; Pannonian = 256). (B) Rate and bearing of climate analogue shifts (n plots Atlantic = 322; Alpine = 266; Boreal = 21; Continental = 2049; Pannonian = 296). Rayleigh's *r* statistic represents a test of uniformity that compares the bearings of shifts to a uniform circular distribution (null hypothesis). Larger values indicate more directed shifts. Asterisk (*) indicates significant deviations from the null hypothesis (*p* < 0.05). Biogeographical boundaries were defined by the European Environment Agency (www.eea.europa.eu).</figcaption></figure>

799

800

801

802

803

804

805

806

807
<figure><img src="image_7.png" /><figcaption>**Figure S5** | Rate and geographical direction of centroid shifts that can be attributed to random noise due to the plot distribution. (A) Rate and geographic direction of species centroids for 1,000 bootstrapped replicates (n species = 266 × 1,000 replicates). (B) Rate and geographic direction of null models where species-specific abundances were randomized across all plots in which the species was observed, and resampled 1,000 times (n species = 266 × 1,000 replicates). In all graphs, the Rayleigh's *r* statistic represents a test of uniformity that compares the bearings of shifts to a uniform circular distribution (null hypothesis). Larger values indicate more directional shifts. Asterisk (*) indicates significant deviations from the null hypothesis (*p* < 0.05). (C) Comparison between species-specific bootstrapped estimates and null models of the absolute rate and the south-north and west-east rate (grey symbols; average across 1,000 replicates). The magnitudes of centroid shifts are significantly different from random noise owing to the spatial distribution of the vegetation plots by 66% (northward shifts were 66% faster than expected based on random movements), 62% (southward shifts), 63% (westward shifts) and 70% (eastward shifts).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure S6** | Framework of the study and flow chart of the data. Species abundance data were collected at the plot-level during a baseline and resurvey period. Environmental change data were collected at the plot-level. Species-specific abundance-weighted centroids were calculated for the baseline survey and the resurvey. The rate of centroid shift was calculated as the absolute rate, the projected south-north and west-east rate, and the northward and eastward distance of colonization and extinction centroids. Species-specific environmental change data were obtained by taking the average environmental change values (realized change between the baseline and resurvey period) in all plots where the species was observed in the baseline period, weighted by its original abundance. Using linear mixed-effects models, biotic centroid shifts were linked to the abundance-weighted environmental change data and all pairwise interactions (indicated by ²). Plant growth form (with five levels: forbs, graminoids, pteridophytes, shrubs and trees) was included as a random effect term in all models.</figcaption></figure>

821

822

823

824

825

826

827

828

829

830

831

832
<figure><img src="image_9.png" /><figcaption>**Figure S7** | Effects of environmental changes on centroid shifts based on the most parsimonious model, including rare species. (A) Results of the mixed effect models (including rare species; *n* species = 596) indicating coefficient estimates and 95% confidence intervals (CI) of the effects of velocity of climate change (*km year*⁻¹, *km north year*⁻¹, *km east year*⁻¹), average nitrogen deposition rate between the baseline survey and resurvey (*kg N ha*⁻¹ *year*⁻¹), and canopy change (% cover increase *year*⁻¹), and all pairwise interactions (indicated by ‘×’) on the absolute rate of centroid shifts (*km north year*⁻¹; negative values indicate southward shifts) and west-east rate (*km east year*⁻¹; negative values indicate westward shifts) extracted from the most parsimonious model (empty rows were not included in the final selected model). All the predictor variables are z-transformed to increase comparability. Model fit is presented as the proportional explained variation by the fixed effect (marginal R², R²<sub>m</sub>) and the proportional explained variation by the fixed and random effects (conditional R², R²<sub>c</sub>). Models accounted for plant growth form as random effect (five levels: forbs, graminoids, pteridophytes, shrubs and trees); (B) Results of the variation partitioning analyses representing the individual contribution of each environmental predictor. Bar plots are proportional to the explained variation by the unique contribution of each fixed effect (expressed as R²<sub>m</sub>). In all graphs, estimates and error bars represent the median value and 2.5 – 97.5 percentiles across 1,000 bootstrap samples.</figcaption></figure>
A

<figure><img src="image_10.png" /></figure>

B

<figure><img src="image_11.png" /></figure>

C

<figure><img src="image_12.png" /><figcaption>**Figure S8** | Effects of environmental changes on the rate of centroid shifts (A), distance of colonization (B) and extinction (C) centroids based on the most parsimonious model. (A) Results of the mixed-effects models (*n* species = 266) indicating coefficient estimates and 95% confidence intervals (CI) of the effects of velocity of climate change (km north year⁻¹, km east year⁻¹), average acidifying deposition rate between the baseline and resurvey (Keq ha⁻¹ year⁻¹), and canopy change (% cover increase year⁻¹), and all pairwise interactions (indicated by ‘×’) on the absolute rate, south-north rate (km north year⁻¹; negative values indicate southward shifts) and west-east rate (km east year⁻¹; negative values indicate westward shifts) of centroids shifts extracted from the most parsimonious model (empty rows were not included in the final selected model). (B, C) Results of the mixed-effects models (*n* species = 202 and 246 for colonization and extinction centroids, respectively) indicating coefficient estimates (95% CI) of the environmental changes on south-north and west-east distance of colonization extinction centroids extracted from the most parsimonious model. In all graphs, estimates and error bars represent the median value and 2.5–97.5 percentiles across 1,000 bootstrap samples. Bar plots are proportional to the variation explained by the unique contribution of each fixed effect (expressed as R²). Model fit is presented as the proportion variation explained by the fixed effect (marginal R², R²ₘ) and the proportion variation explained by the</figcaption></figure>
866 fixed and random effects (conditional R², R²c). Models accounted for plant growth form as random effect

867 (five levels: forbs, graminoids, pteridophytes, shrubs and trees). All the predictor variables were z-

868 transformed to increase comparability.
A

<figure><img src="image_13.png" /></figure>

B

<figure><img src="image_14.png" /></figure>

Figure S9 | Distance and geographic direction of colonization and extinction centroids accounted for rare species. (A) Distance and geographic direction of colonization centroids (including rare species; *n* = 542). (B) Distance and geographic direction of extinction centroids (*n* = 527). Longer distances reflect that colonization or extinction took place further away from the original centroid and that these processes are happening in a preferred direction. Rayleigh's *r* statistic represents a test of uniformity that compares the bearings of shifts to a uniform circular distribution (null hypothesis). Larger values indicate more directed shifts. Asterisk (*) indicates significant deviations from the null hypothesis (*p* < 0.05).
A

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

distance from
centroid (km)

<figure><img src="image_20.png" /></figure>

B

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

Figure S10 | Distance and geographic direction of colonization and extinction centroids for each biogeographic region. (A) Distance and bearing of colonization centroids (n species Alpine = 130; Atlantic = 126; Boreal = 7; Continental = 193; Pannonian = 212). (B) Distance and bearing of extinction centroids (n species Alpine = 140; Atlantic = 125; Boreal = 28; Continental = 211; Pannonian = 250). Rayleigh's *r* statistic represents a test of uniformity that compares the bearings of shifts to a uniform circular distribution (null hypothesis). Larger values indicate more directed shifts. Asterisk (*) indicates significant deviations from the null hypothesis (at the level of p < 0.05). Biogeographical boundaries were defined by the European Environment Agency (www.eea.europa.eu).
<figure><img src="image_27.png" /><figcaption>**Figure S11** | Effects of environmental changes on the distance of colonization and extinction centroids based on the most parsimonious model. (A-B) Results of the mixed-effects models (*n* species = 202) indicating coefficient estimates and 95% confidence intervals (CI) of the effects of velocity of climate change (*km north year*⁻¹, *km east year*⁻¹), average nitrogen deposition rate between the baseline and resurvey (*kg N ha*⁻¹*year*⁻¹), and canopy change (% cover *increase year*⁻¹), and all pairwise interactions (indicated by ‘×’) on northward (*km north*; negative values indicate southward colonization) and eastward distance (*km east*; negative values indicate westward colonization) of colonization centroids extracted from the most parsimonious model (empty rows were not included in the final selected model), and output of the variation partitioning analysis representing the individual contribution of each environmental predictor. The results show that westward colonization was more common in species that have experienced a lower rate of nitrogen deposition across their distribution. (C-D) Results of the mixed-effects models (*n* species = 246) indicating coefficient estimates (95% CI) of the effects of velocity of climate change, average nitrogen</figcaption></figure>
899 deposition rate between the baseline and resurvey, and canopy change, and all pairwise interactions on
900 northward and eastward distance of extinction centroids extracted from the most parsimonious model, and
901 output of the variation partitioning analysis. The results show that species' southward extinction centroids
902 were related to the poleward velocity of climate change and high rates of nitrogen deposition. Eastward
903 extinction centroids were related to a high rate nitrogen deposition and canopy opening. Eastward extinction
904 owing to climate change was more common under canopy opening. In all graphs, estimates and error bars
905 represent the median value and 2.5 - 97.5 percentiles across 1,000 bootstrap samples. Bar plots are
906 proportion to the variation explained by the unique contribution of each fixed effect (expressed as R²). Model fit is presented as the proportional explained variation by the fixed effect (marginal R², R²) and the
907 proportion variation explained by the fixed and random effects (conditional R², R²). Models accounted for
908 plant growth form as random effect (five levels: forbs, graminoids, pteridophytes, shrubs and trees). All the
909 predictor variables were z-transformed to increase comparability. See Fig. S12 for results on the analyses
910 that also included rare species.
911
<figure><img src="image_28.png" /><figcaption>**Figure 12** | Effects of environmental changes on the distance of colonization and extinction centroids based on the most parsimonious model, including rare species. (A-B) Results of the mixed-effect models (including rare species; *n* species = 542) indicating coefficient estimates and 95% confidence intervals (CI) of the effects of velocity of climate change (*km north year⁻¹*, *km east year⁻¹*), average nitrogen deposition rate between the baseline and resurvey (*kg N ha⁻¹ year⁻¹*), and canopy change (% cover increase year⁻¹), and all pairwise interactions (indicated by ‘×’) on northward (*km north*; negative values indicate southward colonization) and eastward distance (*km east*; negative values indicate westward colonization) of colonization centroids extracted from the most parsimonious model (empty rows were not included in the final selected model), and output of the variation partitioning analysis representing the individual contribution of each environmental predictor. (C-D) Results of the mixed-effects models (including rare species; *n* species = 527) indicating coefficient estimates (95% CI) of the effects of velocity of climate change, average nitrogen deposition rate between the baseline and resurvey, and canopy change, and all</figcaption></figure>
pairwise interactions on northward and eastward distance of extinction centroids extracted from the most parsimonious model (empty rows were not included in the final selected model), and output of the variation partitioning analysis. In all graphs, estimates and error bars represent the median value and 2.5 – 97.5 percentiles across 1,000 bootstrap samples. Bar plots are proportional to the explained variation by the unique contribution of each fixed effect (expressed as R²). Model fit is presented as the proportion of explained variation by the fixed effect (marginal R², R²ₘ) and the proportion of explained variation by the fixed and random effects (conditional R², R²c). Models accounted for plant growth form as random effect (five levels: forbs, graminoids, pteridophytes, shrubs and trees); All the predictor variables are z-transformed to increase comparability.
<figure><img src="image_29.png" /><figcaption>Figure S13 | Spatial variation in the nitrogen deposition rate across the baseline centroid positions.</figcaption></figure>

(A) Map of modelled nitrogen (N) deposition rate (sum oxidized and reduced wet and dry deposition expressed in kg N ha⁻¹ year⁻¹; dry deposition accounted for deciduous forest surface) at 0.1° resolution for the reference year 2000. Black dots indicate the abundance-weighted centroid of the 266 most common species in the data set. Grey crosses indicate the distribution of the 2,954 vegetation plots (spatially jittered for clarity). (B) species-specific abundance-weighted nitrogen deposition rates regressed to the degree (°) latitude and longitude of each species' baseline centroid (n species = 266). The regression line represents the model fit of a generalized additive model, with the k-parameter restricted to three to avoid overfitting.

934

935

936

937

938

939

940

941

942
<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>Figure S14 | Effects of nitrogen deposition rates and species nitrogen niche width on colonization centroids based on the most parsimonious model when outlier data points were excluded (n = 3 data points). Results of the mixed-effects model testing for the interaction effect between the average nitrogen deposition rate across each species' distribution (kg N ha⁻¹ year⁻¹) and species Ecological Indicator Value for nitrogen niche width (an index integrating the intra- and inter-regional variability in the nitrogen niche, with higher values indicating more generalist species) on the south-north and west-east colonization centroids. Negative distances indicate southward (A) or westward colonization (B). The effects of nitrogen niche width is plotted with the color gradient. The most parsimonious model structure for west-east colonization did not include the interaction effect anymore, but colonization centroids for generalist species were still often more westward. Colonization centroids of specialist species were either west- or eastward, depending on the experienced nitrogen deposition rate. The most parsimonious model of south-north colonization was an intercept-only model (blue line). Nitrogen generalist species that initially occurred in areas with lower rates of nitrogen deposition moved westward. Colonization in the more specialist species was equally likely westward or eastward, depending on the average rate of nitrogen deposition across their distribution.</figcaption></figure>
<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

Figure S15 | Model output of the effects of acidifying deposition rates (Keq ha⁻¹ year⁻¹) and species acidity (reaction) niche width on colonization centroids based on the most parsimonious model.

Results of the mixed-effects model testing for the interaction effect between the average acidifying deposition rate across each species' distribution (Keq ha⁻¹ year⁻¹) and species Ecological Indicator Value for acidity (reaction) niche width (an index integrating the intra- and inter-regional variability in the reaction niche) on the south-north and west-east colonization centroids. Negative distances indicate southward (A) or westward colonization (B). The acidity niche width was not included in the most parsimonious model structure. Colonization centroids of species were mainly westward for species that experienced lower rates of acidifying deposition rates. The most parsimonious model of south-north colonization was an intercept-only model (blue line).

959
960
961
962
963
964
965
966
967
968
969
<figure><img src="image_34.png" /><figcaption>**Figure S16** | Distribution of the top-soil pH (pH-H₂O) conditions across the studied vegetation plots.

Orange: interpolated top-soil pH (pH-H₂O) conditions by overlaying all out 2,954 studied vegetation plots on top of the gridded soil data layer (SoilGrids.org, 250-meter resolution) available at 250-meter resolution globally. Blue: measured top-soil pH (pH-H₂O) conditions from in-situ soil samples available for 704 vegetation plots across the study area. In-situ soil pH data is available from (58). The median value of top soil pH (pH-H₂O) levels across the plots was between 5.3 (data SoilGrids.org, n plots = 2,954) 5.04 (in-situ soil samples, n plots = 704). The soils across the majority of our plots are relatively well buffered against acidifying deposition due to base cation exchange. Only 0.6 % of all plots (8.8% of the subset of plots with field data) have soil pH levels below the critical pH level of 4.2, a threshold below in which base cations become depleted and levels of Al³⁺ increase in the soil solution (59).</figcaption></figure>
981 **Supplementary Data**

982 **Data S1 | Species list and rates of centroid shifts and environmental changes.** Full list of study species
983 considered, and species-specific values of the absolute rate (km year⁻¹), the south-north rate (km north year⁻¹) and the west-east rate (km east year⁻¹) of centroid shifts; and species-specific experienced rates of
984 atmospheric (nitrogen and acidifying) deposition (kg N year⁻¹ ha⁻¹; K eq year⁻¹ ha⁻¹), forest canopy cover
985 change (% cover increase yr⁻¹) and climate change velocities (km yr⁻¹, km north yr⁻¹, km east yr⁻¹).