Direct and indirect effects of climatic variations on the interannual variability in net ecosystem exchange across terrestrial ecosystems

Abstract

Climatic variables not only directly affect the interannual variability (IAV) in net ecosystem exchange of CO₂ (NEE) but also indirectly drive it by changing the physiological parameters. Identifying these direct and indirect paths can reveal the underlying mechanisms of carbon (C) dynamics. In this study, we applied a path analysis using flux data from 65 sites to quantify the direct and indirect climatic effects on IAV in NEE and to evaluate the potential relationships among the climatic variables and physiological parameters that represent physiology and phenology of ecosystems. We found that the maximum photosynthetic rate was the most important factor for the IAV in gross primary productivity (GPP), which was mainly induced by the variation in vapour pressure deficit. For ecosystem respiration (RE), the most important drivers were GPP and the reference respiratory rate. The biome type regulated the direct and indirect paths, with distinctive differences between forests and non-forests, evergreen needleleaf forests and deciduous broadleaf forests, and between grasslands and croplands. Different paths were also found among wet, moist and dry ecosystems. However, the climatic variables can only partly explain the IAV in physiological parameters, suggesting that the latter may also result from other biotic and disturbance factors. In addition, the climatic variables related to NEE were not necessarily the same as those related to GPP and RE, indicating the emerging difficulty encountered when studying the IAV in NEE. Overall, our results highlight the contribution of certain physiological parameters to the IAV in C fluxes and the importance of biome type and multi-year water conditions, which should receive more attention in future experimental and modelling research.

• net ecosystem exchange
• interannual variability
• climatic variations
• physiological parameters
• direct and indirect effects
• relative importance

To access the supplementary material to this article, please see Supplementary files under ‘Article Tools’.

1. Introduction

During the past century, increased human activities have caused large changes in the global climate system, such as an increase in mean air temperature, a decrease in diurnal temperature range, a redistribution of precipitation and intensification of extreme climate events (IPCC, 2013). Terrestrial ecosystems can not only respond to the climate change, but they also regulate the changing rate of the climate because they absorb 30 % of the anthropogenic CO₂ emissions and determine the interannual fluctuation in the atmospheric CO₂ concentration growth rate (Canadell et al., 2007; Le Quéré et al., 2014; Ahlström et al., 2015). The interannual variability (IAV) in the net ecosystem exchange of CO₂ (NEE) is widely observed in eddy-flux networks, and its drivers have been attributed to a series of climatic and biotic factors (Barr et al., 2007; Pintér et al., 2008; Richardson et al., 2009; Yuan et al., 2009; Dragoni et al., 2011; Humphreys and Lafleur, 2011; Wu et al., 2013). However, model simulations of the IAV in NEE is one of the most urgent and challenging scientific tasks because of the lack of substantial understanding of the intrinsic predictability of IAV in the terrestrial ecosystem carbon (C) cycle (Luo et al., 2015).

The difficulty in predicting the IAV in NEE might result partially from the failure of models to represent the temporal evolution of ecophysiological (e.g. physiological and phenological) factors, which outweigh the climatic factors in most eddy-flux sites at the interannual scale (Shao et al., 2015). Although empirical and semi-empirical research suggests that the ecophysiological properties of ecosystems should be taken into consideration for better simulation of the IAV in C fluxes (Kuzyakov and Gavrichkova, 2010; Migliavacca et al., 2011; Xia et al., 2015), the results from modelling studies indicate that it is still troublesome to predict the temporal changes in the biotic factors or ecophysiological properties of ecosystems (Richardson et al., 2012). Therefore, appropriately representing the dynamics of ecosystem ecophysiology embraces the opportunity to address the issues raised by the IAV in NEE.

The ecophysiological factors that are directly related to the IAV in NEE include photosynthetic and respiratory parameters (hereafter referred as the physiological parameters) that characterize the responses of photosynthesis and respiration to the climatic variations (Shao et al., 2015). Other biotic factors (e.g. biotic disturbance) can only affect the IAV in NEE by changing these parameters. The changes in the physiological parameters can be induced by climatic variations, internal changes in the ecosystem and natural or human disturbances (Richardson et al., 2007; Shao et al., 2015). Among these, the effects of climatic factors on physiological parameters reflect the ecophysiological regulation of ecosystem C cycling by the climatic variations, which occur from the leaf to the ecosystem level (Niu et al., 2012; Wythers et al., 2013; Way and Yamori, 2014; Yamori et al., 2014). At the ecosystem level, a handful of studies have shown that the seasonality of physiological parameters, such as the maximum carboxylation rate (V _cmax, maximum velocity of the carboxylase), maximum photosynthetic rate (A _m , light-saturated photosynthetic rate at the ambient CO₂ concentration), light use efficiency, water use efficiency (WUE), reference respiratory rate at 10 °C (R ₁₀) and temperature sensitivity (Q ₁₀) were related to radiation, temperature and water conditions in some ecosystems (Aber et al., 1996; Barr et al., 2007; Ricciuto et al., 2008; Yu et al., 2008; Ju et al., 2010; Shao et al., 2014; Shi et al., 2014; Peichl et al., 2015). At a larger temporal scale (interannual), the optimum temperature of NEE, the Michaelis constant of gross primary productivity (GPP), WUE and R ₁₀ have also been related to the climatic variations (Bahn et al., 2008; Niu et al., 2011, 2012; Shao et al., 2014).

