Canopy closure and intensifying climate extremes drive understory species loss over 25 years of forest monitoring

Forests cover 31% of the world’s land surface¹, hosting the majority of Earth’s terrestrial biodiversity². Especially in temperate forests, the vascular plants occurring in the understory layer represent a crucial component of their biodiversity³. This component includes more than 80% of the overall plant species richness, showing great sensitivity to environmental changes, and influencing tree regeneration and biogeochemical cycles^4,5. This huge diversity, which is vital for ecosystem functioning and stability, is increasingly facing multiple threats, including climate change, air pollution, and natural or anthropogenic disturbance dynamics (e.g., tree canopy defoliation, management intensification or abandonment)^6,7,8,9, with cascading effects on forests’ capacity to sustain environmental processes underpinning human health and quality of life¹⁰.

Climate warming can drive the thermophilization of plant communities¹¹ and may also lead to the reduction of understory plant species richness¹². Increased soil resource availability – often associated with nitrogen deposition – enhances understory plant species richness by increasing the number of ubiquitous, eutrophic, and acidophilous species¹³. Changes in tree canopy cover are responsible for understory species compositional shifts, with species-poor communities often under dense forests due to low light-availability^14,15,16 (but see Chianucci et al.⁸). Despite the advancement of knowledge regarding the understory’s response to environmental changes, several key issues remain open. First, understory communities are expected to exhibit context-dependent responses to environmental drivers due to inherent site characteristics and the varying influences of multiple local and regional factors (e.g., a unique combination of soil and climatic properties, disturbance history, soil type, canopy density, tree composition)^17,18,19. Indeed, when analyzing forest vegetation shifts after 40 years, Wrońska-Pilarek et al.¹³ found clear patterns only when focusing on single forest types separately, since the magnitude and direction of understory changes varied greatly between forest types. This implies that research outcomes deriving from a given forest type cannot be automatically extended to other forests and highlights the importance of targeted investigations on forests characterized by similar ecological features to evaluate the impact of environmental changes on plant diversity^17,20. Second, a large amount of the research findings come from studies based on space-for-time substitution approaches^21,22 or studies involving the use of semi-permanent plots to resurvey historical locations where vegetation surveys were previously conducted^23,24. Although the resurvey approach with semi-permanent plots provides valuable insights, it is associated with notable limitations, including inconsistencies in sampling timing, relocation errors, and observer-related biases, which can affect data reliability^25,26. Therefore, the use of permanent plots is strongly recommended in repeated surveys (see de Bello et al.²⁷ and references therein). However, even when permanent plots are used, vegetation studies are often based on two or few survey events^9,14,28,29, overlooking the effects of changing environmental conditions^30,31 (but see Bernhardt-Römermann et al.³²). The lack of continuous time series of monitoring vegetation changes fails to account for the intrinsic dynamics of plant communities, including their frequently delayed responses to environmental variations (see Richard et al.³³; Wu et al.³⁴ and references therein), thereby increasing uncertainty in resurvey outcomes²⁵. Additionally, vegetation changes over time can be the result of fluctuation events, that is, non-directional and often transient changes in community composition reflecting natural dynamics rather than consistent trends^33,35.

Vegetation change is most commonly evaluated by alpha diversity metrics, with species richness being by far the most widely applied measure of biodiversity^36,37. Its prevalence reflects both its conceptual clarity and its practical advantages. Species richness provides a direct and intuitive representation of biodiversity change, as local declines in the number of species may indicate local extirpations or failures to maintain viable populations³⁷. Furthermore, species richness is methodologically feasible and robust in long-term monitoring frameworks. Data required to estimate local richness can be collected consistently over time using standardized field protocols, and a wide range of well-developed statistical approaches exists to account for sampling effects and detection uncertainty (see Fletcher et al.³⁷ and references therein). These features make species richness particularly suitable for temporal analyses aimed at detecting directional changes in biodiversity, which require comparable sampling designs through time. The widespread use of species richness, particularly in conservation contexts³⁸, is further reinforced by its integration and central role in international policies³⁹. However, stability in species richness over time can mask even drastic changes in species composition^35,40. Therefore, incorporating other diversity metrics, such as beta diversity^31,35,41, is crucial for a more comprehensive understanding of temporal community dynamics.

