Introduction 

Wood charcoal production is one of the oldest forms of anthropic forest use in most countries of temperate regions ([3]) and the main source of energy since the Iron Age until the 19^th century ([10]). Charcoal production was carried out in the so-called “charcoal kilns” or “charcoal hearths”, which are small terrace-like platforms (35-50 m²) with circular or elliptical shape, commonly established in wind-sheltered forest sites not far from water courses, and connected by ad-hoc footpaths across hill and mountain slopes ([13], [4], [34]). A detailed description of the different production phases and the procedures used for charcoal kiln construction are reported in Landi & Piussi ([34]), Powell ([51]) and Giorgerini ([26]), and are briefly illustrated in Fig. S1 (Supplementary Material).

In northern and central Europe this practice was generally abandoned in the 19th century due to the rapidly increasing and widespread use of coal ([21]), while in most Mediterranean countries the production of wood charcoal increased during the industrial revolution because of the lack of other fuel sources. It generally vanished only in the 1950s, though in some mountain areas it is still in use. As a result, remnants of charcoal platforms are nowadays commonly found in many forest landscapes ([10], [47]). They are characterized by clear alterations of color and texture of the soil material due to the charred woody remains ([43]). In fact, the condensed aromatic structure of wood charcoal allows fragments and particles to persist in soils and other sedimentary records over millennial time-scales ([16], [61]). Anthropogenic deposits of charcoal dating back to the Neolithic period have been documented in Germany and Italy ([19], [38]) and very old samples (> 8000 years BP) originating from wildfires have been found almost unaltered in forest soils of the Italian Alps ([55]).

Several anthracological studies have been carried out to reconstruct the former tree species composition of forests through the identification of the woody taxa used for charcoal production, thus highlighting the usefulness of charcoal platforms in detecting vegetation changes through time ([43], [36], [45], [37]). However, the present-day vegetation of these sites has been rarely investigated in Europe ([59], [14]), and no studies focused on the tree regeneration and forest recolonization processes in charcoal platforms.

Successful recolonization of charcoal platforms by woody species is expected because of the apparently favorable conditions such as the low degree of rockiness, flat terrain and often a higher availability of light ([14]). On the other hand, adult trees are usually absent in most of the platforms, even in those abandoned several decades ago. This suggests that tree regeneration and establishment is hampered by unknown persistent biotic or abiotic factors, which represent a legacy of the past forest use still affecting ecosystem dynamics at the “microhabitat” scale.

In order to test this hypothesis, we analyzed the diversity, density and structural parameters of trees and shrub species in the regeneration layer (< 4 m) of these platforms, and compared them with control sites located in the adjacent forest environment. Differences were then related to variation in light, soil conditions and herb cover measured at the same sites. Using such a paired design and an extensive sampling in three major forest types of central Italy, we could provide the first insights into the long-term effects of the old charcoal platforms on woody species regeneration in the Mediterranean area, and on the possible driving factors.

  Material and methods 

Study area

This study was performed in the forests of Tuscany (central Italy - Fig. 1). The study area includes three main climate and forest types, according to the EEA classification system ([22]): (1) meso-mediterranean evergreen forest dominated by Quercus ilex L. and sclerophyllous shrubs, along the Tyrrhenian coast (hereafter: “sclerophyllous” forests); (2) supra-Mediterranean thermophilous mixed communities dominated by deciduous oaks (Quercus cerris L., Q. petraea (Matt.) Liebl., Q. pubescens Willd.) in the hilly areas of central Tuscany (hereafter: “oak” forests); (3) montane forests with beech (Fagus sylvatica L.), on the Apennine chain (Casentino and Mugello areas) and the volcanic massif of Mt. Amiata (hereafter: “beech” forests). Mean annual rainfall and temperature in this area vary from 780.7 mm and 15.4 °C along the coast to 1.641 mm and 8.7 °C on the Apennines (period 1961-1990 - [53]). The investigated area is characterized by a variety of geolithological formations and soil conditions, but cambisols are the prevalent type according to the Soil Atlas of Europe ([23]).

Fig. 1 - Geographic location of the forest stands and the investigated charcoal platforms in Tuscany (central Italy).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Field sampling design and data collection

