Big Data Palaeoecology (BDP) approach: pollen-inferred landscape change and preindustrial demography
Recently, paleoenvironmental data derived from tree rings or ice cores have been employed to approximate changes in human economic activity related to past epidemics, as well as to warfare and climatic variability45,46. However, none of these proxies is directly related to human demography or provides a basis to estimate variation in Black Death mortality on a regional scale across Europe.
On the contrary, in preindustrial economies, rural labour availability and the spatial scale of cereal cultivation were directly related. An increase in the extent and intensity of cereal cultivation – as reflected in pollen data – would have required not only a predilection and demand for cereals, but greater availability of labour and thus population growth or significant immigration. The maintenance of existing agricultural activity, in turn, would have required relatively stable population levels.47–49 The uniform ~50% mortality postulated for the Black Death across Europe should have resulted in a large and significant decline of cereal cultivation and parallel forest regrowth across Europe, as previously demonstrated for mid-14th c. Sweden25 and singular sites in some regions of Western Europe50. This result agrees with the fact that Black Death mortality could be high among people at productive age, as illustrated for England51,52. Moreover, even in the case of England, a comparatively commercialised and adaptive rural economy in mid-14th-c. Europe, the loss of 50% of population led to a significant decline in the total area under cultivation (as documented by heterogeneous written sources)53. In Italy, another well-developed economy at that time, the expansion of large estates following the Black Death also did not compensate for the general loss of cereal productivity54. This effect, high mortality driving arable contraction, must have been yet more pronounced in more subsistence-oriented and less adaptive economies, with limited surplus production, such as in regions of the Iberian Peninsula, Germany, Sweden, and particularly East-Central Europe. Importantly, palaeoecological evidence for arable contraction may be indicative, to some extent, of not only rural population decline but also urban population decline in the region, as there is evidence in some areas, following the pandemic, of rural-to-urban migration, of country-dwellers repopulating urban centres10. Possibly less common was intraregional rural migration, as marginal lands were abandoned for better quality soils, which were more likely to remain under cultivation25,55.
BDP data collection
Existing online palynological databases (the European Pollen Database (EPD)56 www.europeanpollendatabase.net, and the Czech Quaternary Palynological Database (PALYCZ)57, https://botany.natur.cuni.cz/palycz/), as well as personal contacts of the study authors and a systematic publication search were employed to identify palynological sites in Europe reaching the required chronological and resolution quality for the study of the last millennium. In order to enable statistical analysis, we included only sites clustered in well-defined historical-geographical regions, excluding isolated sites even if the quality of a site’s data was very good. Data of sufficient quality and amount from regions for which the Black Death is well-studied, notably Central and Northern England and the Low Countries, is not presently available; to the best of our knowledge, for each of these regions there currently is not more than a single isolated site50, which does not allow for the application of statistical approaches.
In total, 261 pollen records with the average temporal resolution of 58 years and 14C-age control (or varve chronology), have been collected. The age-depth models of the sequences have been provided by authors in original publications, by the EPD or developed through the Clam package (version 2.3.4) of R software for the purpose of this study. The analytical protocol for pollen extraction and identification is reported in the original publications. The Pollen Sum includes all the terrestrial taxa with some exceptions based on the selection done in the original publications. The full list of sequences, exclusions from the Pollen Sum, age-depth models and full references are reported in Supplementary Data 1.
The taxa list has been normalized by applying the EPD nomenclature. In this respect, the general name Cichorioideae includes Asteraceae subf. Cichorioideae of the EPD and PALYCZ nomenclatures, which primarily refers to the fenestrate pollen of the Cichorieae tribe58. Ericaceae groups Arbutus unedo, Calluna vulgaris, Vaccinium and different Erica pollen types, whereas deciduous Quercus comprehends both Q. robur and Q. cerris pollen types59. Rosaceae refers to both tree and herb species of the family. Finally, Rumex includes R. acetosa type, R. acetosella, R. crispus type, Rumex/Oxyria and Urtica groups U. dioica type and U. pilulifera.
BDP summary pollen indicators
In order to connect changes visible in the pollen data to human demographic trajectories, we assembled four summary pollen indicators that describe specific landscapes related to human activity. They reflect different degrees of demographic pressure on the landscape (cereal cultivation, pastoral activities, which are less-labour intensive than cereal cultivation, abandonment and rewilding) as well as different durations of land abandonment that might have occurred post-Black Death. Our indicators account for the fact that Europe is a continent rich in natural heritage, with a wide range of landscapes and habitats and a remarkable wealth of flora and fauna, shaped by climate, geomorphology and human activity. In order to ensure uniform interpretation of the indicators, we relied on criteria that can be applied to all European landscapes regardless of their local specificity. Cereals and herding are directly related to human activities and are barely influenced by spatial differences. More complex is the succession of natural plants with their ecological behaviour and inter-species competition. For this reason, we relied on existing quantitative indicators of plant ecology.
The Ellenberg L – light availability indicator60 provides a measure of sunlight availability in woodlands and consequently of tree-canopy thickness, reflecting the scale of the natural regeneration of woodland vegetation after cultivation or pasture activities61. Nonetheless, ecological studies have suggested that geographic and climatic variability between different European regions can influence the Ellenberg indicator system62–65. The original indicators were primarily designed for Central Europe58, but several studies developed Ellenberg indicators for other regions, reflecting the specific ecology of the selected taxa (British Isles66; Czech Republic67; Greece68; Italy69; Sweden70). Plants with a L values between 5 and 8 are listed in fast succession indicator, the ones with L values ranging from 1 to 4 are included in slow succession indicator. The result is the following list:
1) Cereals: only cultivated cereals have been included: Avena/Triticum type, Cerealia type, Hordeum type, Secale. 2) Herding includes pastoral indicators linked to the redistribution of human pressure: Artemisia, Cichorioideae, Plantago lanceolata type, Plantago major/media type, Polygonum aviculare type, Rumex, Trifolium type, Urtica, Vicia type. 3) Fast Succession comprises indicators of relatively recent reforestation of cultivated land after abandonment: Alnus, Betula, Corylus, Ericaceae, Fraxinus ornus, Juniperus, Picea, Pinus, Populus, deciduous Quercus, Rosaceae. 4) Slow Succession includes indicators of secondary succession established after several decades of abandonment: Abies, Carpinus betulus, Fagus, Fraxinus, Ostrya/Carpinus orientalis, Quercus ilex type.
In order to validate the indicators overcoming the regional limits of Ellenberg values, a different subdivision has been provided following the Niinemets and Valladares shade tolerance scale for woody species of the Northern Hemisphere71. The subdivision of taxa in the Fast and Slow succession indicators remains the same with only three changes: Fraxinus ornus and Picea move from Fast to Slow succession and Fraxinus from Slow to Fast succession. Extended Data Fig. 1 and 2 show that the two groupings yield the same results, which confirms the reliability of our indicators. There is only one clear exception (Russia), with one more region where smaller-scale diversion occurs for only one indicator, Slow Succession (Norway). The different indicator behaviour results from the different attribution of Picea in our two sets of succession indicators: at high latitude, Picea characterizes the final stage of the ecological succession and hence its different attribution results in different summary indicator values in Russia for the two stages of ecological succession, fast and slow.
Please note our summary indicators are not designed to reflect the entirety of the landscape and reconstruct all of its different components. Rather, they are a means of approximating changes in the landscape related to the types of human activities, and their intensity, as much as they relate to demographic changes in human populations using and inhabiting these landscapes.
BDP analytical methods: statistical and spatial
To control for local specificity, pollen percentages of every taxon from each pollen site were standardised. From the taxa percent in a given year the arithmetic mean calculated for the observations from the period 1250–1450 was subtracted and the result divided by the standard deviation for the 1250–1450 period. Standardized taxa results were assembled for each site into four BDP summary indicators. Since each indicator has different numbers of taxa, the sum of standardized taxa values calculated for a given year and site was divided by the number of taxa in the indicator. For the purposes of replication, this standardised pollen dataset, comprising the four indicators for each sample from each site, is available as Supplementary Data 1.
This dataset has been analysed in two ways, statistically and spatially.
For the statistical approach, standardized regional indices of landscape transformation were created for each region by calculating the average value for all sites within the region, for each of the subperiods analysed in the study (1250–1350 and 1351–1450; 1301-1350 and 1351-1400; 1325-1350 and 1351-1375). Differences between means for each subperiod were measured by the use of the bootstrapping based on 10,000 resamples. The 90% and 95% confidence intervals were estimated with the bias-corrected and accelerated method (BCa)72. These results are visualized in Fig. 3 for the comparison of the sub-periods of 1250–1350 vs 1351–1450, and in Extended Data Fig. 7 and 9 for the comparison of the sub-periods of 1300-1350 vs 1351-1400 and 1325-1350 vs 1351-1375, respectively.
For the spatial approach, we employed the Bayesian model AverageR developed within the Pandora & IsoMemo initiatives (https://pandoraapp.earth/) to map the distribution of pollen indices across Europe. AverageR is a generalized additive model that has been previously described73. It relies on a thin plate regression spline75 to predict new, unseen data using the following model:
Yi=g(longitude,latitude)+εi
where
Yi: independent variable for site i
g(longitude, latitude): spline smoother
εi ~ N(0, σε): error term
The spline smoother can be written as X*β where X is a fixed design matrix and β is the parameter vector. Surface smoothing is controlled by employing a Bayesian smoothing parameter estimated from the data and trades-off bias against variance to make the optimal prediction75. This parameter β is assumed to follow a normal distribution: β ~ N(0, 1 /δ*λ*P), where P is a so-called penalty matrix of the thin plate regression spline, which penalizes second derivatives75. The δ parameter is by default set to 1 but this can be adjusted to suit smoothing needs for each application. In our study δ was set at 0.9 to match the preferred spatial scale of analysis for our dataset (c. 250 to 500km).
AverageR was employed to generate smoothed surfaces for three sets of temporal bins (1250–1350 vs 1351–1450, as well as 1300-1350 vs 1351-1400 and 1325-1350 vs 1351-1375) and for the four BDP indicators (Extended Data Fig. 6, 8 and 10). For the same indicator the difference between the two temporal bins was plotted (Fig. 3; Extended Data Fig. 7 and 9).
Data Availability: All data generated or analysed during this study are included in this published article (as the Supplementary Data).
References – methods
45.
Ljungqvist, F. C. et al. Linking European building activity with plague history. Journal of Archaeological Science 98, 81–92 (2018).
46.
McConnell, J. R. et al. Pervasive Arctic lead pollution suggests substantial growth in medieval silver production modulated by plague, climate, and conflict. PNAS 116, 14910–14915 (2019).
47.
Izdebski, A. A rural economy in transition: Asia Minor from Late Antiquity into the Early Middle Ages. (Taubenschlag Foundation, 2013).
48.
Roberts, N. et al. Europe’s lost forests: a pollen-based synthesis for the last 11,000 years. Scientific Reports 8, 716 (2018).
49.
Bevan, A. et al. The changing face of the Mediterranean – Land cover, demography and environmental change: Introduction and overview. The Holocene 095968361982668 (2019) doi:10.1177/0959683619826688.
50.
Yeloff, D. & Van Geel, B. Abandonment of farmland and vegetation succession following the Eurasian plague pandemic of ad 1347–52. Journal of Biogeography 34, 575–582 (2007).
51.
Razi, Z. Life, marriage, and death in a medieval parish: economy, society, and demography in Halesowen, 1270-1400. (University Press, 1980).
52.
DeWitte, S. N. Age Patterns of Mortality During the Black Death in London, A.D. 1349–1350. J Archaeol Sci 37, 3394–3400 (2010).
53.
Broadberry, S. N. British economic growth, 1270-1870. (University Press, 2015).
54.
Cortonesi, A. L’agricoltura italiana fra XIII e XIV secolo: vecchi e nuovi paesaggi. in Medioevo delle campagne: rapporti di lavoro, politica agraria, protesta contadina. (eds. Cortonesi, A. & Piccinni, G.) 15–56 (Viella, 2011).
55.
Abel, W. Agrarkrisen und Agrarkonjunktur. Eine Geschichte der Land- und Ernährungswirtschaft Mitteleuropas seit dem hohen Mittelalter. (Parey, 1978).
56.
Fyfe, R. et al. The European Pollen Database: past efforts and current activities. Vegetation History and Archaeobotany 18, 417–424 (2009).
57.
Kuneš, P., Abraham, V., Kovářík, O., Kopecký, M., & PALYCZ contributors. Czech Quaternary Palynological Database – PALYCZ: review and basic statistics of the data. Preslia 81, 209–238 (2009).
58.
Florenzano, A., Marignani, M., Rosati, L., Fascetti, S. & Mercuri, A. M. Are Cichorieae an indicator of open habitats and pastoralism in current and past vegetation studies? Plant Biosystems - An International Journal Dealing with all Aspects of Plant Biology 149, 154–165 (2015).
59.
Smit, A. A scanning electron microscopical study of the pollen morphology in the genus Quercus. Acta Botanica Neerlandica 22, 655–665 (1973).
60.
Ellenberg, H. et al. Zeigerwerte von Pflanzen in Mitteleuropa. Scripta Geobotanica 18, 1–258 (1992).
61.
Dzwonko, Z. Assessment of light and soil conditions in ancient and recent woodlands by Ellenberg indicator values. Journal of Applied Ecology 38, 942–951 (2001).
62.
Diekmann, M. & Lawesson, J. E. Shifts in ecological behaviour of herbaceous forest species along a transect from northern central to North Europe. Folia Geobotanica 34, 127–141 (1999).
63.
Gégout, J.-C. & Krizova, E. Comparison of indicator values of forest understory plant species in Western Carpathians (Slovakia) and Vosges Mountains (France). Forest Ecology and Management 182, 1–11 (2003).
64.
Hájková, P., Hájek, M., Apostolova, I., Zelený, D. & Dítě, D. Shifts in the ecological behaviour of plant species between two distant regions: evidence from the base richness gradient in mires. Journal of Biogeography 071103055558002-??? (2007) doi:10.1111/j.1365-2699.2007.01793.x.
65.
Wasof, S. et al. Ecological niche shifts of understorey plants along a latitudinal gradient of temperate forests in north-western Europe. Global Ecology and Biogeography 22, 1130–1140 (2013).
66.
Hill, M., Mountford, J., Roy, D. & Bunce, R. Ellenberg’s indicator values for British plants. ECOFACT Volume 2 Technical Annex. (Institute of Terrestrial Ecology, 1999).
67.
Chytrý, M., Tichý, L., Dřevojan, P., Sádlo, J. & Zelený, D. Ellenberg-type indicator values for the Czech flora. Preslia 90, 83–103 (2018).
68.
Böhling, N., Greuter, W. & Raus, T. Zeigerwerte der Gefäßpflanzen der Südägäis (Griechenland). Indicator values of the vascular plants in the Southern Aegean (Greece). Braun-Blanquetia 32, 1–106 (2002).
69.
Pignatti, S., Menegoni, P. & Pietrosanti, S. Biondicazione attraverso le piante vascolari. Valori di indicazione secondo Ellenberg (Zeigerwerte) per le specie della Flora d’Italia. Braun-Blanquetia 39, 1–97 (2005).
70.
Diekmann, M. Use and improvement of Ellenberg’s indicator values in deciduous forests of the Boreo-nemoral zone in Sweden. Ecography 18, 178–189 (1995).
71.
Niinemets, Ü. & Valladares, F. Tolerance to Shade, Drought, and Waterlogging of Temperate Northern Hemisphere Trees and Shrubs. Ecological Monographs 76, 521–547 (2006).
72.
Hall, P. Theoretical Comparison of Bootstrap Confidence Intervals. The Annals of Statistics 16, 927–953 (1988).
73.
Cubas, M. et al. Latitudinal gradient in dairy production with the introduction of farming in Atlantic Europe. Nature Communications 11, 2036 (2020).
74.
Wood, S. N. Thin plate regression splines. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 65, 95–114 (2003).
75.
Groß, M. Modeling body height in prehistory using a spatio-temporal Bayesian errors-in variables model. AStA Adv Stat Anal 100, 289–311 (2016).