The effects of climate on the IAV in C fluxes mediated by physiological parameters can be considered to be the indirect climatic effects, the strength of which depends on both the correlations between climatic variables and physiological parameters, and the correlations between physiological parameters and C fluxes. If the indirect climatic effects are important contributors to the IAV in C fluxes, the performance of a model can be improved by expressing the physiological parameters as functions of the climatic variables (Ricciuto et al., 2008; Shao et al., 2014). Whether this approach will work or not largely depends on the explicit paths through which the climate directly and indirectly impacts the C fluxes. However, these paths have not been synthetically investigated within eddy-flux networks. To address this issue, we applied a path analysis to 481 site-years of flux data from 65 sites in order to examine the direct and indirect effects of climatic variations on the IAV in NEE and its components (GPP and ecosystem respiration (RE)). Path analysis is a statistical method that is able to quantify the direct and indirect dependencies between or among many variables under consideration (i.e. climatic variables, physiological parameters and C fluxes in this study). The main questions we focused on were: (1) how and to what extent does the climate directly and indirectly affect IAV in C fluxes; and (2) how are these relationships affected by ecosystem characteristics, such as biome types and multi-year water conditions?

2. Materials and methods

2.1. Data sources

The original half-hourly C flux and auxiliary climate data were obtained from AmeriFlux (www.ameriflux.lbl.gov/), CarboEurope (www.europe-fluxdata.eu) and ChinaFLUX (www.chinaflux.org). Only the sites with data for ≥5 yr were selected. In total, there were 481 site-years of data belonging to 65 eddy covariance measurement sites from 1992 to 2010 (Table 1, Supplementary Fig. 1). These sites covered a large range of geographical positions and climatic conditions, which can be broadly classified into evergreen needleleaf forests, deciduous broadleaf forests, evergreen broadleaf forests, mixed forests, shrublands, grasslands and croplands (Table 1). The original data included CO₂ flux (Fc), friction velocity, photosynthetically active radiation (PAR) or global radiation, air temperature (Ta), soil temperature (Ts), annual precipitation (PPT), relative humidity (RH), vapour pressure deficit (VPD) and latent heat flux (LE).

Table 1. Information of study sites