An extensive search for abandoned charcoal platforms was carried out with ad-hoc field trips in representative areas of the three main forest types described above. Charcoal platforms were easily recognized due to the characteristics of the ground surface ([4], [34], [39]). One-hundred fifty-four different platforms were identified, and their position recorded with a GPS device. At each site, altitude, slope, aspect, soil type and tree species composition of the adjacent forest stands were also recorded. Sixty-one out of the above 154 sites were randomly selected among those that were not affected by recent anthropogenic or animal disturbance: (i) 18 “sclerophyllous“ sites in 4 geographically distant forest stands; (ii) 22 “oak” sites in 3 stands; (iii) 21 “beech” sites in 3 stands. In the two areas (Volterra hills and Farma-Merse) both sclerophyllous and oak stands were present (Fig. 1, Tab. 1). The main geographical and environmental variables of these sites are given in Tab. 1. All but two sites are included in protected areas of the Natura 2000 network, Nature Reserves and National Parks. Charcoal production in all sites lasted for centuries and was abandoned at least 60 years ago, while the coppice-with-standards management for firewood production was abandoned later (ca. 45-50 years ago) in both oak and beech forests. In the sclerophyllous stands, this management type was abandoned more recently or it is still in use (Carrari & Selvi, personal observation).

At each selected site, a 3×3 m quadrat (hereafter: “kiln plot”, KP) was established in the center of the charcoal kiln platform. The relatively small size and the general shape of the platforms did not allow the adoption of larger quadrats, which would have included the kiln-forest edge. Light conditions (Photosyntetic Active Radiation, PAR), herbaceous cover and soil factors were measured following the protocol described in Carrari et al. ([14]); observations on the soil profile were also conducted in each quadrat.

Tree regeneration was analyzed in two distinct layers: (I) the understorey vegetation layer (< 1.3 m - hereafter: “UVL”) and the established tree regeneration layer (1.3-4 m - “ETL”). In the UVL, seedlings of all tree and shrub species were identified and their maximum height was recorded. The ground cover percentage of each of these species was quantified as the projection of their crown foliage on the soil surface. According to Pignatti ([48]), we considered as trees all woody angiosperms or gymnosperms belonging to the Raunkiaer’s life forms “Pscap” (scapose phanerophyte) and “Pcaesp” (caespitose phanerophyte). In addition, the number of seedlings (height < 1.3 m) of the dominant species (all belonging to the Fagaceae family) was counted within the quadrat in a sub-sample of five randomly selected platforms for each forest type (15 plots in total - Tab. S4 in Supplementary Material).

In the “ETL”, we recorded the species richness of trees and shrubs and the following structural parameters: number of trees/stools, number of stems per stool, diameter at breast height (dbh, cm) of each stem (> 0.5 cm), and mean height of the stems for each stool.

For each of the 61 kiln plots, the same measurements mentioned above were taken in a control plot (CP) of the same size randomly established in the stand adjacent to the kiln platform at a distance of 10-20 m from the edge of the KP (depending on local topographic and ground conditions); downhill locations were excluded to avoid potential charcoal “contamination” by runoff. This allowed to minimize the variation of the major environmental factors (altitude, parent rock material, slope inclination and slope aspect) between each kiln plot and its respective control plot.

Data analysis

All statistical analyses were performed in R 3.1.2 ([52]). Soil factors, PAR values and herbaceous cover were compared between KP and CP using the non-parametric Mann-Whitney U test, due to the violation of assumptions of normality and homogeneity of the variance (“wilcox.test” function). Differences in woody species composition of the UVL between KP and CP in each forest type were assessed using the non-metric multidimensional scaling (NMDS) based on the cover-weighted Bray-Curtis dissimilarity index (“metaMDS” function of the “vegan” package - [46]). Significance of the differences were tested by permutational multivariate analysis of variance (PERMANOVA) using the “adonis” function with 999 permutations. To verify that such differences were related to the effect of the factor kiln/control (i.e., to dissimilarities in composition between kiln and control plots) and not to a dispersion effect (i.e., dissimilarities within each of the two plot types), we tested for multivariate homogeneity of dispersion using the “betadisper” function, a multivariate analogue of the Levene’s test for homogeneity of variances ([2], [46]).

Density of seedlings of the locally dominant tree species was compared between KP and CP using the non-parametric Mann-Whitney U test. The effects of the charcoal kiln platforms on species richness (SR) at the habitat level (γ-diversity of KP and CP) and plot level (mean SR for KP and CP) in the UVL and in the ETL were also evaluated.

