Climate change reshuffles northern species within their niches

Climate change is a pervasive threat to biodiversity. While range shifts are a known consequence of climate warming contributing to regional community change, less is known about how species’ positions shift within their climatic niches. Furthermore, whether the relative importance of different climatic variables prompting such shifts varies with changing climate remains unclear. Here we analysed four decades of data for 1,478 species of birds, mammals, butterflies, moths, plants and phytoplankton along a 1,200 km high latitudinal gradient. The relative importance of climatic drivers varied non-uniformly with progressing climate change. While species turnover among decades was limited, the relative position of species within their climatic niche shifted substantially. A greater proportion of species responded to climatic change at higher latitudes, where changes were stronger. These diverging climate imprints restructure a full biome, making it difficult to generalize biodiversity responses and raising concerns about ecosystem integrity in the face of accelerating climate change. The authors analyse four decades of distribution data for various taxonomic groups to understand the shift of species within their climatic niches and the changing influences of different climate factors. The diverse and diverging climate imprints raise concerns about future ecosystem integrity.

C limate change impacts are projected to increase, potentially matching land-use change as the main threat to biodiversity [1][2][3] . Current efforts aimed at predicting the consequences of climate change tend to focus on range shifts and on community change following local colonization and extinction of species 4,5 . These predictions are explicitly or implicitly based on the niche concept-that is, species typically tolerate a restricted range of environmental conditions 6 , defining the 'niche space' so that if the environment changes, species either adapt, move or eventually decline and go extinct. Indeed, climate change has prompted extensive shifts in species abundance and distributions 7-10 , potentially leading to increased homogenization of community composition across regions 11 . Less explored is the relative importance of species shifting their position within their climatic niche space over time. Furthermore, as both climate change and niche space involve multiple dimensions, the relative importance of different climatic variables may vary among taxa, across space and over time [12][13][14][15] . Consequently, asymmetric changes in climatic conditions may result in complex responses among species within regional communities. This heterogeneity in both species responses and climatic dynamics hampers our ability to assess the consequences of climate change for the structuring of communities across large spatial scales.
Quantifying how communities respond to ongoing climate change requires long-term spatio-temporally resolved data on phylogenetically diverse arrays of species. Such data are scarce, particularly for northern regions, where climate change is most pronounced 16 . Here we take advantage of a unique dataset to quantify speciesʼ responses to multiple climatic variables over four decades  across an ~1,200 km latitudinal gradient in Finland (Fig. 1). Using occurrence data on 1,478 species of birds, mammals, small rodents, butterflies, moths, forest understory vascular plants and freshwater phytoplankton, we analyse shifts in species' relative niche position over time-that is, whether species occur at the lower end, at the optimum or at the upper end of their climatic niche ( Fig. 2; each dataset was collected and treated independently). As niche dimensions, we explore annual values of mean temperature, total precipitation, duration of snow cover and the North Atlantic Oscillation (NAO) index-a proxy for regional climate variation known to affect ecological processes at various scales 17 .
We first quantified whether these climatic variables have changed at different rates within the latitudinal gradient. We partitioned the study area into three bioclimatic zones, whereby we expect stronger changes in the northernmost zone, specifically regarding temperature and snow cover. Then, to assess how these climatic variables affect species occurrence patterns and whether their relative importance has changed over time or across latitude, we fitted joint species distribution models 18 for each combination of taxonomic group × bioclimatic zone × decade ( Fig. 1 and Methods). Given the stronger changes expected towards the pole 16 , we expect that more species respond to climatic change in the northernmost zone. Furthermore, as most species have geographical ranges expanding further south than the study area, we expect a higher prevalence of responses within the lower and optimum areas of niche space at higher latitudes as climate change progresses (Fig. 2c).
We found strong climatic changes during the study period with uneven rates of change along the latitudinal gradient (Fig. 1c, Extended Data Fig. 1, Supplementary Table S1 and Supplementary Results). Specifically, the northernmost zone experienced stronger increases in temperature and precipitation compared with the middle and southern zones, while the duration of snow cover showed stronger declines both in the northernmost and southernmost zones (Supplementary Table S1). In later decades, average conditions in more northern areas resembled those of earlier decades in middle and southern areas; this pattern was particularly apparent for temperature and snow cover ( Fig. 1c and Extended Data Fig. 1).
The relative position of species within their niche shifted substantially over time (Fig. 3). Across decades, a large proportion of species shifted position between the lower end of niche space (where an increase in the climatic covariate has a positive impact on species occurrence), the niche optimum (that is the bell-shaped area of the curve) or the upper end (where an increase in the climatic covariate has a negative impact on species occurrence; Fig. 3 and Extended Data Fig. 2). Despite there being turnover in species composition within bioclimatic zones (Extended Data Fig. 3), the proportion of species unique to each decade in each region was low, with few exceptions, even when comparing the first and last decades for each taxa (Supplementary Figs. S1 and S2). Furthermore, shifts in occurrence patterns did not lead to strong directional homogenization of communities' composition across space (Extended Data Fig. 3). Specifically, within a given bioclimatic zone, assemblages either remained equally homogeneous in space over decades (for example, mammals), became slightly more (for example, plants) or less (for example, birds, butterflies and moths) homogeneous, or varied in complex ways (for example, small rodents; Extended Data Fig. 3 and Supplementary Fig. S3). Nonetheless, the rate and direction of species replacement over decades differed across zones-the extent of turnover for bird, butterfly and moth assemblages was stronger poleward, while the opposite occurred for phytoplankton. Together, these results strongly suggest that the observed shifts in species responses result from species moving to new domains within their environmental niche.
Across taxonomic groups, the proportion of species responding to the climatic variables increased towards the north but with considerable variation in the direction of those responses, both among climatic variables and among taxa ( Fig. 3 and Extended Data Fig. 2). Temperature emerged as the strongest predictor of variation in species occurrences across all regions, decades and taxa (Supplementary  Table S2). Changes in species' responses to temperature over time were particularly prevalent in the north, mainly occurring towards the lower part of niche space (Fig. 3a). Specifically, over the decades, a growing proportion of bird, mammal and phytoplankton species responded positively to temperature in the northernmost zone (Fig. 3b). In contrast, proportionally more butterfly species occurred at the optimum or upper end of niche space in later decades in the southernmost zone, potentially suggesting that current conditions are getting too warm for these taxa.
Along other niche dimensions, the proportion of species responding to precipitation consistently declined from earlier to later decades, mostly due to fewer species occupying the bell-shaped, optimum part of their niche. In contrast, the proportion of species responding to NAO consistently increased over time (Fig. 3a); this was mostly driven by species occurrences being negatively affected by NAO, particularly bird and moth species (Extended Data Fig. 2). The duration of snow cover emerged as the main predictor of species distributions in one-third of our models and was particularly relevant in explaining plant occurrences (Supplementary  Table S2). Nonetheless, responses to snow cover were more complex, with strong variation across decades and taxa ( Fig. 3a and Extended Data Fig. 2).
The relative importance of the climatic variables as predictors of species occurrences changed markedly through time and differed across zones and taxa. For plants and mammals, all variables tended to explain an increasing proportion of variation in species occurrences over time and poleward for plants, while the opposite occurred for phytoplankton and moth species (Fig. 4). For birds, the explanatory power associated with precipitation and particularly with temperature also increased over time and with latitude ( Fig. 4).

