Relative importance of climatic and edaphic factors as drivers of plant δ¹⁵N along a longitudinal transect

ABSTRACT

Many studies have shown that climatic and edaphic factors influence the variations in nitrogen isotopes (δ¹⁵N) in terrestrial ecosystems. However, the relative importance of co-varying climatic and edaphic controls on plant δ¹⁵N remains somewhat unclear. To address this issue, regional-scale plant (L. chinensis) and soil samples were collected along a longitudinal transect across the temperate grasslands in northern China. The influences of climatic and edaphic factors on the plant δ¹⁵N were disentangled using the variable importance in projection (VIP) approach. We found that the climatic and edaphic variables explained 52.4% and 31.2%, respectively, of the variation in plant δ¹⁵N. The mean annual precipitation (MAP), soil N, and mean annual temperature (MAT) were the most important variables, but the soil water content (SWC), aridity index (AI), and soil clay content were also important variables. The soil pH and soil C: N ratios were identified as unimportant variables. Furthermore, a high-performance model for simulating plant δ¹⁵N values was constructed based on the important variables (VIP > 0.8). Further investigations should focus on the effects of the interaction between the climatic and edaphic variables on the N cycles in the temperate grasslands to provide more reliable predictions of plant δ¹⁵N.

KEYWORDS:

• Relative importance
• nitrogen isotopes
• climatic variables
• edaphic variables
• L. chinensis

Introduction

Temperate grasslands are one of the most widely distributed ecosystems worldwide, and they play important roles in carbon (C) and nitrogen (N) sequestration and climate regulation (Conant et al. 2017; Zhou et al. 2020). Unfortunately, most of the world’s temperate grasslands, especially China’s semi-arid temperate grasslands, are currently being degraded by global warming and excessive human activities (An and Li 2015; Emadodin et al. 2021). Given changes in land use and the global climate, a better understanding of the spatial patterns of the ecosystem functions in semi-arid grasslands is urgently needed to enable sustainable management of these grasslands (Allan et al. 2015). Diverse lines of evidence suggest that ecosystem nitrogen (N) stoichiometry can serve as a focal indicator of ecosystem functions and nutrient limitations (Elser et al. 2007; LeBauer and Treseder 2008; Zhou et al. 2014). Therefore, elucidating the spatial pattern of the N cycling and the natural abundances of N isotopes (δ¹⁵N) in semi-arid grasslands is key for understanding ecosystem functions. It is well known that the variations in δ¹⁵N in ecosystems are controlled by the rates and isotopic compositions of the N inputs and outputs, as well as the internal N transformations occurring in the vegetation-soil systems. In general, in addition to biotic factors, climatic and edaphic factors, including temperature, precipitation, pH, soil C and N contents, and soil texture, can directly or indirectly affect the δ¹⁵N values of plants by affecting the N transformation and release (e.g., ammonification, nitrification, and denitrification) in ecosystems (Robinson 2001; Pardo et al. 2007; Craine et al. 2009; Liu and Wang 2009; Wu, Wang, Chen 2018; Zhao et al. 2020; Luo, Viscarra-Rossel, and Qian 2021). As a result, the δ¹⁵N signal of plant-soil systems has become an effective tool for understanding ecosystem N cycling, and it has been used as a proxy for climate and environment changes (Hobbie, Jumpponen, and Trappe 2005; Ann-Kathrin et al. 2011; Ariz et al. 2015; Craine et al. 2015a). However, the influences of various climatic and edaphic factors on plant δ¹⁵N are not fully understood, which limits the use of plant δ¹⁵N in inferring ecosystem N cycling and environmental changes at the regional and global scales.

To understand the factors controlling plant δ¹⁵N, numerous studies have been conducted to investigate the patterns of plant δ¹⁵N along environmental gradients. In most cases, it has been shown that plant δ¹⁵N values are negatively correlated with the mean annual precipitation (MAP) and positively correlated with the mean annual temperature (MAT) at the regional and global scales (Martinelli et al. 1999; Aranibar et al. 2004; Bai and Houlton 2009; Ma et al. 2012; Peri et al. 2012). Similar temperature/precipitation effects have also been found in studies conducted along altitude gradients (Mannel, Auerswald, and Schnyder 2007; Craine et al. 2009; Nel, Craine, and Cramer 2018). However, some studies have reported that plant δ¹⁵N values do not correlate with MAP and MAT along environmental gradients (Wang et al. 2014; Zhao et al. 2020). This indicates that although the δ¹⁵N -climate correlations for plants seem to be well established at different spatial scales, the climate-driven changes in the δ¹⁵N values of plants remain somewhat uncertain.