Site	Lat (°)	Long (°)	Alt (m)	Eco	n (years)	Ta (°C)	PPT (mm)	References
US-Blo	38.90	−120.63	1315	ENF	9	11.9	1073	Goldstein et al. (2000)
US-Dk3	35.98	−79.09	163	ENF	8	15.0	1050	Stoy et al. (2008)
US-Fmf	35.14	−111.73	2160	ENF	5	9.7	521	Dore et al. (2008)
US-Fuf	35.09	−111.76	2180	ENF	5	8.9	497	Dore et al. (2008)
US-Ho1	45.20	−68.74	60	ENF	13	6.7	840	Richardson et al. (2009)
US-Ho2	45.21	−68.75	91	ENF	11	6.6	842	Xiao et al. (2005)
US-Ho3	45.21	−68.73	61	ENF	5	6.5	877	Davidson et al. (2006)
US-Me2	44.45	−121.56	1253	ENF	9	7.4	441	Thomas et al. (2009)
US-NC2	35.80	−76.67	12	ENF	5	15.8	1212	Noormets et al. (2010)
US-NR1	40.03	−105.55	3050	ENF	10	2.4	663	Blanken et al. (2009)
US-SP3	29.75	−82.16	36	ENF	11	20.1	1034	Powell et al. (2008)
DE-Tha	50.96	13.57	380	ENF	12	7.6	820	Prescher et al. (2010)
DE-Wet	50.45	11.46	703	ENF	7	6.6	967	Rebmann et al. (2010)
ES-ES1	39.35	−0.32	1	ENF	5	17.5	495	Reichstein et al. (2005)
FI-Hyy	61.85	24.29	185	ENF	11	4.8	702	Suni et al. (2003)
FI-Sod	67.36	26.64	189	ENF	7	1.1	425	Aurela (2005)
FR-Pue	43.74	3.60	211	ENF	8	13.6	904	Allard et al. (2008)
IT-Cpz	41.71	12.38	9	ENF	8	15.0	729	Garbulsky et al. (2008)
NL-Loo	52.17	5.74	25	ENF	10	10.3	843	Dolman et al. (2002)
RU-Fyo	56.46	32.92	265	ENF	9	5.2	762	Milyukova et al. (2002)
SE-Nor	60.09	17.48	28	ENF	5	6.4	486	Lagergren et al. (2005)
CN-Qia	26.74	115.06	100	ENF	7	18.3	1139	Zhang et al. (2011)
US-Bar	44.06	−71.29	272	DBF	5	7.4	1330	Jenkins et al. (2007)
US-Dk2	35.97	−79.10	168	DBF	7	15.2	1097	Stoy et al. (2008)
US-Ha1	42.54	−72.17	340	DBF	15	8.5	1153	Urbanski et al. (2007)
US-MMS	39.32	−86.41	275	DBF	5	12.7	1186	Dragoni et al. (2011)
US-MOz	38.74	−92.20	219.4	DBF	5	13.0	1035	Yang et al. (2010)
US-Slt	39.91	−74.60	30	DBF	5	12.4	1032	Hollinger et al. (2010)
US-Syv	46.24	−89.35	540	DBF	5	4.5	696	Desai et al. (2005)
US-UMB	45.56	−84.71	234	DBF	6	7.2	412	Hardiman et al. (2011)
US-WCr	45.91	−90.08	515	DBF	7	5.7	744	Cook et al. (2004)
DE-Hai	51.08	10.45	433	DBF	8	8.3	803	Kutsch et al. (2008)
FR-Hes	48.67	7.07	293	DBF	12	10.4	918	Granier et al. (2000)
IT-Ro1	42.41	11.93	153	DBF	7	15.5	769	Rey et al. (2002)
CN-Din	23.17	112.53	240	EBF	6	20.1	1481	Zhang et al. (2010)
BE-Bra	51.31	4.52	15	MF	9	11.4	792	Carrara et al. (2003)
BE-Vie	50.31	6.00	491	MF	11	8.3	719	Aubinet et al. (2001)
CH-Lae	47.48	8.37	645	MF	5	8.3	605	Eugster et al. (2007)
DK-Sor	55.49	11.64	35	MF	13	8.7	656	Pilegaard et al. (2001)
CN-Cha	42.40	128.10	738	MF	6	3.9	498	Guan et al. (2006)
US-Los	46.08	−89.98	480	SHR	6	3.7	716	Sulman et al. (2009)
MX-Lpa	24.13	−110.44	21	SHR	5	23.4	137	Bell et al. (2012)
US-Aud	31.59	−110.51	1469	GRA	6	15.9	366	Krishnan et al. (2012)
US-Bkg	44.35	−96.84	510	GRA	5	6.3	540	Gilmanov et al. (2005)
US-Dk1	35.97	−79.09	168	GRA	7	15.1	1022	Stoy et al. (2008)
US-FPe	48.31	−105.10	634	GRA	8	5.7	391	Gilmanov et al. (2005)
US-Fwf	35.45	−111.77	2270	GRA	5	8.4	562	Dore et al. (2008)
US-SRM	31.82	−110.87	1116	GRA	7	19.0	324	Scott et al. (2009)
US-Ton	38.43	−120.97	177	GRA	9	16.9	504	Ma et al. (2007)
US-Var	38.41	−120.95	129	GRA	9	15.9	493	Xu and Baldocchi (2004)
US-Wkg	31.74	−109.94	1531	GRA	6	17.2	274	Scott et al. (2010)
AT-Neu	47.12	11.32	1095	GRA	7	6.9	731	Wohlfahrt et al. (2008)
DE-Gri	50.95	13.51	385	GRA	5	7.8	901	Prescher et al. (2010)
DK-Lva	55.68	12.08	15	GRA	5	9.0	655	Gilmanov et al. (2007)
HU-Bug	46.69	19.60	106	GRA	6	10.3	499	Nagy et al. (2007)
NL-Ca1	51.97	4.93	1	GRA	6	11.0	759	Gilmanov et al. (2007)
CN-Dan	29.67	91.33	4333	GRA	5	2.47	470	Xu et al. (2007)
CN-Hab	37.62	101.28	3220	GRA	6	−0.8	484	Kato et al. (2006)
US-ARM	36.61	−97.49	314	CRO	7	15.3	658	Fischer et al. (2007)
US-Bo1	40.01	−88.29	219	CRO	11	11.6	748	Hollinger et al. (2005)
US-Ne1	41.17	−96.48	361	CRO	8	10.7	846	Suyker et al. (2004)
US-Ne2	41.16	−96.47	362	CRO	8	10.2	872	Suyker et al. (2004)
BE-Lon	50.55	4.74	165	CRO	5	11.2	611	Moureaux et al. (2006)
DE-Geb	51.10	10.91	157	CRO	6	9.9	503	Anthoni et al. (2004)
DE-Kli	50.89	13.52	478	CRO	6	7.7	951	Prescher et al. (2010)

Lat, latitude; Long, longitude; Alt, altitude; Eco, ecosystem type; n, data record length; Ta, mean annual air temperature; PPT, annual precipitation; ENF, evergreen needleleaf forest; DBF, deciduous broadleaf forest; EBF, evergreen broadleaf forest; MF, mixed forest; SHR, shrubland; GRA, grassland; CRO, cropland.

2.2. Data processing

To obtain the annual values for the C fluxes, climatic variables and physiological parameters, several data processing procedures were applied. The original half-hourly Fc data were first spike-screened by using the median of absolute deviation as a spike estimator based on the double differenced time series (Papale et al., 2006). Next, the night-time Fc was rejected if the friction velocity was lower than a threshold, which was determined by a 99 % threshold criterion method based on each year's night-time data (Reichstein et al., 2005; Papale et al., 2006). After that, the artificial neural network model with a feed-forward back-propagation algorithm and sigmoid transfer function was applied to partition the NEE into GPP and RE, and to fill the gaps (Shao et al., 2015). The artificial neural network model was trained and validated separately for each year's daytime and night-time data. The input variables of the night-time model were Ts, RH and four seasonal indices, which represented the four seasons (Papale and Valentini, 2003), and the output variable was the night-time RE. The hidden layer had six nodes. The dataset was randomly divided into a training dataset (70 %) and a validation dataset (30 %). This process was repeated 20 times, and the median of 20 simulated outputs was used to fill the gaps in night-time RE and to simulate the daytime RE. The input variables of the daytime model were PAR, Ta, RH and the four seasonal indices, and the output variable was GPP (daytime RE – NEE). Six nodes were set in the hidden layer. The training and validation procedures were repeated 20 times and the median of the simulated GPP was used to fill the GPP gaps. The NEE gaps were filled by the gap-filled RE – GPP. Finally, the gap-filled NEE, GPP and RE, as well as the corresponding climatic variables, were integrated or averaged on an annual basis.