Discussion
Our results underline major reshuffling across the overall biome along a 1,200 km high latitudinal gradient following progressively altered niche space of species. The observed community change is mostly attributable to shifts in climatic conditions relative to species' niches. Quantifying these relative shifts is key to understanding how species respond to ongoing and future climate change, as organisms may be more sensitive at the margins of their niche space and due to different species having different 'starting points' within niche space in relation to a climatic driver 12,13 . Despite complex responses to climate change among taxa, our analysis revealed a clear signal of stronger change towards the pole. Variation in exposure and sensitivity to environmental change is expected to prompt asymmetric trends of change across species and over space and time (for example, varying range or phenological shifts 8,10 ). This is consistent with our findings, as the magnitude and direction of climate-related responses were highly dependent on zone and climatic variables analysed. We further show that the relative importance of the different climatic variables in explaining species occurrences varied across space, time and taxa. This raises concerns about potentially misestimating impacts of climate change on biodiversity, given these may be too context-dependent to allow direct extrapolation, even between consecutive decades.
With accelerating rates of climate change, high latitude communities are thought to become hotspots of change (for example, refs. 19,20 ). Our results are consistent with this expectation in that a growing number of species responded to climatic change. These responses were particularly accentuated at higher latitudes, as northern communities are disproportionately exposed to faster climatic change (for example, ref. 21 ). Additionally, species richness can increase due to an influx of species from lower latitudes, while polar species might lose their habitats altogether 7,8,22 . Nonetheless, only a relatively small fraction of species was unique between decades, even between the first and last decade sampled. This implies that while turnover is a key aspect of ongoing community change, the observed shifts result from species moving to new domains within their niche space, regardless of their 'resident' or 'incoming' status. While some taxa showed stronger rates of turnover in the northernmost zone (birds and moths, in particular), the direction of such trends varied considerably between zones and taxa (see also refs. 23,24 ). Such asynchronous responses are consistent with species emerging as winners or losers, as climate change can either reduce or impose constraints on species fitness, abundance and distributions 8,[11][12][13]25 . Together, these changes translate into altered community structure.
Our findings suggest that with progressing climate change, a growing proportion of species is being reshuffled along their niche space and may increasingly be exposed to more or to less suitable climatic conditions. Compounding the complex dynamics of change for different climatic variables with the expected increase in climate variability 16 is likely to exacerbate species exposure to novel conditions. Specifically, many species appeared to be experiencing a 'thermal release' 13,26 in our analysis, with many responding positively to increasing temperature in later, warmer decades-particularly at higher latitudes. The frequency of species occurring within niche optimum conditions also increased over time. This increase in 'positive' and 'optimum' responses is consistent with species experiencing more favourable conditions and/or with newly suitable areas becoming available to species moving poleward, for instance leading to warm-affiliated species replacing cold-affiliated ones 21,25,27 . Such contrasting responses to climate change among co-occurring assemblages may have far-reaching implications for species interactions and ecosystem functioning 28 .
As the pace of climate change continues to accelerate, there is an urgent need to identify which species and communities may face higher risk of disruption and what the key drivers of community-wide disruption may be 5 . Shifting signatures of climatic variables over time were evident as changes in the proportion of variation explained by each niche dimension. Notably, the imprint of temperature mostly increased over time and did so for taxonomic     groups spanning a range of dispersal abilities, and against a background of complex changes in the other variables. While multiple dimensions underpin climate change 16 , the overriding imprint of warming points to temperature as a reliable predictor of biodiversity change. Indeed, temperature plays a fundamental role determining species metabolism, phenology and distributions, and has been pivotal in biodiversity change assessments and scenarios 5,9,10,12,13,29 . Our taxonomically comprehensive assessment of large-scale climatic drivers adds further support to this role. Beyond temperature, our results highlight the importance of including other variables in analyses of biodiversity change in response to climate change. Despite its importance for high latitude ecosystems, snow cover duration has been largely neglected in climate modelling 30 . The increasing importance of snow cover in explaining plant species occurrences in our study reinforces the need to quantify its effects more systematically, given the potential feedbacks between changes in vegetation at high latitudes and global climate 16,30 . Additionally, the growing proportions of species negatively affected by NAO may point to an increasing role of large-scale regional climatic impacts over time, although this interpretation is challenged by NAO's lower explanatory power compared with local niche dimensions.
We used a state-of-the-art joint species distribution modelling approach to analyse high-resolution community and environmental time series. By explicitly accounting for the multi-variate nature of species assemblages, our approach allows quantification of community-level patterns in how species respond to their environment over space and time, even with sparse data 18 . Yet, our ability to explain species occurrences varied between zones, decades and taxonomic groups. This clearly highlights the challenges in quantifying and forecasting composition of ecological communities 31 . While the evidence for climate change effects accumulates (for example, refs. 7-10,12,15 ), estimated and predicted effects remain highly uncertain among taxa, regions and climate change scenarios [1][2][3][4][5] . Failing to account for the spatio-temporal context of these climatic responses risks strongly misestimating current and future effects of climate change. In addition, future work is needed to determine whether species' niches have expanded, contracted or shifted altogether 32 , and how such responses may eventually translate into changes in abundance 15 or potentially affect species capacity to respond via plasticity or adaptively to climate change. Finally, the observed changes may partly reflect the effects of land-use change 15,24,33 , which may have more direct impacts on species occurrences and community composition. Our approach uses spatial latent variables and thereby accounts for the spatial autocorrelation that may arise from, for example, environmental covariates left unmeasured 18 . Yet, further explicit analyses are needed to understand the combined effects and potential interactions between climate and land-use changes (for example, ref. 3 ). Major challenges thus remain for models and scenarios of biodiversity change in accommodating the potential effects of shifting exposure and relative importance of different climatic variables. Our cross-taxa analysis quantifying detailed species responses across an entire biome provides a first critical step towards understanding variation in biodiversity change, frequently characterized by contrasting trends between metrics, taxa and regions 9,10,23,34,35 . Community changes in our study were mainly driven by species being pushed towards or away from their environmental optima, rather than by major upheaval of regional species pools. Recent reports have identified turnover of species identities as the main signal of global biodiversity change 35 . Our results highlight that even in the absence of strong region-wide colonization or extirpation events, the main signature of climatic change emerged from species being shifted to new domains within their environmental niche. We further revealed that the relative imprints of climatic variables shifted non-uniformly over time and across taxa. This asymmetric restructuring of co-occurring assemblages points to urgent concerns for both species persistence and ecosystems' integrity, as contrasting responses may result in disrupted species interactions and trophic links 28 . Such disruptions appear particularly likely for high latitude biomes, where greater proportions of species responded to climatic change.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-022-01381-x. in response to climate change? Range shifts and community effects in the borderland between forest and tundra.