Temporal variation in community composition (i.e., beta diversity) can be partitioned into two additive components: species turnover and nestedness⁴². This distinction is fundamental because any difference in species compositions between two communities can arise from only two elementary processes: species replacement and species loss (or gain). Consequently, these two components are informative of different underlying ecological mechanisms. Turnover reflects changes in community composition driven by the replacement of species over time, whereby species lost between two time points are replaced by different species. High turnover is therefore expected to be associated with variation in environmental conditions that favour different species at different times, reflecting shifts in competitive hierarchies and immigration–extinction dynamics^22,43. In such a case, community composition changes substantially, while overall species richness may be maintained through compensatory replacement. In contrast, nestedness reflects directional species loss without replacement, indicating progressive community impoverishment. Thus, it arises when communities observed at later time points constitute subsets of past, species-richer communities. This pattern can emerge under changes in environmental conditions imposing ordered disassembly of communities, where specialist or late-successional species are selectively lost^42,43. Such phenomena are often associated with increasing environmental stress and habitat simplification (see Liu et al.⁴⁴; Wang and Yang⁴⁵ and references therein). Accordingly, nestedness is expected to increase when extinction processes dominate over colonization, and to decrease when compensatory dynamics restore previously lost species (recolonization). Disentangling turnover and nestedness is therefore essential for linking compositional change to underlying ecological mechanisms, as similar levels of overall beta diversity may arise either from dynamic species replacement or from directional diversity loss^42,43. Moreover, by analyzing these two components at two different temporal scales (i.e., immediate and cumulative), it is possible to further distinguish between interannual variations and long-term trends in turnover and nestedness. Distinguishing between these two time scales enables inferences of whether observed compositional changes represent transient interannual fluctuations or consistent directional reorganization of forest communities over time (see Di Cavalho et al.⁴⁶). Interannual fluctuations in composition may arise from stochastic demographic variation(i.e., the equilibrium dynamics⁴⁷), conversely, long-term directional changes may reflect persistent shifts in environmental conditions or successional trajectories (see Hillebrand et al.³⁵ and references therein).

In this context, using permanent monitoring plots repeatedly surveyed over a 25-year period (1999–2023), we investigated long-term patterns and drivers (i.e., climate, soil, and forest structure) of forest understory plant species alpha diversity, as well as changes in beta diversity at different time scales. Forest plots were classified into four different biomes⁴⁸: boreal (i.e., alpine coniferous forests of the temperate oro belt⁴⁹), Mediterranean (i.e., sclerophyllous evergreen forests⁴⁹), nemoral beech and nemoral oak (i.e., temperate deciduous forests⁴⁹). Our study utilized the Italian monitoring plots of the ICP Forests network (International Co-operative Program on the Assessment and Monitoring of Air Pollution Effects on Forests, http://icp-forests.net/page/icp-forests-manual), which are distributed across a broad biogeographical gradient, from the Alps to the Mediterranean. This distribution encompasses diverse forest types and a wide variety of vascular plant species⁵⁰, making Italian forests an optimal case study for assessing vegetation dynamics across different environmental conditions. Although temporal case studies from forests of central and northern Europe exist³², southern Europe remains largely overlooked, despite being a hot spot for global changes⁵¹.

We asked the following three questions:

Q1) What are the temporal patterns of understory species richness across different biomes? Does species richness change significantly over a 25-year continuous monitoring period?

Q2) What are the key environmental drivers of understory species richness in different biomes? Do climate, soil, and forest structure influence species richness similarly across the studied biomes?

Q3) To what extent do turnover and nestedness drive temporal trends in understory beta diversity across biomes, and how do these dynamics vary between time scales – i.e., reflecting interannual fluctuations vs. long-term trends?

The outcomes of this research will help understand the ongoing forest diversity dynamics and their drivers, providing valuable insights for science-based conservation policies and forest planning, in alignment with international policies^39,52,53.

