Contrasting isotopic responses of dryland and wetland plants to a century of global anthropogenic changes in nutrient cycling

Anthropogenic emissions of carbon dioxide and reactive nitrogen in various forms disrupt the functioning of ecosystems around the world. In Europe, many valuable habitats, particularly wetlands and semi-natural dry grasslands, are under threat from ongoing eutrophication. However, due to contrasting water regimes, the uptake of anthropogenic nitrogen by plants in these ecosystems is different and is also interrelated with an increase in trophic level in both habitats.

In our study, we measured the δ¹⁵N and δ¹³C values, as well as the total nitrogen content (TN), of 99 pairs of foliar samples collected from 14 species of vascular plants in Poland: seven wetland species and seven dry grassland species. Each pair consisted of a historical sample collected from a herbarium voucher dating from before 1939 (i.e. before artificial fertilisers were widely used in agriculture) and a contemporary sample collected in 2024 from the same species in a similar location.

We performed t tests to determine whether there were significant differences in the means of δ¹⁵N, TN and δ¹³C between samples from the two different habitats. Next, we calculated the differences in δ¹⁵N, TN and δ¹³C between the contemporary and historical samples for each pair. Using the information on date of collection of each herbarium voucher we corrected δ¹³C values taking into consideration changes of atmospheric δ¹³C due to fossil-fuel combustion (so called Suess effect). We then tested whether the difference for each species and habitat type was significantly different from zero using 90 % confidence intervals. Using multiple linear regression, we analysed the relationships between differences in δ¹⁵N and TN over time and the following factors: habitat type, the proportion of farmland in the surrounding landscape, the consumption of synthetic nitrogen fertiliser, and NO_x deposition.

δ¹⁵N and TN values were lower in dry grassland species than in wetland species, in both the contemporary and historical subsets. For dry grassland species, the mean δ¹⁵N value was lower in contemporary samples than in historical ones. For wetland species, the opposite was true. The difference in δ¹⁵N values between pairs of samples was positively related to the amount of farmland in the surrounding landscape. The mean TN was higher in contemporary than in historical wetland plants, but not in dry grassland plants. The mean δ¹³C value, corrected for the Suess effect, was lower in contemporary samples than in historical ones. The mean difference was −0.51 ‰ for dry grassland and −3.85 ‰ for wetland species.

Our study revealed that the century of fossil fuel-derived carbon emissions, increased nitrogen input into the environment, and dominance of artificial fertilisers and combustion-derived nitrogen over biological nitrogen sources have not led to consistent responses across habitats and species. While the isotopic composition of nitrogen and carbon in plant leaves in Central Europe has undoubtedly changed, this change is highly context-dependent. Its magnitude and direction are impacted by the type of habitat and the identity and/or ecology of the species. As expected, man-made alterations appear to be more pronounced in wetland environments than in dryland habitats. Furthermore, the source of disruption may differ between the two habitat types. Specifically, wetlands are exposed to a multitude of anthropogenic nitrogen and carbon sources, whereas dry grasslands seem to be predominantly affected by changes in the composition of the atmosphere.

1 Introduction

Disruption to the global nitrogen and carbon cycles caused by human activities is one of the major threats to the functioning of the biosphere (Hessen et al., 1997; Richardson et al., 2023). The main source of anthropogenic carbon is the combustion of fossil fuels, which produces carbon dioxide (CO₂), one of the main greenhouse gases. However, the high-temperature combustion of fossil fuels in the atmosphere (e.g. in engines of cars and trucks and in electrical power generation plants) also produces nitrogen oxides (NO and NO₂), which play a role in the formation of both smog and acid rain. They also contribute to the atmospheric nitrogen deposition to the ecosystems, mostly after conversion to nitric acid HNO₃ (Seinfeld and Pandis, 2006; Hertel et al., 2012). Fossil fuel combustion began on a massive scale in the Industrial Era (Pirani, 2018). Later, during the 20th century and continuing to the present day, anthropogenic nitrogen input has increased, particularly due to the invention of an efficient method of synthesising ammonium through fixation of atmospheric nitrogen N₂ (Bateman and Kelly, 2007; Bodirsky et al., 2012). This process enabled the large-scale production and use of artificial fertilisers in agriculture, which was one of the foundations of the so-called Green Revolution of the 1960s (Bodirsky et al., 2012). However, excess of anthropogenic, bioavailable nitrogen input into natural ecosystems usually leads to profound changes in their productivity, biodiversity and functioning, and often cannot be easily reversed (Stevens et al., 2004; Dise et al., 2011; van Der Plas et al., 2024). Furthermore, the combined effect of increased CO₂ and reactive nitrogen may have a more negative impact on biodiversity and ecosystem functioning than the effects of each factor acting alone (Bakker et al., 2024; Reich et al., 2024).

An excess of nitrogen input is a well-recognised threat to a wide range of ecosystems that are more susceptible to eutrophication, particularly aquatic habitats (Schindler, 2006), but also species-rich semi-natural grasslands in Europe (Stevens et al., 2004; Bobbink et al., 2010; Dise et al., 2011). Eutrophication in aquatic ecosystems, mainly caused by nitrogen and phosphorus, leads to uncontrolled algal or cyanobacterial blooms. These blooms can result in episodes of hypoxia (low oxygen levels), reduced light penetration, the release of toxic compounds and unpleasant odours, increased water turbidity, and alterations to biocoenoses (Schindler, 2006; Erisman et al., 2008). These changes can include the extinction of underwater plants (Qin et al., 2013) and increased fish mortality (Schindler, 2006). In grassland ecosystems, high nitrogen input may lead to changes in nutrient and water cycling, soil acidification and increased biomass production (Bobbink et al., 2010; Meier et al., 2022). Higher productivity can lead to the dominance of fast-growing, expansive species, while the production of large amounts of litter can inhibit the growth of smaller, less competitive, and often rare species (Meier et al., 2022; Chen et al., 2023). Larger biomass production and litter accumulation also reduce the amount of light reaching ground level, which usually harms germination and seedling survival (Ruprecht et al., 2010; Meier et al., 2022).

