Climate Variability May Delay Post-Fire Recovery of Boreal Forest in Southern Siberia, Russia

: Prolonged dry periods and increased temperatures that result from anthropogenic climate change have been shown to increase the frequency and severity of wildﬁres in the boreal region. There is growing evidence that such changes in ﬁre regime can reduce forest resilience and drive shifts in post-ﬁre plant successional trajectories. The response of post-ﬁre vegetation communities to climate variability is under-studied, despite being a critical phase determining the ultimate successional conclusion. This study investigated the responses of post-ﬁre recruited species to climate change and inter-annual variability at 16 study sites that experienced high-severity ﬁre events, mostly in early 2000, within the Scots pine forest-steppe zone of southeastern Siberia, Russia. These sites were originally dominated by Scots pine, and by 2018, they were recruited by different successional species. Additionally, three mature Scots pine stands were included for comparison. A Bayesian Additive Regression Trees (BART) approach was used to model the relationship between Landsat-derived Normalized Difference Vegetation Index (NDVI) time series, temperature and precipitation in the 15 years after a stand-replacing ﬁre. Using the resulting BART models, together with six projected climate scenarios with increased temperature and enhanced inner-annual precipitation variability, we simulated NDVI at 5-year intervals for 15 years post-ﬁre. Our results show that the BART models performed well, with in-sample Pseudo-R 2 varying from 0.49 to 0.95 for ﬁre-disturbed sites. Increased temperature enhanced greenness across all sites and across all three time periods since ﬁres, exhibiting a positive feedback in a warming environment. Repeatedly dry spring periods reduced NDVI at all the sites and wetter summer periods following such dry springs could not compensate for this, indicating that a prolonged dry spring has a strong impact consistently over the entire early developmental stages from the initial 5 years to 15 years post-ﬁre. Further, young forests showed higher climate sensitivity compared to the mature forest, irrespective of species and projected climatic conditions. Our ﬁndings suggest that a dry spring not only increases ﬁre risk, but also delays recovery of boreal forests in southern Siberia. It also highlights the importance of changing rainfall seasonality as well as total rainfall in a changing climate for post-ﬁre recovery of forest. Contributions: Data analysis,