The total number of resurveys conducted over the study period (1999–2023) was as follows: 51 for the boreal biome, 29 for the Mediterranean biome, 101 for the nemoral Beech biome, and 96 for the nemoral oak biome. Over the 25-year period, in the boreal biome, understory species richness ranged from 27 to 57 species, with a mean of 42 species. In the nemoral beech biome, understory species richness varied between 18 and 33 species, with a mean of 26 species. For the nemoral oak biome, the number of understory species ranged from 29 to 55, with a mean of 43 species. In the Mediterranean biome, understory species richness varied from 16 to 29 species, with a mean of 23 species. The complete values of mean species richness, separated by biome and year, are provided in Table S1.

Trends in understory species richness and their environmental drivers

The species richness of understory vascular plants showed a significant decline from 1999 to 2023 in three out of four biomes (Fig. 1). Specifically, this trend was evident in the boreal (estimate= −0.69, p < 0.01), nemoral oak (estimate= −0.34, p < 0.001), and nemoral beech (estimate= −0.20, p < 0.001) biomes. In contrast, the Mediterranean biome showed no significant trend in the understory species richness (estimate= 0.04, p > 0.7). Aware that the distribution of surveys was not balanced across the years, we repeated the analyses considering only the 1999–2014 period, but the trends of understory species richness remained consistent in all biomes (see Table S2).

Climate and forest structure were the main drivers of these trends (Table 1). In the boreal biome, understory species richness showed a negative relationship with both the tree and shrub cover. In the nemoral beech biome, understory species richness decreased with increasing numbers of consecutive dry days during the growing season and with higher tree cover. A negative effect of soil pH on species richness was observed only in this biome. In the nemoral oak biome, understory species richness declined with increasing precipitation seasonality, the annual frequency of hot days, and higher shrub cover. In the Mediterranean biome, understory species richness decreased with more consecutive dry days during the growing season. However, as previously mentioned, no significant temporal trend in species richness was observed in this biome over the study period.

Temporal beta diversity

Significant temporal trends in the two beta diversity components (i.e., turnover and nestedness) were observed across the analyzed biomes at the considered time scales, highlighting shifts in the composition of understory vascular plant communities over time (Table 2). Specifically, the boreal biome showed a negative temporal trend for turnover at the immediate time scale, while a positive trend was observed for nestedness at the cumulative time scale. In the nemoral beech biome, both turnover and nestedness exhibited positive trends at the cumulative time scale. In the nemoral oak biome, turnover showed a positive trend at the cumulative time scale, whereas nestedness exhibited a positive trend at the immediate time scale. Finally, in the Mediterranean biome, turnover showed a positive trend at the immediate time scale.

Overstory dynamics of canopy closure over time appear to be a key driver of the understory species richness decline in boreal and nemoral forests, as also found by approaches based on space-for-time substitution and semi-permanent plots^54,55,56,57. This is probably related to the natural process of canopy regeneration after the last logging event (open-to-close dynamic¹⁴), given that the permanent plots of the ICP Forests Level II network were not logged during the monitoring period. In many plots, the last logging event occurred 50–60 years before the first sampling year (i.e., 1999)⁵⁸. Our interpretation is confirmed by the significant increase in tree cover over time, particularly for the boreal biome (Table S3). The occurrence of permanent monitoring plots undergoing strong canopy-closure dynamics is not a bias of our study; in fact, this dynamic is mirroring one of the most widespread and relevant phenomena ongoing in European forests, i.e., the abandonment of active forest management and the inherent increase in canopy cover^59,60. Therefore, if the current socio-economic trends are consistent in future decades, we may expect a decline of species richness in nemoral and boreal European forests with the abandonment of forest management practices and related progressive canopy recovery. Opposite patterns should be expected in case of canopy opening due to natural disturbances, tree mortality, or revival of forest management practices^61,62.