Methods
Species Data. We analyse long-term high-quality monitoring data on 1,478 species of birds, mammals, small rodents, butterflies, moths, forest understory vascular plants and freshwater phytoplankton sampled across 6,504 individual sites along an ~1,200 km latitudinal gradient in Finland. Because of differences in sampling methods and in spatial and temporal coverage, each dataset was analysed separately. We note that most species have a distributional range larger than Finland and that for the current purposes, niche space was estimated based on empirical data compiled within Finland alone.
Birds. Bird data have been collected using line transect censuses in Finland since the 1970s 21 . The data are collected yearly based on a one-visit census in which birds are counted along transects with lengths typically 3-6 km. Transects are previously established (that is, with known locations) and not all transects are sampled every year. The census period is June, with observations typically carried between 3:00 and 9:00 am, when the singing activity of birds is highest in dry weather conditions. The observer walks alone at a speed of approximately 1 km h −1 depending on the density of birds along the transect using a map, compass or  36 . The data are curated by the Finnish Museum of Natural History. We used records between 1978 and 2017, including a total of 189 species sampled in 1,105 transects after applying our selection criteria ('Study design and data preparation' below).

Mammals. A systematic monitoring programme of counts of mammal snow tracks was established in 1989 by the Natural Resources Institute Finland (Luke;
Game triangle data 37,38 ). The wildlife triangle scheme is based on a network of 4 + 4 + 4 km triangle-shaped transects (totalling 12 km per triangle) with fixed locations, covering the entire country. The triangles are located in forested areas covering the main forest types and are usually situated in hunting areas with the observations carried out by volunteers (mainly hunters). Around 2,000 triangles have been established and about half of these are counted annually. In the winter count, the transect is walked or skied during one day and all snow tracks of 24 mammal species crossing the count line are recorded, usually from mid-January to mid-March (when snow cover conditions are good). The number of crossings are typically related to the transect length and number of days since last snowfall, when snow tracks have been accumulating. Snowfall can be replaced with a pre-count, where the existing tracks are marked or erased, to be disregarded during the actual count. We used records between 1989 and 2017, including a total of 18 species sampled in 1,958 sites after applying our selection criteria ('Study design and data preparation' below).