Introduction
In Siberia, the length of fire season, the frequency of moderate-to high-severity fires, as well as the magnitude of total burnt areas have all increased markedly over recent decades [1][2][3][4]. Some of the largest increases in fire frequency and severity have been observed close to the southern boundary of the boreal forests (i.e., from 50 • to 55 • N), with regions such as the vast forested areas in the Zabaikal region of southern Siberia burning more than once during the last 20 years [5].
Site-level impacts of fire disturbance on forest structure and composition depend upon fire characteristics in combination with factors such as drainage/soil moisture and fuel load [6,7]. Species such as Scots pine (Pinus sylvestris L.) and Larch Larix spp. can recover from low-severity fire disturbance and may even be promoted by low-to moderate-severity fires in Siberian forests [8][9][10] as this results in an alteration in soil conditions post-surfacefires, which tends to favor pine over other species (pine requiring a mineralized soil for successful regeneration and recruitment [11]). However, high-severity stand-replacing fires can consume all aboveground biomass down to the mineral soil layer (apart from some large woody biomass), which allows various post-fire recovery trajectories with different forest structure and composition, and thus with different consequences for carbon stocks.
In many parts of the boreal zone, there is strong evidence that evergreen conifers are being replaced by deciduous species in taiga forests, as differences in post-fire recolonization strategies exist among these species [12][13][14]. This is likely the case because intensified fires cannot only destroy the aerial seed banks that conifers rely upon for rapid postfire establishment, but also favor faster-maturing deciduous species [15,16]. In the forest-steppe zone of southern Siberia, approaching the southern latitudinal limit of boreal forest, conifer trees can also be replaced by steppe vegetation [5,17]. Recovery, regeneration and ultimately maintained presence of taiga forest following fire disturbance in the southern boreal zone is, therefore, a critical consequence of a warming climate.
A recent study highlighted an increased vulnerability of forests to climate-induced regeneration failures following wildfires [18]. This is particularly the case for post-fire recruited species in their early developmental stages when ecosystem environmental conditions are harsh. Previous studies suggest that if tree species fail to establish in the first five to ten years post-fire, it is highly likely to result in a permanent forest loss [16,19,20]. Young forest stands are more vulnerable to climate variability than mature forests, largely because of immature root systems that can further induce hydraulic failure over prolonged drought periods [21,22].
In boreal forests across the northern mid-to high-latitudes, including Siberian forests, increased air temperature, more frequent prolonged dry periods and climatic extremes have already been recorded or are predicted based upon empirical observations [23][24][25]. Warmer and drier conditions particularly facilitate fire disturbance in coniferous forests and the boreal biome [26]. Furthermore, intensified fire activity might even amplify interannual climate variability during both spring and summer by inducing changes in such parameters as surface albedo and the hydrological cycle [27][28][29]. This can result in greater difficulties in predicting vegetation recovery trajectories post-fire due to the changing current and future temperature and precipitation regimes, as post-fire recruitment success appears to respond positively to the presence of winter snow cover in the first several years after a burn [17].
Temperature has increased over the entirety of Siberia over recent decades, also notably in 2020 (January to May 2020 being >5 • C above the 1951-1980 average [30]), and it is likely to be more than 2 • C above reference under all IPCC CMIP scenarios over the remainder of the twenty-first century [31]. Unlike rising temperature, precipitation regimes are highly variable, dependent upon geographical location. For example, in the forest-tundra ecotone of northeastern Siberia, annual precipitation has changed little over recent decades [32]. This is also the case for some regions of southern Siberia over the period from 1982 to 2017, according to TerraClimate data [33]. However, in some parts of southern Siberia, an increase in annual precipitation (of ca. 14 mm per decade) has been recorded over the period 1935-2013 [34]. In addition to changes in the annual precipitation budget, studies have also revealed that summer precipitation input has varied greatly in Siberia, suggesting a strong inter-annual variability in precipitation regime [10,35,36]. Future estimates of precipitation have even greater uncertainty, and one possible change presented in the Coupled Model Intercomparison Project phase 5 (CMIP5) modeled precipitation is that colder months are likely to experience larger fractional precipitation increases than warmer months, which can result in greater variation in seasonal precipitation [36]. Understanding the responses of post-fire recruited species to climate change and intra-and inter-annual variability is therefore crucial for assessing the interactive effects of fire and climate upon forest development. However, to date, such studies are generally lacking.
The aim of the present study was to determine the post-fire responses of various vegetation communities, which were following different successional pathways, to altered climatic conditions in a forest-steppe zone of southeastern Siberia. To achieve this, a combination of remote sensing data, field observational data and machine learning algorithms under a Bayesian framework were employed.
Landsat-derived Normalized Difference Vegetation Index (NDVI) time series data were used to characterize the condition of vegetation communities pre-and post-fire. NDVI, as a proxy for vegetation productivity and coverage [37], has been widely used in studies relating to fire disturbance and post-fire recovery [38][39][40][41][42][43]. For example, Shvetsov et al. [39] aimed to evaluate the performance of SWIR-based satellite metrics (e.g., the normalized burn ratio time series) in assessing post-fire vegetation dynamics in the Zabaikal region, and to model the short-term vegetation recovery. Their research findings showed that fire frequency and severity, as well as surface temperature anomalies, are affecting vegetation recovery. In addition, they encouraged future studies to investigate the effects of moisture and precipitation on post-fire vegetation recovery in this region. Such a derived NDVI metric has been proven to be a good indicator of vegetation productivity in Siberian forests [44]. We selected a range of sites that have recently been burnt by a stand-replacing fire mostly in the early 2000s. All these sites were dominated by Scots pine, but now are following different trajectories post-fire with different recruited species components. The approximately 15-to 20-year recovery period since the last stand-replacing fire also matches the fire return interval of the study region, as reported by Kukavskaya et al. [5]. Given the timeframes of the majority of fire events (most in the early 2000s), Landsat data were our primary data sources. Despite the fact that a daily MODIS-derived NDVI record (MCD43A4) has been freely available since 2000, the MODIS spatial resolution (500 m) is much coarser compared to the Landsat images (30 m). Further, we proposed a range of potential climatic scenarios for this study, which include increases in air temperature and annual total precipitation input, as well as shifts in seasonal precipitation distribution.
Our first hypothesis was that forest age would modulate the sensitivity to climate variability in southern Siberia. It has been suggested that some of the recruited species may resist certain levels of unfavorable climatic conditions in the forest-steppe zones of Siberia [45]. For example, Scots pine is considered to be of exceptionally high tolerance to water stress compared to other tree and grass species present in the region [46][47][48]. Therefore, our second hypothesis was that sites following a coniferous recovery pathway would respond differently to climate variability than those following a mixed coniferous/deciduous and/or grassland pathways. Thirdly, we hypothesized that total annual precipitation input as well as shifts in seasonal precipitation distribution during growing seasons (April-September) would influence recovery of regenerating sites. In other words, precipitation distribution (both timing and seasonal amounts) is important-not just the total annual input in this region. Our study sites in the region were on west-and southwest-facing slopes locate tween 537 and 1035 m a.s.l., receiving generally less than 400 mm precipitation ann (see Figure 2; data from one of the weather stations in the area, ID: GHCND:RSM0003 GPS location: 51.35° N, 110.47° E, accessed from https://www.noaa.gov/, (accessed December 2019)). Meteorological data for the region are sparse, and this station (pin angle in Figure 1) was chosen as central to our spread of sites (for further details, als Supplementary Figure S1, showing median, maximum and minimum annual preci tion over the period 1980-2018 for the six stations available in the region). Our study sites in the region were on west-and southwest-facing slopes located between 537 and 1035 m a.s.l., receiving generally less than 400 mm precipitation annually (see Figure 2; data from one of the weather stations in the area, ID: GHCND:RSM00030844, GPS location: 51.35 • N, 110.47 • E, accessed from https://www.noaa.gov/, (accessed on: 1 December 2019)). Meteorological data for the region are sparse, and this station (pink triangle in Figure 1) was chosen as central to our spread of sites (for further details, also see Supplementary Figure S1, showing median, maximum and minimum annual precipitation over the period 1980-2018 for the six stations available in the region). ) and proposed climatic scenarios which included shifts in precipitation distribution (tri and square) and changes in total precipitation input (spring and summer seasons spanning April to September; plus sign and square with cross). Shaded area in grey shows the maximum minimum of monthly total precipitation (all forms, including snow) over the observational pe Figure 3 outlines the integrated processes followed within this study, details of of which can be found in the following sections. In the forest-steppe of the Zabaikal and Republic of Buryatia regions, the tree v tation is dominated by Scots pine (Pinus sylvestris L.), with sporadic small patches of a (Populus spp.), birch (Betula pendula Roth) and larch (Larix gmelinii (Rupr.) Rupr. ex K and proposed climatic scenarios which included shifts in precipitation distribution (triangle and square) and changes in total precipitation input (spring and summer seasons spanning from April to September; plus sign and square with cross). Shaded area in grey shows the maximum and minimum of monthly total precipitation (all forms, including snow) over the observational period. Figure 3 outlines the integrated processes followed within this study, details of each of which can be found in the following sections. ) and proposed climatic scenarios which included shifts in precipitation distribution (triangle and square) and changes in total precipitation input (spring and summer seasons spanning from April to September; plus sign and square with cross). Shaded area in grey shows the maximum and minimum of monthly total precipitation (all forms, including snow) over the observational period. Figure 3 outlines the integrated processes followed within this study, details of each of which can be found in the following sections. In the forest-steppe of the Zabaikal and Republic of Buryatia regions, the tree vegetation is dominated by Scots pine (Pinus sylvestris L.), with sporadic small patches of aspen (Populus spp.), birch (Betula pendula Roth) and larch (Larix gmelinii (Rupr.) Rupr. ex Kuzen and L. sibirica Ledeb.) [5]. The understory is mainly rhododendron (Rhododendron spp.), rose (Rosa spp.), diamond willow (Salix bebbiana Sarg.) and spirea (Spiraea spp.). Leaf area In the forest-steppe of the Zabaikal and Republic of Buryatia regions, the tree vegetation is dominated by Scots pine (Pinus sylvestris L.), with sporadic small patches of aspen Remote Sens. 2021, 13, 2247 6 of 20 (Populus spp.), birch (Betula pendula Roth) and larch (Larix gmelinii (Rupr.) Rupr. ex Kuzen and L. sibirica Ledeb.) [5]. The understory is mainly rhododendron (Rhododendron spp.), rose (Rosa spp.), diamond willow (Salix bebbiana Sarg.) and spirea (Spiraea spp.). Leaf area index (LAI) varied between 0.29 to 2.28 at our study sites, calculated from digital photography (following Macfarlane et al. [50]). Soils of the top 30 cm depth were collected during 2018 field sampling and further characterized as light, medium or sandy loam. Soil bulk density (calculated as the dry weight of soil divided by its volume) at 30 cm depth ranged from 1.16 to 1.57 g cm −3 (Table 1). Soil pH and EC were measured in a 1:5 soil:deionized water suspension, and soils were non-saline (EC 1:5 < 0.4 dS m −1 ) and had pH 1:5 ranging from 6.0 to 7.1. Table 1. Site characteristics in a boreal forest-steppe zone of southern Siberia (data collected in 2018) for sites that were dominated by grass (G), Scots pine (SP), a mixture of species (Mix) and mature Scots pine (C), respectively. The dash sign (-) indicates that data are not available (in the case of site G3, as a result of extreme weather preventing soil sampling).

Site Selection
A high-severity fire can change an original pine-dominated stand into a mixed species forest stand that comprises one or more of the following species: aspen, birch, larch and pine, or occasionally into a grass-dominant ecosystem. Therefore, we determined that our sites had been burnt only once in the past 30 years and were already following various successional pathways to the present day. Thus, we included sites that were either dominated by grass ("G"), Scots pine ("SP") or a mixture of species ("Mix"). These fire-disturbed sites represented major successional pathways.
In this study, a Scots pine-dominated site was delineated from a mixed stand group by the absolute and relative quantity of pine present in the stand. In addition, we included three mature Scots pine stands (referred to as "C") that have not been burnt since 1987 as a comparator. Overall, we had 19 sites comprising 7 grass-dominated sites, 5 Scots pine recruitment sites, 4 mixed recruitment sites and 3 mature pine sites.

Fire History Data
The fire history of each site was verified by using records from the local forestry authority in southern Siberia and from remotely sensed images.
To briefly illustrate the latter assessment approach, we usedCMIP5 data, in combination with a time lapse movie for each study site showing R-G-B scenes as well as near-infrared (NIR)-R-G and shortwave infrared (SIR)-NIR-R false-color composites, with all possible Landsat satellite images (Landsat-5 TM, Landsat-7 ETM+ and Landsat 8 OLI) collected from 1987 to August 2019 (see example shown in Figure 4). Specifically, Landsat-5 TM images from 1987 to 2011, Landsat-7 ETM+ images from 1999 to 2019 and Landsat 8 OLI images from 2013 to 2019 were used for the time lapse movie at each site. A minimum of two Landsat images per month were collected. We then identified fire events at each site manually by two experienced researchers independently viewing the entirety of each frame of the 'animation footage'. The assessment reports submitted by the two experienced independent researchers were then compared to determine the final fire histories. Due to differences in the Landsat series of satellites' timelines, a minimum of two Landsat images per month were collected. Therefore, over the past 32 years (from 1987 to 2019), at least 768 images (=32 × 12 × 2) were used to assess stand-replacing fire events for each site.

Fire History Data
The fire history of each site was verified by using records from the local forestry authority in southern Siberia and from remotely sensed images.
The historical fire data from the local forestry authority were initially collected by local foresters and firefighters and then submitted to the regional Federal Forest Service. Such information can be viewed through their official websites (the Zabaikal regional Federal Forest Service office, Chita, Russia: http://лесслужба.забайкальскийкрай.рф/ (accessed on: 1 December 2019), and Buryatia regional Federal Forest Service office, Ulan-Ude, Russia: https://egov-buryatia.ru/ralh/ (accessed on: 1 December 2019)).
To briefly illustrate the latter assessment approach, we usedCMIP5 data, in combination with a time lapse movie for each study site showing R-G-B scenes as well as nearinfrared (NIR)-R-G and shortwave infrared (SIR)-NIR-R false-color composites, with all possible Landsat satellite images (Landsat-5 TM, Landsat-7 ETM+ and Landsat 8 OLI) collected from 1987 to August 2019 (see example shown in Figure 4). Specifically, Landsat-5 TM images from 1987 to 2011, Landsat-7 ETM+ images from 1999 to 2019 and Landsat 8 OLI images from 2013 to 2019 were used for the time lapse movie at each site. A minimum of two Landsat images per month were collected. We then identified fire events at each site manually by two experienced researchers independently viewing the entirety of each frame of the 'animation footage'. The assessment reports submitted by the two experienced independent researchers were then compared to determine the final fire histories. Due to differences in the Landsat series of satellites' timelines, a minimum of two Landsat images per month were collected. Therefore, over the past 32 years (from 1987 to 2019), at least 768 images (=32 × 12 × 2) were used to assess stand-replacing fire events for each site.  NGR is a near-infrared, red, green false-color image. RGB is the true-color image and SNR is the shortwave infrared, near-infrared, red false-color image. The blue box is a 1 × 1 km area around the site.

Field Data
To quantify post-fire vegetation recovery, a 15-meter radius sampling plot was established at each study site in 2018. All stems of regenerating trees within the sampling area were counted. Diameter at the base, diameter at breast height (at 1.3 m; only if seedlings or saplings exceeded this 1.3 m threshold), as well as height of trees were also recorded using a manual clinometer ('Clinomaster', Sisteco, Finland).

Landsat-Derived NDVI Time Series Data
Landsat satellite images at each site were used to assess fire disturbance, post-fire vegetation recovery and climate sensitivity analysis. These images were obtained via Google Earth Engine. The Landsat-5 Thematic Mapper (TM) sensor (from 1987 to 2011) and Landsat-7 Enhanced Thematic Mapper Plus (ETM+) sensor (from 1999 to 2017), deriving the Normalized Difference Vegetation Index (NDVI) (every 16 days at 30 m spatial resolution), were used. Gaps in the time series data were filled using a spline interpolation approach via the "na.spline" function, from the "greenbrown" package [51] in R [52]. On average, the percentage of gaps in the Landsat-5-derived NDVI time series from 1987 to 2011 was 16%. Similarly, the percentage of data gaps in the Landsat-7-derived NDVI time series from 1999 to 2017 was 24%. We conducted a general linear regression analysis to detect the relationship between annual maximum of NDVI and time since fire at each site.

Climate Sensitivity Analysis
To determine the response of NDVI to climate variability, we used a Bayesian Additive Regression Trees (BART) model-a machine learning approach-in the 'bartMachine' package [53] in R. BART is a non-parametric Bayesian regression approach which uses dimensionally adaptive random basis elements [54]. It relies on an underlying Bayesian probability model rather than a pure algorithm and provides uncertainty estimates of model outputs [53]. Readers are directed to paper [53] which provides R-based codes for both regression and classification analysis with associated detailed explanations of their respective applications.
Monthly air temperature and precipitation data time series were extracted from a gridded climate product (TerraClimate) which has a high-spatial resolution of ca. 4 km (1/24 • ) [33]. For each study site, climate data from 1987 to 2017 were used. Monthly maximum of Landsat-5-and Landsat-7-derived NDVI from each site and associated monthly air temperature and precipitation were used to build the algorithm in BART. To investigate whether vegetation at different phases of recovery would respond to climate differently, we separated NDVI datasets by time since fire (i.e., first 1-5, 6-10 and 11-15 years post-fire). As Landsat-5 time series datasets from 1987 to 2011 were used, we designated data collected from 1987 until fire year as the 'pre-fire group'. As mature stands did not experience any fire disturbance since 1987, we artificially divided the NDVI data recorded from 1987 to 1999 and from 2000 to 2017 as pre-and post-fire datasets to match up fire history of the burnt sites.
We built BART algorithms for each site and different time since fire (i.e., pre-fire, first 1-5, 6-10 and 11-15 years post-fire) respectively, using 1000 iterations with 50 trees, excluding a 250-iteration burnt-in period. NDVI time series data collected at each site over each period were treated as an individual training dataset. In BART, to assess model performance, in-sample (training data) root mean squared error (RMSE), Pseudo-R 2 and out-of-sample (testing data) Pseudo-R 2 and RMSE were computed. To conduct out-of-sample model assessment, a 10-fold cross-validation with 500 iterations was also used. A 95% credible interval was used to predict the uncertainty estimate from the BART posterior probability.
To quantify the responses of NDVI per site to climate variability, we introduced a set of six climatic scenarios, one baseline scenario, that reflects the typical climate in our study region, and five scenarios with modified temperature or precipitation distribution. All BART models were run with the same climate scenarios, that is climate did not vary from site to site but rather the site-to-site differences are encapsulated in the BART models. This approach was chosen so that responses depend only upon site characteristics and not varying climate between sites, which thus allows us to compare responses directly. None of the climate scenarios were unusual for the region. We can therefore use them for each site without the risk of forcing a model with a climate pattern that is very different from the climate that was used for model calibrations, which could clearly lead to inappropriate artefacts. As baseline climate, we adopted the long-term (1999-2017) monthly means of a meteorological station located centrally within our sites (ID: GHCND:RSM00030844-see earlier), with an annual mean air temperature of −2.1 • C and total annual precipitation of 364 mm. The average annual precipitation input is very comparable to other weather stations in the study region.
The climate scenarios included raising temperature by 2 • C and changes in precipitation distribution and total annual input, respectively. To simulate shifts in precipitation distribution in this strongly continental climatic region, we reduced or increased late spring to early summer rain input (April-June) by 60% and then adjusted to the late summerearly autumn period (July-September) accordingly, such that overall annual totals did not change (Figure 2). A ±60% shift in spring-early summer precipitation was chosen to simulate either a prolonged dry or moist spring-early summer condition, respectively. These adjustments were well within the range of historic minimum and maximum precipitation input (1999-2017) (Figure 2), over at least the past two decades and likely to continue in the future. In addition, we decreased or increased 30% of rain input during the entire growing season from April to September to simulate changes in total annual precipitation input conditions.
To summarize, our climate sensitivity analysis therefore included six climatic scenarios in total: (1) no changes in both air temperature and precipitation (with baseline data), (2) a 2 • C increase in annual mean air temperature with average precipitation, (3) 60% spring precipitation shift to summer while the annual totals and air temperature remained at long-term averages, (4) 60% increase in spring precipitation input while the annual totals and air temperature remained at long-term averages, (5) 30% reductions in both spring and summer precipitation input with average temperature and (6) 30% increases in both spring and summer precipitation input with average temperature. Note: we did not shift 60% summer rainfall input into spring months. This is because with an increased amount of rainwater, the total amount of precipitation in spring would exceed its historic records ( Figure 2).
Using the resulting BART models (for TM dataset-58 algorithms, for ETM+ dataset-47 algorithms), together with projected climate, we simulated NDVI at five-year intervals for 15 years post-fire. Vegetation response to climate variability, expressed as relative change (%), was determined by subtracting NDVI simulated under baseline climate from NDVI simulated under altered climatic conditions and then divided by NDVI simulated under baseline climate.
As our study sites present the majority of successional pathways in the study region, we then continued to investigate if the relative change of simulated NDVI across all the burnt sites would respond to proposed climatic scenarios differently at different times since fire respectively (i.e., first 1-5, 6-10 and 11-15 years post-fire, as previously used by Barrett et al. [17]), and if recruited species would induce changes in the relative change to climate variability. To answer the questions, a general linear mixed model was used to analyze the effects of altered climatic conditions (five scenarios: (1) increased air temperature, (2) dry spring and wet summer, (3) wet spring and dry summer, (4) dry spring and dry summer, as well as (5) wet spring and wet summer), time (1-5, 6-10 and 11-15 years post-fire), recruited species (i.e., grass, recruited Scots pine and mixed species) and their interactions as fixed effects upon the relative change of simulated NDVI. The analysis was performed with the 'lmer()' function, fitted into the 'lme4' package [55], and p-values were obtained by using the 'lmerTest' package [56] in R. Significance was set at p ≤ 0.05. A pairwise post hoc Tukey test was also conducted to test differences in combinations of previous rainfall patterns and time, by using the 'multcompView' package [57] in R.

Fire History and Tree Regeneration Post-Fire
The majority of our studied sites (13 out of the 16) were burnt in/around 2000 (i.e., 1996-2003) ( Table 1). In each case, the fire events were recorded in the spring and early summer seasons. In addition, soil texture analysis revealed that soils at our study site are all free-draining, being very sandy in general, with similar bulk densities (similar level of compaction).
Both TM-and ETM+-derived NDVI time series data showed a substantial decline post-fire across all sites, as illustrated by the examples in Figure 5. Furthermore, only sites that were dominated by Scots pine post-fire showed significant positive trends with time (Landsat-5 dataset: slope = 0.0068, F = 13.97, p < 0.001; Landsat-7 dataset: slope = 0.0074, F = 22.12, p < 0.001). However, sites that were dominated by grass, mixed species, as well as mature Scots pine, did not show a trend over approximately 20 years (p > 0.05).
wise post hoc Tukey test was also conducted to test differences in combinations of previous rainfall patterns and time, by using the 'multcompView' package [57] in R.

Fire History and Tree Regeneration Post-Fire
The majority of our studied sites (13 out of the 16) were burnt in/around 2000 (i.e., 1996-2003) ( Table 1). In each case, the fire events were recorded in the spring and early summer seasons. In addition, soil texture analysis revealed that soils at our study site are all free-draining, being very sandy in general, with similar bulk densities (similar level of compaction).
Both TM-and ETM+-derived NDVI time series data showed a substantial decline post-fire across all sites, as illustrated by the examples in Figure 5. Furthermore, only sites that were dominated by Scots pine post-fire showed significant positive trends with time (Landsat-5 dataset: slope = 0.0068, F = 13.97, p < 0.001; Landsat-7 dataset: slope = 0.0074, F = 22.12, p < 0.001). However, sites that were dominated by grass, mixed species, as well as mature Scots pine, did not show a trend over approximately 20 years (p > 0.05). Our field observations, collected in 2018, showed that post-fire tree recruitment varied markedly across sites ( Figure 6). Grass-dominated sites generally failed to recruit any seedlings post-fire or recruited very few seedlings per hectare. Abundantly recruiting Scots pine sites had an average stocking density of 8900 ± 1772 (n = 5; range 3750 to 13,250) stems ha −1 for the five sites, respectively. Mixed species stands generally had a greater number of regeneration seedlings compared to Scots pine-dominated sites, with an average stocking density of 18,389 ± 5137 (n = 4; range 33,250 to 10,500) stems ha −1 , respectively. For comparison, the mature Scots pine stands had an average stocking density of 1137 stems ha −1 . Average diameter of base (DB) and height of regenerated aspen, birch, larch and Scots pine varied substantially between sites (p < 0.001) (Supplementary Figure S2). Our field observations, collected in 2018, showed that post-fire tree recruitment varied markedly across sites ( Figure 6). Grass-dominated sites generally failed to recruit any seedlings post-fire or recruited very few seedlings per hectare. Abundantly recruiting Scots pine sites had an average stocking density of 8900 ± 1772 (n = 5; range 3750 to 13,250) stems ha −1 for the five sites, respectively. Mixed species stands generally had a greater number of regeneration seedlings compared to Scots pine-dominated sites, with an average stocking density of 18,389 ± 5137 (n = 4; range 33,250 to 10,500) stems ha −1 , respectively. For comparison, the mature Scots pine stands had an average stocking density of 1137 stems ha −1 . Average diameter of base (DB) and height of regenerated aspen, birch, larch and Scots pine varied substantially between sites (p < 0.001) (Supplementary Figure S2).
Remote Sens. 2021, 13, x FOR PEER REVIEW Figure 6. Post-fire stocking density of aspen, birch, larch and Scots pine from all the burnt were undertaking different successional pathways used in the study. Sites: grass (G), Scots p a mixture of species (Mix), respectively.

Response of NDVI to Temperature, Precipitation and Post-Fire Climate Sensitivity
The BART models performed well, particularly with the datasets that were c for post-fire events (Figure 7, Supplementary Tables S1 and S2). On average, in Pseudo-R 2 was 0.81 and 0.85 for the burnt sites at the 1-5 years and 6-10 years p TM time series dataset, whilst it was 0.79, 0.86 and 0.83 for the burnt sites at the 1-6-10 years and 11-15 years post-fire respectively, with the ETM+ time series datas of-sample Pseudo-R 2 was 0.64 and 0.71 for the burnt sites at the 1-5 years and 6-1 post-fire respectively, with the TM time series dataset, whilst it was 0.61, 0.72 and the burnt sites with the ETM+ dataset during the same time periods as mentioned Further, simulations with TM time series, recorded from 1987 to the year of fire, ge showed exceptionally low Pseudo-R 2 (ca. 0.40 for in-sample and 0.22 for out-ofcompared to the post-fire values. On the other-hand, in-sample Pseudo-R 2 of the Scots pine stand was low and remained stable throughout the study period (0.32 and 0.51 for ETM+). In-sample and out-of-sample RSME varied between 0.02 a across all the sites and time periods since fire (Supplementary Tables S1 and S2).

Response of NDVI to Temperature, Precipitation and Post-Fire Climate Sensitivity
The BART models performed well, particularly with the datasets that were collected for post-fire events (Figure 7, Supplementary Tables S1 and S2). On average, in-sample Pseudo-R 2 was 0.81 and 0.85 for the burnt sites at the 1-5 years and 6-10 years post-fire TM time series dataset, whilst it was 0.79, 0.86 and 0.83 for the burnt sites at the 1-5 years, 6-10 years and 11-15 years post-fire respectively, with the ETM+ time series dataset. Outof-sample Pseudo-R 2 was 0.64 and 0.71 for the burnt sites at the 1-5 years and 6-10 years post-fire respectively, with the TM time series dataset, whilst it was 0.61, 0.72 and 0.69 for the burnt sites with the ETM+ dataset during the same time periods as mentioned above. Further, simulations with TM time series, recorded from 1987 to the year of fire, generally showed exceptionally low Pseudo-R 2 (ca. 0.40 for in-sample and 0.22 for out-of-sample) compared to the post-fire values. On the other-hand, in-sample Pseudo-R 2 of the mature Scots pine stand was low and remained stable throughout the study period (0.32 for TM and 0.51 for ETM+). In-sample and out-of-sample RSME varied between 0.02 and 0.15 across all the sites and time periods since fire (Supplementary Tables S1 and S2).
OR PEER REVIEW 12 of 19 Figure 7. In-sample (a, b) and out-of-sample (c, d) coefficient of determination (Pseudo-R 2 ) for the relationship between measured and simulated NDVI by using Bayesian Additive regression Trees models (mean ± standard error). Models were performed for each site at different times since fire respectively (four periods: prior to fire event =, from 1987 to the fire year (with TM-derived time series), first 1-5 years, 6-10 years and more than 10 years post-fire (with ETM+-derived time series)). Recovery Pseudo-R 2 was calculated based on sites that were dominated by grass, Scots pines and a mixture of species (yellow bars), and Mature Pseudo-R 2 were from mature Scots pines only (grey bars).
Predicting NDVI under the different climatic scenarios did not depend on whether the models were calibrated with Landsat-5 TM or Landsat-7 ETM+ data-predictions from the Landsat-5 and 7 calibrated models were similar (p = 0.33). The modeling outputs from the models calibrated with the different Landsat data were then combined for further analysis. Relative changes of simulated NDVI under different climatic scenarios varied significantly across all the burnt sites (p < 0.001) (Figure 8). However, sites that were undergoing different successional pathways responded similarly to each climatic condition (Supplementary Figure S3). NDVI was enhanced under the increased air temperature scenario over all the three time periods since fire. Shifts in precipitation distribution, but with total annual input remaining at its long-term average, resulted in contrasting NDVI responses: a drier spring (April-June) followed by a wetter summer (July-September) reduces NDVI, indicating that prolonged drought in spring would significantly reduce greenness of all sites; in contrast, a drier summer to early autumn period following a wetter spring to early summer period could enhance greenness across all sites. A 30% reduction in rainfall input during both spring and summer months would reduce NDVI, whilst a 30% increase in precipitation would slightly enhance NDVI.
NDVI was significantly enhanced over the initial 5 years from burn compared to 6-10 and 11-15 years post-fire under the increased air temperature and wetter spring, drier In-sample (a,b) and out-of-sample (c,d) coefficient of determination (Pseudo-R 2 ) for the relationship between measured and simulated NDVI by using Bayesian Additive regression Trees models (mean ± standard error). Models were performed for each site at different times since fire respectively (four periods: prior to fire event =, from 1987 to the fire year (with TM-derived time series), first 1-5 years, 6-10 years and more than 10 years post-fire (with ETM+-derived time series)). Recovery Pseudo-R 2 was calculated based on sites that were dominated by grass, Scots pines and a mixture of species (yellow bars), and Mature Pseudo-R 2 were from mature Scots pines only (grey bars).
Predicting NDVI under the different climatic scenarios did not depend on whether the models were calibrated with Landsat-5 TM or Landsat-7 ETM+ data-predictions from the Landsat-5 and 7 calibrated models were similar (p = 0.33). The modeling outputs from the models calibrated with the different Landsat data were then combined for further analysis. Relative changes of simulated NDVI under different climatic scenarios varied significantly across all the burnt sites (p < 0.001) ( Figure 8). However, sites that were undergoing different successional pathways responded similarly to each climatic condition (Supplementary Figure S3). NDVI was enhanced under the increased air temperature scenario over all the three time periods since fire. Shifts in precipitation distribution, but with total annual input remaining at its long-term average, resulted in contrasting NDVI responses: a drier spring (April-June) followed by a wetter summer (July-September) reduces NDVI, indicating that prolonged drought in spring would significantly reduce greenness of all sites; in contrast, a drier summer to early autumn period following a wetter spring to early summer period could enhance greenness across all sites. A 30% reduction in rainfall input during both spring and summer months would reduce NDVI, whilst a 30% increase in precipitation would slightly enhance NDVI. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 19 Figure 8. Relative changes of simulated NDVI (mean ± standard error) under (a) increased (by 2 °C) monthly air temperature with average precipitation, (b) a 60% reduction in spring precipitation, whilst the annual total remained the same (dry spring, wet summer), (c) a 60% increase in spring precipitation, whilst the annual total remained the same (wet spring, dry summer), (d) 30% reductions in both spring and summer precipitation (dry spring, dry summer) and (e) 30% increases in both spring and summer precipitation (wet spring, wet summer) in grass-dominated, regenerating Scots pine, mixed and mature Scots pine stands.

Current Climate Sensitivity
The present study has shown that post-fire recruited vegetation was highly sensitive to climate variables in a forest-steppe zone of southeastern Siberia when compared to mature forest stands, as indicated by the greater coefficient of determination (R 2 ) values from the BART models. The high out-of-sample R 2 (>0.61) suggests that the larger proportion of the variance in NDVI can be explained by changes in air temperature and precipitation. This supports our first hypothesis and is unsurprising given that seedlings in their initial years after germination (i.e., 0-5 and 6-10 years) are known to be extremely sensitive to small variations in climate, while mature forests are much more resistant to unfavorable climatic conditions, because of more extensive and developed root systems that can use water conservatively [45,58]. The results of this study are consistent with a recent study that suggests that climate sensitivity reduces with increasing forest age in boreal-temperate forests in North America [59], and with those of Anenkhonov et al. [46] who showed that young Scots pines are more vulnerable under climate change compared to mature individuals in the foothills of the Khamar-Daban and Dzhidinskiy mountain ridges, south of Lake Baikal in SE Russia.
Scots pine and mixed species sites showed very little difference in climate sensitivity at every time interval (p > 0.05) (Figure 8). Therefore, our second hypothesis that post-fire stands with different dominant species would respond to climate variability differently has not been confirmed. This is surprising given that previous studies showed that Scots pine stands are of exceptionally high tolerance to water stress compared to other tree species in the study region [46][47][48]. There are a number of possible reasons for this. The first is that the climatic scenarios we proposed all fall within the ranges of approximately 20year historic minimum and maximum temperature and precipitation (Figure 2), suggesting that these climatic conditions can reoccur anytime now and in the near future. Thus, these species may be adapted to these conditions already. Whilst there are differences in the climate envelopes of the species, the study sites all fall within the overlap of these Relative changes of simulated NDVI (mean ± standard error) under (a) increased (by 2 • C) monthly air temperature with average precipitation, (b) a 60% reduction in spring precipitation, whilst the annual total remained the same (dry spring, wet summer), (c) a 60% increase in spring precipitation, whilst the annual total remained the same (wet spring, dry summer), (d) 30% reductions in both spring and summer precipitation (dry spring, dry summer) and (e) 30% increases in both spring and summer precipitation (wet spring, wet summer) in grass-dominated, regenerating Scots pine, mixed and mature Scots pine stands. NDVI was significantly enhanced over the initial 5 years from burn compared to 6-10 and 11-15 years post-fire under the increased air temperature and wetter spring, drier summer scenarios (p < 0.01; Figure 7). The combined drier spring and wetter summer scenario induced a comparable reduction in the relative change of simulated NDVI over the three time periods since fire (p > 0.05). Under the increased total precipitation input over the growing season scenario, NDVI was enhanced more from 6-10 and from 11-15 years post-fire than over the first 5 years (p < 0.001). In contrast, NDVI was reduced by a greater extent from 6-10 and from 11-15 years post-fire than over the initial 5 years from burn (p < 0.01) when the growing season gets drier.

Current Climate Sensitivity
The present study has shown that post-fire recruited vegetation was highly sensitive to climate variables in a forest-steppe zone of southeastern Siberia when compared to mature forest stands, as indicated by the greater coefficient of determination (R 2 ) values from the BART models. The high out-of-sample R 2 (>0.61) suggests that the larger proportion of the variance in NDVI can be explained by changes in air temperature and precipitation. This supports our first hypothesis and is unsurprising given that seedlings in their initial years after germination (i.e., 0-5 and 6-10 years) are known to be extremely sensitive to small variations in climate, while mature forests are much more resistant to unfavorable climatic conditions, because of more extensive and developed root systems that can use water conservatively [45,58]. The results of this study are consistent with a recent study that suggests that climate sensitivity reduces with increasing forest age in boreal-temperate forests in North America [59], and with those of Anenkhonov et al. [46] who showed that young Scots pines are more vulnerable under climate change compared to mature individuals in the foothills of the Khamar-Daban and Dzhidinskiy mountain ridges, south of Lake Baikal in SE Russia.
Scots pine and mixed species sites showed very little difference in climate sensitivity at every time interval (p > 0.05) (Figure 8). Therefore, our second hypothesis that post-fire stands with different dominant species would respond to climate variability differently has not been confirmed. This is surprising given that previous studies showed that Scots pine stands are of exceptionally high tolerance to water stress compared to other tree species in the study region [46][47][48]. There are a number of possible reasons for this. The first is that the climatic scenarios we proposed all fall within the ranges of approximately 20-year historic minimum and maximum temperature and precipitation (Figure 2), suggesting that these climatic conditions can reoccur anytime now and in the near future. Thus, these species may be adapted to these conditions already. Whilst there are differences in the climate envelopes of the species, the study sites all fall within the overlap of these envelopes. The second possible explanation for the similarity between recruitment trajectories (Figure 8) is that immediately after a stand-replacing fire, tree saplings are often surrounded by grass and/or shrub species (as seen from the other newly burnt sites during field campaigns). This could explain the similarity of both the pine and mixed trajectories to the grass trajectory, particularly for the first 5 years, but not the similarity of the pine and mixed to each other after more than 50 years (as in our mature stands). The third possible explanation for the similarity of trajectories is that the spatial and temporal resolution of the Landsat data is insufficient to detect the differences. Further research using a higher resolution remote sensing dataset and/or more extensive field observations is needed to determine which, if any, of these reasons is the cause behind this surprising result.

Recovery under Changes in Climatology
Under current CO 2 emissions (IPCC RCP 4.5, 6.0 and 8.5), the world is on track to reach a 2 • C increase in global mean air temperature by 2046-2065, and boreal forest is no exception [31]. Figure 8 shows that a 2 • C increase in air temperature could enhance greenness and thus speed up recovery of all sites post-fire in forest-steppe zones of southern Siberia, assuming no change in precipitation regime with an average annual precipitation input of 364 mm. This finding is consistent with the results of other field [60], satellite [61,62] and modeling [63] based studies in boreal forest and steppe ecosystems, in which it was found that a warming climate can promote regeneration success and biomass increases in the boreal forest biome. In northern mid-to high-latitudes (>50 • N), rising temperature is the dominant driver of the positive greening trends [64] because the onset of photosynthesis in the spring in boreal forests is triggered by air temperature [65,66]. Thus, a warmer spring and autumn can extend the entire period of photosynthetic activity in evergreen species. For deciduous species and grasses specifically, warmer air temperature can also result in an earlier leaf expansion, as well as potentially delaying autumn senescence, which would also lead to a longer length of growing season [67]. Interestingly, the modeled NDVI at the burnt sites responded more positively to rising temperature over the initial 5 years from burn compared to 6-10 and 11-15 years post-fire, suggesting that plants are likely to be more susceptible in the earlier stages following fire disturbance.
Unsurprisingly, the simulations in this study show that a higher amount of precipitation input during the growing season (April-September) resulted in an increase in NDVI, while a 30% reduction in precipitation input over the growing season reduced NDVI. Our study region generally lacks any substantial snow input during the winter season (Figure 3), and thus the total amount of precipitation input during the growing season ultimately determines the overall greenness and vegetation productivity. Even winter precipitation can be an important factor in the survival of seedlings post-fire [17]. Here, Barrett et al. [17] aimed to distinguish among areas in southern Siberia that recruit trees abundantly post-fire, intermediate recruitment and recruitment failure (AR, IR and RF, respectively) using remotely sensed data and information from field observations. They used a long time series of remotely sensed indices of greenness, moisture, snow cover and fire severity from Landsat and phenological characteristics derived from MODIS indices to compare these recruitment types prior to and post-burning events. They found that the earliest separability of AR and RF sites using remotely sensed indices occurs in the winter months, 3-4 years post-fire.
In a dry environment, with annual precipitation input less than 400 mm, it is expected that a higher precipitation input will cause greater vegetation coverage and productivity [68]. This result does raise concerns about the future sustainability of forests in this southern boreal-steppe region, with a recent study showing that Scots pine stands are experiencing post-fire recruitment failure related to moisture limitations [17].

Recovery under Changes in Climate Seasonality
While changes in total precipitation are important, our results suggest that shifts in precipitation distribution (annual precipitation remains the same, but with shifts in seasonality) can also lead to contrasting recovery outcomes (Figure 8). For example, increased spring rainfall input with a lower summer rainfall can speed up the recovery of vegetation post-fire, with an increase in modeled NDVI of 1.2%. This is comparable to the 1.3% increase that was estimated under a 30% increase in total rainfall. This importance of seasonality is reinforced by our finding that a shift in precipitation from spring to summer resulted in a reduction in the model NDVI of ca. 5.2% in the initial 5 years post-fire. This reduction in NDVI was 3.8% for the 6-10 years and 4.5% for the 11-15 years post-fire periods, highlighting the strong impact of spring precipitation at all vegetation developmental stages post-fire, with drier springs resulting in slow recovery and development of vegetation from both forest-and grass-based ecosystems.
It has been suggested that in Siberian forests, trees are more likely to accumulate carbon mainly in the production of leaves and needles at the beginning of the growing season (April-May), and only subsequently (June-July) do they begin to store carbohydrates in the stem wood, whilst continuing the allocation to green parts [69,70]. Furthermore, Mencuccini and Grace [71] have shown a greater ratio of leaf to sapwood area in relatively 'moist' compared to 'dry' Scots pine stands, suggesting that trees can invest more effort in producing photosynthetic tissues when environmental conditions become wetter. In this manner, a dry spring would slow the development of foliage, which can further reduce leaf area, shorten length of growing season and ultimately result in a reduction in vegetation coverage and productivity by delaying leaf unfolding and by reducing net photosynthesis. Several studies carried out in Eastern and central Siberia have confirmed that the maximum daily net photosynthesis was strongly influenced by spring soil water availability [72] and that radial growth is positively affected by precipitation just before, and at the beginning of, the growing season [73,74]. In the present study, enhanced summer precipitation, from July to September only, did not necessarily compensate for the overall reduction in NDVI. This is likely since high evapotranspiration demand during summer can potentially offset extensive summer rainfall input [75,76]. Additionally, plants may not be capable of using immediately available water for gaining canopy coverage and biomass, as Berner et al. [32] revealed that tree growth positively correlated with prior summer moisture availability. Instead, extensive spring precipitation is highly likely to be used by plants in summer. This is because an increase in precipitation intensity can push rain water deeper into the soil, particularly into a sandy soil such as those of the Zabaikal region and the Republic of Buryatia [48]. Subsequently, the stored soil water formed at depth in spring can be used progressively in summer months. In our study region, 60% of total spring rainfall input is approximately equivalent to the magnitude of 30% of summer rainfall input (Figure 2). A dry or wet spring with a wet summer respectively, as shown in Figure 7, induced contrasting NDVI responses. The finding of this study, that seasonality can have a significant impact, which is even greater than changes in annual precipitation input, alerts us to the fact that future climate and vegetation model-based studies need to fully consider the role of precipitation seasonality in determining vegetation coverage and productivity.
Machine learning algorithms have been increasingly used in remote sensing timeseries data over recent years [77,78]. The BART model is particularly advantageous when comparing with other machine learning regression algorithms (e.g., random forests and boosted regression trees), mainly because BART relies on an underlying Bayesian probability model rather than a pure algorithm and thus provides uncertainty estimates of model outputs [53]. In this study, the BART model generally performed well in both detecting current climate sensitivity and projecting future changes under varying climatic scenarios, which can be used in similar regions and cases.

Conclusions
Intensified fire activity can induce different successional pathways in a forest-steppe zone of southern Siberia. The present study suggests that inter-annual climate variability affects post-fire vegetation recovery in the region, irrespective of any difference in recruited species. In addition, a prolonged dry spring has a strong impact consistently over the entire early developmental stages, from the initial 5 years to 15 years post-fire. Higher spring precipitation input is likely to induce an increase in photosynthetic tissues and possibly earlier leaf unfolding, which in turn increases greenness and productivity. Our study highlights the importance of precipitation seasonality in post-fire boreal forest recovery, which warrants greater focus in future climate and vegetation model-based studies. Further, Landsat-derived NDVI has shown a great sensitivity to changes in temperature and precipitation, particularly post-fire disturbance, which can be considered as a good indicator for further investigating the interactive effect of climate and fire disturbance on forest recovery at a broader scale. In addition, latest generation sensors beyond Landsat (such as Sentinel 2 of the Copernicus series), although restricted in image run years, offer both better spatial and spectral resolution going forward. This offers the future possibilities of further significant improvements to the methodologies reported here.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13122247/s1, Figure S1: (a) Locations of available meteorological station data in south eastern Siberia within the study region, (b) boxplot showing median, 25%, 75% quartiles, maximum and minimum (whiskers) and outlier (dot) of annual precipitation (mm) for the six stations indicated in (a). Four stations ("AGINSKOE", "CHITA", "HILOK" and "KRASNYJ") from west to east recorded comparable annual precipitation input over the course of the past 40 years (from 1980 to 2018) (Oneway ANOVA test: F3,152 = 0.81, p = 0.489). This indicates that the majority of our study region has similar annual precipitation input. Precipitation is lower (ca. 250 mm) to the extreme east of the set, but still highly variable, as shown by maximum and minimum values. Figure S2: Diameter at base (a) and height (b) of regenerated trees at our sites that were dominated by grass (G), Scots pine (SP) or a mixture of species (Mix). Figure S3: Relative changes of simulated NDVI under (a) increased monthly air temperature (by 2 • C) with average precipitation, (b) a 60% reduction in spring precipitation, whilst the annual total remained the same, (c) a 60% increase in spring precipitation, whilst the annual total remained the same, (d) 30% reductions in both spring and summer precipitation and (e) 30% increases in both spring and summer precipitation at each individual site. Spring season in this study is defined as being from April to June, and summer season is from July to September for sites that were dominated by grass (G), Scots pine (SP), a mixture of species (Mix) and mature Scots pine (C), respectively. Table S1: In-sample coefficient of determination (Pseudo-R 2 ) and root mean square error (RMSE) for the relationship between measured and simulated NDVI by using Bayesian Additive Regression Trees models. Models were performed for each site at different times since fire respectively (four periods: prior to fire event-from 1987 to the fire year (with TM-derived time series), first 1-5 years, 6-10 years and more than 10 years post-fire (with ETM+-derived time series)) for sites that were dominated by grass (G), Scots pine (SP), a mixture of species (Mix) and mature scots pine (C), respectively. Table S2: Out-of-sample coefficient of determination (Pseudo-R 2 ) and root mean square error (RMSE) for the relationship between measured and simulated NDVI by using Bayesian Additive Regression Trees models. Models were performed for each site at different times since fire respectively (four periods: prior to fire event-from 1987 to the fire year (with TM-derived time series), first 1-5 years, 6-10 years and more than 10 years post-fire (with ETM+-derived time series)) for sites that were dominated by grass (G), Scots pine (SP), a mixture of species (Mix) and mature scots pine (C), respectively.