Modelling the impact of climate change on the occurrence of frost damage in Sitka spruce ( Picea sitchensis ) in Great Britain

Climate change is predicted to increase temperature and seasonal temperature variance in Great Britain (GB). Sitka spruce ( Picea sitchensis (Bong.) Carr) is the most important tree species used in commercial plantations throughout Europe and GB. Frosts that occur outside the winter dormancy period can negatively affect trees, since they happen after dehardening. Damage can be especially severe at bud burst, before emerging needles mature and form protective barriers. Here, we modelled the impact of climate change on frost sensitivity in Sitka spruce with temperature data from ﬁve climate projections. The UKCP09 climate model HadRm3 uses emission scenario SRESA1B for the years 2020–2099. The global and downscaled versions of the UKCP18 HadGem3 model use the emissions scenario RCP 8.5. The global model CMCC-CM uses the RCP 4.5 and RCP 8.5 emissions scenarios. The predictions based on these models were compared with results from gridded historical data for the period 1960–2015. Three indicators that assessed the frost sensitivity of Sitka spruce were explored: the total number of frosts between the onset of dehardening and the end of summer, which use three different temperature thresholds (Index 1 0 ◦ C , 1 –3 ◦ C , 1 –5 ◦ C ); the total number of frosts after bud burst (Index 2); the number of days with minimum temperatures below the resistance level (backlashes) during the hardening–dehardening period (September–August) (Index 3). The indices were validated with historical data for frost damage across GB, and Index 1 –3 ◦ C , Index 1 –5 ◦ C and Index 3 were shown to be signiﬁcantly correlated. The frequency of all frosts and backlashes is expected to decrease with climate change, especially under higher emissions scenarios. Post-bud burst frosts have been historically very rare in GB and remain so with climate change. Downscaled regional climate models detect geographic variability within GB and improve prediction of overall trends in frost damage in comparison to global climate change models for GB.


Sitka spruce in Great Britain
Sitka spruce (Picea sitchensis) is the most important commercial conifer species in Great Britain (GB). In the UK, Sitka spruce plantations cover ∼665 000 hectares and represent 25 per cent of the total woodland area (Forestry Commission, 2019). Observations of Sitka spruce in its native North America have indicated that genotypes from southern origins, such as Washington, experience higher levels of spring frost damage than northern origins such as QCI (Haida Gwaii/Queen Charlotte Island) (Forestry Commission, 1957). As QCI origin was suitable to local conditions in the UK, it was deemed unnecessary to identify trees of greater frost resistance from more northern origins, such as Alaska (Forestry Commission, 1957) and QCI was selected as the provenance of Sitka spruce for importation to the UK. Subsequent timber quality improvements were made by the breeding programmes of British tree nurseries (Lee, 1999).

Geographic regions of GB
The climate of the island of GB is classified as temperate oceanic (Peel et al., 2007) with temperature and precipitation regimes that differ across geographic regions. To account for these differences, GB has been divided into four climatically similar regions of provenance that divide GB north to south and east to west ( Figure S1), as defined in the Forest Reproductive Material Regulations in 1977 for a range of native species (Herbert et al., 1999). It is important to note that the regions of provenance are defined for native species at current and historical conditions, which may Impact of climate change on the occurrence of frost damage in Sitka spruce change with climate change. Seed zones are subdivisions of the regions of provenance, divided according to geoclimatic variation and natural boundaries.

Freezing damage
In plants, subfreezing temperatures can cause the destruction of cell walls and membranes due to ice formation (Holopainen and Holopainen, 1988), and protein denaturation. As perennial species, trees have evolved mechanisms to survive periods of cold temperatures during the winter season, such as protection of meristematic tissue by the formation of buds (Kuprian et al., 2014;Sakai, 1978) or cold hardening of sensitive tissues.
Cold hardening is triggered by shortening photoperiods and augmented by temperatures below 5 • C, environmental factors that also lead to the termination of bud development (Cannell et al., 1990). A study in white spruce (Picea glauca) found that a reduction in photoperiod is sufficient to start the hardening process, whereas cold temperatures alone do not induce hardening (Hamilton et al., 2016).
Cold hardiness is metabolically costly and there is little growth during the frost hardy period. This is why as soon as temperatures begin to rise in spring, trees lose these adaptations and deharden quickly, so that the threshold temperature at which frost damage could occur increases by up to 2 K day −1 (Neilson et al., 1972;Repo, 1992). During the process of dehardening the frost hardiness of bud tissue decreases to become the most sensitive tissue in trees, more so than needles to frost (Aitken and Adams, 1997;Anekonda et al., 2000). Buds also lose their scales during the period before bud burst, which confer protection during bud burst (Alfaro et al., 2000;Cannell and Sheppard, 1982;Cannell et al., 1990;Man et al., 2016).
Frost hardiness can be quantified using the median lethal time, LT 50 , that is defined as the temperature at which 50 per cent of trees subjected to a temperature die (Zhang and Willison, 1987). Spring frost damage is commonly reported after bud burst and during the growing stages, whereas damage to trees by frost in winter or before bud burst, is rare (Man et al., 2016), because the hardening process confers resistance to temperatures below the lowest temperature that typically occurs during this period (Sakai and Larcher, 1987).