In addition to the influences of climatic factors, nonclimatic factors, such as plant types (species), edaphic variables, and land use changes (e.g., agriculture and grazing), may significantly affect plant δ¹⁵N (Yang et al. 2013; Craine et al. 2015b; Lv, Xu, and Gao 2021). For instance, in a natural ecosystem, the δ¹⁵N values differ considerably between plant species, i.e., by more than 3.2‰, e.g., between N₂-fixing and non-N₂-fixing plants, mycorrhizal and non-mycorrhizal plants, different photosynthetic plants, and plants of different genotypes (Hobbie, Jumpponen, and Trappe 2005; Kahmen, Wanek, and Buchman 2008). This means that if the δ¹⁵N values of multiple plant species are selected as a whole, the relationships between the plant δ¹⁵N and climatic and edaphic factors may be obscured, and incorrect results may even be obtained. Besides that, several researchers have also suggested that plant δ¹⁵N may be affected by edaphic factors, such as the soil N concentration, pH, and other soil variables (Menge et al. 2011; Yang et al. 2013; Craine et al. 2015a) because these variables can play important roles in controlling the N-cycling processes through their influences on microbial activity in terrestrial ecosystems (Butterbach-Bahl and Gundersen 2011; Tang et al. 2018). For example, soil water can promote mineralization and microbial nitrification when the soil is unsaturated and aerobic (Butterbach-Bahl and Gundersen 2011). The pH of the soil can affect internal N-cycling processes via its effects on microbial nitrification and denitrification, and it can control N loss by affecting NH₃ volatilization (Peri et al. 2012; Mayor et al. 2014). The soil N content may regulate foliar δ¹⁵N values at different spatial scales (Craine et al. 2009). Furthermore, the soil texture has been demonstrated to influence N turnover by indirectly affecting other physicochemical properties of the soil (Kramer et al. 2003; Menge et al. 2011). Overall, previous studies conducted at different scales were very comprehensive and have greatly increased our understanding of the controls and patterns of plant δ¹⁵N in natural ecosystems. Nonetheless, our understanding of the patterns and controls on the δ¹⁵N values of plants is still limited to a certain extent due to the following two reasons. One limitation is our inadequate understanding of the extent to which individual soil variables regulate plants δ¹⁵N values because these soil variables jointly control the input and output of soil N and the overall isotopic fractionation through microbial decomposition (Craine et al. 2015b; Luo, Viscarra-Rossel, and Qian 2021). Additionally, previous studies have mainly focused on plant δ¹⁵N-climate relationships or the relationships between soil δ¹⁵N values and climatic variables (Wu, Wang, and Chen 2018; Zhao et al. 2020), and therefore, our understanding of the effects of edaphic factors on plant δ¹⁵N variations is less complete than that of the effects of climatic factors on plant δ¹⁵N variations along environmental gradients. Another limitation is that it is not very clear how co-varying climatic and edaphic factors independently influence the plant δ¹⁵N signatures of natural ecosystems. Luo et al. (2015) reported that the climate not only directly but also indirectly (via its effects on soil variables) significantly affects plant δ¹⁵N values. Although many studies have been conducted to explore the impacts of various factors on plant δ¹⁵N at different scales, the relative effect of each variable on plant δ¹⁵N has seldom been separated from the overall effects of environmental variables. Due to the potential co-varying effects of climatic and edaphic factors on N isotopes and the complexity of N cycling in ecosystems (Luo, Viscarra-Rossel, and Qian 2021), we know little about the relative contributions of the various variables to the variations in plant δ¹⁵N values along environmental gradients. For example, the increase in plant δ¹⁵N with increasing MAT may also be caused by decreasing precipitation, and yet the degree to which variables can explain the variations in plant δ¹⁵N values is still uncertain due to the covariation of temperature and aridity. This is likely one of the major causes of the relatively low reliability of interpretations of ecosystem N cycles and environmental change based on plant δ¹⁵N values reported in previous studies. Consequently, it is imperative that we gain a better understanding of the factors controlling plant δ¹⁵N. In addition, most previous studies primarily used relatively simple bivariate regression or ordinary multiple regression to determine how climatic and edaphic variables influence plant δ¹⁵N values at different spatial scales. However, such analyses cannot separate the influences of climatic and edaphic variables on plant δ¹⁵N values on the regional scale due to the mutual correlations among the explanatory variables.

The goal of this study was to disentangle the relative importance of co-varying climate and soil factors in controlling the variance of plant δ¹⁵N at the regional scale. Based on a literature review, we hypothesize that distinct factors control the plant δ¹⁵N signature. (1) Climatic factors such as the MAP, MAT, and aridity index (AI) can exert stronger effects on plant δ¹⁵N than edaphic factors because the climate significantly affects the plant δ¹⁵N directly and indirectly by affecting soil variables. For example, soil N mineralization and nitrification are largely controlled by temperature and precipitation. In addition, (2) the model for predicting plant δ¹⁵N based on important explanatory variables has a higher reliability than that based on all of the explanatory variables since the reliability of the model is largely related to the relative equilibrium between the degrees of the importance of the independent variables. In this study, a framework of how climatic and edaphic factors influence plant δ¹⁵N was developed (Figure 1). This framework is based on core theoretical assumptions. To test the above hypotheses, we collected plant (Leymus chinensis, a dominant perennial species with wide adaptability) and soil samples along a longitude transect across the temperate grasslands in northern China, and then, we determined the δ¹⁵N values and soil physicochemical properties of these samples. The variable importance in projection (VIP) method based on partial least squares (PLS) regression was used to analyze the contributions of the climatic and edaphic factors to the plant δ¹⁵N.

Figure 1. Conceptual framework of the climatic and edaphic factors influencing plant δ¹⁵N. The solid and dashed black arrows indicate positive and negative effects, respectively.

Materials and methods

Transect description

This study was conducted along an east-west transect across the temperate grasslands in Inner Mongolia, northern China (Figure 2). The longitude of the transect ranged from 112°48′ to 121°58′E, and the latitude varied between 42°19′ and 43°98′N. The transect was predominantly characterized by a semi-arid and sub-humid continental monsoon climate, with a MAT of 1.71–7.10°C and a MAP of 154–446 mm. Along the entire transect, from east to west, the main vegetation types were temperate meadow grassland, typical temperate grassland, and temperate desert grassland, in which L. chinensis communities were widely and continuously distributed (Yang et al. 2007). The soils related to these three types of grassland were predominantly chernozem, chestnut soil, and eolian sand soil, which had the same substrate age and belonged to the Kastanozem group, according to the Food and Agriculture Organization soil classification system (Cheng et al. 2009; Gerasimova 2010; Luo et al. 2013). The pH values of the 0–20 cm soil layer varied between 6.2 and 8.5. With such climatic and edaphic gradients, this transect provides an ideal environment for exploring the variations in the plant δ¹⁵N values and their influencing factors.