Small rodents.
Since the 1960s, the Natural Resources Institute Finland (Luke) carries out an inventory of vole species to support forest management planning (small rodents data 39 ). Data were collected biannually in spring (mid-March to mid-June) and autumn (mid-August to mid-October) in forested and field habitats in 34 locations throughout Finland. Individual trapping sites within locations were often not constant over the study period, primarily due to changes in land use. New sites were selected to be as close and as similar as possible to earlier sites regarding habitat characteristics. Trapping was conducted in two habitat types in almost all locations: spruce (Picea abies) or mountain birch (Betula pubescens) forests and in open grassland habitats, primarily old agricultural land no longer in use. We used records between 1978 and 2017, including a total of 15 species sampled in 19 sites after applying our selection criteria ('Study design and data preparation' below).

Butterflies.
We combined two similar surveys of butterflies conducted in agricultural landscapes in Finland. The first is a butterfly monitoring network based on volunteer transects initiated by the Finnish Environment Institute (SYKE) 40 . Between 1999 and 2017, the network included 101 sites, with an average of 47 sites recorded annually (range 30-59). At each site, the walking route (transect) is kept constant from year to year and walked repeatedly during the summer. Along each transect, the number of individuals for each species is recorded from a 5 × 5 × 5 m 3 cube ahead of the observer 41 . The transects are monitored by volunteer butterfly enthusiasts with high species identification skills, who are asked to conduct a minimum of seven annual visits per transect, approximately once a fortnight from late May to late August. Weekly counts are recommended and are usually carried out on nearly half of the transects. The sampling period is typically no longer than 16 weeks and less than ten weeks in the northernmost transects. The second survey type spans between 2001 and 2014 and consists of 68 standardized transects of 1 km length in southern Finland and the Åland islands. These transects were sampled by researchers with a constant sampling frequency of seven counts per summer 42 . The median transect length of the combined data is 1.95 km (mean = 2.41 km). We used records between 1999 and 2017, including a total of 68 species sampled in 98 transects after applying our selection criteria ('Study design and data preparation' below). 43,44 , which is coordinated by the Finnish Environment Institute (SYKE). Nocturnal moths were sampled using 'Jalas' light traps 45 equipped with 125 W Hg vapour or 160 W mixed-light bulbs, located mainly in forested areas across Finland. Traps are located in the same location from year to year and are usually emptied weekly. Sampling occurred every night from early spring to late autumn (usually between April and October). Sampling effort (that is, trapping period) was constant across years for each trap, but given that sampling aimed to cover the entire moth activity period at each location, the trapping period was longer in more southern traps. Volunteers empty the traps and identify the specimens 43 , with a variable number of traps being sampled per year. The taxonomic skills of the volunteer lepidopterists were typically excellent, and data quality control and cross checking was carried out by the monitoring coordinators 43,44,46 . The data used here consist of species records collected from 65 traps with at least eight years of sampling. We used records between 1993 and 2017, including a total of 615 species after applying our selection criteria ('Study design and data preparation' below).