The physiological (photosynthetic and respiratory) parameters used in this study were the maximum photosynthetic rate (A _m ), apparent quantum yield (α), WUE, R ₁₀ and Q ₁₀. A _m , α, R ₁₀ and Q ₁₀ were derived from the Michaelis–Menten and Q ₁₀ equations (Falge et al., 2001; Richardson and Hollinger, 2005):1 2

These equations were parameterised on the basis of 30-d moving windows with a step of 1 d, and only the GPP and RE partitioned from the NEE data were used. The estimated parameters with coefficients of variation larger than 0.5 were excluded (Reichstein et al., 2005) and the median of a certain parameter in a year was considered to be the annual value of the parameter. The WUE is the ratio of GPP to evapotranspiration (ET, calculated from LE [ET (mm/day)=LE (W/m²)/28.94]) (Yuan et al., 2009).

2.3. Path analysis

The dataset used in this study had previously been analysed by Shao et al. (2015) but with different methods and aims. The additive model in Shao et al. (2015) aimed to quantify the relative contribution of climatic and biotic factors to IAV in C fluxes, but it could not identify the paths through which the climate directly and indirectly affected the C fluxes. Therefore, in order to elucidate the direct and indirect effects of climatic variables on the IAV in C fluxes, we applied path analysis, which is an extension of regression model and can deal with multiple casual relationships between or among correlated variables (Shipley, 2004), to the dataset in Shao et al. (2015). Theoretically, path analysis can incorporate nonlinear interactions between climatic variables and C fluxes by using nonlinear functional forms. However, in our study, no nonlinear relationships between variables were found at the interannual scale, though they might exist at shorter temporal scales. The standardised path coefficient (ρ, the path coefficient when all the variables were standardised before the analysis), which is an analogy of the correlation coefficient, was used to quantify the effect size of one variable on another.

Shao et al. (2015) suggested that the correlations between climatic and biotic factors might depend on biome types and multi-year water conditions. Therefore, the study sites were grouped into evergreen needleleaf forests, deciduous broadleaf forests, grasslands and croplands according to their plant function type. In addition, we used the water balance index (WBI), which is the difference between evapotranspiration and precipitation (Law et al., 2002), to group the sites into wet (WBI below −500 mm), moist (−500 mm<WBI<0 mm) and dry (WBI>0 mm) ecosystems according to their multi-year water conditions. In order to exclude the influence of different magnitudes of the variables among ecosystems, the annual climatic variables, physiological parameters and C fluxes in each site were first standardised by subtracting the mean and then dividing it by the standard deviation of the multi-year data. Path analysis was then applied to each group separately.

For each ecosystem group, two path models were conducted: one for NEE and another for GPP and RE. The use of two models instead of one is because even though some variables are critical to GPP and RE, they may not necessarily be important to NEE. The potential paths in the models mainly included three aspects: the direct effects of climatic variables on C fluxes, the direct effects of climatic variables on physiological parameters and the direct effects of physiological parameters on C fluxes. The first of these three groups were the direct climatic effects; the latter two together characterised the indirect climatic effects on the IAV in C fluxes. In addition, we also considered the regulation of GPP on RE (Supplementary Tables 1 and 2).

The potential paths were based on the following considerations (Supplementary Tables 1 and 2). (1) For the causal paths from climate to C fluxes, five climatic variables (PAR, Ta, Ts, VPD and WBI) were considered to affect NEE and GPP directly, whereas Ta, Ts, VPD and WBI affect RE (precipitation was not considered as an explanatory variable because its effects were actually mediated by water availability, which was the opposite of WBI). (2) The IAV in photosynthetic parameters (A _m , α and WUE) may be induced by changes in radiation (PAR), temperature (Ta, Ts) and water conditions (VPD, WBI), whereas the IAV in respiratory parameters (R ₁₀ and Q ₁₀) might be driven by temperature and water conditions. (3) The photosynthetic and respiratory parameters may influence GPP and RE, respectively, and both kinds of physiological parameter may be important to NEE. (4) The correlation between GPP and RE was found within and across ecosystems (Chen et al., 2015), and may be caused by the substrate supply of photosynthesis to both autotrophic and heterotrophic respiration (Schimel et al., 1994). Therefore, we also assumed that GPP might influence R ₁₀. However, autotrophic and heterotrophic respiration might be driven by different factors and correlated differently with GPP, which might result in relationships other than that between GPP and R ₁₀. Thus, the control of GPP on RE was also included in the path model. (5) Water availability in autumn and winter might influence the plants’ phenology in the next spring and thus was important to the C fluxes in the next year (Gordo and Sanz, 2005). Therefore, the lag effects of autumn–winter WBI on next year's C fluxes and physiological parameters were considered.

