To investigate the impact of irrigation and aridity on crop carbon isotope values, studies have relied on outdoor growing experiments that correlate broad categories of water availability with δ¹³C and carbon isotope discrimination (Δ¹³C) values in aggregate samples (Wallace et al., 2013; Styring et al., 2016a). Crop carbon stable isotope values are primarily influenced by photosynthetic pathway and water availability (Araus et al., 1997a; Cappers and Neef, 2012; Farquhar et al., 1989; Tieszen, 1991; Vogel, 1993), although field proximity to forests where light availability may be lower would also impact crop δ¹³C values (Van Der Merwe, 1982; van der Merwe and Medina, 1991; Bonafini et al., 2013). Internal stomatal conductance of CO₂ also further modifies the carbon isotope discrimination of plants, with water-efficient flora typically exhibiting lower stomatal conductance and/or higher photosynthetic capacity that consequently impart a ¹³C-enrichment in plant tissues (Ma et al., 2021). This is reflected in the carbon isotopic composition of wheat and barley grown under the same conditions, with barley δ¹³C and Δ¹³C values ~1‰ higher than wheat (Flohr et al., 2019; Styring et al., 2016a; Wallace et al., 2013).

Modern aggregate samples show limited carbon isotope variation in crop plants grown in the same field, whereas more significant differences could exist between crops grown in different fields under similar watering conditions (Ferrio et al., 2005; Wallace et al., 2013; Jones et al., 2021). Research establishing the initial frameworks for assessing water availability drew on 168 sample sets, each consisting of ~25 to 50 aggregated grains of wheat and barley grown in arid and semi-arid fields in Span and Syria, that were collected from a mix of rainfed and irrigated water regimes. Aggregated samples of barley harvested in Spain in 2007 and 2008 demonstrated overall low variation in mean Δ¹³C values. In unirrigated (n = 2) and moderately irrigated (n = 3) barley from a single annual harvest, these differences were extremely limited (2008: 18.9 ± 0.2‰ vs. 18.7 ± 0.5‰), but fully irrigated fields (2007: n = 6, 2008: n = 3) exhibited higher mean Δ¹³C values as well as variation (2007: 18.8 ± 0.6‰, 2008: 19.6 ± 1.1‰; Wallace et al., 2013). This higher variation in fully irrigated wheat mean Δ¹³C values in different years (2007: n = 5, 2008: n = 3) was not as pronounced (2007: 17.81 ± 0.81‰, 2008: 18.14 ± 0.56‰; Wallace et al., 2013). Similar ranges of limited intra-field variation were also observed in aggregate samples (n_aggregate = 36, each of 50 grains) in plants subjected to three distinct regions with different watering in semi-arid fields in Morocco (Styring et al., 2016a). Two of the regions, Rainfed North and Rainfed South, reflect the overall water availability and precipitation exhibiting average carbon isotope values of −27.6 ± 0.5 ‰ (703 mm) and −23.6 ± 0.6‰ (272 mm), respectively. The irrigated and flood-cultivated Oasis crops exhibited a mean δ¹³C value of −26.3 ± 0.6‰ with similar watering conditions to Rainfed North. There is extremely limited overall variation between aggregate samples within each region of ± 0.6‰ (Styring et al., 2016a).

Notably, the Δ¹³C values of aggregate barley (n = 10 grains) grown on unirrigated rainfed fields located in Jordan showed greater carbon isotope differences between harvests carried in separate years but from the same field (n = 3 per year, 2005–2006: 15.4 ± 0.1‰, 2006–2007: 16.5 ± 0.2‰, 2007–8: 17.3 ± 0.1‰; Flohr et al., 2019). In addition to these temporal differences, Jones et al. (2021) documented spatial differences in Δ¹³C values, up to 2‰, between barley grown in different fields but under similar watering conditions across Northwest India. For 124 aggregate samples (n_grain = 20–30 grains each) from uncontrolled flooded (n_sample = 101, 18.0 ± 1.5‰), sprinkler-irrigated (n_sample = 19, 16.8 ± 1.4‰), and rainfed fields (n_sample = 4, 17.2 ± 2.8‰), this study found that while water availability was the primary driver of carbon isotope variation, the variation between fields with the same watering conditions could be significant (Jones et al., 2021).

Single-grain isotope analysis has found limited intra-panicle and intra-field carbon isotope variation (Heaton et al., 2009). In two bread wheat (Triticum aestivum ssp. vulgare) plants grown in the same uncontrolled field in Nottingham (UK; n_{grain =} 18), single-grain δ¹³C values ranged between −26.5 and −27.5‰. In corresponding aggregate samples from the same fields, each consisting of ~300 uncarbonized grains from 6 panicles, yielded δ¹³C standard deviations between ± 0.37 and ± 0.83‰ (Heaton et al., 2009). These studies have found consistent limited ranges of δ¹³C and Δ¹³C variation found in aggregate and single-grain samples of grains grown in the same field and year (Araus et al., 1997b; Styring et al., 2016a; Flohr et al., 2019; Wallace et al., 2013; Araus et al., 2003). However, more notable variation in carbon isotope values is possible between inter-annual harvests from the same fields, as well as spatial variation independent of watering conditions (Wallace et al., 2013; Flohr et al., 2019; Jones et al., 2021). These patterns of carbon isotope variation indicate that information of watering practices may be lost in aggregate samples if grains originate from different fields or years.