The reduction of species richness is not necessarily detrimental to the conservation of forest understory diversity. Indeed, species richness can be a misleading diversity metric for conservation purposes since it does not consider species identity and functioning^63,64,65. Boch et al.⁶⁶ found that high plant species richness indicates forest disturbances rather than conservation status. Furthermore, Hilmers et al.⁶⁷ observed a U-shaped response of species richness to forest succession, with peak diversity values in early and late stages. We therefore do not exclude a future increase in species richness if the monitored forests are left undisturbed.

We speculate that light-demanding species are progressively filtered out by the increase in canopy cover, while shade-tolerant species are facilitated, as shown by many researchers^60,68. Indeed, this is supported by the analysis of Ellenberg indicator values, which shows a decreased mean light value over time (Table S4). Shade-tolerant species are often forest specialists, sometimes of conservation interest, while light-demanding species are often generalists or typical of open ecosystems⁶⁹. Our speculation is supported by data observation, which reveals that some non-forest species (sensu Heinken et al.⁷⁰) – like Holcus lanatus L. and Rhinanthus minor L. – were present in the nemoral biome plots during the early years of monitoring, but later disappeared. However, this aspect requires a further and more specific evaluation⁷¹.

Our results showed a significant effect of climate – particularly longer dry periods, higher frequency of hot days, and greater precipitation seasonality – in reducing species richness in the nemoral biomes. This suggests that the occurrence of extreme climatic events has a stronger effect than variations in total precipitation or mean temperature. Indeed, several studies have highlighted the predominant role of climatic extremes in shaping biodiversity responses^72,73,74. However, it is important to note that our analysis focused on macroclimatic variables without considering the potential buffering effect of canopy cover. In fact, canopy cover may favour more stable microclimatic conditions (i.e., the climate experienced by the understory layer), as higher canopy cover is expected to generate a larger climate buffering capacity, thus increasing the difference between the macro- and the microclimate^19,75.

Our findings point toward a potential shift in forest understory communities toward xerophilization, consistent with findings by Wrońska-Pilarek et al.¹³, and thermophilization, as reported by Zellweger et al.⁷⁶ and Braziunas et al.¹². This is further supported by the analysis of Ellenberg indicator values, which revealed a decrease in moisture values and an increase in temperature values over time in the nemoral oak biome (Table S4).

Soil variables showed a limited effect on understory species richness despite previous evidence⁷⁷. However, our results indicate a negative effect of soil pH on species richness in nemoral beech forests. Probably, these forests have achieved a balanced cycle of uptake and release of basic cations (BC; i.e., K, Mg, Ca), so that the soil buffer capacity tends to improve, thanks to the high levels of atmospheric BC deposition⁷⁸. Changes in soil pH are known to influence both plant species richness and functioning, but with context-dependent effects^21,79. However, we must acknowledge that our soil data were not available for all sampling years of the vegetation surveys, which may have constrained the assessment of the overall contribution of soil predictors to species richness.

Contrary to boreal and nemoral forests, Mediterranean forests did not show significant changes in understory species richness over time. This may depend on one or a combination of the following aspects. First, the understory of Mediterranean forests is largely composed of woody species, which are usually more persistent across years than herbaceous species⁸⁰. Second, the canopy cover of the sampled Mediterranean forests remained relatively stable throughout the monitoring period (see Table S5), suggesting that the understory species pool has not been subject to filtering effects from changes in canopy dynamics.

Another consideration concerns the amount of variance explained by the models. We found that the fixed factors explained less than 6.9% of the overall variance in all models. This suggests that a larger portion of the variance depends on the specific features of the study plots (i.e., the random factors), further confirming that local environmental conditions and disturbance history of each plot play a crucial role in determining diversity patterns^81,82,83.

The complex temporal dynamics in species turnover and nestedness across biomes emphasize the multifaceted nature of beta diversity and its scale-dependent patterns. Examining both immediate and cumulative temporal scales, we disentangled key processes shaping vegetation communities over time. In general, both turnover and nestedness contributed to changes in beta diversity at both temporal scales, although with distinct trends and biome-specific exceptions.