The full path models contained too many paths compared to the sample size, which might have caused an overparameterisation problem or may have even failed to derive any statistical tests. Therefore, we simplified the model structure by excluding the paths where the corresponding bivariate correlation coefficients (r) were larger than −0.1 and smaller than 0.1. A path model with the remaining paths was conducted, and the standardised path coefficients (ρ) and the corresponding p-values were calculated via maximum likelihood estimation (Beaujean, 2014). To obtain the most parsimonious path model for each ecosystem group, the non-significant paths were excluded recursively by dropping the path with the largest non-significant p-value each time and then re-parameterising the new model. This procedure was repeated until all the paths were significant (p<0.05). In order to make the results among ecosystem groups comparable, we constructed final path models that contained all the significant paths in the most parsimonious models of all the ecosystem groups to be compared. This treatment helped to eliminate the uncertainty caused by variable selection. The path analysis was achieved by using the package lavaan in R software (Rosseel, 2012; R Core Team, 2014).

3. Results

3.1. Direct and indirect climatic effects in different biomes

In all the biomes being studied (evergreen needleleaf forests, deciduous broadleaf forests, grasslands and croplands), the maximum photosynthetic rate (A _m ) was the most important factor for the IAV in C fluxes (the standardised path coefficient (ρ) as an analogy of the correlation coefficient fell between −0.32 and −0.68 for NEE and 0.72–0.85 for GPP), which was mainly induced by the variation in VPD (Figs. 1 and 2). For the IAV in RE, the most important drivers were GPP (ρ=0.47–0.67) and reference respiratory rate (R ₁₀, ρ=0.24–0.29 except for croplands; Fig. 2). Though we also considered the lag effects of autumn–winter WBI on the physiological parameters and C fluxes in the next year (Supplementary Table 3), they were not selected in the final path models.

Fig. 1 The direct and indirect effects of climatic variation on the interannual variability (IAV) in net ecosystem exchange (NEE) for evergreen needleleaf forests (a), deciduous broadleaf forests (b), grasslands (c) and croplands (d). Green and red lines represent positive and negative effects, respectively. Solid and dashed lines represent significant (p<0.05) and non-significant (p≥0.05) effects, respectively. The values beside the paths are the standardised (0–1) path coefficients (ρ), which are only shown for significant effects. PAR, photosynthetically active radiation; Ta, air temperature; Ts, soil temperature; VPD, vapour pressure deficit; WBI, water balance index; A _m, maximum photosynthetic rate; α, apparent quantum yield; WUE, water use efficiency; R ₁₀, reference respiratory rate at 10 °C; Q ₁₀, temperature sensitivity.

Fig. 2 The direct and indirect effects of climatic variation on in the interannual variability (IAV) in gross primary productivity (GPP) and ecosystem respiration (RE) for evergreen needleleaf forests (a), deciduous broadleaf forests (b), grasslands (c) and croplands (d). Green and red lines represent positive and negative effects, respectively. Solid and dashed lines represent significant (p<0.05) and non-significant (p≥0.05) effects, respectively. The values beside the paths are the standardised (0–1) path coefficients (ρ), which are only shown for significant effects. PAR, photosynthetically active radiation; Ta, air temperature; Ts, soil temperature; VPD, vapour pressure deficit; WBI, water balance index; A _m, maximum photosynthetic rate; α, apparent quantum yield; WUE, water use efficiency; R ₁₀, reference respiratory rate at 10 °C; Q ₁₀, temperature sensitivity.

There were also different direct and indirect paths from climate to C fluxes among biomes. For example, a positive correlation between temperature and apparent quantum yield (α), as well as a negative correlation between temperature sensitivity (Q ₁₀) and RE, was only found in forests, while the controls of GPP on reference respiratory rate (R ₁₀) only appeared in non-forests (Fig. 2). Compared with deciduous broadleaf forests, the evergreen needleleaf forests showed direct climatic effects on C fluxes (PAR on GPP, and Ts on RE), while there were more significant indirect climatic effects in deciduous broadleaf forests than in evergreen needleleaf forests (Fig. 2a and b). Although PAR influenced the maximum photosynthetic rate in both grasslands (ρ=−0.39) and croplands (ρ=0.34) within non-forests, it did so in opposite directions. In addition, the indirect climatic effects on RE in grasslands were not found in croplands (Fig. 2c and d).

3.2. Direct and indirect climatic effects in ecosystems with different aridity

The direct and indirect effects of climate on IAV in C fluxes were also different among ecosystems with different multi-year water conditions (Fig. 3). Overall, the moist ecosystems had more complicated direct and indirect paths than wet and dry ecosystems. Compared with wet ecosystems, the IAV in GPP of moist and dry ecosystems can be induced by apparent quantum yield (α, an indicator of light use efficiency, ρ=0.08 and 0.17, respectively) and WUE (ρ=0.13 and 0.24, respectively), and the IAV in RE can be induced by reference respiratory rate (R ₁₀, ρ=0.31 and 0.25, respectively). On the other hand, direct effects of WBI on RE were found in moist (ρ=−0.10) and dry (ρ=−0.18), but not wet ecosystems.

