Multivariate Conditional Granger Causality Analysis for Lagged Response of Soil Respiration in a Temperate Forest

Ecological multivariate systems offer a suitable data set on which to apply recent advances in information theory and causality detection. These systems are driven by the interplay of various environmental factors: meteorological and hydrological forcing, which are often correlated with each other at different time lags; and biological factors, primary producers and decomposers with both autonomous and coupled dynamics. Here, using conditional spectral Granger causality, we quantify directional causalities in a complex atmosphere-plant-soil system involving the carbon cycle. Granger causality is a statistical approach, originating in econometrics, used to identify the presence of linear causal interactions between time series of data, based on prediction theory. We first test to see if there was a significant difference in the causal structure among two treatments where carbon allocation to roots was interrupted by girdling. We then expanded the analysis, introducing radiation and soil moisture. The results showed a complex pattern of multilevel interactions, with some of these interactions depending upon the number of variables in the OPEN ACCESS Entropy 2013, 15 4267 system. However, no significant differences emerged in the causal structure of above and below ground carbon cycle among the two treatments.


Introduction
The seasonal rate of change in atmospheric carbon dioxide (CO 2 ) above temperate regions is largely driven by the net sum of two opposing processes: carbon (C) uptake through Gross Primary Production (GPP) and carbon release from the vegetation and soil through plant and microbial respiration.The magnitude of the annual ecosystem-driven change to the global soil C pool is such that even a small increase in soil respiration (R s ) rates could result in a net reduction in net terrestrial C sink strength and an increase in atmospheric CO 2 concentration [1][2][3].
The soil ecosystem comprises a multitude of organisms acting over a broad range of spatial and temporal scales, from heterotrophic soil fauna to plant roots, all embedded in a complex and highly variable physical environment [4].This variation is both spatial, through differences in soil microbial communities and root structural characteristics, and temporal, through fluctuations in energy and water dynamics and through disturbances.
A principal challenge to understanding biogeochemical cycles is identifying the causal agents of process rates, which requires separation of endogenous dynamics from the effects of the timedependent, and often correlated, forcing variables.Any coupling among ecological variables induced by external drivers is difficult to distinguish from the endogenous correlation structures generated by the core ecological system.Major complications arise from the inter-correlation of external forcing variables, indirect effects of one variable over another, and feedbacks.For example, air and soil temperatures are strongly correlated to incoming shortwave radiation.Part of this energy reaching the surface, stems and leaves, is transformed into sensible heat, which warms the air above, and into soil heat flux, which warms the soil below.Sensible and latent heat fluxes are correlated, as both are constrained by available energy (net radiation minus soil heat flux) [5].Soil moisture can limit plant metabolic rates and microbial activity and therefore controls, with other drivers, the reaction rates of many biogeochemical processes in the soil.Hence, it is not surprising that in ecological time series analysis, one of the main challenges is to determine how one state (or rate) variable relates to or "causes" the magnitude or dynamics of another.The problem is naturally compounded when both the affecting and affected variables are subjected to multiple periodic and external forcings, and when state changes lag behind causal agents.The conventional approach to studying the relationships between affected variables and the forcing that drives them employs correlation [6].However, the existence of correlation itself does not necessarily demonstrate causation.First, due to the symmetric nature of correlation it cannot distinguish between the causal roles of the affected and affecting variable.Second, correlation cannot distinguish whether both variables are directly interacting or if the appearance of correlation is a result of both variables being forced by a common third variable.
Identification of causality, based on the argument that the response cannot precede the perturbation, has been at the core of a number of disciplines of prediction [7].Among the more popular measures of causality are the so-called Granger causality (hereafter referred to as G-causality)-the main subject of this study-and transfer entropy.For Gaussian variables it has been demonstrated that these two approaches are equivalent [8], and this result has been extended to other common types of distributions [9].
The concept of G-causality originated in econometrics [10,11] but is now being used in a wide range of different fields, among them the analysis of climatic and environmental time series [12][13][14][15].Compared to several information-theory driven methods that are used in ecohydrology [16,17] the main advantage of G-causality is that it has an explicit spectral and non-parametric formulation [18,19].Ecological processes, especially those connected with the exchanges of mass and energy between biosphere and atmosphere, often display strong periodicities due to diurnal and seasonal cycles [20,21], and other stochastic oscillations due to large climatic events, such as front passages and atmospheric pressure fluctuations.Hence, it is convenient to have a measure of causality in the frequency domain.Furthermore, when the autoregressive order of the data and time lags among variables are not known a priori, non-parametric approaches eliminate the need of autoregressive modeling [19].When applying G-causality analysis to ecological processes such as soil biological active emissions Detto et al. [15] and Hatala et al. [22] showed that this spectral approach presents a decisive advantage when compared with the time based scheme.With this method, quasi-periodic forcing can be identified, localized within a specific and limited number of frequencies (e.g.seasonal, diurnal) and removed via conditional analysis.
Extensive logging and fire in the late 19th century through the early 20th century reduced most of the forest in the region to an early successional forest [23].As the early successional trees of this forest approach maturity, ecological succession is driving an intermediate disturbance process in this region, characterized by increased mortality of aspen (Populus) and birch (Betula) trees.To hasten this process of mortality, we have manipulated a 33 ha plot by stem girdling all early successional (aspen and birch) trees in the Forest Accelerated Succession Experiment (FASET) [24].A flux tower was installed in the treatment plot, and an existing flux tower in a nearby untreated plot serves as a control.
Several dynamic chambers were deployed in the FASET and control sites for continuous measurements of CO 2 effluxes.Understanding the underlying variability of R s responses to disturbance and canopy structural change is confounded by the covariance at multiple frequencies of the biophysical factors hypothesized to be the most important drivers of R s -soil temperature, soil moisture, radiation, and GPP [25,26].
We used this long-term, large-scale experiment to provide a test case for the application of Granger-causality analysis in determining the driving variables of R s dynamics, and particularly its outcome-R s in a forest of the Great Lakes region.To understand the motivations of using this particular statistical approach, we next frame the problem within the relevant ecological context.We describe a multivariate approach for testing the causal structure of relationships between R s and its meteorological and ecological drivers, and whether they differ between an undisturbed site and a site following prescribed disturbance in the form of stem girdling.We apply spectral G-causality analysis in a layered multi-variant approach for testing of the significance of causal relationships between one variable and another, given all other co-occurring environmental conditions, thereby choosing which variables, among an available set, are important predictors.

