Introduction

Sandstone regions represent unique landscapes within central Europe and are characterized by steep geomorphological gradients with deep valleys, rocks and flat plateaux. Sandstones set up fine-grained mosaics with a high diversity of habitats (Härtel et al. 2007). Patches with contrasting disturbance regimes, microclimates and light conditions enable the coexistence of species with opposite ecological demands (Gutzerová and Herben 2001). They also provide conditions for the deposition of various palaeoecological archives, through which we can study this species co-existence at large time scales. The fine-grained mosaics of conditions have played a key role in the vegetation history of central Europe, because they have allowed interglacial flora to survive during glacial periods (Kuneš and Jankovská 2000; Svoboda et al. 2018) as well as glacial tree species during the Holocene (Pokorný et al. 2017; Svoboda et al. 2018); they also became refugia for human societies for the whole of prehistory (Ptáková et al. 2020) even though they were never intensively cultivated and therefore remained permanently forested.

The geomorphological relief strongly determines the ecological conditions and mosaics of present forest habitats. This ecological factor also represents an advantage for the interpretation of a pollen record that has usually originated from the humid bottom of one of the gorges. Pollen studies from about a decade ago have already brought the first kind of spatially explicit interpretation, even for the rocks and plateaux, thanks to the steep ecological gradient and known ecological preferences of the species (Kuneš and Jankovská 2000; Pokorný and Kuneš 2005; Kuneš et al. 2007). Wet valleys hosted spruce during the middle Holocene and fir and beech during the Late Holocene, while plateaux and rock edges remained during these periods covered by oak and pine, respectively. Those pollen interpretations have recently been matched by the charcoal record from the soil (Novák et al. 2019). The advantage of having both charcoal proxies together, the charcoals from the soil, and an archaeological context, is the possibility to work on a finer spatial scale then ever before (Nelle et al. 2010; Novák et al. 2018). During the last three decades, palaeoecological exploration in Czech sandstone areas has produced tens of records for multiple proxies (Novák et al. 2012, 2019; Svoboda et al. 2018). Here, we build on these data, and use them in a spatially explicit model, providing vegetation estimates independent of the ecological demands of the species. Moreover, we have revealed the ecological demands of tree taxa using palaeoecological reconstruction.

The development of tools for spatially explicit reconstructions from pollen records has advanced significantly and without exception reveal vegetation estimates using pollen productivity estimates and assumptions about pollen dispersal. We applied the Extended Downscaling Approach—EDA (Theuerkauf and Couwenberg 2017), which represents the latest reconstruction algorithm (see comparison with earlier models below).

Palaeoecological research is usually based on the combination of several methods in a multi-proxy approach, where proxies complement each other, for example, different sides of a wetness gradient. The pollen record mainly originates from wet peatbogs and also from a source area around the site. All pollen transport between the source plant and the core is approximated by a Langrangian dispersal function (Kuparinen et al. 2007). Charcoals separated from the soil or soil sediments originate from at least one past fire event and are basically omnipresent (Nelle et al. 2013); thus pedoanthracology can also complement the pollen from drier parts of the area (Touflan et al. 2010; Robin et al. 2013). The charcoal catchment area is influenced by the topography and in general provides information on the vegetation over a very local spatial scale (5 m—Ohlson and Tryterud 2000) much smaller than pollen. Their comparison is possible by making assumptions about pollen and charcoal productivity, and about the dispersal-taphonomic vectors which bring the plant remains from a specific source area. Considering those assumptions explicitly brings interpretation of both proxies together and it means a change from a complementation of the proxies to a validation of the proxies.

The weakness inherent with soil charcoals is the discontinuous character of the record; however, charcoals from archaeological contexts can balance this. The processes involved between wood collection and charcoals in the sediment are still inevitably complex and affected by human behaviour (Théry-Parisot et al. 2010). However, generally, the collection of common firewood is not strongly selective. This assumption is based on the availability of the wood within a short distance from the settlement (Principle of least effort; Shackleton and Prins 1992).

Aims:

  • Provide EDA reconstructions based on pollen and charcoals from soil and archaeological contexts.

  • Compare scenarios according to the input data, and validate and evaluate the best scenario.

  • Link the past composition from reconstruction maps with the topography and obtain long-term autecological preferences of the trees concerned.