A mixed model approach was adopted in order to examine all tree regeneration parameters (seedling density of dominant species, SR in UVL and ETL and structural parameters) as response variables to the following predictor variables: forest type (levels: sclerophyllous, oak, beech), environmental variables (altitude, parent rock, slope inclination and slope aspect), soil and light factors (N, C, pH and PAR) and herbaceous cover. First, predictor variables were tested for collinearity using Pearson’s and Spearman’s correlation analysis to obtain quantitative and ordinal correlation coefficients, respectively. The association between nominal variables or between a nominal and a continuous variable was tested by chi-square or Student’s t-test, respectively. As a result, parent rock material, altitude and pH (associated with forest type) were excluded from further analyses.

The variable “forest stand” was included in the model as random factor to account for possible effects of environmental and/ or management-related factors on the tree regeneration variables specific to each forest area.

Tree seedling density and SR in the understory were initially analyzed using the following model (eqn. 1):

where y is the response variable, α is the intercept, β is the parameter to be estimated for each factor, FT is the forest type, HC is the herb cover, PAR is the photosynthetic active radiation measured, C is the soil carbon content, N is the soil nitrogen content, SA is the slope aspect, SI is the slope inclination, 1|FS indicates the forest stand as random factor and ε is the error. The model was tested using the “glmer” function in R with a Poisson error distribution and the “loglink” function of the “lme4” package ([6]). For SR in the ETL, HC was excluded from the model as not relevant for woody species growth in this layer; light was also excluded since PAR values were measured at 1 m above the forest floor.

Starting from these full models, we looked for model parsimony following the approach of Zuur et al. ([62]). First, we deleted the random variation across sites while keeping the fixed effect term, using the “glm” function from the “stats” package in R, with a Poisson error distribution, log link and parameter estimation via maximum likelihood. The model yielding the lowest value of the Akaike’s Information Criterion (AIC - [1]) was considered as the best fitting to the data. Once an optimal random structure was found, we searched for the optimal fixed effect structure by comparing the AIC of models with the same random effect structure, but a different fixed effect structure (here parameter values were maximum likelihood estimates). In case of over-dispersion, the standard errors were corrected using a quasi-GLM model ([62]). Accordingly to Zuur et al. ([62]) models were validated looking at response residuals (observed minus fitted values - Fig. S2 and Fig. S3 in Supplementary Material).

Using a similar mixed model approach we tested the effects of the relevant predictors and their interactions on the structural parameters (number of trees/stools, number of stems per plot and their basal area and mean height in the ETL) using the following model (eqn. 2):

The number of trees/stools and the number of stems were fitted using the “glmer” function with a Poisson error distribution, as for the species richness models. Basal area and mean height were fitted allowing for random variation across “forest stands” using the “lmer” with a Gaussian error distribution. In this case, fixed effect models were tested using the “gls” from the “nlme” package ([49]) with a Gaussian error distribution, and parameter estimation was calculated with a restricted maximum likelihood. Basal area was preliminary log-transformed to match the requirements of normality and homoscedasticity.

For the optimal models selected for each variable, we calculated the R-squared (R2), i.e., the fraction of the total variation in the response variable accounted for by the model. For models with fixed effects only, the adjusted R2 of the linear model was reported; for models that also contained random effects, a conditional R2 was calculated using the “MuMIn” package ([7]), according to Nakagawa & Schielzeth ([44]).

  Results 

Soil, light and herbaceous cover

All kiln plots (KP) were characterized by a well distinct layer, 10-46 cm thick, of black soil formed by the accumulation of charcoal, with fragments of medium-large size. This layer was sharply delimited from the underlying mineral soil layer, whose color depended mainly on the parent rock material. The total carbon content in the charcoal layer of KP was almost two-fold higher than in CP (control plots) in sclerophyllous and oak forests, and significantly higher in KP of beech forests, though to a lower extent (Tab. 2). The same trend was generally recorded for N, but such differences were significant for oak forests only. Regarding pH, higher values were found in KP of the beech forest, while no significant differences were detected for the other two forest types. Both light availability (PAR) and herbaceous cover (HC) were higher in KP compared with CP in oak forests, while no significant differences could be found in the other two forest types (Tab. 2).

