Holocene vegetation dynamics of the Eastern Mediterranean region: Old controversies addressed by a new analysis

We reconstruct vegetation changes since 12 ky in the Eastern Mediterranean to examine four features of the regional vegetation history that are controversial: the extent of non‐analogue vegetation assemblages in the transition from the Late Glacial to the early Holocene, the synchroneity of postglacial forest expansion, the geographical extent of temperate deciduous forest during the mid‐Holocene and the timing and trigger for the re‐establishment of drought‐tolerant vegetation during the late Holocene.


| INTRODUC TI ON
The Eastern Mediterranean-Black Sea Caspian Corridor (EMBSeCBIO) region (33°-49° N, 20°-60° E) provides opportunities to examine Holocene vegetation changes and how these have been shaped by climate and human activities, because it is characterised by large climate and elevation gradients and has a long history of human occupation stretching back to the introduction of Neolithic farming in the early Holocene.Several works have summarised vegetation changes for parts of the region or specific time periods (e.g.Bottema et al., 1994;Bozilova & Tonkov, 1995;Connor & Sagona, 2007;Elenga et al., 2000;Magyari et al., 2019;Prentice et al., 1996;Tarasov et al., 1998;Zeist & Bottema, 1991).
However, the limited spatio-temporal coverage of these analyses means that there is no comprehensive overview of the vegetation history of the EMBSeCBIO region and several features of the Holocene record are still a matter of debate, including the extent of non-analogue vegetation in the transition from the Late Glacial to the early Holocene, the timing of postglacial afforestation across the region, the geographical extent of temperate deciduous forest during the mid-Holocene and the exact timing of the expansion of drought-tolerant vegetation during the late Holocene.
The large number of pollen records from individual sites across the region compiled by the EMBSeCBIO project (Cordova et al., 2009;Marinova et al., 2018) allows a more comprehensive analysis of Holocene vegetation changes to be made.
Pollen records from the northern extratropics indicate that some vegetation assemblages during the Late Glacial and early Holocene combined species that are not present together today, so-called non-analogue vegetation types (Magyari et al., 2014(Magyari et al., , 2018(Magyari et al., , 2019;;Williams & Jackson, 2007;Zanon et al., 2018).The existence of non-analogue vegetation types in the EMBSeCBIO region has not been explored in a systematic way.However, it has been argued that some areas of the Carpathians (Feurdean et al., 2007), western Georgia (Connor & Kvavadze, 2009) and the southwestern Black Sea (Bottema et al., 1994) were refugia for tree taxa, including species of Fagus, Ulmus and Acer, during the glacial period.The pollen percentages of these tree taxa are much higher in some early Holocene samples than in samples from vegetation where they occur today, but are consistently associated with high abundance of grasses (Bottema et al., 1994;Feurdean et al., 2007;Robles et al., 2022).
There has been no systematic evaluation of whether these anomalous values exceed observed variability across the full environmental range of the biomes in which these tree species occur today, and therefore whether these samples genuinely represent non-analogue vegetation.
The postglacial timing of forest spread in the EMBSeCBIO region is also controversial.It has been suggested that forest development in the western part of the region began in the early Holocene, while further east the increase in forests occurred 2-5 millennia after the beginning of the Holocene (Djamali et al., 2010;Messager et al., 2017).This delay has been attributed to dry conditions during the early Holocene in the east (Djamali et al., 2010;Joannin et al., 2014;Messager et al., 2013;Wright et al., 2003).However, δ 18 O isotopes from lake carbonates and speleothems and leaf wax δ 13 C isotopes from the region have been interpreted as showing humid conditions since the beginning of the Holocene (Bar-Matthews et al., 2003;Bliedtner et al., 2020;Cheng et al., 2015;Eastwood et al., 2007;Fleitmann et al., 2009;Göktürk et al., 2011;Roberts et al., 2008;Stevens et al., 2001).
Several hypotheses have been proposed to explain the apparent delay in forest development, including the time lag for tree migration from glacial refugia (Connor & Kvavadze, 2009;Messager et al., 2017), the impact of fire (Roberts, 2002;Turner et al., 2010), the fact that increases in precipitation were concentrated outside the growing season (Brayshaw et al., 2011;Dean et al., 2015;Göktürk et al., 2011), a negative feedback from the freshwater lake which occupied the Black Sea basin before it was joined to the Mediterranean (Göktürk et al., 2011;Messager et al., 2017) and potential anthropogenic modification of the landscape (Asouti & Kabukcu, 2014).However, the extent of dry conditions and the inhibition of forest growth in the eastern part of the EMBSeCBIO region are still unclear, making it difficult to evaluate these hypotheses.
Early analyses of vegetation patterns in the mid-Holocene suggested that temperate deciduous forest was more extensive in southern Europe at 6 ky (cal.yr BP) than today, replacing cool mixed and conifer forest in the mountains of southern Europe and Mediterranean vegetation in the lowlands (Huntley & Prentice, 1993;Peyron et al., 1998;Prentice et al., 1996).The replacement of xerophytic vegetation by temperate deciduous forest implies an increase in plant available moisture during the growing season, while winters remained cooler than today.Collins et al. (2012) also showed an expansion of temperate deciduous forest at 6 ky compared to present but argued that the extent of this expansion was overestimated in earlier reconstructions which considered all sites with a mid-Holocene record and did not measure the expansion relative to reconstructions of modern day vegetation at the same site.According to the Collins et al. (2012) reconstruction, drought-tolerant, open vegetation was still present in southern Europe, including in the Eastern Mediterranean region which is the focus of this study, implying the continuation of a pronounced seasonality in the precipitation regime.However, none of these vegetation reconstructions consider possible within-biome variability of temperate deciduous forest and all of them either focus only on the 6 ky window or have limited coverage of the EMBSeCBIO region.Thus, the question of the extent of temperate deciduous forests during the middle Holocene is still unanswered.
It has been suggested that there has been a general increase in drought-tolerant vegetation, reflecting a trend towards more arid conditions, since the mid-Holocene in the Mediterranean region (Jalut et al., 1997(Jalut et al., , 2009;;Magny et al., 2002;Roberts, Brayshaw, et al., 2011;Roberts, Eastwood, et al., 2011;Roberts et al., 2001Roberts et al., , 2004Roberts et al., , 2008)).However, there is debate about whether this occurred simultaneously across the region or indeed was ubiquitous (Herzschuh et al., 2022).Isolating a long-term trend is further complicated by the fact that multi-centennial oscillations in precipitation have been inferred from pollen and speleothem records in the Eastern Mediterranean, but these oscillations are neither simultaneous nor seen everywhere (Bar-Matthews & Ayalon, 2011;Bini et al., 2019;Burstyn et al., 2019;Göktürk et al., 2011;Kaniewski et al., 2018).Furthermore, there is some debate about whether the increase in drought-tolerant vegetation was a consequence of changes in climate or the effect of changes in land use on runoff (Jacobson et al., 2021;Schilman et al., 2001).In this article, we use a recently developed method to reconstruct vegetation from pollen samples which takes account of within-biome variability in taxon abundance and also identifies assemblages that are not typical of any modern biome (Cruz-Silva et al., 2022).This allows us to address whether existing controversies are in part a reflection of the use of subjective approaches to vegetation reconstruction.We address four questions: (1) the extent of non-analogue vegetation during the postglacial period; (2) the synchroneity of forest expansion during the early Holocene; (3) the extent of temperate deciduous forest in the mid-Holocene and (4) the expansion of drought-tolerant vegetation during the late Holocene.

| MATERIAL S AND ME THODS
Vegetation reconstructions for fossil pollen sites in the EMBSeCBIO region were made using the technique described by Cruz-Silva et al. ( 2022) (Appendix S1), using 5765 modern pollen samples from Europe, the Middle East and northern Eurasia derived from the EMSeCBIO and SPECIAL databases (Harrison, 2019;Harrison et al., 2021), which were assigned to biomes based on potential natural vegetation reconstructions (Hengl et al., 2018) (Figure 1) using a previously tested search window of 20 × 20 km [close to the theoretical pollen source area for small (<1 km 2 ) basins accordingly to Prentice (1988) and Sugita (1993)] to characterise the within-biome means and standard deviations of the abundances of individual taxa (see details of the methods in Cruz-Silva et al., 2022).This technique is an improved method of vegetation reconstruction compared to the biomisation approach that was previously applied to the EMBSeCBIO region (Marinova et al., 2018) because it allows the identification of non-analogue vegetation types, i.e. assemblages that consist of associations of taxa with abundances very different from those seen in any modern biome.The new method produces a substantial increase in the prediction accuracy for individual biomes compared to biomisation and a reduction in the number of misassigned samples (Cruz-Silva et al., 2022).We have made a further evaluation of the method to determine how accurate it is in predicting the vegetation represented by core top samples (age <150 cal.yrs.BP) from fossil pollen cores from the EMBSeCBIO region.We discuss the application of this method to fossil records from the EMBSeCBIO region and then describe the analyses of the fossil reconstructions to address the four questions posed earlier.

| Application to fossil records
Fossil pollen data were obtained from the EMBSeCBIO database (Harrison et al., 2021) they come from extinct lakes or stratigraphic sequences), (c) where there is a major hiatus after the youngest radiocarbon date, (d) where more than half of the radiocarbon dates were rejected by the original authors, (e) where more than half of the ages are based on pollen correlation with other radiocarbon-dated records and (f) marine records or cores from very large lakes (>500 km 2 ).We kept records where the original publication indicates a major hiatus but where there are radiocarbon dates above this hiatus allowing an age model to be constructed for the post-hiatus part of the record.As a result of this filtering, vegetation reconstructions were made for 122 continental records from basins ranging in size from <0.01 to 500 km 2 .New age-depth models were produced for these 122 records using the IntCal20 calibration curve (Reimer et al., 2020) and the 'rbacon' R package (Blaauw et al., 2021) in the 'AgeR' R package (Villegas-Diaz et al., 2021).This package provides an optimum model for each record, based on the lowest quantified area between prior and posterior accumulation rate distribution curves (see Harrison et al., 2022).The selected model was checked manually and through visual inspection to ensure that the final age models represented the date information accurately and did not manifest abrupt shifts in accumulation rates or changes at the dated depth.The records have a mean length of 6700 years (Figure 1) and a mean resolution of 265 years (Figure S1).We calculated dissimilarity and similarity scores (Appendix S1) for each fossil sample with respect to each biome in turn.

| Non-analogue vegetation types
We applied the per-biome threshold values determined by Cruz-Silva et al. (2022) to detect whether a given sample exceeded the threshold for assignment to a given biome.Samples that exceeded this threshold for all biomes were considered to represent non-analogue vegetation.Since there is some uncertainty in the matches between fossil and modern samples, we used a 5% threshold to distinguish false positives, and only consider samples that exceed this threshold as actual non-analogues.We quantified the proportion of sites with at least one non-analogue sample in 300-year windows with 50% overlap to produce a time series of the occurrence of non-analogues through time.The size of the window was chosen to approximate the mean resolution of the records (265 years).Since the time series are expressed as a percentage of total records in each window, we examined the number of pollen records available in each window to verify that changes were not an artefact of data availability.We also checked that individual samples were not depauperate and the counts sufficiently large to ensure an adequate representation of the vegetation in the pollen assemblage.Previous research has shown that high-elevation records may be contaminated by upward transport of pollen from lower elevation sites (Takahara et al., 2000), and this might provide an explanation for apparently non-analogue samples.To test this, we created separate time series for high (>1500 m a.s.l.), medium (between 500 and 1500 m a.s.l.) and low (<500 m a.s.l.) elevation sites.

| Timing of the early Holocene forest expansion
We used terrestrial sites spanning at least 3000 years during the period between 7 and 11.6 ky, with at least 13 pollen samples during this interval, and where the oldest pollen sample had an age of at least 11.2 ky to analyse the timing of forest expansion.Temperate deciduous malacophyll broadleaf forest (TEDE) and cool mixed evergreen needleleaf and deciduous broadleaf forest (CMIX) were considered as forest; samples classified as belonging to other biomes were considered together since they represent less moisturedemanding vegetation types.The non-normalised similarity scores for each pollen sample were summed separately for both forest and other vegetation and then expressed as a percentage.The resulting curve was smoothed using the mean value in windows of 300 years with 50% overlap.This compositing approach has been used with other types of data, for example charcoal and speleothem data (Parker et al., 2021;Power et al., 2010), to minimise the effect of age uncertainty and differences in sampling resolution between records and to emphasise the common signal across sites.Break point analysis (Zeileis et al., 2003) was applied on the time series for each individual site to obtain the optimal number (smallest magnitude of residuals) of break points in the trends in the sequence (score ~ age) using the R package 'strucchange'.Linear regressions were performed between the identified break points in each sequence.The earliest point indicating a change from zero or negative slope (no change or decrease in forest) to a positive slope (increase in forest) was taken as the start of forest expansion.The mean and standard deviation between all estimated inception points was calculated to obtain the approximate age of the start of forest expansion at each site.We explored whether there were differences in the timing of initial forest expansion as a function of latitude, longitude and elevation.

| Timing of mid-Holocene expansion of temperate deciduous forest
We examined all continental records that wholly or partially covered the mid-Holocene, defined here as the interval from 7 to 4 ky.We determined the predicted biome in successive non-overlapping 300-year windows centred on 7.0, 6.5, 6.0, 5.5, 5.0, 4.5 and 4.0 ky for each record, where the width of the window allows for potential chronological uncertainties.In cases where more than one biome was predicted during a time window at a given site, we used the most common biome across all samples (following Bigelow et al., 2003) or, in the case of a tie, the biome with the highest sum of similarity scores.In order to express the changes of TEDE through time, we filtered the sites predicted as TEDE and expressed as the percentage of TEDE per time window relative to the total number of sites.To express the changes of TEDE through space as maps, in each time-window, we categorised the sites predicted as TEDE into one of the following categories: 'Remains as TEDE' for records where TEDE had been present in the previous interval and persisted, 'change to TEDE' for records which were not TEDE in the previous interval and became TEDE and 'TEDE' for sites which were TEDE in the selected time window but had no record in the previous time window.

| Evaluation of the expansion of drought-tolerant vegetation in the late Holocene
We examined the vegetation trends from 6 ky onwards and through the late Holocene by creating a composite time series of the number of sites classified as forest, where forest includes samples allocated to TEDE and CMIX, in 300-year windows with 50% overlap.As the intention was to compare the vegetation time series with evidence for changes in both human population density and speleothem oxygen isotope records of hydroclimate, we focused these analyses on the region where there are speleothem isotope records covering at least 3000 years from 6 ky to the present.Additionally, we assessed the robustness of the relationships between changes in forest cover, human population and climate by examining these trends at a sub-regional scale.We focus on the three regions where there is sufficient archaeological data to construct reliable summed probability distribution (SPD) curves: Greece, Anatolia and the Levant.
The SPD of archaeological radiocarbon dates has been widely used as a measure of changes in human population density (see Rick, 1987;Robinson et al., 2019), including in the Eastern Mediterranean (Palmisano, Lawrence, et al., 2021;Weiberg, Bevan, et al., 2019).We extracted data from the Mediterranean basin region from the p3k14c global database of archaeological radiocarbon dates (Bird et al., 2022).An SPD curve covering the interval from 6 to 2.5 ky was produced using 5979 calibrated radiocarbon dates (using the Intcal20 calibration curve) from 620 sites, without normalising and using the 'rcarbon' R package (Bevan et al., 2022).This approach cannot be used for intervals younger than 2.5 ky because dating of archaeological records in the more recent period tends to rely on typo-chronological schemes defined by short-lived pottery types and coins rather than radiocarbon (Palmisano, Bevan, et al., 2021;Roberts et al., 2019;Weiberg, Hughes, et al., 2019).Settlement data from archaeological surveys have been used to reconstruct populations changes in the interval younger than 2.5 ky (e.g.Palmisano et al., 2017;Roberts et al., 2019), but are only available for a few limited sub-regions, for example for the Peloponnese and Macedonia (Weiberg et al., 2016) in Greece, the Sagalassos basin and the Burdur province in south western Anatolia (De Cupere et al., 2017;Dusar et al., 2012) and thus cannot be used to construct a regional composite.Intervals of higher and lower than average human population density according to the regional SPD curve between 6 and 2.5 ka were identified using an exponential null model of population growth, where the periods of growth (and decline) were those intervals in the SPD that were significantly different from the fitted model (Timpson et al., 2014).An exponential null model was used because, theoretically, this accounts for constant homogeneous taphonomic loss of archaeological sites without having to resort to any other correction bias (Timpson et al., 2014).We also compared the regional SPD with SPD curves constructed for sub-regions, Greece, Anatolia and Levant, using the same dataset.
Oxygen isotope records from speleothems are generally interpreted as indicating changes in moisture in the region above the cave site (Bar-Matthews et al., 2003;Burstyn et al., 2019;Cheng et al., 2015;Fleitmann et al., 2009).We extracted five oxygen isotope records in the EMBSeCBIO region from the SISAL v2.2 database (Comas-Bru et al., 2020, b).The SISAL database contains multiple age models for each record; we preferentially used models built with Bchron (which is the most common model type in the database), followed by Bacon or the originally provided ages.

| Prediction accuracy of the modern core-top reconstructions
The quantitative evaluation of the ability of the training set to predict modern core-top pollen samples from the EMBSeCBIO region showed a balanced accuracy of 54% when only the dominant biome within a 21 km search window was considered (Appendix S1; Table S1).As in the evaluation using the modern samples from the region (Cruz-Silva et al., 2022), the accuracy increased substantially when both the dominant biome and the sub-dominant biome were considered, reaching 69% (Table S1).This level of accuracy is not as good as the values obtained using modern pollen samples (with ages <50 years) from the region (64% and 76% respectively), most likely reflecting the wider time interval (c.150 years) covered by the core-top samples.Nevertheless, the technique correctly predicted modern vegetation patterns in the region including grasslands in the southeastern Caucasus, Mesopotamia and the Levant, temperate forest in the Balkans and warm temperate shrublands and forest in the Mediterranean lowlands.As in the case of the modern evaluation, misallocations occurred between closely related biomes (e.g.TEDE and CMIX) and were frequently in the Carpathian region (Figure 2) where such misallocations likely reflect pollen transport up elevational gradients.

| Non-analogue vegetation types
The percentage of records with non-analogue samples in the earliest part of the record, between c. 12.3 and 11.8 ky, does not exceed the 5% threshold for false positives (Figure 3a).There is a rapid increase in the number of records with non-analogue samples after c. 11.8 ky, and the highest values of the entire record occur between c. 11.5 and 9.5 ky (Figure S2).There is a gradual decrease in the percentage of records with non-analogue samples until c. 6.0 ky, after which the values do not exceed the 5% threshold for false positives.There is no obvious relationship between changes in the number of sites with non-analogue samples and the total number of sites per time window, and there are no windows with less than 20 records (Figure 3a), so it is unlikely that the changes in the representation of non-analogues are an artefact.
Most (~40%) of the records with non-analogue samples during the interval characterised by the maximum occurrence of non-analogues are from sites above 1500 m a.s.l.(Figure S3).The records are concentrated in the Carpathian region (Figure 3b), where there are seven sites with non-analogue samples.There is a single record west of Georgia, and one from the mountains south of the Black Sea.Nonanalogue samples are generally characterised by atypical abundances of specific pollen types rather than simply the presence or absence of specific pollen types.The taxonomic composition of samples with no modern analogue, aggregated for the Carpathian region (Figure S4), showed that the most abundant taxon is Ulmus (>20%), but other tree species including Pinus, Alnus, Picea, Betula, Fraxinus and Corylus were present and herbaceous taxa including Cyperaceae, Artemisia, Oxyria/Rumex and Ericaceae were also abundant.The taxonomic composition of samples from the Gagra site, eastern Georgia, showed that the most abundant taxa are Abies (>20%), followed by

| Post-glacial forest expansion
Comparison of the similarity scores for forest versus other vegetation shows a transition from a predominantly open or drought-tolerant Late Glacial vegetation to dominance of moisture-demanding forests during the early Holocene (Figures S2 and S5).However, the nature of this transition varies geographically.At some sites in the Carpathians (e.g.Kismohos, Avrig 1, Mohos1, Sterogoiu), the Rila and Rhodope mountains (e.g.Dry Lake 2, Kupena) and eastern Caucasus (e.g.Didachara), the increase in the forest was large and it became the dominant vegetation type in a matter of centuries.Other sites in the Aegean region (e.g.Gramousti, Orestias G25, Lake Maliq S1) show a more gradual and fluctuating increase in forest, while sites in Zagros (e.g.Lake Zeribar) show only a limited and slower increase.Considering all sites, the mean date of the start of forest increase is 10.64 ± 0.65 ky (Figure 4).The individualistic response of each site is consistent with the fact that there is no discernible difference in the timing of the transition from open vegetation to forest either with latitude, longitude or elevation across the sites (Figure 4).

| Temperate deciduous forest
There is no obvious expansion of TEDE at 6.0 ± 0.15 ky compared to the previous intervals (Figure 5), but rather changes in the geographic distribution of these forests (Figure 6).The distribution of TEDE at 7.0 ± 0.15 ky was restricted to montane areas in the Carpathians, Zagros, Balkans and Caucasus, and the Euxine region around the Black Sea (Figure S2; Figure 6).TEDE had disappeared from the Carpathian region in the next time window but was present at 6.0 ± 0.2 ky, and also appeared at sites in the Euxine region and in the Aegean region.The maximum expansion of TEDE occurs after 6.0 ky (Figure 5;

| Late Holocene forest dynamics
The maximum in forest cover occurred between c. 4.5 and 3.6 ky (Figure 7c; Figure S2), when more than 50% of the sites are reconstructed as forest.There is an abrupt decline in forest cover from 50% to 30%, starting at c. 3.3 and reaching a minimum at 2.6 ky.
There is no obvious coherence between the changes in forest cover and population changes as evidenced by the SPD curve.
The interval of highest forest cover overlaps with a period of high human population density; the decrease in human population density around 3.5 ky occurs before the abrupt decrease in forest sites at 3.3 ky.Although the brief increase in forest cover around 2.8 ky corresponds to an interval of declining population, the subsequent interval of low forest cover is not associated with an increase in population, though this may reflect uncertainties in the SPD.However, the lack of coherency at a regional scale does not preclude the possibility that there were anthropogenic impacts on vegetation at more local scales.
The changes in the number of forest sites does show a close resemblance to intervals of wetter/drier conditions as inferred from the oxygen isotope records (Figure 7).The records from Jeita, Mavri Trypa, Skala Marion and Sofular show wetter conditions when the percentage of forest sites is high, particularly between 4.5 and 3.3.
They also show a sharp increase in aridity near 3 ky (though this occurs somewhat later at c. 2.8 ky in Skala Marion) that persists until at least c. 2.3 ky.All of the speleothem records show somewhat wetter conditions corresponding to the forest expansion between c.

and 1.3 ky, though this is most marked at Mavri Trypa and Skala
Marion.The Sofular record shows increasing aridity during the last millennium, and this would be consistent with the cessation of speleothem formation at Jeita, Mavri Trypa and Skala Marion, but is inconsistent with the relatively high number of forest sites.However, the Kocain record does not show a strong increase in aridity in the recent 1000 years, suggesting that there may be considerable spatial heterogeneity in the nature of late Holocene moisture changes.
Further analyses of radiocarbon data for specific sub-regions (Figure 8) reveal distinct human population density trajectories.
Population levels were low in Greece during the mid-Holocene and maximum population levels were reached around 3 ky.In contrast, the Levant and Anatolia exhibited high population levels during the mid-Holocene, followed by a general population decline.There are limited numbers of pollen records (Figure S8) for some of these sub-regions, with an average of only 7.5 records per time window for Anatolia sites, and an average of only three sites for the Levant, which may affect the reliability of the composite forest cover curves.Nevertheless, the forest curves for Greece and Anatolia show a closer resemblance to changes in moisture, as shown by speleothem isotope records, than to human population changes.In both regions, increases in forest cover correspond to both increases and decreases in population; there is a similar lack of synchroneity with population changes in times of decreasing forest cover.The forest cover for the Levant does not align well with either humidity or population changes, but this likely reflects the very limited number of records from the region.

| DISCUSS ION
The Late Glacial in the EMBSeCBIO region was characterised by open vegetation analogous to modern cold and dry affinity assemblages.However, a large number of early Holocene records represent vegetation types that are not present in the region today.The maximum number of records with non-analogue samples occurred from c. 11.5 to 9.5 ky.Samples from records from c. 6.0 ky onwards were all attributed to existing biome types (Figure 3a).This pattern is consistent with the analysis of vegetation changes across Europe using a modern analogue approach (Zanon et al., 2018), which also found an increase in the number of non-analogue records from the start of the Holocene, reaching a peak between 9.0 and 6.0 ky.This pattern contrasts with records from e.g.North America, where nonanalogue vegetation is more typical of the Late Glacial than the early Holocene (Williams & Jackson, 2007).
Almost all of the sites with non-analogue samples come from the Carpathians (Figure 3b).This region has been identified as a glacial refugium for a number of coniferous and deciduous trees.
High-resolution pollen records and macrofossils have shown the occurrence of Pinus, Juniperus, Betula, Salix and Picea in lowland areas and the Romanian Carpathians prior to 14.7 ky (Feurdean et al., 2007(Feurdean et al., , 2010;;Pató et al., 2020;Tanţău et al., 2014).Plant genetic data also point to this region as a glacial refugium (Magri et al., 2006;Petit et al., 2003), while genetic evidence of canopy forest snails also supports the existence of temperate trees in the Carpathians during the Late Glacial (Juřičková et al., 2019).
Evidence for a very rapid expansion of Larix, Pinus and Ulmus in the early Holocene strongly supports the idea of a glacial refugium for these species in Romania (Feurdean et al., 2007).The extremely high values (>20%) of Ulmus in the non-analogue samples in the early Holocene (Figure S4), when this taxon does not exceed a mean of 0.4% in spectra from forest biomes where it is found today (Cruz-Silva et al., 2022), support the idea of rapid expansion from refugial conditions.Similarly, the lowlands of eastern Georgia on the edge of the Black Sea, where the Gagra site is located, has also been previously identified as a Late Glacial refugium where Abies was the dominant taxon (Connor & Kvavadze, 2009).This is consistent with the fact that Abies was the dominant taxon at this site (Figure S4), though Ulmus is also present in relatively high abundance (3.5%).
Quantitative reconstructions of Last Glacial Maximum climates (Davis et al., 2022) indicate somewhat moderate cooling in regions to the west and south of the Black Sea accompanied by increased growing season precipitation, which could favour the persistence of moisture-demanding broad-leaved trees along with conifers.
The timing of forest expansion after the Younger Dryas in response to increasing temperature and humidity occurred between 10.64 ± 0.65 ky (Figure 4).There is no evidence of systematic differences in timing between the western and eastern parts of the  (Feurdean et al., 2007(Feurdean et al., , 2010;;Finsinger et al., 2017;Huntley & Prentice, 1993;Jamrichová et al., 2017;Prentice, 1985;Tantau et al., 2003;Tanţău et al., 2011Tanţău et al., , 2014)).There is, however, a tendency for sites in the western part of the region and high-elevation sites to show more abrupt increases in forest cover than those in the eastern part of the region.This, coupled with the fact that forest density in the eastern sites never reached the levels observed in the western sites, may underpin the perception that forest expansion was delayed in the eastern part of the EMBSecBIO region.
Temperate deciduous forests were present in Anatolia and southern Greece by 6 ky (Figure 6) and temperate deciduous trees were an important component of the vegetation across much of the region at that time.This is consistent with earlier findings (e.g.Huntley & Prentice, 1993;Peyron et al., 2017;Prentice et al., 1996Prentice et al., , 2000;;Roberts et al., 2004) and with quantitative climate reconstructions showing that most of the Mediterranean region was characterised by lower temperatures and greater plant-available moisture during the mid-Holocene (e.g.Bartlein et al., 2011;Cheddadi et al., 1997;Davis et al., 2003).It is also consistent with speleothem oxygen isotope records from the Middle East, which show higher precipitation than present during the mid-Holocene (Burstyn et al., 2019).
However, our analyses suggest that the expansion of TEDE occurred later than this, since both the 5.5 and 5.0 ky time windows registered a higher abundance of TEDE than at 6 ky (Figure 6).The focus on 6 ky in earlier studies was motivated by the fact that this was an interval chosen for palaeoclimate modelling (Prentice et al., 2000), being a compromise between the maximum changes in northern hemisphere insolation and the relict presence of northern hemisphere ice sheets (Joussaume et al., 1999), rather than the idea that this might represent the maximum expression of climate-driven vegetation changes.The increase in TEDE after 6 ky is consistent with warm and wet conditions inferred from marine and terrestrial biological and biogeochemical proxies in sediment cores from North and SE Aegean and Levant Seas (Triantaphyllou et al., 2014).
Late Holocene forest cover exhibited centennial-scale variations rather than a systematic decline.The abundance of forest sites shows millennial scale variability with peaks between 4.5 and 2.3 ky, at 2.7 ky, and 2.3 and 1 ky.The early interval of high forest cover is coincident with a wetter interval starting at 4.8 ky in the Jeita record and 4.5 ky in the Skala Marion record and lasting until 2.9 ky (Cheng et al., 2015).The δ 13 C record from Sofular also shows increased humidity between c. 4.5 and 3 ky (Fleitmann et al., 2009;Göktürk et al., 2011).The gradual increase in forest cover after 2.3 ky is also consistent with the speleothem evidence of increased humidity.The cessation of speleothem formation at Mavri Trypa, Skala Marion and Jeita at varying times after 2 ky could reflect aridification, but this is not supported by the Sofular and Kocain records which show a continuation of relatively moist but fluctuating conditions, consistent with the moderate fluctuations in forest cover.
Intervals of high forest cover occur during intervals of both high and low population density.At a regional scale, population dynamics and forest cover showed limited coherence, emphasising the importance of considering climate influences on forest changes.Even at sub-regional scale, there is no relationship between population density and forest cover in Greece and Anatolia.In regions where there are sufficient pollen data to construct a reliable composite curve, the fluctuations in late Holocene forest cover seem to follow centennial-scale patterns that are more tied to variations in climate.An improvement in the coverage of vegetation records for other parts of the EMBSeCBIO region is required to test the relationship between human activities and vegetation more thoroughly.It is not possible to use the SPD 'dates as data' approach to reconstruct changes in population density after 2.5 ky in the EMBSeCBIO region and thus we are unable to evaluate the role of humans in the fluctuations of forest cover during the most recent interval.However, the abrupt drop in population at c. 3 ky is robust and has been documented previously (Palmisano, Lawrence, et al., 2021;Weiberg, Bevan, et al., 2019;Weiberg, Hughes, et al., 2019) and is thought to correspond to a period of political collapse and social transformation (Knapp & Manning, 2016).While the short-lived increase in forest cover at 2.9 ky might be associated with this abrupt decline, the subsequent interval of reduced forest cover cannot be explained as a consequence of human activities since population levels remained low during this time.Thus, while there is scope for more thorough investigations of the role of human activities on the landscape at a local scale, the evidence currently available suggests that late Holocene vegetation changes were more likely to have been driven by changes in climate than human activities.
F I G U R E 6 Dynamics of temperate deciduous forest during the mid-Holocene in the Eastern Mediterranean-Black Sea Caspian corridor (EMBSeCBIO) region.The plots show sites where temperate deciduous malacophyll broadleaf forest (TEDE) was inferred as the dominant biome within specific 300-year windows, where the colours show whether these were records where TEDE had been present in the previous interval and persisted (remains as TEDE), records which were not TEDE in the previous interval and became TEDE (change to TEDE) or records which were TEDE in the selected time window but had no record in the previous time window (TEDE).Sites where the predicted biome is not TEDE are shown as open circles.The biome reconstructions for these sites are shown in Figure S2.
We have used a new approach taking account of within-biome variability in pollen assemblages to make vegetation reconstructions (Cruz-Silva et al., 2022), and shown that this provides more reliable reconstructions than the classic biomisation method.Nevertheless, it is not perfect since there are still some mismatches between the predicted and observed vegetation.These mismatches may, in part, reflect the diversity of site types or differences in the size of the basin since both of these factors can influence the pollen source area and the quality of pollen preservation (Davis, 2000;Prentice, 1985).However, the rigorous filtering applied to the records, for example by removing large basins and records that showed obvious indications of poor preservation, should have minimised the impact of site differences on the vegetation reconstructions.Even in the early Holocene there are at least 20 records from the region after this filtering.Our focus on the regional picture of vegetation changes to address specific controversies should also minimise the impact of potential uncertainties in the reconstructions.
We have largely focused on vegetation changes that reflect changing plant-available moisture and the seasonality of precipitation, and have used speleothem data as an independent source of information on moisture changes through the Holocene.Although quantitative reconstructions of climate variables have been made at some sites from the EMBSeCBIO region (e.g.Davis et al., 2003;Herzschuh et al., 2022;Robles et al., 2022), they either use a PFTbased modern-analogue approach or local calibrations, which may provide reasonable reconstruction statistics but tends to underestimate climate variability.Furthermore, these reconstructions focus on mean annual precipitation rather than plant-available moisture.
It would be useful to develop robust reconstructions of bioclimatic variables for the EMBSeCBIO region in order to be able to determine how far the observed vegetation changes reflect changes in climate factors influencing plant growth.F I G U R E 8 Comparison of reconstructed changes in vegetation in over the past 6000 years with evidence for changes in human population and palaeoclimate at sub-regional scale (Anatolia, Greece and Levant).The small coloured dots on the map indicate sites of archaeological radiocarbon dates (Bird et al., 2022), the coloured dots with black border on the map indicates sites of pollen records and the triangles on the map show location of oxygen isotope records (δ 18 O) from speleothems. of forest increase between sites suggests multiple factors, including relative location of refugia, plant dispersal abilities, competition and disturbance regimes, could have influenced the response to early Holocene climate change.Temperate deciduous forest was not more extensive than today at 6 ky; the maximum mid-Holocene expansion occurred at c. 5.5 and 5.0 ky.Late Holocene forest cover exhibited centennial-scale variations rather than a systematic decline, at both regional and sub-regional scale.Fluctuations in forest cover during this interval are broadly consistent with speleothem records of changing moisture availability than aggregate regional changes in human population.However, more well-dated pollen cores and archaeological evidence on human population density are necessary to establish definitive conclusions about the interaction of climate, vegetation and people at sub-regional scales.

ACK N O WLE D G E M ENTS
, which contains records derived from publicaccess databases and the original authors.The database contains 187 Holocene records.Vegetation reconstructions were made after filtering to remove records (a) with no radiocarbon dating, (b) where the age of the uppermost pollen sample was unknown (e.g. because F I G U R E 1 Distribution of modern and fossil pollen sites (WGS84 Mercator projection).(a) Distribution of modern pollen entities from the SPECIAL modern pollen dataset (SMPDS) and the Eastern Mediterranean-Black Sea Caspian corridor (EMBSeCBIO) database.The colours represent the biome at each site derived from the Hengl et al. (2018) reconstruction of potential natural vegetation.The biome codes are: CENF, cold evergreen needleleaf forest; CMIX, cool mixed evergreen needleleaf and deciduous broadleaf forest; DESE, desert; ENWD, evergreen needleleaf woodland; GRAM, graminoids with forbs; TEDE, temperate deciduous malacophyll broadleaf forest; TUND, tundra; WTFS, warm-temperate evergreen needleleaf and sclerophyll broadleaf forest; XSHB, xeric shrubland.(b) Distribution of the fossil pollen records from the EMBSeCBIO region that completely or partially span the last 12 ky.The colour represents the time length of the record.
Alnus, Carpinus, Pinus, Fagus, Ulmus, and non-arboreal taxa including Amaranthaceae, Artemisia and Asteraceae.In general, non-analogue samples from the Carpathians and eastern Georgia are characterised by assemblages dominated by moisture-demanding trees, but taxa indicative of open landscapes were present in high abundance.On the other hand, non-analogue samples from the Yenicaga site, southwest of the Black Sea, have a different taxonomic composition, being dominated by Equisetum (>20%), followed by Poaceae, Cyperaceae, Pinus, Abies and deciduous Quercus.Although the high abundances of tree taxa suggest that the assemblages represent in situ vegetation, nevertheless the co-occurrence of moisture-demanding trees with high abundance of non-arboreal taxa raises questions about whether these samples could represent contributions from nearby wetland or alluvial settings.

F
Comparison of (a) observed and (b) reconstructed biomes from modern samples in the Eastern Mediterranean-Black Sea Caspian corridor (EMBSeCBIO) region (WGS84 Mercator projection).The observed biomes are derived from the Hengl et al. (2018) reconstruction of potential natural vegetation.The biome codes are: CMIX, cool mixed evergreen needleleaf and deciduous broadleaf forest; DESE, desert; ENWD, evergreen needleleaf woodland; GRAM, graminoids with forbs; TEDE, temperate deciduous malacophyll broadleaf forest; WTFS, warmtemperate evergreen needleleaf and sclerophyll broadleaf forest.F I G U R E 3 (a) Proportion of records from the Eastern Mediterranean-Black Sea Caspian corridor (EMBSeCBIO) region having at least one sample identified as non-analogue in a 300-year time-window over the past 12,000 years, where the time windows were constructed with 50% overlap.The red line indicates the 5% threshold to separate false positives (values below the threshold).The blue highlighted area corresponds to the period of maximum occurrence of non-analogues.The grey histogram shows the number of records through time in 300-year time window with 50% overlap.(b) The location of sites with more than 5% non-analogues samples during the interval between 11,500 to 9500 yr BP (red circles); sites with less than 5% non-analogues samples are shown as grey circles (WGS84 Mercator projection).The size of the circles represents the elevation of the site.
Figure S6), with additional sites appearing in the Balkans and in Greece in the 5.5 ky time window and several additional sites showing a shift to TEDE in the 5.0 ky time window.Subsequently, TEDE disappeared from lowland sites around the Aegean and on the Greek peninsula, though TEDE persisted or expanded in the Balkans and the northwestern Caucasus.
2.3 ky to reach levels of around 45% between c. 2.3 and 1 ky.Although there are fluctuations in the number of forest sites after 1 ky, the number remains relatively high until the present.Thus, there is no explicit evidence of a systematic decrease in forest over the late Holocene; the pattern F I G U R E 4 The timing of initial forest expansion during the early Holocene in the Eastern Mediterranean-Black Sea Caspian corridor (EMBSeCBIO) region.The sites are organised by latitude (from north to south) and colour coded to show longitude.The size of the circles shows the elevation of the sites.The black line at 10.64 ky indicates the mean value of forest inception, while the grey bar shows the standard deviation of ±0.65 ky.F I G U R E 5 Proportion of records reconstructed as temperate deciduous malacophyll broadleaf forest (TEDE) in a 300-year window centred on a specific time slice in the mid-Holocene, here defined as 7 to 4 ky.rather points to centennial-scale variations in forest abundance and corresponding increases in more open vegetation, as represented by samples reconstructed as tundra (TUND), desert (DESE), graminoids et al.EMBSeCBIO region, or between high and low elevation sites.For instance, at 10.5 ± 0.15 ky (Figure S2), sites located at both easternmost and westernmost longitudes, spanning from the Caucasus to the Balkans, underwent a transition from predominantly open biomes, primarily GRAM, to forested biomes including CMIX and TEDE.This spatial distribution aligns with the broader shift from open or drought-tolerant Late Glacial vegetation to moisture-demanding forests during the early Holocene.Differences in the speed of forest expansion, and hence the timing of peak forest cover, indicate a somewhat individualistic response at different sites.This is consistent with previous studies which have shown that vegetation changes were of variable amplitude and timing and have explained this in terms of interactions of external and internal drivers including glacial refugia, dispersal abilities, migration patterns, climate, soil development, competition and disturbance regimes Vegetation with no modern analogue occurred from c. 11.8 to c. 6 ky, mostly in the Carpathians, with a maximum number of records between 11.5 and 9.5 ky.The composition of non-analogue samples and the location of the records suggest that these vegetation types arose because of rapid expansion of individual species from Late Glacial refugia.There is no discernible difference in the timing of forest expansion in the early Holocene either with latitude, longitude or elevation across the sites.Differences in the exact timing and speed F I G U R E 7 Comparison of reconstructed changes in vegetation in over the past 6000 years with evidence for changes in human population and palaeoclimate in the Eastern Mediterranean-Black Sea Caspian corridor (EMBSeCBIO) region.The map (a) shows sites used in the reconstruction of biome changes.Changes in human population density (b) are shown by the summed probability distribution (SPD) of radiocarbon dates on archaeological material.The biome plot (c) shows the proportion of records in the region characterised by forest in a 300-year time-window through the past 6000 years, where the time windows were constructed with 50% overlap.Changes in palaeoclimate are inferred from the oxygen isotope records (δ 18 O) from individual speleothem records (d-h), where more negative values are taken to indicate wetter conditions.
ECS and SPH acknowledge support from the European ResearchCouncil-funded project GC 2.0 (Global Change 2.0: Unlocking the past for a clearer future; grant number 694481).ICP acknowledges support from the ERC under the European Research Council under the European Union Horizon 2020 research and innovation programme (grant agreement no: 787203 REALM).No permits were needed to carry out this work.CO N FLI C T O F I NTE R E S T S TATE M E NTNone of the authors have a conflict of interest to disclose.