Methods

Study area

Sandstone landscapes with their characteristic castellated form are found in the Bohemian Cretaceous Basin, in the north of the Czech Republic. We include four sub-regions: Bohemian Switzerland, the Doksy region, Bohemian Paradise and the Broumov region. For the purposes of visualization, we have placed within each sub-region one example plot of 2.5 × 2.5 km (Fig. 1). In all these sub-regions, basaltic hills are a reminder of tertiary volcanism. Soils in the sandstone areas are mostly acidic; however, the presence of sandstones with higher carbonate contents implies less acidic soils in the Doksy region and Bohemian Paradise. Moreover, all subregions were covered by loess during the cold stages of the Quaternary (Valečka and Šebesta 2013). Shallow lakes that formed during the Full and Late Glacial became terrestrialised during the Early Holocene, with the accumulation of peat also on elevated plateaux.

Fig. 1
figure 1

Map of study area and the four sub-regions: a Bohemian Switzerland, b Doksy region, c Bohemian Paradise, and d Broumov region. The four example plots 2.5 × 2.5 km (red) are shown in Fig. 2

Full size image

Pollen analysis and target taxa

We used 11 pollen sequences (Supplementary Information 1 Table 1) from the Czech pollen quaternary database (PALYCZ; Kuneš et al. 2009). Analysis was focused on seven dominant taxa for the sandstone landscapes: Pinus, Picea, Quercus, Corylus, Betula, Fagus and Abies. The pollen taxa include all species from the genus except for Pinus cembra-type and Betula nana-type, which were excluded from the dataset. The pollen dataset consisted of 735 original samples, which were aggregated in time windows and produced 173 samples with minimal pollen sum 100 grains.

Charcoal analysis

Charcoal analysis was performed in soil profiles and archaeological sites in rock shelters (Supplementary Information 1 Table 1). At all sites charcoal fragments were extracted using a water flotation technique followed by wet sieving with a 250-μm mesh size (Carcaillet and Thinon 1996). Charcoal fragments were manually selected from the sieved fractions, and charcoal analysis was performed only on fragments of the largest fraction (> 1 mm). Charcoals were identified using a standard identification key (Schweingruber 1990) and charcoal reference collection material. The taxonomical identification was performed under a stereomicroscope and a reflection light microscope with 200–500 × magnification. Determination of charcoal fragments was generally feasible to the genus level.

Charcoal analyses in rock shelters

The anthracological research on 20 sites was carried out in the context of parallel archaeological investigations. Charcoal samples from rock shelter profiles were collected from sections cut into 5-cm-thick layers down to the bedrock.

Charcoal analysis in soil profiles

Our study is focused on 27 soil profiles which were dug and sampled for charcoal analyses. Profiles were excavated in three basic geomorphological positions: valley, plateau, and rocky edges. All profiles were taken in flat terrain with 0–5° inclination.

Soil profiles were dug down to the bedrock and each sample contained 10 L of fine-grained soil. Soil samples were taken from 10-cm-thick layers. In samples with a low number of fragments, all fragments were analysed. In samples with abundant charcoal fragments, only a selection of the 250 fragments were analysed.

Chronological control

Age-depth models of 58 pollen and charcoal records were constructed using a total of 199 radiocarbon dates. For three of the pollen profiles, additionally 210Pb dating was available. The age-depth models of 13 charcoal sequences from archaeological sites were also constrained by the dated pottery found. The surface (zero depth) was considered to be the present in most of the models. Age-depth models were constructed by Bacon (Blaauw and Christen 2011) or classical modelling (Blaauw 2010), depending on the number of datings (for details, see Supplementary Information 1 Table 1). Palaeobotanical counts of all proxies were binned into nine time-windows of 1000 years, spanning from 400 to 9400 cal BP.

GIS analysis

