No impact of nitrogen fertilization on carbon sequestration in a temperate Pinus densiflora forest

Setting

This study was conducted in approximately 40-year-old naturally regenerated P. densiflora stands in Wola National Experimental Forest in Gyeongnam province in South Korea (35°12′ N, 128°10′ E; Table 1). The productivity of this forest is low, with a dominant tree height of 10 m at 20 years of age. Over the last 10 years, the mean annual precipitation was 1490 mm, of which one third fell during summer (July–August), and the mean temperature was 13.1 °C. The vegetation growing season generally lasts for approximately 200 days, extending from early April to October. The soil texture is a silt loam originating from sandstone and shale (clay 13.0 ± 0.8%, silt 44.1 ± 1.3%, sand 42.9 ± 1.6%; n = 18). The given texture results in volumetric water contents at 13.4 ± 0.7% (m³ m⁻³) at permanent wilting point (1500 kPa) and 40.7 ± 1.2% at field capacity (10 kPa)⁵⁵. The understory is covered with Lespedeza spp., Quercus variabilis, Q. serrata, Smilax china, and Lindera glauca.

In 2010, we selected two adjacent P. densiflora stands approximately 100 m apart from each other (180 m and 195 m above sea level, on slopes of 15° and 33°, both stands face south). Following a completely randomized design, we established nine plots (10 × 10 m² with a 5 m untreated buffer) within each stand, of which three were randomly assigned to annual NPK fertilization, three to PK fertilization, and the rest to a control treatment without fertilization. The fertilizer, composed of urea, fused superphosphate and potassium chloride (N₃P₄K₁) or P₄K₁ was added manually by deposition on the forest floor for 3 years in April 2011, April 2012, and March 2013. Over these 3 years, the NPK plots received 33.9 g N, 45 g P, and 11.1 g K m⁻², while the PK plots received 45 g P and 11.1 g K m⁻².

Tree and stand measurements

The standing biomass of trees was estimated using a combination of site-specific allometric equations based on destructive harvesting⁵⁶ and repeated measurements of the dimensions of all trees in each plot (5–18 trees plot⁻¹). The stem diameter at 1.2 m (D) was measured for all trees in each plot for which D was ≥ 6 cm. Selecting a representative tree in size for each plot within the 4 × 4 m² center of the plot, we measured the tree height (H) and crown base for the representative trees. Measurements were performed in April and September 2011, September 2012–2014, and October 2021. We observed no effect of fertilization on the relationship between D and H or between D and crown base, so we assumed no effect on the allometric functions for foliage or branch biomass. A dendrometer band (Series 5 Manual Band, Forestry Suppliers Inc., Jackson, MS, USA) was installed on 18 representative trees (one per plot) to monitor radial growth monthly.

Three 0.25 m² circular litter traps were installed 60 cm above the forest floor in each plot in April 2011. Litter was collected at 3-month intervals between June 2011 and March 2015. The litter from each trap was transported to the laboratory and then oven-dried at 65 °C for 48 h. All dried samples were separated into needles, bark, cones, branches, and miscellaneous components, and weighed separately.

In September 2014, we estimated the biomass of understory vegetation, separately for woody plants and herbaceous plants. All woody plants < 1.2 m tall in each plot were harvested, while herbaceous plants were harvested within three 1 × 1 m² subplots within each plot. The samples were separated into stems, branches, and leaves, and then oven-dried at 65 °C and weighed.

The leaf area index (LAI) was estimated for each treatment as the product of the specific leaf area (cm² g⁻¹) and the sum of the standing leaf biomass and annual foliage litterfall⁵⁷. As most litterfall occurred after summer, adding the foliage litterfall to the total foliage mass is necessary to accurately estimate the LAI during the summer season. The estimated LAI and recorded annual photosynthetically active radiation (PAR) were used to estimate the canopy-intercepted PAR using the Beer-Lambert light extinction function [intercepted PAR = PAR·(1–e^−k·LAI)] with an extinction coefficient (k) of 0.52⁵⁸.

Estimation of stocks and fluxes of soil carbon and nitrogen

Annual fluxes of C and N were estimated from the measured monthly fluxes. Biomass stocks were estimated annually, while soil stocks were estimated only at the end of the investigation.

Soil C fluxes were characterized by determining the total soil CO₂ efflux and heterotrophic respiration, whereas soil N fluxes were characterized by determining the net mineralization and nitrification rates using an in situ incubation method. Total soil CO₂ efflux was separated into root- and heterotrophic respiration using a trench method⁵⁹. A root exclusion barrier (50 cm diameter, 30 cm depth, PVC) was installed in each plot, 4 weeks before the treatments were initiated. Two soil collars (10 cm diameter, 5 cm depth, PVC) were installed inside the rooting barrier to measure heterotrophic respiration and another two were installed outside the rooting barrier to measure heterotrophic + root respiration. Vegetation was removed within the rooting barriers but litterfall was kept. Soil CO₂ efflux was measured once a month using an infrared gas analyzer (Model EGM-4, PP-Systems, Hitchin, Hertfordshire, UK) equipped with a flow-through closed chamber (Model SRC-2). Measurements were performed between 10:00 and 13:00 over the study period. Soil temperature was also measured at a depth of 8 cm near the soil CO₂ efflux collar using a digital soil temperature probe (K-type, Summit SDT 200, Seoul, Korea).