It is estimated that in 2020, the total amount of nitrogen produced globally by humans was 226 Tg N yr⁻¹, while natural nitrogen fixation ranged from 58 to 128 Tg N yr⁻¹ (Galloway et al., 2021). It is estimated that global consumption of artificial nitrogen fertilisers increased nearly tenfold from 11.3 Tg N yr⁻¹ in 1961 to 108 Tg N yr⁻¹ in 2022 (Lu and Tian, 2017; FAO, 2024). The highest use of nitrogen-based fertilisers in agriculture was among the European Union's main agricultural producers, particularly in France (1.7 million t), Poland and Germany (1.0 million t each) (Eurostat, 2025). In Poland in 2020, the average amount of pure nitrogen used in artificial fertilisers was 70.4 kg per hectare of agricultural land (GUS, 2021). As a result, Poland can be classified as one of the countries with the highest artificial fertiliser consumption per hectare of agricultural land, with fertilisation levels exceeding the EU average by about one third (Piwowar, 2022). Additionally, around 22 kg of nitrogen per hectare of Polish agricultural land is added in the form of organic fertilisers both solid and liquid ones (GUS, 2026). It is estimated that only about 30 %–50 % of the nitrogen from the fertilisers applied is taken up or assimilated by crops, with the remainder being lost through leaching, surface runoff and gaseous emissions (Lehnert et al., 2018). In addition to fertilisation, important anthropogenic sources of nitrogen include high-temperature combustion (Seinfeld and Pandis, 2006; Hertel et al., 2012) and wastewater from households, industry, and agriculture. The latter is a source of mostly organic nitrogen that impacts aquatic environments (Chen et al., 2011).

The physiology and evolutionary strategy of nitrogen uptake by plants, such as through mycorrhizae or the presence of nitrogen-fixing bacteria, influences the ratio of stable nitrogen isotopes, i.e. ${}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}$ (δ¹⁵N). Usually, plants discriminate against ¹⁵N, resulting in lower δ¹⁵N values in plant tissue than in their N source (Dawson et al., 2002). But the main driver of variation in δ¹⁵N in plants is the type and availability of nitrogen sources and their isotopic composition (Evans, 2001). Therefore, stable isotopes are a useful tool to trace anthropogenic nitrogen input (Kendall et al., 2007; Hastings et al., 2013). There are three major sources of human-derived nitrogen with distinct δ¹⁵N signatures: (1) Synthetic fertiliser: widely used since the 1950s, produced via the Haber–Bosch process, has δ¹⁵N values close to air N₂ or slightly negative (δ¹⁵N between −4 ‰–0 ‰) (Bateman and Kelly, 2007). (2) Organic nitrogen from untreated wastewater and natural fertilizers (also in the form of atmospheric deposition of agricultural NH₃). Volatilisation and microbial transformations of this type of nitrogen come with strong ¹⁴N loss, giving the remaining nitrogen high δ¹⁵N values (approx. 5 ‰–20 ‰) (Bateman and Kelly, 2007). (3) Atmospheric deposition from fossil fuel combustion in power plants and engines. The δ¹⁵N of reactive nitrogen from fossil fuel combustion is generally negative but can be variable due to atmospheric transformations (Elliott et al., 2019) or soil microbial processing (Shi et al., 2014).

Nitrogen content (relative to other elements) in plant tissues is a species-specific trait, with species with competitive and ruderal strategies often being more N-rich and species which are stress-tolerant and have strongly lignified tissues being N-poor (Tilman and Wedin, 1991; Tian et al., 2019). However, it has been proven that the relative content of nitrogen in relation to other macronutrients, but also the nitrogen concentrations in the tissues, could be dependent on nitrogen availability in the environment (Palpurina et al., 2019). Moreover, under elevated CO₂, the nitrogen concentration in plant tissues tends to decrease (Cotrufo et al., 1998), which can be a result of multiple mechanisms, including the simple dilution effect but also changes in root architecture (Taub and Wang, 2008).

Stable carbon isotopes (δ¹³C) in plant biomass can be used to assess long-term changes in carbon cycling, photosynthetic conditions, and environmental stressors affecting plant life (Cernusak et al., 2013). The carbon isotopic signatures of terrestrial plants reflect CO₂ values, so studying historical samples, including plant biomass, provides evidence for the Suess effect: the progressive depletion of atmospheric ¹³C due to anthropogenic fossil fuel emissions. This is because the carbon returned to the atmosphere through the combustion of fossil fuels is depleted in both ¹³C and ¹⁴C compared to atmospheric carbon dioxide (Keeling, 1979). To disentangle the effect of fossil-fuel combustion derived variation of δ¹³C from other possible factors a correction for the Suess effect has to be applied (Dombrosky, 2020). On a much smaller scale, changes in CO₂ concentration also affect ¹³C discrimination by plants (Keeling et al., 2017). Photosynthesis in general is a process strongly discriminating against ¹³C, but the magnitude of this discrimination depends on plant physiology and may be altered by stomatal activity: long periods of closed stomata (higher water-use efficiency) mean that the plant is using a limited CO₂ pool, resulting in lesser discrimination and higher δ¹³C values (Seibt et al., 2008). Aquatic plants can use dissolved inorganic carbon (DIC) for photosynthesis; in fully submerged plants, it is the only carbon source. As opposed to atmospheric CO₂, δ¹³C of DIC (and aquatic plants) undergo various site-dependent changes in both directions, depending on the source of carbon, type of water body (Pronin, 2024), water chemistry – mainly pH (Thomson et al., 2001) and abundance of DIC (Schulte et al., 2011; Chappuis et al., 2017). This variety of factors determines the final value of δ¹³C, making the practical use of carbon isotopes as a marker of eutrophication in wetland ecosystems comparable to the use of nitrogen isotopes (Wang et al., 2015; Heredia et al., 2022).

