Species-specific climate–growth interactions determine tree species dynamics in mixed Central European mountain forests

Increasing growing season temperatures and the seasonal redistribution of precipitation due to climate change have recently been recorded across the globe. Simultaneously, increases of severe droughts and windstorm frequency have also been documented. However, the impacts of climate change on tree growth performance and fitness might largely differ among coexisting species. Consequently, ongoing temperature increases could lead to extensive changes in tree species compositions in many forest biomes including temperate mountain forests. In this study we used an extensive dataset of 2824 cored trees of three species from two sites, and parameterized a purely climate driven process-based model (Vaganov–Shaskin) to simulate the growth dynamics and climatic limitations of coexisting Picea abies, Fagus sylvatica and Abies alba in two of the oldest mountain forest reserves in Central Europe (the Boubín and Žofín Primeval Forests). We assumed that the species composition reflects climatic growth limitations, and considered between-site differences in mean temperature due to elevation as a model of future climate change effects on mountain forests. Our results show a complexity of site- and species-specific responses of Central European forests to climate change. Over the last 70 years, the proportion of F. sylvatica in Central European natural forests has increased at the expense of conifers. During the investigated period, we observed an increase in the growth rates of the studied species mainly at the higher elevation site, while for the lower elevation site there was increasing intensity of moisture limitation. Despite being the most moisture-limited species, P. abies showed the highest simulated growth rates. In contrast, A. alba was the least moisture limited of all considered species. Given its recent proportion in the forest species composition and intermediate drought resistance, we anticipate the future expansion of F. sylvatica in Central European mountain forests.


