C₃ plant isotopic variability in a boreal mixed woodland: implications for bison and other herbivores

Study location: Elk Island National Park, Alberta

Elk Island National Park (EINP) is a ~200 km² protected area located ~40 km east of Edmonton. The park is located within Canada’s southern boreal plains ecozone, an area of transition between semi-arid prairie and wetter boreal forest (Nicholson, 1995). Topographically, the park is part of the Beaver Hills region, an area of knob-and-kettle terrain with abundant lakes and wetlands. Vegetation within the park is a patchy mosaic of aspen parkland, boreal mixed woodland, grassy/shrub meadows, marshes, and lacustrine areas (Fig. 1) (Best & Bork, 2004; Holsworth, 1960; Nicholson, 1995). All identified plant taxa in the park utilize C₃ photosynthesis (Hanna Schoenberg, personal communication, May 18, 2021). EINP’s mean annual temperature was 1.7 °C and mean annual precipitation was 460 mm between 1951 and 1980, but both temperature and precipitation have been increasing due to climate change (https://climatedata.ca/). EINP typically experiences moderate summers and cold, dry, windy winters. Temperatures range from average lows of −18 °C in January to average highs of 23 °C in July (weather-atlas.com). Peak summer rains occur in July (mean of 112 mm precipitation) and snowfall reaches a high of 206 mm in March (weather-atlas.com). Spring blooms typically begin to appear in May and the growing season lasts from approximately mid-May to mid-September.

EINP is home to several large ungulate species, including moose (Alces alces), elk/wapiti (Cervus canadensis), white-tailed deer (Odocoileus virginianus), mule deer (Odocoileus hemionus), plains bison (Bison bison bison), and wood bison (Bison bison athabascae) (Telfer & Cairns, 1986). For many decades, EINP has been a source for genetically-pure disease-free bison that have been introduced to conservation herds across the continent (Markewicz, 2017). The plains and wood bison areas are separate: plains bison range freely within the fully-fenced northern portion of the park and wood bison range freely within the separate, fully-fenced southern portion of the park (Fig. 1). Bison in both areas have access to the same types of habitat and vegetation.

Sample collection and preparation

Plant samples were collected with the permission of Parks Canada (Research and Collection Permit EI-2016-21863). Grass, sedge, forb, shrub, and horsetail (Equisetum spp.) samples were collected on June 27–29, 2016 (n = 133) and November 6, 2016 (n = 32) from dry open areas (e.g., meadows, hill slopes), dry closed-canopy environments, and wet areas (shorelines of ponds or lakes) (Figs. 1, 2). Site categorizations were based on observations at the time of sampling rather than on generalized vegetation maps, because wetlands can be ephemeral. Sampling sites were selected based on recent sightings of bison and physical evidence of bison (e.g., dung, wallows, hoofprints) in the area. To mimic bison foraging patterns, only terrestrial above-ground plant parts were collected. For the same reason, graminoids were prioritized for collection and sampled relative to their abundance at each location. Since only two sampling sites were wetlands and sedge cover was sparse, the sample size for sedges is low. Plants were identified to genus or species with reference to Johnson et al. (1995).

All samples were air-dried and ground to a fine powder with a Wig-L-Bug device prior to isotopic analysis. Most of the samples (n = 131) were homogenized into ‘whole plant’ samples such as might be consumed by a relatively indiscriminate herbivore, including varying proportions of leaves, stems, seeds, and/or flowers (Table 2). To test for isotopic differences among plant parts, leaves and seeds/flowers were analyzed separately for a subset of samples (n = 34). Grass leaves are wrapped around stems before diverging as a separate blade, making stems and leaves difficult if not impossible to separate in bulk samples. Grass flowers are complex structures that include a rachis and many tiny pedicels which are likewise difficult or impossible to separate from the floret. As a result, grass leaf and seed/flower samples include variable proportions of these other tissues as well.

Carbon and nitrogen isotope measurements

Carbon and nitrogen isotope values (δ¹³C, δ¹⁵N) and carbon and nitrogen contents (dry weight %C, %N) were obtained using an Elementar VarioMicro Cube elemental analyzer coupled with an Isoprime isotope-ratio mass spectrometer in continuous-flow mode. Carbon and nitrogen isotope values were obtained during the same run by combusting approximately 1 mg of sample and using a high level of dilution to reduce the carbon dioxide gas peaks. Nitrogen isotope results from samples with N₂ peaks <1 nA were excluded unless duplicate analyses exhibited similar reproducibility to samples with larger gas peaks. Carbon isotope values of the low-nitrogen samples were retained since the carbon peaks were large enough to produce reliable results. The samples with low nitrogen-gas peaks are those lacking δ¹⁵N values in Table 2.

δ¹³C values were calibrated to VPDB and δ¹⁵N values were calibrated to AIR using USGS-40 and USGS-41 or 41a (accepted δ¹³C values of −26.39, +37.63 and +36.55‰ and accepted δ¹⁵N values of −4.52, +47.57 and +47.55‰, respectively). Sample replicates (minimum 10% of samples in each run) and internal check standards of methionine, amaranth, and red lentil (long-term mean δ¹³C of −28.60, −13.59, −26.12‰; long-term mean δ¹⁵N of −5.04, +2.94 and −1.09‰, respectively) were used to monitor measurement uncertainty. Uncertainty measures were calculated following the method of Szpak, Metcalfe & Macdonald (2017). For δ¹³C, precision u(R_w) was 0.11‰, accuracy (u(bias)) was 0.09‰, and total analytical uncertainty (u_c) was 0.14‰. For δ¹⁵N, precision was 0.23‰, accuracy was 0.23‰, and total analytical uncertainty was 0.33‰.

Statistical analyses

Statistical analyses were conducted using Excel for Office 365 and PAST (PAleontological STatistics) 4.03. Shapiro–Wilk W tests were used to assess the normality of distributions. Levene tests were used to evaluate the homogeneity of variance. Normally distributed datasets (carbon isotope values) were compared using Student’s t-tests (two independent samples), paired-sample t-tests (two paired samples), or one-way ANOVA F-tests with Tukey’s post-hoc comparisons (three or more independent samples). Non-normally distributed datasets (nitrogen isotope values) were compared using Mann–Whitney U tests (two independent samples), Wilcoxon sign-rank tests (two paired samples) or Kruskal–Wallis H tests with Dunn-Bonferroni post-hoc comparisons (three or more independent samples). Alpha was set to 0.05 for all statistical comparisons. In the text below, means are reported with standard deviations, unless noted otherwise.

Whole sample

Plant δ¹³C values ranged from −32.6 to −24.9‰, with a mean and standard deviation of −28.5 ± 1.5‰ (Tables 2, 3). Plant δ¹⁵N values ranged from −3.9 to +9.9‰, with a mean and standard deviation of +0.4 ± 2.7‰. The shape of the distribution was normal for δ¹³C (Shapiro−Wilk W = 0.99, n = 165, p = 0.7; skewness = −0.05) and positively skewed for δ¹⁵N (Shapiro–Wilk W = 0.92, n = 141, p < 0.001; skewness = 1.14) (Fig. 3).

Plant types

The mean δ¹³C values of grasses, sedges, shrubs, forbs and horsetail were within 1.9‰ of one another (Table 3), and an ANOVA showed no statistically significant differences among the groups (F(4,160) = 1.3, p = 0.28). With horsetail removed (because of its small sample size), there were still no significant differences in δ¹³C among grasses, sedges, shrubs, and forbs (F(3,158) = 1.0, p = 0.39) There was a significant difference among the δ¹⁵N values of plant types (H(4) = 12.9, p = 0.01), but the Dunn–Bonferroni test suggested that only the horsetail-forb comparison was significant (p = 0.03). With horsetails removed there was no statistically significant difference among grasses, sedges, shrubs, and forbs (H(3) = 7.0, p = 0.07), and their medians were within 2.3‰ of one another. Although the median grass δ¹⁵N value did not significantly differ from that of any other group, grasses had the greatest variability of any plant type, and grass samples had both the highest (>+5.1‰) and lowest (<−2.3‰) individual plant δ¹⁵N values (Table 3, Fig. 4). A Levene’s test from medians (i.e., Brown–Forsythe test) indicated that the difference in the variability of δ¹⁵N among plant types was statistically significant (p = 0.01).

Habitats

Plant growth habitat had a significant effect on the carbon isotope compositions of plants (F(2,162) = 48.8, p < 0.001). The differences among all three groups were statistically significant, with the highest δ¹³C values in open areas (−27.9 ± 1.2‰, n = 108), intermediate values in wet areas (−28.9 ± 1.4‰, n = 12) and the lowest values in closed-canopy areas (−30.0 ± 1.1‰, n = 45) (Table 3, Fig. 5). Growth habitat also affected δ¹⁵N values (H(2) = 7.7, p = 0.02), with higher δ¹⁵N values in wet habitats (+2.6 ± 2.7‰, n = 10) compared to those in either open areas (+0.1 ± 2.4‰, n = 94) or closed canopy areas (+0.5 ± 3.1‰, n = 37). Although wet areas had higher mean (and median) δ¹⁵N values than the open or closed-canopy areas, the latter two habitat types hosted the plants with the highest individual δ¹⁵N measurements (Fig. 5). As mentioned previously, these extreme δ¹⁵N values were all from grass samples. There was a positive skew in the δ¹⁵N values of plants from open environments (W = 0.9, n = 94, p < 0.001) and closed environments (W = 0.9, n = 37, p < 0.001).

Plant parts

Carbon isotope compositions of leaves were on average 1.2‰ lower than those of seeds/flowers from the same plants (paired samples t = 7.8, df = 33, p < 0.001) (Table 3). Furthermore, the great majority of individual plant samples had lower leaf than seed/flower δ¹³C values, with seed/flower minus leaf differences (Δ¹³C_seed–leaf) of individual plants ranging from −0.5 to +3.1‰ (Fig. 6). The lowest mean and individual δ¹³C values were obtained from leaves in closed habitats, and the highest mean δ¹³C from seeds in open habitats (Fig. 6).

Nitrogen isotope compositions of leaves were 0.5‰ lower on average than those of seeds/flowers from the same plants (Table 3), but the difference was not statistically significant (Wilcoxon W = 250, df = 27, p = 0.06). Individual plants had highly variable seed-minus-leaf differences (Δ¹⁵N_seed–leaf), ranging from −2.4 to +2.9‰ (Fig. 7).

Seasonal changes

Seasonal shifts in plant δ¹³C and δ¹⁵N occurred between early summer (late June) and mid fall (early November) (Table 3, Fig. 7). Plant δ¹³C increased slightly during fall, both for the whole dataset (t(163) = 2.1, p = 0.04, mean difference of 0.6‰) and when only locations sampled in both seasons were included (t(62) = 2.2, p = 0.03; mean difference of 1.0‰). Plant δ¹⁵N also increased during fall, both for the whole dataset (U = 582.5, p = 0.003; mean difference of 2.5‰) and when only samples from matched locations were compared (U = 145, p = 0.02; mean difference of 2.0‰). Plant nitrogen contents (%N) also significantly decreased from summer to fall (whole sample: U = 265, p < 0.001; mean difference of 0.8% matched locations: U = 71, p < 0.001; mean difference of 0.7%) (Fig. 7). The true seasonal decrease in plant nitrogen content is likely greater than this value implies, since proportionally more fall plant samples were excluded due to their small gas peaks (see Table 3).

Plant isotopic distributions

The distribution of plant δ¹³C values was normal. The EINP whole-sample mean δ¹³C of −28.5‰ is somewhat lower than the modern global mean C₃ plant δ¹³C value of −27.0‰ determined by Kohn (2010). This can be attributed to two main factors: (1) the δ¹³C of atmospheric CO₂ during our sample collection (in 2016) was significantly lower than Kohn’s (2010) normalized value of −8.0‰ because of the ongoing effects of fossil fuel burning (Long et al., 2005), and (2) Kohn’s (2010) study excluded understory plants with δ¹³C values below −31.5‰, whereas the present study did not.

Distributions of plant nitrogen isotope compositions were positively skewed. Skewness of isotopic distributions is seldom explicitly evaluated, and isotopic data presentations that facilitate visual examination of skewness (e.g., frequency histograms, box-and-whisker plots) are relatively rare, so it is difficult to determine how common skewed plant nitrogen isotope distributions may be. Metcalfe & Mead (2019) observed a negatively skewed δ¹⁵N distribution for Pleistocene plants. Funck et al. (2020: Supplemental Material) provide box-plots that appear to illustrate positively skewed modern grass δ¹⁵N and negatively skewed modern herb δ¹⁵N distributions, but they did not explicitly evaluate skewness. The other plant isotopic studies reviewed here neither evaluated skewness nor presented data in forms that make it easy for readers to evaluate. Determining the shape of a distribution is often overlooked but testing for normality is a critical first step before utilizing parametric statistical methods, at least when sample sizes are small (which is typical in most archaeological and paleontological studies) (Ghasemi & Zahediasl, 2012). Failing to recognize skewed isotopic distributions can result in inappropriate data reporting (e.g., use of means and standard deviations rather than medians and interquartile ranges) and use of statistical tests whose assumptions are not met (i.e., parametric statistical tests), potentially producing invalid results and leading to erroneous interpretations. Assessing the skewness of dietary components (and other characteristics of data distribution) is also critical for studies using stable isotope mixing models, which typically assume normal distributions and require dietary inputs of means and standard deviations (Cheung & Szpak, 2020).

It is possible that skewed plant δ¹⁵N distributions could help explain the strong nitrogen isotopic niche partitioning that has been observed among herbivores in some ecosystems. In particular, mammoths tend to have significantly higher δ¹⁵N values than co-existing herbivores, which is related to a dietary (rather than physiological) difference (Schwartz-Narbonne et al., 2015). In the present study, grasses had the greatest variability in δ¹⁵N of any plant taxon and all of the most positive δ¹⁵N values in the skewed tail of the distribution (i.e., values >5.1‰) were from grasses (Table 3, Fig. 4). Grasses are the predominant food of mammoths, but also of bison, who do not have enriched δ¹⁵N values. If variables could be identified that predict which grass specimens within a given ecosystem have high δ¹⁵N values (i.e., taxa, parts, growth-stages, growth habitats), then it might be possible to determine if mammoths were likely to have been selecting such grasses (for example, by employing different feeding strategies or preferring different microhabitats). In general, a herbivore preferentially selecting plants from the skewed ‘tail’ of an isotopic distribution would be predicted to occupy a distinct isotopic niche relative to herbivores that are randomly selecting plants from throughout the distribution. This would also be true of herbivores selecting plants whose δ-values fall within the tails of a normal distribution, but a skewed plant distribution would be predicted to result in greater herbivore isotopic niche differentiation due to the more extreme values of outliers in the skewed tail of the distribution.

Plant types

The overlapping δ¹³C and δ¹⁵N values of grasses, sedges, forbs, and shrubs in EINP highlights the importance of understanding local plant variability when interpreting herbivore isotopic compositions. Previous research has established some generalities about isotopic differences among primary producers. For example, lichens often have higher δ¹³C values than terrestrial plants (e.g., Brooks et al., 1997; Teeri, 1981), woody gymnosperms generally have higher δ¹³C values than woody angiosperms (Hare & Lavergne, 2021), and aquatic plants tend to have higher δ¹⁵N values than terrestrial plants (e.g., Kielland, 2001; Plint, Longstaffe & Zazula, 2019. However, comparisons of differences among plant types at local levels can produce disparate results (e.g., Drucker et al., 2010: Fig. 4), which is perhaps not surprising when one considers the complex range of environmental factors that affect δ¹³C and δ¹⁵N, as well as the fact that researchers select different plant groups for study and even categorize them differently (Table 1). Compilations of isotopic data from plants growing in various habitats (i.e., global or regional datasets) can obscure the effects of microhabitats (e.g., degree of canopy cover, altitude, aridity, etc.), which may be more important variables for interpreting herbivore isotopic compositions than plant type. Studies that compare herbivore isotopic compositions in ancient C₃ ecosystems to a plant baseline organized by plant type (e.g., Schwartz-Narbonne et al., 2019; Schwartz-Narbonne et al., 2021) presuppose that type is the most important predictor of a plant’s isotopic compositions. An alternative approach is to put equal or greater emphasis on major environmental factors that influence plant isotopic compositions, such as the canopy effect (e.g., Drucker et al., 2008; Hofman-Kamińska et al., 2018) or ecosystem changes (e.g., Drucker et al., 2011; Metcalfe & Longstaffe, 2014).

The plant-type data in the present study highlight the importance of ecosystem-specific contexts. In particular, it is not appropriate to assume that grasses, sedges, shrubs and forbs will have consistent relative isotopic differences in disparate environments and temporal intervals. Consequently, similar isotopic niches among animals does not necessarily mean that animals ate the same things, or that “one species could fulfill another’s ecological role” (Schwartz-Narbonne et al., 2019: 1). Rather, isotopic niche overlap could simply indicate that there are minimal isotopic differences among the disparate plants consumed by herbivores in that environment. Likewise, minimal isotopic variations in serially-sampled animal tissues does not necessarily suggest that animals had highly specialized diets with minimal variation. On the contrary, minimal seasonal isotopic variations in herbivore tissues could occur even when animals undertake significant seasonal changes in diet if there are no significant differences among plant types in that area. Given these complexities, the key to being able to make meaningful interpretations of herbivore isotopic compositions is to have a good understanding of which isotopic baselines and variables are most important for a particular study and to seriously consider alternative interpretations based on the various factors that can influence isotopic systems. In general, isotopic niches are far from equivalent to dietary niches or dietary specializations.

Plant parts and habitats: carbon isotopes

Lower plant δ¹³C values in EINP closed habitats compared to open habitats (~2‰ on average) is consistent with the well-known canopy effect, in which understory plants have significantly lower δ¹³C values than plants that make up the canopy or emergent layers, or plants that grow in open areas (e.g., Bonafini et al., 2013; Chevillat et al., 2005; Drucker et al., 2008; Van Der Merwe & Medina, 1991). The lower δ¹³C values in EINP leaves relative to seeds/flowers (~1‰ on average) is likewise in agreement with the 1–3‰ difference that has been reported in many other studies (e.g., Badeck et al., 2005; Ghashghaie & Badeck, 2014; Metcalfe & Mead, 2019).

The EINP plant isotopic data suggest that among herbivores, a combined effect of plant-part and habitat-selection could result in significant carbon isotope niche partitioning within C₃ environments, with the largest differences between animals consuming seedy/flowery plants in open environments (higher δ¹³C) and those selecting seedless/flowerless plants in closed environments (lower δ¹³C). This offers an alternative to assuming that animal niche partitioning in C₃ environments is due to differing proportions of grass vs browse or consumption of different plant taxa. Many previous studies have used herbivore δ¹³C to infer the ‘openness’ of utilized habitats (e.g., Bocherens et al., 2015; Doppler et al., 2017; Drucker et al., 2003; Drucker et al., 2011), but few have considered the additional isotopic effects of plant-part differences, such as the decrease in leaf δ¹³C than occurs as the leaf expands (Vogado et al., 2020) or differences among seedier versus seedless plant parts (but see Guiry et al., 2020 for an exception). The effects of ‘seedy’ vegetation on herbivore isotopic compositions deserves further study, since there may also be differential digestibility among seeds and leaves that influences their incorporation into herbivore tissues.

In general, herbivore feeding specializations go beyond selection of particular plant forms, species and habitats to include specialization on particular plant parts and growth stages. These differential feeding strategies might have particularly pronounced isotopic effects in an environment like the mammoth steppe, where co-existing grazers likely consumed different parts of the same plants. For example, elephantids rip out tall (potentially seedy) bunches of grasses by grabbing them with their trunks, whereas bison break off short (probably less seedy) grasses and tall/mid-level new growth with their tongues and teeth (Guthrie, 1982). On the mammoth steppe, bison tended to have higher δ¹³C values than mammoths in a range of locations and temporal intervals (e.g., Bocherens, 2015). Higher δ¹³C values in a taxon that consumes shorter grasses is the opposite of what would be expected if ‘seediness’ was a factor in isotopic niche differentiation. However, the higher δ¹³C values of bison could result from bison consuming a larger proportion of short, newly-grown leaves, which tend to have higher δ¹³C values than older mature leaves (Vogado et al., 2020). Regardless of what drives isotopic niche differentiation on the mammoth steppe, the results of the present study suggest that in some environments, habitat and plant-part selection could have greater isotopic effects on herbivore isotopic compositions than selection of different plant taxa.

Plant habitat: nitrogen isotopes

EINP plants from the wet habitat tended to have higher δ¹⁵N values than plants from the dry (open or closed-canopy) environments. Although this contrasts with the general trend towards higher δ¹⁵N values in drier locations that is often observed on regional and global scales (Craine et al., 2009; Handley et al., 1999; Wang et al., 2014), it is consistent with the higher plant δ¹⁵N values often observed in aquatic systems relative to terrestrial systems (e.g., Cloern, Canuel & Harris, 2002; Kielland, 2001; Plint, Longstaffe & Zazula, 2019. It is possible (and perhaps likely) that EINP terrestrial plants growing in seasonally wet areas obtained some nitrogen from aquatic sources, leading to higher δ¹⁵N values. It is also possible that herbivore dung is frequently deposited in wetland areas when animals come to drink, contributing ¹⁵N-enriched nitrogen to the wetland system and mimicking the established effects of manuring on plant δ¹⁵N (e.g., Bogaard et al., 2007; Szpak et al., 2014). It is important to note that the sample size available for EINP wetland habitats was small, so the reliability of this habitat difference should be re-examined in future studies. Nevertheless, in combination with previous studies that clearly show higher δ¹⁵N values among aquatic plants, these results suggest caution for archaeologists and paleoecologists who interpret higher herbivore δ¹⁵N as indicators of increased aridity. An alternative explanation (among others) for high herbivore δ¹⁵N values could be the consumption of plants growing in or near nutrient-rich wetlands.

Seasonal changes in plant isotopic compositions

A summer-to-fall (late June to early November) increase in both δ¹³C and δ¹⁵N (by ~1 and 2‰, respectively) was observed in EINP plants. This could be due to a combination of factors, including changes in the biochemical compositions of tissues, changes in source C and N isotopic compositions, remobilization of nutrients into roots for winter, and early decomposition. The direction and magnitude of seasonal isotopic changes in plants may vary among environments and locations. For example, Karlsson, Eckstein & Weih (2000) found that the δ¹⁵N values of most Subarctic plants in northern Sweden increased between the snowmelt (May) and mid-June, but decreased in August and September, with a range in seasonal variation of 2.1 to 5.3‰. On the other hand, the timing of key seasonal changes (e.g., temperature increases and decreases) varies considerably among locations and makes seasonal generalizations challenging.

Reconstructing ecosystem-specific seasonal changes in plant δ¹³C and δ¹⁵N could help researchers interpret serial-sampling studies of herbivore isotopic compositions, which may vary due to seasonal changes in diet, physiology, and/or isotopic variations in plants. Seasonal changes in the diets of a range of herbivores have been studied within C₃-dominated ecosystems, and these changes are often relatively small in magnitude (~2 to 3‰ or less). For example, Funck et al. (2020) observed temporal changes in sectioned wood bison (Bison bison athabascae) hair δ¹³C and δ¹⁵N that they attributed to nutritional stress. Julien et al. (2012) serially-sampled steppe bison (Bison priscus) teeth and interpreted small winter increases in δ¹³C as an indication of lichen consumption. Metcalfe & Longstaffe (2014) identified different seasonal patterns in the tooth enamel of mastodons (Mammut americanum) that lived in the same geographical area during different time periods, which they suggested were the result of major vegetational shifts. Kielland (2001) serially-sampled Alaskan moose (Alces alces) hooves and interpreted variations of about 2–3‰ as evidence for seasonal changes in diet. Plant isotopic values and their variability underlie the interpretations of all these studies.

Bison generally consume graminoids year-round but may seasonally switch between grasses and sedges, and/or consume forbs and woody plants when graminoids are not available (Gogan et al., 2010). The minimal isotopic differences among plant taxa in EINP suggests that these seasonal shifts in bison foraging strategies might not be recorded in the isotopic compositions of incrementally growing bison tissues such as teeth or hair. However, based on the EINP seasonal plant data, one might predict that seasonal isotopic shifts in the plants themselves could be recorded in serially-sampled bison tissues. Generalizing to other environments, researchers should be aware that seasonal changes in herbivore isotopic compositions do not necessarily indicate changes in foraging strategies, but can result from isotopic changes within the plants themselves.