Hot summers ahead? Multi-decadal spring season warming precedes sudden summer temperature rise in pre-anthropogenic climate change

ABSTRACT Waning annual seasonality is documented in an up to one-month advance in spring onset since the 1980’s in northern latitudes, perturbing ecosystems and socio-economic performance. Summer temperatures, in contrast, have been rising only recently, indicating an offset in seasonal warming. The limited time span of this observational data makes the asynchronous pattern difficult to quantify, hindering projections of intra-annual dynamics. We explore temporal phase relations of seasonal warming over the Late Pleniglacial/Bølling and the Younger Dryas/Holocene climate transitions that preceded present anthropogenic warming. We determine past spring onset and thermal properties from dwarf birch paleo-phenology. Reconstructed spring warming led maximum summer warming by about a century during both transitions. Long-term reconstruction of intra-annual temperature regimes provides the perspective required for seasonal response analysis. Our results document that multi-decadal spring season warming precedes sudden summer temperature rise also during natural climate change. The rapidity of present seasonality changes, however, is unprecedented.


Introduction
Annual seasonality is diminishing in northern latitudes through rapid winter warming, opposed to more moderate summer warming (NASA GISS data 2019). Northern vegetation has responded with a significant advance in spring onset (defined by bud-burst and vegetation greening) exceeding a month in some regions since the 1980's (Menzel et al. 2006;Xu et al. 2013;Wang et al. 2015). Prolonged vegetation greening contributes to climate change through altered albedo and hydrological pathways; exemplified starkly by the Arctic amplification phenomenon (Parmesan et al. 2012;Chae et al. 2015;Overland et al. 2016). In addition, earlier spring onset causes various phenological asynchronicities between plants and insects, with cascading idosynchratic effects through terrestrial trophic ecosystems (Primack et al. 2009;Richardson et al. 2013;Thackeray et al. 2016;Kharouba et al. 2018). Recently, higher summer temperatures have been observed globally, particularly in mid to high latitudes (Hansen & Sato 2016), suggesting a temporal offset between warmer springs and warmer summers in the ongoing climate change (Fig. 1).
Hotter summers are forecasted to further exacerbate negative effects on ecosystems and human societies, including extreme weather events and wildfires, loss of habitat, decreased food security, spread of pests and disease, as well as sharp increase in mortality by exceeding human thermoregulatory capacity (Sherwood & Huber 2010;Polgar & Primack 2011;Hansen et al. 2012;IPCC 2014;Mora et al. 2017). These rising summer temperatures are a robust trenda "new normal"and not a mere product of weather variability inherent to the climate system (Coumou & Rahmstorf 2012;Trenberth et al. 2015;Hansen & Sato 2016;Lewis et al. 2017). The ecological and socio-economic consequences of this intimate phenologyclimate link are far-reaching, but extremely difficult to quantify or predict due to the severely limited time span of observational data. To place these emerging asynchronous spatial and temporal patterning of seasonality changes into a sufficiently broad climatological context, the analysis of long-term climate proxy records is essential.
Here we investigate two late-glacial natural warming episodes that each occurred within a few decades, thus resembling the ongoing global warming in rapidity and spatial expression (Steffensen et al. 2008;Heiri et al. 2014). These are the transitions from the Younger Dryas to the Holocene (YD/H:~11.800 years before present = 11.8kyr BP) and the Late Pleniglacial to the Bølling (LP/B:~14.8kyr BP). For these previous, fully completed, warming episodes we present high-resolution paleo-phenological records derived from image-analysis of sub-fossil leaves of the arctic shrub Betula nana (dwarf birch). The leaf cuticle morphology records growing season thermal conditions expressed as Growing Degree Days (GDD 5 : the cumulative seasonal sum of degree Celsius (°C) above a threshold value of 5°C for subarctic biomes) during the growth period and correlates to budburst dates indicative for spring onset (Wagner-Cremer et al. 2010). In combination with independently reconstructed summer temperatures, seasonal warming patterns are quantified, their phase relations deciphered and compared to the ongoing climate development. So far, winter vs summer season-resolved paleo-temperature data from single-source material are available only from Greenland ice-core records (Vinther et al. 2010) and from rare palynological studies . Such records, however, reconstruct the individual temperature dynamics without synthesising intraannual seasonal dynamics. Increasingly applied, dataconstrained spatio-temporal paleoclimate models stress the importance of diverging seasonal warming during the last deglaciation (Buziert et al. 2014(Buziert et al. , 2018. Model comparisons to temperature reconstructions based on proxy-data often implement summer temperatures only, and suffer from inconsistencies in regional analysis of seasonal temperature changes (Zhang et al. 2016(Zhang et al. , 2017Schenk et al. 2018). We here present long-term data such as required for full-extent seasonal response analyses in the terrestrial realm. These data include morphological measurements of plant leaf cells, compared with new and previously published pollen and chironomid reords. Importantly, our seasonal reconstructions are based on multi-proxy analysis from single sediment sections, overcoming numerous error sources relating to age determination, chronology and between-site variability. By explicitly quantifying spring-onset related phenological traits such as bud-burst date, we directly compare past and present vegetation response to seasonal climate change and capture the early season dynamics.

Localities and leaf databases
The first Betula nana image databases used in this study derives from a Late Glacial paleo-lake sedimentary sequence of the Hässeldala Port (HÄ) locality located in southern Sweden (Fig. 1). Previous studies have reconstructed the stratigraphy and chronology through pollen, tephra and 14 C data analyses, as well as the paleo-climate and environment with multiple biogeochemical proxies including a chironomid-based temperature record and climate model simulations (Steinthorsdottir et al. 2013(Steinthorsdottir et al. , 2014Muschitiello et al. 2015;Wohlfarth et al. 2017). The HÄ database encompasses the period 12.5-11.75kyr BP and thus includes the abrupt warming at the end of the Younger Dryas.
The second Betula nana database derives from leaves preserved in 15.2-14.4kyr BP old sediments from Schleinsee (SCH) in southern Germany ( Fig. 1), spanning the termination of the last glacial and the onset of the interstadial (i.e. the Bølling-Allerød warm phase). The chronology for the SCH sediment section was obtained by 14 C data analyses and a bulk carbonate stable oxygen isotope (δ 18 O) record that is tied to the Greenland NGRIP δ 18 O record (Wagner-Cremer & Lotter 2011). The July temperature record has been reconstructed based on pollen from the same sediment section (Fig. S2).

Micromorphological analysis of fossil Betula nana leaf samples
Microscopic analysis of micromorphological epidermal cell parameters has been performed on all available leaf remains ranging from a total number of 2-12 subfossil leaf fragments in each sampled horizon of 1 cm core section at both sites totalling 101 leaf fragments deriving from HÄ and 53 from SCH. Sample availability strongly depends on abundance and cuticle preservation in the sediments. Analysis of the leaf-area margin was avoided to minimize errors induced by the intrinsic variability of the leaves; species determination is based on species-specific stomatal size (Wagner et al. 2000). Micromorphological epidermal cell characteristics are subsequently converted to GDD5 and bud-burst dates using transfer functions deduced from modern training sets as described below. Legend: Temperature changes expressed as degree Celsius/100 year trends based on observational data using the CRU TS 4.00 data-set via the KNMI explorer (climexp.knmi.nl). White stars indicate locations of the sites Hässeldala (HÄ), Sweden and the Schleinsee (SCH), Germany.

Undulation index (UI)
Ultimate leaf size is a function of initial epidermal cell division and successive lateral epidermal cell expansion leading to the final leaf size at the end of the growing season (Wagner-Cremer et al. 2010). The waxy cuticles of leaves reflect epidermal cell properties as size and shape and are resistant to decomposition and diagenetic processes during burial in sediments. Epidermal cell areas (CA [μm 2 ]) and epidermal cell circumferences (CC [μm]) of minimal 30 cells per leaf fragment, after which standard deviations become constant, were analyzed using computer aided image analysis of microscopic photographs taken from modern and fossil leaf material using either transmitted light or fluorescence microscopy (Fig. 2).
From CA and CC values, the epidermal cell-wall undulation index: was calculated for each leaf fragment, following Kürschner (1997). The UI measures the deviation between the circumference of a cell (CC) with a given cell area (CA) and the circumference of a circle with the same area.

GDD 5 and bud-burst date inference models
Routinely applied leaf cuticle morphology analysis of the (sub-) arctic canopy species Betula nana from the phenological station at Kevo, Finland, has shown that the thermal properties of the growing season (April-September, AMJJAS), expressed as GDD 5 , produce a quantifiable imprint in leaf morphology (Wagner-Cremer et al. 2010). The UI of Betula nana in the modern training set is linearly correlated to the range of 450-1050°C GDD5 that occurred during the past 20 years of continuous monitoring (Wagner-Cremer et al. 2010) (Fig. 2;  Fig. S1). Moreover, a significant correlation between UI and the date of bud-burst allows also for concomitant estimates of the timing of phenological spring onset (Wagner-Cremer et al. 2010) (Fig. S1).
The mean UI values subsequently determined from fossil leaves are used to infer GDD 5 and bud-burst dates by applying the species-specific prediction models based on modern training sets (Wagner-Cremer & Lotter 2011) (Table 1).

Reconstruction of spring onset and seasonal temperature dynamics
The HÄ reconstruction ( Fig. 3A; Fig. S2), beginning at 12.5kyr BP, records GDD 5 of 650-680°C and reaches minimum values of 640°C between 12.2kyr BP and 12.1kyr BP, after which GDD 5 increase to maximum values of 800°C at 11.8kyr BP. The chironomid-inferred July temperature over the record at HÄ shows minor temperature fluctuations centred around 9.2°C for the interval from 12.4kyr BP to 11.95kyr BP, where a rapid increase in July temperature cumulates in maximum summer temperatures of 11.4°C at 11.75kyr BP (from Muschitiello et al. 2015). Summer temperature reconstructions are commonly provided as July temperatures, but here we present a direct comparison between our spring onset data and maximum summer temperatures by providing bud-burst data as days before mid-July (as maximum temperature tiepoint, which we refer to from here on as maximum summer warmth: MSW), to avoid potential differences in season timing between the Late Quaternary and today. From 12.5kyr BP to 12.25kyr BP spring onset occurs~43 days before MSW, then retreats to a minimum of 39 days before MSW until 12.1kyr BP, after which spring onset advances to a maximum of 53 days before MSW at 11.8kyr BP. Besides the clear dynamics of the growing season length, we determine a distinct 130 year offset between spring onset and MSW towards the end of the Younger Dryasthat is, it takes MSW~130 years to catch up to the early spring onset (see Fig. 3A).
The SCH dataset ( Fig. 3B; Fig. S2) begins at 15.12kyr BP where Late Pleniglacial low GDD 5 values of 665°C are recorded, followed by a successive increase to~850°C prevailing from 15kyr BP to 14.75kyr BP. A sudden jump to maximum GDD 5 of 1000°C is completed within less than a century where high GDD 5 prevail until 13.2kyr BP. SCH summer temperatures are estimated from pollen analysis revealing relatively constant low MSW values of~10.5°C between 15.2kyr BP to 14.7kyr BP, from where a continues increase leads up to MSW of 13.3°C at 14.4kyr BP. Spring onset at SCH occurs 33 days before MSW at 15.2kyr BP, and a two-step advance to~40 days before MSW between 15kyr BP and 14.8kyr BP from where a rapid increase leads up to a maximal spring advance to 50 days before MSW at 14.75kyr BP, prevailing for the remaining record. Importantly, also here we detect a distinct 110 year offset between spring onset and MSW towards the end of the Late Pleniglacial, with MSW reaching a maximum~110 years later than early spring onset. Note that the overall much higher GDD 5 for SCH compared to HÄ reflects the latitudinal temperature gradient between the two sites.

Discussion
For both studied sections, a spring advance of 1-1.5 days/ decade is calculated for the intervals of most rapidly changing bud-burst dates. At SCH the advance seems to be slightly more rapid than at HÄ, however, a more detailed analysis requires a higher temporal resolution than is available for our sections. This 1-1.5 days/decade spring advance is significantly slower than the modern spring greening advance, shown by remote sensing data-based calculations to be between 3 to more than 12 days/decade since the 1980's, in the Eurasian arctic to boreal realm (Menzel et al. 2006;Xu et al. 2013). Betula nana, however, has been shown to emerge late in the modern plant community bud-burst succession (Post et al. 2016), which can explain the discrepancy with satellite monitoring of the entire vegetation in an area. Our results thus either pick up a potentially weakened signal, or illustrate the unprecedented pace of the Anthropogenic climate change in the northern latitudes.
Our paleo-records indicate an offset between spring and summer warming during natural phases of climate change, where rapid spring warming leads by approximately a century. This divergence in seasonal response to climatic triggers constrains modelled asynchronous responses determined for the warming episodes during the last deglaciation (Buizert 2015). Such temporal and spatial response estimates suggest that over Greenland, intensive winter warming exceeds the ultimate summer and, thus, mean annual temperature shifts (Buziert et al. 2018;Salonen et al. 2018). Asynchronous warming over Greenland is attributed to reorganization of Atlantic meridional overturning circulation, freshwater forcing and atmospheric CO 2 dynamics (Buizert 2015;Buziert et al. 2018). Data-constrained simulations of July temperatures over the B-YD time interval show that relatively high summer temperatures prevail during the YD over large parts of the northern hemisphere, where the main temperature change is explained by intensification of Table 1. GDD 5 and bud-burst inference models. Inference models for GDD 5 and bud-burst dates based on double log transformation of modern training sets of Betula nana UI values from Kevo and climate station data as well as phenological records from the same site. GDD 5 inference model GDD 5 = 10 (2.4232+6.0284*log(UIfossil)) , r 2 = 0.681, p = <0.001, RMSE = 62 GDD 5 Bud-burst inference model  A. Green line -GDD 5 reconstruction based on fossil Betula nana leaf fragments from the Hässeldala site (Sweden) with RMSE as grey envelope. Red line: MSW (July temperature (°C) reconstructed by chironomid analysis from the same sediment core, 3-point average of data (Muschitiello et al. 2015). Green shaded area is the onset of vegetation greening deduced from reconstructed Betula nana bud-burst dates as number of days before MSW for comparison to proxy-based summer temperature. Betula nana leaf pictures visualise the relative size changes in foliage under varying thermal conditions. Stippled lines delineate onset of GDD 5 increase and initial increase in July temperatures in the proxy records. B. GDD 5 reconstruction based on fossil Betula nana leaf fragments from the Schleinsee site (Germany) with RMSE as grey envelope. Red line: MSW (July temperature (°C) reconstructed by pollen analysis) from the same sediment core. Green shaded area and stippled lines show onset of greening and offset between GDD 5 accumulation and July temperature (°C) for the Late Pleniglacial -Bølling transition quantified as in panel A.
continentality and, thus winter to spring cooling (Schenk et al. 2018). Our seasonally resolved temperature reconstructions underpin the relevance of winter and successive spring warming in annual temperature dynamics also during rapid warming phases. Uniquely, we are able to quantify the extremely dynamic seasonal warming outside of Greenland, where ice-core based modelling efforts may exceed their prediction limits. Such records from the terrestrial realm are thus a prerequisite to document the spatio-temporal phase relations of seasonal temperature over the large geographical realm undergoing rapid climate change, after oceanic and/or atmospheric reorganizations. Decoupling of seasonal temperatures during other geological periods which may serve as potential analogues for present climate change, such as the Eemian interglacial, has recently been highlighted by proxy-based temperature reconstructions, where the mixed imprints of insolation and oceanic forcing regulates seasonality expression in the high northern Latitudes (Salonen et al. 2018). The link to potential CO 2 forcing during the Eemian, however is low. For the more recent climate episodes studied here, CO 2 forcing needs to be considered (Steinthorsdottir et al. 2013(Steinthorsdottir et al. , 2014. The CO 2 dynamics during past rapid climate change as well as warmer-than-present episodes are related, however, to natural feedback mechanisms whereas ongoing warming is clearly related to the anthropogenic CO 2 increase (IPCC 2014) and potentially develops alternative feedbacks. In the context of the ongoing climate change, where close to 40 years of spring warming has already occurred, assuming that the Arctic behaves similarly during the current warming to the two warming transitions we have shown for the late glacial, then we can expect that the rise in summer temperature will occur within 60 years. Given the pace of the present warming however, with at least twice the rate of spring onset advance measured in days/decade, as well as recent documented shifts in mid-low latitudes towards hotter summer temperatures (Hansen & Sato 2016), it is well possible that this threshold will be reached much sooner.