Spring onset and seasonality patterns during the Late Glacial period in the eastern Baltic region

. Spring onset is an important phenological obser-vation that is sensitive to modern climate change and can be traced back in geological time. The Late Glacial ( ∼ 14 500– 11 700 cal yr BP) spring onset and growing season (growing degree days) dynamics in the eastern Baltic region were reconstructed using the micro-phenological approach based on the dwarf birch ( Betula nana ) subfossil leaf cuticles. The presented study sites, Lake Lielais Svetinu (eastern Latvia) and Lake Kosilase (central Estonia), are located ∼ 200 km apart in the region affected by the south-eastern sector of the Scandinavian Ice Sheet. During the Late Glacial period the region and its biota were inﬂuenced by the retreating glacier and the different stages of the Baltic Ice Lake. The plant macrofossil data conﬁrm that the study sites were in different vegetation zones (arctic-to-boreal) during the Late Glacial period. The dynamics of the estimated length of the growing season and spring onset, combined with the regional collection of plant macrofossil records, suggest the importance of local settings to species migration. During the Late Glacial warming period (Bølling–Allerød), a notable spring warming and longer growing season was calculated based on micro-phenology, but the treeline did not extend beyond central Estonia. The comparison of pollen- and chironomid-inferred past temperature estimations with spring onset, growing degree days, and plant macrofossil data shows coherent patterns during the cooler Older Dryas and warmer Bølling–Allerød periods, while suggesting more complicated climate dynamics and possible warmer episodes during the Younger Dryas cold reversal.


Introduction
Changing seasonality belongs to the most eminent characteristics of ongoing global climate change. Phenological observations are sensitive and easily obtainable indicators of biospheric changes in response to climate change (Peñuelas and Filella, 2001;Badeck et al., 2004;Cleland et al., 2007). The earlier unfolding of leaves has been observed since the mid-20th century all over the Northern Hemisphere (Ahas et al., 2002;Badeck et al., 2004 and references therein;Menzel et al., 2006;Jeong et al., 2011;Montgomery et al., 2020). The spring phenology is influenced by pre-season temperature that is described as temperature before the spring phenological date, winter temperature, and accumulated precipitation (Wang et al., 2015). These seasonality changes are the subject of intensive discussions, given the strong impact of the earlier spring onset on many sectors, including agriculture, ecosystem stability, and vegetation dynamics (e.g. Buermann et al., 2018;Menzel et al., 2020).
Palaeoecological studies have the potential to contribute to a better understanding of spatio-temporal seasonality dynamics by reconstructing seasonal temperature changes during the phases of natural climate change that occurred in the recent past. While past summer and winter temperatures have been reconstructed more often, the dynamics of the spring season is still underrepresented in the temperature records available so far. The micro-phenological proxy based on the dwarf birch (Betula nana) subfossil leaf cuticle features (Wagner-Cremer et al., 2010) offers a means to shed light on past seasonality patterns, the dynamics of spring onset, and the amount of growing degree days (GDDs). This proxy makes use of the direct correlation between the lateral epidermal cell expansion that regulates cell size and shape during the maturation of leaves (Wagner-Cremer et al., 2010;Ercan et al., 2020). The longer the growth period available to plants and higher the accumulated GDDs, the larger and increasingly undulated the epidermal cells grow. Quantified as the undulation index (UI) in long-term monitoring studies and validated in free-field growth experiments (Ercan et al., 2021), this proxy has demonstrated asynchronous spring versus summer temperature dynamics during the warming from the Late Pleniglacial to the Bølling-Allerød period and during the transition from the Younger Dryas to the Holocene (Wagner-Cremer and Lotter, 2011;Steinthorsdottir and Wagner-Cremer, 2019). During both warming phases, spring onset and GDD accumulation rose more than a century before the summer temperatures started to increase (Steinthorsdottir and Wagner-Cremer, 2019). A systematic application of this proxy to B. nana leaf-bearing lake and peat sequences thus has the potential to add information on the expression of spring season warming over a large geographical range.
In the present study, we focus on the seasonality changes in the eastern Baltic region by compiling multiple temperature and vegetation reconstructions based on a variety of biological proxies. The Baltic region was covered by the south-eastern sector of the last Scandinavian Ice Sheet (SIS) (Kalm, 2012). During the past decade, the Late Glacial environmental history of this region has been intensely studied (e.g. Stančikaitė et al., 2009;Heikkilä et al., 2009;Veski et al., 2012;Druzhinina et al., 2020;Šeirienė et al., 2020). The Late Glacial changeable climatic conditions (Rasmussen et al., 2014), past summer and winter temperature reconstructions , and associated arctic-to-boreal vegetation dynamics Veski et al., 2012;Amon et al., 2016) and tree-line presence  provide an excellent case to study the seasonal dynamics of rapid natural climate warming and cooling episodes. Here we apply the micro-phenological UI proxy to two sites that hold uncommonly continuous and well-preserved B. nana subfossil leaf fragments in order to determine spring onset dates and the thermal properties of the growing season. The first indication that the UI proxy is also applicable in a more continental settings has been provided in experimental studies on B. nana relict stands in Poland (Ercan et al., 2021) and is applied here to a subfossil leaf record from the sediment core of the continental (eastern Baltic) region for the first time. Compared to the already available alternative temperature proxies for the same sites, this approach uniquely enables the first detailed analysis of seasonality changes during the Late Glacial period in the eastern Baltic region.