On each day when the soil CO₂ efflux was measured, two soil cores were collected from each plot at a depth of 5 cm using a 100 cm³ core soil sampler. One sample was placed in a plastic bag and the other was returned to the soil and incubated to estimate the net N mineralization rate. We thus collected two soil samples from each plot every month, one fresh and one incubated. Both samples were transported to the laboratory and their fresh weight was measured. A 10 g portion of the fresh soil sample was oven-dried for 48 h at 105 °C to quantify the soil’s gravimetric water content and the rest was kept to determine the concentrations of nitrate and ammonium and the soil pH (measured using a 1:5 soil water suspension with an ion-selective glass electrode; Istec Model pH-220L, Seoul, Korea). The bulk density of each soil sample was calculated from its gravimetric water content and fresh weight.

Ammonium and nitrate were extracted from soil samples (5 g) using 50 ml of a 2 M KCl solution in a mechanical vacuum extractor (Model 24VE, SampleTek, Science Hill, KY, USA). The resulting solutions were then immediately placed in a cooler at 4 °C for storage. The ammonium and nitrate concentrations of the solutions were determined using an auto analyzer (AQ2 Discrete Analyzer, Southampton, UK). The mineral N concentration was measured throughout the period when fertilizer was applied, from April 2011 to April 2014. The net rates of ammonification and nitrification were estimated based on the differences in the ammonium and nitrate contents of the soil, respectively, between before and after the incubation. The sum of these two rates was taken as the net mineralization rate. If the net mineralization rate was negative, it was regarded as a net immobilization rate.

Additional soil samples were collected from four randomly selected points in each plot at depths of 0–15 cm both 4 weeks before starting the treatments and at the end of the investigation (2014). Samples were pooled within a plot and their concentrations of C, N, and available P were determined to capture the soil’s initial condition and the treatment response. C and N were analyzed with an elemental analyzer (vario Macro cube, Elementar Analysensysteme GmbH, Germany). The available P concentration was determined by extraction using NH₄F and HCl solutions⁶⁰ and analyzed using a UV spectrophotometer (Jenway 6505, Staffordshire, UK). C and N stocks were estimated by measuring the concentrations of both elements in a 100 cm³ core and dividing by the bulk density of the soil in the 0–5 cm layer.

Data analysis

The trees’ NPP was estimated as the product of the annual biomass increment, woody components of litterfall, and the C content of each biomass component⁵⁶. Foliage NPP was estimated as the sum of the annual leaf litterfall and the annual foliage biomass increment, which in turn was estimated using the relevant allometric equation. Total soil CO₂ efflux and heterotrophic respiration were modelled as exponential functions of the soil temperature and moisture content in each treatment:

Here, R_{S ijkl} is the lth observation of the total soil CO₂ efflux or heterotrophic respiration in replicated plot i (i = 1–3) of treatment j (F, j = control, PK, NPK fertilizations) within stand k (S, k = 1, 2). T_S is the soil temperature (°C) and T_S × F is the interaction term. The bs are coefficients to be estimated and ε represents the final residuals. After finding no effect of soil moisture content on the temperature sensitivity of R_S (b₁) or base R_S (b₀) (p = 0.3 for total CO₂ efflux and 0.5 for heterotrophic respiration; Table S3), only the soil temperature and plot-stand factors were included in the final model. The models included both terms separately for each year to enable propagation of spatiotemporal variation in the annual estimates within a treatment. Because we only measured the soil temperature manually once a month, hourly estimates of the soil temperature were obtained by deriving a relationship between the soil temperature and the corresponding air temperature (Fig. 2a; R² = 0.954, relationship unaffected by treatment, p = 0.99). Using the above model and the air temperature, we estimated the annual soil CO₂ efflux and heterotrophic respiration.

The effects of the fertilization treatments on above- and belowground variables were determined using a mixed effect ANOVA model to account for the completely randomized design used in the adjacent two stands. The stand term was treated as a random effect and the fertilization treatments were treated as fixed effects.

Here, Response_ijkl is a response variable in replicated plot i of treatment j within stand k in year l (l = 1–4).  µ is the grand population mean, F is the fixed effect of the treatment, Y is the fixed effect of the year as a dummy, and S is the random effect associated with stand. For the monthly soil N mineralization and nitrification models, we replaced year Y with month (l = 1–36).

All statistical analyses and parameter estimations for the devised models were performed using R (v. 4.1.3): the lm function was used to devise linear relationships and the lme function of the nlme package was used to determine the coefficients of the mixed effect ANOVA model (Eq. 2). Means were compared and separated using Tukey’s test with a significance threshold of p < 0.05, which was performed using the lsmeans and cld functions of the lsmeans package. The final residuals (ε_ijkl) were normally distributed, which was confirmed visually by plotting them against predictions.

Ethical approval

The present study conducted in the Wola National Experimental Forest on P. densiflora trees, including the collection of plant/soil materials, complied with local (Gyeongsangnam-Do) and national (South Korea) guidelines and legislation, with the permission from the administering institute, National Institute of Forest Science.

Source: Ecology - nature.com