In biomes where understory species richness decreases over time, both cumulative turnover and nestedness are observed. This suggests consistent increasing trends of beta diversity patterns during the monitoring period. In detail, the boreal biome showed an increasing trend in cumulative nestedness, indicating a progressive process of species filtering, leading to communities that are increasingly subsets of their original state^42,84. In contrast, no evidence in the cumulative turnover was observed, while a decrease in immediate turnover was shown, suggesting that the rate of species replacement is stable considering the long-term perspective and is reducing considering consecutive years. These findings might be related to a stabilization process where nestedness is of overwhelming importance in the long-term perspective and is potentially associated with environmental-induced constraints, such as the reduction of light availability due to higher canopy cover.

Nemoral biomes exhibited increasing trends in cumulative turnover, suggesting that species replacement is the most important process behind understory community changes over time. Many researchers confirm that species turnover is a common response to substantial environmental changes from the initial conditions in terms of forest structure and resource availability⁶². The nemoral beech biome was also characterized by cumulative nestedness, revealing that these forests are experiencing a noticeable and complex reshaping of the understory communities over time.

The Mediterranean biome, characterized by understory species richness stability over time, showed a compositional dynamic marked by increasing turnover acting at the immediate time scale only. This pattern suggests a “carousel model”⁸⁵, where species enter and exit the plots without leading to a cumulative trend over time. Such a mechanism reflects a dynamic equilibrium where the community maintains richness but is subject to interannual compositional fluctuations probably related to stochastic events or interannual climate variability. This interpretation is further supported by recent findings from Jandt et al.⁴⁰ and Bonari et al.³¹, who demonstrated that stable species richness can conceal underlying trends of biodiversity reorganization, often resulting from species-specific responses to environmental changes.

In conclusion, in this paper, we used permanent monitoring plots spanning a wide environmental gradient to assess 25-year trends and drivers of forest understory plant species richness. We found that forests belonging to boreal and nemoral biomes (i.e., coniferous alpine forests and temperate deciduous forests, respectively) appear more sensitive to global changes than Mediterranean forests. Specifically, boreal and nemoral forests are experiencing a decline in understory species richness, associated with increased canopy cover—partly a legacy of past forest management practices—as well as with greater precipitation seasonality, longer dry periods, and heat waves. In contrast, Mediterranean forests, which are generally adapted to recurrent droughts and characterized by a more stable canopy cover, showed stable levels of species richness over time.

Beta diversity analyses revealed that the processes behind species richness patterns vary by biome and temporal scale. In biomes with declining richness, both turnover and nestedness operate over the long term, suggesting deterministic processes such as environmental change and progressive species filtering. In contrast, Mediterranean forests exhibited interannual fluctuations without long-term trends in beta diversity, indicating a dynamic equilibrium likely driven by stochastic events or interannual climate variability.

Overall, we highlight the importance of long-term monitoring using permanent plots repeatedly sampled over time to uncover changes in diversity, as well as the processes and drivers behind these changes at different temporal scales. Despite the relatively high sampling effort and the costs associated with monitoring infrastructure, these types of data represent an unprecedented resource for understanding the site-specific effects of global change. For instance, continued long-term monitoring will be essential to assess how future canopy stabilization in boreal and nemoral biomes will interact with climate change in shaping understory dynamics. Additionally, continued monitoring can assess whether Mediterranean forests will experience shifts in understory species richness if tipping points related to climate change are reached.

Our findings also highlight the importance of developing adaptive forest management strategies considering site-specific responses to environmental change.

Study Area and Vegetation Sampling

The study considers 31 permanent monitoring plots established under the ICP Forests Level II framework in Italy. These plots represent the intensive forest monitoring program and are key for providing insight into cause-and-effect relationships between the condition of forest ecosystems and various stress factors such as air pollution and drought (Fig. 2).

These ICP Forests Level II plots are located within environmentally homogeneous forest areas. Each plot has a dimension of 50 m × 50 m, subdivided into 25 sub-plots of 10 m × 10 m, of which 12 permanent sub-plots were selected in a chessboard pattern as sampling units for vegetation monitoring over time⁸⁶. All the plots are fenced, preventing disturbances from human activities (e.g., trampling and logging), as well as from medium- to large-sized mammals.