For the classification of the sandstone landscape, we used spatial K means clustering into eight spatially similar classes using the Hill-Climbing method (Rubin 1967). For clustering, we selected nine topographic indices selected to capture the landscape morphology, landscape level water movement and solar energy balance (Table 1). These indices were derived from the 5th generation of the digital terrain model of the Czech Republic, provided by the State Administration of Land Surveying and Cadastre, Land Survey Office, Czech Republic, https://www.cuzk.cz. To reduce fine-scale variation irrelevant for our landscape classification, we scaled up from the original 2-m-cell size to the 10-m-cell size using b-spline interpolation (Lee et al. 1997). All topographic analyses were carried out in SAGA GIS 6.4.0 (Conrad et al. 2015). Prior to hydrological analyses we removed sinks using breaching (Lindsay 2014). Ground water level for calculating the vertical distance to the channel network was interpolated from the national (water) channel network DIBAVOD (https://www.dibavod.cz/17/geodatabaze-dibavod.html).

Table 1 List of indices used for the classification of the sandstone landscape
Full size table

EDA and vegetation modelling

We applied the Extended Downscaling Approach – EDA (Theuerkauf and Couwenberg 2017), because it deals with many sites of a similar size, e.g. forest hollows in our case, whereas the first alternative Landscape Reconstruction Algorithm – LRA (Sugita 2007a, b) requires the combination of small sites and a large lake, which is unavailable for this area. EDA estimates are based on pollen records and data on the spatial pattern of landscape classes around sites, whereas LRA estimates are revealed only from pollen records and their variation at all sites. The vegetation pattern in LRA is assumed to be patchy – an even vegetation mosaic within which the sites are placed randomly (Sugita 1994). This can hardly be fulfilled in a real landscape; peatbogs are usually surrounded with specific species, for example, those tolerating a high water level in the soil. The EDA reveals just one scenario, the most optimal one for landscape units, and thus EDA was also preferred against the second alternative – Multiple Scenario Approach (MSA; Bunting and Middleton 2009), which reveals many scenarios.

We used the Extended Downscaling Approach (EDA; Theuerkauf and Couwenberg 2017) to generate likely vegetation composition within each landscape class. The EDA is a forward modelling approach, which assumes vegetation composition within each landscape class. Based on this assumption it searches for the most likely vegetation composition within each landscape class. To this end, the EDA uses an iterative approach. In each iteration, first the landscape pattern around each pollen site is translated into simulated vegetation using the given vegetation composition within each landscape class. Secondly, pollen deposition at each pollen site is calculated using the vegetation pattern, pollen productivity estimates and a pollen dispersal model. Thirdly, the distance between modelled pollen and empiric pollen composition at each pollen site and for each time slice is calculated as the squared-chord-distance.

Iterations start with a random vegetation composition within each landscape class. Using the DEoptim optimizer (Mullen et al. 2011), vegetation composition is than adjusted within each iteration until the distance between simulated and empiric pollen composition is minimal, i.e. until simulated pollen deposition is most similar to the empiric pollen composition for each site and time slice.

In the present study, we included for each pollen record the environmental pattern within a 4 km distance from each pollen site. To this end, we extracted the cover of the eight landscape classes in consecutive rings with a radius of 5, 10, 15, 20, 40, 60, 80, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000 and finally 4000 m. Productivity estimates (Supplementary Information 1 Table 2) were taken from Andersen (1970). Fall speed of pollen was adopted from Soepboer et al. (2007). The dispersal model was set to the Lagrangian stochastic model from Kuparinen et al. (2007), and adjusted to unstable atmospheric conditions (‘lsm unstable’). Distance-weighting considering pollen dispersal calculates with all consecutive rings without the need for explicit setting of a relevant source area of pollen.

When applying EDA to the charcoal records, it had to be slightly adjusted. We now assumed that charcoal composition reflected vegetation composition and hence compared simulated vegetation composition and charcoal composition. We assumed that in the case of the pedoanthracological (soil) charcoal samples, charcoal composition reflected vegetation composition within a 5 m distance from the sample site (Ohlson and Tryterud 2000), and in the case of the archaeological sites (rock shelters) we set the source area charcoal samples to a distance of 200 m, considering the daily need for firewood (Théry-Parisot et al. 2010).

Moreover, we added a charcoalization-factor, which is a parameter analogous to pollen productivity estimates. For a given species, it expresses the relative ratio between the quantity of anthracomass in the charcoal record and the forest biomass of the corresponding species. For the present study, this factor was set to 1 for all taxa, since no empirical values are yet available. We did not apply any distance weighting.