Bud burst in Sitka spruce
Temperature is well-known to be the main abiotic driver of bud burst (Aitken & Hannerz, 2001;Cannell & Smith, 1983). Bud burst requires a genetically determined period of chilling, followed by a period of warmth that is accumulated over time, known as Growing Degree Days (GDD) (Cannell et al., 1990). Warmer temperatures seem to advance bud burst dates, although that advance is limited by an inverse relationship between the thermal time needed for bud burst and the chill days requirement (Bigler and Bugmann, 2018;Cannell and Sheppard, 1982;Fu et al., 2015;Morin et al., 2009Morin et al., , 2010Murray et al., 1989). Higher temperatures seem to decrease chilling more than they increase thermal time, which reduces the effect of warming on bud burst hastening, resulting in a nonlinear effect (Bigler and Bugmann, 2018;Fu et al., 2015;Morin et al., 2009Morin et al., , 2010. Early models of changes in phenology were based on daily mean temperatures, but recent studies have found that models that account for maximum temperatures, T max , and minimum temperatures, T min , can improve accuracy (Fu et al., 2016) through the contrasting effects of chilling and warming. Recently, minimum temperatures were shown to contribute towards the chilling requirement for bud burst, whilst maximum temperatures provide the warming to advance bud burst (Meng et al., 2019;Meng et al., 2020). Understanding the contrasting effect of accumulated heating and chilling on bud burst has refined model accuracy.
The recent work of Peaucelle et al. (2019a) showed that GDD and chilling explain 30 per cent of spatial variance of leaf unfolding, while the inclusion of water and light availability (intensity and photoperiod) combined explained another 10 per cent of spatial variance. Indeed, a study that examined the phenology of birch (Betula pendula Roth.), beech (Fagus sylvatica L.) and oak (Quercus robur L.) trees, found that maximum temperatures had a greater effect on birch and beech than on oak trees (Fu et al., 2016). Thus, species-specific phenological models need to be developed to account for this added complexity.
In general, species with a low chilling requirement will have a longer bud burst to bud set period mediated by global warming (Ma et al., 2019;Morin et al., 2009;Murray et al., 1989). However, within species plasticity results in some populations possessing both chilling sensitive and insensitive characteristics (Morin et al., 2009). Plants that are adapted to cold climates generally have a longer chilling requirement, otherwise any short period of warmth in winter could prematurely deharden tissues and result in tissue death. Global warming is likely to result in a slight increase in the length of the bud burst-bud set period, but could also increase the cold hardiness threshold temperature in spring, thus increasing the likelihood of frosts with temperatures below the threshold (Richardson et al., 2018). However, this reduced frost tolerance may be compensated by a reduction in the total number of frosts (Pohl et al., 2019), and by an increased period between bud burst and damaging frosts (Bigler and Bugmann, 2018).
According to Cannell and Smith (1983), 92 per cent of the variation in bud burst dates in Sitka spruce trees is explained by the number of GDD needed for bud burst when exposed to a certain number of chill days (equation 1). As the number of days to bud burst decreases, so does cold hardiness, reaching the highest threshold temperature around the time of bud burst (Aitken and Hannerz, 2001). In the work of Cannell and Smith (1983) Scottish forests comprised of QCI provenance Sitka spruce were used to calculate the bud burst date, and the GDD model produced is still considered very accurate for current climate conditions (Cannell et al., 1985), and more complex models are frequently found to be less accurate (Olsson and Jönsson, 2015).

Climate change effects on frost tolerance
In an early study, Cannell and Smith (1983) reported that spring frosts that are potentially damaging to Sitka spruce have a return time of 3-5 years in upland regions. As trees are most sensitive to frosts during the first 3-4 years after planting, an increase in the return time of frosts could allow some young plantations to completely avoid frost damage. If frosts that cause severe damage to young plantations only happen every 7-8 years, some Forestry tree plantations could go for their entire establishment period without damages caused by frost.
Climate change will affect the conditions under which frosts occur through the increased likelihood of warmer winter, autumn and spring periods, which may impact bud set, bud burst, winter hardiness and spring hardiness (Guak et al., 1998;Murray et al., 1989). Although warmer temperatures hasten metabolic processes, the annual total duration of primary production may be unchanged, due to earlier seed maturation and senescence (Reyes-Fox et al., 2014). Thus, trees tend to start bud formation earlier in autumn (Repo, 1992;Richardson et al., 2018), and bud burst starts later in spring, due to the reduction in chilling during the winter (Bigler and Bugmann, 2018;Morin et al., 2009Morin et al., , 2010Murray et al., 1994).

Estimates of frost occurrence and impacts on vegetation
Estimates of the frequency and harmfulness of spring frosts can be achieved in two ways: (1) by observing the frequency of total frosts; or (2) by linking frost occurrence to the expected phenological or physiological state of plant tissues and estimating the harmfulness of the frost based on that state. Modelling frost damage of coniferous tree species requires models that account for the damage to mature needles that have not been shed, as well as the tenderness of newly emerging tissues. Models suited for annual crops that simply time emergence to specific climatic events would not be suitable because of differences in the underlying biology. Subsequently, linking the expected frost damage to the phenological state, e.g. bud burst or leaf unfolding, seems to be the most common method of estimating frost damage for trees (Jeong et al., 2018;Jönsson et al., 2004;Ma et al., 2019;Meier et al., 2018;Morin and Chuine, 2014;Vitasse and Rebetez, 2018;Zohner et al., 2020).
Global circulation models (GCM) can estimate the occurrence of frost events under predicted climates. Thus, modelled climate data have often been used as the first approximation of predicted frost damage (Jönsson et al., 2004;Woldendorp et al., 2008). Recently, Jeong et al. (2018) have increased GCM complexity by the inclusion of additional microclimatic factors (e.g. humidity or snow cover) in an attempt to refine predictions. However, Ambroise et al. (2020) exemplified the complexity of including such factors in models without considering the associated plant interactions by showing that decreased snow cover could expose fine root tissue to damaging low temperatures that may result in stem and needle dehydration due to loss of hydraulic function.
Models that are also parameterised with the physiological state of the plant can take into account the interactions between frosts and physiologically mediated frost hardiness to improve accuracy (Jönsson et al., 2004;Kellomäki et al., 1995;Rochette et al., 2004;Woldendorp et al., 2008).

Aims and objectives
The overarching aim of this study was to estimate the occurrence of frosts and potential frost damage to Sitka spruce trees under future climate change scenarios. Our specific objectives were to: (1) approximate future trends of frost occurrence using three global and two regional climate change models; (2) utilize three indicators to determine how predicted future climate alters the occurrence and interaction of frosts with the physiological and phenological state of trees throughout the spring; (3) examine trends throughout GB to see whether there are regional differences in the expected changes to frost damage; and (4) explore how the resolution of different climate models influences model prediction precision. Figure 1A) measures the number of days with daily minimum temperatures below a threshold of 0 • C, −3 • C and − 5 • C following a five-day consecutive period of daily mean temperatures above 5 • C after the 1 January (late winter to early spring), which is assumed to be the trigger that starts the dehardening process, to the end of August (Jönsson et al., 2004). The purpose of the index is to provide an estimate of the total number of frosts that occur between late winter and late summer, as a first approximation to the effects of climate change. The different temperature thresholds were selected for different reasons. First, that 0 • C is the freezing point of water and is a natural threshold to use. Second, −3 • C was chosen as the threshold where trees that had undergone supercooling to achieve baseline frost hardiness would likely be damaged. And, third, −5 • C was chosen to determine how frequently temperatures below the natural frost hardiness level occurred. After validating Index 1, the threshold that was correlated with historical frost damage was used for Index 2.