A couple of studies have already used historical plant material stored in herbaria to track changes in the isotopic and chemical composition of plant biomass since the 19th century (Peñuelas and Filella, 2001; McLauchlan et al., 2010; Alldred et al., 2024; Kühn et al., 2025) or even the 18th century (Peñuelas and Matamala, 1990). Some of these studies have identified similar trends; for example, a decrease in δ¹⁵N was observed in foliar samples during the 20th century, as demonstrated for prairie plants in Kansas (McLauchlan et al., 2010) and Mediterranean terrestrial plants in north and east Spain (Peñuelas and Filella, 2001). However, McLauchlan et al. (2010) indicated that the main reason for this trend was an increase in atmospheric CO₂ levels (and thus stronger N-limitation) thought, Peñuelas and Filella (2001) attributed it to higher anthropogenic nitrogen input into ecosystems. Furthermore, McLauchlan et al. (2010) found that leaf nitrogen content decreased over the 20th century, while Kühn et al. (2025) observed an opposite trend in vascular plants from mesic grasslands and crop fields in Germany, a country with a similar, to our study area, intensive and crop-based agricultural sector (Eurostat, 2025). It is also worth noting that most of the aforementioned studies focused on terrestrial habitats, while recent study by Alldred et al. (2024) also found herbarium vouchers of macroalgae to be a very promising source of information in studying past biogeochemical trends in water environments.

To date, no studies have attempted to compare at the same time different habitat types that may be exposed to varying levels of anthropogenic nitrogen. Similarly, no studies have analysed the correlation between the magnitude of change and local land use factors. Considering the above discrepancies and gaps in knowledge, we formulated the following hypotheses to verify in our study:

1. Wetland plants and dry grassland plants exhibit different patterns of shifts in stable isotope composition.

2. Effects of eutrophication are more visible in wetlands, especially in submerged species, which show a drop in δ¹³C and an increase in δ¹⁵N.

3. TN in plant leaves increased over the last century.

4. The change in TN and δ¹⁵N is correlated with the proportion of agricultural areas' consumption of artificial nitrogen fertilisers in agriculture, and the atmospheric nitrogen deposition in the modern landscape.

5. δ¹³C of dry grassland plants decreased by about 1.5 ‰ over the last century due to the Suess effect.

2 Methods

2.1 Selection of species and specimens for analyses

For our study, we selected species that had to be represented by well-preserved and well-developed herbarium vouchers in the Herbarium of the Faculty of Biology, University of Warsaw (WA). Those herbarium vouchers should be collected before 1939, so before the artificial fertilisers became widespread in Central Europe (Merta-Staszczak and Zwyrzykowska-Wodzińska, 2017; Tian et al., 2022), from at least six locations that could be linked to contemporary occurrences of the species. We considered a species to be currently present at a given location (i.e. appropriate to the location of the herbarium specimen): (i) if we had personally observed the species there in the last decade, (ii) if we could identify patches of potential habitat of the species on ortophotographs (available at https://www.geoportal.gov.pl/, last access: 15 September 2024), or (iii) if we found observations of the species reported in recent botanical literature or in iNaturalist (https://www.inaturalist.org/, last access: 15 September 2024).

To reduce the impact of various potential nitrogen sources on the chemical composition of plant biomass, we excluded species that fix nitrogen via symbiotic bacteria. We chose plant species from taxa which are known to be non-mycorrhizal (Caryophyllaceae, Brassicaceae, Hydrocharitacae, Potamogetonacea, Nymphaea alba) or weakly associated with fungi (especially in wetland habitats) and if so, mostly of AM mycorrhizal type through which they obtain rather phosphorus than nitrogen (Alismataceae, genus Carex) (George et al., 1995; Miller et al., 1999; Cornwell et al., 2001; Unger et al., 2016). All of our species were C₃ plants, as patterns of excessive nitrogen acquisition and their impact on the isotopic composition of plant leaves seem to be more uniform among plants with the C₃ photosynthetic type (Yao et al., 2011). Our selection included relatively common species in two habitat types: dry grasslands and wetlands (springs, lakes and rivers) (see Table 1). The stability of moisture conditions in these habitats (often deficient in grasslands and always abundant in wetlands) helped to reduce the impact of variation in water availability at the time of material collection on the results. Finally, we selected 14 vascular plant taxa, two of which were aggregations of two closely related species (Table 1).

Table 1 Taxon characteristics and number of paired samples collected from historical and appropriate contemporary specimens. Species nomenclature follows World Flora Online (2025).

Figure 1 Localities of the sampled populations of dry grassland species (on the left) and wetland species (on the right). Overlapping locations were randomly dispersed to increase visibility (spatial data source: Bohn et al., 2001–2003).

2.2 Study area

The spatial arrangement of our study was predominantly driven by the origin of the herbarium specimens of the selected species. As Warsaw was part of the Russian Empire until the First World War, most of the materials stored in the herbarium from before that period were collected within its former borders, except in the vicinity of the Tatra Mountains, which was a popular destination among botanists at that time. We excluded vouchers collected outside the current Polish borders. To avoid oversampling certain popular areas, we did not include in our study samples of the same species collected less than a few kilometres apart. Finally, majority of the selected pairs originated from the eastern half of Poland (Fig. 1). According to the Köppen-Geiger climate classification, the region predominantly has a warm temperate, fully humid climate with a warm summer (Cfb) (Kottek et al., 2006). The mean annual temperature of the sampled locations varied from 6.35 to 9.05 °C (mean 8.26 °C), while annual precipitation ranged from 496 to 1187 mm (mean 632 mm). The mean temperature of the coldest month (January) ranged from −3.75 to −0.15 °C (mean −2.01 °C), while the mean temperature of the warmest month (July) ranged from 16.15 to 19.45 °C (mean 18.67 °C) in the period 1981–2010 (Karger et al., 2017).

2.3 Sampling of plant material

Each herbarium voucher sample consisted of two to five leaf fragments collected from one or several individuals, depending on the availability of material and the state of preservation of the specimens. According to information obtained from the head of the WA Herbarium, the specimens stored there were not prepared or preserved using nitrogen-containing substances. We attempted to obtain contemporary samples in the same locality as the historical one. The accuracy of relocation was related mostly to the level of detail provided in the herbarium voucher description. If the original location description was precise (it allowed for site relocation with uncertainty less than 1 km), we searched for the species in the same location, and if we were unable to find it there, then we searched in a radius of a maximum of 10 km from the original site. However, if the description was very general (e.g. pointing to a large town), we looked for the focal species in the broader vicinity of that town in suitable habitats. The total number of samples collected was 198, comprising 99 historical samples collected between 1855 and 1938 (median 1905) and 99 contemporary samples that were paired with their respective locations. All contemporary samples were collected between April (for the early-flowering species as Arabis hirsuta and Erysimum spp.) and September (for some of the aquatic plants) of 2024. Each contemporary sample, similarly to the historical ones, consisted of two to five leaf fragments collected from several individuals per locality (Fig. 1). We selected fully developed, green and healthy plants that were mostly in the flowering or early fruiting phase. Additionally, we collected and deposited two to three specimens from most of the currently sampled populations as herbarium vouchers in the Herbarium WA for future reference.

2.4 Laboratory analyses

Fresh leaf samples were dried in 40 °C for 24 h within 1 d after collecting. In our analysis we did not perform pre-treatment of samples with acid to remove carbonates as it can affect both total nitrogen and nitrogen isotopic composition, as well as δ¹³C values of organic carbon (Brodie et al., 2011; Schlacher and Connolly, 2014). Most of our samples were collected from the parts of the plants growing above the water and we wanted to follow a standardized protocol for all of the samples. In our dataset we had only two submerged species which could contain CaCO₃ on their surface: Stuckenia pectinata and Elodea canadensis. To investigate the possible presence of surface carbonates and their impact on our results we performed an additional check for those species, which details are described in Appendix A. The conclusions of this test were that there was no measurable effect of carbonates on the results for S. pectinata, while the possible effect on E. canadensis was minor and had no impact on the overall results or our reasoning (see Appendix A). Each sample (both fresh and historical) was homogenized by crushing in a pestle and mortar or cut into small pieces using scissors. Two subsamples were weighted into tin capsules – approx. 5 mg for quantitative analysis and 3 mg for stable isotope analysis. Total carbon and nitrogen were measured using Thermo Flash 2000 Elemental Analyzer with TCD detector. Stable carbon and nitrogen were measured in a Thermo Delta V Plus Isotope Ratio Mass Spectrometer coupled with a Flash 2000 Elemental Analyser. Stable isotope ratios are reported as δ (‰) in line with the equation:

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}_{\mathrm{sample}}={\displaystyle \frac{R_{\mathrm{sample}}}{\left(R_{\mathrm{standard}}\right)-\mathrm{1}}}\times \mathrm{1000},\end{array}$$