Seedlings density

In general, seedling density of the dominant tree species (Fagaceae) increased with light availability, decreased with herb cover, and was affected by forest type and by its interaction with the above variables. Moreover, density decreased on steeper slopes and increased on south-facing slopes. No significant differences between KP and CP were found, with the exception of Quercus pubescens in sclerophyllous forests (p = 0.04 - Tab. 3). However, total C content had a weak negative effect on seedling density, while the latter was positively influenced by the N content (+1.68). The random effect of forest stands was not included in the selected model (Tab. 4), which explained 57% of the total variation.

Composition and species richness

Concerning the composition of woody species in the UVL, NMDS analysis yielded different results for the three forest types. Significant differences between KP and CP were found in oak forests (p_perm = 0.008), with a similar level of compositional variation within these two plot types (p_disp = 0.541 - Fig. 2). Beech and sclerophyllous understorey showed no differences in species composition between the two plot types.

Fig. 2 - Scatter plot from NMDS based on Bray Curtis dissimilarity index of woody species in the understorey (UVL) of oak forests. Differences in woody species composition of the understory of kiln (KP) and control plots (CP) were significant. p_perm indicates the combined significance of the location and dispersion effect, based on PERMANOVA with 999 permutations; p_disp indicates the significance of the dispersion effect.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

In total, 45 woody species (26 trees and 19 shrubs) were recorded in the UVL of the 122 (= 61 × 2) plots. In this layer, γ-diversity was always higher in KP compared to CP for all three forest types (22 vs. 14 in sclerophyllous forests; 28 vs. 22 in oak forests; 13 vs. 9 in beech forests - Fig. 3a). In sclerophyllous forests, nine species were unique to KP, while only one to CP; in oak forests, eight species were unique to KP and two to CP; in beech forests, four species were unique to KP (Acer platanoides, Castanea sativa) and one to CP (Prunus avium - Tab. S1 in Supplementary Material). Species richness in the UVL was significantly higher in KP compared with CP in oak forests; the same trend was also recorded for the other two forest types, though without significant differences (Fig. 3a).

Fig. 3 - Differences in plot-level species richness between kiln (KP) and control (CP) plots in the UVL (A) and ETL (B) layers of each forest type. Error bars represents the standard error. Total species richness (γ) values of each plot and forest type are reported above each bar. (***): p<0.001; (ns): p>0.05.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The ETL included a total of eight tree and nine shrub species. In contrast to the UVL, γ-diversity was larger in CP than in KP (Fig. 3b). In the KP of sclerophyllous forests, this layer included four species (present in 16.7% of the plots - Tab. S1 in Supplementary Material), compared to seven species in the total sample of CP (Fig. 3b). In oak forests, no species were recorded in KP, while eight were present in this layer in 72.7% of the CP plots (Fig. 3b, Tab. S2). No tree species were found in the ETL of KP in beech forests (Fig. 3b), while beech was present in 14.3% of the CP (Tab. S3). Consequently, species richness at the plot-level in ETL was always significantly higher in CP compared with KP in all forest types (Fig. 3b).

Mixed model analysis showed that the random variation across the examined forest stands was not relevant for species richness in the two layers, while the fixed effect of forest type was always significant, explaining 66% of the total variation in the UVL and a lower proportion (26%) in ETL (Tab. 4). In the UVL, C content affected positively SR, while N had a negative effect. Moreover, light and herb cover also increased SR, though to a varying extent depending on the forest type (Tab. 4). In the ETL, SR increase was mainly due to N (+3.1) and slope inclination (+0.73), with differences among forest types, while its decrease was mainly associated with C content (Tab. 4).

Structural parameters

Regarding the structural variables in the ETL, large differences were found between the three forest types. As expected, the density of trees and shrubs in the sclerophyllous and oak stands was considerably higher than in the beech stands (Fig. 4). Remarkable differences were found comparing KP with CP. No stools or stems were present in KP of oak and beech forests (Fig. 4a, Fig. 4b). In sclerophyll forests, though a few species were present, the mean basal area of KP was always much lower than in CP (Fig. 4c). The height of these species was also significantly lower in KP (Fig. 4d).