Climate change indices
Bud burst dates were calculated using equation 1, which estimates the thermal time, measured in GDD with a base temperature of 5 • C, that plants require for bud burst. For bud burst in Sitka spruce the chill days were defined in this study as days with mean temperatures below 5 • C that occurred between the 1 November and bud burst. The model uses a floor value of 85 for the chill days, setting it at that number even when actual chill days are below that, and a maximum level of 180 (Cannell et al., 1985). The thermal time is defined as GDD(T base > 5 • C) between 1 February and bud burst and is dependent on the chilling experienced by the trees.
It was found that for Sitka spruce populations of QCI provenance (the most common provenance used in commercial forestry in GB), average values of the coefficients were a = 67.4, b = 4401.8, c = −0.042 (Cannell et al., 1985). To account for the variation that could be found among improved QCI origin trees coefficients at both extremes were taken, from a maximum of a = 100, b = 7500, to a minimum of a = 50, b = 2500.
Chill days are calculated first between the 1 November and the 31 May. An approximate bud burst date is then calculated from this value which is used to recalculate the value for the number of chill days that more accurately approaches reality, since chill days that occur after bud burst should not be counted towards the chill days total. After the GDD needed for bud burst are calculated using equation 1, GDD accumulated from the 1 Impact of climate change on the occurrence of frost damage in Sitka spruce January is calculated, and the date at which the calculated value is achieved is considered the bud burst date.
Index 2 ( Figure 1B) measures the frequency of daily minimum temperatures below −3 • C after bud burst (spring to summer period), calculated by equation 1, using the mean values of the coefficients (Cannell et al., 1985). Index 2 max does the same for the maximum value of the coefficients (late bud burst populations), whereas Index 2 min uses the minimum value of the coefficients (early bud burst date populations). These indices were used to evaluate whether frosts occurring after bud burst are a significant factor in expected historical and predicted future frost damage, for populations with differing sensitivities to chilling and warming.
Index 3 measures the number of days with a daily minimum temperature below the hardiness level during the hardening-dehardening period. The hardening-dehardening period is defined as the period between the 1 September and the 31 August of the next year. These dates were chosen because the autumn equinox, the 22 September, was considered the starting point of the hardening process in autumn. The daily hardiness threshold was a function of daily mean temperatures, the rate of hardening, and the rate of dehardening, limiting the highest hardiness temperature to −3 • C (Cannell and Sheppard, 1982), and the lowest hardiness temperature to −40 • C (Cannell et al., 1990). This index gives a more detailed approximation of the expected levels of damage with climate change, depending on the effects of both the warming temperatures and the occurrence of frosts.
The rate of hardening and dehardening used for Index 3 were described by a sigmoidal function of the daily mean temperatures (equation 2). Where d defines the lower asymptote, or the minimum rate of hardening-dehardening, and a is the distance between the upper and lower asymptote, which when combined with d describes the maximum rate of hardening-dehardening, and b and c define the shape of the curve, with c defining the centre of the sigmoidal curve when c = T mean . Parameter b defines whether the sigmoid slope is decreasing or increasing, depending on whether it is positive or negative.
Rate of (de)hardening = a Since a decreasing photoperiod affects hardening rates even at warm temperatures, for the days between the autumn equinox and winter solstice, or the winter dormancy period, a minimum value of 0.1 K day −1 was established for hardening. For the rest of the year, the minimum hardening rate was set to 0 K day −1 . The maximum rate of hardening, 1.15 K day −1 , was estimated to be reached at −10 • C (Cannell et al., 1990), whereas the minimum rate of hardening was reached at 5 • C. The coefficients for the curve were calculated by simulating a dataset describing a symmetrical sigmoidal curve with the parameters contained in Table 1. For the hardening rate during the winter dormancy period, values of a = 1.05, b = 0.96, c = 2.5 and d = 0.1, were estimated for equation 2, whereas for the rest of the year, the values were a = 1, b = 0.96, c = 2.5 and d = 0.
The maximum rate of dehardening in Sitka spruce was not found in the available literature, so a value of 1.25 K day −1 was used as a proxy, taken from the average of published dehardening rates (Repo, 1992) for Norway spruce (Picea abies; 6 K week −1 ), and Scots pine (Pinus Sylvestris; 11 K week −1 ). The minimum rate of dehardening, 0 K day −1 , was assumed to be reached at 5 • C, whereas the maximum rate was reached at 15 • C. For the dehardening rate, values of a = 5, b = −1.69, • C = 10.07 and d = 0, were estimated for coefficients of the sigmoid curve.
All the indices were calculated for all years using R 3.6.3 (R Foundation for Statistical Computing and R Core Team, 2019) with the RStudio IDE 1.3.959 (RStudio Team, 2019). The ncdf4 1.17 R package (Pierce, 2019) was used to open and manipulate the minimum and mean temperature ncdf format datafiles. The raster 3.1-5  and rgdal 1.4-8  R packages were used to manipulate the mean and minimum daily temperature data, extracted from the ncdf datafile and organized into a uniform grid of GB. The dplyr 1.0 R package (Wickham et al., 2020a) was used to calculate the indices for each grid coordinate. The PMCMRplus 1.4.4 R package (Pohlert, 2020) was used for the Bonferroni-Dunn post-hoc test after the Kruskal-Wallis H test.
The indices did not follow a normal distribution. The distribution of each index is the sum of different non-correlated non-normal distributions in different zones, thus the overlap of different distributions across different zones was not normal and could not be transformed. The non-parametric Kruskal-Wallis H test was used to compare the indices by different factors, with a post-hoc Bonferroni-Dunn test used for pairwise comparisons. A model that uses a higher emissions scenario, the UKCP18, was also selected. The UKCP18 Regional Projections uses a 12 × 12 km grid over the UK for 1980-2080 for the dynamically downscaled