The ICP Forests program involves vegetation sampling in permanent plots (at least) every five years, following a standard and harmonized protocol, including crew training and intercalibration exercises⁸⁷. Sampling followed quality assurance procedures⁸⁸. All surveys were carried out during the summer season, when vegetation reached maximum development, collecting species occurrences and abundances (percent cover) according to different forest layers (herb, shrub, and tree layers). Here, we considered only occurrences of vascular plants, including both herbaceous and woody plants up to 0.5 m in height (hereafter, “understory”). Species nomenclature followed Pignatti⁸⁹.

We grouped the 31 plots into homogeneous biomes according to the classification of Mucina et al.⁴⁸, as species response is context-dependent. The classification resulted in 7 plots belonging to the boreal biome (i.e., alpine coniferous forests of the temperate oro belt⁴⁹), 4 plots to the Mediterranean biome (i.e., sclerophyllous evergreen forests⁴⁹), and 20 plots to the nemoral biome (i.e., temperate deciduous forests⁴⁹). Due to the high number of plots within the nemoral biome and the significant differences in terms of climate and vegetation of this biome⁹⁰, we further refined the classification by introducing a finer level⁴⁸. As a result, we identified 10 plots as nemoral oak forest biome and 10 plots as nemoral beech forest biome, reflecting warm temperate and cool temperate climates, respectively (Fig. 2)⁹⁰.

Explanatory variables

Based on existing ecological knowledge on forest response to environmental conditions^21,56 and data availability, we selected, for each specific plots, a set of predictors at different spatial scales, including climate, soil, and forest structure parameters (Table 3).

