Front. Environ. Archaeol. Frontiers in Environmental Archaeology Front. Environ. Archaeol. 2813-432X Frontiers Media S.A. 10.3389/fearc.2024.1351411 Environmental Archaeology Original Research Carbon isotope values of hazelnut shells: a new proxy for canopy density Styring Amy K. 1 * Jirdén Elin 2 Lagerås Per 3 Larsson Mikael 4 Sjöström Arne 4 Ljung Karl 2 * 1School of Archaeology, University of Oxford, Oxford, United Kingdom 2Department of Geology, Lund University, Lund, Sweden 3The Archaeologists, National Historical Museums, Lund, Sweden 4Department of Archaeology and Ancient History, Lund University, Lund, Sweden

Edited by: Xinyi Liu, Washington University in St. Louis, United States

Reviewed by: Yu Dong, Shandong University, China

Nicasio T. Jiménez-Morillo, University of Evora, Portugal

*Correspondence: Karl Ljung karl.ljung@geol.lu.se Amy K. Styring amy.styring@arch.ox.ac.uk
29 02 2024 2024 3 1351411 06 12 2023 09 01 2024 Copyright © 2024 Styring, Jirdén, Lagerås, Larsson, Sjöström and Ljung. 2024 Styring, Jirdén, Lagerås, Larsson, Sjöström and Ljung

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Hazel (Corylus avellana) has been abundant in the vegetation of northern and central Europe since the early Holocene and has provided food and materials for humans ever since. Here we use stable carbon isotope (δ13C) values of hazelnut shells to infer woodland openness based on the premise of the “canopy effect”. It is well established that plants growing in dense, shaded forests have lower carbon isotope (δ13C) values than plants growing in open areas. By measuring δ13C values in hazelnuts collected from trees growing in different levels of light intensity, we show that the canopy effect is preserved in hazelnuts and that their δ13C values can be used to infer woodland openness in the past. We apply the method to hazelnuts recovered from sites dated to between the Mesolithic and Iron Age (c. 7000 BCE−1000 CE) in southern Sweden. Our results show that the nuts dated to the Mesolithic were harvested from hazels growing in a range of closed to open settings while nuts from subsequent periods were harvested from progressively more open environments. Given the abundance of hazelnuts recovered from many archaeological contexts, this method has the potential to reconstruct the microhabitats exploited by humans in the past and explore the impact of humans on their environment.

archaeology environmental reconstruction δ13C values Mesolithic Sweden section-at-acceptance Archaeological Isotope Analysis