We prepared the following three scenarios of EDA reconstructions according to input data. The scenarios are: P-EDA – pollen, S-EDA – soil charcoals, A-EDA – archaeological charcoals.

Statistical testing

Each EDA scenario, with an input record covering all eight classes and all nine time-windows, has compositional data for 72 spatiotemporal cells. We compared the three scenarios in pairs (A-P, S-P, A-S) by two ways. Quantitative comparison meant the calculation of dissimilarity coefficients using SQchord distance (Simpson 2007) between vegetation compositions of each pair of cells from two compared scenarios.

Qualitative testing consisted in the comparison of presences and absences of taxa in the P, S and A scenarios. We converted EDA percentages to presences and absences for individual taxa, thus we obtained a matrix with 72 spatiotemporal cells for each taxon. Representation between 4 and 100% of EDA composition was considered as a presence, less than 4% as an absence. Comparison of two present-absent matrices led to matching cells (taxa is present or absent in both matrices) and mismatching cells (taxa is present only in one scenario). We randomized both presence-absence matrices and recorded the number of mismatching cells. We repeated the randomization 1000 times and ordered the result by the number of mismatching cells. The number of mismatching cells from real data was compared against the number of mismatching cells produced from randomization. Our p-value is this order of random-based mismatches, which corresponds to real mismatches, and expresses a probability of random pattern, when comparing pattern of presences in both spatiotemporal matrices.

Based on the result of these comparisons we selected one scenario with all time-windows for which we calculated a multivariate analysis, in order to get the relationship between species and nine topographic factors (Table 1), age, fire activity, archaeological occupancy and climate (see their overall trends in Supplementary Information 1 Fig. 1). Fire activity was calculated by Bobek et al. (2019) for the same sandstone area as the one considered in the present paper. It was based on eight microcharcoal sequences and, moreover, three of them were extracted from the same cores as the pollen in this study. In order to obtain a fire trend within a millennial time window, the original z-score with a resolution of 100 years was regressed against reversed age. The positive slope, which we obtained from the regression, shows increasing fire activity within the time window and a negative slope decreasing. Archaeological occupancy is a model of human activity based on archaeological data with a resolution of grid cell 33.3 km2 and a time window of 500 years (Kolář et al. 2022). We characterize the climate by the mean annual temperature and annual precipitation extracted from a global climate model with a grid cell 1 km2 and a time window of 100 years (CHELSA; Karger et al. 2021). Firstly, for both grids we calculated a mean value from grid cells overlapping with our sites—42 grid cells for climate and 17 grid cells for archaeological occupancy. Secondly, we averaged the original temporal resolution to our millennial time windows. Using the compositional data in the best EDA scenario, we calculated Nonmetric Multidimensional Scaling and subsequently projected the fourteen environmental variables (Oksanen et al. 2020).

Results

Sandstone relief and its sampling by palaeoecological methods

The landscape classification reveals three rarer (< 10% each of the study area) and five more abundant landscape classes (Supplementary Information 1 Fig. 2a). The first group, classes 1–3, has classes with higher ruggedness and are mainly found in sandstone landscapes with well-developed sandstone phenomena. These extreme habitats are: 1—southern slopes in valleys with the highest diurnal anisotropic index, 2—northern slopes and bottoms of valleys with the lowest topographic and lowest mass balance indices; and 3 – rocks and edges with the highest topographic and vector ruggedness indices. Classes 1 to 3 are quite rare in the Doksy region (Fig. 2b). The second group, classes 4 to 8, represent smoother parts of the landscape. Considering their topographic position index, class 4 has zero representing flat topography, the rest are concave (classes 6 and 8) with negative values, or convex (classes 5 and 7) with positive values at low (classes 5 and 6) or high (classes 7 and 8) elevations (Supplementary Information 1 Table 3).

Fig. 2
figure 2

Maps of landscape classes in example plots in the four sub-regions a Bohemian Switzerland, b Doksy region, c Bohemian Paradise, and d Broumov region

Full size image

The pollen and archaeological charcoal records cover all the 72 spatial–temporal cells (Supplementary Information 1 Fig. 2b). However, though soil charcoal records are well present in class 7, in the youngest time-window, they are rare in classes 4, 6, 8; for some time-windows they are also present in classes 5 and 7. In total, soil charcoal records cover 54 of the 72 spatial–temporal cells.