Figure 2. Distribution of the sampling sites across the temperate grasslands in northern China.

Sample collection

The 40 undisturbed sampling sites were selected 500–1000 m away from major roads along the transect across the temperate grasslands (Figure 2). The longitude, latitude, and altitude of each sampling site were determined using a Global Positioning System (Thales, Arlington, USA). The detailed information about the sampling sites is summarized in Table 1. During the months of July and August in 2017, samples of a widely distributed plant (L. chinensis) and soil samples were collected from 40 sites along this transect. At each site, five 1 m × 1 m quadrats within an area of 10 m × 10 m were randomly established for soil and plant sampling. Within each quadrat, the aboveground parts of healthy L. chinensis individuals were harvested, and then, all of the harvested aboveground parts from each site were combined into one sample. All of the plant samples were air-dried and brought back to the laboratory for further treatment. During the sampling of the plants, three soil cores (2.5 cm diameter) were collected from each quadrat to a depth of 20 cm. A total of 15 soil cores were obtained from each site and were mixed thoroughly to create one composite sample. After removing the fine roots and other coarse materials from the soil samples, each composite sample was sieved through a 2.0 mm screen and divided into two parts. One part was stored in a thermal insulation box (at 4°C) until its physicochemical properties could be analyzed, and the other part was oven-dried at 65°C in preparation for isotopic analysis. Additionally, the soil water content (SWC) of the 0–20 cm depth interval was determined in each quadrat using a Digital Soil Moisture Meter (Soil tester-300), and then, the mean value of the five quadrats in each site was calculated.

Table 1. The general information about the sampling sites along the studied transect.

Sample sites	Longitude (°E)	Latitude (°N)	Altitude (m)	MAT (°C)	MAP (mm)	Main soil type	Main community type	Average plant δ^15N (‰)
Tongliao1*	121.77	43.63	202	6.81	387	A	L. chinensis + Chloris virgata	−0.42
Tongliao2*	121.98	43.61	207	5.92	446	A	L. chinensis + Chloris virgata	−0.29
Tongliao3	121.96	43.60	205	5.90	442	A	L. chinensis + Chloris virgata	−0.41
Kailu1*	121.95	43.73	158	6.62	375	A	L. chinensis meadow	1.76
Kailu2*	121.99	43.57	184	7.10	396	A	L. chinensis meadow	−0.57
Kailu3	122.58	43.74	164	6.76	408	A	L. chinensis meadow	−0.36
Arhorqinqi1*	120.56	43.74	342	6.06	359	B	L. chinensis + C. virgata	−0.11
Arhorqinqi2*	120.75	43.66	281	6.49	358	B	L. chinensis + C. virgata	2.15
Arhorqinqi3	120.64	43.63	285	6.53	358	B	L. chinensis + C. virgata	0.63
Bairinzuoqi1*	119.36	43.88	486	5.43	339	B	S. baicalensis + L. chinensis	−0.24
Bairinzuoqi2*	119.40	43.89	472	5.26	353	B	S. baicalensis + L. chinensis	−0.34
Bairinzuoqi3	119.23	43.58	477	5.31	384	B	S. baicalensis + L. chinensis	0.87
Bairinyouqi1*	119.99	43.40	698	4.14	362	B	S. baicalensis + L. chinensis	0.54
Bairinyouqi2*	119.62	43.45	520	4.98	334	B	S. baicalensis + L. chinensis	1.50
Bairinyouqi3	119.71	43.44	453	5.34	329	B	S. baicalensis + L. chinensis	2.33
Linxi1*	118.19	43.61	758	4.19	341	B	L. chinensis + herbage	1.43
Linxi2*	118.54	43.53	648	4.86	339	B	L. chinensis + herbage	1.28
Linxi3	118.71	43.58	664	4.68	347	B	L. chinensis + herbage	1.20
Kesketonqi1*	117.71	43.47	1346	2.62	365	B	S. grandis + L. chinensis	0.98
Kesketonqi2*	117.84	43.39	1159	2.65	351	B	S. grandis + L. chinensis	0.42
Kesketonqi3	117.86	43.53	1050	3.42	347	B	S. grandis + L. chinensis	1.02
Xilin east1*	117.09	43.31	1267	1.75	398	B	S. grandis + L. chinensis	2.10
Xilin east2*	117.32	43.30	1269	1.71	408	B	S. grandis + L. chinensis	1.29
Xi Ujimqin1*	117.82	43.87	1168	2.23	230	C	S. baicalensis + L. chinensis	2.07
Xi Ujimqin2*	117.09	43.85	1047	2.87	219	C	S. baicalensis + L. chinensis	2.36
Xi Ujimqin3	117.64	43.95	1076	2.48	244	C	S. baicalensis + L. chinensis	2.11
Xilinhot1*	116.40	43.77	1099	2.28	324	B	S. grandis + L. chinensis	1.36
Xilinhot2*	116.54	43.71	1200	1.74	350	B	S. grandis + L. chinensis	2.67
Xilinhot3	116.60	43.76	1150	1.97	342	B	S. grandis + L. chinensis	1.79
Abagqi1*	115.01	43.98	1120	2.13	265	B	S. krilovii + L. chinensis	1.92
Abagqi2*	115.65	43.93	1133	2.02	295	B	S. krilovii + L. chinensis	1.42
Abagqi3	115.28	43.93	1164	1.91	286	B	S. krilovii + L. chinensis	2.06
Baiqi1*	115.13	42.24	1406	2.06	363	B	S. krilovii + L. chinensis	1.43
Baiqi2*	115.12	42.12	1405	2.03	364	B	S. krilovii + L. chinensis	1.21
Sonidzuoqi1*	114.09	43.85	1047	2.37	219	C	S. baicalensis + L. chinensis	2.87
Sonidzuoqi2*	114.47	43.92	1022	2.57	228	C	S. baicalensis + L. chinensis	2.65
Sonidzuoqi3	114.82	43.98	1131	2.09	259	C	S. baicalensis + L. chinensis	2.71
Erlanhot1*	112.53	43.67	998	3.64	154	C	S. baicalensis + L. chinensis	2.72
Erlanhot2*	112.73	43.70	987	3.63	159	C	S. baicalensis + L. chinensis	2.39
Erlanhot3	112.80	43.73	960	3.74	156	C	S. baicalensis + L. chinensis	2.42