香京julia种子在线播放

    1. <form id=HxFbUHhlv><nobr id=HxFbUHhlv></nobr></form>
      <address id=HxFbUHhlv><nobr id=HxFbUHhlv><nobr id=HxFbUHhlv></nobr></nobr></address>

      1 Introduction

      Palynological studies show that hazel (Corylus avellana L.) has been continuously present in northern and central Europe since the early Holocene (e.g., Huntley and Birks, 1983; Giesecke et al., 2011). Moreover, hazelnut shells are frequently recovered from archaeological sites, particularly those dated to the Mesolithic and Neolithic (Zvelebil, 1994; Perry, 1999; Regnell, 2012; Bishop et al., 2014), demonstrating their important role in human subsistence. Although their actual importance in past diets is still debated because nutshells are more resistant to charring and other taphonomic processes compared to other plants (Jones, 2000; Bishop, 2019), they have the potential to provide an important source of additional information about the past environments in which hazel trees grew and which humans exploited. Hazelnut shells are also a preferred sample for radiocarbon dating because their carbon was fixed within a single year (Bayliss and Marshall, 2022, p. 36), meaning that any environmental information they do yield can be directly dated.

      Stable carbon isotope (δ13C) values of plants are influenced by the ratio of leaf intercellular (ci) to ambient (ca) carbon dioxide concentrations (Farquhar et al., 1982). In C3 plants, this ratio is strongly affected by stomatal conductance and photosynthetic activity, which are in turn affected by water availability and light intensity (Farquhar et al., 1989). When soil moisture levels decrease, stomatal conductance decreases, decreasing the ci and resulting in less negative δ13C values. Numerous studies have observed a negative relationship between C3 plant δ13C values and mean annual precipitation, which flattens off above about 1,000 mm per year (summarized in Kohn, 2010). Light intensity is an additional factor that has been found to influence plant δ13C values, with plants growing in lower light levels having lower rates of photosynthesis, leading to higher ci that results in lower δ13C values (Ehleringer et al., 1986, 1987). Lower δ13C values of plants growing under dense canopy have also been attributed to the recycling of 13C-depleted carbon dioxide, produced by decomposing leaf litter on the forest floor (Vogel, 1978). The observation that plants growing in denser forests have lower δ13C values has been dubbed the “canopy effect” and has been invoked in archaeological studies to attribute variation in wild herbivore bone collagen δ13C values to variability in the density of forest cover (e.g., Bocherens et al., 1999; Cerling and Harris, 1999; Drucker et al., 2008). The δ13C values of archaeological plant remains are rarely considered as a proxy for canopy density, however.

      Hazel grows as a multi-stemmed bush in the understorey of deciduous and coniferous forests, but also as isolated trees in wood-pastures. It is generally found in temperate climates with moderate/high annual rainfall (800–1,000 mm; Fideghelli and De Salvador, 2009). There have been numerous studies on the effect of irrigation, fertilizer, and genotype on the characteristics (e.g., protein content, fatty acid profile, morphology) and yield of commercially grown hazelnuts (e.g., Girona et al., 1994; Yao and Mehlenbacher, 2000; Akçin and Bostan, 2018; Külahçilar et al., 2018), but most of these have been in the Mediterranean region where the majority of hazelnuts are now grown commercially. Studies have found that hazel nut yields are sensitive to water stress (Girona et al., 1994; Mingeau et al., 1994; Tombesi and Rosati, 1997), and that the nuts are a strong sink for water during ripening (Bregaglio et al., 2016). It has been posited that their sensitivity to water stress is due to their low capacity for stomatal control and relatively shallow root system (0.4–0.5 m in depth; Portarena et al., 2022). These studies, however, have all been carried out in southern Europe, where rainfall levels are relatively low and water availability is a limiting factor in plant growth. Indeed, a global review of the δ13C values of the leaves and wood tissues of conifers, which also have shallow root systems, concluded that water availability only has a significant effect in seasonally dry climates, where evaporation is greater than precipitation (Warren et al., 2001).

      Corylus avellana has been found to be highly adaptable to different light conditions (Catoni et al., 2015), allowing it to grow in both sunny and shaded environments. It seems possible, therefore, that in northern Europe, hazelnut shell δ13C values will be most closely correlated with light levels during the time of fruit formation. According to a study in Denmark, closest to our study region, the fruit (which includes the kernel and nutshell) starts to form in May but does not begin to grow appreciably in size until fertilization takes place at the beginning of July. By the end of July, the shell is fully developed and begins to lignify at the pointed apex. Once the shell is formed, the nut begins to grow quickly and reaches full size in about one month. The fruit then remains on the tree for another month with no further external changes occurring, and during September it will detach itself from the involucre and is ready to harvest. The nut grows from the pointed apex toward the rounded base, known as the hilum (Hagerup, 1942). It can therefore be posited that the carbon routed to the growing nutshell will reflect growing conditions in July.

      In this study, we investigate whether the δ13C values of hazelnut shells vary with light levels, governed by canopy density. We collected hazelnuts from trees growing in varying light levels at three locations in southern Sweden (Figure 1). We determine the intra-shell variability in δ13C values, intra-tree variability in δ13C values and the relationship between nutshell δ13C values and light level, measured as the leaf area index (LAI). We collected hazelnuts from the same locations in two subsequent years to preliminarily explore a potential effect of water availability (due to differences in precipitation amount and relative humidity between the years; Table 1) on hazelnut shell δ13C values. This was not the focus of this study, however, and future work needs to be carried out to more thoroughly assess the relative effects of light level and water availability on the δ13C values of hazelnuts in this region. We then determine the δ13C values of hazelnut shells from archaeological sites in southern Sweden, dated to between the Mesolithic and Iron Age (c. 7000 BCE−1000 CE), as a pilot study to gain an initial insight into how the environment—in terms of canopy density—in which people collected hazelnuts changed through time.

      Map of southern Sweden showing the location of the modern hazelnut sampling sites, archaeological hazelnut sampling sites, and meteorological stations providing precipitation, relative humidity, cloud cover, hours of sunshine, and solar radiation data. Sites: 1. Dörröd, 2. Holma, 3. Linnebjer, 4. Fiskeby, 5. Ingelstad 6:1, 6. Ystad VA schaktöverv, 7. Uppåkra regional center, 8. Uppåkra 2:14, 9. Uppråkra 2:25, 10. Hjärup 9:8, 11. Hjärup 7:1, 22:1, 12. Flädie Tankställe Damm, 13. Borstahusen SU 2020, 14. Stångby 5:28, 15. Slabälta 1, 16. Ringsjöholm, 17. Rönneholm 10:3, 18. Rönneholm 6:1, 19. Vomb, 20. Hörby A, 21. Lund Sol.

      Meteorological observations for 2021 and 2022, when hazelnuts were sampled in this study.

      Parameter Units Meteorological station Latitude Longitude 2021 2022
      Jan–Aug May–Aug Jan–Aug May–Aug
      Total precipitation mm Vomb 55.66 13.53 447 259.5 417.9 202.2
      Average relative humidity % Hörby A 55.86 13.67 79 76 75 73
      Average cloud cover % Hörby A 55.86 13.67 66 62 55 54
      Total hours of sunshine hours Lund Sol 55.71 13.21 1485 937 1690 1018
      Total solar radiation kWh/m2 Lund Sol 55.71 13.21 880 631 955 665
      2 Materials and methods 2.1 Field sampling methods

      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).

      2.2 Sampling of archaeological hazelnut shells

      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.

      Information on modern and archaeological sites from which hazelnuts were sampled for this study.

      No. Site Latitude Longitude Site type Archaeological period Date range
      1 Dörröd 55.60 13.46 Modern sampling site 2021–2022
      2 Holma 55.95 13.55 Modern sampling site 2021–2022
      3 Linnebjer 55.73 13.30 Modern sampling site 2021–2022
      4 Fiskeby 58.6 16.11 Archaeological residential Iron Age 500 BC–AD 400
      Bronze Age 1100–1 BC
      5 Ingelstad 6:1 56.75 14.92 Archaeological residential Iron Age 500–1 BC
      6 Ystad VA schaktöverv 55.43 13.82 Archaeological residential Iron Age AD 800–1050
      7 Uppåkra regional center 55.66 13.16 Archaeological residential Iron Age 45 BC–AD 600
      8 Uppåkra 2:14 55.67 13.16 Archaeological residential Iron Age AD 0-250
      9 Uppåkra 2:25 55.67 13.17 Archaeological residential Iron Age AD 100-400
      10 Hjärup 9:8 55.66 13.13 Archaeological residential Iron Age AD 900–1000
      11 Hjärup 7:1, 22:1 55.67 13.14 Archaeological residential Iron Age AD 1–1050
      12 Flädie Tankställe Damm 55.73 13.07 Archaeological residential Bronze Age 1700–1100 BC
      Neolithic 4000–2400 BC
      13 Borstahusen SU 2020 55.9 12.81 Archaeological residential Bronze Age 1100–520 BC
      14 Stångby 5:28 55.75 13.2 Archaeological residential Neolithic 2400–1700 BC
      15 Slabälta 1 55.94 13.44 Archaeological short-lived campsite Mesolithic 7040–6600 BC
      16 Ringsjöholm 55.9 13.43 Archaeological residential Mesolithic 7070–5840 BC
      17 Rönneholm 10:3 55.93 13.42 Archaeological residential Mesolithic 6020–5660 BC
      18 Rönneholm 6:1 55.93 13.42 Archaeological residential Mesolithic 5728–5462 BC

      More details can be found in Supplementary Table 2.

      2.3 Laboratory analysis

      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 δ13C values compared with portions of the same powdered samples that had not been treated. There was no systematic change in the δ13C values with acid treatment and there was a maximum difference in the δ13C 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*13C −35.05 ± 0.02‰; University of Indiana) and an internal seal bone collagen reference material (δ13C −12.54 ± 0.13‰). Measurement uncertainty was monitored using two reference materials: alanine (DL alanine, δ13C −27.18 ± 0.16‰; internal reference material from Sigma Aldrich) and spruce (δ13C −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.

      2.4 Data analysis

      The Δ13C values of modern hazelnut shells were calculated from the measured δ13C values (δ13Cplant) and an average δ13C value of atmospheric CO213Cair) 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).

      Δ13C = δ13Cairδ13Cplant1+ δ13Cplant/1000

      The Δ13C values of archaeological hazelnut shells were calculated from the measured δ13C values (δ13Cplant) and a δ13Cair value approximated by the AIRCO2_LOESS system (Ferrio et al., 2005). The δ13Cair 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 δ13C 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 δ13C 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 δ13C 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.

      Examples of hazels growing in (A) open, (B) semi-open, and (C) closed canopy environments.

      3 Results 3.1 Variability in δ<sup>13</sup>C values within single hazelnut shells

      Figure 3 shows variation in δ13C 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 δ13C value than the apex in all cases.

      Intra-nutshell variation in δ13C values.

      3.2 Variability in hazelnut shell δ<sup>13</sup>C values within single trees

      Figure 4 shows boxplots of the variation in δ13C 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 δ13C 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 δ13C values varies between ± 0.4 and 1.9‰ (mean = 1.1‰, n = 9).

      Intra-tree variation in modern nutshell δ13C values. Boxes represent the quartiles, the bold line represents the median and whiskers represent 1.5 × the interquartile range. Trees are color-coded by the density of canopy (open, semi-open, closed) in which they were growing. These categories were assigned by binning the measured LAI values.

      3.3 The effect of canopy density on hazelnut shell δ<sup>13</sup>C values

      Figure 5 shows the relationship between hazelnut shell δ13C values and the measured LAI at each sub-site for hazelnuts collected in 2021 and 2022. The δ13C 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 δ13C 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.

      Modern nutshell δ13C values plotted against leaf area index (LAI). Symbols differ by sampling year and are color-coded by site. The dotted line represents a fitted linear model relating δ13C values and leaf area index. The gray shading represents the 95% confidence interval of this relationship.

      Figure 6 shows boxplots of the variation in δ13C 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 δ13C 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 δ13C values nested by location) found a significant difference in the nutshell δ13C values among the three categories of canopy density (F = 19.93, p < 0.0001). The least squares mean of nutshell δ13C 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‰).

      Comparison of modern nutshell δ13C values among trees categorized as growing in open, semi-open and closed environments. Boxes represent the quartiles, the bold line represents the median and whiskers represent 1.5 × the interquartile range.

      Although the nested analysis of variance finds a significant difference in the nutshell δ13C values among the three categories of canopy density, we wanted to test the reliability of this model to impute (predict) canopy density from measured δ13C values when the true category of canopy density is unknown (i.e., for archaeological samples). To do this, we split the modern hazelnut δ13C data into a training block and a testing block and fit the regression on the training data. We use the linear regression between δ13C 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 δ13C 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 Δ13C values, which will be used when applying it to archaeological data.

      Confusion matrix for prediction of canopy density category using modern hazelnut δ13C values split into a training and test dataset.

      Predicted
      Open Semi-open Closed
      Open 69 21 10
      True Semi-open 35 19 46
      Closed 15 21 65

      Numbers are percentages and sum to 100 across each row (with exceptions due to rounding). Numbers in bold are the percentage of hazelnut shells assigned to the correct canopy density category.

      3.4 Change in archaeological hazelnut shell δ<sup>13</sup>C values between the Mesolithic and Iron age in southern Sweden

      Figure 7 shows boxplots of the Δ13C 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 δ13C values were converted to Δ13C values (see Section 2.4) to account for fluctuations in the δ13C value of atmospheric CO2 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 Δ13C value (as in Section 3.3, but with Δ13C instead of δ13C value). We then inverted this relationship to impute the missing canopy density value in the archaeological data from their measured Δ13C value. The imputed category of canopy density is the category that minimizes the difference between the fitted and observed Δ13C values for each category. The individual datapoints for archaeological nutshells in Figure 7 are color-coded by their imputed canopy density category.

      Comparison of modern and Mesolithic nutshell Δ13C values from different growing environments (modern) and sites (archaeological). Boxes represent the quartiles, the bold line represents the median and whiskers represent 1.5 × the interquartile range. Datapoints are color-coded by the canopy density category. For Mesolithic nutshells this category was interpolated from the measured Δ13C values.

      The nutshell Δ13C 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 Δ13C values varies between ± 0.6 and 1.1‰ (mean = 0.8‰, n = 4). The Δ13C 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 Δ13C 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 Δ13C 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 Δ13C 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 Δ13C 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 Δ13C 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 Δ13C values among the three periods [H(2) = 7.32, p = 0.026]. Jonckheere's test revealed a significant trend in the data: a decrease in median nutshell Δ13C value between the Mesolithic and Iron Age, J = 284, p = 0.014. Comparisons of the mean ranks between periods showed that nutshell Δ13C 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 Δ13C 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.

      Comparison of modern and archaeological nutshell Δ13C values from different growing environments (modern) and time periods (archaeological). Boxes represent the quartiles, the bold line represents the median and whiskers represent 1.5 × the interquartile range. Datapoints are color-coded by the canopy density category. For archaeological nutshells this category was interpolated from the measured Δ13C values.

      4 Discussion 4.1 Variability in hazelnut shell δ<sup>13</sup>C values within single nuts and trees

      The mean 0.5‰ range in δ13C values within a single hazelnut shell (Figure 3) is small compared to the mean 3.1‰ variation in nutshell δ13C values from a single tree (Figure 4) and the total range in hazelnut nutshell δ13C 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 δ13C results.

      The mean 3.1‰ intra-tree range in nutshell δ13C values (Figure 4) is relatively large and is actually the same as the mean intra-tree range of Corylus avellana pollen δ13C 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 13C; or (ii) different nuts are formed at slightly different times when conditions that affect 13C-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 δ13C 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 δ13C values and height. The intra-tree leaf δ13C values had a range of c. 3‰ (Schleser, 1992: Figure 1), but this range is likely to have been higher if the δ13C 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 δ13C 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 δ13C 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 δ13C 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 δ13C 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 δ13C values between two sites/phases of 3.1‰, which is the difference between δ13C 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 δ13C values between two sites/phases of 1.2‰, which is the difference between δ13C 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.

      4.2 Hazelnut shell δ<sup>13</sup>C values as a proxy for canopy density?

      There is a significant negative correlation between LAI (as a proxy for canopy density) and modern hazelnut shell δ13C values (Figure 5), indicating that hazelnut shell δ13C values can be used to infer the density of the canopy in the environment in which they were collected. The conditional R2 value of 0.52 indicates that 52% of the variation in nutshell δ13C 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 δ13C 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 δ13C 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 δ13C values and therefore we can only be confident of this relationship for regions with similar climate to modern southern Sweden.

      Studies on tree leaf δ13C 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 δ13C 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 δ13C values growing in differing levels of shade. They found that differences in light intensity had a stronger influence on the δ13C value of the grass tissue than temperature or rainfall. It therefore seems likely that the variation in hazelnut shell δ13C values with LAI are due to differences in the δ13C value of carbon fixed by leaves that discriminate against 13C to varying degrees according to the rate of photosynthesis that is influenced by light levels.

      Studies of plant δ13C values growing in differing densities of canopy have determined a depletion in 13C 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 δ13C 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 δ13C values of hazelnut shells from these three different settings differ significantly means that hazelnut δ13C 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 δ13C 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 δ13C 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 δ13C values can be interpreted as a shift toward more open or more closed environments because the relationship between nutshell δ13C values and LAI/canopy density is strong and significant.

      4.3 Reconstructing changes in hazel growing environments through time

      The Δ13C 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 Δ13C 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 Δ13C values of nutshells from Slabälta 1 could therefore possibly be linked to more open vegetation earlier in the Mesolithic.

      The decreasing Δ13C 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 13C. 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 δ13C values. Until further studies comprehensively assess the effect of water availability on hazelnut shell δ13C 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 δ13C values.

      If canopy density was the main factor driving the changes in hazelnut shell δ13C 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 Δ13C 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.

      5 Conclusions

      The δ13C values of hazelnut shells have been found to vary with canopy density. Stable carbon isotope analysis of hazelnut shells from archaeological sites can therefore be used to provide direct insights into the nature of wooded environments in which people were collecting food resources in the past. This can complement the more general understanding of vegetation cover and density provided by pollen analysis and δ13C values of wild faunal bone collagen that integrate vegetation over a wider area. On Mesolithic and Neolithic sites, where hazelnut shells are abundant, nutshells from discrete palimpsests can help to characterize the microhabitats that were exploited by different groups. In later periods, nutshells with δ13C values indicative of growing in closed canopy environments can identify sites that were in the proximity of more closed woodlands.

      Data availability statement

      The original contributions presented in the study are included in the article/Supplementary material and the R script used for statistical analyses is in https://github.com/AmyStyring/Hazelnut-project. Further inquiries can be directed to the corresponding authors.

      Author contributions

      AKS: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Supervision, Visualization, Writing—original draft, Writing—review & editing. EJ: Data curation, Investigation, Writing—review & editing. PL: Resources, Writing—review & editing. ML: Resources, Writing—review & editing. AS: Resources, Writing—review & editing. KL: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing—original draft, Writing—review & editing.

      Funding

      The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study was funded by the Ebbe Kocks Stiftelse (Grant No. EK202t-0021).

      We are grateful to Amy Holguin for sharing unpublished isotopic data from the experimental charring of hazelnuts and to Peter Ditchfield for technical assistance on the EA-IRMS. We also thank the two reviewers for their insightful comments that have improved the manuscript.

      Conflict of interest

      The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

      The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

      Publisher's note

      All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

      Supplementary material

      The Supplementary Material for this article can be found online at: /articles/10.3389/fearc.2024.1351411/full#supplementary-material

      References Akçin Y. Bostan S. Z. (2018). Variations in some traits affecting yield in different hazelnut cultivars. Acta Hortic. 1226, 149152. 10.17660/ActaHortic.2018.1226.21 Asner G. P. Scurlock J. M. O. Hicke A. J. (2003). Global synthesis of leaf area index observations: implications for ecological and remote sensing studies. Glob. Ecol. Biogeogr. 12, 191205. 10.1046/j.1466-822X.2003.00026.x Bayliss A. Marshall P. (2022). Radiocarbon Dating and Chronological Modelling: Guidelines and Best Practice. London: Historic England. Berglund B. E. (1991). Landscape reconstruction in south Sweden for the past 6000 years. Proc. Br. Acad. 77, 2537. Berglund B. E. (2003). Human impact and climate changes—synchronous events and a causal link? Quat. Int. Environ. Variabil. East West Eur. 105, 712. 10.1016/S1040-6182(02)00144-1 Berglund B. E. Gaillard M.-J. Björkman L. Persson T. (2008). Long-term changes in floristic diversity in southern Sweden: palynological richness, vegetation dynamics and land-use. Veget. Hist. Archaeobot. 17, 573583. 10.1007/s00334-007-0094-x Berry S. C. Varney G. T. Flanagan L. B. (1997). Leaf δ13C in Pinus resinosa trees and understory plants: variation associated with light and CO2 gradients. Oecologia 109, 499506. 10.1007/s00442005011028307333 Bishop R. R. (2019). Experiments on the effects of charring on hazelnuts and their representation in the archaeological record. J. Archaeol. Sci. Rep. 26:101839. 10.1016/j.jasrep.2019.05.004 Bishop R. R. Church M. J. Rowley-Conwy P. A. (2014). Seeds, fruits and nuts in the Scottish Mesolithic. Proc. Soc. Antiq. Scot. 143, 971. 10.9750/PSAS.143.9.71 Bocherens H. Billiou D. Mariotti A. Patou-Mathis M. Otte M. Bonjean D. . (1999). Palaeoenvironmental and palaeodietary implications of isotopic biogeochemistry of last interglacial neanderthal and mammal bones in Scladina Cave (Belgium). J. Archaeol. Sci. 26, 599607. 10.1006/jasc.1998.0377 Bonafini M. Pellegrini M. Ditchfield P. Pollard A. M. (2013). Investigation of the ‘canopy effect' in the isotope ecology of temperate woodlands. J. Archaeol. Sci. 40, 39263935. 10.1016/j.jas.2013.03.028 Bregaglio S. Orlando F. Forni E. De Gregorio T. Falzoi S. Boni C. . (2016). Development and evaluation of new modelling solutions to simulate hazelnut (Corylus avellana L.) growth and development. Ecol. Modell. 329, 8699. 10.1016/j.ecolmodel.2016.03.006 Broadmeadow M. S. J. Griffiths H. (1993). “Carbon isotope discrimination and the coupling of CO2 fluxes within forest canopies,” in Stable Isotopes and Plant Carbon-Water Relations, eds J. R. Ehleringer, A. E. Hall, and G. D. Farquhar (San Diego, CA: Academic Press), 109129. Buchmann N. Kao W.-Y. Ehleringer J. (1997). Influence of stand structure on carbon-13 of vegetation, soils, and canopy air within deciduous and evergreen forests in Utah, United States. Oecologia 110, 109119. 10.1007/s00442005013928307459 Catoni R. Granata M. U. Sartori F. Varone L. Gratani L. (2015). Corylus avellana responsiveness to light variations: morphological, anatomical, and physiological leaf trait plasticity. Photosynthetica 53, 3546. 10.1007/s11099-015-0078-5 Cerling T. E. Harris J. M. (1999). Carbon isotope fractionation between diet and bioapatite in ungulate mammals and implications for ecological and paleoecological studies. Oecologia 120, 347363. 10.1007/s00442005086828308012 Drucker D. G. Bridault A. Hobson K. A. Szuma E. Bocherens H. (2008). Can carbon-13 in large herbivores reflect the canopy effect in temperate and boreal ecosystems? Evidence from modern and ancient ungulates. Palaeogeogr. Palaeoclimatol. Palaeoecol. 266, 6982. 10.1016/j.palaeo.2008.03.020 Ehleringer J. R. Field C. B. Lin Z. Kuo C. (1986). Leaf carbon isotope and mineral composition in subtropical plants along an irradiance cline. Oecologia 70, 520526. 10.1007/BF0037989828311493 Ehleringer J. R. Lin Z. F. Field C. B. Sun G. C. Kuo C. Y. (1987). Leaf carbon isotope ratios of plants from a subtropical monsoon forest. Oecologia 72, 109114. 10.1007/BF0038505328312905 Farquhar G. D. Ehleringer J. R. Hubick K. T. (1989). Carbon isotope discrimination and photosynthesis. Annu. Rev. Plant Physiol. Plant Mol. Biol. 40, 503537. 10.1146/annurev.pp.40.060189.002443 Farquhar G. D. O'Leary M. H. Berry J. A. (1982). On the relationship between carbon isotope discrimination and the intercellular carbon dioxide concentration in leaves. Aust. J. Plant Physiol. 9, 121137. 10.1071/PP9820121 Ferrio J. P. Araus J. L. Buxó R. Voltas J. Bort J. (2005). Water management practices and climate in ancient agriculture: inferences from the stable isotope composition of archaeobotanical remains. Veget. Hist. Archaeobot. 14, 510517. 10.1007/s00334-005-0062-2 Fideghelli C. De Salvador F. R. (2009). World hazelnut situation and perspectives. Acta Horticult. 845, 3952. 10.17660/ActaHortic.2009.845.2 France R. (1996). Carbon isotope ratios in logged and unlogged boreal forests: examination of the potential for determining wildlife habitat use. Environ. Manage. 20, 249255. 10.1007/BF01204009 Gaillard M.-J. (1984). “Water-level changes, climate and human impact: a palaeohydrological study of Krageholm lake (Scania, southern Sweden),” in Climatic Changes on a Yearly to Millennial Basis: Geological, Historical and Instrumental Records, eds N.-A. Mörner, and W. Karlén (Dordrecht: Springer Netherlands), 147154. Giesecke T. Bennett K. D. Birks H. J. B. Bjune A. E. Bozilova E. Feurdean A. . (2011). The pace of Holocene vegetation change – testing for synchronous developments. Quat. Sci. Rev. 30, 28052814. 10.1016/j.quascirev.2011.06.014 Girona J. Cohen M. Mata M. Marsal J. Miravete C. (1994). Physiological, growth and yield responses of hazelnut (Corylus avellana L.) to different irrigation regimes. Acta Horticult. 351, 463472. 10.17660/ActaHortic.1994.351.50 Graham H. V. Patzkowsky M. E. Wing S. L. Parker G. G. Fogel M. L. Freeman K. H. (2014). Isotopic characteristics of canopies in simulated leaf assemblages. Geochim. Cosmochim. Acta 144, 8295. 10.1016/j.gca.2014.08.032 Hagerup O. (1942). The Morphology and Biology of the Corylus-Fruit. I kommission hos E. Copenhagen: Munksgaard. Hellman S. (2007). Validating and testing the landscape reconstruction algorithm in Southern Sweden - towards quantitative reconstruction of past vegetation (Ph.D. dissertation). Linnaeus University, Kalmar, Sweden. Holguin A. Charles M. Mithen S. Schulting R. (2022). In a nutshell: using structural and chemical changes to establish the charring conditions of archaeological hazelnut shells. J. Archaeol. Sci. 144, 105623. 10.1016/j.jas.2022.105623 Huntley B. Birks H. J. B. (1983). An Atlas of Past and Present Pollen Maps for Europe : 0-13,000 Years Ago. Cambridge: Cambridge University Press. Jones G. (2000). “Evaluating the importance of cultivation and collection in Neolithic Britain,” in Plants in Neolithic Britain and Beyond, ed A. S. Fairbairn (Oxford: Oxbow Books), 7984. Kohn M. J. (2010). Carbon isotope compositions of terrestrial C3 plants as indicators of (paleo)ecology and (paleo)climate. Proc. Nat. Acad. Sci. U. S. A. 107, 1969119695. 10.1073/pnas.100493310721041671 Külahçilar A. Tonkaz T. Bostan S. Z. (2018). Effect of irrigation regimes by mini sprinkler on yield and pomological traits in ‘Tombul' hazelnut. Acta Hortic. 1226, 301308. 10.17660/ActaHortic.2018.1226.45 Lagerås P. Fredh E. D. (2020). Long-term development of landscape openness and arable land use in an agricultural region of southern Sweden: the potential of REVEALS estimates using pollen records from wells. Veget. Hist. Archaeobot. 29, 113124. 10.1007/s00334-019-00756-z Larsson L. Sjöström A. (2010). Mesolithic research in the bog Rönneholms mosse, southern Sweden. Mesolithic Miscellany 21, 29. Licursi D. Antonetti C. Fulignati S. Vitolo S. Puccini M. Ribechini E. . (2017). In-depth characterization of valuable char obtained from hydrothermal conversion of hazelnut shells to levulinic acid. Bioresour. Technol. 244, 880888. 10.1016/j.biortech.2017.08.01228847076 Mingeau M. Ameglio T. Pons B. Rousseau P. (1994). Effects of water stress on development growth and yield of hazelnut trees. Acta Horticult. 351, 305314. 10.17660/ActaHortic.1994.351.33 Müller C. Hethke M. Riedel F. Helle G. (2020). Inter- and intra-tree variability of carbon and oxygen stable isotope ratios of modern pollen from nine European tree species. PLoS ONE 15, e0234315. 10.1371/journal.pone.023431532516353 O'Dwyer R. Marquer L. Trondman A.-K. Jönsson A. M. (2021). Spatially continuous land-cover reconstructions through the holocene in Southern Sweden. Ecosystems 24, 14501467. 10.1007/s10021-020-00594-5 Perry D. (1999). Vegetative tissues from mesolithic sites in the Northern Netherlands. Curr. Anthropol. 40, 231237. 10.1086/200008 Portarena S. Gavrichkova O. Brugnoli E. Battistelli A. Proietti S. Moscatello S. . (2022). Carbon allocation strategies and water uptake in young grafted and own-rooted hazelnut (Corylus avellana L.) cultivars. Tree Physiol. 42, 939957. 10.1093/treephys/tpab16434875099 Regnell M. (2012). Plant subsistence and environment at the Mesolithic site Tågerup, southern Sweden: new insights on the “Nut Age.” Veget. Hist. Archaeobot. 21, 116. 10.1007/s00334-011-0299-x Roberts J. Hopkins R. Morecroft M. (1999). Towards a predictive description of forest canopies from litter properties. Funct. Ecol. 13, 265272. 10.1046/j.1365-2435.1999.00312.x Schleser G. H. (1992). (delta)(13)C pattern in a forest tree as an indicator of carbon transfer in trees. Ecology 73, 19221925. 10.2307/1940045 Sjöström A. (1997). Ringsjöholm: a Boreal–Early Atlantic Settlement in Central Scania, Sweden. Lund Archaeol. Rev. 3, 520. Sjöström A. (2004). Rönneholm 6-10, 12, 14 och 15: Arkeologisk undersökning av ett mesolitiskt boplatskomplex i Rönneholms mosse: Hassle 32:18, Stehag socken, Eslövs kommun, Skåne, Rapporter från Institutionen för arkeologi och antikens historia, Lunds universitet. Lund: Institutionen för Arkeologi och Antikens Historia, Lunds Universitet. Sjöström A. (2013). Slabälta 1 - en boplats från sen maglemosetid vid Ageröds mosse: Arkeologisk undersökning 2012: Munkarp 2:3, Munkarp socken, Höörs kommun, Skåne, Rapporter från Institutionen för arkeologi och antikens historia, Lunds universitet. Lund: Institutionen för Arkeologi och Antikens Historia, Lunds Universitet. Tombesi A. Rosati A. (1997). Hazelnut response to water levels in relation to productive cycle. Acta Horticult. 445, 269-278. 10.17660/ActaHortic.1997.445.36 van der Merwe N. J. Medina E. (1991). The canopy effect, carbon isotope ratios and foodwebs in amazonia. J. Archaeol. Sci. 18, 249259. 10.1016/0305-4403(91)90064-V Vogel J. C. (1978). Recycling of carbon in a forest environment. Oecol. Plant. 13, 8994. Warren C. R. McGrath J. F. Adams M. A. (2001). Water availability and carbon isotope discrimination in conifers. Oecologia 127, 476486. 10.1007/s00442000060928547484 Wastegård S. (2022). The Holocene of Sweden – a review. GFF 144, 126149. 10.1080/11035897.2022.2086290 Yao Q. Mehlenbacher S. A. (2000). Heritability, variance components and correlation of morphological and phenological traits in hazelnut. Plant Breed. 119, 369381. 10.1046/j.1439-0523.2000.00524.x Zvelebil M. (1994). Plant use in the mesolithic and its role in the transition to farming. Proc. Prehist. Soc. 60, 3574. 10.1017/S0079497X00003388
      ‘Oh, my dear Thomas, you haven’t heard the terrible news then?’ she said. ‘I thought you would be sure to have seen it placarded somewhere. Alice went straight to her room, and I haven’t seen her since, though I repeatedly knocked at the door, which she has locked on the inside, and I’m sure it’s most unnatural of her not to let her own mother comfort her. It all happened in a moment: I have always said those great motor-cars shouldn’t be allowed to career about the streets, especially when they are all paved with cobbles as they are at Easton Haven, which are{331} so slippery when it’s wet. He slipped, and it went over him in a moment.’ My thanks were few and awkward, for there still hung to the missive a basting thread, and it was as warm as a nestling bird. I bent low--everybody was emotional in those days--kissed the fragrant thing, thrust it into my bosom, and blushed worse than Camille. "What, the Corner House victim? Is that really a fact?" "My dear child, I don't look upon it in that light at all. The child gave our picturesque friend a certain distinction--'My husband is dead, and this is my only child,' and all that sort of thing. It pays in society." leave them on the steps of a foundling asylum in order to insure [See larger version] Interoffice guff says you're planning definite moves on your own, J. O., and against some opposition. Is the Colonel so poor or so grasping—or what? Albert could not speak, for he felt as if his brains and teeth were rattling about inside his head. The rest of[Pg 188] the family hunched together by the door, the boys gaping idiotically, the girls in tears. "Now you're married." The host was called in, and unlocked a drawer in which they were deposited. The galleyman, with visible reluctance, arrayed himself in the garments, and he was observed to shudder more than once during the investiture of the dead man's apparel. HoME香京julia种子在线播放 ENTER NUMBET 0016www.ifdi.com.cn
      lxchain.com.cn
      livecn.com.cn
      etxyse.com.cn
      www.lianlao.org.cn
      www.ieuhqx.com.cn
      www.suyin.net.cn
      qclvyou.net.cn
      www.xfchain.com.cn
      wzpetq.com.cn
      处女被大鸡巴操 强奸乱伦小说图片 俄罗斯美女爱爱图 调教强奸学生 亚洲女的穴 夜来香图片大全 美女性强奸电影 手机版色中阁 男性人体艺术素描图 16p成人 欧美性爱360 电影区 亚洲电影 欧美电影 经典三级 偷拍自拍 动漫电影 乱伦电影 变态另类 全部电 类似狠狠鲁的网站 黑吊操白逼图片 韩国黄片种子下载 操逼逼逼逼逼 人妻 小说 p 偷拍10幼女自慰 极品淫水很多 黄色做i爱 日本女人人体电影快播看 大福国小 我爱肏屄美女 mmcrwcom 欧美多人性交图片 肥臀乱伦老头舔阴帝 d09a4343000019c5 西欧人体艺术b xxoo激情短片 未成年人的 插泰国人夭图片 第770弾み1 24p 日本美女性 交动态 eee色播 yantasythunder 操无毛少女屄 亚洲图片你懂的女人 鸡巴插姨娘 特级黄 色大片播 左耳影音先锋 冢本友希全集 日本人体艺术绿色 我爱被舔逼 内射 幼 美阴图 喷水妹子高潮迭起 和后妈 操逼 美女吞鸡巴 鸭个自慰 中国女裸名单 操逼肥臀出水换妻 色站裸体义术 中国行上的漏毛美女叫什么 亚洲妹性交图 欧美美女人裸体人艺照 成人色妹妹直播 WWW_JXCT_COM r日本女人性淫乱 大胆人艺体艺图片 女同接吻av 碰碰哥免费自拍打炮 艳舞写真duppid1 88电影街拍视频 日本自拍做爱qvod 实拍美女性爱组图 少女高清av 浙江真实乱伦迅雷 台湾luanlunxiaoshuo 洛克王国宠物排行榜 皇瑟电影yy频道大全 红孩儿连连看 阴毛摄影 大胆美女写真人体艺术摄影 和风骚三个媳妇在家做爱 性爱办公室高清 18p2p木耳 大波撸影音 大鸡巴插嫩穴小说 一剧不超两个黑人 阿姨诱惑我快播 幼香阁千叶县小学生 少女妇女被狗强奸 曰人体妹妹 十二岁性感幼女 超级乱伦qvod 97爱蜜桃ccc336 日本淫妇阴液 av海量资源999 凤凰影视成仁 辰溪四中艳照门照片 先锋模特裸体展示影片 成人片免费看 自拍百度云 肥白老妇女 女爱人体图片 妈妈一女穴 星野美夏 日本少女dachidu 妹子私处人体图片 yinmindahuitang 舔无毛逼影片快播 田莹疑的裸体照片 三级电影影音先锋02222 妻子被外国老头操 观月雏乃泥鳅 韩国成人偷拍自拍图片 强奸5一9岁幼女小说 汤姆影院av图片 妹妹人艺体图 美女大驱 和女友做爱图片自拍p 绫川まどか在线先锋 那么嫩的逼很少见了 小女孩做爱 处女好逼连连看图图 性感美女在家做爱 近距离抽插骚逼逼 黑屌肏金毛屄 日韩av美少女 看喝尿尿小姐日逼色色色网图片 欧美肛交新视频 美女吃逼逼 av30线上免费 伊人在线三级经典 新视觉影院t6090影院 最新淫色电影网址 天龙影院远古手机版 搞老太影院 插进美女的大屁股里 私人影院加盟费用 www258dd 求一部电影里面有一个二猛哥 深肛交 日本萌妹子人体艺术写真图片 插入屄眼 美女的木奶 中文字幕黄色网址影视先锋 九号女神裸 和骚人妻偷情 和潘晓婷做爱 国模大尺度蜜桃 欧美大逼50p 西西人体成人 李宗瑞继母做爱原图物处理 nianhuawang 男鸡巴的视屏 � 97免费色伦电影 好色网成人 大姨子先锋 淫荡巨乳美女教师妈妈 性nuexiaoshuo WWW36YYYCOM 长春继续给力进屋就操小女儿套干破内射对白淫荡 农夫激情社区 日韩无码bt 欧美美女手掰嫩穴图片 日本援交偷拍自拍 入侵者日本在线播放 亚洲白虎偷拍自拍 常州高见泽日屄 寂寞少妇自卫视频 人体露逼图片 多毛外国老太 变态乱轮手机在线 淫荡妈妈和儿子操逼 伦理片大奶少女 看片神器最新登入地址sqvheqi345com账号群 麻美学姐无头 圣诞老人射小妞和强奸小妞动话片 亚洲AV女老师 先锋影音欧美成人资源 33344iucoom zV天堂电影网 宾馆美女打炮视频 色五月丁香五月magnet 嫂子淫乱小说 张歆艺的老公 吃奶男人视频在线播放 欧美色图男女乱伦 avtt2014ccvom 性插色欲香影院 青青草撸死你青青草 99热久久第一时间 激情套图卡通动漫 幼女裸聊做爱口交 日本女人被强奸乱伦 草榴社区快播 2kkk正在播放兽骑 啊不要人家小穴都湿了 www猎奇影视 A片www245vvcomwwwchnrwhmhzcn 搜索宜春院av wwwsee78co 逼奶鸡巴插 好吊日AV在线视频19gancom 熟女伦乱图片小说 日本免费av无码片在线开苞 鲁大妈撸到爆 裸聊官网 德国熟女xxx 新不夜城论坛首页手机 女虐男网址 男女做爱视频华为网盘 激情午夜天亚洲色图 内裤哥mangent 吉沢明歩制服丝袜WWWHHH710COM 屌逼在线试看 人体艺体阿娇艳照 推荐一个可以免费看片的网站如果被QQ拦截请复制链接在其它浏览器打开xxxyyy5comintr2a2cb551573a2b2e 欧美360精品粉红鲍鱼 教师调教第一页 聚美屋精品图 中韩淫乱群交 俄罗斯撸撸片 把鸡巴插进小姨子的阴道 干干AV成人网 aolasoohpnbcn www84ytom 高清大量潮喷www27dyycom 宝贝开心成人 freefronvideos人母 嫩穴成人网gggg29com 逼着舅妈给我口交肛交彩漫画 欧美色色aV88wwwgangguanscom 老太太操逼自拍视频 777亚洲手机在线播放 有没有夫妻3p小说 色列漫画淫女 午间色站导航 欧美成人处女色大图 童颜巨乳亚洲综合 桃色性欲草 色眯眯射逼 无码中文字幕塞外青楼这是一个 狂日美女老师人妻 爱碰网官网 亚洲图片雅蠛蝶 快播35怎么搜片 2000XXXX电影 新谷露性家庭影院 深深候dvd播放 幼齿用英语怎么说 不雅伦理无需播放器 国外淫荡图片 国外网站幼幼嫩网址 成年人就去色色视频快播 我鲁日日鲁老老老我爱 caoshaonvbi 人体艺术avav 性感性色导航 韩国黄色哥来嫖网站 成人网站美逼 淫荡熟妇自拍 欧美色惰图片 北京空姐透明照 狼堡免费av视频 www776eom 亚洲无码av欧美天堂网男人天堂 欧美激情爆操 a片kk266co 色尼姑成人极速在线视频 国语家庭系列 蒋雯雯 越南伦理 色CC伦理影院手机版 99jbbcom 大鸡巴舅妈 国产偷拍自拍淫荡对话视频 少妇春梦射精 开心激动网 自拍偷牌成人 色桃隐 撸狗网性交视频 淫荡的三位老师 伦理电影wwwqiuxia6commqiuxia6com 怡春院分站 丝袜超短裙露脸迅雷下载 色制服电影院 97超碰好吊色男人 yy6080理论在线宅男日韩福利大全 大嫂丝袜 500人群交手机在线 5sav 偷拍熟女吧 口述我和妹妹的欲望 50p电脑版 wwwavtttcon 3p3com 伦理无码片在线看 欧美成人电影图片岛国性爱伦理电影 先锋影音AV成人欧美 我爱好色 淫电影网 WWW19MMCOM 玛丽罗斯3d同人动画h在线看 动漫女孩裸体 超级丝袜美腿乱伦 1919gogo欣赏 大色逼淫色 www就是撸 激情文学网好骚 A级黄片免费 xedd5com 国内的b是黑的 快播美国成年人片黄 av高跟丝袜视频 上原保奈美巨乳女教师在线观看 校园春色都市激情fefegancom 偷窥自拍XXOO 搜索看马操美女 人本女优视频 日日吧淫淫 人妻巨乳影院 美国女子性爱学校 大肥屁股重口味 啪啪啪啊啊啊不要 操碰 japanfreevideoshome国产 亚州淫荡老熟女人体 伦奸毛片免费在线看 天天影视se 樱桃做爱视频 亚卅av在线视频 x奸小说下载 亚洲色图图片在线 217av天堂网 东方在线撸撸-百度 幼幼丝袜集 灰姑娘的姐姐 青青草在线视频观看对华 86papa路con 亚洲1AV 综合图片2区亚洲 美国美女大逼电影 010插插av成人网站 www色comwww821kxwcom 播乐子成人网免费视频在线观看 大炮撸在线影院 ,www4KkKcom 野花鲁最近30部 wwwCC213wapwww2233ww2download 三客优最新地址 母亲让儿子爽的无码视频 全国黄色片子 欧美色图美国十次 超碰在线直播 性感妖娆操 亚洲肉感熟女色图 a片A毛片管看视频 8vaa褋芯屑 333kk 川岛和津实视频 在线母子乱伦对白 妹妹肥逼五月 亚洲美女自拍 老婆在我面前小说 韩国空姐堪比情趣内衣 干小姐综合 淫妻色五月 添骚穴 WM62COM 23456影视播放器 成人午夜剧场 尼姑福利网 AV区亚洲AV欧美AV512qucomwwwc5508com 经典欧美骚妇 震动棒露出 日韩丝袜美臀巨乳在线 av无限吧看 就去干少妇 色艺无间正面是哪集 校园春色我和老师做爱 漫画夜色 天海丽白色吊带 黄色淫荡性虐小说 午夜高清播放器 文20岁女性荫道口图片 热国产热无码热有码 2015小明发布看看算你色 百度云播影视 美女肏屄屄乱轮小说 家族舔阴AV影片 邪恶在线av有码 父女之交 关于处女破处的三级片 极品护士91在线 欧美虐待女人视频的网站 享受老太太的丝袜 aaazhibuo 8dfvodcom成人 真实自拍足交 群交男女猛插逼 妓女爱爱动态 lin35com是什么网站 abp159 亚洲色图偷拍自拍乱伦熟女抠逼自慰 朝国三级篇 淫三国幻想 免费的av小电影网站 日本阿v视频免费按摩师 av750c0m 黄色片操一下 巨乳少女车震在线观看 操逼 免费 囗述情感一乱伦岳母和女婿 WWW_FAMITSU_COM 偷拍中国少妇在公车被操视频 花也真衣论理电影 大鸡鸡插p洞 新片欧美十八岁美少 进击的巨人神thunderftp 西方美女15p 深圳哪里易找到老女人玩视频 在线成人有声小说 365rrr 女尿图片 我和淫荡的小姨做爱 � 做爱技术体照 淫妇性爱 大学生私拍b 第四射狠狠射小说 色中色成人av社区 和小姨子乱伦肛交 wwwppp62com 俄罗斯巨乳人体艺术 骚逼阿娇 汤芳人体图片大胆 大胆人体艺术bb私处 性感大胸骚货 哪个网站幼女的片多 日本美女本子把 色 五月天 婷婷 快播 美女 美穴艺术 色百合电影导航 大鸡巴用力 孙悟空操美少女战士 狠狠撸美女手掰穴图片 古代女子与兽类交 沙耶香套图 激情成人网区 暴风影音av播放 动漫女孩怎么插第3个 mmmpp44 黑木麻衣无码ed2k 淫荡学姐少妇 乱伦操少女屄 高中性爱故事 骚妹妹爱爱图网 韩国模特剪长发 大鸡巴把我逼日了 中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片 大胆女人下体艺术图片 789sss 影音先锋在线国内情侣野外性事自拍普通话对白 群撸图库 闪现君打阿乐 ady 小说 插入表妹嫩穴小说 推荐成人资源 网络播放器 成人台 149大胆人体艺术 大屌图片 骚美女成人av 春暖花开春色性吧 女亭婷五月 我上了同桌的姐姐 恋夜秀场主播自慰视频 yzppp 屄茎 操屄女图 美女鲍鱼大特写 淫乱的日本人妻山口玲子 偷拍射精图 性感美女人体艺木图片 种马小说完本 免费电影院 骑士福利导航导航网站 骚老婆足交 国产性爱一级电影 欧美免费成人花花性都 欧美大肥妞性爱视频 家庭乱伦网站快播 偷拍自拍国产毛片 金发美女也用大吊来开包 缔D杏那 yentiyishu人体艺术ytys WWWUUKKMCOM 女人露奶 � 苍井空露逼 老荡妇高跟丝袜足交 偷偷和女友的朋友做爱迅雷 做爱七十二尺 朱丹人体合成 麻腾由纪妃 帅哥撸播种子图 鸡巴插逼动态图片 羙国十次啦中文 WWW137AVCOM 神斗片欧美版华语 有气质女人人休艺术 由美老师放屁电影 欧美女人肉肏图片 白虎种子快播 国产自拍90后女孩 美女在床上疯狂嫩b 饭岛爱最后之作 幼幼强奸摸奶 色97成人动漫 两性性爱打鸡巴插逼 新视觉影院4080青苹果影院 嗯好爽插死我了 阴口艺术照 李宗瑞电影qvod38 爆操舅母 亚洲色图七七影院 被大鸡巴操菊花 怡红院肿么了 成人极品影院删除 欧美性爱大图色图强奸乱 欧美女子与狗随便性交 苍井空的bt种子无码 熟女乱伦长篇小说 大色虫 兽交幼女影音先锋播放 44aad be0ca93900121f9b 先锋天耗ばさ无码 欧毛毛女三级黄色片图 干女人黑木耳照 日本美女少妇嫩逼人体艺术 sesechangchang 色屄屄网 久久撸app下载 色图色噜 美女鸡巴大奶 好吊日在线视频在线观看 透明丝袜脚偷拍自拍 中山怡红院菜单 wcwwwcom下载 骑嫂子 亚洲大色妣 成人故事365ahnet 丝袜家庭教mp4 幼交肛交 妹妹撸撸大妈 日本毛爽 caoprom超碰在email 关于中国古代偷窥的黄片 第一会所老熟女下载 wwwhuangsecome 狼人干综合新地址HD播放 变态儿子强奸乱伦图 强奸电影名字 2wwwer37com 日本毛片基地一亚洲AVmzddcxcn 暗黑圣经仙桃影院 37tpcocn 持月真由xfplay 好吊日在线视频三级网 我爱背入李丽珍 电影师傅床戏在线观看 96插妹妹sexsex88com 豪放家庭在线播放 桃花宝典极夜著豆瓜网 安卓系统播放神器 美美网丝袜诱惑 人人干全免费视频xulawyercn av无插件一本道 全国色五月 操逼电影小说网 good在线wwwyuyuelvcom www18avmmd 撸波波影视无插件 伊人幼女成人电影 会看射的图片 小明插看看 全裸美女扒开粉嫩b 国人自拍性交网站 萝莉白丝足交本子 七草ちとせ巨乳视频 摇摇晃晃的成人电影 兰桂坊成社人区小说www68kqcom 舔阴论坛 久撸客一撸客色国内外成人激情在线 明星门 欧美大胆嫩肉穴爽大片 www牛逼插 性吧星云 少妇性奴的屁眼 人体艺术大胆mscbaidu1imgcn 最新久久色色成人版 l女同在线 小泽玛利亚高潮图片搜索 女性裸b图 肛交bt种子 最热门有声小说 人间添春色 春色猜谜字 樱井莉亚钢管舞视频 小泽玛利亚直美6p 能用的h网 还能看的h网 bl动漫h网 开心五月激 东京热401 男色女色第四色酒色网 怎么下载黄色小说 黄色小说小栽 和谐图城 乐乐影院 色哥导航 特色导航 依依社区 爱窝窝在线 色狼谷成人 91porn 包要你射电影 色色3A丝袜 丝袜妹妹淫网 爱色导航(荐) 好男人激情影院 坏哥哥 第七色 色久久 人格分裂 急先锋 撸撸射中文网 第一会所综合社区 91影院老师机 东方成人激情 怼莪影院吹潮 老鸭窝伊人无码不卡无码一本道 av女柳晶电影 91天生爱风流作品 深爱激情小说私房婷婷网 擼奶av 567pao 里番3d一家人野外 上原在线电影 水岛津实透明丝袜 1314酒色 网旧网俺也去 0855影院 在线无码私人影院 搜索 国产自拍 神马dy888午夜伦理达达兔 农民工黄晓婷 日韩裸体黑丝御姐 屈臣氏的燕窝面膜怎么样つぼみ晶エリーの早漏チ○ポ强化合宿 老熟女人性视频 影音先锋 三上悠亚ol 妹妹影院福利片 hhhhhhhhsxo 午夜天堂热的国产 强奸剧场 全裸香蕉视频无码 亚欧伦理视频 秋霞为什么给封了 日本在线视频空天使 日韩成人aⅴ在线 日本日屌日屄导航视频 在线福利视频 日本推油无码av magnet 在线免费视频 樱井梨吮东 日本一本道在线无码DVD 日本性感诱惑美女做爱阴道流水视频 日本一级av 汤姆avtom在线视频 台湾佬中文娱乐线20 阿v播播下载 橙色影院 奴隶少女护士cg视频 汤姆在线影院无码 偷拍宾馆 业面紧急生级访问 色和尚有线 厕所偷拍一族 av女l 公交色狼优酷视频 裸体视频AV 人与兽肉肉网 董美香ol 花井美纱链接 magnet 西瓜影音 亚洲 自拍 日韩女优欧美激情偷拍自拍 亚洲成年人免费视频 荷兰免费成人电影 深喉呕吐XXⅩX 操石榴在线视频 天天色成人免费视频 314hu四虎 涩久免费视频在线观看 成人电影迅雷下载 能看见整个奶子的香蕉影院 水菜丽百度影音 gwaz079百度云 噜死你们资源站 主播走光视频合集迅雷下载 thumbzilla jappen 精品Av 古川伊织star598在线 假面女皇vip在线视频播放 国产自拍迷情校园 啪啪啪公寓漫画 日本阿AV 黄色手机电影 欧美在线Av影院 华裔电击女神91在线 亚洲欧美专区 1日本1000部免费视频 开放90后 波多野结衣 东方 影院av 页面升级紧急访问每天正常更新 4438Xchengeren 老炮色 a k福利电影 色欲影视色天天视频 高老庄aV 259LUXU-683 magnet 手机在线电影 国产区 欧美激情人人操网 国产 偷拍 直播 日韩 国内外激情在线视频网给 站长统计一本道人妻 光棍影院被封 紫竹铃取汁 ftp 狂插空姐嫩 xfplay 丈夫面前 穿靴子伪街 XXOO视频在线免费 大香蕉道久在线播放 电棒漏电嗨过头 充气娃能看下毛和洞吗 夫妻牲交 福利云点墦 yukun瑟妃 疯狂交换女友 国产自拍26页 腐女资源 百度云 日本DVD高清无码视频 偷拍,自拍AV伦理电影 A片小视频福利站。 大奶肥婆自拍偷拍图片 交配伊甸园 超碰在线视频自拍偷拍国产 小热巴91大神 rctd 045 类似于A片 超美大奶大学生美女直播被男友操 男友问 你的衣服怎么脱掉的 亚洲女与黑人群交视频一 在线黄涩 木内美保步兵番号 鸡巴插入欧美美女的b舒服 激情在线国产自拍日韩欧美 国语福利小视频在线观看 作爱小视颍 潮喷合集丝袜无码mp4 做爱的无码高清视频 牛牛精品 伊aⅤ在线观看 savk12 哥哥搞在线播放 在线电一本道影 一级谍片 250pp亚洲情艺中心,88 欧美一本道九色在线一 wwwseavbacom色av吧 cos美女在线 欧美17,18ⅹⅹⅹ视频 自拍嫩逼 小电影在线观看网站 筱田优 贼 水电工 5358x视频 日本69式视频有码 b雪福利导航 韩国女主播19tvclub在线 操逼清晰视频 丝袜美女国产视频网址导航 水菜丽颜射房间 台湾妹中文娱乐网 风吟岛视频 口交 伦理 日本熟妇色五十路免费视频 A级片互舔 川村真矢Av在线观看 亚洲日韩av 色和尚国产自拍 sea8 mp4 aV天堂2018手机在线 免费版国产偷拍a在线播放 狠狠 婷婷 丁香 小视频福利在线观看平台 思妍白衣小仙女被邻居强上 萝莉自拍有水 4484新视觉 永久发布页 977成人影视在线观看 小清新影院在线观 小鸟酱后丝后入百度云 旋风魅影四级 香蕉影院小黄片免费看 性爱直播磁力链接 小骚逼第一色影院 性交流的视频 小雪小视频bd 小视频TV禁看视频 迷奸AV在线看 nba直播 任你在干线 汤姆影院在线视频国产 624u在线播放 成人 一级a做爰片就在线看狐狸视频 小香蕉AV视频 www182、com 腿模简小育 学生做爱视频 秘密搜查官 快播 成人福利网午夜 一级黄色夫妻录像片 直接看的gav久久播放器 国产自拍400首页 sm老爹影院 谁知道隔壁老王网址在线 综合网 123西瓜影音 米奇丁香 人人澡人人漠大学生 色久悠 夜色视频你今天寂寞了吗? 菲菲影视城美国 被抄的影院 变态另类 欧美 成人 国产偷拍自拍在线小说 不用下载安装就能看的吃男人鸡巴视频 插屄视频 大贯杏里播放 wwwhhh50 233若菜奈央 伦理片天海翼秘密搜查官 大香蕉在线万色屋视频 那种漫画小说你懂的 祥仔电影合集一区 那里可以看澳门皇冠酒店a片 色自啪 亚洲aV电影天堂 谷露影院ar toupaizaixian sexbj。com 毕业生 zaixian mianfei 朝桐光视频 成人短视频在线直接观看 陈美霖 沈阳音乐学院 导航女 www26yjjcom 1大尺度视频 开平虐女视频 菅野雪松协和影视在线视频 华人play在线视频bbb 鸡吧操屄视频 多啪啪免费视频 悠草影院 金兰策划网 (969) 橘佑金短视频 国内一极刺激自拍片 日本制服番号大全magnet 成人动漫母系 电脑怎么清理内存 黄色福利1000 dy88午夜 偷拍中学生洗澡磁力链接 花椒相机福利美女视频 站长推荐磁力下载 mp4 三洞轮流插视频 玉兔miki热舞视频 夜生活小视频 爆乳人妖小视频 国内网红主播自拍福利迅雷下载 不用app的裸裸体美女操逼视频 变态SM影片在线观看 草溜影院元气吧 - 百度 - 百度 波推全套视频 国产双飞集合ftp 日本在线AV网 笔国毛片 神马影院女主播是我的邻居 影音资源 激情乱伦电影 799pao 亚洲第一色第一影院 av视频大香蕉 老梁故事汇希斯莱杰 水中人体磁力链接 下载 大香蕉黄片免费看 济南谭崔 避开屏蔽的岛a片 草破福利 要看大鸡巴操小骚逼的人的视频 黑丝少妇影音先锋 欧美巨乳熟女磁力链接 美国黄网站色大全 伦蕉在线久播 极品女厕沟 激情五月bd韩国电影 混血美女自摸和男友激情啪啪自拍诱人呻吟福利视频 人人摸人人妻做人人看 44kknn 娸娸原网 伊人欧美 恋夜影院视频列表安卓青青 57k影院 如果电话亭 avi 插爆骚女精品自拍 青青草在线免费视频1769TV 令人惹火的邻家美眉 影音先锋 真人妹子被捅动态图 男人女人做完爱视频15 表姐合租两人共处一室晚上她竟爬上了我的床 性爱教学视频 北条麻妃bd在线播放版 国产老师和师生 magnet wwwcctv1024 女神自慰 ftp 女同性恋做激情视频 欧美大胆露阴视频 欧美无码影视 好女色在线观看 后入肥臀18p 百度影视屏福利 厕所超碰视频 强奸mp magnet 欧美妹aⅴ免费线上看 2016年妞干网视频 5手机在线福利 超在线最视频 800av:cOm magnet 欧美性爱免播放器在线播放 91大款肥汤的性感美乳90后邻家美眉趴着窗台后入啪啪 秋霞日本毛片网站 cheng ren 在线视频 上原亚衣肛门无码解禁影音先锋 美脚家庭教师在线播放 尤酷伦理片 熟女性生活视频在线观看 欧美av在线播放喷潮 194avav 凤凰AV成人 - 百度 kbb9999 AV片AV在线AV无码 爱爱视频高清免费观看 黄色男女操b视频 观看 18AV清纯视频在线播放平台 成人性爱视频久久操 女性真人生殖系统双性人视频 下身插入b射精视频 明星潜规测视频 mp4 免賛a片直播绪 国内 自己 偷拍 在线 国内真实偷拍 手机在线 国产主播户外勾在线 三桥杏奈高清无码迅雷下载 2五福电影院凸凹频频 男主拿鱼打女主,高宝宝 色哥午夜影院 川村まや痴汉 草溜影院费全过程免费 淫小弟影院在线视频 laohantuiche 啪啪啪喷潮XXOO视频 青娱乐成人国产 蓝沢润 一本道 亚洲青涩中文欧美 神马影院线理论 米娅卡莉法的av 在线福利65535 欧美粉色在线 欧美性受群交视频1在线播放 极品喷奶熟妇在线播放 变态另类无码福利影院92 天津小姐被偷拍 磁力下载 台湾三级电髟全部 丝袜美腿偷拍自拍 偷拍女生性行为图 妻子的乱伦 白虎少妇 肏婶骚屄 外国大妈会阴照片 美少女操屄图片 妹妹自慰11p 操老熟女的b 361美女人体 360电影院樱桃 爱色妹妹亚洲色图 性交卖淫姿势高清图片一级 欧美一黑对二白 大色网无毛一线天 射小妹网站 寂寞穴 西西人体模特苍井空 操的大白逼吧 骚穴让我操 拉好友干女朋友3p