where R is the isotopic ratio (${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ or ${}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}$), and standard is an international reference value: atmospheric air for δ¹⁵N and Pee Dee Belemnite (PDB) for δ¹³C. International standards were measured for calibration and measurement precision, which were <0.1 ‰ for δ¹³C and < 0.2 ‰ for δ¹⁵N. We use four standards for the calibration curve (δ¹³C and δ¹⁵N values respectively in brackets): USGS65 (−23.29 ‰, 20.68 ‰), USGS40 (−4.52 ‰, −26.39 ‰), USGS41a (47.55 ‰, 36.55 ‰) and Urea#3a (5.89 ‰, 42.05 ‰) from Schimmelmann Research and IAEA-600 (−27.77 ‰, 1.02 ‰) for recalculation and precision assessment. A calibration curve is prepared at the beginning of each day of measurements, along with 3 IAEA-600 standards. Measurements only started if the results of all three samples do not differ to the certified values by 0.2 ‰ for δ¹³C and 0.1 ‰ for δ¹⁵N. Daily, about 40–60 samples were analysed. Each 10 samples a IAEA-600 standard was measured, and sample results were only included, if the standard result falls within the abovementioned threshold. In order to facilitate comparison of the samples, we corrected the values of the historical specimens for the Suess effect (Dombrosky, 2020). The obtained values of all measurements, including the uncorrected δ¹³C values, are available in Table S1 in the Supplement. We performed the analyses in the Laboratory of Biogeochemistry and Environmental Conservation at the University of Warsaw Biological and Chemical Research Centre (CNBCh UW).

2.5 Anthropogenic nitrogen input and land use characteristics

To assess the current levels of anthropogenic nitrogen input at the sample collection sites, we obtained data on the consumption of mineral nitrogen fertilisers per hectare of agricultural land within given municipalities in Poland from the 2020 General Agricultural Census (GUS, 2021). Moreover, we used European air quality data for the 2023 dataset (European Union, 2025) as a source of data on the NO_x atmospheric deposition in the sampling locations. The dataset provides NO_x concentrations for vegetation-related indicators (annual averages) at a 2 km² grid. It combines air quality monitoring data, chemical transport model results, and other supplementary data using a regression-interpolation-merging mapping methodology (European Union, 2025).

Using GIS tools (QGIS Development Team, 2024), we calculated the proportion of agricultural land within 10 km buffers around the contemporary specimens collection sites. We used CORINE Land Cover (CLC, 2018) as a source of land-use cover data, which has a resolution of 100 m. Within each buffer size, we calculated the sum of all subcategories of Level 1 “2. Agricultural area” occurring in our study area: 2.1.1 Non-irrigated arable land, 2.2.2 Fruit trees and berry plantations, 2.3.1 Pastures, 2.4.2 Complex cultivation patterns, and 2.4.3 Land principally occupied by agriculture, with significant areas of natural vegetation. We then divided the sum by the total area of the respective buffer.

2.6 Data analysis