Spatial differences
The gridded data for each year was grouped into provenance regions in order to establish differences between biologically relevant geographical areas across GB (Forestry Commission, 2012). The matrix was then transformed into a data frame with the addition of coordinates for each grid point. The sp 1.4-2 R package (Pebesma et al., 2019) was used to extract the polygon coordinates for each provenance region from the GB Forestry Commission (FC) provenance region shapefile, overlay the data frame and group all points by zone. A Kruskal-Wallis H test was conducted with indices as dependants and the zones as factors.
The level of spatial autocorrelation was measured by calculating Moran's I, with a Markov chain of 10 000. The shapefiles 0.7 (Stabler, 2013)

Validation
Data on the level of tissue damage to 41 commercial forests across GB by frosts in the spring of 2015 was obtained from Tilhill forestry. Individual trees were scored on a scale of 0-100 points scale by evaluating both the presence of leader damage (the most harmful frost damage for commercial forests, thus assigned 50 points if present) and the level of observed frost damage among lateral branches. The final score was a sum of the overall percentage of damage and the presence/absence of leader damage (Table S1).
In order to calculate the coefficient of correlation between the indices and the observed values of frost damage, a correlation matrix was calculated using the Hmisc 4.4-1 R package (Harrell, 2020), giving the value of the Pearson's r correlation coefficient and the P-value for the correlation. Impact of climate change on the occurrence of frost damage in Sitka spruce