Comparison of EDA reconstructions

The three EDA reconstructions show prominent changes through time in all landscape classes (Fig. 3). Class 1 constantly showed a low dissimilarity between P-EDA and S-EDA throughout the whole study period (Fig. 4). Both these strikingly similar reconstructions showed a coniferous succession of Pinus-Picea-Abies. Class 2 was occupied by Pinus during the period 8900–6900 cal BP according to S-EDA and P-EDA, and by Picea during 5900–4900 cal BP according to all three scenarios. For the last three-time windows, S-EDA lacked Fagus and, together with A-EDA, rather indicated Abies and Picea. Correspondingly to P-EDA, Fagus appeared in A-EDA in the last time window. The lowest dissimilarity in classes 3 and 4 was found between P-EDA and A-EDA due to the appearance of Corylus and Fagus, which were lacking in S-EDA. Class 5 was co-dominated by Pinus and Quercus according to all three scenarios, while S-EDA and P-EDA indicated Picea in the period 6900–3900 cal BP. The lowest dissimilarity in classes 6 and 8 is hard to estimate, because of the sparse S-EDA record. Class 7 showed very low dissimilarity between S-EDA and A-EDA due to a dominance of Pinus and Quercus; both taxa appeared in P-EDA together with Abies, Corylus, Betula and Fagus. In total, the comparison between P-EDA and A-EDA had the lowest dissimilarity in 24 spatiotemporal cells, A-EDA and S-EDA lowest in 16, and S-EDA and P-EDA in 14 cells. The dissimilarity coefficient suggests that A-EDA is the most similar scenario to P-EDA. This was even visible in the qualitative comparison of presences and absences of species in individual classes, as, for example, class 3. It contained Betula, Corylus and Fagus in A-EDA and P-EDA (darkgreen in Fig. 5a), whereas S-EDA lacked them (orange in Fig. 5b), and the dissimilarity of P-A is lower than that of P-S.

Fig. 3
figure 3

Mean EDA estimates based on a P—pollen, b S—soil charcoals and c A—archaeological charcoals. For error estimates, see Supplementary Information 1 Figs. 3, 4 and 5

Full size image
Fig. 4
figure 4

Dissimilarity between EDA reconstructions in landscape classes (lines). Dots show the lowest dissimilarity from the three comparisons

Full size image
Fig. 5
figure 5

Comparison between presences ( +) and absences (−) of taxa in EDA reconstructions showing matching (dark and light green) and non-matching (orange and red) cells. Proportions higher than 4% were considered as a presence in all landscape classes (x axis)

Full size image

However, in general, the qualitative comparison indicated a different pattern between P-EDA and A-EDA. The latter, A-EDA, predicted Fagus, Abies, Betula and Corylus to more spatial–temporal cells than S-EDA (dark-green or red in Fig. 5a or orange in Fig. 5c). The taxa Abies and Fagus distinctly appeared in older time-windows for A-EDA (red in Fig. 5a), whereas P-EDA and S-EDA predicted these taxa in the younger time-windows (see Abies in class 1; darkgreen in Fig. 5b and orange in Fig. 5a).

Statistical testing of spatial temporal pattern in A-EDA and P-EDA showed only one significant taxa (Corylus). The comparing of A-EDA and S-EDA provided one significance-edging taxa (Betula), whereas comparing S-EDA and P-EDA resulted in three significant (Abies, Fagus and Picea) and one significance-edging (Quercus) taxa (Fig. 6). We thus selected the pollen and soil charcoal record as the best input data for EDA reconstructions, but for our subsequent analysis (see next section) we considered only P-EDA due to the discontinuous character of the soil charcoal record (for further details, see Discussion).

Fig. 6
figure 6

Ranges of p-values that express a probability that a comparison of two spatiotemporal matrices (Fig. 5) corresponds to a random pattern. Ranges below alpha level (0.05, vertical dashed line) show scenarios with significantly similar patterns

Full size image