We performed all the statistical analyses in R ver. 4.4.2 (R Core Team, 2024), and visualisation of the results using the “ggplot2” R-package (Wickham, 2016). We used a two-sample, two-sided t test with unequal variances to see whether there were differences in the mean δ¹⁵N between samples of dry grassland and wetland plants in both the historical and contemporary specimen groups. Prior to other statistical tests, we calculated the differences in δ¹⁵N, total nitrogen content (TN) and δ¹³C by subtracting the measurements obtained from historical samples from those obtained from contemporary samples belonging to the same pair. We then used 90 % confidence intervals (CIs), which equals to one sided 95 % confidence of lower or higher mean, to test whether the differences were statistically significant from zero within each of the species studied, as well as within groups of species representing the two studied habitat types (dry grassland species vs. wetland species). The calculation of the CIs was done using one sample t test implemented in the “ci” function from the R-package “gmodels” (Warnes et al., 2024). To test whether differences in δ¹⁵N and TN over time were influenced by the habitat type, proportion of farmland in the surrounding area, the quantity of synthetic nitrogen fertilisers used, and the amount of NO_x deposited, we employed multiple linear regression. Additionally, we performed a linear regression model of δ¹³C against δ¹⁵N, incorporating all the wetland plant samples studied.

Figure 2 δ¹⁵N (a) and TN (b) in historical and contemporary foliar samples from dry grassland (n=51) and wetland (n=48) plants. The asterisks represent the significance level based on a two-sample t test assuming unequal variances (** p<0.01, *** p<0.001). The boxes indicate the upper and lower quartiles, the black rhombus indicates the mean, and the line within a box indicates the median, whiskers extend from the hinges to the largest/smallest values within 1.5 times the Interquartile Range (IQR) and the points represent the samples with values outside of this range.

3 Results

3.1 Nitrogen isotopes and content across habitats and species

The δ¹⁵N value was significantly lower in the foliar samples of dry grassland species than in wetland plant species, in both groups of contemporary and historical samples (see Fig. 2a). In contemporary samples, the difference between the means of plants representing the two habitat types was 3.83 ‰ (dry grassland species: mean ${\mathit{\delta }}^{\mathrm{15}}\mathrm{N}=-\mathrm{1.91}$ ‰; wetland species: mean δ¹⁵N=5.73 ‰), whereas in historical samples, the difference was 3.22 ‰ (dry grassland species: mean δ¹⁵N=1.24 ‰; wetland species: mean δ¹⁵N=4.45 ‰) (Fig. 2a, for more details see Tables B1 and B4). The value of TN was also significantly lower in the foliar samples of dry grassland species than in wetland plant species, in both contemporary and historical groups of samples (Fig. 2b). In contemporary samples, the difference between the means of TN in plants representing the two habitats was 1.06 % (dry grassland species: mean TN = 2.72 %; wetland species: mean TN = 3.78 %), whereas in historical samples the difference was 0.55 % (dry grassland species: mean TN = 2.92 %; wetland species: mean TN = 3.47 %) (Fig. 2b, Tables B2 and B4).

Figure 3 Mean differences between the δ¹⁵N value and TN of contemporary samples and historical samples across the species and habitat types. Each observation is the difference between a contemporary and a historical sample paired with a location. The thick lines represent the 90 % confidence interval (CI) for the population mean: black for significant difference and grey for non-significant one.

In the dry grassland species, the mean δ¹⁵N value was significantly lower in the contemporary samples than in the historical ones, with a mean difference of −3.14 ‰ (Fig. 3a). The mean difference in δ¹⁵N values was lower for all species representing this habitat type, but was significant for all except one (Carex flacca) (Fig. 3a, Tables B1 and B4). In the wetland species, the mean δ¹⁵N value was significantly higher in the contemporary samples than in the historical ones, with a mean difference of 1.28 ‰. However, the difference was not significant in any of the single species representing this habitat (Fig. 3a, Table B4). We did not confirm a difference in the mean TN between contemporary and historical samples of all dry grassland species pooled together (Fig. 3b). While most species in this habitat type had a lower mean TN in contemporary samples than historical ones, this difference was only significant in Carex praecox. Two other species showed an opposite trend, but only one of these (Silene otites) had a significantly higher mean TN in contemporary samples than in historical ones (Fig. 3b, Tables B2 and B4). The mean TN was significantly higher in contemporary than in historical samples of wetland species, with a mean difference of 0.31 % (Fig. 3b, Table B4). Significantly higher mean TN was also confirmed in three species that represent this habitat type: Cardamine amara, Elodea canadensis, and Stuckenia pectinata (see Fig. 3b and Table B4). The relative mean difference of TN followed similar patterns to absolute values as shown in Fig. C1.

Table 2 The results of the multiple linear regression models where the habitat type, the share of farmland within a 10 km radius, consumption of synthetic N fertilisers (average use in kg ha⁻¹ per municipality), and NO_x deposition were used as predictors of the difference in both δ¹⁵N and TN between contemporary and historical foliar samples paired by location. Significant p values (p<0.05) are marked in bold.

3.2 Change of differences against nitrogen input and land cover

The difference in δ¹⁵N between contemporary and historical samples paired by location was significantly higher in wetland plants than in dry grasslands and positively related to the share of farmland in a 10 km radius around the location (Table 2, for visualisation see Fig. C2). NO_x deposition was also positively related with the difference in δ¹⁵N, although at a lower significance level (Table 2). Among the predictors of the difference in TN between contemporary and historical samples paired by location only the habitat type was significant, indicating higher increase in TN in case of wetland plants than in dry grassland plants (Table 2).

Figure 4 Mean differences between the δ¹³C of contemporary samples and historical samples across the species and habitat types. Each observation represents the difference between contemporary and historical sample values (δ¹³C of contemporary sample – δ¹³C of historical sample) paired with a location. The thick lines are 90 % CI for the population mean: black for significant difference and grey for non-significant one. Please note that prior to the statistical analyses the historical values were corrected for the Suess effect.

3.3 Carbon isotopes across habitats and species

When samples from each habitat type were pooled, the mean δ¹³C value, corrected for the Suess effect, was significantly lower in contemporary samples than in historical ones. The mean difference was −0.51 ‰ for dry grassland and −3.85 ‰ for wetland species (Fig. 4). A significantly lower mean difference in δ¹³C, corrected for the Suess effect, was confirmed for Carex praecox among the dry grassland species and for the following five wetland species: Alisma spp., Cardamine amara, Elodea canadensis, Sagittaria sagittifolia and Stuckenia pectinata. The largest differences in the corrected δ¹³C values were found between contemporary and historical samples of E. canadensis and S. pectinata, with means of −8.38 ‰ and −7.55 ‰, respectively (Fig. 4, Table B4). We also found a negative correlation between the δ¹³C and δ¹⁵N values in these two species, but this was only significant for S. pectinata (R²=0.41, p<0.01; see Fig. C3).