Note: MAT and MAP are the mean annual temperature and mean annual precipitation, respectively. A denotes chernozem; B denotes chestnut soil; and C denotes eolian sand soil.

Laboratory measurements

The dried plant and soil samples were ground into a fine powder using a ball mill (NM200, Retsch, Haan, Germany), and then, the samples were stored in plastic bags. Following the method described by Harris, Horwath and Kesswl (2001), the soil samples were washed using 150 ml of 0.5 M HCL to remove any carbonate before the soil organic carbon (SOC) content was analyzed. Another balled soil subsample without HCL treatment was used to analyze the total N (TN) concentration and N isotope ratios. The SOC and TN concentrations of the soil samples were determined using the Walkley-Black modified acid-dichromate ferrous sulfate titration method and the micro-Kjeldahl digestion method coupled with colorimetric determination, respectively. The soil C/N ratio was calculated as the quotient of the SOC and TN content. Approximately 3 mg of plant samples and 65 mg of soil samples were loaded into a capsule and analyzed to determine their δ¹⁵N values using a Finnigan Delta^PlusXP (Thermo Scientific, Waltham, MA, USA) coupled to an automatic elemental analyzer (Flash EA1112, Thermo Finnigan, Milan, Italy). The natural abundances of ¹⁵N and ¹⁴N (δ¹⁵N) are reported in per mil (‰) relative to atmospheric N₂:

(1) ${\delta }^{15}N\left(0/00\right)=\left(R_{\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}}/R_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{d}}-1\right)\times 1000,$(1)

where R_sample is the (¹⁵N/¹⁴N) ratio of the sample, and R_standard is the ¹⁵N/¹⁴N of atmospheric N₂, which was used as a standard. The standard deviation of the repeated measurements was ±0.15‰ for the isotopic analysis.

In addition, the pH value of a subsample soil (10 g) with a dry soil-water ratio of 1:2 was measured using a pH meter (HI-9125, Hanna Instruments Inc., Woonsocket, RI). The clay component (< 2 μm) of another subsample was separated out using the ultrasonic energy method. The results of the particle size analysis are reported as the percentage by weight of the oven-dried soil.

Climate data for sample sites

The data for three climatic variables, including the MAT, MAP, and AI (the ratio of the potential evapotranspiration to the precipitation), for each sampling site, were obtained from the China Meteorological Data Service Center, Resource and Environment Science and Data Center, and from local weather stations. These climate data were the averages of observation data collected during a 30 year period. The data extraction was conducted using the spatial analysis tool in ArcGIS 10.2.

Statistical analyses

In this study, the VIP method was used to measure the degrees of importance of the climatic and edaphic controls on the plant δ¹⁵N and to identify the important controlling factors. The VIP method is a variable importance measurement technique based on PLS analysis. It combines multiple linear regression, principal component analysis, and correlation analysis, and it can eliminate potential collinearity among variables (Kano and Fujiwara 2013). It can be used to identify and screen independent variables with the most explanatory power for explaining the dependent variable based on their VIP scores and the absolute regression coefficient values of the PLS model. It also can reflect the relative importance of the independent variable to the dependent variable and can express the contributions of each independent variable to the dependent variable. The VIP value of the jth independent variable can be calculated as follows:

(2) $VI{P}_J=\sqrt{\frac{k}{{\sum }_{i=1}^h{r}^2\left(y;t_1,\cdots ,t_h\right)}}{\sum }_{i=1}^h{r}^2\left(y;t_i\right)w_{ij}^2$(2)

where k is the number of independent variables; and h is the total number of components; r is the correlation coefficient of the dependent variable (y) and the principal component ($t_i$) extracted from the relevant independent variables, which indicates the ability of $t_i$ to explain y; and $w_{ij}$ is the weight of the jth independent variable in the ith principal component, which is usually determined according to the default settings of the software. As the ability of x_j to explain y is conveyed through the principal component t_i, the stronger the ability of t_i to explain y, the greater the interpretability of x_j to y, and the larger the VIP value. Thus, the VIP values were monitored as threshold indicators to determine the relative contributions of the individual independent variables. The results of previous studies suggest that the VIP threshold defining the most influential variables is between 0.83 and 1.21 (Mkhabela, Bullock, and Sapirstein 2018). In this study, if VIP $\ge$ 1.2, the variable was classified as one of the most important variables. If 0.8 $\le$ VIP < 1.2, the variable was considered to be an important variable. These two types of variables should be retained in the PLS model. If VIP < 0.8, the variable was classified as unimportant and was excluded from the model. In addition, the relative effect was calculated based on the contribution (i.e., the VIP value) of each variable to the plant δ¹⁵N and its corresponding weight (i.e., the absolute value of the regression coefficient of each variable estimated using the PLS method). Then, the relative effect of each variable was scaled so that the sum of the influences of all of the variables was equal to 100. Furthermore, we also calculated the 95% confidence interval as the 2.5% and 97.5% quantiles of the relative influences via bootstrapping estimates, which serves as a measure of the uncertainty of the degrees of importance of the variables.