Fig. 3 The direct and indirect effects of climatic variation on the interannual variability (IAV) in net ecosystem exchange (NEE, a–c), gross primary productivity (GPP) and ecosystem respiration (RE, d–f) for wet (water balance index, or WBI, below −500 mm; a, d), moist (−500 mm<WBI<0 mm; b, e), and dry (WBI>0 mm; c, f) ecosystems. Note that a larger WBI indicates a drier environment. Green and red lines represent positive and negative effects, respectively. Solid and dashed lines represent significant (p<0.05) and non-significant (p≥0.05) effects, respectively. The values beside the paths are the standardised (0–1) path coefficients (ρ), which are only shown for significant effects. PAR, photosynthetically active radiation; Ta, air temperature; Ts, soil temperature; VPD, vapour pressure deficit; WBI, water balance index; A _m, maximum photosynthetic rate; α, apparent quantum yield; WUE, water use efficiency; R ₁₀, reference respiratory rate at 10 °C; Q ₁₀, temperature sensitivity.

3.3. Model performance

The results derived from the NEE model were generally consistent with those from the GPP and RE model, but with some discrepancies (Figs. (1–3)). For example, the complicated indirect climatic effects in grasslands seen in the GPP and RE model (Fig. 2c) were largely simplified in the NEE model (Fig. 1c). In the dry ecosystems, Ts emerged to directly influence NEE, though it did not affect GPP or RE (Fig. 3e and f). Overall, the path models were better for explaining the IAV in GPP (r ²=0.54–0.81) and RE (r ²=0.36–0.60) than that in NEE (r ²=0.17–0.61; Figs. (1–3)). Among ecosystem groups, the model performed best in grasslands and croplands (the r ² values of C fluxes were all above 0.5) and worst in evergreen needleleaf forests (the r ² values of NEE were below 0.2). The climatic variables explained the variations in the physiological parameters best in grasslands, with the r ² of maximum photosynthetic rate being 0.31. In other ecosystems, climatic variations could only induce very small levels of the IAV in the physiological parameters (Figs. (1–3)).

Because the physiological parameters were not directly measured but derived from equations, the potential correlations among them might cause some uncertainty. Although significant correlations were found in many pairs of parameters, the coefficients were generally very small, except for those between maximum photosynthetic rate and WUE in grasslands and croplands and dry ecosystems (r>0.4, Table 2). In addition, there was no significant correlation between reference respiratory rate and temperature sensitivity (Table 2).

Table 2. Correlation coefficients (r) of the physiological parameters in different ecosystem types

Ecosystem group	r(A m , α)	r(A m , WUE)	r(α, WUE)	r(R 10, Q 10)
ENF	0.21*	0.10	0.11	0.14
DBF	0.10	0.11	0.28*	0.07
GRA	0.29*	0.41*	0.27*	0.07
CRO	0.02	0.56*	0.10	0.07
Wet	0.05	0.11	0.02	0.10
Moist	0.15*	0.29*	0.19*	0.10
Dry	0.09	0.49*	0.09	0.14

ENF, evergreen needleleaf forests; DBF, deciduous broadleaf forests; GRA, grasslands; CRO, croplands; Wet, ecosystems with a water balance index (WBI) below −500 mm; Moist, ecosystems with −500 mm<WBI<0 mm; Dry, ecosystems with WBI>0 mm; A _m, maximum photosynthetic rate; α, apparent quantum yield; WUE, water use efficiency; R ₁₀, reference respiratory rate at 10 °C; Q ₁₀, temperature sensitivity.

*p<0.05.

4. Discussion

4.1. Regulation of biome type

The climatic variations can drive the IAV in C fluxes directly and indirectly by changing the physiological parameters, with both shared and unique paths among ecosystem types. In all the biome types being studied (i.e. evergreen needleleaf forests, deciduous broadleaf forests, grasslands and croplands), the importance of maximum photosynthetic rate (A _m ) was highlighted (Figs. 1 and 2), a result that has also been found in a tundra and a fen (Humphreys and Lafleur, 2011). Xia et al. (2015) also showed that the seasonal maximal photosynthetic capacity was the best indicator of IAV in GPP for all the biomes. The control of VPD on A _m at the interannual scale (Figs. 1 and 2) reflected the control of water stress on stomatal conductance (Medlyn et al., 2001). For the IAV in RE, there was no consistent driver across all the biome types except for the strong control of GPP (Fig. 2). The tight correlation between GPP and RE might be derived from the influence of GPP on the components of RE: autotrophic and heterotrophic respiration. The relationship between GPP and autotrophic respiration was mainly caused by the proportional changes in growth respiration to GPP (Amthor, 2000). The heterotrophic respiration can be correlated with GPP because the microbes prefer to decompose the younger organic matter derived from more recent photosynthate (Schimel et al., 1994).