Fig. 4 - Differences in structural variables of the ETL layer (1.30-4 m) between kiln (KP) and control plots (CP) for each forest type. (A): mean number of stools per plot; (B): mean number of stems per plot; (C) mean basal area (G, cm²); (D) mean height of tree stems (m). Error bars represent the standard error. (***): p<0.001; (ns): p>0.05.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Looking at the effect of the measured variables on such parameters, FT (forest type), C and N were the main predictors. In particular, C content was the main factor negatively affecting the number of stools and stems, especially in oak forests. Contrastingly, N content was strongly and positively related with number of stools and stems, both alone and in the interaction with forest type, while it was not significant for the other two variables. Moreover, slope inclination and aspect had other minor effects on structural variables (Tab. 5). The models for the number of stools and stems explained up to the 42% and 44% of the variation, respectively, while the models for the mean basal area and height explained not more than 30% of the total variation. The best model for number of stools and stems included the random factor of the forest stand location (Tab. 5).

  Discussion 

The results of this study, obtained by analyzing three distinct forest types, demonstrated that former charcoal platforms represent a peculiar microhabitat for the regeneration of woody species.

The effect of kiln platforms on the γ-diversity in the understorey layer was clearly positive, especially in oak forests, where richness was also enhanced at the plot-level. In the latter forest type, kiln platforms also affected the composition in woody species, in line with previous studies on the understorey in mixed oak stands of the same area ([14]), as well as in beech forests of Germany ([59]). Interestingly, no selection effect seems to occur at this early stage of regeneration, since the species recorded have different ecological requirements and functional traits, especially in thermophilous oak forests ([15]). Regeneration of pioneer shrub species such as Cytisus scoparius, Crataegus monogyna and Prunus spinosa occurred in the understorey mixed with those of early- and late-successional trees, such as Fraxinus ornus and Quercus ilex, respectively, depending on the forest type.

Recent findings have proven that kiln platforms are microhabitats which favor species diversity in the understory ([14]). The results of our models revealed that various factors significantly affect woody species composition in the understory, among which herb cover. Facilitation mechanisms among plants ([12]) may explain the positive influence of herbaceous cover on tree species richness in platforms of oak forests. On the other hand, herb cover negatively influenced the seedling density of dominant trees, likey due to the higher competitive ability of herbaceous plants for soil nutrients ([40]). Moreover, other important factors affecting woody species richness on the platforms are the higher C content in the soil and the increased light availability, the latter being a major driver of understory diversity and productivity in European forests ([5], [24]). Indeed, charcoal platforms may initially act as small canopy gaps, increasing the light available on the forest ground and positively affecting the early regeneration stages ([50], [8]). It is well documented that canopy gaps can provide recruitment opportunities for tree seedlings and thus increase the number of species, promoting forest diversity ([11]). In beech forests, the soil pH was higher in our platforms than in control plots, confirming previous findings in other kiln sites in North America and Europe ([42], [41], [60], [59], [20]). This is likely another positive factor for species richness, also based on evidence from deciduous woodlands in central Europe ([17], [5]). On the contrary, the higher N availability resulted to negatively affect woody species richness, since most plants are outcompeted by few species able to better exploit this resource ([18]).

Regarding the effects of kiln platforms on the established regeneration layer (1.3-4 m), the response of woody species was generally largely negative. Remarkably, a similar detrimental influence on the woody species was found in the old charcoal platforms of south-eastern Pennsylvania ([42]) and the Appalachian mountains in north America ([60]). In our study, both total (γ-diversity) and plot-level species richness were generally lower on the platforms than in the control plots for both oak and sclerophyllous forests, clearly indicating less suitable conditions for trees and shrubs at later stages of their development.

A similar negative effect was detected for most structural variables. Although seedling density in the understorey was similar between platforms and control plots, the number of stools and stems per plot above 1.3 m (ETL) was much lower in platforms than in the adjacent stands of sclerophyllous forests, or even stools and stems are completely lacking in oak and beech forests.

Overall, these findings suggest that such contrasting “kiln platform effect” is not simply a “forest gap effect”, but is driven by peculiar factors that are likely associated with the long-lasting production of wood charcoal at the same sites. According to previous works ([42], [41], [14]), these include abiotic conditions, most likely in the soil, and possibly other components of biotic nature. The results of our models allow to exclude the increased soil nitrogen content as a negative factor for species richness, as well as the number of stems and stools. Instead, there is evidence for a general negative effect of soil carbon content on species richness and most structural parameters, with differences among forest types. Hence, there is a contrasting response of woody species to increased soil C (positive in the UVL and negative in ETL) and N content (negative in the UVL and positive in ETL).