Additionally, to explore the relationships between the plant δ¹⁵N and the important climatic and soil variables, two PLS models were employed to model the plant δ¹⁵N as a function of the climatic and edaphic variables based on the same set of testing data. In the first model (PLS1), all of the climate and soil variables (excluding soil δ¹⁵N) mentioned earlier were used as the input, and the plant δ¹⁵N was the output. In the second model (PLS2), only the important variables with VIP values of > 0.8 were used as predictors of the plant δ¹⁵N. In addition, five accuracy indexes, namely, the cumulative interpretation rate in the X direction (R²X), the cumulative interpretation rate in the Y direction (R²Y), the cumulative prediction rate of the model (Q²), the mean square error (MSE), and the mean magnitude of the relative error (MMRE), were used as the prediction performance parameters. To accomplish this, the data samples were divided into two groups. The screening group contained 28 samples (marked with an asterisk in Table 1) and was used to determine the importance of the variable and to screen the important variables. The testing group contained 12 samples (without an asterisk in Table 1) and was mainly used to test the availability of the established prediction model. All of the statistical analyses were performed using the SIMCA-P 14.1 software (Sartorius-Umetric, Umeå, Sweden).

Results

Changes in plant δ¹⁵N and climatic and edaphic factors along the longitudinal transect

From west to east along the longitudinal transect, the plant δ¹⁵N, soil N concentration, soil clay content, and soil δ¹⁵N decreased (Figures 3(a,d)); while the SWC, AI, MAP, and MAT increased (Figures 3(e,h)). However, the soil C: N ratio and pH value did not exhibit an obvious trend along the transect (Figures 3(i,j)). The linear regression and bivariate correlation analysis revealed that the plant δ¹⁵N was significantly negatively correlated with the MAP, MAT, AI, and SWC (Figures 4(a,d)), and it was significantly positively correlated with the soil N content, soil clay content, and soil δ¹⁵N (Figures 4(e,g)). In comparison, the plant δ¹⁵N was not correlated with the soil C: N ratio and soil pH (Figures 4(h,i)). In addition, there was a negative linear correlation between the soil pH and MAP along the longitudinal transect, but this correlation was not significant (Figure 4(j)).

Figure 3. Variations in the plant δ¹⁵N and climatic and edaphic factors along the longitudinal transect. *, **, and *** denote significance levels of P < 0.05, 0.01, and 0.001, respectively.

Figure 4. Correlations between the plant δ¹⁵N and the climatic and edaphic factors along the longitudinal transect. *, **, and *** denote significance levels of P < 0.05, 0.01, and 0.001, respectively.

Importance of climatic and edaphic variables

The PLS regression revealed that the three-latent variables explained 83.6% of the variations in the plant δ¹⁵N values, and 86.5% of the variability in the predictor variables. The VIP values (Table 2) calculated based on the definition of the importance of the variable projection (Eq. (2(2) $VI{P}_J=\sqrt{\frac{k}{{\sum }_{i=1}^h{r}^2\left(y;t_1,\cdots ,t_h\right)}}{\sum }_{i=1}^h{r}^2\left(y;t_i\right)w_{ij}^2$(2) )) indicate that of the climatic factors, the MAP made the greatest contribution to the variability of the plant δ¹⁵N (VIP = 1.34), closely followed by the MAT (VIP = 1.23). Of the edaphic factors, the soil N content contributed the most significantly to the plant δ¹⁵N variation (VIP = 1.27), followed by the SWC (VIP = 1.03). According to the aforementioned criteria for determining the importance of a variable in explaining the variations in the plant δ¹⁵N in the temperate grasslands in northern China, the MAP, soil N concentration, and MAT were classified as the most important variables (VIP $\ge$ 1.2); while the SWC, AI, and soil clay content were classified as important variables (0.8 $\le$ VIP < 1.2). The remaining two variables (soil pH and soil C: N ratio) were classified as unimportant variables (VIP < 0.8). For the same data, using one climate principal component that was selected according to the minimum requirement criterion (Worley and Powers 2016) as the predictor, the PLS model could explain 52.4% of the variations in the plant δ¹⁵N. When the first two soil principal components were included, an additional 31.2% of the variance of the plant δ¹⁵N could be explained (Figure 5), of which the first principal component of the soil variables alone explained 27.6% of the variance of the plant δ¹⁵N. It should be noted that the fitted model in this study underestimated the high plant δ¹⁵N values and overestimated the low plant δ¹⁵N values (Figure 5). This bias was probably due to the interferences of the unimportant independent variables in the model.

Figure 5. The performance of the PLS model considering the soil variables (green circles) and not considering the soil variables (red circles) when explaining the variations in the plant δ¹⁵N values. The black dashed line is the 1:1 line, and the green and red lines are the trend lines with and without the soil variables as predictors, respectively.

Table 2. Variable importance in projection (VIP) values and regression coefficients (RCs) for plant δ¹⁵N based on the PLS-VIP method.

Variable	VIP	Standard error	RC	Standard error	Variable	VIP	Standard error	RC	Standard error
MAP	1.34	0.345	−0.318***	0.071	SWC	1.02	0.341	−0.153*	0.019
MAT	1.23	0.361	−0.247**	0.081	Soil pH	0.36	0.137	0.091	0.027
AI	0.88	0.285	0.054	0.020	Soil clay	0.87	0.107	0.056	0.021
Soil N	1.27	0.264	0.307**	0.081	Soil C:N	0.57	0.136	−0.026	0.043

Note: MAT, MAP, AI, and SWC are the mean annual temperature, mean annual precipitation, aridity index, and soil water content, respectively. *, **, and *** denote statistically significant correlations at the 0.1, 0.05, and 0.01 levels, respectively.