Site Description
Our study site is located at the University of Michigan Biological Station in northern Michigan (45°35.5'N 84°43' W).The forest is dominated by maturing early successional tree species, relatively even-aged Populus grandidentata Michx.(bigtooth aspen), Populus tremuloides Michx.(trembling aspen), and Betula papyrifera Marsh.(paper birch).Other canopy tree species include Fagus grandifolia Ehrh.(American beech), Acer saccharum Marsh.(sugar maple), Acer rubrum L. (red maple), Pinus strobus L. (white pine) and Quercus rubra (red oak).The early successional aspen and birch species are beginning to senesce and will continue to do so over the next 50 years [27].Total basal area density is 30.83m 2 /ha in the treatment plot and 36.87 m 2 /ha in the control with aspen and birch composing 40.0% of the total basal area in the treatment plot and 53.1% of the control plot.On 20 April-3 May 2008 we initiated an experiment in which we stem girdled all aspen and birch trees in a 33 ha treatment plot, totaling approximately 6700 trees.Girdling removes the bark and underlying phloem in a strip around a tree, preventing the translocation of photosynthate from leaves to roots.An unmanipulated control plot with similar site and soil characteristics is located nearby and within the footprint of an eddy-flux tower (the US-UMB Ameriflux tower, [24]).The mean annual temperature is 5.5 °C (1942-2003).The mean annual precipitation is 817 mm (including 294 cm snowfall).The site is snow covered in most years from November to April.The study area encompasses 140 ha on a high outwash plain and adjacent gently sloping moraine.Nearly all soils in the study area are well to excessively well drained Haplorthods of the Rubicon, Blue Lake, or Cheboygan series.About 53% of the fine-root mass is located within the upper 20 cm of the soil profile.Forest floor C mass is 5-15 Mg C/ha, and the mineral soil is 95% sand and 5% silt, with pH 4.5-5.5 in water.Soils are of low fertility, with total N capital to 40 cm depth of 2000 kg/ha, an average in situ net N-mineralization rate of 42 kg N/(ha yr) and < 2% net nitrification [24,28].