Crop nitrogen isotopic composition is driven by plant uptake of bioavailable nitrogen through either nitrate (NO3-) or ammonium (NH4+; Denk et al., 2017), both of which are affected by manuring and soil moisture (Craine et al., 2015; Szpak et al., 2017). Grain δ¹⁵N values reflect nutrient availability during grain-filling, and practices such as manuring enrich soil ¹⁵N through the direct addition of ammonia to soils, driving floral nitrogen isotope values up, sometimes considerably (up to 10‰; Fraser et al., 2011; Kanstrup et al., 2011; Szpak, 2014; Craine et al., 2015). In addition to marked ¹⁵N enrichment in grain values visible in cultivars grown under manured conditions (Bogaard et al., 2007; Chadwick et al., 2000; Aguilera et al., 2008; Fraser et al., 2013; Styring et al., 2016a), there can be differences in the amplitude of nitrogen isotopic variation for grains from manured plots compared to those from unmanured plots, but this appears to be dependent in part on regional environmental conditions. In controlled experimental plots located in temperate environments in Hertfordshire (UK), aggregate unmanured and manured samples yield widely different δ¹⁵N values but similar ranges of variation (δ¹⁵N_manured = 1.9 ± 1.4‰, range = 0.6 to 4.0‰; δ¹⁵N_manured 7.4 ± 1.1‰; Bogaard et al., 2007; Fraser et al., 2013). In contrast, aggregate samples from Morocco exhibited dramatic differences in nitrogen isotope values between manured and unmanured fields, but low nitrogen isotope variation within both unmanured fields (Rainfed North δ¹⁵N_mean 0.8 ± 0.6‰) and irrigated manured fields (Oasis δ¹⁵N_mean 14.0 ± 1.0‰), but wide variation in rainfed manured fields (Rainfed South δ¹⁵N_{mean =} 7.8 ± 2.1‰).

Analyses of single grains also found slightly greater variation in manured versus unmanured grains (Bogaard et al., 2007; Larsson et al., 2019). In plants from the same fields, single-grain intra-panicle variation in unmanured wheat δ¹⁵N values was limited to ± 0.1‰ (range = −0.5 to 0.5‰) while manured grains varied by ± 1.3‰ (range = 5 to 7.5‰; Bogaard et al., 2007). A similar pattern was identified in manured and unmanured experimental plots in Borgeby, Sweden, where intra-panicle variation in δ¹⁵N values measured from the single grains of unmanured 2-row hulled barley (Hordeum vulgare ssp. distichon) averaged 5.4 ± 0.6‰ (range = 4.1 to 6.2‰), but in manured grains 8.9 ± 1.6‰ (range = 6.4 to 11.8‰; Larsson et al., 2019). The variation observed between individual grains from the same fields (ca. 1–2‰), and the relatively high variation present between some aggregate sample wheat and barley δ¹⁵N mean values, suggests that important information on fertilization practices may be heavily obscured in these pooled samples.

Archaeobotanical assemblages are in most cases made up of plant material preserved through charring. There has been much research documenting the impact of carbonization on grain morphology and taphonomic processes that impact preservation (Boardman and Jones, 1990; Charles et al., 2015; Hillman et al., 1993; Pearsall, 2015; van der Veen, 2007). With increasing interest in the use of stable isotope analysis to understand ancient agricultural systems, this work has extended to explore how charring influences original in vivo carbon and nitrogen isotope values through experimental research on modern grains (DeNiro and Hastorf, 1985; Aguilera et al., 2008; Kanstrup et al., 2012; Fraser et al., 2013; Nitsch et al., 2015; Stroud et al., 2023a).

Grains exposed to temperatures above 200°C undergo Maillard reactions that volatilize carbon (C-) and nitrogen (N-) containing compounds as seed starches convert to dextrin (Pazola and Cieslak, 1979; Styring et al., 2013). Grain %C and %N rise by ca. 20 and 2.5%, respectively, while grains lose 20 to 40% of their mass, proportional to temperature (Czimczik et al., 2002; Braadbaart et al., 2004; Kanstrup et al., 2012). Between 230 and 300°C, these reactions lead to some loss of ¹⁴N and a corresponding increase in δ¹⁵N values (Bogaard et al., 2007; Kanstrup et al., 2012; Styring et al., 2013), while δ¹³C values remain largely unchanged (Aguilera et al., 2008; Fraser et al., 2011; Nitsch et al., 2015). At temperatures above 300°C, gross deformation in grain morphology occurs, with mass loss over 50% (Boardman and Jones, 1990; Braadbaart et al., 2004; Charles et al., 2015). This is accompanied by significant deviations from in vivo isotope values, with observed shifts of ≥1‰ in δ¹³C and ≥2‰ in δ¹⁵N values (Czimczik et al., 2002; Kanstrup et al., 2012; Fraser et al., 2013; Stroud et al., 2023a).

Many charring studies have compared the mean carbon and nitrogen isotope values between aggregate charred and uncharred samples. In grains charred between 200°C and 300°C, a mean increase of +0.6‰ in the δ¹⁵N values of 5 samples each composed of 15 grains was found when compared with uncharred aggregate samples (Kanstrup et al., 2011, 2012). Fraser et al. (2013) charred aggregate samples of 25 to 50 grains at 230°C from plots in Syria, Germany, and the United Kingdom; charred grains yielded higher δ¹⁵N values up to +0.8‰ compared to uncharred samples but no notable differences in δ¹³C values (Bogaard et al., 2007; Fraser et al., 2011).

Other research has assessed charring effects on grain isotope values using statistical correlations or linear regression models to predict isotopic offsets for archaeobotanical grain isotope samples. Aguilera et al. (2008) compared samples of five aggregate wheat or barley grains charred at 250°C with uncharred samples, testing correlations between both to suggest a δ¹⁵N offset of +0.68‰ to compensate for charring impacts. Other charring experiments on modern grains applied multiple linear regression models to isotope data sets to predict how time and temperature impact grain δ¹³C and δ¹⁵N values (Nitsch et al., 2015; Stroud et al., 2023a). In unmanured samples of eight taxa, including bread wheat (Triticum aestivum) and hulled barley (Hordeum vulgare), samples of 10 grains were fired between 230 and 300°C between 4 and 24 h and compared with uncharred samples. For every 15°C above 200°C, δ¹³C values increased by +0.05‰ and for δ¹⁵N by 0.12‰, whereas every 4 h of charring resulted in a 0.016‰ increase δ¹³C, and 0.04‰ for δ¹⁵N values (Nitsch et al., 2015; Stroud et al., 2023a). Nitsch et al. (2015) predicted overall isotopic offsets of +0.31‰ in δ¹⁵N, and−0.11‰ in δ¹³C for grains carbonized between 230 and 260°C. Stroud et al. (2023a), expanding and reanalyzing this data, argued for offsets of +0.32‰ in δ¹⁵N, and −0.16‰ in δ¹³C to be subtracted from archaeological wheat and barley grain isotope values (Stroud et al., 2023b). Notably, these valuable charring studies were based on grains randomly selected from multiple plants from a single field, carbonized and aggregated into a single sample, then compared with a corresponding aggregate sample of randomly selected sample of uncharred grains. To date, no analyses have compared the isotopic variation of single grains—either charred or uncharred, within a single panicle.