However, there were some distinctive differences between forests and non-forests (Fig. 2). For example, in forests, the warmer environment stimulated GPP by promoting the apparent quantum yield (α), which is physiologically defined as the ratio of photosynthesis to the absorbed PAR and is equivalent to the concept of light use efficiency, which has been applied in light use efficiency models (Johnson and Barber, 2003). This might be because the enzyme activity or leaf area index was enhanced in warmer environments (Shi et al., 2014; Gitelson and Gamon, 2015). The negative correlation between temperature sensitivity (Q ₁₀) and RE found in forests seemed to contradict common sense. However, previous studies confirmed that Q ₁₀ generally declined with the increase in temperature (Tjoelker et al., 2001), which, in turn, indicated a negative relationship between Q ₁₀ and RE because both autotrophic and heterotrophic respiration were positively correlated with temperature.

Compared with forests, the non-forest ecosystems showed the control of GPP on reference respiratory rate (R ₁₀). This might be caused by the different transportation rates of photosynthetic production between forests and non-forests. The carbohydrate synthesised in leaves usually takes a couple of days to arrive at the roots in trees but the duration is shorter than 1 d in herbs (Högberg et al., 2001; Högberg and Read, 2006; Kuzyakov and Gavrichkova, 2010). In addition, the maintenance respiration might account for the larger proportion of autotrophic respiration in trees compared with herbs, which will weaken the relationships between photosynthesis and respiration (DeLucia et al., 2007).

Previous studies have suggested that the deciduous broadleaf forests are more sensitive to climatic variations than evergreen needleleaf forests (Welp et al., 2007; Yuan et al., 2009). Our results showed that the C fluxes in deciduous broadleaf forests were more sensitive to indirect climatic effects than those in evergreen needleleaf forests, while there were direct climatic effects on GPP and RE in evergreen needleleaf forests but not in deciduous broadleaf forests (Fig. 2a and b). These differences might be due to the different leaf life spans, phenology and the environment being adapted to. Compared to deciduous broadleaf forests, evergreen needleleaf forests are usually with longer foliage longevity, lower specific leaf area and are adapted to unfavourable environments, which make the nutrient (e.g. nitrogen and phosphorus) concentrations and thus the photosynthetic parameters more constant in a changing environment (Bloom et al., 1985; Givnish, 2002). On the other hand, the foliage nutrient concentrations would be more variable and correlated with climatic variables in deciduous broadleaf forests, as the leaf development is largely controlled by climate (Polgar and Primack, 2011). Likewise, the annual litter fall with varying litter quality (e.g. C: N) in deciduous broadleaf forests might contribute to the indirect climatic effects on RE.

The most significant difference between grasslands and croplands was that the effect of PAR on maximum photosynthetic rate (A _m ) was negative in the former but positive in the latter (Fig. 2c and d). In grasslands, the relatively low productivity and high water stress resulted in excess energy when the radiation increased, damaging the photosynthetic system and reducing A _m (Schulze et al., 2005). In contrast, high productivity and irrigation allowed the croplands to enhance A _m when the radiation increased. Another difference was that the indirect climatic effects on RE that occurred in grasslands was not found in croplands, which might be because the management (i.e. irrigation, fertilisation) in croplands hampered the influence of climate.

4.2. Regulation of multi-year water conditions

Our previous study suggested that the multi-year water conditions might be the most important regulator of the biotic responses to climatic variation (Shao et al., 2015). Current study also showed that the direct and indirect effects of climate on IAV in C fluxes were different among ecosystem groups with different aridity (Fig. 3). For example, WUE and apparent quantum yield (α, an indicator of light use efficiency) had significant effects on GPP in moist and dry, but not wet ecosystems. The discrepancy in importance of WUE among different ecosystem groups might be due to the different partitioning of evapotranspiration into evaporation and transpiration. The WUE in wet ecosystems was generally larger than that in other ecosystems due to the stronger stimulation of GPP compared to evapotranspiration when water availability increased (Xiao et al., 2013; Xue et al., 2015; Zhu et al., 2015). Higher WUE means that relatively more water is lost through transpiration, which is tightly coupled with photosynthesis. Therefore, the variation of WUE in wet ecosystems is smaller than that in other ecosystems, which ultimately results in the non-significant relationship between WUE and GPP. Likewise, light use efficiency also reaches its maximum when there is no water stress (Yuan et al., 2014), and this results in a non-significant relationship between light use efficiency and GPP. Interestingly, the water condition (i.e. VPD) was important to IAV in C fluxes, even for wet ecosystems (Fig. 3b). Due to flux partitioning and gap-filling, VPD can be related to C fluxes because it was influenced by air temperature and relative humidity. These are input variables of the artificial neural network model, and our results showed that it had little direct effect on C fluxes. This suggests that the effects of VPD on GPP in wet ecosystems were not a by-product of data processing, but a reflection of the underlying causal links. This might be due to a seasonal dry environment or an extreme drought year occurring in wet ecosystems, which can highlight the importance of water conditions (Tian et al., 2000; Lewis et al., 2011).

The moist ecosystems had more complicated indirect paths than wet and dry ecosystems, which might be because the former had a longer time lag (8–10 months) between the water balance and C fluxes compared with the latter ones (Vicente-Serrano et al., 2013). This lag effect might reflect the influence of autumn–winter precipitation stored in soil on the C fluxes of the following year. In order to investigate this issue, we added the effects of autumn–winter WBI on next year's physiological parameters and C fluxes. Although some significant correlations were found (Supplementary Table 3), none of these lag effects were selected in the final path models, suggesting that the lagged climatic effects were much weaker than the concurrent ones. Other studies found even longer lags (2–4 yr) between water conditions and C fluxes (Dunn et al., 2007; Ito, 2011). However, it is difficult to evaluate these long-term lags because of insufficient data.

