Whether nitrogen (N) availability will limit plant growth and removal of atmospheric CO₂ this century is controversial. Studies have suggested that N could progressively limit plant growth, as trees and soils accumulate N in slowly cycling biomass pools in response to increases in carbon sequestration. However, a question remains over the longer-term (decadal to century) feedbacks between climate, CO₂ and plant N uptake. The symbiosis between plants and microbes can help plants with mycorrhizal N uptake or biological N2 fixation – the pathway through which N can be rapidly brought into ecosystems and thereby partially or completely alleviate N limitation on plant productivity. Here we present results for plant N isotope composition (δ¹⁵N) in a peat core that dates to 15000 cal. yr BP to ascertain ecosystem-level N cycling responses to rising atmospheric CO₂ concentrations in the past. We found that an increase in atmospheric CO₂ concentration happened with a decrease in δ¹⁵N values of both Sphagnum moss and Ericaceae over this time period when constrained for climatic factors. A modern experiment demonstrated that δ¹⁵N of Sphagnum mosses decreased with increasing N2 fixation rates. These findings suggested that N2 fixation in Sphagnum moss by symbiosis with cyanobacteria and N uptake in Ericaceae by symbiosis with mycorrhizal fungi both likely increased with rising atmospheric CO₂ concentrations, highlighting a longer-term feedback mechanism whereby N constraints on terrestrial carbon storage can be overcome. 

1. Study sites and sampling monoliths

The peat core (42.227° N, 126.525° E) was collected in Hani, a valley peatland situated in the middle of the Longgang Mountain Chain which is part of the western flank of the Changbai Mountains, Northeast China. The present surface vegetation community mainly comprises of Betula ovalifolia, Ledum palustre, Phragmites australis, Vaccinium uliginosum, Carex limosa, S. palustre, S. magellanicum, S. cuspidatum, S. subsecundum, S. imbricatum, Polytrichum strictum and Aulacomnium palustre. The mean annual temperature and precipitation range from 0.8 to 5.4 ℃ and from 490 to 1300 mm, respectively from local instrumental records (Jingyu County) during 1955–2019 AD. The length of this core was 8.77 m. The peat core was described and cut with stainless steel at 1 cm intervals in the field. All of the samples were packed in plastic bags and stored at 4 ℃ in a refrigerator for further analyses.

2. Radiocarbon dating

The chronology of the peat core was based on the ¹⁴C accelerator mass spectrometer (AMS) analyses of plant residual stems and leaves subjecting to an acid–alkali–acid treatment (Table S1). The graphite samples were finished in our ¹⁴C dating preparation lab and then sent to the NTUAMS laboratory at the National Taiwan University (NTU) for AMS measurement with a HVE 1.0 MV Tandetron Model 4110 BO-Accelerator Mass Spectrometer (AMS). The measured ¹⁴C ages were converted to calibrated calendar ages (a BP = years before 1950 AD) with the IntCal13 calibration curve (Reimer et al., 2013) using the CALIB 8.10 program and expressed as cal. yr BP with 2σ ranges. The calendar age was calibrated using the Bacon age-depth model (Figure S1).

3. Plant macrofossil preparation procedures

Sub-samples (1 cm3) for plant macrofossils were analyzed using standard method (Mauquoy et al., 2010). Plant residues were washed with deionized water after alkaline solution treatment and then examined under a stereomicroscope at 10–20 magnification and macrofossils were identified. Moss leaves and small seeds were examined at high magnifications (×100–400). Five plant groups: Ericaceae, Carex, Equisetum fluviatile (E. fluviatile), Sphagnum (S. palustre) and Drepanocladus (washed only with deionized water) were picked for further isotope analyses because of their high abundance (> 10 mg for each sample) along the whole core and their residues were stored in a 5 mL centrifuge tube in a 4 ℃ refrigerator. Tilia 2.0.41 (Grimm, 1992) was used to calculate plant macrofossil percentages and to plot the diagram.

4. Stable nitrogen isotope analyses

Each residue sample of the five selected groups was dried at 50 ℃ for 48 h, homogenized in a ball mill, and analyzed for N concentration and isotopic abundance in an elemental analyzer coupled online with an isotope ratio mass spectrometer (DELTA V Advantage; ThermoFisher, Germany).

N isotopic abundance was expressed as δ15N and was calculated as follows: δ (‰) = [(R_sample/R_standard)-1)], where R is the ¹⁵N/1⁴N ratio of the sample and the correspondent standard, air N₂. To check the system stability, standard materials with known δ¹⁵N values were analyzed between every ten samples. Standard deviations of laboratory standards for δ¹⁵N were less than 0.15 ‰.

5. Analysis of bulk peat and calculation of nitrogen accumulation rate (NAR)