As is shown in Figure 6, without considering the other factors, the overall relative effect of the climatic variables on the plant δ¹⁵N was significantly larger than that of the edaphic variables (54.72% vs. 45.29%). For a single variable, we found that the relative effect of the MAP was the largest among all of the variables, accounting for 29.99% of the overall influence. The soil N content was the second most important, accounting for 27.47% of the total influence of all of the variables, followed by the MAT (21.38%), SWC (11.01%), and soil clay content (3.44%). The relative influences of the remaining three variables (i.e., AI, soil pH, and soil C: N ratio) were marginal compared to those of the five variables mentioned above, and together these three variables accounted for less than 7.0% of the overall influence of all of the climatic and edaphic variables.

Figure 6. The relative influences of the climatic and edaphic variables on the plant δ¹⁵N.

Verification of the relationships between the plant δ¹⁵N and the important climate and soil variables

Figure 7 and Table 3 show that model PLS2 had a better prediction than model PLS1. The predicted values of model PLS2 were closer to the observed values than those of model PLS1, with cumulative explained variances of R²X = 0.832, R²Y = 0.954, and Q² = 0.872. In addition, we also found that the MSE and MMRE of model PLS2 were lower than those of model PLS1 (Table 3). These results demonstrate that model PLS2, which was constructed using the important independent variables, is more accurate than the conventional model PLS1 in terms of their prediction performances.

Figure 7. Actual and predicted values of the testing samples based on different models (unit: ‰).

Table 3. Comparison of prediction accuracies of models PLS1 and PLS2 for testing samples.

Predicted model	Number of explanatory variables	R2X (cum)	R2Y (cum)	Q2 (cum)	MSE	MMRE (%)
PLS1	8	0.764	0.862	0.825	6.78	13.92
PLS2	6	0.832	0.954	0.872	3.58	7.15

Discussion

Importance of climatic controls on plant δ¹⁵N

Many studies have explored the influences of various climatic factors on the δ¹⁵N values of terrestrial plants on different spatial scales (Amundson et al. 2003; Cheng et al. 2009; Craine et al. 2015a). However, few scholars have paid particular attention to the relative importance of the effect of a climatic factor on the plant δ¹⁵N along an environmental gradient. In this study, the relative influences of the climatic variables accounted for approximately 54.72% of the total influence of all of the relevant variables on the plant δ¹⁵N, which demonstrates that the climatic factors exerted a significant control on the plant δ¹⁵N. This result supports the results of previous regional-scale studies, i.e., climate variables can explain much of the variations in plant δ¹⁵N (Yang et al. 2013; Craine et al. 2015b). In particular, the MAP was the most important climatic determinant of the plant δ¹⁵N in the temperate grasslands in northern China because it had the largest relative impact (Figure 6). The strong negative effect of precipitation on the plant δ¹⁵N values (Table 2) along the longitudinal transect is consistent with reports on global and other regional scales (Schulze et al. 1999; Brenner et al. 2001; Swap et al. 2004; Cheng et al. 2009; Craine et al. 2009; Liu et al. 2010; Ma et al. 2012). The above relationship can be explained by the indirect effect of precipitation on the soil δ¹⁵N and plant root development. First, precipitation can enrich soil N with ¹⁵N in temperate grassland ecosystems by influencing NH₃ volatilization. Generally, precipitation will affect the leaching of alkaline cations in the soil, so it controls the pH value of the topsoil. According to our analysis, there was a negative correlation between the soil pH and the MAP along the transect (Figure 4(j)). Due to this, soils in arid and semi-arid ecosystems usually have high pH values. It is well known that a higher pH value in the topsoil layer can significantly accelerate the volatilization of NH₃, resulting in high soil δ¹⁵N values (Handley et al. 1999; Murphy and Bowman 2009; Craine et al. 2015b). Because of the dependence of plant δ¹⁵N on soil δ¹⁵N (Figure 4(g)), the δ¹⁵N values of the plants decreased as the precipitation increased along the entire transect. Second, plant species in arid and semi-arid environments have deeper roots that can obtain deeper soil N resources, which are usually enriched in ¹⁵N.