4.3. Uncertainty and implications

Our study showed how climatic variations may directly and indirectly regulate the IAV in C fluxes in different biomes and ecosystems with different water conditions. However, there might be some uncertain patterns in the results. Overall, the explanatory ability of the path models was best for the IAV in GPP, followed by those in RE and NEE (Figs. (1–3)). The difficulty of simulating IAV in RE and NEE might be due to the complex responses of both fluxes to climatic variations. For example, RE consists of autotrophic and heterotrophic respiration and each of these components includes multiple processes with different responses to environmental changes (Kuzyakov and Gavrichkova, 2010; Thornley, 2011; Zhou et al., 2014). Likewise, NEE is the difference between RE and GPP, whose responses to climate are quite different (Anderson-Teixeira et al., 2011; Shi et al., 2014).

Another source of uncertainty might be the autocorrelations between variables because they were not independently observed but were derived from a few observed variables (e.g. NEE and climatic variables). The correlation between GPP and RE was thought to be spurious because these were partitioned from NEE data (Vickers et al., 2009). However, recent studies suggested that this spurious correlation was relatively small (Lasslop et al., 2010; Baldocchi et al., 2015). The physiological parameters were also not independent because they were derived from the Michaelis–Menten and Q ₁₀ equations. Fortunately, these correlations were not strong enough to introduce large uncertainty into the results (Table 2). Strictly, the physiological parameters might be just by-products of parameterisation process, which would lead to spurious relations in the path models. In this situation, there should be stronger relationships between climatic variables and C fluxes than those between climate and physiological parameters. However, we found that climatic variables were more highly correlated with physiological parameters rather than C fluxes, indicating that there was no obvious spurious correlation problem. Furthermore, it is conceptually reasonable to treat the physiological parameters as drivers rather than by-products of parameterisation because the physiological parameters determine the capacities of photosynthesis and respiration, which are also widely adopted in ecosystem models. In addition, although some relationships between climatic variables, physiological parameters and C fluxes could be expected (e.g. the relationship between PAR and photosynthetic parameters, and that between temperature and RE) due to the data processing procedure, they did not appear in all the ecosystem groups, indicating that the methodological uncertainty was not the main source of the correlations between variables.

Nevertheless, some implications and suggestions can be made for future studies. First, the weak correlations between climatic variables and physiological parameters may reflect the fact that NEE is not a biochemical process, so it is difficult to characterise its biotic drivers, or may indicate that the internal changes of ecosystems, and natural or human disturbance might be the primary drivers of the IAV in physiological parameters. These drivers may include the morphology, biochemistry and phenology of leaves, soil properties, nutrient availability, vegetation characteristics and human management (Buchmann and Schulze, 1999; Cook et al., 2008; de Beeck et al., 2010; Zhou et al., 2013; Madani et al., 2014; Duchemin et al., 2015; Gitelson and Gamon, 2015; Yan et al., 2015; Zhu et al., 2015). However, great caution should be applied when extrapolating findings from the seasonal scale up to the interannual scale, because the relationships between climate and physiological parameters might change across temporal scales. For example, Tan et al. (2015) found that VPD dominated WUE in a tropical rainforest at the instantaneous scale but became the second important factor at the seasonal scale and was not correlated to WUE at the interannual scale.

Second, the patterns of the direct and indirect climatic effects on the IAV in C fluxes differed among ecosystem types. However, identifying the different underlying mechanisms of C cycle dynamics and evaluating their importance to the IAV in C fluxes will require quantitative data describing the important features of the different ecosystems. For example, tree growth might be an important aspect for forest ecosystems and leaf phenology might regulate photosynthesis in deciduous broadleaf forests (Desai et al., 2005; Dragoni et al., 2011). Species composition may shift in the grasslands across several years (Kahmen et al., 2002), whereas human management regimens are critical to croplands (Lobell et al., 2006). In order to incorporate this information into C cycle models, long-term observations should be carried out and appropriate procedures that can better weigh the importance of these data should be developed in future studies.

5. Acknowledgements

This study was financially supported by the National Natural Science Foundation of China (Grant No. 31370489), the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning and the ‘Thousand Young Talents’ Program in China. The authors also thank Marc Aubinet, Dennis Baldocchi, Christian Bernhofer, Gil Bohrer, Paul Bolstad, Nina Buchmann, Alexander Cernusca, Reinhart Ceulemans, Kenneth Clark, Ankur Desai, Allen Goldstein, Andre Granier, David Hollinger, Gabriel Katul, Alexander Knohl, Werner Kutsch, Beverly Law, Anders Lindroth, Timothy Martin, Tilden Meyers, Eddy Moors, William Munger, Kimberly Novick, Ram Oren, Serge Rambal, Corinna Rebmann, Andrew Richardson, Maria-Jose Sanz, Russ Scott, Lise Soerensen, Riccardo Valentini and Timo Vesala for their generous permission to access the flux data.

Notes

To access the supplementary material to this article, please see Supplementary files under ‘Article Tools’.