The P-EDA results started with a dominance of Pinus (classes 1, 5 and 6) and Corylus (classes 3 and 7) in the oldest time-window (8900 cal BP). For the period 7900–3900 cal BP, results showed a dominance of Picea in classes 1, 2, 5, 6, 8, which starts earlier in class 1 and 8 ([8900-]7900 cal BP) than in class 2 and 6 (~ 5900 cal BP). The youngest time-windows, 3900–900 cal BP, are co-dominated by Fagus and Abies (classes 2, 7 and 8) or Fagus and Betula (classes 3 and 4). Classes 1 and 6 are composed of Abies, Pinus and Picea, and class 5 of Quercus, Pinus and Abies (Fig. 3a).

Relationship of trees to topographic factors

NMDS (nonmetric multidimensional scaling) of the P-EDA ordinated Abies and Fagus on the left, Quercus and Pinus in the centre, Corylus on the right, Picea at the top and Betula at the bottom. Age decreased from right (Corylus) to left (Fagus) and elevation decreased from left and right towards the centre, both variables making the first main gradient in the data. Higher abundance of Fagus correlates with higher mean annual temperature, annual precipitation, archaeological occupancy, but lower age and decreasing trend of fire activity. Picea, Quercus, Pinus and Betula tended to occupy classes with a lower elevation than Corylus, Abies and Fagus. Topographic factors were correlated with each other, following two main patterns. The first pattern was similar to the topographic position index (TPI250) or convergence index (CI). Both indices decreased from the right (high at Corylus) to the top or top left (low at Picea and Abies). The second pattern was similar to the vector ruggedness measure (VRM). This index was high on the right and left (Corylus and Fagus) and decreased towards the top centre or bottom centre (Picea and Betula). Figure 7 and Supplementary Information 1 Fig. 6 show the order of all taxa according to any of the topographic factors. For example, based on vertical distance to channel network (VCN), the closest taxon to the groundwater was Picea; intermediate distances were held by Abies, Betula, Pinus and Quercus. Fagus grew far from the channel network and Corylus even further.

Fig. 7
figure 7

NMDS of pollen-based scenario displaying samples (points), species (blue) and environmental variables (red). For the rest of the environmental variables, see Supplementary Information 1 Fig. 6

Full size image

Discussion

Evaluation of the best scenario and EDA model

The A-EDA results provide a lower dissimilarity to P-EDA than to S-EDA, but the comparison of presences and absences shows only one significant taxon (Corylus). The surroundings of human settlements were affected by the collection of firewood, forest grazing, fire, direct management or their abandonment; for example, stands of Corylus were managed by Mesolithic populations (Ptáková et al. 2020). People thus created site-specific types of vegetation, including early-succession stages, which were an ecological opportunity for spreading taxa. Abies and Fagus appear in the earlier time-windows of A-EDA, and it then becomes detectable by other methods. This phenomenon has already been documented by the finding of Late-Holocene trees (Carpinus, Abies and Fagus) in the Neolithic charcoal record (Novák et al. 2021). This would indicate that we can compensate for the lack of broadleaved taxa in soil charcoals by the use of archaeological charcoals. Additionally, the sediments under rock shelters are preserved against erosion, and bioturbation, and they also have the possibility of archaeological dating. Archaeological charcoals provide valuable information of a qualitative character, but soil charcoals are here shown to be more suitable for the present quantitative approach.

The S-EDA results systematically underestimate the abundance of Betula, Corylus and Fagus, but, in the case of Fagus, it provides a significantly similar pattern to P-EDA. The low estimates are produced by the lack of empirical values for charcoalization factors, which could potentially correct them. Our results show that our ad-hoc setting of putting the charcoalization factor to 1 for all species is not accurate, because the characteristics of wildfires also control the composition of the charcoal assemblage. For example, low anthracomass indicates that mainly dead wood burnt during the wildfires (Novák et al. 2019) and that such wildfires are not able to ignite the large stems of living trees, especially Fagus, known for the low inflammability of its stands (Bobek et al. 2019). The lower representation of Fagus in our soil charcoal record would support its low (< 1) factor of charcoalization. Corylus is also under-represented in the charcoal record, possibly due to its thin and soft wood. At the opposite end of this gradient would lie fire-prone species, e.g. Pinus or Betula (Adámek et al. 2018), which might have factors of charcoalization higher than 1.