In addition to the MAP, we found that the MAT was also an important climate variable controlling the δ¹⁵N values of the plant. However, it should be noted that in this study, the MAT had a negative impact on the plant δ¹⁵N (Table 2) along the transect, which is contrary to previous findings that the plant δ¹⁵N is positively correlated with the MAT at the global scale (Martinelli et al. 1999; Amundson et al. 2003; Craine et al. 2009; Kang et al. 2011). Moreover, several studies have also shown that the plant δ¹⁵N is not linearly correlated with temperature (Wang et al. 2014; Craine et al. 2015b; Chen 2018). These conflicting results could be due to three reasons. First, as was previously above, in the temperate grassland ecosystems in northern China, precipitation is not only the most important limiting factor that determines plant production, N absorption, and assimilation, but it is also the most important factor controlling the variations in the plants δ¹⁵N values (Chen and Wang 2000; Wang et al. 2014). Although increasing the temperature can accelerate ¹⁴N loss via nitrobacteria and result in positive soil and plant δ¹⁵N values, due to the synchronization of rain and heat in the eastern part of the transect (Figures 3(g)–(h)), the precipitation can contribute more negative δ¹⁵N values to the plants in this region, which may conceal the positive correlation between the plant δ¹⁵N and temperature to a certain extent (Liu and Wang 2009). Moreover, the decrease in the MAT along the western part of the transect (Figure 3(h)) further intensifies the negative correlation between the plant δ¹⁵N and temperature along the entire transect. Consequently, the coupled effect of these two factors may lead to a negative correlation between the plant δ¹⁵N and MAT. Second, the negative correlation between the plant δ¹⁵N and temperature is likely related to the differences in the N deposition in the eastern and western parts of the transect. Yang et al. (2007) reported that the eastern part of the transect has had a higher amount of N deposition than the western part in the past few decades. For example, the annual N deposition rate in Tongliao in the eastern part of the transect can be as high as 11 kg N ha⁻¹ year⁻¹, while the annual N deposition rate in Xilinhot in the western part of the transect is only 5 kg N ha⁻¹ year⁻¹. This intensified N deposition is likely to shift the N cycle into an unsteady state and to significantly change the pattern of the ¹⁵N abundance in the grassland ecosystem since the δ¹⁵N values of deposited N can vary from −10‰ to 5‰ (Handley et al. 1999). It has been reported that an increase in temperature can inhibit the absorption of N by L. chinensis in temperate grasslands and can significantly reduce the foliar N concentration, thus leading to a negative correlation between the plant δ¹⁵N and N deposition (Li 2015). Since the N deposition and temperature exhibited similar decreasing trends along the east-west transect in this study, the plant δ¹⁵N along this N deposition gradient may decrease with increasing MAT. An alternative explanation for the decrease in the plant δ¹⁵N with increasing temperature is that response of some parameters (e.g., nitrification and N loss) to elevated N deposition may be temporarily suppressed when the sampling sites along the N deposition gradient are not N saturated (Pardo et al. 2007). Although the exact reasons for the obvious negative correlation between the plant δ¹⁵N and MAT along the transect, which is not consistent with the results of other studies, are yet to be explored, the difference in the N deposition is likely one of the causes. Third, the relatively narrow temperature range (MAT: 1.71–7.10°C) along this transect, compared to those in previous studies, may cause this inconsistency. Shifts in the N cycling from organic or ammonium-dominated status to nitrate-dominated status have been hypothesized as a potential control of the plant and soil δ¹⁵N values (Amundson et al. 2003; Sun et al. 2010). Accordingly, it is possible that the relatively narrow range of the MAT along the studied transect is insufficient to cause a large shift in the N dynamics. This possibility is supported by the relationship between the variations in the plant δ¹⁵N and the MAP being stronger than that between the plant δ¹⁵N and the MAT (Figures 4(a)–(b)).

Moreover, the inconsistency between the above results may also be related to the difference in the ¹⁵N abundances of the different plant species (Swap et al. 2004; Liu and Wang 2010). As far as we know, non-environmental variables, especially plant species, can substantially influence an ecosystem’s ¹⁵N abundance, and the differences between the δ¹⁵N values of plant species can be greater than 3.20‰ (Martinelli et al. 1999; Kahmen, Wanek, and Buchman 2008). The differences in the plant δ¹⁵N values of different plant species have been confirmed in the northeastern United States (Pardo et al. 2007) and the agro-pasture area in northern China (Chen 2018; Liu et al. 2018; Zhang, Chen, and Chen 2020). For example, in the farming-pasture ecotone in northern China, the δ¹⁵N of Asiatic plantain has been found to be positively correlated with temperature, while the δ¹⁵N of Artemisia absinthium has been found to be significantly negatively correlated with temperature (Chen 2018; Liu et al. 2018). In this study, we only sampled one main plant species that grew along the entire transect, which largely separated the effect of the plant species on the foliar ¹⁵N abundance from the influences of the climate and soil factors. However, in some previous regional-scale studies, the N isotopes of many plant species were studied as an integral part of the investigation (Bai and Houlton 2009, 2009; Ma et al. 2012). Therefore, the impact of interspecies differences may modify the unified and positive correlation between the plant δ¹⁵N and temperature and may even result in a negative correlation (Pardo et al. 2007; Chen 2018; Liu et al. 2018).

For the climate variable AI, although its VIP value was greater than 0.8, its relative impact was only 3.35% (Figure 6), so it was less important in determining the plant δ¹⁵N. This may be due to the narrow range of AI values, that is, the AI only varied from 1.23 in the eastern part of the transect to 2.79 in the western part of the transect.

Importance of edaphic controls on plant δ¹⁵N

Soil factors may directly or indirectly affect dynamic N processes (e.g., the N input-output balance and N transformation) by influencing microbial activity and the accessibility of N to microbes and thereby the plant δ¹⁵N values (Menge et al. 2011; Wang et al. 2014). In this study, the overall relative effect of the edaphic variables on the plant δ¹⁵N was 45.29%, and the first two soil principal components explained 31.2% of the variance of the plant δ¹⁵N. This illustrates that the edaphic factors should be taken into account in explaining the variations in plant δ¹⁵N values. For example, the SWC, which was the second most important contributor of the edaphic factors (Table 2) and accounted for 11.01% of the total explanatory variance (Figure 6), exerted a significant effect on the soil aeration conditions, which affect the rates of N mineralization and NH₃ volatilization, as well as the microbial decomposition of organic matter and thereby the soil and plant δ¹⁵N (Butterbach-Bahl and Gundersen 2011). We argue that the negative δ¹⁵N-SWC relationship was caused by two reasons. First, in this regional-scale study, the utilization efficiency of the soil N by plants may have increased with increasing soil moisture along the transect, which resulted in more ¹⁵N-depleted soil, thus causing the soil δ¹⁵N values to be lower. Because there was a close positive correlation between the plant δ¹⁵N and soil δ¹⁵N (Figure 4(g)), the plant δ¹⁵N was negatively correlated with the SWC. However, in the relatively dry sites where the vegetation was thin, the net loss of N through NH₃ volatilization was far more likely to occur, which led to enrichment of the soil N in ¹⁵N and an increase in the plant δ¹⁵N (Frank, Evans, and Tracy 2004; Chang et al. 2009). Second, the soil δ¹⁵N often increases due to decomposition of soil organic matter. As the SWC increases, the incomplete decomposition of plant litter may bring more enriched ¹⁴N organic matter into the soil N pool, thus causing more negative plant δ¹⁵N values (Bai and Houlton 2009; Liu and Wang 2009).