Each peat sample was weighed before and after drying at 105 ℃ until constant weight was reached. Dry bulk density (DBD) was calculated from the dry weight and sample volume. The dried samples were milled and homogenized for C, N content and C/N ratio analysis. Nitrogen accumulation rate (NAR) was calculated by the following equation based on raw data including peat accumulation rates (PAR), DBD and N contents. NAR = PAR * DBD * N (%) * 100. Peat accumulation rates (PAR) (cm yr^-1) for each profile were calculated for peat layers between the dated depths of the peat profile (Tolonen and Turunen, 1996).

6. CO2 enrichment experiment of living Sphagnum and calculations of N₂ fixation rates

In early September 2020, the dominant Sphagnum species (S. fallax, S. magellanicum and S. fuscum) were collected from both hollows and hummocks in Dongfanghong peatland, Northeast China (42.185°N, 128.312°E). This experiment had two periods: the treatment period (one month) and the assay period (48 h). First, we incubated the Sphagnum in lab for a month under two treatments (ambient CO₂ and elevated CO₂ with 600 ppm) with seven replicates. Plants were incubated under prevailing light conditions (7200 lx, 400–700 nm, day length of 16 h) and in the dark (0 lx, night length of 8h) and at 25 ℃ for daytime and 20 ℃ for night in incubators. Distilled water was added into each sample every two days to prevent overdrying and maintain the original water level.

After treatment period, Sphagnum mosses were incubated for assay of N₂-fixation rates. Samples were incubated in 120 mL glass vials and 10 mL of water was added into the vials to keep the samples moist. A volume of 20 mL of air was removed, and 20 mL of 15N2 tracer gas [98% (vol/vol)] were injected into the vials to reach headspace enrichment levels of 16% (vol/vol) for ₁₅N₂. Incubations were terminated after 48 h by opening the vials, emptying water, and freezing them at -20 ℃. Then the samples were dried at 50 ℃ to a constant mass and analyzed for their bulk δ¹⁵N values and N contents using a stable isotope ratio mass spectrometer (DELTA V Advantage; ThermoFisher, Germany). Fe and Mo (below detection limit) contents in the tissue of the incubated Sphagnum samples were determined using sector field ICP-MS (7500C; ThermoFisher, Germany) employing well established analytical procedures (Kempter et al., 2017). Powdered samples were dissolved in a microwave autoclave (CEM, Mars6) using 3 mL high purity HNO₃ and 0.1 mL HBF₄. The certified plant material (BW0000-2016) was used for quality control.

• N₂ fixation rate was calculated according to Leppänen et al. (2013):
• N uptake (nmol g^-1 h^-1) = (1/100) * (%N/100) * [atom%_sample - atom%control/MW(N₂)] * [(100/atom%_headspace)/48] * 10⁹ (Eqn. 1); where % N is the N percent of the dried sample, and MW(N2) is molecular weight for N2 (28.0134462);
• atom%headspace = f_N2 * (atom%_N2) - atom%_15N2) + atom%_15N2(Eqn. 2); where f_N2 is the percentage for N2 in air or water, atom%_N2 = 0.36630, and atom%_15N2 = 98;
• f_N2 = V_N2/(V_N2 + V_15N2)  (Eqn. 3); where V is the volume of the corresponding gas.

7. Data of atmospheric CO₂ concentrations and other hydroclimate factors 

Continuous records of the atmospheric CO₂ concentrations were from published ice core records and recent measurements on firn air and atmospheric samples spanning the interval from the penultimate glacial maximum (～156 kyr BP) to the beginning of the year 2016 AD (Köhler et al., 2017). The newest published age scales were used for the ice core data and atmospheric CO₂ concentrations represented an integrated global signal in this paper. We used raw data instead of splined-smoothed calculations. Meteorological data of mean annual temperature (MAT) had two sources depending on time. From 2010–1950 AD, the hydroclimate factors were from local instrumental records (Jingyu County). Before 1950 AD, the hydroclimate data were from previous publications: MAT was based on air temperature derived from brGDGT distributions with calibration in the same peatland as our sampling site of the peat core (Zheng et al., 2017). DWT was reconstructed by a transfer function (WATOL-inv) using plant macrofossil of the peat core (Yang et al., 2022). All the climatic factors were matched to plant δ¹⁵N data within plus and minus five years. 

8. Statistical analyses

Normality assumption was tested with Shapiro-Wilk´s test and homogeneity of variances was assessed using Bartlett´s test. Treatment effects were tested using one-way analysis of variance (ANOVA) and a Tukey HSD post-hoc test. The level of significance (α) was set at 0.05. In addition, bivariate relationships between plant δ¹⁵N and climatic factors were estimated using zero-order correlations (same as Pearson’s correlations) and partial correlations by controlling for a single edaphic variable (i.e. DWT, MAT, TN, CO2). The detailed significance level of the partial correlation matrix was shown in Table S2. All statistical analyses were conducted using SPSS 19.0 (SPSS Inc., Chicago, IL, USA).