The EDA reconstruction tool was designed for the pollen record, thus it corrects most biases of pollen analysis, except for the regional pollen component originating from larger distances (Sugita 1994). Our sites are small forest hollows, and thus this regional component is small; so we hope it can be neglected. An alternative tool of stand-scale palynology, the Landscape Reconstruction Algorithm (LRA) (Sugita 2007b), which removes the regional pollen component from the local reconstruction, strongly depends on reliable REVEALS estimates from large lakes (Sugita 2007a). When lakes are not available in the study area, EDA is more appropriate than LRA. This can be illustrated by the following comparison with our previous study: dissimilarity coefficients between pollen- and charcoal- based EDA reconstructions are much lower than between the pollen-based LOVE model and proportions in soil-charcoal sequences (Fig. 7 in Abraham et al. 2017; Fig. 4).

Vegetation history and present communities in sandstone landscapes

For southern slopes in valleys and wide valleys at low elevations, classes 1 and 6, the results show a coniferous tree succession—Pinus (8900–6900 cal BP)-Picea (7900–3900 cal BP)-Abies (2900–900 cal BP). Class 6 also has some low proportions of broadleaved trees, but the coniferous character of both classes and low admixture of beech links the reconstructed forest from the period 2900–900 cal BP to the present community of species-poor fir forests. A concave relief is a corresponding characteristic for this vegetation type (Boublík 2013). Interestingly, this forest type presently appears primarily in the southern Czech Republic but is missing in the north.

For northern slopes and the valley bottoms, class 2, the results also suggest the presence of coniferous forests, yet with an admixture of broadleaved taxa. Fagus substitutes Quercus in the last three time-windows and forms fir-beech forest.

As regards rocks, rocky edges and non-rocky edges, classes 3 and 4, reconstruction mainly suggests the presence of Corylus, Quercus, Betula and Pinus, whereas Picea has low abundance. Despite the compositional changes during the Holocene from Corylus- to Betula-dominated forest, all species present are highly light-demanding, which corresponds to the ecological conditions on the rocks. The reconstructed forest community of both classes partly overlaps with the present communities of acidophilous pine forests (Zelený 2013); however, the presence of Corylus indicates a further overlap with communities found on deeper and nutrient-rich soils (Kovanda 1990).

For the plateaux and a convex relief at low elevations, class 5, we obtained the following tree succession—Pinus, Quercus and Corylus (8900–7900 cal BP)-Picea (7900–2900 cal BP)-Pinus and Quercus (2900–900 cal BP). The composition and the habitat of the most-recent reconstructed vegetation is likely similar to modern-day acidophilous oak woodlands (Roleček 2013) or boreo-continental pine woodlands (Chytrý 2013). The latter community is supposedly a relict element of the hemiboreal taiga (Novák et al. 2012) occupying either flat plateaux, here class 5, or the basiphilous tops of rocks, here class 3. According to our results, class 3 was never colonized by Picea, whereas Pinus woodland within class 5 evolved on sites which were previously occupied by Picea for several millennia. Class 5 in the Doksy region (Fig. 2) overlaps mainly with arenic podzols or arenic cambisols (Tomášek 2003) and we can support the possible role of Picea in their genesis (Olsson and Melkerud 1989; Augusto et al. 2002).

For both classes 7 and 8 at high elevations, tree succession starts and ends in quite a similar manner—Corylus and Betula (8900–6900 cal BP) and Abies, Fagus (3900–900 cal BP), the only difference being the co-dominant admixture that came up especially at 5900–4900 cal BP. Hill tops and other areas of convex relief at high elevations, class 7, are co-dominated by Quercus and Pinus, whereas the concave relief of high elevations, class 8, is co-dominated by Picea. Our reconstructed vegetation from the last time-windows matches present-day fir-beech forests (see animation in Supplementary Information 2).

Tree autecology based on vegetation history