Experimental studies conducted on modern uncharred and charred grains have repeatedly demonstrated the viability of correlating carbon and nitrogen stable isotope ratios measured in seeds to cultivation practices and environmental conditions (Ferrio et al., 2005; Bogaard et al., 2013; Wallace et al., 2013; Nitsch et al., 2015; Flohr et al., 2019; Jones et al., 2021). However, most experimental research has either assessed isotopic variability between fields or examined the impact of charring on aggregated samples using modern grains grown under known conditions. Archaeobotanical grains complicate the use of aggregate samples as they cannot be assumed to derive from a single field or set of growing conditions.

The clearest argument making the case for aggregating archaeological grains derives from work applying a multiple linear regression model calculating potential charring offsets (Nitsch et al., 2015; Stroud et al., 2023a). These studies used this model to define the range of carbon and nitrogen isotope variation within a single modern growing context (i.e., a field), and define adequate sample sizes for archaeobotanical stable isotope analysis (Nitsch et al., 2015; Stroud et al., 2023a). The analysis, based on 70 aggregate samples consisting of 10 grains representing eight taxa from unmanured fields, calculated a residual standard error (SE) of ca. 0.25‰ in δ¹³C and 0.5‰ in δ¹⁵N. From this calculation, this work argued that within a 95% confidence interval, any expected variation in grain δ¹³C and δ¹⁵N values from a single field would be ~± 0.5‰ in δ¹³C and ± 1.0‰ in δ¹⁵N (1.96 × SE; Nitsch et al., 2015; Stroud et al., 2023a). Nitsch et al. (2015) hypothesized that decreasing the number of grains homogenized together within a single aggregate sample would in turn increase the standard error, increasing the variation and uncertainty of a single measurement, thereby rendering individual grain isotope values too variable to be interpreted (Stroud et al., 2023a). From these findings drawing on aggregate samples from modern growing contexts, Nitsch et al. (2015) recommended the inclusion of ten (10) archaeological grains from a single archaeological context to create a single aggregated sample. This recommendation was based on their observation that uncertainty at the 95% confidence level decreased with the inclusion of 10 specimens grown under uniform conditions; increasing the sample size beyond this point would unnecessarily consume additional archaeobotanical specimens without further reducing the inherent isotopic variability within an aggregate sample (Nitsch et al., 2015; Stroud et al., 2023b). However, averaged isotope values from aggregated samples from a single field grown in a single season may not be representative of archaeobotanical assemblages of charred seeds that potentially span multiple unknown temporal and spatial scales. This underscores the necessity of single-grain isotope analysis to capture the isotopic variation present in past cultivation practices.

We performed carbon (δ¹³C) and nitrogen (δ¹⁵N) isotope analyses on single grains from the same panicle to better define the range of intra-panicle isotopic variation for unmanured crops growing in a well-watered C₃ environment and assess how charring might impact that range in variation. Four panicles of naked barley (Hordeum vulgare var. nudum) were selected from legacy collections generated by growing experiments conducted by Institute of Prehistoric and Protohistoric Archaeology, University of Kiel at Steinzeitpark Dithmarschen, Albersdorf, Germany (Figure 1). The ca. 5 m × 5 m plot lies on the grounds of the Steinzeitpark Dithmarschen, with sandy clay soils that have not been manured or received fertilizer since at least 2005, when the park was founded (Burbaum et al., 2019; Beuker, 2020). The barley was sown by hand in April 2017, left to grow with no intervention until harvest, and received no manure or additional water beyond environmental precipitation of ca. 315 mm (Deutscher Wetterdienst., 2024), with mature ripe plants harvested in August 2017. Each panicle was assigned a letter (A, B, D, E). Grains on each panicle were sequentially numbered in ascending order, from the base to the apex, recording their relative position (Figure 2). All grains from two panicles, D (n = 33) and E (n = 30), were individually measured for carbon and nitrogen isotopes in order to assess intra-panicle isotopic variation within and between single inflorescences. These uncarbonized grains were ground individually into powder using an agate mortar and pestle and then weighed (2.5 mg) into tin capsules for carbon and nitrogen isotope ratio mass spectrometry.

Plant remains recovered from archaeological contexts are typically preserved through exposure to fire (Gallagher, 2015). Studies of charring impacts on grain stable isotope values using aggregate samples have found minimal increases to grains charred between 230 and 300°C (Fraser et al., 2013; Vaiglova et al., 2014). To assess how charring potentially impacts the range of intra-panicle isotopic variation, we measured individual grains subjected to a range of charring conditions. Four unaltered grains each from panicles A and B were selected for mass spectrometry and set aside. Additional 26 grains from panicles A and B were individually carbonized in both anoxic and oxidizing conditions at 230°C (n = 8), 300°C (n = 16), and 400°C (n = 4) in a muffle furnace (Table 1, Supplementary material 1). Grains carbonized in anoxic conditions were individually wrapped in aluminum foil, with each foil packet then placed together in an enclosed ceramic crucible. Grains carbonized in oxidizing conditions were placed individually in an open ceramic crucible. All uncarbonized and carbonized grains were ground individually into powder using an agate mortar and pestle and then weighed (2.5 mg) into tin capsules for carbon and nitrogen isotope ratio mass spectrometry.