Soil CO 2 Efflux
Automated chambers were used to continuously measure soil CO 2 efflux in the control and treatment plots.Air was circulated from the chamber to an infrared gas analyzer (IRGA, LI-6252, Li-Cor, Lincoln, NE, USA) and then returned to the chamber.At the treatment and control sites, eight closed system automated chambers were randomly distributed on the forest floor near the base of each flux tower.Each chamber consisted of an aluminum frame with a pneumatic actuator to control the opening and closing of the lid and a clear plexiglass cylinder (radius 14.75 cm) with a small fan near the top of the chamber to homogenize the chamber air as samples were taken (Figure 1).Within each chamber a manifold returned sampled air to the bottom of the chamber.Air from the chamber headspace was sampled through a polyethylene coated aluminum tube and pumped to the IRGA for analysis.The system was controlled by a data logger (CR23X, Campbell Scientific Inc., Logan, UT, USA) and each chamber was sampled at 10 Hz for seven minutes each hour, so that each of the eight chambers in each array was sampled once every hour.Measurements started at both sites in 18 April 2008.

Figure 1.
(A) Automated chambers closed when pressure was applied to the pneumatic actuator which created a seal between the lid and the chamber.(B) Samples were collected from the black tube near the top center of the chamber, and analyzed air was returned via the perforated manifold at the base of the chamber.(C) Respiration was determined by analyzing the timeseries of concentration accumulation in the chambers-raw CO 2 concentration data (blue line) were graphed for 330 s while CO 2 accumulated in the closed automated chambers.Each accumulation curve was fit individually with a non-linear curve (black line) in MATLAB [29] (see methods).Coefficients were computed from the curve and exported to yield R s .
The program fit a non-linear curve (Equation ( 1)) to those data points which yielded coefficients from which R s was calculated: (1) When C c (t) is the concentration of CO 2 in ppm at time t (s), C s is the soil CO 2 concentration (ppm), C 0 is the CO 2 concentration of the air sample (ppm) at time t = 0, and k is the concentration saturation rate (1/s).The initial rate of change of carbon concentration in the chamber immediately after closing the lid can be solved as: As this rate of CO 2 accumulation is true at the instant of closing the lid it is also characteristic to the rate of emission from the soil immediately before closing the lid, an thus used to calculate the diffusional flux of CO 2 from the soil (D c ), which is equivalent to soil respiration, R s (µmol/m 2 /s): (3 where S is the chamber surface area (0.068 m 2 ), V is the chamber volume (0.0144 m 3 ) and ρ is the molar density of air (41.6 mol/m 3 ).To solve for C s and k an estimate of C 0 is needed.To obtain C 0 a minimizing function was employed to find the lowest 20 s moving average within the first 70 data points.Using C 0 , t, and C c data, a non-linear least squares fit with least absolute residuals returned estimates of C s and k. Figure 1C shows examples of raw CO 2 data and the fitted curve for one hourly 7 min period of data from which coefficients were calculated.C 0 , C s , and k values for each seven minute period for each site, day, hour, and year and chamber combination were calculated.
Periods for which the estimated C 0 value was less than 380 ppm, the current ambient air CO 2 concentration, were removed.Estimates where the R 2 value for the fitted line was below 0.9 were removed.Negative D c values and values above 15 µmol/m 2 /s were removed as outliers.Fluxes above 15 µmol/m 2 /s were determined to be outliers because they did not fall within the range of values reasonable for our site [30]; outliers made up less than 1% of the data.The complete dataset is illustrated in Figure 2

Other Measurements
The following micrometeorological data were collected in both treatment and control plots: air temperature (T a ) and humidity from an aspirated and shielded temperature/humidity probe (model Hygroclip S3, Rotronic, Bassersdorf, Switzerland); above canopy total incoming and diffuse photosynthetic active radiation (PAR) using a pyranometer (model BF2, Delta-T Services, Cambridge, UK); Net radiation (RAD n ), short and long wave downwelling and upwelling radiation streams using a 4-channel radiometer (model CNR1, Kipp & Zonen, Delft, The Netherlands); precipitation above the canopy using a tipping bucket (model TR-525-M, Texas Electronics, Dallas, TX, USA).Sensors were deployed at 46 m above ground on a tower at the control plot and at 32 m on a tower at the treatment plot.Soil temperature (T s ), was measured at depths of 2, 7.5, 20, 50 and 100 cm below ground at each tower site (Type E Thermocouples, Campbell Scientific, Inc.); Soil water content (SWC) was measured at two points near each tower site and at two depths (7.5, 20 cm) at each point using a soil-moisture sensor (models CS615-L and CS616-L, Campbell Scientific, Inc.).For each site, soil moisture and temperature series were computed taking the spatial average along all depths and location.Above-canopy and soil observations were logged every minute and averaged to half-hour periods.
The net flux of CO 2 between atmosphere and forest (net ecosystem exchange, NEE), the surface sensible (H) and latent (LE) heat fluxes were measured using the flux-covariance approach.Wind velocity and temperature fluctuations were sampled at 10 Hz using three-dimensional sonic anemometers (model CSAT3, Campbell Scientific, Inc.); CO 2 and water vapor concentrations were sampled at 10 Hz using closed-path infrared gas analyzers (IRGAs) (model LI-7000, Li-Coer Bioscience).Measurement outliers and measurements marked as faulty by the sensors' quality control variables were removed (despiked).Wind data were processed using a 3-D coordinate rotation assuming half-hourly mean vertical wind velocity is 0 and rotating the horizontal wind components toward the mean wind direction [31].Lag time between the anemometer and scalar concentrations was calculated using the maximal covariance method [32].The Schotanus correction [33,34] and conversion to "real" temperature [35] were applied to sonic anemometer measurements of temperature.Wind and concentrations observations were compiled into half-hour block averages of net fluxes following the AmeriFlux recommendations [36].Water vapor and CO2 concentrations were adjusted using the Webb, Pearman, and Leuning correction [37] in a modified form derived by Detto and Katul [38] as a correction for the 10 Hz time series of the scalar.Observations were divided into seasons according to the carbon flux phenology in the site [39] and only observations during the growing season were used here.Data were filtered based on a seasonal frictional velocity (u*) threshold criterion [40,41], with a prescribed maximum u* threshold of 0.35 m/s [42], typical growing season u* threshold values were between 0.29 and 0.35.We used a bilinear periodic method to fill gaps in temperature, moisture, humidity, and radiation observations and assumed that CO 2 fluxes during nighttime were driven entirely by ecosystem respiration (R e ).R e was calculated using site-specific empirical formulas developed for our site [30], which relate nighttime NEE to soil moisture and temperature [43].We used this R e equation to gap-fill NEE during the nighttime.Gaps in GPP were filled using the mean of 100 neural network simulations [44].Gap-filled GPP was added to R e to provide gap-filled NEE.Since the FASET-plot tower sampled a footprint that was sometimes larger than the 33 ha of contiguous experimental girdling area, we used a 2-D footprint model [45] that we modified to automatically integrate the flux-source probability over the treatment area and thus provide an index for the treatment footprint probability in each 30 min block average period.We then used the probabilistic flux footprint climatology [46] to scale our conclusions to fluxes originating only from the FASET plot.More details about the flux data analysis in our site are listed in [24,39].Time series are shown in Figure 3.

Granger Causality
Consider two discrete time random variables X and Y with autoregressive representation in the form: where  and  are white noise prediction errors with covariance matrix: a, b, c and d are coefficients describing the linear interactions between the variables, with subscript j indicating time lags.When the above equation is compared with a univariate model , and when the multivariate model outperforms the univariate case, (e.g., Y is said to have a causal effect on X (and similarly for the effect of X on Y).This is the statistical interpretation of causality proposed by Granger [11] and is commonly referred to as Granger or G-causality.G-causality is a measure of coupling with explicit time directionality.For this reason, it is based on prediction errors rather than on linear interactions among coefficients.Traditionally, it is expressed as the logarithmic ratio between the residual variance of the bivariate and the null model (i.e., the univariate case).
Let's consider a multivariate system of k + 2 stochastic variables (X, Y, Z 1 , …, Z k ) which admits spectral representation and factorization in the form [47]: S and U are spectral matrices respectively, of the complete system and of a system which does not include the variable Y, whose causality is tested, * indicates matrix adjoint.
The conditional G-causality of Y on X given Z 1 , Z 2 , …, Z k is computed as [48]: where: are the corrected transfer function matrices that separate the pure directional interactions [18].The rotation matrices P are normalization matrices needed to recast the multivariate systems in the canonical form (with uncorrelated errors).
If the variables X and Y are not interacting directly, the numerator and denominator of Equation ( 7) are equal and ( )  0

Estimation of G-Causality
There are two approaches to estimate G-causality: parametric and non-parametric.The former fits the parameters of the autoregressive models and computes the spectral matrices.This is usually done with least-square fit (e.g., Yule-Walker estimator) upon prescribing the order of the model a priori.For short series exhibiting apparent periodicities, superimposed on long-range memory, this approach may be problematic, because higher order autoregressive models are required to reconstruct the observed dynamics.The estimation error of an autoregressive model is known to increase with order and order misspecifications may produce spurious causality effects [49].Moreover, for the multivariate case, further problems arise when the spectral estimates of different AR models (Equation ( 14)) are not identical due to practical estimation errors.This requires implementing a partition matrix improvement [48].
Alternatively, the spectral matrix S() can be estimated directly from the data using classic routines (multitaper, wavelet transform or analogous spectral methods).The factorization theorem of a spectral matrix [50] ensures that, under fairly general conditions, any spectral matrix, can be decomposed, or factorized, into S() = H()  H * ().Iterative algorithms have been proposed to evaluate H and  from S. Here, we used the algorithm formulated by Wilson [51].Since S, H and  can be inferred from the time series, the method does not require any assumption regarding the autoregressive order of the model representing the data and the estimation of G() is considered non-parametric [19].
The robustness of this statistic must be tested against the null hypothesis that the two series do not exhibit any causal interactions at a particular frequency.In general, the analytical solution for the causality structure of the system is not known a priori, necessitating the use of Monte Carlo simulation methods and surrogate data [52].In our analysis, the data were divided into blocks of 60 days and spectral matrices were calculated as an ensemble.A significant interaction is found when both bivariate and conditional G-causalities are greater than the 99th percentile of the results from the randomly permutated blocks.

Direct Connectivity Diagram
Cause-and-effect diagrams can be built based on the outcome of the multivariate G-causality analysis.Each variable is represented as a node of a circular network and significant causal interactions are depicted with arrows pointing in the direction of the influenced variable and its width proportional to the intensity of the interaction.To construct such a diagram we follow two rules: (1) Only direct causal interactions are depicted.This implies that if the bivariate interaction is significant, but the multivariate is not, an arrow is not drawn.However, if the bivariate G-causality is not significant, no further action is taken.(2) If both directional interactions are detected among two nodes, only the greater is retained.
The spectral G-causality allows creation of connectivity diagrams at different frequencies or averaged in specific bands.

Results
Automated soil respiration measurements in the control plot were consistent with those of long-term point measurements [49].Findings from the present study demonstrate that the treatment site had lower R s than the control site for all years (Figure 3).
Bivariate and conditional G-causalities on GPP and T s over R s showed a large effect of temperature and a much smaller effect of GPP, especially at low frequencies (Figure 4).At daily frequency, and in minor extent at the corresponding half day sub-harmonic, the effect of GPP on R s appears to be important, although it is largely reduced when T s is included in the analysis.There are no apparent differences between the two sites in response to temperature and/or GPP.This hypothesis was tested computing the difference . Confidence intervals obtained upon bootstrapping showed that the differences were not significant (Figure 5).
The analyses were more complicated when other variables were progressively included.To synthetize the results we aggregated the G-causalities in either long time scales (0.2 < frequency < 0.7625 day −1 ) and short time scales (0.7625 < frequency < 1.7 day −1 ).Given the fact that the results in both sites were very similar, we also aggregated the data without differentiating among sites.The connectivity graphs (Figure 5) depict only significant direct interactions, i.e., not mediated by other variables.We consider three cases.In the first case only GPP, R s and T s are included.This diagram supports the hypothesis of GPP control over R s at long time scales, though at short time scales a significant interaction cannot be detected (Figure 6A).Note the biologically improbable interaction of GPP over T s especially at short time scales.
In the second case RAD n also is included in the analysis (Figure 6B).Introducing this variable changes the causality pattern slightly.The causal effect of GPP over R s can still be detected at long scales, although considerably reduced.The causality of GPP over T s , however, disappears and RAD n influences all variables at all time-scales, except for GPP at short scale.This makes sense, as RAD n controls all energy processes and is closely coupled with biological processes such as GPP.The disappearance of the influence of GPP on T s , which was previously detected, indicated that GPP influences T s indirectly through an unmeasured variable not included in the 3-variable analysis, and most probably RAD n .The lack of G-causality of RAD n over GPP at short scales can be explained by the fact that radiation influences GPP almost instantaneously, while G-causality only detects lagged responses at the time-scale of the measurements, here-one hour.
In the third case, soil water content was introduced (SWC) (Figure 6C).At short time scale the scheme of interactions does not change as SWC does not act at daily or higher frequencies.However, at long time scales, a complex pattern of casual interactions emerges.RAD n and SWC largely control all other variables and the influence of GPP over R s is lost.Interestingly, SWC is not influenced by any of the other variables, i.e., SWC is Granger autonomous while RAD n is influenced by SWC.

Discussion
Soil respiration rates and the contribution of this flux to R s varies considerably among forest types, though the fundamental variables regulating carbon cycling process rates are generally conserved across ecosystems [53].In the mixed deciduous forests of Michigan, USA, respiration from soils is the largest component of ecosystem respiration, comprising approximately 70% of the total flux of CO 2 [30].Ecological transitions and regional climatic changes can affect Rs.Carbon sequestration by an ecosystem is determined as a small difference between two large fluxes: respiration, which emits CO2 to the atmosphere and gross primary productivity (GPP) that takes up CO 2 from the atmosphere through the process of photosynthesis.The small difference between these two opposing fluxes makes the net CO 2 sequestration rate highly sensitive to any change in either respiration of GPP.
These climatic changes have motivated the scientific community to quantify R s and to better understand the factors controlling it, particularly during and after disturbance [54].Primary climatic controls on R s are temperature, solar radiation, and precipitation [55].Soil respiration increases as soil temperatures rise, because microbial community abundance and activity both respond positively to warming [27,56].Drought and extremely wet conditions limit microbial activity and decrease respiration [57].Over periods of hours to days, photosynthetic activity also controls carbohydrate supply from leaves to roots resulting in increased C availability to mycorrhizal fungi and microorganisms in the root environment (rhizosphere), thereby increasing R s [58][59][60][61].
Forest disturbance can change the rates of photosynthetic activity through a reduction in live leaf area and thus GPP, and resulting in changes the soil's energy budget as temperature and humidity shift in response to reduced shading and water uptake by trees.Disturbance also changes canopy structure and the sub-canopy environment in ways that may affect R s , modifying light and nutrient use efficiency [62], roughness length, canopy-atmosphere coupling, and increasing turbulence mixing above the soil [63,64].In both 2008 and 2009 respiration was lower in the treatment site than in the control.Due to the lack of respiration measurements prior to the girdling treatment we cannot rule out that this difference is a result of the treatment.
However, the difference in respiration between the sites was rather constant over the period from 3 months to a year and a half past the treatment, despite stem girdling and eventual and slow mortality (trees where fully alive in summer 2008, and mostly dead in 2009) of canopy dominant trees that occupied roughly 30% of the canopy in 2008.Additionally, the dynamics and causality structure of the drivers of respiration in both sites was not different.Given these two indications, it seems unlikely that that difference is a result of the treatment.This was not because of an ineffective treatment.In summer 2009 the aspen and birch trees were visibly defoliated.However our study indicates resilience of ecosystem functioning to moderate disturbance, and particularly the soil ecosystem which maintained high respiration rates despite a reduction of photosynthates from girdled trees.In part, this resilience is due to compensatory growth and increased productivity of the remaining un-girdled trees from increased light and nitrogen availability [65].Alternatively, soil respiration may have declined in the treatment very abruptly following the 2008 girdling, shortly before our chamber measurement started.Several studies document an almost immediate decline in R s that reaches its lowest point within days following stem girdling [54,66], suggesting that the exhaustion of stored labile carbohydrates progressively limits root and rhizosphere respiration.
Our G-causality analysis suggested that GPP did not directly control R s in either treatment or control site, and that RAD n was a primary environmental constraint on R s through its effects on T s at short time scales.At longer scale, both SWC and RAD n controlled R s .This finding is similar to other studies suggesting no diurnal relationship between GPP and R s [67] and in contrast to others showing that GPP or photosynthesis constrains R s over diurnal timescales [15,58,59,68].Respiration from our sites is principally limited by soil moisture and radiation (heat to the soil), suggesting that RAD n , which is remotely sensed with high accuracy, may serve as a useful predictor of R s .In another site, the Duke Forest, where a significant link between GPP and R s was established [15], moisture and temperature were less limiting factors.
Furthermore, under natural conditions, there is a substantial attenuation of the diurnal amplitude of photosynthesis along the phloem path length, primarily influenced by vegetation height and phloem conductivity [69].At this site, trees were quite tall (around 25 m).Unfortunately, information about phloem conductivity was not available.
Several studies have shown that recently produced photosynthate is a key driver of root respiration and thus, girdling will typically decrease soil respiration [54,66].Species-specific responses to girdling can be highly variable [70] and in some cases, no effect of girdling on soil respiration was reported [71].Studies of respiration in pest-infected forests, where the trees were girdled by insects, found mixed results from no effect whatsoever, to a large decline [72,73].
Other methods for testing GPP and R s interactions were reviewed [60].Among non-manipulative approaches, methods based on time series analyses, such as cross-correlation, pulse-response and Fourier decomposition, were suggested to detect lagged response.However, we note that only G-causality is able to provide a directionality of the interactions, which is the real signature of a causal mechanism.Furthermore, the multivariate approach investigated in this study showed that including more biophysical factors may reveal indirect or apparent interactions that, in lower dimension systems, may appear as direct and significant.However, the proposed method has also some limitations.Long, continuous and multiple time series are required to construct robust spectral matrices and reduce noise.Although G-causality is less sensitive to random Gaussian noise, other non-Gaussian source of error (e.g., spikes) may cause unpredictable effects.Any additional variable that contains such errors may reduce the statistical power to the extent that weak interactions may not be detectable.

Figure 2 .
Figure 2. Daily soil respiratory carbon loss (R s ) calculated from the 1 hour soil respiration measurements (points) between Control (blue), and Treatment (red) plots in 2008, 2009, and 2010 growing seasons.

Figure 3 .
Figure 3.Time series of the four variables measured at control (AmeriFlux in blue) and treatment (FASET in green) sites.

Figure 4 .
Figure 4. Bivariate (dashed lines) and conditional (thick lines) G-causalities of GPP and T s over R s .Top panels represent control (AmeriFlux) and bottom panels, treatment (FASET).The thick lines represent the desired goal of the analysis because they isolate direct effects, while the dashed lines indicate the apparent connections which are due to a combination of direct and indirect effects.

Figure 5 .
Figure 5. Difference of conditional G-causalities of GPP given T s over R s between control (AmeriFlux) and treatment (FASET).Dashed lines are the 95th quantiles obtained boostrapping the 60 days blocks of the pulled dataset and divide them in two random groups.

Figure 6 .
Figure 6.Connectivity diagram based on the significant direct G-causalities using (A) three, (B) four or (C) five multivariate systems.The widths and values on the arrows indicate average G-causalities for long time scale (0.2 < frequency < 0.7625 day −1 ) and short time scales (0.7625 < frequency < 1.7 day −1 ).