Study area
The topography of the eastern Baltic has been largely shaped by Pleistocene glaciations, where Weichselian glaciation (SIS) contributed to present-day topography. The region is currently situated in the hemiboreal vegetation zone, within the transitional zone between the boreal and nemoral forest zones of Europe. Estonia and Latvia are located between 56 and 59.5 • N on the eastern coast of the Baltic Sea in a transitional area from maritime to continental climate, characterised by a west-east gradient in the continentality of the climate. Climatic conditions are mainly determined by high cyclonic activity and the prevalence of westerlies (Jaagus et al., 2010). The mean annual temperature in Estonia varies between 4.1 and 6.5 • C (Riigi Ilmateenistus, 2021). The mean annual temperature in Latvia is 7.2 • C (Krauklis and Draveniece, 2004). The snow cover in Estonia persists from December to late March (Jaagus, 1997) and in Latvia from December to March-April (Draveniece, 2009;Krauklis and Draveniece, 2004).
The vegetation types during the Late Glacial period in the eastern Baltic region spanned from arctic tundra in North Estonia to mixed forests in Latvia . The vegetation dynamics are influenced by Late Glacial hemispheric climate fluctuations (Rasmussen et al., 2014), as well as by local factors . In the Late Glacial period (14 500-11 700 cal yr BP), an important palaeogeographical feature of the deglaciated eastern Baltic region was the formation of ice lakes (Vassiljev and Saarse, 2013;Rosentau et al., 2009;Amon et al., 2014). The vegetation communities during the Late Glacial period in the study region, as described by plant macrofossil and pollen records, spanned from pioneer snow-patch tundra to mixed boreal forests (Amon et al., , 2016. The present study explores the rich subfossil dwarf birch (B. nana) leaf records from two sites, namely Lake Kosilase and Lake Lielais Svetinu (L. Svetinu) ∼ 200 km southeast of Lake Kosilase (Fig. 1).

Method
The subfossil leaf material was extracted from the sediments during the plant macrofossil analyses. The sediments were wet-sieved and the remaining material was identified under a stereo-microscope and a light microscope. The subfossil B. nana leaves were stored in distilled water at 4 • C until a further analysis with the fluorescence microscope. Leaves with preserved cuticle cells were photographed under the microscope (Leica) using fluorescent light at 200× magnifi- cation. Epidermal cell properties were measured on digital images using the software ImageJ. From the measured cell area (CA, µm 2 ) and cell circumference (CC, µm), the undulation index (UI, dimensionless) was calculated according to Kürschner (1997). The UI was determined for a minimum of three epidermal pavement cells on at least three leaf fragments per sample.
Budburst dates are commonly used to determine spring onset. Spring onset dates are subsequently expressed as dayof-year (DoY) before midsummer warmth (MSW), referring to the number of days before 15 July (DoY: 196), which is set as the date for maximum summer warmth in order to facilitate a comparison with summer temperature proxies such as the chironomid-based July temperature reconstructions (Steinthorsdottir and Wagner-Cremer, 2019). This procedure enables the quantification of the time period between spring onset and maximum summer temperatures by calculating the DoY data rather than providing numbers for individual months.
From Lake Kosilase a 30 cm long silty sediment section was analysed and two samples were selected for 14 C dates. The chironomid-based temperature reconstruction from Lake Nakri was published previously (Heiri et al., 2014), as were the pollen-based temperature reconstructions  and the plant macrofossil data  from Lake L. Svetinu.

Results
Radiocarbon dating results for Lake Kosilase and Lake L. Svetinu are given in Table 1 and Fig. 2. The radiocarbon dates were calibrated using Oxcal 4.4 (Bronk Ramsey, 2009) and the IntCal20 calibration curve (Reimer et al., 2020).
The mean UI from Lake L. Svetinu subfossil leaf record varies between 1.17 and 1.25 around a total average of 1.2, with an average standard error of 0.03. The UI values of Lake Kosilase subfossil dwarf birch leaves range from mean 1.17 to 1.23 with a total average of 1.21 and a standard error of 0.04 (Fig. 3). The estimated GDD5 ranges from 674 to 1009 ( Fig. 4). The spring onset calculations were based on inferred budburst dates recorded by the fossil leaves and revealed spring onset dates between 21 May and 6 June, equivalent to DoY 146-157 and translating to 39-55 d before MSW (Fig. 4).

Discussion
The combination of (palaeo-)botanical and phenological records reveals the dynamics of the vegetation, growing season duration, and spring onset during the Late Glacial period in the study area (Ercan et al., 2021). The captured spring onset signal is probably regional. We have a temporal overlap of the two spring onset records from the study sites ∼ 200 km apart. The sample density in the overlapping part is not high, but the general trend (earlier spring onset around 13 900-14 000 cal yr BP, later spring onset around 13 800 cal yr BP, and earlier spring onset again around 13 600 cal yr BP) observed in both datasets suggests comparable patterns of the spring onset signal over the entire study region.
In the present study, we present spring onset and GDD5 records based on the UI analysis of B. nana leaves from two study localities ∼ 200 km apart (Fig. 1). Both lakes are currently surrounded by hemiboreal vegetation, but during the Late Glacial period the floral situation was much more dynamic. During that time, Lake L. Svetinu and Lake Kosilase were located in vegetation zones that oscillated northwards and southwards, following changing climatic conditions and species migration. The pioneer vegetation phases at both localities were characterised by tundra dominated by treeless environments, B. nana, and Dryas octopetala. During the Allerød, warming permitted the migration of vari- ous tree species to eastern Latvia, mainly from the southern and south-eastern regions. The palaeobotanical record of Lake L. Svetinu consists of the macro-remains of birch, pine, and aspen, suggesting the occurrence of mixed forest  that was one of the most diverse vegetation communities during the Late Glacial period in the region . However, the proximity of the re-treating glacier and areas submerged by the Baltic Ice Lake hindered the northward migration of trees, and the regional Late Glacial maximum northern treeline was likely located in southern Estonia . This implies that even during the warmer phases of the Late Glacial, the surroundings of Lake Kosilase most likely remained treeless and were mainly covered by tundra vegetation. During the Younger   and Kosilase (KOS), pollen-based May-June-July-August (MJJA) and December-January-February (DJF) temperature estimations of Lake L. Svetinu , chironomid-based July temperature estimation of Lake Nakri (location 1 in Fig. 1 Dryas cooling (∼ 11 650-12 850 cal yr BP) the tree remains disappeared from the records of Lake L. Svetinu, suggesting the re-establishment of tundra vegetation and a southward regression of the regional treeline. Although such vegetation shifts are recognised in many sedimentary records, the driving role of growing season dynamics has so far hardly been considered, given the lack of temperature records for the early season. Here we have the unique opportunity to synthesise vegetation dynamics and link it to climatic changes through local macro-remain analysis, GDD5 and spring onset reconstructions from B. nana leaf remains, and more regional temperature signals for winter and summer deduced from pollen and chironomid data.
The length of the growing season is correlated with the UI of both B. nana (Wagner-Cremer et al., 2010) and Betula pubescens subsp. czerepanovii (Ercan et al., 2020). The estimated length of the Late Glacial growing season in the Baltic region varies between ∼ 670 and 1000 cumulative growing degree days (GDD5) (Fig. 4). It is comparable with current northern and mid-Finnish conditions (Ercan et al., 2020). The GDD5 estimations from the microphenological fossil dataset in Germany (Schleinsee) stop at 14 400 cal yr BP, when the present dataset from the Baltic region starts (Steinthorsdottir and Wagner-Cremer, 2019). The GDD5 estimations from Hässeldala, Sweden, (Steinthorsdottir and Wagner-Cremer, 2019) start at 12 500 cal yr BP and span until the start of the Holocene; they are in general lower than the Baltic GDD5 estimations; however, both display an interval of higher GDD5 values during the Younger Dryas cooling period. The Holocene onset is more pronounced in the Hässeldala dataset than in the Baltic (Latvia).
The shifting of vegetation zones and migration of the species is evident from the studied B. nana record -no B. nana leaf remains are recognised in the L. Svetinu core between 13 500 and 12 500 cal yr BP. Despite the absence of B. nana leaves during this period, the continuous plant macrofossil records document the persisting local flora and provide evidence for vegetation dynamics. The dominant plant species substituting the dwarf birch during the Allerød (∼ 13 250-13 900 cal yr BP) in Lake L. Svetinu were tree birch (Betula pendula) and pine (Pinus sylvestris) (Fig. 4). The occurrence of these two species provides additional information on the spring conditions through their contemporaneous growth requirements. Phenological observations of B. nana budburst dates from the Kevo research station in Finnish Lapland suggest that the number of days from spring onset to MSW varies between 25 and 56 d, with a mean of 42 d. A modern dataset for B. pendula from several study sites in Finnish Lapland indicates that the number of days from spring onset to MSW ranges from 40 to 62 d; average spring onset date is 51 d before MSW (Pudas et al., 2008). The modern observations of the budburst dates of P. sylvestris in the Inari region in Finnish Lapland (Salminen and Jalkanen, 2015) suggest the spring onset range of P. sylvestris to be from 50 to 75 d before MSW. Compari-son of the modern phenological characteristics to the plant macro-remain data thus indicates that the major shift in vegetation composition was related to spring onset dynamics, which might have taken place earlier in the warmer periods of the Late Glacial.
The notable shift towards earlier spring onset at ∼ 13 300-13 400 cal yr BP (at the end of GI-1c or Allerød) evident from Lake Kosilase coincides with the findings of pine macrofossils in the sediments of Lake L. Svetinu (Fig. 4). It suggests an ameliorating, warmer environment in this region. The micro-phenology of dwarf birch growing in the tundra or at the verge of the treeline (Kosilase) parallels the warming signal that supported the formation of mixed boreal forests in eastern Latvia, ∼ 200 km southwards (L. Svetinu).
The regional spring warming at ∼ 13 300-13 400 cal yr BP and its effect on past vegetation may be similar to the modern tundra greening process and earlier spring dates reported by modern phenological observations Montgomery et al., 2020 and references therein). The warmer climate, including spring, supported the northward shift of the treeline: pine was present in Latvia, and at ∼ 13 300-13 000 cal yr BP the scarce tree-birch macrofossils were found in the sediments of Lake Nakri (∼ 100 km southeast of Lake Kosilase, Amon et al., 2012). The plant macrofossil records suggest that the Late Glacial treeline probably did not reach beyond central Estonia (Amon et al., , 2016, although a single poorly preserved, probably tree-birch or hybrid, macrofossil was found in the Lake Kosilase sediments. At the same time, spring in central Estonia was already taking place much earlier as confirmed by the UI of the B. nana leaves and could have been suitable for tree species growth. Thus, the limiting factor for the advancement of the treeline may have been related to other local conditions in the partially submerged mosaic landscape, the proximity of the melting or retreating glacier, and the maximum summer temperatures. The modern distribution of the treeline-forming tree birch in Scandinavia follows the 13.2 • C isotherm for summer temperature (Odland, 1996) and the 10-12 • C summer temperature isotherm for treeline formation (CAVM Team, 2003). The chironomid-based July temperature estimates around Lake Nakri for this period were ∼ 12.7-12.8 • C (Heiri et al., 2014), and the pollenbased reconstructed mean May-June-July-August temperature range was ∼ 12.2-12.3 • C, which are close to the modern tree-line limits. Assuming that the summer temperature was one of the limiting factors for the northwards migration of the treeline, the results resemble the vegetation response under modern climate conditions. The spring onset thus occurs before a significant rise in the summer temperature. The observed warmer climatic conditions were subsequently disrupted by the onset of the Younger Dryas cold reversal. A similar set -warming of the spring preceding the warming of the summer temperatures -is observed by Steinthorsdottir and Wagner-Cremer (2019).

Late Glacial climate fluctuations and seasonality signal
Pollen-based temperature estimations for May-June-July-August (MJJA) and December-January-February (DJF) are available for Lake L. Svetinu for the Late Glacial period  and provide an indication of the winter and summer conditions prevailing in the studied region. A minor offset between the pollen record and the microphenological dataset is likely related to the different sampling strategies. The pollen samples were taken as 1 cm thick slices after some interval while the plant macrofossil samples from where the dwarf birch leaves were extracted were 5 cm thick continuous samples. Additional chironomid-based July temperature estimates from Lake Nakri (Fig. 4) ∼ 90 km to the north are included to complement the temperature data (Heiri et al., 2014).
In the oldest part of the record, between 14 400 and ∼ 14 000 cal yr BP, assigned to GI-1e (Bølling, ∼ 14 050-14 650 cal yr BP), all proxies (pollen-based DJF and MJJA temperatures, GDD5, and spring onset estimates; see Fig. 4) follow a comparable pattern of cool moderate temperatures during all seasons. Short pulses of cold conditions occurred between 14 300 and 14 400 cal yr BP and 13 800 and 13 900 cal yr BP. Since these coolings were short, the temporal resolution of our chronology is not suitable to determine the offsets in the seasonality changes between different proxies. The later cooling is interpreted to reflect the GI-1d (Older Dryas) phase (∼ 13 900-14 050 cal yr BP) (Rasmussen et al., 2014). A gradual decrease in spring thermal conditions has been observed earlier in a GDD5 spring onset reconstruction from the Schleinsee section in Southern Germany, where the early season and spring cooling during the Bølling culminates in the Older Dryas, after which the conditions improve again towards the Allerød (Wagner-Cremer and Lotter, 2011). Although the Schleinsee site is located much further to the south, the comparable records point towards a dynamic seasonality pattern during this period.
The following warm period (Allerød) is characterised by a temperature increase occurring in all proxies. The summerrelated proxies reveal peaks at ∼ 13 600-13 700 cal yr BP with +12 • C pollen MJJA temperature and +13.7 • C chironomid-inferred July temperature. This initial peak is followed by a small decrease in the pollen-based MJJA and GDD5 spring onset data that is not picked up by the chironomid record. During the Allerød, B. nana disappears from the record while B. pendula and P. sylvestris become common in the macro-remain record. This succession provides additional evidence for stable, mild conditions. During the Younger Dryas, the signals captured by different proxies diverge. At ∼ 12 700 cal yr BP, a rapid temperature decline from ∼ 13 to ∼ 10 • C occurred according to chironomid-based reconstructions, but there is a notable variability between the individual samples ( Fig. 4). At the same time (∼ 12 650 cal yr BP) the pollen-based winter (DJF) temperatures decrease as well, reaching to minimum −20 • C at ∼ 12 000-12 100 cal yr BP. The chironomid-based July temperatures (+8.6 • C) reached minimum (+8.6 • C) around 12 500 cal yr BP. This cooling is not entirely paralleled by pollen-based summer (MJJA) temperature reconstructions, which remains at the previous (Allerød) level, even indicating a small rise between −12 100 and 12 500 cal yr BP, and only decreasing from 12 100 cal yr BP until the onset of the Holocene. Spring season conditions follow the declining winter temperatures rather than the summer signal. The GDD5 values decrease from the ∼ 12 400 cal yr BP peak values of above 1000 to low ∼ 720 values at 12 250 cal yr BP. Unfortunately, we do not have good data coverage for this time interval and thus cannot provide a detailed analysis on the validity of the peak values or the exact timing of the cooling. From 12 200 cal yr BP, the GDD5 values stabilise on very low values, varying between 700 and 800. The results also indicate very late spring onset during this phase, where the growing season did not start earlier than 39-43 d (in June) before MSW. The general pattern during the Younger Dryas in terms of seasonality changes very closely resembles the results obtained for the Hässeldala, where GDD5 values ranging from 650 to 700 were reconstructed, and spring onset occurred ∼ 40 d before MSW (Steinthordottir and Wagner-Cremer, 2019).
The strong variability in several proxies during the Younger Dryas cold reversal (Fig. 4) points to a notably dynamic temperature system. There is an ongoing debate about warm summers on the northern latitudes during the Younger Dryas cooling (Björck et al., 2002;Schenck et al., 2018). The pronounced variability on chironomid-based dataset, duration of the growing season and spring onset data, and rather stable pollen-based MJJA temperatures may indicate short, occasional phases with possible early springs and warmer summer temperatures. However, it does not substantially confirm persistent warm summer conditions during the Younger Dryas, as the main trend in chironomids, pollenbased DJF temperatures, and spring onset and GDD5 records indicate a cooling period. Global climate simulations for the Bølling, Allerød, and Younger Dryas (Schenk et al., 2018) show that for the Baltic region the month of May is highly important, with colder May temperatures during the Younger Dryas than during the Allerød period. This is in broad accordance with the trend in the spring onset record but punctuated by occasional short phases with earlier budburst in the first part of the Younger Dryas. Modern studies suggest that the ongoing lengthening of the vegetation period is also due to the earlier budburst and spring onset (Jeong et al., 2011). The relevance of the early season temperatures requires proxies that can distinguish between the summer or winter temperatures and spring temperatures.

Conclusions
The present study is one of the rare examples of the application of the micro-phenological proxy to the Late Glacial fossil leaf record. By applying a novel micro-phenological proxy and plant macrofossil analysis to the multi-proxy approach, we can constrain early season and spring temperature dynamics rather than traditionally inferred winter or summer temperatures.
The combination of spring onset and GDD5 reconstructions and estimates and plant macrofossil records indicates that the local tree-line dynamics are linked to local rather than regional environmental conditions. The spring onset in central Estonia during the Allerød period was earlier than during the Older Dryas cooling and was comparable to the timing in the modern forest species. However, the Late Glacial plant macrofossil data of northern Estonia displays the lack of any tree remains. Therefore, we could reconstruct a similar situation with modern climate change when the spring onset takes place earlier, but the summer temperatures are not yet rising. The explanation for the eastern Baltic Late Glacial situation could be the proximity of the SIS and the Baltic Ice Lake.
The Late Glacial seasonality dynamics are evaluated by comparing spring onset estimates and pollen-based DJF and MJJA temperature reconstructions in combination with chironomid-based July temperatures. During the GI-1d and the subsequent warming period (Bølling-Allerød) all proxies agree, while the divergence during the YD suggests a dynamic seasonality pattern with possible occasional warm episodes.
Author contributions. LA is the main author and contributed to conceptualisation, investigation, and writing. FWC contributed to conceptualisation, methodology, and writing. JV contributed to chronology. SV contributed to past vegetation and climate reconstructions. LA prepared the manuscript with contributions from all co-authors.
Competing interests. The contact author has declared that none of the authors has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Financial support. This research has been supported by the Eesti Teadusagentuur (grant nos. MOBTP140 and PRG323).
Review statement. This paper was edited by Nathalie Combourieu Nebout and reviewed by two anonymous referees.