Climatic conditions for each study plot and year of vegetation sampling were obtained from the E-OBS dataset provided by the Copernicus Climate Change Service (https://surfobs.climate.copernicus.eu). This dataset is available on regular latitude-longitude grids with a spatial resolution of 0.1 degrees (approximately 8 km²) and daily temporal resolution. Using mean daily temperature (°C) and total daily precipitation (mm) data, we computed the following climatic variables: mean annual temperature, mean temperature of the growing season, total annual precipitation, total precipitation of the growing season, and precipitation seasonality. Precipitation seasonality was expressed as the coefficient of variation (CV) and calculated as σ/μ × 100, where σ and μ, respectively, represent the standard deviation and the mean of the total monthly precipitation. This metric reflects the intra-annual (seasonal) variability in precipitation amount, indicating how stable or variable precipitation is across the twelve months of the year. In addition, temperature seasonality was calculated as the standard deviation (SD) of monthly mean temperatures, based on daily maximum temperature (°C) data. Similarly, this metric describes the intra-annual variability in temperature, indicating how stable or variable thermal conditions are throughout the year. Additionally, the following climatic indices were calculated: the number of annual consecutive dry days (CDD_y), the number of consecutive dry days during the growing season (CDD_gs), and the percentage of hot days during the year (TX90p). The Consecutive Dry Days (CDD) index is calculated by identifying the maximum number of consecutive days with daily precipitation less than 1 mm within a given period (https://etccdi.pacificclimate.org/list_27_indices.shtml). The TX90p index represents the percentage of days in a given year during which the daily maximum temperature exceeds the 90th percentile of daily maximum temperatures of the base period (1970–2000 in this study)⁹¹. The TX90p index was calculated using the ‘climdex.tx90p’ function from the ‘Climdex.pcic’ package. These climatic indices are designed to monitor extreme climate events, such as heat waves and droughts⁹². Widely used in previous studies, these indices are also straightforward to interpret^33,92,93. They focus on capturing patterns of extreme conditions that typically occur multiple times a year, rather than rare, high-impact events that may only occur once a decade⁹⁴.

Soil solutions were collected bi-weekly in the field at the study plots, following the standard procedures of the ICP Forests manual (http://icp-forests.net/page/icp-forests-manual). Tension lysimeters were installed in mineral soil to continuously collect gravitational and weakly retained soil water at specific tension levels (i.e., down to −60 kPa), allowing for ongoing monitoring of soil solutions without causing soil disturbance. This enables tracking nutrient availability and pollutants. Specifically, we employed the following variables from the topsoil (at a depth of 20 cm): pH, K, NH₄⁺, NO₃^–, and SO₄^2-. The bi-weekly measurements were aggregated into annual means. Correlation between variables was evaluated using the ‘cor’ function from the R package, and no intercorrelations were detected among the selected variables (Table S6). Soil data were not available for the Mediterranean biome.

Forest stand characteristics, also measured in the field at each plot during the same years as the vegetation surveys, followed ICP Forest procedures (http://icp-forests.net/page/icp-forests-manual). Specifically, data were collected on tree and shrub layer cover (in percentage) as well as average tree defoliation percentage, as these factors influence the sub-canopy microclimate, directly affecting the understory layer^76,95. Tree cover and tree defoliation are fundamental features for assessing forest conditions. Tree defoliation refers to the reduction in the foliage density of a tree’s crown compared to a reference tree assumed to exhibit full foliage density⁹⁶ and is considered an indicator of tree and forest health⁹⁷. In contrast, tree cover refers to the percentage of forest area occupied by the vertical projection of tree crowns^98,99, representing the actual extent of overstory cover¹⁰⁰.

Data analysis

We prepared two datasets for analysis: a species matrix and a predictor matrix. The species matrix comprised understory species occurrences at each plot and in each year, along with the total number of sampling units per plot. Since we aimed to analyze each biome separately, we created four distinct species matrices coupled with four predictor matrices, one for each forest biome. Species richness (S) was calculated from the species matrix by summing the number of species observed across all sampling units within each plot for each sampled year. The predictor matrix comprised all environmental variables for each sampling plot and year. Predictors were standardized using the ‘scale()’ function from the ‘base’ R package, which transforms each variable to have a mean of zero and a standard deviation of one. This standardization ensures comparability among variables with differing units, enhances model stability, and allows for comparable coefficients in statistical analyses.

For the statistical analysis, we fitted linear mixed-effects models (LMMs), using the ‘lme4’ package. The analysis consisted of three major parts, each addressing one of the study’s objectives: the first two focused on species richness as a metric of alpha diversity, while the third focused on Sørensen dissimilarity index as a metric of beta diversity.

Specifically, the first group of analyses (i) tested variation in the species richness over time in each biome. For this purpose, we ran models of the form “S ~ Year + (1|SiteID) + (1 | N.Su)”, where ‘S’ is species richness and ‘N.Su’ is the number of sampling units. The analyses were conducted at the level of individual plots (SiteID), which serve as the statistical units of analysis. This approach was necessary to account for repeated measurements collected over time within each plot and adjust for an unequal number of sampling units across surveys.

The second group of analyses (ii) tested the effects of environmental variables (listed in Table 3) on species richness. For this task, we used LMMs of the form “S ~ Environmental variables + (1|SiteID) + (1 | N.Su)”. Separate models were fitted for each group of environmental variables (climate, soil, forest structure) within each biome, to maximize data retention. In fact, not all environmental variables were measured during every survey or in every plot, and therefore, including all predictors in a single model would have excluded all cases with missing values, substantially reducing sample size and statistical power. Multicollinearity was assessed for each model using the variance inflation factor (VIF) from the ‘car’ package. Thus, none of the fixed factors included in our final models showed intercorrelation issues (VIF < 5)¹⁰¹. Then, we subjected the linear mixed models built with all variables of each group to model selection. Specifically, we applied a bidirectional stepwise selection procedure (both forward and backward) based on the Akaike Information Criterion (AIC)^102,103,104. This approach is based on an algorithm that adds or removes variables depending on whether they decrease the AIC, in order to obtain the model with the lowest AIC value. Model selection was performed with the ‘step’ function in ‘lmerTest’ package¹⁰⁵.

Model assumptions (homoscedasticity, residual independence, and normal distribution of residuals) were assessed through visual inspection of histograms and residual plots, including Q-Q plots and fitted values versus residuals plots^106,107,108. Additionally, spatial and temporal autocorrelation of model residuals were evaluated using Moran’s I and Mantel tests, respectively (Table S7 and S8)^56,109. Temporal autocorrelation was computed as the Pearson correlation between the distance matrix of model residuals and the temporal distance matrix based on sampling year, using 9.999 permutations with the ‘mantel’ function in ‘vegan’ package¹¹⁰. Spatial autocorrelation was evaluated on plot-level residuals using geographical coordinates and distance-based neighbour lists with the ‘moran.test’ function in ‘spdep’ package¹¹¹. Only the boreal biome showed temporal autocorrelation, and therefore, we incorporated an autoregressive correlation structure (corAR1)¹¹². This model was fitted using the ‘lme’ function in the ‘nlme’ package, as ‘lme4’ package does not support correlation structures¹¹³. P-values for the variables in the model were obtained using the ‘summary’ function from the ‘lmerTest’ package¹⁰⁵.

For the third group of analyses (iii), we investigated the temporal trends in beta diversity¹¹⁴ within each biome, focusing on its two primary components: turnover and nestedness⁴². Beta diversity was quantified using the Sørensen dissimilarity index, which is one of the most widely used measures due to its dependence on the proportion of species shared between two communities⁴². The Sørensen dissimilarity index (βsor) is formulated as:

where a is the number of species common to both plots, b is the number of species occurring only in the first plot, and c is the number of species occurring only in the second plot.

To describe species turnover, we used the Simpson pairwise dissimilarity index (βsim), calculated as:

where a, b, and c are the same variables defined above. Since βsor and βsim are equal in the absence of nestedness, their difference represents the nestedness component of beta diversity. Accordingly, the nestedness-resultant dissimilarity (βnes) was derived and formulated as:

Therefore, two distance matrices were computed to account for turnover and nestedness resultant components using the ‘betapart’ package¹¹⁵. Since beta diversity analysis is based on pairwise calculations and multiple sampling events were conducted over time at each plot, the analysis was performed across two temporal scales: immediate and cumulative⁴⁶. For the immediate time scale, beta diversity was calculated between consecutive sampled years (e.g., t₁–t₂, t₂–t₃, …), reflecting interannual changes in community composition. On the cumulative time scale, beta diversity was calculated between the first sampled year and each subsequent sample taken over time (e.g., t₁–t₂, t₁–t₃, t₁–t_n), capturing cumulative (i.e., long-term) changes in community composition. Temporal trends in turnover and nestedness were investigated by building separate LMMs for each biome, with distinct models for each combination of beta diversity component and temporal scale.

For immediate time scale models, beta-diversity values between consecutive years were used as the response variable, and a new variable identifying each pairwise comparison was included as a fixed factor. The ‘SiteID’ and ‘temporal distance’ (e.g., the time interval between consecutive samples t₁–t₂) were treated as random factors to account for the hierarchical structure of the sampling design and heterogeneous repeated measures over time. For example, if three repeated samplings were present for a plot (t₁–t₂–t₃), temporal distance represented the interval between t₁–t₂, and t₂–t₃. For cumulative time scale models, ‘temporal distance’ (i.e., the number of years elapsed since the first sampling) was included as the fixed factor, while ‘SiteID’ was treated as a random factor. In this case, if three repeated samplings were present for a plot (t₁–t₂–t₃), the temporal distance represented the interval between t₁–t₂, and t₁–t₃. This approach allows for assessing whether compositional changes are continuous (e.g., indicating a linear relationship between cumulative turnover and temporal distance) or whether assemblages revert to earlier states⁴⁶. When the LMM did not meet the required assumptions, the dependent variable was transformed using a square-root function and/or a Generalized Linear Mixed Model (GLMM) was employed. This was the case for the boreal biome, where a GLMM with a beta distribution and a logit link function was fitted to investigate temporal trends in nestedness at the cumulative time scale. Nestedness values were adjusted to avoid 0 and 1 values in order to meet the assumptions of the beta distribution, and the analysis was performed using the ‘glmmTMB’ package.

All statistical analyses and visualizations were performed in the R statistical programming language¹¹⁶ (version 4.2.2).