Understory vascular plants of forests.
Understory vegetation was surveyed on a systematic network of 1,700 sites established on mineral soil in forested land between 1985 and 1986 (as part of the 8th Finnish National Forest Inventory 47 ). This network consists of clusters, each including four sites with 400 m intervals. The clusters were located 16 km apart from each other in southern Finland, and 24 km and 32 km apart in northern Finland along east-west and north-south axes, respectively. All 1,700 sites were resurveyed in 1995, and a subset of 443 sites were resurveyed in 2006. The spatial extent of sampling was comparable across surveys covering the whole country. In all three surveys, vascular plant species (dwarf shrubs, herbs, ferns and graminoids, including also tree and shrub seedlings and saplings up to 50 cm tall) were identified and species' cover (0.1-100%) was visually estimated; this was based on three to six permanent square-shaped sampling plots of 2 m 2 , located 5 m apart from each other within each site. The data are curated by the Natural Resources Institute Finland (Luke). For this analysis, we selected all sites with four sampling plots. The average of species cover across these four sampling plots is used as an estimate of species abundance at each site. This included occurrence records from 1,518, 1,673 and 443 sites in years 1985, 1995 and 2006, respectively. After applying the selection criteria ('Study design and data preparation' below), the data included a total of 109 species sampled in 1,712 sites.
Phytoplankton. The National Finnish Phytoplankton Monitoring Database maintained by the Finnish Environment Institute (SYKE; open data portal http:// www.syke.fi/en-US/Open_information) comprises nationwide phytoplankton community data of lake surface water samples. We used data collected in the late summer months with samples taken during early July to late August, reflecting the peak productivity season of lake phytoplankton communities. To ensure consistent sampling methodology, we included only data between the years 1977 and 2017. All phytoplankton samples were preserved with acid Lugol's solution and analysed using the standard Utermöhl technique 48 . We used records between 1978 and 2017, including a total of 464 species sampled in 1,544 sites after applying our selection criteria ('Study design and data preparation' below).

Study design and data preparation.
Before running the joint species distribution models (below), we converted abundance data into presence records. Each site was assigned to one of the four bioclimatic zones in Finland 49 -from south to north: hemiboreal (HB), southern boreal (SB), middle boreal (MB) and northern boreal (NB). We combined the two southernmost regions by pooling the occurrence records to obtain a better distribution and number of samples given the much smaller extent of the HB zone. Each occurrence record was also assigned to a different decade based on the year of sampling: decade 1 (1978-1987), decade 2 (1988-1997), decade 3 (1998-2007) and decade 4 (2008-2017) (Fig. 1). Scarce records before 1978 were excluded. Splitting the data into discrete zone and decade subsets allowed us to use independent (and computationally manageable) data to jointly model species responses within each taxonomic group and to disentangle any contrasting imprints of climatic changes between regions and periods. Each subset therefore covered a wide range of climatic conditions for each taxon, with the majority including ten years of data. While the number of sites varied over decades, zones and taxa, the change was not systematic, neither over time nor across taxonomic groups-that is, there was no consistent pattern of more sampling in later decades, and for each group, there could be more sampling sites in a given zone in a later or in an earlier decade (Supplementary Table S2). In addition, the frequency distribution of the pairwise distances between all sites remained similar across decades in each zone ( Supplementary Fig. S4), suggesting no changes in the spatial aggregation of sites over time. Furthermore, our analyses are model based and thus explicitly account for the number and distribution of sampling sites, while making inference on both environmental covariates and spatial random effects 50 . Thus, the variation in the number and distribution of sampling sites affects the uncertainty on trend estimates, rather than affecting the estimates themselves.
Species were included if they had a minimum of ten occurrences in each zone × decade combination. Finally, due to the difficulty of obtaining reliable model estimates from very sparse data, we set two additional criteria, including only data subsets with at least 20 samples and at least six species; this meant that despite data being available, some subsets were not included in the analyses (for instance, small rodents in NB).
Environmental data. For each site across the different taxonomic datasets and for each year sampled, we extracted values of daily mean temperature, daily precipitation sum and daily snow depth from the Finnish Meteorological Institute (https://etsin.fairdata.fi/datasets/fmi?keys=Finnish%20Meteorological%20 Insitute&terms=organization_name_en.keyword&p=1&sort=best; first accessed in April 2019 and updated in May 2020). These datasets are part of the Finnish Meteorological Institute ClimGrid, which is a gridded daily climatology dataset of Finland, with a spatial resolution of 10 × 10 km 51 . From these, we calculated values of annual mean temperature, total precipitation and number of days with snow cover. We extracted annual NAO values from the Climate Analysis Section, National Center for Atmospheric Research (NCAR), Boulder, Colorado, United States 52 . The NAO index is calculated based on surface sea level pressure difference between the Subtropical (Azores) High and the Subpolar (Iceland) Low, with a high index (NAO +) indicating cool summers and mild and wet winters, whereas low values (NAO −) indicate cold dry winters 53 . We explored whether other climatic variables might also be relevant for our analysis. Specifically, we additionally calculated January temperature, July temperature, temperature range, mean temperature standard deviation, temperature seasonality, growing degree days (over 5 °C), summer precipitation, precipitation range and mean snow depth. For calculating these additional variables, we extracted daily maximum and minimum temperature data from the ClimGrid data. We evaluated the correlation patterns among these variables and found they were highly correlated, particularly with annual values (Supplementary Fig. S5). As such, we included annual mean and summed values in our models because the relevance of the more detailed variables is likely to vary among taxa. Using annual values also facilitates comparisons with other studies and climate scenarios and allows overcoming issues regarding the overlap between some variables and seasonal sampling of the different species surveys.
Quantifying climatic change patterns. We used two approaches to quantify changes in the different climatic variables. First, we used linear regressions with an interaction term between decade and bioclimatic zone to test whether changes over time differed between the different zones in our analysis framework. Second, for a more spatio-temporally resolved assessment of changing patterns we used the k-means clustering method to characterize regions of common climatic profiles for each variable considering all available data over space and time. We set k = 4, resulting in four groups of respectively similar variable conditions. Subsequently we calculated the mean and standard deviation of each resulting cluster for all variables and highlighted the average climate regimes over the whole latitudinal gradient and decades.
Joint species distribution modelling. We used a joint species distribution modelling framework to (1) determine how the different climatic variables affected species occurrence patterns and (2) assess whether their relative importance in structuring assemblages has changed over time or across latitude. We fitted separate spatially explicit models for each combination of taxonomic group × bioclimatic zone × decade, yielding a total of 63 models. We modelled the probability of species presence in response to temperature, precipitation and snow cover with quadratic and to NAO with a linear function. Because we model presence data, we used probit regression models. Towards (1), we used the species responses to the climatic variables to quantify the proportion of species at the lower end of their niche (that is, occurrences increasing along the climatic gradient), at the upper end of their niche (that is, occurrences decreasing along the gradient) or at the optimum of their niche (that is, occurrences peaking within the gradient; Fig. 2b; 'Scoring species' position within niche domains' below) within each bioclimatic zone and decade. Towards (2), we compared the proportion of explained variance attributed to each variable and examined whether their relative contribution shifted through time and/or space.
For each taxon × bioclimatic zone × decade combination, we fitted latent-variable joint species distribution models using the Hierarchical Modelling of Species Communities (HMSC) framework. HMSC is a multi-variate Bayesian generalized linear mixed-effect model framework, which allows joint modelling of the responses of entire species assemblages and explicit modelling of spatial and temporal autocorrelation 18,54,55 . We used spatially structured latent variables which were originally proposed by Ovaskainen et al. 56 . and later expanded to big spatial data by Tikhonov et al. 57 . We fitted the models with the 'Hmsc' v 3.0.9 package 55 in R 58 with a probit link function and assuming the default prior distributions. As fixed effects, we included the climatic variables described above, estimating a second-order polynomial term for all covariates except for NAO, for which we estimated a linear term only. To account for variation in other (unmeasured) environmental variables and potential year-to-year variation not captured by the climatic covariates, we included the random effects of site and year, respectively. All models had the same structure for all the taxon × zone × decade subsets, except for the understory vegetation data for which we did not include the covariate NAO nor the random effect of year because these data were collected only in three individual years, corresponding to a single year per decade in our analytical framework (Fig. 1). Finally, due to computational bottlenecks for large data subsets, some model runs failed to complete with available resources; when this was the case, we randomly subsetted 1,000 records before re-fitting the models (specifically, bird and phytoplankton subsets for the last decade in SB and six mammal subsets for MB and SB in the second, third and fourth decades). We performed posterior sampling using four Markov Chain Monte Carlo chains, each collecting 250 samples, yielding a total of 1,000 samples. We used a thinning interval of 100 and excluded the first 12,500 iterations as burn-in, only sampling the subsequent 25,000 iterations per chain. For phytoplankton in the southernmost region in decade 4, we used a thinning of 10 due to computational constraints due to large site and species numbers. To evaluate Markov Chain Monte Carlo convergence, we examined the distribution of the potential scale reduction factor over the parameters related to the fixed effects and the random effects (equivalent to the Gelman-Rubin statistic 59 ). We assessed model fit via the Area Under the Curve (AUC) statistic 60 and model discriminatory power was quantified by Tjur's R 2 , which is recommended as a standard measure of discriminatory power for binary outcomes 61 .
To quantify shifts in the explanatory power associated with each covariate, we assessed variance component estimates, that is, the relative explanatory power of each environmental covariate in the HMSC models 18,54 . We estimated how the relative importance of the covariates in explaining species occurrences varied over time by fitting linear regression models to the species variance component estimate values as a function of decade (using the function 'lm' in R) and then compared these changes across zones for the different taxonomic groups (Fig. 4). These model comparisons were carried out after weighing the variance component values by each model's ability to explain species occurrence patterns (that is, discriminatory power quantified using Tjur's R 2 values).

Scoring species' position within niche domains.
To analyse whether a species occurred at the lower end, at the optimum or at the upper end of its climatic niche within a particular bioclimatic zone and decade, we assessed the species' responses to each of the climatic variables as follows. First, we classified a species as non-responsive to a specific climate variable within the measured range of that variable if the posterior distribution of the corresponding beta parameter estimates included zero with a probability of more than 10% (corresponding to having less than 90% posterior probability for the response). The non-zero responses were then classified as positive, negative or 'bell-shaped' based on the sign of the derivative of the response over the observed environmental gradient. A positive response corresponds to a species being at the lower end of its niche, 'bell-shaped' response corresponds to a species being at the optimum of its niche, and a negative response corresponds to a species being at the upper end of its niche. In cases where the derivative is positive/negative at both ends of the environmental gradient, responses were classified as either positive or negative, respectively. Cases where the derivative changed from positive to negative required subsequent evaluation. More specifically, if the derivative was positive or negative over less than 20% of the gradient, we classified the response as negative or positive, respectively. If the derivative was positive or negative over more than 80% of the gradient, we classified the response as positive or negative, respectively. Finally, if the derivative was positive or negative over at most 60% of the environmental gradient, we classified the response as bell-shaped. We evaluated whether this threshold affected the overall results by implementing the same classification using two other criteria for the derivatives being positive or negative: over less than 10% and more than 90% of the environmental gradient, and over less than 30% and more than 70% of the gradient. This showed that our classification of species' responses was robust to these choices ( Supplementary Fig. S6). This classification procedure did not apply to the beta parameters for NAO because we did not include a polynomial term for this covariate, as explained above. Thus, we obtained the number of species for each taxon × bioclimatic zone × decade model that showed responses to the different covariates and calculated the proportion of these species relative to the total number of species in each model ( Fig. 3 and Extended Data Fig. 2).
Comparing overall species composition between decades. We compiled the species list present in each decade and each zone for the different taxonomic groups and compared these lists between consecutive decades, that is, comparing decades 1 and 2, 2 and 3, and 3 and 4. For each comparison, we noted how many species were present in both decades ('shared' = A), were present only in the first decade in the comparison ('unique to earlier decade' = C) or were present only in the last decade in the comparison ('unique to later decade' = B). We did this exercise for all taxa in all zones, plotting the sum of 'shared' and 'unique species in each decade' (Supplementary Fig. S2) and all the consecutive decade comparisons for each taxon (Supplementary Fig. S1). To quantify these patterns over a larger temporal extent, we implemented the same procedure but only comparing the first and last decades sampled for each taxon (Supplementary Figs. S1 and S2; note that for butterflies, these analyses are identical, because this taxon was sampled only in decades 3 and 4).
Quantifying community dissimilarity. To assess how community composition changed over space and time, we calculated overall dissimilarity among all the sampling units within a given zone and decade-that is, we quantified variation in composition among sites within a given spatio-temporal extent, regardless of their location. Dissimilarity indices range between 0 and 1, representing cases where all or no species are shared between sites, respectively. We used the same occurrence matrices that were analysed with the HMSC models, that is, the raw species data matrices. We used the function 'beta.sample' in the 'betapart' package v 1.5.2 62,63 to calculate total dissimilarity (Sørensen index), which can be additively decomposed into the turnover (Simpson index) and nestedness components 64 . 'beta.sample' randomly selects a specified number of sites to generate distributions of the multiple-site dissimilarity measures. This is important because the number of sites affects the estimated compositional change values. For each taxonomic group, we first determined the minimum number of sites among the different zone × decade combinations, which was used to define the number of sites to be randomly sampled from the original occurrence matrix, performing this subsampling 1,000 times. We then plotted the mean and standard deviation of these distributions to compare compositional change for each taxonomic group across the different zones and decades. We focus on the turnover metric (that is, species replacement among sites independent of changes in species richness; Extended Data Fig. 3), as it was systematically the main component of dissimilarity except for the small rodent data where the nestedness component had a relatively higher contribution to total dissimilarity. We show the results for the three metrics in Supplementary Fig. S3.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data used in the current analyses are available in GitHub (https://github. com/benweigel/RECChange) and in an online archive at Zenodo 65 . Moth monitoring data is available through the Finnish Biodiversity Information Facility, FinBIF (https://laji.fi/en/observation/list?sourceId=KE.1501). Phytoplankton data was obtained from the open access data service of Finnish Environment Institute (SYKE). Raw data is available from data owners on reasonable request. For the Natural Resources Institute Finland (Luke) any data requests should be sent to kirjaamo@luke.fi. The climatic variables data is available from the Finnish Meteorological Institute (https://etsin.fairdata.fi/datasets/ fmi?keys=Finnish%20Meteorological%20Insitute&terms=organization_name_ en.keyword&p=1&sort=best). NAO values were exported from the National Center for Atmospheric Research (https://climatedataguide.ucar.edu/climate-data/ hurrell-north-atlantic-oscillation-nao-index-pc-based).

Code Availability
Code to reproduce the analysis is available in GitHub (https://github.com/ benweigel/RECChange) and in an online archive at Zenodo 65 .