The topography classification and vegetation composition of the classes already allows interpretation of ecological preferences visually; however, our NMDS (nonmetric multidimensional scaling) results show them more explicitly. Abies and Fagus are dominant in the temperate ecosystem of central Europe with very similar ecological demands (Ellenberg 1988). Fagus occupies both classes at higher elevations (classes 7 and 8), but is almost absent at lower elevations, in contrast to Abies (classes 5 and 6), due to insufficient precipitation (Koblížek 1990). Both taxa, together with Picea, occupy rather concave habitats (negative TPI250) compared to other taxa and they cover the whole gradient of topographic ruggedness and vertical distances to the stream channel network. Fagus (high VRM and VCN) does not tolerate waterlogged soils, and also late spring frosts, and thus prefers hill slopes rather than the bottom of valley basins (San-Miguel-Ayanz et al. 2016). Picea (lowest TPI250 and VCN) is able to avoid the water table with a very shallow root system (Farjon 2017) and the bottoms of sandstone valleys provide cooler climatic conditions for this mountain species (Härtel et al. 2007). The rest of the taxa from flat and convex habitats (positive TPI250) also partly cover the gradient of vertical distances to the channel network. Corylus (high VCN) does not tolerate waterlogged soils (Savill 2013), whereas trees with lower VCN can grow on peatbogs (Pinus; Skalická 1988) or tolerate periodic flooding (Quercus; San-Miguel-Ayanz et al. 2016). The positive relationship between Fagus abundance and mean annual temperature, annual precipitation and archaeological occupancy is visible even from the stratigraphic plot (Supplementary Material 1 Fig. 1), all these variables increase towards the present. However, the spread of Late Holocene taxa like Fagus or Abies was triggered by human disturbance (Giesecke et al. 2017) rather than climate, which played a role during the Early Holocene with the spread of Corylus (Giesecke et al. 2011). Abandonment during the Early Bronze Age (Kolář et al. 2022) fits with the spread of Fagus, Abies (Supplementary Material 1 Fig. 1) and Carpinus in the lowlands (Ralska-Jasiewiczowa et al. 2003). The temporal pattern of fire activity is less straightforward; however, NMDS results show a decreasing trend of fire activity with an increase of Fagus abundance. The inhabiting effect of Fagus-stands on the spread of wildfires has been proposed by Bobek et al. (2019).

Strengths, limitations and future strategies

We must admit that our climatic and archaeological data have coarser spatial resolution (1 and 33 km2) than the topographic factors (100 m2), thus we could only study the relationship of taxa to the overall trend through the Holocene. On the other hand, more precise climatic data for the last nine millennia will always contain a global climatic model downscaled to the topography. Both datasets are already included in the analysis. Archaeological findings are referenced to parish, which is more precise than the grid cell; however, the topography in sandstone areas is even finer than the parish, and thus human impact in the past with corresponding spatial precision must be estimated using some model.

Our taxa selection covers the Late Holocene well; however, for a better insight into the Middle Holocene it would merit the inclusion of Tilia, Ulmus or Fraxinus. For the EDA analysis, we only consider the landscape classification based on a digital elevation model (DEM) and its derivatives. In this respect, we neglect patterns in geology and soil types (cf. Theuerkauf et al. 2014). This can be considered a limitation for ca. 5% of the study area, where volcanic outcrops penetrate the sandstone, such as Mužský hill in the Bohemian paradise. On the other hand, soil properties depend mainly on the strong topographic gradient and overlap with our landscape classification, for example, class 7 with loess deposits on elevated plateaux (Valečka and Šebesta 2013), or class 5 with arenic podzols as mentioned above. This underlines our classification of the Digital Elevation Model and its derivatives as an appropriate approach (Man et al. 2014) for EDA modelling in sandstone landscapes.

We have shown that a spatially-explicit modelling of vegetation pattern using the charcoal record is theoretically possible and the results prove to be reasonable even with ad-hoc parameters. Nevertheless, some empirical estimates would be needed for future reconstructions. The only empirical value we have adopted from the flat boreal landscape is the circular source radius for soil charcoals (Ohlson and Tryterud 2000). However, any steep relief would require the additional testing of different shapes considering the watershed, since the main taphonomic vector is water (Clark 1988).

Vegetation modelling needs to be supported by calibration studies quantifying the relationship between the present vegetation prior to a wildfire, and the resulting charcoal assemblages produced (e.g., Scott 2010). We would also like to suggest that EDA-based vegetation models would obtain more complete results if the landscape classes were known before the site selection for soil charcoals. One can then stratify at least the spatial coverage, albeit the temporal coverage might still remain discontinuous.