The high C content in the platform soil is mainly due to the large proportion of charred fragments, as recently demonstrated for a hearth in the Eastern Alps ([20]). Thus, the negative effect on the established regeneration layer are likely associated with the charcoal accumulation in the soil, rather than to its carbon content per se. Consequences of charcoal accumulation on tree growth are still poorly known. A positive influence has been reported by authors who support the “biochar” practice to promote forest restoration ([54]). Indeed, previous studies suggested that the higher availability of some macro- and micronutrients in these soils could exert a positive effects on microbial communities, as in the case of the south American Terra Preta soils ([28], [27], [35], [31]) and charcoal hearths in the Alpine region ([20]). In this context, our results are in line with the positive effect of charcoal on nutrient availability (especially due to the reduced leaching of nitrogen) reported in several agricultural biochar studies ([9], [32], [56]).

On the other hand, less favourable effects of chacoal on plant growth have been reported in the literature, such as an altered nutrition as a consequence of the poor development of mycorrhizal fungi. For example, experimental “biochar” studies on agricultural crops showed its negative effect on arbuscular mycorrhizal fungi ([58]), and a decrease in both bio-available organic C and N in their ectomycorrhizal system ([57]). Further, Gaur & Adholeya ([25]) found that the biochar media hampered the P uptake by host plants, indicating that in some cases charcoal may reduce the formation of mycorrhiza by decreasing the nutrient availability or creating unfavorable nutrient ratios in the soil ([57]).

Another consequence of charcoal addition is the alteration of the soil water retention. Although charcoal is known to positively affect soil water retention due to its porosity, the real availability of charcoal-adsorbed water to plants still needs to be assessed ([30]).

Finally, we cannot exclude that other factors not considered in this study could patly explain the “kiln platform effect”, such as the repeated soil heating on the platforms. It is well documented that slash pile burning in forests can cause cascading effects on soil structure (e.g., porosity) and chemical properties (e.g., pH, nutrient availability) and has a negative impact on the populations of arbuscular mycorrhizal fungi ([33]) and microbial communities ([29]). In addition, Mikan & Abrams ([42]) suggested the presence of potentially harmful soluble salts in the soil of charcoal platforms that may cause physiological drought. In our study, the presence of drought-tolerant species in the overstorey of the charcoal sites (e.g., Fraxinus ornus, Arbutus unedo) and the lack of some more mesophytic species found in the overstorey of the controls (e.g., Acer campestre, Carpinus betulus) may lend circumstantial support to this hypothesis.

  Conclusions 

The extensive networks of old charcoal platforms in Mediterranean forests provide a natural experimental setting to investigate the legacy effects of an ancient form of forest use on tree regeneration dynamics. In this work, we demonstrated a general positive influence of kiln platforms on woody species at their first stage of regeneration in three different forest types. On the other hand, we found that the further growth and development of trees is negatively affected by modified soil and environmental conditions that prevent their access and establishment into the higher layers, thus causing a substantial lack of forest recolonization. This effect hinders or at least slows down the recolonization of even the oldest kiln platforms, unlike in forest gaps of different origin. Whether charcoal accumulation per se has a negative impact on tree growth remains an open issue which should be investigated on different species with an experimental approach.

  List of abbreviations 

KP: charcoal kiln plot; CP: control plot; UVL: understorey vegetation layer; ETL: established tree regeneration layer; SR: species richness.

  Acknowledgements 

Authors wish to thank Giacomo Certini (University of Firenze, Italy), Giovanni Mastrolonardo (University of Bologna, Italy) and Davide Ascoli (University of Torino, Italy) for discussion and comments that improved the first version of the manuscript; Giovanni Quilghini (Corpo Forestale dello Stato, UTB Pratovecchio, AR, Italy) for allowing field sampling in the Badia Prataglia forest, and Andrea Barzagli for his help during the field work. This study was supported by funds from the Italian Ministry of University and Scientific Research and the University of Firenze to FS; the research of KV is supported by the ERC Consolidator Grant 614839 - PASTFORWARD.

  References