4 Discussion

4.1 Different patterns of shifts in stable isotope composition between dry grasslands and wetlands

The first hypothesis was confirmed for δ¹⁵N, as well as δ¹³C. In dry grassland species, we found a significant decrease in δ¹⁵N values by 3.14 ‰ on average in contemporary plant leaf leaves compared to those from the pre-Green Revolution period. This is consistent with the results obtained by other authors: McLauchlan et al. (2010) detected a long-term decrease in δ¹⁵N in leaves beginning in 1940 – in their study, δ¹⁵N values in leaves decreased on average by 4.2 ‰ (change of mean value from 3.4 ‰ to −0.8 ‰) between 1940 and 2008. Similarly, Peñuelas and Filella (2001) detected a systematic decrease in mean δ¹⁵N values between 1920–1930 and 1985–1995 by 3.32 ‰ (from 2.26 ‰ to −1.06 ‰). Peñuelas and Matamala (1990) also found a downward trend in δ¹⁵N in herbarium material. These declining trends in δ¹⁵N values over time have also been detected in other environmental records such as lake sediments (Holtgrieve et al., 2011), tree rings (Poulson et al., 1995; Bukata and Kyser, 2007) and ice cores (Freyer et al., 1996; Hastings et al., 2009). McLauchlan et al. (2010) related the decreasing trend in foliar δ¹⁵N of plants from the North American prairie with increasing concentration of CO₂ in the atmosphere, leading to stronger N-limitation of primary productivity. The decrease in nitrogen availability is usually associated with a stronger decrease in δ¹⁵N values in leaves (Caldararu et al., 2022). However, when we compared changes in the total concentration of nitrogen in leaf tissues, our study did not confirm the trend of increasing nitrogen limitation resulting from the carbon fertilisation of ecosystems observed by McLauchlan et al. (2010). The decrease in the δ¹⁵N values of dry grassland plants observed in our study can be explained by several mechanisms, which are not mutually exclusive. Firstly, the cessation of free-range grazing and area of pastures in Poland over the last few decades (Radkowska and Radkowski, 2023; Sucholas et al., 2025) has caused livestock and their ¹⁵N-enriched faeces (Bateman and Kelly, 2007) to disappear from dry grasslands. In fact, the vast majority of our samples were collected from dry grasslands that, in recent years, were only mown or were undergoing secondary succession after abandonment of grazing. Secondly, and probably the most important driver, is the atmospheric deposition that has become an important source of nitrogen depleted in ¹⁵N (Hastings et al., 2009). Such reasoning was also proposed by Peñuelas and Filella (2001) to explain the declining trend in foliar δ¹⁵N of plants from dry habitats in Spain. They suggested that in dry conditions, the nitrogen pool is more dependent on atmospheric deposition than in the case of more humid and cold conditions, which promote microbial activity and organic nitrogen mineralisation. This effect may be exuberated in the recent decades by increasing temperatures which may further decrease soil water availability for microbial activity in dry grasslands, which are generally water-stressed ecosystems (Samaniego et al., 2018). Our results may indicate that these vulnerable habitats may be receiving a particularly high input of anthropogenic nitrogen via atmospheric deposition, which can reach areas that are remote from towns and farms (Holtgrieve et al., 2011).

In contrast, the δ¹⁵N values of wetland species had increased slightly compared to the pre-Green Revolution period. While confirming the first hypothesis, this also supports our second hypothesis. Such results could be explained by the dominance of nitrogen pollution derived from organic matter (with higher δ¹⁵N values) in aquatic habitats. This is due to human sewage entering surface water and groundwater (Bateman and Kelly, 2007), which can sometimes result in extremely high δ¹⁵N levels (Alldred et al., 2024). Despite significant progress in wastewater treatment in Poland in recent decades, environmental monitoring of rivers and lakes shows high levels of organic waste contamination (Szołdrowska and Smol, 2025). Our δ¹³C results for wetland plants (discussed below) support this explanation, as they also suggest that increased water pollution by allochtonous organic matter is probably the main driving force behind the observed changes in wetlands. In comparison with the pre-Green Revolution era, there is a possibility that the level of organic (nitrogen and carbon) pollution may be rising due to factors other than the quantity of organic fertilisers. Such factors may include the method and distribution of manure fertilisation (which was less intensive and more scattered in the past than in modern agriculture), uncontrolled sewage inputs from intensively developing settlements and tourist infrastructure close to rivers and lakes, and the runoff of peat from ploughed organic soils or drained peatlands (a fairly common phenomenon in Poland) (Krüger et al., 2015; Heredia et al., 2022). We also cannot exclude a short term trend related to disruption of artificial fertiliser markets following Russian invasion of Ukraine in 2022, which would push farmers to use more natural fertilisers in the years directly preceding our sampling. Indeed the use of artificial fertilisers decreased by 11 % in 2022 compared to 2020, and then by a further 6 % in 2023 (GUS, 2026). However, there is no data to indicate whether this was due to lower fertilisation levels or the replacement of artificial fertilisers with natural ones. In the regions studied, the ratio of nitrogen provided by artificial and natural fertilisers varies from nearly equal proportions (with slightly more organic fertiliser) in the Podlaskie region, to less than a third provided as natural fertiliser in the Lubelskie region (GUS, 2026). This pattern is probably even more complicated at smaller spatial scales, with sub-regions specialising in animal production, where manure is available and widely used, and locations where crop production dominates and artificial fertilisers are used more frequently. Alternatively, the observed pattern in wetland species could be attributed to increased runoff of artificial nitrogen fertilisers into water systems, leading to elevated nitrogen concentrations and resulting in higher rates of organic nitrogen mineralisation in water, followed by increased rates of nitrification, leading to elevated δ¹⁵N values (Diebel and Zanden, 2009). This is particularly relevant given that the use of organic fertilisers in European agriculture has remained at similar levels to those observed before the Green Revolution (Tian et al., 2022). In dry habitats, the complex microbial processes leading to ¹⁵N enrichment, as previously described for wetlands, are hampered due to water deficits; thus, the inorganic nitrogen resulting from atmospheric deposition (¹⁵N-depleted) may serve more directly as a source of that macroelement for plants.