Compared with the SWC, the soil N content was a more important factor controlling the plant δ¹⁵N (Table 2), and its relative influence explained approximately 27.47% of the variations in the plant δ¹⁵N (Figure 6). This indicates that the plant δ¹⁵N was more sensitive to changes in the soil N content than to changes in the other soil-related variables along the transect. The reason for this is that the plant δ¹⁵N is closely related to the TN content of the leaves, which is closely related to the soil N content (Craine et al. 2009; Yang et al. 2013). The positive and indirect effect of the soil N content on the plant δ¹⁵N may be related to climate-induced changes in the N availability for plant growth. In this study, precipitation was considered to be the most important factor regulating plant growth and the soil nutrient status (Hobbie, Macko, and Williams 2000; Falkengren-Grerup et al. 2004; Lü et al. 2012). Due to the decrease in precipitation from east to west along the transect (Figure 3(g)), the microbial activity and soil N availability may have been greater in the warm and wet sites than in the cold and dry sites (Luo et al. 2015). Under high N availability conditions, isotopically depleted N is preferentially lost from the soil through NH₃ volatilization, denitrification, and leaching, which leads to ¹⁵N enrichment in the soil and a subsequent increase in the plant δ¹⁵N.

The plant δ¹⁵N slightly increased as the soil clay content increased along the transect (Figure 4(f)). Although the soil clay content was an important variable (VIP > 0.8), it was less important in determining the plant δ¹⁵N, and its relative effect was only 3.44% (Figure 6). The increase in the plant δ¹⁵N with increasing clay content may have been driven by its indirect effects on the soil moisture and oxygen concentration, and the subsequent effects on the loss of gaseous N from the soil. Butterbach-Bahl and Gundersen (2011) reported that the loss of gaseous N from fine-textured soils can be expected to be higher due to the more frequent stimulation of the predominantly anaerobic process of denitrification. Subsequently, this large gaseous N loss may lead to a larger δ¹⁵N value in fine-textured soils because of the strong potential effects on ¹⁵N/¹⁴N fractionation.

The soil pH and C/N ratio had insignificant effects on the plant ¹⁵N abundance compared with the other edaphic variables along the transect. Our analyses revealed that the soil pH and C/N ratio together only explained approximately 3.37% of the total variations in the plant δ¹⁵N (Figure 6), and their VIP values were less than 0.8 (Table 2), which indicates that the soil pH and C/N ratio did not exert significant effects on the plant δ¹⁵N. This phenomenon may be due to the relatively uniform soil pH values and the fluctuating soil C/N ratios along the transect (Figures 3(i)–(j)). Nevertheless, we cannot speculate that the lack of obvious trends in the plant δ¹⁵N values with the soil pH and C/N ratio (Figures 4(h)–(i)) indicates that there was no significant variation in the N cycling because the soil pH values and C/N ratios along the transect were high enough for substantial gaseous N loss to occur via NH₃ volatilization (Butterbach-Bahl and Gundersen 2011; Yang et al. 2013). Therefore, the soil C/N ratio and soil pH are still important drivers of the N cycle at the regional scale.

Limitations and conclusions

Nonetheless, there are still some limitations. First, some of the soil properties, especially the chemical variables such as the soil pH, may actively respond to external environmental changes. Treating such variables as constant may lead to a biased estimation of the variable’s importance if the variable exhibits marked temporal variation. Second, the transect selected in this study was relatively narrow (approximately 1200 km), and whether some relationships between the plant δ¹⁵N and explanatory variables are universal needs to be verified further. Third, although the edaphic factors seem to be individual controls on the plant δ¹⁵N, the climate may have an effect on those edaphic factors and thus on the plant δ¹⁵N. Therefore, to better understand the environmental factors controlling the plant δ¹⁵N dynamics, future research should expand the transect range and collect data from more sites to address the above discussed limitations at the regional scale. In addition, the effects of the interactions between the climatic and edaphic variables on the ecological N cycles in temperate grasslands should be investigated. This will be important for eliminating the uncertainty in the impacts of the environmental factors on the plant δ¹⁵N and providing more robust predictions based on the roles of the various factors in the plant δ¹⁵N dynamics.

In summary, we investigated the relative degrees of importance of the climatic and edaphic factors regulating the plant δ¹⁵N along a transect across the temperate grasslands in northern China. The results of this study have important implications for understanding the mechanisms by which the plant δ¹⁵N changes. We concluded that compared with the edaphic factors, the climatic factors investigated in this study explained more of the variation in the plant δ¹⁵N along the transect. In terms of the climatic factors, the MAP exerted a first-order control on the variability of the plant δ¹⁵N, followed by the MAT and AI. For the edaphic factors, the soil N content exerted the most significant control on the plant ¹⁵N abundance, followed by the SWC, soil clay content, soil pH, and C: N ratio. These findings provide important information about the responses of the ecosystem’s N cycle to climatic and edaphic changes in the temperate grasslands in northern China. Moreover, a high-performance model for predicting the plant δ¹⁵N was created using the important independent variables (VIP > 0.8).

Highlights

• The relative importance of climatic and edaphic controls on plant δ¹⁵N was systematically studied across the temperate grasslands in northern China.

• The influences of climatic and edaphic factors on the plant δ¹⁵N were disentangled using the VIP approach.

• Compared with the edaphic factors, the climatic factors explained more of the variation in the plant δ¹⁵N along the transect.

• The important information about the responses of the ecosystem’s N cycle to climatic and edaphic changes was provided in the temperate grasslands in northern China.

Declaration of conflict interest statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the Natural Science Foundation of Hunan Province (grant number 2020JJ4296). The authors would like to thank the editors and the anonymous reviewers for their valuable comments and suggestions on this paper. We also thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This research was financed by the Natural Science Foundation of Hunan Province (No. 2020JJ4296).