Introduction
The composition of Central European mountain forests has changed plastically over the Holocene, reflecting the Quaternary climate cycle (e.g. Tinner and Lotter 2006, Fletcher et al 2010, Bobek et al 2019. Birch-pine forests predominating in the mountain areas at the beginning of the Holocene gradually shifted to predominantly Picea abies in the mid-Holocene. Since that time, a gradual expansion of Fagus sylvatica and Abies alba has been observed, causing significant changes in the forest structure and disturbance dynamics (Bobek et al 2018). During the last century, however, the dynamics of temperate forests have accelerated unprecedentedly, with an even faster expansion of F. sylvatica and decline of A. alba (e.g. Vrška 2008, Elling et al 2009). These changes have been happening mostly due climate change and associated direct and indirect human impacts (Elling et al 2009, Seidl et al 2011. According to climate scenarios, Central Europe is one of the most quickly warming parts of the world (Lewis et al 2019), with temperatures at the end of the 21st century expected to be more than 4 • C higher than the 1961-1990normal (Fischer and Schär 2010, IPCC 2014. Therefore, ongoing increases of growing season temperatures, related to higher evapotranspiration and soil moisture deficits, will bring long periods of drought that will recur at higher frequencies in low and mid-elevation forests (Briffa et al 2009, Spinoni et al 2015. In addition to direct climatic stress, Central European forests might also be affected indirectly, e.g. through an increasing frequency of both biotic and abiotic disturbances (Elling et al 2009, Seidl et al 2011. However, in contrast to the limiting effects of climate, temperature increases might stimulate an earlier start and extension of the growing season in cold forests (Menzel et al 2006, Schwartz et al 2006, Jin et al 2019. These environmental changes will most likely cause a response of tree species distribution patterns (Martín-Benito et al 2010, Buras and Menzel 2019), and may also have consequences on the productivity of Central European forests (Ciais et al 2005, Boisvenue andRunning 2006, Trotsiuk et al 2020). However, it is a question to what extent the changing growing conditions will affect the fitness of different tree-species and regions under future climate conditions. While the ecological niches and current distribution areas of the dominant tree species of Central European forest (P. abies, F. sylvatica and A. alba) largely overlap, their ecological demands differ in numerous aspects (e.g. light demand; Bolte et al 2007, Tinner et al 2013, Pretzsch et al 2014. These differences are expressed in primeval forests by changing species distributions in relation to local ecological gradients, such as soil moisture or the amount of nutrients in soils (Daněk et al 2019). The ability of trees to survive or expand is a result of complex and interacting factors including ecological demands (e.g. Daněk et al 2019), competition traits (e.g. Bošeĺa et al 2013, susceptibility to disturbances (e.g. Kašpar et al 2020), game pressure (e.g. Diaci et al 2011), air pollution (e.g. Kolář et al 2015, Cienciala et al 2017 and climate (e.g. Rossi et al 2016). Due to the large amount of interacting environmental factors, it is necessary to compare samples from multiple environmental gradients or to apply complex numerical models to quantify the importance of individual factors in tree growth performance.
Traditional dendrochronology offers a set of tools for retrospective analyses of tree growth using treering parameters with annual (Büntgen et al 2006) or sub-annual precision (von Arx and Carrer 2014). However, the response of tree growth to climate conditions often has a non-linear character (Rossi et al 2007, 2008, Wilmking et al 2020. Therefore, tools capable of simulating this non-linear response of tree growth to climate conditions are necessary for successful predictions of tree growth and disentangling climate-growth limitations (Jevšenak and Levanič 2016). Climate-driven process-based models of tree ring formation with an intermediate level of complexity are ideal for this purpose (Vaganov et al 2006, Tolwinski-Ward et al 2011. Such process-based models have been used to simulate inter-annual growth variability (Shishov et al 2016), climate-growth limitations (Tumajer et al 2017), years with abrupt growth changes (Tychkov et al 2019) or cambial phenology and kinetics (Tumajer et al 2021). Because of the fine temporal resolutions of simulated growth rates (days), it is possible to use process-based models to quantify between-species differences in ecological demands within one site, or between-site differences along ecological gradients.
In this study, we used tree ring width chronologies of the three main Central European forest species from two of the oldest forest reserves in Central Europe (Boubín approx. 1020 m a.s.l.; and Žofín approx. 780 m a.s.l.) to quantify the spatial-and species-specific differences in climate-growth limitations. Both forests are now undergoing pronounced changes in tree species composition (Šamonil and Vrška 2008), and understanding the perspectives of individual species under current and future climate conditions is of key importance. To address this question, we used a Vaganov-Shashkin process-based model (Vaganov et al 2006) to describe climatic controls on the daily growth rates of the main coexisting species. We assumed that recent climate limitations of trees in the lower elevation Žofín could shed light on future scenarios for the higher elevation Boubín with ongoing climate change. The main goals of our study were to (a) determine the influence of the main climatic factors in observed changes in tree species composition, (b) quantify the species-specific and site-specific differences in growth rates and climatic limitations, and (c) elucidate trends in climate growth limitations over the studied period due to climate change.

Study site
This study was conducted in two of the oldest forest reserves in Central Europe, Žofín (N48 • 40 ′ , E14 • 42 ′ , mean elevation 780 m a.s.l.-further abbreviated as ZOF) and Boubín (N48 • 58 ′ , E13 • 48 ′ , mean elevation 1020 m a.s.l.-further abbreviated as BOU; figure 1). The forest reserves have been strictly protected since 1838 and 1858, respectively. The mean annual temperatures and precipitation are 6.2 • C and 866 mm at ZOF, and 4.9 • C and 1067 mm at BOU, respectively (Tolasz et al 2007). In ZOF, F. sylvatica is currently the dominant species, with an additional high proportion of P. abies, while in BOU the proportion of P. abies and F. sylvatica is approximately equal (figures S1 and S2 (available online at stacks.iop.org/ERL/16/034039/mmedia)). The proportion of A. alba at both sites, based on recent tree inventories, is below 5% (figures S1 and S2).

Studied tree species and their ecological demands
P. abies is an intermediately shade-tolerant species with a shallow root system (Schmid andKazda 2002, Puhe 2003), which makes it highly susceptible to drought stress (Mäkinen et al 2002, 2003, Hartl-Meier et al 2014, Zang et al 2014, Ponocná et al 2016, Tumajer et al 2017, Vitali et al 2017. In addition to direct negative effects of drought stress on P. abies growth, its fitness can also decline due to indirect drought-pronounced factors, such as bark beetle outbreaks (Jakoby et al 2019). P. abies occurs in natural mixed forests from ca. 600 m a.s.l., and its proportion gradually increases with elevation (Pretzsch et al 2014). In contrast to P. abies, F. sylvatica and A. alba are highly shade-tolerant tree species, capable of surviving under closed canopy for a long time (Collet et al 2001, Lichtenthaler et al 2007, Bošeĺa et al 2013. Even though F. sylvatica has competitive dominance on sites with moderate acidity and moisture (Bolte et al 2007), A. alba can effectively overgrow F. sylvatica inside small gaps or under dense canopy cover (Čater and Levanič 2013). The resistance of both species to drought is higher than that of P. abies

Data sampling and processing
Both living and recently dead trees were cored at both study sites, using a randomized schematic design of tree selection (see details in Šamonil et al 2013). Both studied areas were overlaid by a square grid with a plot size of 44.25 m (dimensions derived from the Czech National Inventory, see www.uhul.cz). In ZOF, at each plot we cored six exposed (canopy trees) tree individuals that were closest to the plot center, and within the recent gap occupied by thin young trees we cored nine tree individuals. In BOU we cored one currently exposed tree nearest the plot center plus three additional exposed trees closest to the circles of radii 5, 10, and 15 m around the center. By this approach we obtained a core series of young trees recently exposed in existing younger gaps as well as a core series of old adult canopy trees. One core was extracted from each tree 0.5-1.3 m above the ground. Only those cores with an estimated pith offset of less than 3 cm (Applequist 1958) and without the presence of rot were taken for further analysis; other cores were discarded. In ZOF, 3020 cores were extracted between 2007 and 2010 (Šamonil et al 2013), and of those 2110 cores were used for further analysis (table 1). In BOU, 1378 cores were extracted from 2012 to 2018 (Kašpar et al 2020), 714 of which were used for further analysis (table 1).
Both studied forests are characterized by high soil variability. Hydromorphic soils (gleysols, stagnosols, histosols or fluvisols) are present in 35%-40% of the forest, with the remainder of both forests composed of dryer terrestrial soils (cambisols, podzols, and leptosols). Despite the fact that both are mixed forests (figure S2), P. abies more often occurs on hydromorphic soils than both coexisting tree species (figure S2, table S2). Both sites are also characterized by intermediate soil acidity.
Samples were prepared using standard dendrochronological procedures (Stokes and Smiley 1996). First, all cores were dried at room temperature, then fixed in slices and sanded by fine sandpaper. Subsequently, tree-ring widths were measured on a positioning table using PAST4 software (SCIEM©). Samples taken after 2018 were scanned by a highresolution scanner (EPSON LA2400) and measured using WinDendro software (Régent instru-ments©). Cross-dating of all cores was done in PAST4 (SCIEM©) and later checked using COFECHA (Grissino-Mayer 2001). Both measuring systems yield consistent data (Maes et al 2017).
Our dataset included both young and old trees ( figure S3). Therefore, a regional standardization curve approach was used to estimate growth trends (RCS; figure S4), and then to detrend each raw treering width series (Briffa et al 1992). All steps of this process were handled in R (R Core Team 2019), using the package 'dplR' (Bunn 2008(Bunn , 2010. Finally, we built standard and residual chronologies (figure S5) for each site and species using biweight robust averaging of individual series (six chronologies in total). For each chronology, we calculated the moving interseries correlation (rbar) for 31 year long intervals with 30 year overlaps. Moving correlations were also calculated between all chronologies and between chronologies of individual tree species. The inputs of this model include site latitude (to estimate intra-annual patterns of day length), daily climatic data (mean temperatures and precipitation totals), and residual tree-ring chronology (used as a target to maximize coherence between simulated and observed site chronologies). Residual chronologies were preferred for the VS model calibration, because the model does not consider overwintering effects and autocorrelation in tree-ring width chronologies (Vaganov et al 2006). Daily mean temperatures and precipitation totals were downloaded for the 1950-2018 period as gridded data at 0.1 • resolution from the E-OBS database (Cornes et al 2018). From those data we calculated monthly mean temperatures and precipitation totals. Values of the Palmer drought severity index (PDSI) were downloaded from the CRU database (Barichivich et al 2020).

Tree growth modeling
The model first determines partial growth responses to temperature [Gr T (d,y)] and soil moisture [Gr M (d,y)] for each day 'd' and year 'y' . These are based on a pair of non-linear response functions converting temperatures and moisture (calculated from precipitation by an incorporated hydrological model; Thornthwaite and Mather 1955) into growth responses on a scale from 0 to 1 (figures S6(a)-(c)). The shape of both response functions is determined by four parameters-minimum and maximum temperature/soil moisture enabling growth (T1 and T4/M1 and M4), and minimum and maximum temperature/soil moisture of optimal growth conditions (T2 and T3/M2 and M3). The lower of both partial growth responses (following Liebig's law of minimum) is then multiplied by the daily day length relative to the day length of the summer solstice [Gr S (d)] to produce daily integral growth rates [Gr(d,y)] (figure S6(e)). Annual sums of daily integral growth rates then represent tree-ring width estimates over individual calendar years (figure S6(g)). For a full description of the model algorithm see Vaganov et al (2006).
In addition to simulating annual tree-ring widths, we also extracted key phenological dates from VS model outputs. We defined the growing season onset as the day of the year when the simulated cumulative growth (Gr) exceeded 2.5% of the total cumulative growth of a given year (i.e. when 2.5% of the total simulated tree-ring width was reached). Analogously, the growing season ended on the day of the year when simulated cumulative growth exceeded 97.5% of the cumulative growth of the given year. This approach is less sensitive to outliers and provides more robust results than the default approach implemented in VS (Shishov et al 2016). A similar approach to determining phenological dates has also been used on xylogenesis (e.g. Rathgeber et al 2018) and band dendrometer data (e.g. van der Maaten et al 2018).
The approach mentioned above describes mean simulated growth. However, the growth of trees might significantly differ in exceptionally dry or wet years (Au et al 2020). To assess intra-annual patterns of climatic factors controlling cambial activity, we averaged Gr T (d,y), Gr M (d,y) and Gr(d,y) for each calendar day (d) over 1950-2018 and plotted them together into a single chart. To visualize the intraannual intensity of moisture limitations during dry years we constructed the same chart based only on 10% of the driest and wettest years (He et al 2017). To highlight possible legacy effects of dry years, we then built a similar chart for the years following dry and wet years.
The VS model is controlled by 17 parameters. To iteratively estimate the set of parameters resulting in the highest coherence between simulated and observed chronologies we used VS oscilloscope (Shishov et al 2016). Based on different ecological demands of the studied species, we expected that model parameters would largely differ between species, but marginally between sites (Tumajer et al 2017, He et al 2018. Indeed, we first calibrated the model for each tree species in BOU, because we expected that since it is located closer to ecological limits a stronger climate control on tree growth might be expected (Ponocná et al 2016). After finding optimal sets of parameters for BOU, we applied them to ZOF and modified parameter values to find their optima for this site (table S1). The coherence between simulated and observed chronologies was quantified using the Pearson correlation coefficient, RMSE and synchronicity (Shishov et al 2016).

Statistical analysis
To quantify the difference between the real (observed) growth of all studied tree species, the mean treering widths of each tree species during 1950-2018 at both sites were bootstrapped to obtain 95% confidence intervals of the tree growth. In addition, we evaluated a linear climate-growth relationship computing the correlation of monthly climatic data from the previous year March to the current year October (mean temperatures, precipitation sums and drought expressed as PDSI index; Barichivich et al 2020) with standard chronology. Growth trends of individual trees were evaluated using the Mann-Kendal test.
From the results of the VS model we bootstrapped the mean growing season length, the mean number of days when tree growth was climatically unlimited (Gr T = Gr M = 1), and the number of days limited by temperature (Gr M > Gr T > 0) and moisture (Gr M > Gr T > 0). Because three relatively different values were calculated, the set of all values in each exact bootstrap always came from the same years. For subsequent comparisons we used values of the 0.025th, 0.5th and 0.975th percentiles, again representing 95% confidence intervals and medians. The trends in the dates of the beginning and end of the growing season, the growing season duration, and in the cumulative growth of individual tree species at both sites were evaluated using linear regression.

Observed growth trends and linear climate-growth relationships
Considering all tree-ring data (covering approximately from 1800 to the present), P. abies had significantly wider tree-rings than the other two species. A. alba grew faster than F. sylvatica in BOU but not in ZOF (table 1). Since the 1950s, P. abies expressed significantly higher growth rates than both other tree species in ZOF, while F. sylvatica had significantly higher growth rates than A. alba in BOU. Since 1990, differences in growth rates among species have been insignificant regardless of site (table 1). This was confirmed by an analysis of growth trends showing mostly positive slopes in tree-ring widths for all the tree species in the calibration period. Specifically, 49.5% and 54.6% of P. abies, 59.4% and 64.7% of F. sylvatica and 67.9% and 75% of A. alba individuals showed positive trend in BOU and ŽOF respectively since 1950. Considering entire length of all series, P. abies showed mostly negative trends, F. sylvatica mostly positive trends and A. alba showed positive growth trends in ZOF and negative trends in BOU.
Site standard chronologies showed significant positive correlations with growing season temperatures, and in the case of F. sylvatica with summer temperatures ( figure S7). Overall, calculated correlations were generally higher in BOU than in ZOF ( figure  S4). For P. abies the strongest correlations occurred with March temperature (figures S7(a) and (d)). Site chronologies of F. sylvatica correlated strongly with temperature in July and August in ZOF ( figure S7(b)), and in August and September in BOU ( figure S7(e)). The growth of A. alba was significantly correlated with temperatures throughout the entire growing season in BOU ( figure S7(c)), while in ZOF correlations were significant only for April temperatures ( figure S7(f)). In contrast to temperatures, correlations with precipitations were mostly insignificant (figure S7).

Model parameterization and simulated tree growth dynamics
The highest sensitivity of simulated chronologies was found for the variation in parameters of minimal temperature for growth (T1) and temperature defining the minimal threshold for optimal growth (T2), followed by cumulative temperatures for the beginning of growth (tbeg), and the minimum soil moisture and lower range of optimal soil moisture conditions (M1 and M2). The final set of parameters for the VS model are given in table S1. Correlations between observed and simulated chronologies were generally higher for BOU (figures 2(a)-(c)) than for ZOF (figures 2(d)-(f)), except for P. abies. Statistically significant correlations were observed between both curves (VS model simulations and observed residual chronologies) for all tree species and sites, except F. sylvatica in ZOF (figure 2). The high-frequency variability of simulated chronologies was visually well-synchronized to the variability in observed standard and residual chronologies ( figure S5).
The simulated growth rates in the peak growing season (June and July) reflected variability in temperature conditions, because soil moisture dropped below optimal levels throughout the summer. The highest growth rates (Gr > 0.6) occurred on average from the beginning of May (P. abies) or the beginning of June (F. sylvatica) till the end of September, overlapping the period with the lowest temperature limitation (figures S8 and S9). Mean moisture limitation was relatively low at both sites (Gr M > 0.95 in most of the year; figure S8). Nevertheless, a weak moisture limitation was apparent for some species/sites from mid-June to mid-October, with a peak in August, when Gr M occasionally fell below 0.8 ( figure  S8). The overall moisture limitation was stronger in ZOF than in BOU. The most sensitive tree species to drought was P. abies (figure S8). Moisture limitations for A. alba and F. sylvatica were considerably lower, but there were also between site differences. In BOU, F. sylvatica was basically unlimited by moisture, while in ZOF, A. alba was less limited by moisture during 10% of the driest years (figure S10). The total annual growth (Gr) of P. abies was reduced on average about 11% in BOU and about 5% in ZOF. By contrast, F. sylvatica at both sites and A. alba in ZOF showed opposite results in dry years, with higher growth rates than average (from 3% to 5%). Higher growth in the 10% of the wettest years was observed for P. abies and F. sylvatica in ZOF (2.7% and 1.8%, respectively). In BOU, the growth of F. sylvatica and A. alba was reduced during the wettest years (about 5.8% and 5.4%, respectively). Long periods of drought also affected P. abies in the following growing season, with the growth reduced by 5% in BOU and 4% in ZOF (figure S11). Also, F. sylvatica in ZOF was still affected by the drought in the subsequent year, when simulated growth was on average about 8% lower. In the years following the 10% of the wettest years, only P. abies at both sites had higher simulated growth (about 2% in BOU and 4% in ZOF).
The simulated growing season in BOU was shorter than the growing season in ZOF (with a later start and earlier end (figures 3 and S8)). The mean simulated growing season duration was longest for P. abies (156.9 d in BOU and 166.2 d in ZOF). The F. sylvatica growing season lasted 119.4 d in BOU and 122.1 d in ZOF, and for A. alba 139.5 d in BOU and 144.2 d in ZOF. Statistically significant differences in mean growing season duration were documented between P. abies and F. sylvatica and between A. alba and F. sylvatica. In contrast, growing season duration did not significantly differ between P. abies and A. alba. Between-site differences were statistically significant for all species except F. sylvatica.
The linear regression showed an increasing duration of the growing season at both sites. The slopes of the linear trend in growing season duration spanned from 1.2 to 1.5 d decade −1 , depending on the species and site (figure 3(a)). For BOU, the extension of the growing season was 1.5, 1.2 and 1.4 d decade −1 for P. abies, F. sylvatica and A. alba, respectively. In ZOF the extension of the growing season was 1.5 d decade −1 (A. alba and P. abies) and 1.4 d decade −1 (F. sylvatica). This trend was due to both the earlier onset and later cessation of the growing season (figure S12).

Simulated growth trends and climate-growth limitations
At both sites, P. abies has the highest simulated growth rates from 1950 to 2018 (i.e. the long-term mean of the annual sum of Gr), followed by A. alba and F. sylvatica. The differences for all pairs of species were statistically significant. The simulated growth rates of P. abies and A. alba were significantly higher in ZOF, but for F. sylvatica this difference between sites was statistically insignificant. The highest simulated mean monthly growth rates were for P. abies (0.94 at both sites), followed by A. alba (0.86 and 0.89) and the lowest for F. sylvatica (0. 80 and 0.83 in BOU and ZOF,respectively). The highest growth rates were always observed in July (figures 4(a) and S8), except for P. abies that had the highest growth rates already in June. All the species showed lower growth rates between the 1960s and 1980s, followed by increasing growth rates in the more recent decades (figures S4 and S5).
Simulated data for all tree species in BOU showed higher increase in annual growth rates than in ZOF ( figure 3(b)). The increase in annual growth rates of F. sylvatica in BOU was almost triple compared to ZOF (+0.30% and +0.11% per year of mean simulated annual growth rates during entire simulated period). The slopes of increases in growth rates of A. alba were +0.19% in BOU and +0.13% in ZOF and of P. abies +0.08% in BOU and +0.02% in ZOF.
There was lower proportion of growing season days with temperature limitation in ZOF compared to BOU (figure 4). This difference was statistically significant for F. sylvatica (on average 71.0% days in BOU and 51.0% days in ZOF; figures 4(c) and (f)) and A. alba (on average 56.7% days in BOU and 51.2% days in ZOF; figures 4(d) and (g)), but insignificant for P. abies (on average 28.1% days in BOU and 25.9% days in ZOF;figures 4(b) and (e)). On the other hand, the growth of trees in ZOF was more limited by soil moisture (figure 4). Differences in the number of moisture limited days were statistically significant for P. abies (on average 7.5% days in BOU and 23.6% days in ZOF; figures 4(b) and (e)) and F. sylvatica (no limitation in BOU and on average 18.9% days in ZOF; figures 4(c) and (f)), but insignificant for A. alba (3.1% in BOU and 0.3% in ZOF;figures 4(d) and (g)). Therefore, the numbers of days with optimal growth conditions (days without any climate limitation; figure 4) were significantly higher in BOU for P. abies (on average 64.5% days in BOU and 50.7% days in ZOF;figures 4(b) and (e)), significantly lower in BOU for A. alba (on average 40.2% days in BOU and 48.5% days in ZOF; figures 4(d) and (g)), and insignificantly different for F. sylvatica (on average about 29.0% days in BOU and 30.0% days in ZOF; figures 4(c) and (f)).

Discussion
Throughout the Holocene, changes in the species composition of Central European mixed forests were driven by climate oscillations (e.g. Tinner andLotter 2006, Fletcher et al 2010). The results of our study clearly indicate that climate-growth interactions have also markedly determined forest dynamics during the last century, when unprecedented climate change has been observed (Luterbacher et al 2004). Although an increasing synchronicity in growth across the globe suggests the existence and increasing importance of a common growth-limiting climatic factor (Manzanedo et al 2020), our simulated prominent between-species differences in intra-annual climatic limitations highlight the contrasting responses of the main coexisting tree species of Central European mountain forests to climate change. This will have a direct impact on Central European old-growth forests, leading to significant changes in composition and structure. For instance, the simulations revealed significant moisture limitation for P. abies, while F. sylvatica proved to be drought-resistant but most limited by temperatures. Such contrasting responses of tree growth to climate on an intra-annual scale resulted in large differences in the simulated response of total annual growth to climate change. Indeed, highly positive growth trends were simulated for F. sylvatica in the higher-elevated BOU but marginal trends were found for P. abies at both sites. In addition, our results correspond to changes seen in the forest structure over the last century, and especially to processes occurring during the most recent decades. Accordingly, we presume that changes in climategrowth limitations are the driving factor of ongoing forest structure alterations in Central European forests (Fuhrer et al 2006) and will continue shaping them in the near future ( Buras and Menzel 2019).

Model parameterization and simulated tree growth dynamics
The spatial distribution of dominant climatic controls on tree growth and simulated VS parameters are in accordance with findings of studies focused on wood formation in both mountain and lowland forests (e.g. Giagli et al 2015, Treml et al 2015, Rossi et al 2016, Antonucci et al 2019. For instance, the critical temperature for the growth reactivation (parameter T1) was the most sensitive parameter for all investigated tree-species in the higher elevation BOU. On the contrary, minimal moisture limitation was more important in the lower elevation ZOF, mainly for P. abies. This is in line with results of previous studies Simulated chronologies were well correlated with observed residual chronologies for both highfrequency variability and decadal trends (figures 2 and S5). Overall, the reliability of the simulations was higher in BOU (except for P. abies) due to stronger temperature controls on tree growth near species ecological distribution limits (Breitenmoser et al 2014). This is also seen in the highest coherence between simulated and observed chronologies the relatively colder periods of the 1970s and 1980s, or in case of P. abies in ZOF, in the most recent period of drought stress (figure 2). The lower correlation between simulated and observed growth chronologies in ZOF is probably a consequence of a mixed climatic signal typical for intermediate elevation stands (Ponocná et al 2016, Kolář et al 2017, Tumajer et al 2017. By contrast, the growth of P. abies in ZOF is largely limited by soil moisture, indicating that P. abies is closer to its ecological distribution limit here and not in BOU. The overall poorer performance of the model for F. sylvatica can be explained by the lower environmental sensitivity of this species (Bošeĺa et al 2013), the significant effect of masting years distorting the climatic signal of site chronologies (Hacket-Pain et al 2015), and by the fact that the VS model was not originally designed for deciduous tree species (Vaganov et al 2006). The lower simulated growth of A. alba could partly be explained by the inclusion of a period with significant air pollution into the calibration data.
Despite the fact that the VS model successfully reproduced high-frequency variability and intermediate-frequency trends, the simulated growth intensity did not fully correspond to actual tree-ring widths. P. abies had the highest simulated growth regardless of site, followed by A. alba and F. sylvatica. This corresponds to the situation observed since the mid-19th century but does reflect growth rates measured during the model calibration period from 1950 to 2018 (table 1). After 1990, observed differences in growth rates among species were insignificant.

The consequences of different species-specific climate growth trends in simulated growth
Despite the observed effects of moisture limitation on growth at both sites, the increase in growing season temperatures had a dominant positive effect on the annual growth of all studied tree species. Indeed, simulated tree-ring widths of F. sylvatica and A. alba were significantly narrower in wet years due to lower temperatures associated with rainy weather. This highlights the complexity of climate change effects on forest ecosystems (Trotsiuk et al 2020), which is significantly modulated by species and site conditions. One such effect is clear when evaluating trends in simulated annual growth, with trees in the higherelevation BOU increasing more than in ZOF. On the one hand, increasing temperatures lead to an extension of the growing season ( Boisvenue and Running 2006), but on the other hand lead to a higher level of summer drought stress due to increased evapotranspiration (Briffa et al 2009). The consequences of these processes might differ between species, as shown by contrasting patterns of climate-growth interactions (Kolář et al 2017, Vitali et al 2017, Bošel'a et al 2019. In accordance with the between-species differences in drought resistance reported by these studies, our simulation showed a much lower increase in the annual growth of P. abies over the last seven decades compared to both coexisting tree species. As expected, P. abies, a common treeline tree species (Kašpar and Treml 2018), showed the lowest temperature limitation (figures 4(a) and (d)). On the other hand, it experienced the highest moisture sensitivity, as had already been demonstrated across large latitudinal and altitudinal gradients (e.g. Mäkinen et al 2002, 2003, Giagli et al 2016, Ponocná et al 2016. The reduction of P. abies simulated growth was apparent not only in the dry year, but through legacy effects also in the subsequent year. The growth reduction in the following year accounted about 5% of long-term mean tree-ring width. This might be a consequence of shallow root system of P. abies (Dobbertin 2002, Schmid and Kazda 2002) and low hydraulic resistance to negative water potential, making it prone to drought-induced embolism (Charra-Vaskou et al 2012). Therefore, it is currently being competitively displaced at both study sites by other species into sites with relatively hydromorphic soils. The results of our simulations were obtained based on a large proportion of samples originating from trees on hydromorphic soils. Approximately half of the P. abies used in this study already showed negative growth trends in simulated period (and more than a half when considering entire tree-ring series). Indeed, even more extensive changes of tree species composition may be expected to the detriment of P. abies, especially on sites with dry soils. In addition, negative trends in the growth of P. abies might be accelerated due to drought-pronounced indirect stress factors, such as bark beetle outbreaks (Jakoby et al 2019).
Of the studied tree species, F. sylvatica showed the highest temperature limitation, which is because both sites are relatively close to the altitudinal limit of its distribution (Breitenmoser et al 2014). However, simulated growth also revealed a relatively severe moisture limitation of F. sylvatica in ZOF (figures 4(b) and (e)). The growth of F. sylvatica was even more affected in the year following the dry season than in the case of P. abies. This is in accordance with Pretzsch et al (2013), who suggests a limited stomatal control of transpiration during days with negative water potential in F. sylvatica. An associated high intensity of vessel cavitation (Hacke and Sauter 1995) might be a legacy effect on growth in subsequent growing seasons. The delayed response to drought compared to spruce may also be due to a deeper root system. F. sylvatica relies on soil moisture from deeper layers of the soil, which resaturate slowly after strong drought spells due to the delayed infiltration of rain water. Cavitation in broadleaves can also be induced by spring frosts at high altitudes ( Hacke and Sauter 1996). Increasing frequencies of spring frost events might occur as a side effect of an earlier start of the growing season (Schwartz et al 2006). Nevertheless, simulated growth rates show the highest increasing trends for this tree species, mainly due to decreasing temperature limitation (Di Filippo et al 2007). Therefore, at our study sites a positive effect of growing season extension (Prislan et al 2019) and lower temperature limitation will likely outweigh the negative effects of drought stress for F. sylvatica. In addition, the frequency of years with high precipitation and low summer temperatures, with negative effects on growth, might be expected to decline due to climate change as well.
In contrast, A. alba showed the highest drought resistance and an intermediate temperature limitation (figures 4(c) and (f)). This is likely related to its ability to regulate transpiration and the relatively high minimum water potential (even lower than F. sylvatica ;Nourtier et al 2014). Considering mainly its resistance to drought, our results showed a sufficient ability of A. alba to resist stress related to ongoing climate change and increasing drought, making it a well-adapted tree species for future Central European mountain forests (Fuhrer et al 2006, Büntgen et al 2014, Gazol et al 2015. Despite this, however, A. alba has the lowest proportion of the studied species in current forests due to pollution-related dieback culminating in the 1980s (Elling et al 2009, Bošel'a et al 2016) and preferential game browsing preventing its regeneration (Diaci et al 2011). However, due its climatic suitability and endangerment by human activities, forestry management procedures should be aimed at supporting its regeneration in Central European managed forests (e.g. using fencing).

Effects of climate growth limitations on future forest composition
Changes in the intensity of moisture and temperature limitations generally follow elevation gradients (Tumajer et al 2017). Indeed, simulated temperature limitations were higher in BOU, though the number of days with optimal growing conditions was lower in ZOF due to higher moisture limitations. Simulated growth rates increased more in BOU than in ZOF (figure 3), suggesting a significant increase in the productivity of the higher elevation forest (Tumajer et al 2017, Trotsiuk et al 2020. The only discrepancy in this trend was the simulated growth of A. alba, which showed a stronger growth limitation at higher elevations due to soil water over-saturation and subsequent root hypoxia. Climate change together with large scale disturbances will affect the regeneration of the climatically most-limited P. abies. Its regeneration on nonhydromorphic soils will be restricted by drought and possibly also by the competitive advantage of drought-resistant F. sylvatica and A. alba (table S2). However, growing under such conditions forces the development of a shallow root system (Puhe 2003), making trees prone to future windstorms. A. alba has a similar pattern of climategrowth limitations as F. sylvatica (Zang et al 2014), but is better adapted to drought conditions (Zang et al 2014, Hoffmann et al 2018 and has recently shown high growth rates (Bošeĺa et al 2014, Büntgen et al 2014. Nevertheless, its limited proportion in the composition of current forests will prevent it becoming a dominant tree species in the near future ( Buras and Menzel 2019). Despite this, in the context of our simulations as well as actual observations of the forest structure at both sites (Šebková et al 2011(Šebková et al , Král et al 2014, the eventual expansion and dominance of F. sylvatica might be expected (Janík et al 2016, Buras and Menzel 2019). In particular, this species might benefit from its high adaptability to future climate conditions (Di Filippo et al 2007, Hoffmann et al 2018, Bošel'a et al 2019, as is reflected by the larger increase in annual growth compared with P. abies and mostly positive growth trends of individual trees. This trend may also largely explain the increase in the population of F. sylvatica across Central Europe. Anticipated species composition changes in Central European mountain forests might have significant consequences for other environmental processes in forest ecosystems, e.g. disturbances (Bobek et al 2018(Bobek et al , 2019 or soil development (Daněk et al 2019). If a locality approaches the ecological limit of a tree species (e.g. moisture limitation), the future spatial distribution of this species will better reflect the local habitat conditions. Under such a scenario, geomorphological conditions may become the dominant factor controlling the tree species distribution (Dyderski and Pawlik 2020).

Conclusions
Contrasting patterns of climatic limitations exist among the main tree species of Central European mixed mountain forests. We found species-specific differences in the response to ongoing climate change, which are modulated by local climate conditions over the elevation gradient. Forests at higher elevations are expected to become more productive due to an increase in temperature, leading to increasing summer growth rates and a longer growing season, while the intensity of moisture limitations is expected to remain marginal. Trees at lower elevations, however, will experience insufficient moisture and a limitation of summer growth rates due to drought. This moisture limitation is expected to be the driving factor in species composition changes. Specifically, a decline of drought-sensitive P. abies in lower elevations at the expense of tree species more suitable to future climate conditions can be expected. Our results suggest that A. alba is well adapted for the future climate, but expansion in natural forests will probably be slower because of its limited proportion in current forest stands. In contrast, the expansion and dominance of F. sylvatica can be expected due to the large number of juveniles, its sufficient level of drought resistance, projected reductions of temperature growth limitations, and the recent ongoing positive trend in growth rate.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.