It is worth noting that there was variation in the δ¹⁵N difference among the wetland species studied. The increase in δ¹⁵N was particularly high in the submerged species Stuckenia pectinata, which would agree with our second hypothesis, but the temporal change for this species was not statistically significant. Generally, the variation in δ¹⁵N difference within this group was greater than in dry grassland species, making statistical inference more challenging (given the relatively small number of observations). Overall, changes to wetland habitats in Poland over the last century seem to be relatively significant not only due to ongoing pollution, but also due to widespread drainage of wetlands for agricultural purposes (Janssen et al., 2016), increasingly frequent droughts leading to increased concentrations of pollutants (Szołdrowska and Smol, 2025) and rise of average water temperature following the climate change and other anthropogenic impacts (Tomalski et al., 2025). Moreover, the destruction of the wetland buffer zones due to intensification of agriculture may impede denitrification process in the catchment leading to more organic-derived nitrogen reaching directly water bodies, where our submerged plant grew (Walton et al., 2020). In many cases, finding current specimens of wetland species in the locations indicated in the herbarium vouchers was challenging or even impossible due to the disappearance of wetlands, profound changes in their trophic state, or changes to the surrounding landscape. Consequently, the sources of nitrogen and their isotopic signature may evolve in wetlands more complexly over time than in dry grasslands, making it difficult to observe general trends in this group.

4.2 Changes of total nitrogen content in plants over the last century

We found that the change in the total nitrogen content of the foliar samples varied greatly depending on the species and habitat. The increase predicted by our third hypothesis was confirmed only for wetland species. In the case of dry grassland species we observed no general trend of TN changes over time. At the species level one species, Carex praecox, showed statistically significant decline as observed by McLauchlan et al. (2010) who attributed TN decrease observed in prairie plants in Kansas to CO₂ fertilisation effects and increase of water use-efficiency under elevated CO₂ (especially important in dry habitats where primary productivity could be limited by drought) which could lead to increase in primary productivity and finally to productivity limitation by nitrogen. On the other hand, in our study, Silene otites showed the significantly opposite pattern, with an increase in TN, despite usually occupying habitats with extreme water deficits. Dianthus carthusianorum also had a slightly increased TN (although not significant), which may suggest a family-specific evolutionary strategy, as both species belong to the Caryophyllaceae family. On the other hand we also did not observe a general increase in foliar nitrogen content as Kühn et al. (2025). This may mean that for most of the species additional input of anthropogenic nitrogen have compensated for the increased demand for this nutrient related to the rising carbon dioxide levels but has not exceeded it as happened in probably even more overfertilised habitats studied by Kühn et al. (2025).

In the case of wetland plants there was an overall visible increase of total nitrogen content in foliar samples when contemporary plants were compared with herbaria vouchers from the pre-Green Revolution period. This trend was driven by strong and statistically significant responses in Cardamine amara, Elodea canadensis, and Stuckenia pectinata where we observed relative increases in foliar total nitrogen concentrations of more than 20 % (see Fig. C1). Those species were not standing out in the case of δ¹⁵N difference, but their δ¹³C difference was significantly bigger than expected from the global ambient δ¹³C trends. Both could point to water bodies being polluted by sewage and organic fertilisers, bringing both ¹⁵N-enriched nitrogen and ¹³C-enriched organic carbon, which, after mineralisation, is turned into DIC, which is utilised by those, mostly submerged, plants.

4.3 Factors influencing the patterns of change in chemical composition of plants over time

We confirmed our fourth hypothesis only for the trends in δ¹⁵N. The decline of δ¹⁵N in foliar samples was larger in areas with a lower share of farmland in the surrounding landscape. There are a number of studies that found positive correlations between the anthropogenic nitrogen load and δ¹⁵N values in diverse elements of the environment. For example, Xue et al. (2009) showed a significant positive relationship between river NO₃-δ¹⁵N values and the proportion of agricultural and urban land use. Dailer et al. (2010) obtained low δ¹⁵N values in algae collected from areas with low anthropogenic influence. An increase in δ¹⁵N related to agricultural and urban areas was found in macrophytes from coastal (Kohzu et al., 2008) and river ecosystems (Pronin et al., 2025). Alternatively, our results can be attributed to the increased use of fertiliser in agriculture today compared to the past (Tian et al., 2022), regardless of whether they are organic or artificial. Intensively used agricultural land, as a result of the use of nitrogen fertilisers, manure and other agricultural practices, is characterised by increased nitrogen availability and a higher rate of nitrogen circulation, resulting in greater losses of this element through processes such as denitrification, ammonia volatilisation and nitrate leaching (Craine et al., 2009; Hall et al., 2015; Gerschlauer et al., 2019). As a result, the remaining nitrogen pool is gradually enriched with ¹⁵N, which may be reflected in high δ¹⁵N values in leaves (Gubsch et al., 2011).

4.4 Changes in δ¹³C over time deviate from the Suess effect