Validation
The percentage value of frost damage was correlated with the number of frosts below −3 • C, Index 1 -3 • C (Pearson's r = 0.47, Pvalue < 0.01), and with the number of frosts below −5 • C, Index 1 -5 • C (Pearson's r = 0.34, P-value < 0.05), for the forests evaluated. A correlation was also found between the percentage of damage and the number of backlashes, Index 3 (Pearson's r = 0.23, P-value < 0.05).
No correlation was found between the percentage of damage and the number of frosts below 0 • C, (Index 1 0 • C ), or frosts occurring after bud burst (Index 2).

Overall changes with climate change
The number of frosts between the onset of dehardening and late summer for all three reference temperatures of Index 1 (0 • C, −3 • C and −5 • C), was significantly reduced with climate change in all climate change scenarios (P < 0.001), when compared with the average for 1960-2015, with the exception of the CMCC RCP 4.5 model for both the −3 • C and −5 • C threshold temperatures ( Figure 2). This follows a general trend during the historical period, with reductions in the occurrence of frosts between 1960 and 2015 (Table S2).
For the CMCC RCP 4.5 predicted scenario, the frequency of frosts <−3 • C and <−5 • C increased, and the change was significantly different from the historical average for the −5 • C threshold temperature (P < 0.001). All predicted climate change scenarios were significantly different between each other (P < 0.001). The scenarios with higher CO 2 levels showed a higher decrease in the incidence of frosts. However, there were differences between different models that used the same emissions scenarios; the CMCC model showed a lower decrease of frost incidence than the UKCP18 model ( Figure 2). There was also a difference in the scope of the modelled area; the global HadGEM3 model showed a higher decrease than the downscaled version.
The number of late season frosts that occur after bud burst (Index 2) decreased significantly in all scenarios, for both average, early, and late bud burst dates (Index 2, 2 min , 2 max ), when compared with the average for 1960-2015. There was, however, an exception for the UKCP09 model, based on the SRESA1B emissions scenario, which, for early bud burst dates, showed an increase in the number for frosts after bud burst ( Figure S2).
There were no significant differences in Index 2 min between all climate change models used. For Index 2, only UKCP09 and UKCP18 downscaled outputs remained significantly different (P < 0.001). For Index 2 max , the UKCP09 model differed significantly from the rest (P < 0.001), whereas the other four showed no significant differences.
The number of frosts after bud burst was very low, both historically and with changes due to climate change. In fact, the average number of post-bud burst frosts <−3 • C per year, Index 2, was <0.05, for all provenance regions and periods (Table S3). The average return time for a post-bud burst frost was found to extend to the entire length of the period evaluated (55 years for the historical period, 80 years for the 2020-2100 period).
Bud burst dates were significantly hastened (P < 0.001) in all climate models except the CMCC RCP 4.5 model, which showed The average number of frosts <0, −3 and −5 • C between the onset of dehardening and late summer (Index 1) for the historical  and the predictions for 2020-2100 made by five climate change datasets (CMCC model's data for RCP4.5, the UKCP09 datasets for model HadRM3, the CMCC model's data for RCP 8.5, the UKCP18 dataset for both the downscaled and global HadGEM3 model). Values given for all four provenance regions of GB (10,20,30,40), and the average for the entire GB. Note that Y axes in all three figures have different scales. a delay in bud burst by 3 days (Table S4). The UKCP18 model hastened bud burst dates by >50 days for both the global and downscaled model resolution. UKCP09 also hastened the bud burst date by 17 days, 2 weeks more than the CMCC RCP 8.5. All differences between models were significant (P < 0.001). This follows a similar trend between 1960 and 2015, with a hastening of the bud burst date between 1960-1979 to 1980-1999 and 1980-1999 to 2000-2015 (Table S5).
The number of chill days (T mean < 5 • C) during winter decreased significantly with all climate change scenarios (P < 0.001), except Forestry Figure 3 Average number of backlashes, days during which the frost tolerance of a tree is predicted to be higher than the minimum temperatures achieved that day (Index 3). Values given for the historical  and the predictions for 2020-2100 made by five climate change models (CMCC model's data for RCP4.5, the UKCP09 datasets for model HadRM3, the CMCC model's data for RCP 8.5, the UKCP18 dataset for both the downscaled and global HadGEM3 model). Values given for all four provenance regions of GB (10,20,30,40), and the average for the entire GB.
for the CMCC RCP 4.5 model, by a range of 1-15 days (Table S4). The CMCC model showed the smallest decrease, and the higher CO 2 emissions scenario showed a higher effect on chill days (Table S4). All differences between different climate change models were significant (P < 0.001). Chill days were also decreased in the 20-year periods between 1960 and 2015 (Table S5).
The number of backlashes (Index 3) significantly decreased with all climate change scenarios (P < 0.001), except for the UKCP09 model, where backlashes increased by 70 per cent ( Figure 3). All the models that used higher CO 2 emissions scenarios showed a decrease in backlashes, with the largest decrease in higher CO 2 emission scenarios ( Figure 3). The CMCC model estimated a greater decrease in backlashes than the UKCP18 model, and the downscaled version of UKCP18 showed a smaller change in backlashes than the global model ( Figure 3). All differences between climate change model results were significant (P < 0.001). This follows the historical trend, where the number of backlashes decreases between 1960-1979 and 1980-1999, although the decrease slows between 1980-1999 and 2000-2015 (Table S6).

Geographical differences within GB
All values showed a significantly (P < 0.01) high level of spatial autocorrelation for the historical  dataset (Table S7). All indices, except Index 2 min , were significantly spatially autocorrelated for UKCP09, as were all indices, except for Index 2 max for UKCP18. For the global model of UKCP18, fewer indices were significantly spatially autocorrelated, whereas for the CMCC global model, both the RCP 4.5 and RCP 8.5 had fewer significantly spatially autocorrelated indices and the autocorrelation was weaker.
The climatic gradient of temperatures did not change with climate change, with the coldest provenance region being the northeast (zone 20), followed by the northwest (10), southwest (30) and southeast (40) of GB (Table S4).
However, geographical differences in the number of frosts below −3 • C (Index 1 -3 • C ) between the onset of dehardening and late summer between provenance regions generally reduced with climate change, especially with the higher CO 2 emissions scenario ( Figure 4). Geographic differences were more visible in the downscaled model when compared to the global UKCP18. Geographic differences in the incidence of frosts below −5 • C (Index 1 -5 • C ) were non-existent historically and were predicted to remain so with climate change.
The number of frosts <−3 • C (Index 1 -3 • C ) was reduced more in the southern regions, where the historical average number of frosts was lower, than in the northern regions ( Figure 4). There are also longitudinal differences in the rate of change. The northeast showed a higher rate of decrease of frost incidence than the northwest, whereas in the south of GB the southwest (provenance region 40, Figure S2) showed a higher rate of decrease than the southeast (provenance region 30, Figure S2).
The provenance regions were significantly different from each other in the number of frosts <−3 • C (Index 1 -3 • C ) they experienced, and the differences remained significant (P < 0.001) with climate change. The CMCC global model showed fewer pairs of provenance regions remaining significantly different (P < 0.001) from each other, with no differences between provenance regions 20, 30 and 40 for the RCP 8.5, and no differences between provenance regions 30-40 for the RCP 4.5.
For the historic period 1960-2015 for GB the geographical variation in the number of frosts after bud burst, whether early, average or late (Index 2 min , 2, and 2 max ), was low but significant (P < 0.001). Whereas for all climate change models, no significant geographical differences in the number of frosts after bud burst were observable ( Figure S2).  The geographical variation in the number of backlashes, Index 3, was found to increase with climate change, at least for the downscaled models ( Figure 5). Geographic differences were much less observable with the global models.
The same pattern of change was observed for Index 1 and Index 3, with differences in the rate of change across longitudinal and latitudinal clines (Figures 2 and 3).
The hastening of bud burst occurred more markedly in the south than in the north, although differences in the rate of bud burst change between east and west were not apparent. Geographical differences in the date of bud burst remained significantly different (P < 0.001) both on the longitudinal and latitudinal cline (Table S5).

Validation of indices
Frost damage observed by commercial foresters in the field was most often correlated with a threshold temperature of −3 • C. This suggests that the natural hardiness of Sitka spruce trees to daily minimum temperatures between 0 • C and −3 • C is sufficient to prevent damage, with the bifurcation threshold where damage occurs being close to −3 • C. The weak correlation between frost damage and the number of frosts below −5 • C in British commercial forests, despite the increased severity of the frosts, can be explained by the rarity of these events, as GB has a mild, temperate climate. Thus, it can be expected that frosts cause severe damage more frequently due to unusual timing than due to their severity.
The correlation of tissue damage with the number of times minimum temperatures drop below the natural frost hardiness of a tree (i.e. backlashes; Index 3) suggests that Index 3 can be used to predict expected levels of damage due to frosts. Whilst we have shown that Index 3 has utility in frost prediction, it should be noted that the data used for validation was based on frost damage observations made in a single year and that correlation of frost damage observation with indices does not necessarily indicate causation. For example, a single severe frost event could have caused all the observed damage, although this is unlikely, because forest managers were unable to ascribe the damages to a single event. This suggests that an additional measure of frost severity could be useful in making predictions.

Effects of climate change
Climate change, with warmer winter and spring periods, will change the conditions under which frosts occur in GB. Our modelling shows that predicted climate change will result in a decrease in the number of backlashes (Index 3) during the entire year thereby decreasing the chance of damage to Sitka spruce (Figure 3). Less frost damage is likely despite the delay of the Forestry onset of hardening in autumn, which could decrease maximum winter frost hardiness thresholds by as much as 7 K (Guak et al., 1998).
The total number of frosts (Index 1 -3 • C ) between the onset of dehardening and late summer is also predicted to decrease with climate change. However, the number of frosts matters less than their timing. Our data showed that climate changedriven increases in spring daily mean temperatures (average temperatures above 5 • C) will advance dehardening and bud burst, resulting in conifers that are at their lowest level of frost hardiness earlier in the year (Cannell and Sheppard, 1982;Cannell and Smith, 1983;Ma et al., 2019;Morin et al., 2009Morin et al., , 2010. Frosts that occur before trees have dehardened do not affect trees severely, as harm occurs when air temperatures are below the frost hardiness level of the tree tissues. Therefore, the timing of the frosts is very important, and frosts occurring after a tree has dehardened or during bud burst are expected to be much more harmful (although rarer) than those occurring when a tree is hardened.
Despite the overall significant decrease in chill days during the winter with climate change, the simultaneous significant increase in thermal time, as measured by GDD (Table S4), is expected to lead to earlier bud burst dates for Sitka spruce across GB, as calculated by equation 1 (Table S4). Bud burst is expected to occur 2-57 days earlier, according to the downscaled climate models. It should be noted that the error of bud burst dates is about 12 days on the national scale (Peaucelle et al., 2019b). Cannell and Smith's (1983) bud burst model is based exclusively on GDD and ignores the difference between the effects of daily maximum and minimum temperatures. Despite this omission the GDD based model seems to accurately describe bud burst dates for the 1960-2015 historical period used in this study (Olsson and Jönsson, 2015). However, our study found that accounting for the different effects of maximum and minimum temperatures would lead to a model too specific to a site to be useful at larger resolutions, where generic models tend to work best (Peaucelle et al., 2019b). Thus, the generic model works best for the scale that is being used, especially considering it is calibrated to the most common local origin of Sitka spruce (QCI).
Historically, the number of post-bud burst frosts precited using the Cannell and Smith (1983) GDD model is low throughout the entire 1960-2015 period, with extremely long return times between post-bud burst frosts (Table S3). Our GDD model is specific to Sitka spruce varieties introduced to GB. It is not locally calibrated to each seed zone as: (I) Sitka spruce is an introduced species; (2) most commercial Sitka spruce varieties are of the same origin; and (3) local adaptations have not occurred yet due to the limited time scale. However, due to the extremely low value for Index 2 of post-bud burst frosts (an average occurrence of 0.05 post-bud burst frosts a year), it is expected that changes in phenology would not affect the amount of frost damage in commercial forests in GB (Table S3). Even considering early bud burst varieties of Sitka spruce (Index 2 min ), the number of postbud burst frosts remains extremely low both historically and with climate change.
Thus, despite the inaccuracies that result from the lack of calibration of the phenological models to local adaptations, the size of the effect means that, for Sitka spruce in GB, climate change is unlikely to change the number of post-bud burst frosts by orders of magnitude to a level where it has consequences for commercial plantations. Occurrences of post-bud burst frosts (Index 2) with climate change were significantly (P < 0.001) different to the number of historical frosts, but, at the UK scale, the extremely low frequency of frosts means that they are unlikely to cause noteworthy damage to tree form and productivity. Thus, whilst the damages inflicted on Sitka spruce plantations by postbud burst frosts may be significant when they do occur, the rarity of modelled return times (55-80 years; Table S3) suggests that post-bud burst frosts are unlikely to be a factor that should be considered in relation to commercial forestry. Although the rotation length of a Sitka spruce forest in GB is typically 35-45 years, frosts are most damaging in the first 3-4 years of the establishment of a young plantation, after which damage is unlikely to affect tree form and timber quality. Due to the rarity of such frosts the length of the return time calculated in this study is limited to the length of the period studied and is likely to be even longer than predicted.
Phenology-based predictions of frost damage with climate change are very common with many authors considering the frequency of frosts after bud burst for indigenous tree species (Jönsson et al., 2004;Ma et al., 2019;Morin and Chuine, 2014) or for all vegetation (Rigby and Porporato, 2008). It should be noted that post-bud burst damage only affects some tree species (Ma et al., 2019), and our study shows Sitka spruce is not affected. Since phenology is such a common concern for breeders when it comes to frost damage, it is useful to highlight that efforts to alter Sitka spruce phenology are unlikely to reduce frost damage in plantations throughout GB.
Locations where late spring frosts are more common, such as North America, are typically populated by tree species that are late-leafing, so these should suffer from less frost damage, whereas Europe and Asia have native species that quickly react to warm temperatures, and these will thus in the future become more vulnerable (Zohner et al., 2020). This suggests that Sitka spruce, which is of North American origin, is likely to benefit from climatic change when compared to native species in GB in terms of frost tolerance.
Bud burst date is negatively correlated with injury from frost, which means an earlier bud burst date leads to more damage (Aitken and Adams, 1997;Fløistad and Granhus, 2010). However, the number of frosts that occur after bud burst (Index 2) is significantly reduced with climate change (Table S2). This means that, despite bud burst happening earlier, the levels of harm from frosts do not increase as trees are most sensitive in the period preceding and following bud burst when the number of frosts happening during this most critical period is expected to decrease.
These anticipated reductions in frost harm will depend on whether the warming caused by climate change is uniform. A warmer climate during the winter and spring without a reduction of frosts would increase the chances of potential frost damage experienced by trees. Despite this, the onset of plant growth following dehardening does not prevent rehardening, and fluctuations in temperature can result in cycles of hardening and dehardening (Leinonen et al., 1997). Plants cycling through the different physiological states of frost hardiness have also been shown to have higher levels of frost hardiness than plants grown Impact of climate change on the occurrence of frost damage in Sitka spruce in consistent temperatures, as long as the average temperature is the same (Leinonen et al., 1997).
There is an overall decrease in the occurrence of Index 3, which calculates the number of days with minimum temperatures below the frost hardiness level achieved during daily hardening-dehardening ( Figure 3, Table S3), following the historical trend (Table S6). However, it should be noted that this reduction is based on the absence of sudden local fluctuations that cannot be predicted by big-scale climate models.
An observational study in northeastern Ontario showed that warming in that area was not as uniform as was predicted by previous modelling studies. There was an increase in extreme events, like a sudden −8 • C frost which was preceded and followed by daily maximum temperatures up to 30 • C (Man et al., 2009). Therefore, we cannot rule out wider fluctuations in temperatures without the general reduction of frosts in some areas in GB where overall increases in temperature lead to an increase in sensitivity, which would mean that in those areas, the levels of harm could increase.
Changes that accompany warming, such as an increase in dry seasons (Gazol et al., 2019;Vitasse et al., 2019), may also counteract the overall beneficial effects of warming. Thus, drought years will harm the growth of trees due to the reduction of xylem conductivity (Vitasse et al., 2019). Frosts can also induce xylem embolism (Gazol et al., 2019), in addition to harming needles.

Geographical differences in GB
All these climate improvements are happening to a greater degree in the north of GB relative to the south, which would facilitate an increased productivity of plantations in the north, through lower risk of frost during establishment of young plantations and increased productivity and yield due to warmer growing periods. Decreases in subfreezing temperatures are expected to lead to increases in yield of only 10-15 per cent for Sitka spruce in GB (Waring, 2000). This is because the main limitations to growth in Sitka spruce plantation in GB seems to be solar radiation and soil fertility more than the length of the bud burst to bud set period. In fact, even winter growth is limited by solar radiation rather than temperatures for Sitka spruce in GB (Waring, 2000).
However, these decreases in the predicted number of harmful frosts depend on the reliability of climate models. Some studies show that models over-estimate mean temperatures, while better predicting increases in maximum temperatures (Pohl et al., 2019). There are also issues in prediction certainty because model predictions are based on the data from close meteorological stations, but local microclimatic conditions (e.g. frost hollows) may create conditions for frost damage. Thus, estimated temperatures are a proxy for real temperatures in certain locations, which are unknown (Keenan et al., 2020). This means that even if an area is predicted to have decreased frost events, local conditions may contribute to higher than expected levels of harm from frosts (Man et al., 2009).
Autocorrelation measures the level to which the value of a variable depends on its neighbouring values. If the autocorrelation is significant, the spatial distribution of high values will be spatially clustered with other high values, contrary to what would occur if the spatial processes were random. This means that the underlying geography significantly correlates with the values obtained. Thus, different regions of GB will have different spatial trends based on their geography. Global models that are not downscaled have much lower levels of spatial autocorrelation. This would indicate that low resolution global models are unable to account for underlying geography of smaller regions, thus ignoring important regional differences. The UKCP18 Regional Projections model used in this study uses RCP 8.5, and the UKCP09 Met Office HadRM3-PPE UK run for the HadRM3 model uses the SRESA1B scenario. Improvements to the modelling accuracy of this study could be achieved if datasets of Tmin and Tmax generated from local scale (grid size ≤ 25 km) models with different RCPs were available to estimate future effects of climate change.

Conclusion
The overall effect of climate change on Sitka spruce in GB is predicted to result in earlier bud burst, a decrease in the occurrence of frosts, and a decrease in the occurrence of frost events during a vulnerable period for plant tissues, known as backlashes. The lengthening of the bud burst to bud set metabolically active period due to earlier bud burst and the decreased risk of harm from frosts, are likely to increase the productivity of tree plantations, improve the quality of wood products, and could reduce establishment and management costs incurred by the commercial forest industry.

Supplementary data
Supplementary data are available at Forestry online, in both graphical and table format. Figure S1 shows the map of GB's seed zones. Figure S2 shows the absence of geographical patterns for Index 2. The data used to validate indices are presented in Table S1. Average values for the 20-year periods between 1960 and 2015 for Index 1 0 • C , 1 -3 • C , 1 -5 • C for each seed zone are presented in Table S2. Average values of all indices for each climate model are presented in Table S3. Average estimated bud burst dates for each provenance region and climate model are given in Table S4. Average estimated bud burst dates for each seed zone and each 20-year period between 1960 and 2015 are given in Table S4. Average number of backlashes (Index 3) for each seed zone during the 21-year periods between 1960 and 2015 are given in Table S6. Values on the level of autocorrelation between neighbouring seed zones are given in Table S7.