Modern hazelnuts were collected on 8th and 10th September 2021 and 16th and 17th September 2022 at three sites in Scania, southern Sweden: Linnebjer nature reserve, Dörröd nature reserve and Holma (Figure 1; Supplementary Table 1). Nuts were collected in between five and ten locations (sub-sites) at each site from individual trees or from the ground below trees. A minimum of three nuts were collected at each sub-site. Where nuts were collected directly from the tree, or there was only one hazel tree in the location, we could use the nuts to establish intra-tree variability. Leaf Area Index (LAI) was measured with a LAI 2000 Plant Canopy Analyzer (PCA, Li-Cor, Lincoln, NE, USA) at each sub-site to determine the density of the canopy. Sample codes comprise a site identifier (L = Linnebjer, D = Dörröd, H = Holma), year (absent = 2021, 22 = 2022), subsite/location identifier (1–10) and nut identifier (1–10). Three of the nuts were cut into three, perpendicular to the axis of growth (see Figure 3), to assess intra-nut variability and were labeled as apex (u), mid (m) and base (o) sections. Precipitation, relative humidity, cloud cover, sunshine time and solar radiation for 1st January to 31st August and 1st May to 31st August 2021 and 2022 were extracted from data provided by the SMHI, Swedish Meteorological and Hydrological Institute (https://www.smhi.se/data/meteorologi/ladda-ner-meteorologiska-observationer; Table 1).

Ten hazelnut shell fragments each from the sites of Slabälta 1, Rönneholm 6:1, Rönneholm 10:3 and Ringsjöholm (Figure 1) were selected for isotope analysis. These shell fragments were preserved by waterlogging and had minimal signs of charring. The sites have been radiocarbon dated to the Mesolithic, between c. 7000 and 5500 BCE, spanning the late Maglemose and Kongemose cultures (Table 2; Sjöström, 1997, 2004, 2013; Larsson and Sjöström, 2010). Charred hazelnut shell fragments were also selected from sites in southern Sweden dated to between the Neolithic and Iron Age (Figure 1; Table 2). The cross-sections of nutshell fragments from Uppåkra, Uppåkra 2:14, Uppåkra 2:25, and Hjärup 9:8 resemble those of modern hazelnut shells experimentally charred for 2 h at less than 320°C (see Holguin et al., 2022) when viewed under a light microscope (Supplementary Figure 1). Unfortunately, none of the other charred nutshells were examined prior to analysis, but their external appearance was similar to those whose cross-sections were examined.

Modern hazelnut shells were freeze-dried and crushed using a ball mill. Archaeological nutshell fragments were crushed to a powder, either in an agate mortar and pestle or some of the Mesolithic nutshells with a ball mill. A subset of the crushed nutshells was analyzed using fourier transform infrared spectroscopy with attenuated total reflectance (FTIR-ATR). There are no peaks corresponding to carbonate, nitrate or humic acid contamination (Supplementary Datasheet 1). The charred archaeological nutshells were treated with 10 mL 0.5 M hydrochloric acid at 70°C for 30–60 min, then rinsed in distilled water three times before freeze-drying. Eight (out of forty) archaeological waterlogged hazelnut shells were treated with acid and their δ¹³C values compared with portions of the same powdered samples that had not been treated. There was no systematic change in the δ¹³C values with acid treatment and there was a maximum difference in the δ¹³C value between treated and untreated samples of 0.29‰ (mean = 0.18‰). This is a similar difference to that observed between duplicate carbon isotope measurements on the same homogenized sample of modern hazelnut shells (maximum difference 0.17‰; mean = 0.07‰). It was therefore decided not to proceed with further acid treatment of the waterlogged nutshells.

Approximately 1 mg of each homogenized sample was weighed into tin capsules for carbon isotopic analysis on a Sercon 20–22 isotope ratio mass spectrometer coupled to a Sercon GSL elemental analyser at the Research Laboratory for Archaeology and the History of Art, University of Oxford, UK. Stable carbon isotope values were calibrated to the VPDB scale using Caffeine-2^* (δ¹³C −35.05 ± 0.02‰; University of Indiana) and an internal seal bone collagen reference material (δ¹³C −12.54 ± 0.13‰). Measurement uncertainty was monitored using two reference materials: alanine (DL alanine, δ¹³C −27.18 ± 0.16‰; internal reference material from Sigma Aldrich) and spruce (δ¹³C −25.44 ± 0.02‰; Elemental Microanalysis). Precision [u(Rw)] was determined to be ±0.15‰, accuracy or systematic error [u(bias)] was ±0.19‰ and the total analytical uncertainty was estimated to be ±0.24‰. Raw and normalized isotope data of samples and reference materials are in Supplementary Table 3. The R script for calculating accuracy and precision using data in Supplementary Table 3 is in https://github.com/AmyStyring/Hazelnut-project.

The Δ¹³C values of modern hazelnut shells were calculated from the measured δ¹³C values (δ¹³C_plant) and an average δ¹³C value of atmospheric CO₂ (δ¹³C_air) determined from air sampled at weekly intervals during 2020 and 2021 (from Global Monitoring Laboratory https://gml.noaa.gov/dv/iadv/, Pallas-Sammaltunturi station, Finland), using the equation below from Farquhar et al. (1982).

The Δ¹³C values of archaeological hazelnut shells were calculated from the measured δ¹³C values (δ¹³C_plant) and a δ¹³C_air value approximated by the AIRCO2_LOESS system (Ferrio et al., 2005). The δ¹³C_air are in Supplementary Tables 1, 2. LAI values were split into three equal bins. Nuts collected from sub-sites with an LAI < 2.65 were assigned to the “open” category of canopy density; those with an LAI of 2.65–5.31 to the “semi-open” category; and those with an LAI > 5.31 to the “closed” category. These designations largely matched those that were subjectively assigned in the field (see Figure 2 for examples of these categories). The minimum LAI for the closed canopy category also aligns with the global mean LAI for temperate broadleaved forests (mean LAI = 5.1; Asner et al., 2003) and so could be considered to be the cut-off between understorey shrubs growing under woodland canopy and those either growing in the open or on the border of woodlands. Since nutshells charred at 320°C for 4 h have been found to have δ¹³C values that are on average 0.51 ± 0.20‰ higher than their uncharred counterparts (Holguin et al., in prep), we subtract 0.51‰ from the δ¹³C values of charred nutshells prior to plotting in figures and prior to any statistical analyses in order to make them comparable to the uncharred modern and waterlogged Mesolithic nutshells. There have been no studies of the effect of waterlogging on plant stable isotope values and so we make no adjustment to the δ¹³C values of the waterlogged hazelnut shells. Statistical analyses were performed in R v.4.2.2 and the R script used is in https://github.com/AmyStyring/Hazelnut-project.

Figure 3 shows variation in δ¹³C values within single hazelnut shells. The values vary by up to 0.6‰ within shells (mean = 0.5‰, n = 3). There is no consistent directional change in values from the base to apex of the hazelnut shells, although the base does have a lower δ¹³C value than the apex in all cases.

Figure 4 shows boxplots of the variation in δ¹³C values of hazelnut shells sampled from single trees. The nutshell values differ by up to 4.9‰ within trees (mean = 3.1‰, n = 9) and the standard deviations in δ¹³C values within a single tree vary between 0.5 and 2.2‰ (mean = 1.3‰, n = 9). Calculated from the standard error, the 95% confidence interval (CI) of intra-tree variability in hazelnut shell δ¹³C values varies between ± 0.4 and 1.9‰ (mean = 1.1‰, n = 9).

Figure 5 shows the relationship between hazelnut shell δ¹³C values and the measured LAI at each sub-site for hazelnuts collected in 2021 and 2022. The δ¹³C values (n = 192) are normally distributed (W = 0.99, p-value = 0.22). A linear mixed-effects model, including location (sub-site) as a random effect, finds a significant negative relationship between nutshell δ¹³C values and LAI (Beta = −0.47, SE = 0.07, t = −6.37, p < 0.0001). The model produced a reasonable fit, with a conditional RGLMM2 of 0.52 [calculated using r.squaredGLMM (mixed-effects model) in R]. Date was not found to have a significant effect on the relationship and the fit of the model was not improved by adding date as an additional coefficient. The fit of the model was also not improved by nesting location within site.

Figure 6 shows boxplots of the variation in δ¹³C values of hazelnut shells sampled from trees growing in three different densities of canopy: open, semi-open and closed (see Figure 2 for reference). These categories were assigned using the measured LAI values (see Section 2.4). The δ¹³C values are normally distributed within categories and Levene's test for equality of variance found that the variances were not statistically different among categories [F_{(2, 189)} = 0.32, p = 0.73]. Nested analysis of variance (with nutshell δ¹³C values nested by location) found a significant difference in the nutshell δ¹³C values among the three categories of canopy density (F = 19.93, p < 0.0001). The least squares mean of nutshell δ¹³C values growing in the open is −26.7‰ (95% CI −27.3 to −26.0‰); of those growing under semi-open canopy is −28.6‰ (95% CI −29.3 to −27.9‰); and of those growing under closed canopy is −29.8‰ (95% CI −30.6 to −29.0‰).

Although the nested analysis of variance finds a significant difference in the nutshell δ¹³C values among the three categories of canopy density, we wanted to test the reliability of this model to impute (predict) canopy density from measured δ¹³C values when the true category of canopy density is unknown (i.e., for archaeological samples). To do this, we split the modern hazelnut δ¹³C data into a training block and a testing block and fit the regression on the training data. We use the linear regression between δ¹³C value and category of canopy density with a random effect of location to fit the training data but set the random effect to zero in the prediction step as location identity will be unknown in archaeological data. The imputed canopy density is the category that minimizes the difference between the fitted and observed δ¹³C values. We split the data at random in the ratio 8:1 training to test and computed the prediction score. The prediction score is the proportion of correctly imputed canopy density levels. We repeated this 1,000 times with independent splits and report the average prediction score (Table 3). The average success rate for the model is 53%, which means that it performs better than random predictions, which would have a success rate of 1/3. The model correctly assigns nutshells to open and closed categories in 69 and 65% of cases, respectively, but performs poorly in correctly assigning nutshells to the semi-open category (only 19% correct assignments). The model performs the same with nutshell Δ¹³C values, which will be used when applying it to archaeological data.

Figure 7 shows boxplots of the Δ¹³C values of hazelnut shells from the four archaeological sites dated to the Mesolithic in southern Sweden. The archaeological sites are in chronologically ascending order (Slabälta 1 is the oldest and Rönneholm 6:1 and Rönneholm 10:3 are roughly contemporaneous). The measured δ¹³C values were converted to Δ¹³C values (see Section 2.4) to account for fluctuations in the δ¹³C value of atmospheric CO₂ over time. We carried out single imputation to assign a category of canopy density (open, semi-open, or closed) to each of the archaeological nutshell fragments. First, we used the modern hazelnut data to model the relationship between openness and Δ¹³C value (as in Section 3.3, but with Δ¹³C instead of δ¹³C value). We then inverted this relationship to impute the missing canopy density value in the archaeological data from their measured Δ¹³C value. The imputed category of canopy density is the category that minimizes the difference between the fitted and observed Δ¹³C values for each category. The individual datapoints for archaeological nutshells in Figure 7 are color-coded by their imputed canopy density category.

The nutshell Δ¹³C values differ by up to 5.8‰ within sites (mean = 4.0‰, n = 4). This is greater than the variation observed within single trees at the modern study sites. Calculated from the standard error, the 95% confidence interval (CI) of intra-site variability in hazelnut shell Δ¹³C values varies between ± 0.6 and 1.1‰ (mean = 0.8‰, n = 4). The Δ¹³C values are normally distributed within sites and Levene's test for equality of variance found that the variances were not statistically different among sites [F_{(3, 36)} = 1.00, p = 0.40]. Analysis of variance found a significant difference in the nutshell Δ¹³C values among the four sites [F_{(3, 36)} = 4.94, p < 0.05, ω = 0.48]. Post hoc comparisons using the Tukey HSD test indicated that the mean nutshell Δ¹³C value at Slabälta 1 (mean = 19.2‰, SD = 1.8) is significantly lower than that at Rönneholm 10:3 (mean = 21.2‰, SD = 1.3). While only one of the nutshells from Slabälta 1 is assigned to the closed canopy category, six (out of ten) of the nutshells from Rönneholm 10:3 are assigned to closed canopy.

Figure 8 shows boxplots of the Δ¹³C values of hazelnut shells from archaeological sites dated to the Mesolithic, Neolithic, Bronze Age and Iron Age in southern Sweden. We carried out single imputation to assign a category of canopy density (open, semi-open, or closed) to each of these archaeological nutshell fragments. The individual datapoints for archaeological nutshells in Figure 8 are color-coded by the imputed canopy density value. The standard deviations in nutshell Δ¹³C values within time periods are similar (SD = 1.4–1.5‰), apart from the Neolithic, which is only represented by two samples. Further statistical comparisons do not include the Neolithic samples because they are too few. The Δ¹³C values are not normally distributed within the Mesolithic but Levene's test for equality of variance found that the variances were not statistically different among periods [F_{(2, 58)} = 0.56, p = 0.58]. A Kruskal-Wallis test found a significant difference in the nutshell Δ¹³C values among the three periods [H₍₂₎ = 7.32, p = 0.026]. Jonckheere's test revealed a significant trend in the data: a decrease in median nutshell Δ¹³C value between the Mesolithic and Iron Age, J = 284, p = 0.014. Comparisons of the mean ranks between periods showed that nutshell Δ¹³C values are not significantly different between the Mesolithic and Bronze Age (difference = 3.3) or between the Bronze Age and Iron Age (difference = 16.4), but there is a significant difference between the Mesolithic and Iron Age nutshell Δ¹³C values (difference = 13.2). While 12 (out of 40, or 30%) of the Mesolithic nutshells are assigned to the closed canopy category, only two (out of 18) dated to the Iron Age are assigned to closed canopy.

The mean 0.5‰ range in δ¹³C values within a single hazelnut shell (Figure 3) is small compared to the mean 3.1‰ variation in nutshell δ¹³C values from a single tree (Figure 4) and the total range in hazelnut nutshell δ¹³C values determined across three sites in southern Sweden (11.2‰; Figure 5). It is therefore unlikely that sampling fragments from random parts of a nutshell will bias δ¹³C results.

The mean 3.1‰ intra-tree range in nutshell δ¹³C values (Figure 4) is relatively large and is actually the same as the mean intra-tree range of Corylus avellana pollen δ¹³C values determined in another study (3.1‰, n = 4; Müller et al., 2020). This suggests that the carbon used to synthesize the lignin and structural carbohydrates (cellulose, hemicellulose) that make up the majority of nutshell (Licursi et al., 2017) is either (i) sourced from different (groups of) leaves that experience different conditions, which influence their stomatal conductance and photosynthetic activity and therefore the level of discrimination against ¹³C; or (ii) different nuts are formed at slightly different times when conditions that affect ¹³C-discrimination differ. Given that nuts tend to develop at a similar time within the year (Hagerup, 1942), it seems most likely that carbon being derived from different leaves accounts for the variation in nutshell δ¹³C values. A study of leaves sampled at different heights (150 leaves were homogenized for each height) within a single beech tree (Fagus sylvatica) found a positive correlation between leaf δ¹³C values and height. The intra-tree leaf δ¹³C values had a range of c. 3‰ (Schleser, 1992: Figure 1), but this range is likely to have been higher if the δ¹³C values of individual leaves had been determined. Hazelnut shells must be a sink for carbon from neighboring leaves as opposed to receiving carbon from a mix of sources, which would lead to much more homogenous hazelnut δ¹³C values from within a single tree.

In this study we were unable to pick nuts from designated parts of the tree canopy to explore the relationship between nutshell δ¹³C values and position on the branch, height above the ground and irradiance received by the surrounding leaves, but this could be investigated in the future. Our results indicate that single nutshell fragment δ¹³C values are only 95% likely to encompass the true tree mean if their confidence range is extended by ±1.1‰ (mean 95% CI of intra-tree variability). Differences in nutshell δ¹³C values of less than 2.2‰ therefore do not necessarily reflect differences in growing condition and it is recommended that multiple nutshells are sampled from the same context/site/period to adequately characterize the growing conditions. Indeed, to identify a difference in mean nutshell δ¹³C values between two sites/phases of 3.1‰, which is the difference between δ¹³C values of modern hazelnuts growing in the open and under closed canopy (Section 3.3), given that the pooled standard deviation of nutshells from a single tree is 1.3‰, four samples from each site/phase are required. To identify a difference in mean nutshell δ¹³C values between two sites/phases of 1.2‰, which is the difference between δ¹³C values of modern hazelnuts growing under closed and semi-open canopy (Section 3.3), 20 samples from each site/phase are required. This is with a statistical power of 80% and a significance level of α = 0.05.

There is a significant negative correlation between LAI (as a proxy for canopy density) and modern hazelnut shell δ¹³C values (Figure 5), indicating that hazelnut shell δ¹³C values can be used to infer the density of the canopy in the environment in which they were collected. The conditional R² value of 0.52 indicates that 52% of the variation in nutshell δ¹³C values can be accounted for by canopy density as opposed to other variables. There were slightly lower levels of precipitation and higher levels of solar radiation in 2022 compared to 2021 (Table 1), which could have led to higher nutshell δ¹³C values due to lower water availability, but these meteorological differences (categorized by the date variable) were not found to have a significant effect on the fit of the linear model regressing nutshell δ¹³C values and LAI. It cannot be ruled out that more extreme changes in precipitation and solar radiation levels could have a significant effect on nutshell δ¹³C values and therefore we can only be confident of this relationship for regions with similar climate to modern southern Sweden.

Studies on tree leaf δ¹³C values have found that they vary with height in the canopy. This variation has been attributed predominantly to light intensity, with leaves in the upper canopy that receive more light having higher δ¹³C values (Berry et al., 1997; Buchmann et al., 1997; Roberts et al., 1999; Graham et al., 2014). A study by Bonafini et al. (2013) investigated the variation in grass δ¹³C values growing in differing levels of shade. They found that differences in light intensity had a stronger influence on the δ¹³C value of the grass tissue than temperature or rainfall. It therefore seems likely that the variation in hazelnut shell δ¹³C values with LAI are due to differences in the δ¹³C value of carbon fixed by leaves that discriminate against ¹³C to varying degrees according to the rate of photosynthesis that is influenced by light levels.

Studies of plant δ¹³C values growing in differing densities of canopy have determined a depletion in ¹³C of between 2 and 5‰ between plants growing in open and closed environments (van der Merwe and Medina, 1991; Broadmeadow and Griffiths, 1993; France, 1996; Bonafini et al., 2013). In this study the difference between the mean δ¹³C values of hazelnut shells growing in the open and those growing under closed canopy is 3.1‰ (−26.7 minus−29.8‰), and between the highest and lowest 95% CI values of hazelnut shells growing in the open and those growing under closed canopy is 4.6‰ (−26.0 minus−30.6‰). Our values are therefore consistent with the results of previous studies.

We have grouped the measured LAI values into three bins, which correspond to open, semi-open and closed canopy environments. These categories are useful because they represent different environments that can be more easily visualized than a gradient of LAI values. Figure 2 shows examples of the three categories. The “open” category can broadly be thought of as trees growing in wood-pastures, with open space around each tree. This is how hazels are cultivated today (e.g., Portarena et al., 2022). The “semi-open” category could encompass hazels growing on the margins of woodlands, with open land on one side, or those growing in woodland clearings, with no canopy directly above them. The “closed” category represents hazels growing as understorey trees/shrubs in woodland. These three settings have different implications for the types of environments in which hazels were growing and in which humans were collecting them. The fact that the δ¹³C values of hazelnut shells from these three different settings differ significantly means that hazelnut δ¹³C values are a promising index to differentiate between them.

We tested the reliability of the nested analysis of variance model to impute canopy density from measured δ¹³C values (Table 3). The model correctly categorizes hazelnut shells into open and closed categories in over 65% of cases but does not perform well in correctly categorizing nutshells into the semi-open category. It is therefore recommended that studies that seek to characterize past hazel growing environments using nutshell δ¹³C values focus on the relative proportions of nutshells from open and closed environments and refrain from interpreting those that are categorized as having grown in semi-open conditions. Nonetheless, general trends toward higher or lower δ¹³C values can be interpreted as a shift toward more open or more closed environments because the relationship between nutshell δ¹³C values and LAI/canopy density is strong and significant.

The Δ¹³C values of hazelnut shells from Mesolithic archaeological sites within 5 km of one another in southern Sweden are found to vary significantly (Figure 7). There is an even split between nutshells from open and closed environments, suggesting that people had access to and exploited a range of wooded and more open environments for their food resources. This aligns with pollen analysis that indicates that southern Sweden was dominated by natural broad-leaved forests with smaller areas of open land throughout the Mesolithic (O'Dwyer et al., 2021). It isn't clear why the Δ¹³C values of hazelnut shells from Slabälta 1 and Rönneholm 10:3 should differ significantly from one another. Slabälta 1 (7000–6640 BCE) was occupied at least one thousand years before the Rönneholm sites (5728–5462 BCE) and regional pollen-based land cover reconstructions of southern Sweden show that the pine and hazel-dominated vegetation of the early Mesolithic (c. 8500–7500 BCE) was replaced by broad-leaved forest with more abundant Quercus and Tilia from around 7000 BCE (Gaillard, 1984; Berglund et al., 2008; O'Dwyer et al., 2021). The lower Δ¹³C values of nutshells from Slabälta 1 could therefore possibly be linked to more open vegetation earlier in the Mesolithic.

The decreasing Δ¹³C values of hazelnut shells from archaeological sites in southern Sweden dated to between the Mesolithic and Iron Age (c. 7000 BCE to 1000 CE) are consistent with hazelnuts being gathered from increasingly open environments (Figure 8). It is possible that this trend could be due to decreasing water availability over this period, leading to decreased stomatal conductance and decreased discrimination against ¹³C. This is contrary to what is known about the past climate in Sweden, however, with a thermal maximum between c. 6000 and 3000 BCE and a shift to colder and wetter conditions from c. 2000 BCE (Wastegård, 2022), making it unlikely that water availability is the main factor influencing the nutshell δ¹³C values. Until further studies comprehensively assess the effect of water availability on hazelnut shell δ¹³C values in this region, however, we cannot rule out the possibility that changes in relative humidity through time had some effect on the archaeological hazelnut shell δ¹³C values.

If canopy density was the main factor driving the changes in hazelnut shell δ¹³C values, it seems that by the Iron Age, the vast majority of nuts were being gathered from outside woodlands, perhaps from open wood-pastures. This is consistent with pollen records, which indicate an opening up of the landscape between the Mesolithic and Iron Age, with tree cover decreasing from about 80% in the middle Mesolithic to 30% by the pre-Roman Iron Age (Berglund, 1991; Hellman, 2007; O'Dwyer et al., 2021). The degree of openness did not increase uniformly across southern Sweden, however, with more open land along the coasts in the Bronze Age (Berglund, 2003) and local catchment pollen cores (like wells) showing high variability in tree cover (Lagerås and Fredh, 2020). Pollen records indicate a sudden shift toward a more open landscape in the middle of the Bronze Age (Berglund, 1991; Hellman, 2007; O'Dwyer et al., 2021), and so increasing the number of nutshell carbon isotope measurements from the Bronze Age should be a priority to explore whether there was a shift in the environment in which hazelnuts were gathered at this time. The number of nutshells from Neolithic sites should also be increased. Overall, the Δ¹³C values of archaeological nutshell fragments provide a new insight into the environments in which hazelnuts were gathered in the past that complement pollen-based reconstructions of overall vegetation cover.