We expand our exploration of single-grain isotope analysis by examining the degree to which carbon and nitrogen stable isotope ratios derived from an aggregated pool of single grains encompass the range of isotopic variation exhibited by those single grains. This is done through the stable isotope analysis of aggregated and single-grain sample sets all obtained from the same primary archaeological context. Charred samples of hulled barley (Hordeum vulgare) were collected from a burned storage bin made of wattle and daub discovered by the Harappa Archaeological Project (HARP) in the eroded upper levels of Trench 42, on Mound AB at the Indus site of Harappa, Pakistan (Figure 1, ca. 1900–1700 BCE; Meadow and Kenoyer, 2001). This charred bin and associated wash layers contained thousands of charred barley grains (Meadow and Kenoyer, 2001). Cemetery H ceramics associated with the context and two samples of charred barley that have been radiocarbon dated confirm that these samples date to Periods 4/5 ca. 1900–1700 BCE (Meadow and Kenoyer, 2001).

Covering over 150 ha in size, Harappa was a large urban center that may have supported a population of over ca. 40,000 people (Kenoyer, 2008; Kenoyer and Meadow, 2016). Harappa was linked with other Indus cities by trade and a regionally integrated material culture that included similar ceramic traditions and systems of weights (Kenoyer, 1997; Possehl, 2002), set within a mosaic of dense rural settlements spanning a wide range of environmental conditions (Weber et al., 2010a; Petrie and Bates, 2017). Harappa is located adjacent to a paleochannel of the Ravi River and today receives ca. 250 mm annual precipitation with ca. 90 mm falling between November to April and ca. 160 mm between May to October (Pendall et al., 1990). Like other Indus sites, Harappa agricultural systems were structured by the South Asian monsoon (Bates, 2022; Clift and d'Alpoim Guedes, 2021; Weber, 2003; Fuller, 2006). Traditional farming practices for the region take advantage of seasonal precipitation that support two cropping seasons: the summer (kharif) season when crops such as millets and grams are grown during the June to October monsoonal season, and the winter (rabi) seasons when barley, wheat, and lentils are grown between November to April (Cappers et al., 2016; Miller, 2006). Harappan agriculture relied heavily on Rabi cropping, suggested by the high ubiquity, counts, and weight of charred barley, followed by winter grown wheat and lentil grains found in the archaeobotanical assemblage (James et al., 2024; Weber, 2003; Weber et al., 2010a).

Large-grained domesticated barley, the primary crop cultivated throughout the duration of occupation at Harappa, was intensively exploited (James et al., 2024), although millets and grams increased in importance by 2,600 BCE (Weber et al., 2010a,b). Floodplain cultivation or some form of localized irrigation potentially provided barley fields with an additional water source required to meet barley watering requirements (Miller, 2006, 2015; You, 2019). Crop production at Harappa may have been further enhanced through the use of zebu cattle and water buffalo (Meadow, 1996; Patel and Meadow, 2017), which served as a source of manure (Lancelotti, 2018; James et al., 2024), as well as traction, suggested by pathologies on their skeletal extremities (Miller, 2004), and terracotta models of cattle pulling carts and plows (Kenoyer, 2004). Secondary products formed a substantial component of Harappa agricultural systems, with cattle, sheep, and goats also husbanded for their meat and milk (Meadow, 1996; Miller, 2003; Patel and Meadow, 2017).

A total of 46 carbonized grains were selected from the Trench 42 burnt storage feature. Grains were randomly selected and placed into two sample groups, the first representing single grains (S) each individually analyzed for carbon and nitrogen isotope ratios, and the second of aggregated samples (A), representing multiple seeds mechanically homogenized for subsequent isotope analyses. To assess if the range of isotopic variation measured from randomly selected single grains (S1) matched a paired sample of aggregated seeds (A1), 20 grains were collected from the Trench 42 barley.

Ten (S1; grain n = 10) of these grains were each individually ground into a fine powder using an agate mortar with ~2.5 mg of the resulting powder weighed into a tin capsule for single-grain measurements. The remaining 10 grains (total number of grains in sample A1, grain n = 10) were collectively homogenized into a single aggregate sample A1 (total number of grains in sample A1, n = 10). The powder from A1 was measured six times, and weighed into tin capsules (A1; n = 6).

Carbon (δ¹³C) and nitrogen (δ¹⁵N) isotope values and range of isotopic variation exhibited by an aggregated sample, representing multiple barley seeds, were then compared to the isotopic composition of single grains also used to create the aggregate sample sets. Twenty-six single grains (S2, n = 26) were individually homogenized with ~2.5 mg of each grain's powder separately weighed into tin capsules for single-grain stable isotope measurements. The remaining powders from the S1 individual grains were retained and homogenized together into the aggregate sample A2 (total number of grains in sample, n = 26). The powder from A2 was measured six times, with each sample weighed into tin capsules (A2; n = 6).

All archaeological powders were assessed by Attenuated Total Reflectance—Fourier-transform infrared spectroscopy (ATR-FTIR) for contamination and diagenetic alteration (Styring et al., 2013; Vaiglova et al., 2014). All ATR-FTIR data was preprocessed within OPUS software ver 8.1, utilizing the software's built-in baseline corrections and vector normalization. Existing screening protocols rely on observable wavelength peaks indicating high concentrations of contaminating nitrate (1,085 cm⁻¹, 1,450 cm⁻¹, 3,300 cm⁻¹), carbonate (720 cm⁻¹, 870 cm⁻¹), or humic contamination (1,010 cm⁻¹, 1,080 cm⁻¹; Vaiglova et al., 2014; Brinkkemper et al., 2018). No peaks indicating contamination were found, so no pretreatment was conducted.

Carbon and nitrogen isotope analyses were undertaken at the Archaeology Stable Isotope Laboratory (ASIL), Institute for Prehistoric and Protohistoric Archaeology, University of Kiel. Samples were measured using an isoprime visION continuous flow isotope ratio mass spectrometer coupled to a vario PYRO cube elemental analyzer (Elementar Analysesysteme GmbH, Langenselbold, Germany). Stable carbon and nitrogen isotope values were calibrated relative to the VPDB (δ¹³C) and AIR (δ¹⁵N) scales using the glutamic acid standards USGS40 (δ¹³C −26.39 ± 0.04‰, δ¹⁵N −4.52 ± 0.06‰) and USGS41a (δ¹³C 36.55 ± 0.08‰, δ¹⁵N 47.55 ± 0.15‰; Qi et al., 2003, 2016). Samples were measured in 14 analytical runs. Measurement uncertainty was monitored using wheat and millet internal standards. The isotope values reported here for internal standards represent long-term averages calibrated to VPDB and AIR with USGS40 and USGS41a (Supplementary material 1). Precision [u(R_w)] was ± 0.07‰ for δ¹³C and ± 0.16‰ for δ¹⁵N on the basis of repeated measurements of calibration standards, check standards, and sample replicates. Accuracy [u_(bias)] was determined to be ± 0.12‰ for δ¹³C and ± 0.13‰ for δ¹⁵N on the basis of the difference between the observed and known δ values of the check standards and the long-term standard deviations of these check standards. The total analytical uncertainty was ± 0.19‰ for δ¹³C and ± 0.11 for δ¹⁵N (after Szpak et al., 2017).

Carbon isotopes in archaeological grains are most often expressed in analysis as carbon isotope discrimination (Δ¹³C) to accommodate for differences in past and modern atmospheric carbon dioxide (CO₂) concentrations and atmospheric δ¹³C (Farquhar et al., 1989; Araus et al., 1997a; Wallace et al., 2015; Rosen et al., 2019). Here, all δ¹³C values measured from archaeological seeds are presented as Δ¹³C calculated using an estimated δ¹³C_air of 6.4‰ (Farquhar et al., 1989; Francey et al., 1999; Ferrio et al., 2005). This calculation indicates the carbon isotope values of the archaeological seed in relation to the source of atmospheric CO₂ i.e., Δ reflects the difference in δ¹³C values between the air and plant tissue (Farquhar et al., 1989).

Uncarbonized grains (n = 71) from the four panicles of modern, experimentally grown barley exhibited low variation in δ¹³C averaging −28.9 ± 0.3‰ (range = −29.5 to −28.5‰; Table 1, Figure 3). The number of grains sampled per panicle influenced the isotopic variation expressed within each panicle. Panicle A (n = 4) yielded on average a δ¹³C value of −28.7 ± 0.2‰ (range = −28.8 to −28.4‰), Panicle B −29.1 ± 0.3‰ (n = 4, range = −29.5 to −28.8‰), Panicle D (n = 33) −28.9 ± 0.2‰ (n = 33, range = −29.5 to −28.5‰), and Panicle E −28.9 ± 0.2‰ (n = 30, range = −29.4 to −28.3‰).

For nitrogen isotopes, uncarbonized grains yielded an average δ¹⁵N value of 0.4 ± 0.6‰ (range = −1.2 to 1.3‰). Panicle A exhibited a mean δ¹⁵N value of 1.2 ± 0.1‰ (range = 1.2 to 1.3‰), Panicle B 1.0 ± 0.3‰ (range = 0.6 to 1.0‰), Panicle D 0.7 ± 0.3‰ (range = −0.6 to +1.2‰), and Panicle E −0.2 ± 0.4‰ (range = −1.2 to 0.4‰). No differences in δ¹³C or δ¹⁵N values based on the position of each grain on its respective panicle were observed (Figure 3).

Panicle A grains carbonized at 230°C (n = 6) averaging −28.5 ± 0.3‰ in δ¹³C (range = −28.9‰ to −28.1‰) were similar in their isotopic composition as uncarbonized grains (Table 1, Figure 4). Panicle B grains carbonized at 230°C (n = 8) also yielded similar carbon isotope values as uncarbonized grains from the same panicle, with an average of −29.0 ± 0.1‰ (−29.4 to −28.6‰). Panicle B grains carbonized at 300°C (n = 8) also yielded similar carbon isotope values as uncarbonized grains averaging −29.0 ± 0.1‰ in δ¹³C (range = −29.7 to −28.5‰). In nitrogen isotopes, grains from Panicle B carbonized at 230°C (n = 8), seeds averaged 1.2 ± 0.2‰ in δ¹⁵N (range =1.0 to 1.4‰), similar to uncharred grains from the same panicle. Panicle B seeds carbonized at 300°C (n = 8) exhibited on average 0.5‰ in δ¹⁵N relative to uncarbonized seeds from the same panicle with an average of 1.5 ± 0.2‰ (range = 1.2 to 1.8‰).

Grains carbonized at 400°C (n = 4) from Panicle A exhibited strongly different carbon and nitrogen isotope values with two specimens ashed during the firing process. The grains (n = 2) yielded an average δ¹³C of −26.0 ± 0.6‰ (range = −26.6 to −25.4‰) and an average δ¹⁵N of 4.0 ± 1.2‰ (range = 3.1 to 4.8‰).

Steinzeitpark barley from Panicles A and B carbonized at 230°C and 300°C showed overlapping ranges of carbon and nitrogen isotope variation with the corresponding uncarbonized grains from those panicles. Inter-panicle differences in δ¹³C were low, but more pronounced in δ¹⁵N regardless of charring condition; this is particularly clear in the diverging nitrogen values of Panicle E (Table 1, Figure 4).

Single-grain barley from Harappa Trench 42 yielded extremely wide variation in single-grain isotope values but notable overlap with isotope values measured from aggregate samples (Table 2, Figure 5). Across the 46 individually sampled grains, one sample failed due to low %N (Supplementary material 1). Single grain measurements (n = 35) averaged 18.3 ± 0.8‰ in Δ¹³C (range = 16.7 to 20‰) and 5.2 ± 0.6‰ in δ¹⁵N. Two grains yielded outstandingly high δ¹⁵N values of 16.0 and 16.4‰, respectively; removal of these outliers results for an average δ¹⁵N value of 4.6 ± 1.9‰ for single grains (range = 0.6 to 8.1‰; Table 1, Figure 5).

S1 grains show average Δ¹³C values of 18.3 ± 0.8‰ (n = 9; range =16.7 to 19.4‰) and δ¹⁵N of 6.8 ± 4‰. The corresponding aggregate sample A1 (grain n = 10) shows a mean Δ¹³C value of 18.8 ± 0.1‰, and δ¹⁵N of 4.5 ± 0.2‰, enriched on average 2.3‰ in ¹⁵N relative to the S1 seeds. The average carbon discrimination values of S2 and A2 are effectively identical, with slight differences between their average nitrogen isotope values. S2 (n = 26) grains average Δ¹³C value is 19.0 ± 0.7‰ (range =17.4 to 20‰) and S2 mean δ¹⁵N values are 4.6 ± 2.9‰ (range = 0.6 to 16.4‰). The corresponding aggregate sample A2 (grain n = 26) yielded a mean Δ¹³C of 19.0 ± 0.03‰ and 4.3 ± 0.2‰ in δ¹⁵N. There are no significant differences between the single grain or aggregate samples in either Δ¹³C [Kruskal-Wallis: H (1, 43) = 7.61, p = 0.643] or δ¹⁵N values [Kruskal-Wallis: H (1, 43) = 3.59, p = 0.742] (Table 2, Figure 6).

The Steinzeitpark Dithmarschen (SD) barley exhibited low variation in both δ¹³C and δ¹⁵N values as expected for C₃ crops grown in a temperate environment with ample water and no manuring inputs (Figures 2, 3). For C₃ grasses such as wheat and barley, grain δ¹³C values reflect water availability during crop grain-filling, whether mediated through precipitation, environmental humidity, soil moisture content, or irrigation (Araus et al., 1997a; Cappers and Neef, 2012; Farquhar et al., 1989; Tieszen, 1991; Vogel, 1993). The average δ¹³C value of Steinzeitpark uncarbonized grains (−28.9 ± 0.3‰) corresponds to crops grown in temperate conditions with ample water (“well-watered” = ≤ −26‰; Kohn, 2010; Wallace et al., 2013). There is also low carbon isotope variation (ca. 0.5‰) between inter-panicle mean δ¹³C values. Panicle D and E grains exhibit slightly greater ranges in δ¹³C values than uncarbonized grains from A and B, which is likely due to differences in the sample sizes of uncarbonized grains (Table 1).

The SD single-grain δ¹³C values display extremely low variation across all four panicles. The limited carbon isotope variation displayed in naked barley is consistent with previous studies examining single-grain intra-panicle variation in temperate grown wheat, which show ca. ± 0.5‰ variation in grains from a single panicle and grains between panicles (Heaton et al., 2009). This homogeneity is also consistent with studies showing ca. ± 0.5‰ δ¹³C variation between aggregate samples comprised of grains from crop plants from the same field (Nitsch et al., 2015; Wallace et al., 2015; Jones et al., 2021).

Crop nitrogen stable isotope values are influenced by soil conditions, aridity, and manuring intensity. The frameworks used for interpreting crop δ¹⁵N values derive from growing experiments where known amounts of manure were applied to fields in arid, semi-arid, and temperate environments δ¹⁵N values in England, Germany, Morocco, and Syria. Rather than correlating specific amounts of manure to grain nitrogen values, these frameworks rely on estimated benchmarks that gauge the relationship between grain δ¹⁵N values and manuring practices. Unmanured δ¹⁵N values fall below ≤ 2.5‰, moderately manured crops range between 2.5 and 6.0‰, and heavily manured crops are ~≥6.0‰ (Bogaard et al., 2007; Fraser et al., 2011, 2013; Kanstrup et al., 2014; Styring et al., 2016a).

The absence of fertilizer in the SD experimental plot in the form of manure or other amendments means bioavailable nitrogen in the low-nutrient sandy soils are driving plant tissue δ¹⁵N values (Craine et al., 2015; Szpak, 2014; Denk et al., 2017; Larsson et al., 2024). The low average δ¹⁵N value of 0.4 ± 0.6‰ for all uncarbonized grain from the SD plot reflects the lack of direct soil amendments through manuring and low bioavailable soil N (unmanured δ¹⁵N = ≤ 2.5‰; Fraser et al., 2011; Styring et al., 2016a; Larsson et al., 2019). This may be due to increased nitrate leaching in the sandy soils, and result in soil nitrogen depleted in ¹⁵N (Craine et al., 2015; Takebayashi et al., 2010), possibly leading to the low nitrogen isotope values seen here. The uncarbonized grain mean δ¹⁵N values from Panicles A and B exhibit greater similarity (A, n = 4, 1.2 ± 0.1‰; B n = 4, 1.0 ± 0.4‰) to each other than with Panicles D (n = 30, 0.7 ± 1.8‰) or E (n = 33, 0.2 ± 0.4‰) from which a larger number of grains were measured. Grain δ¹⁵N values from Panicle D and E display slightly more variation compared to A and B. In particular, Panicle E displays the lowest mean δ¹⁵N values, with 16 grains falling below −0.2‰ (Figure 4). The lower nitrogen isotope values of Panicle D and E may be due to slight spatial variation in soil conditions across the cultivation plot, with D and E root systems reflecting particularly low bioavailable nitrogen accessible by each plants root system (Amundson et al., 2003; Craine et al., 2009, 2015). The nitrogen isotope values displayed by Panicle D and E barley notably overlap with archaeological barley δ¹⁵N values from the Funnel Beaker Neolithic site Oldenburg LA77, located ca. 100 km away, with δ¹⁵N values of 0.1‰ (#264 grain n = 10, δ¹⁵N 0.1‰, #10284 grain n = 15, δ¹⁵N 0.6‰; Filipović et al., 2019).

These results underscore the minimal variation in carbon and nitrogen isotope ratios among grains from the same field. Aggregate samples of grain from arid rainfed fields in Jordan, Syria, and Spain have demonstrated an ± 1‰ difference in aggregate grain δ¹³C values from the same field but harvested in different years, attributed to annual differences in precipitation in arid environments (Flohr et al., 2011, 2019; Wallace et al., 2013). These studies highlight that intra-seasonal variation in precipitation, especially in semi-arid environments, in conjunction with anthropogenic watering practices, significantly influence carbon isotopic variation in charred seeds. In addition, they demonstrate how archaeological grains from contexts that were formed over multiple years, if combined into an aggregate sample, may yield a single carbon isotope ratio that erases key information of annual precipitation or water availability. Nitrogen isotope variation between unmanured grains harvested from single fields in England, Germany, Sweden, and Morocco averages ca. ± 1‰, and increases to ca. ± 2‰ in both manured single-grain (Bogaard et al., 2007; Larsson et al., 2019), and manured aggregate samples (Kanstrup et al., 2011; Fraser et al., 2013).

Archaeological studies using single grains have interpreted differences in δ¹⁵N values of greater than ≥2‰ between depositional contexts as either grains from crop plants that experienced different growing conditions in different fields (Lightfoot and Stevens, 2012; Larsson et al., 2019), or disparate access to and application of manure (Larsson et al., 2019). The consistent intra-panicle and inter-panicle isotopic differences found in grains from the same field, as shown here and in other studies, demonstrate the interpretive promise of single-grain isotope analysis in archaeological contexts.

The minimal impact of charring on the carbon and nitrogen isotopic composition of modern seeds from SD confirms previous research demonstrating that charring between 200 and 300°C leaves grains morphology intact and results in a slight ca. 0.1‰ enrichment in ¹³C and ca. 0.3‰ in ¹⁵N (Fraser et al., 2013; Charles et al., 2015; Nitsch et al., 2015; Stroud et al., 2023a). Grains carbonized at 230 and 300°C from Panicles A and B exhibit nominally more variation in their nitrogen isotopic compositions than uncarbonized grains (Table 1, Figure 4). There were slight incremental increases in the mean δ¹⁵N values of grains carbonized at 230 and 300°C in both Panicle A (230°C, n = 4, 1.0 ± 0.2‰; 300°C n = 8, 1.5 ± 0.2‰) and Panicle B (230°C, n = 8, 1.2 ± 0.2‰; 300°C, n = 8, 1.5 ± 0.2‰), but no changes in the overall variation (Table 1, Figure 4). In the four grains from Panicle A carbonized at 400°C, only two yielded enough material unrecognizable as grains, for analysis, and these specimens displayed clear deviations from the mean δ¹³C values (n = 2, 26.0 ± 0.8‰) and δ¹⁵N (n = 2, 4.0 ± 1.2‰). Critically, while there are indications of slightly increased nitrogen isotope variation between grains at higher temperatures, these differences are less than the limited variation displayed between different panicles (Table 1).

The aggregate and single-grain data from Trench 42 provides two slightly different interpretations of cultivation practices at Harappa ca. 1900 BCE. The high mean aggregate Δ¹³C values (A1 18.8 ± 1.5‰; A2 18.9 ± 1.5‰) for Harappa seeds fall within expectations for well-watered barley (≥ 18.5‰; Jones et al., 2021; Styring et al., 2016a; Wallace et al., 2013). Single-grain Δ¹³C values show only moderate variation (n = 35, 18.8 ± 0.8‰, range = 16.7 to 20.0‰), with most falling within the approximate range corresponding to medium-to-well-watered barley (Δ¹³C ≥ 17‰), with the exception of a single grain falling into the narrow “medium-to-low” (Δ¹³C 16 to 17‰) watering threshold (Wallace et al., 2013; Flohr et al., 2019; Jones et al., 2021). These values are similar to δ¹³C values of aggregated modern barley from Northwest India grown in flooded fields (n = 101, 18.0 ± 1.5‰), as opposed to rainfed (n = 4, 17.2 ± 2.8‰) or sprinkler-irrigated fields (n = 19, 16.8 ± 1.4‰; Jones et al., 2021). The high mean Δ¹³C values exhibited by single-grain samples are accompanied by relatively limited carbon isotopic variation, indicating that not only were many of the grains from well-watered fields, but also that this water availability was consistently managed (i.e., Styring et al., 2016a). This suggests that at Harappa, water access during the winter growing season was not restricted.

The aggregate mean δ¹⁵N values (A1: 4.5 ± 0.2‰; A2: 4.3 ± 0.2‰) fall within the benchmarks of low to medium levels of manuring (2.5 to 6.0‰; Bogaard et al., 2007; Fraser et al., 2011; Styring et al., 2016a), implying that the agricultural soils around Harappa received some manuring amendments, but not in large quantities. In contrast, single-grain δ¹⁵N values exhibit wide variation (n = 35, 5.2 ± 3.3‰, 0.55 to 16.38‰). Grain nitrogen isotope values are distributed across unmanured (n = 5, ≤ 2.5‰), low-to-medium (n = 19, 2.5 to 6.0‰), and heavily-manured (n = 11 ≥ 6.0‰) thresholds (Figure 5, Supplementary material; Fraser et al., 2011, 2013; Larsson et al., 2024; Styring et al., 2016a). While many of the grains span low to moderate manuring conditions, there is also the presence of some grains indicating no manuring, whereas a substantial number of grains indicate highly manured fields. This suggests either the ranges of δ¹⁵N variation from a single field are greater than might be expected, or that this storage feature held grains from multiple fields subjected to varying manuring practices or growing conditions.

The relationships between crop stable carbon and nitrogen values and agricultural organization hinges on the correlation of manuring and watering inputs to varying labor inputs and land use (Bogaard, 2004; Styring et al., 2017a; Bogaard et al., 2019). The interpretative framework for archaeobotanical stable isotope data guiding this correlation relies on models of labor- and land-limited agriculture derived from ethnographic studies of traditional smallholder farmers in Africa and Europe (Netting, 1993; Halstead, 2014). Labor-limited intensive practices are focused on small households dependent solely on scarce household labor for field management (Bogaard, 2004, 2005; Halstead, 2014). The expectation is that an increase in production under labor-limited conditions would entail high watering and high manuring inputs concentrated into small plots by nuclear households (Bogaard, 2005). These agricultural practices would be reflected in the isotopic composition of cultivated plants visible as high crop Δ¹³C and δ¹⁵N values (Wallace et al., 2013; Styring et al., 2016a). Land-limited extensive practices are associated with the use of plow agriculture, increasing available labor through traction to bring more land under cultivation, but the attendant spatial extensive land use means there are fewer watering, manuring, and labor inputs per field (Styring et al., 2017a). These less intensive applications of water and manuring would translate to a decrease in Δ¹³C and δ¹⁵N values (Styring et al., 2016a).

This framework traces connections between carbon and nitrogen isotope values that may indicate changing agricultural production encompassing a shift from intensive to extensive practices (Bogaard et al., 2013, 2019; Styring et al., 2017b,a; Yang et al., 2022). To date, exploration of land vs. labor-limited agriculture in ancient cultivation systems has relied on isotope values derived from aggregate grain samples. As archaeological specimens could span multiple growing seasons and thus different growing conditions (Riehl, 2020), aggregate samples containing large numbers of seeds attempt to ensure that samples encompass the entire range of variation of growing conditions in a field (Nitsch et al., 2015; Bogaard et al., 2016). Thus, aggregate studies attempt to endeavor to render comparison viable by averaging the isotopic information of all potential cultivation practices or growing conditions represented at a given site (Kanstrup et al., 2012; Styring et al., 2017a).

The Trench 42 aggregate samples (A1 and A2) mean Δ¹³C values correspond to those expected for moderately to well-watered crops (Jones et al., 2021), while the mean δ¹⁵N values fall within the benchmarks of medium to potentially low levels of manuring (Bogaard et al., 2007; Fraser et al., 2011; Styring et al., 2016a). Within the framework of labor- and land-limited agriculture, in conjunction with the widespread evidence of traction and plow agriculture at Harappa, this data suggests limitations on available land as opposed to available labor. It may be that labor was abundant and the cultivation of additional fields gave rise to land becoming more scarce and increasingly valuable (Bogaard, 2005). Therefore, within this framework, the isotope values from aggregate samples broadly correspond to evidence for land-limited extensive agricultural practices at Harappa (Bogaard, 2004; Bogaard et al., 2019).

The overlap between the mean Δ¹³C and δ¹⁵N values of single-grain S1 and S2 and aggregates A1 and A2 indicate that an aggregated sample pool captures a representative signature of cultivation. However, the variation displayed in S2 compared to A2 illustrates how isotope values measured from aggregate samples do not fully describe the full range of isotopic variation expressed in archaeological grains. Along these lines, the carbon and nitrogen isotopic variation shown by individual barley grains (S1 and S2) from Trench 42 potentially complicates an interpretation of land-limited production (Figures 5, 6). The mean Δ¹³C and δ¹⁵N values produced by S1 and S2 are similar to A1 and A2, indicating that the aggregate samples capture an average isotopic variation from a primary context (Figure 5). However, the variation in all single-grain δ¹⁵N values from all Trench 42 barley (5.2 ± 3.3‰, range = 0.6 to 16.4‰; Figure 6), is greater than observed intra-field variation in modern manured barley single-grain δ¹⁵N values (8.9 ± 1.64‰, range = 6.4 to 11.8‰; Larsson et al., 2019). At the same time, the limited variation in single-grain Δ¹³C values indicates the consistent watering of cultivars by farmers, or environmental conditions leading to well-watered barley (Δ¹³C 18.8 ± 0.8‰).

This pattern of constrained Δ¹³C values and wide δ¹⁵N variation at Harappa might be explained in two ways. First, grains from moderately manured fields may express a wider variation in δ¹⁵N values than so far documented in experimental plots and archaeological studies, as the irregular distribution of manure and patchy release of nitrogen across a field may be responsible for increased δ¹⁵N variation in manured plants (Larsson et al., 2019; Styring et al., 2024). If this is the case, S1 and S2 could represent a single well-watered field, but with varying application of soil amendments leading to wide nitrogen isotope variation in measured grains. Alternatively, this pattern in barley cultivation might indicate active water management which was accompanied with both low- and high-input manuring. Wide ranges of carbon and nitrogen isotope variation resulting from different field conditions might be expected from extensive agricultural practices that bring more land under cultivation (Styring et al., 2017a).

These results underscore the need to further evaluate the isotopic variation expressed in single grains from crops grown in intensive and extensive agricultural systems In particular, additional research exploring single-grain δ¹⁵N values for both low manure input fields and intensively manured fields is needed to better characterize how these differences in practice might impact the amplitude of nitrogen isotopic variation observed in charred seeds. What is clear is that single-grain isotope analysis increases the visibility of environmental and anthropogenic sources of uncertainty embedded in archaeobotanical isotope data and has additional implications for sample selection that should be investigated further (Calin-Jageman and Cumming, 2019). It potentially allows for the detection of smaller isotopic effect sizes (i.e., meaningful variation) that might be lost in an aggregate sample. In contrast, aggregate samples, while containing large numbers of grains, effectively create a single average with fewer statistical degrees of freedom, in turn increasing potential uncertainty, and lesser statistical precision. The use of single grains could in turn allow better quantification and interpretation of this uncertainty (Styring et al., 2024).