We did not fully confirm our fifth hypothesis. After applying the correction for the Suess effect, we were left with a small (0.51 ‰) but significant decrease in δ¹³C since historical times in dry grassland foliar samples. This is contrary to the expected reduced stomatal opening time (Seibt et al., 2008) caused by water shortage, and is hard to explain. Perhaps our results should be interpreted as signals of increased water-use efficiency in plants due to increased ambient CO₂ concentrations, leading to more efficient CO₂ uptake and thus lower water loss, which reduces drought stress and lowers δ¹³C values in plants (Medlyn et al., 2001), even though climate change is leading to more severe droughts in terms of absolute water availability, as mentioned above. The significant decrease in δ¹³C in wetland plants agrees with our second hypothesis. Wetland plants can show such a response to eutrophication due to two mechanisms: (1) increased nutrient availability causes higher productivity that causes plants to increase stomatal conductance (in not fully submerged plants), which leads to higher fractionation (Wang et al., 2015). The parameter strongly correlated with δ¹³C in most studies with such results (Chappuis et al., 2017; Wrosz et al., 2025) is total phosphorus in water. We do not have such values for our study, but we can safely assume that P has increased due to anthropogenic input. (2) An oversupply of organic matter (with “terrestrial” δ¹³C values) from sewage and organic fertiliser into wetlands causes intensive bacterial mineralisation. At some point, most DIC can originate from this process, thus has very low δ¹³C values (Hellings et al., 1999). Heredia et al. (2022) attributed the decrease in δ¹³C of submerged plants in an eutrophication gradient to allochtonous organic matter decomposition. Similar findings – surprisingly low δ¹³C values – came from the study on urban rivers in NW Poland (Pronin et al., 2025). The simultaneous response on δ¹³C and δ¹⁵N (Fig. C2) confirms that eutrophication is the main driver for these changes and an important factor affecting wetland ecosystems at a large scale.

5 Conclusions

The present study revealed that, whilst the isotopic composition of nitrogen and carbon in foliar samples in Central Europe has changed over the last century, this change is very context-dependent. The magnitude and direction of the observed change are impacted by both the habitat type and the species identity and/or ecology. This may be regarded as a caveat for future researchers, indicating that studies on the stable isotope composition of plants should consider both their taxonomy and ecology. Nevertheless, the anthropogenic global change has left an imprint on the isotopic composition of the flora from the habitats under study. As would be expected, the alterations appear to be more pronounced in wetland environments than in dryland habitats. It is evident that both habitat types are among the most threatened in Poland (Perzanowska and Korzeniak, 2022) and the entire European Union (Janssen et al., 2016). However, it appears that human activity has resulted in more significant alterations to biogeochemical processes in wetlands compared to dry grasslands. While the latter have also diminished in size, the functionality of the remnants appears to be less affected. Moreover, the source of this disruption may vary between the two habitat types. Specifically, wetlands are exposed to a multitude of anthropogenic nitrogen and carbon sources, while dry grasslands are predominantly affected by alterations in atmospheric composition.

Appendix A: Check of the possible impact of carbonate contamination of Stuckenia pectinata and Elodea canadensis samples as well as of the acid treatment for their δ¹³C and δ¹⁵N values as well as total nitrogen (TN) content

To exclude the possible effect of surface carbonate contamination of the two submerged species of our dataset Stuckenia pectinata and Elodea canadensis for their δ¹³C values we employed a modified approach suggested by Kennedy et al. (2005). Due to the lack of visible CO₂ release, we could not use this as an indicator of complete carbonate removal. This enabled us to assume a very low carbonate content, so we used a weaker acid (0.3 M) to minimise the loss of organic matter.

We weighed four aliquots from each homogenised sample into silver capsules. We then added 100 mL of 0.3 M HCl directly to two of the aliquots in the capsules. After 1 h, we dried the samples at 40 °C overnight. This process was repeated once more. We then closed all the capsules, and analysed one pair of acid-treated and control samples for total carbon and nitrogen, and the second pair for δ¹³C and δ¹⁵N.

The results of this test show no effect of acid treatment on δ¹⁵N and total nitrogen (TN) content, nor on δ¹³C for S. pectinata. The only statistically significant difference was observed in E. canadensis. However, the mean difference in δ¹³C between treated and untreated samples was 0.44 ‰ (90 % confidence interval: 0.141 to 0.739), whereas the mean difference observed in our study between historical and contemporary samples was −8.38 ‰ (90 % confidence interval: −15.22 to −1.55). Thus, the possible error (including its uncertainty and the uncertainty of the results) is no greater than the difference observed between contemporary and historical samples. Moreover, as we treat contemporary and historical samples of plants inhabiting similar habitats in the same way, systematic methodological error would not lead to the detection of false patterns.

Appendix B

Table B1 Summary statistics of δ¹⁵N (values given in ‰) of the foliar samples of 14 species selected for the study, from both contemporary specimens collected in 2024 and historical specimens in the form of herbarium vouchers collected before 1939.

Table B2 Summary statistics of total nitrogen content (TN, values given in %) in the foliar samples of 14 species selected for the study from both contemporary specimens collected in 2024 and historical specimens in the form of herbarium vouchers collected before 1939.

Table B3 Summary statistics of δ¹³C (values given in ‰) of the foliar samples of 14 species selected for the study, from both contemporary specimens collected in 2024 and historical specimens in the form of herbarium vouchers collected before 1939. Please note that the historical values were corrected for the Suess effect.

Table B4 Mean differences between the δ¹⁵N, TN, and δ¹³C of contemporary samples and historical samples across the species and habitat types. The lower and upper 90 % confidence intervals for means were calculated using a one-sample t test. Please note that the historical values were corrected for the Suess effect prior to the calculations.

Appendix C

Figure C1 Relative mean differences between the TN of contemporary samples and historical samples across the species and habitat types. Each observation is the difference between a contemporary and a historical sample paired with a location. The relative difference was calculated according to a formula: (TN contemporary sample − TN historical sample) $/$ TN historical sample ⋅ 100. The thick lines represent the 90 % confidence interval (CI) for the population mean: black for significant difference and grey for non-significant one.

Figure C2 The relationship between the proportion of farmland within a 10 km radius and the difference in δ¹⁵N between contemporary and historical foliar samples, which were paired by location. The black line represents the predicted value of the dependent variable for each value of the main independent variable, while accounting for the influence of the other independent variables. Partial R² and p value refer to share of farmland in a 10 km radius. The dashed lines represent the regression lines for the given habitat type with their equations given with respective colours.

Figure C3 The relationships between δ¹⁵N and δ¹³C in all the foliar samples studied (i.e. historical and contemporary samples combined) from seven wetland species. The blue line shows the fitted linear regression model, and the grey area shows the 95 % confidence interval.