Detecting the Differences in Responses of Stomatal Conductance to Moisture Stresses between Deciduous Shrubs and Artemisia Subshrubs

Shrubs and subshrubs can tolerate wider ranges of moisture stresses in both soil and air than other plant life forms, and thus represent greater nonlinearity and uncertainty in ecosystem physiology. The objectives of this paper are to model shrub/subshrub stomatal conductance by synthesizing the field leaf gas exchanges data of 24 species in China, in order to detect the differences between deciduous shrubs and Artemisia subshrubs in their responses of stomatal conductance to changes in the moisture stresses. We revised a model of stomatal conductance by incorporating the tradeoff between xylem hydraulic efficiency and cavitation loss risk. We then fit the model at the three hierarchical levels: global (pooling all data as a single group), three functional groups (deciduous non-legume shrubs, deciduous legume shrubs, and subshrubs in Artemisia genus), and individual observations (species × sites). Bayesian inference with Markov Chain Monte Carlo method was applied to obtain the model parameters at the three levels. We found that the model at the level of functional groups is a significant improvement over that at the global level, indicating the significant differences in the stomatal behavior among the three functional groups. The differences in tolerance and sensitivities to changes in moisture stresses are the most evident between the shrubs and the subshrubs: The two shrub groups can tolerate much higher soil water stress than the subshrubs. The analysis at the observation level is also a significant improvement over that at the functional group level, indicating great variations within each group. Our analysis offered a clue for the equivocal issue of shrub encroachment into grasslands: While the invasion by the shrubs may be irreversible, the dominance of subshrubs, due to their lower resistance and tolerance to moisture stresses, may be put down by appropriate grassland management.


Introduction
Dynamic global vegetation models (DGVMs) and regional ecosystem models have been making projections of structures and functions in response to changes in climate and anthropogenic activities [1,2]. However, there existed great disagreements among the models, largely due to the uncertainties in physiological parameters that control the carbon assimilation into and emission out of ecosystems [3]. Ecosystem model parameters are often derived from plant functional traits, which have been recognized as important links to ecosystem functions, structures, and adaptations [4,5]. Relationships among maximum leaf photosynthesis, maximum stomatal conductance, specific leaf area, leaf life span, leaf size, and leaf nitrogen have been found from the global leaf traits dataset by means of meta-analyses [6,7,8,9]. Shrubland ecosystem physiology in the sub-humid, semiarid, and arid regions involves great uncertainty because of the superior tolerance to moisture stresses, morphological plasticity, and adaptability [10].
Shrubs have been reported to invade grasslands in many places over the world [10,11,12,13,14,15]. Many causal factors have been hypothesized to trigger the processes of shrub encroachment. Among these hypotheses, increased drought frequency and shifted rainfall seasonality/intensity have been considered as major drivers [16,17]. These hypotheses are based on a common assumption that shrubs are more tolerant to drought than grasses. The parameterized stomatal model used in the Patch Arid Land Simulator (PALS) [18,19] showed that the shrub stomata can tolerate much more severe soil water stress than the grass stomata. When soil moisture is ample, the grasses showed greater stomata conductance than the shrubs.
Shrubs tend to have deep roots, high ratio of leaf to sapwood areas, low vertical shading, and strong morphological plasticity to adapt to variations in climate and soils [20]. The same shrub species can be phreatophyte in sandy soils, but xerophyte in heavy clayey soils to maintain active when soil water potential lower than 25 MPa [21,22].
Models of stomatal conductance are important tools to quantify leaf production processes of shrubs, since leaf photosynthesis is colimited by stomatal conductance and biochemical carboxylation [23]. Many stomatal models exist in the literature, with different mechanistic/empirical assumptions and treatments. Ball et al. [24] proposed a simple empirical model of stomatal conductance as a function of net photosynthesis rates and relative humidity on leaf surface. Another model was developed to count the composite effects of net photosynthesis, vapor pressure deficit, and CO 2 pressure on stomatal conductance [25]. More mechanistic model of stomatal conductance that considered the transient response of stomatal conductance to changes in driving variables was developed by Buckley et al. [26,27] Incorporating stomatal conductance with hydrological structure and function of tree canopy showed complex interactions between crown structure and hydrological function [28].
A semi-mechanistic model was developed by Gao et al. [29] to calculate stomatal conductance as a function of soil water potential, vapor pressure deficit in air, intercellular CO 2 concentration, and light irradiance. The model was revised [30]to consider the cavitation loss of xylem conductance and hydrological capacitance, the model considered the cavitation loss of xylem hydraulic conductivity linearly dependent on xylem water potential, so that where K soil{to{leaf is the apparent soil-to-leaf conductance, g p is the maximum conductance, and f is the parameter signifying the linear dependence of K soil{to{leaf on y x , the xylem water potential. However, recent prevailing literatures on functional anatomy of plant xylem of various life forms revealed that the loss of xylem hydraulic conductivity due to cavitation and embolism is significantly related to the vessel diameters [31,32,33,34,35,36,37,38], despite of some equivocal results with insignificant relationship from a few studies [39,40]. According to the prevailing literatures, vessels with larger diameter provide greater hydraulic conductivity, however, larger vessels are more prone to cavitation damage by negative pressure inside. Therefore, there exists a tradeoff between the hydraulic conductivity and its sensitivity to cavitation loss. The results of these experiments also indicated that the whole-plant (soil-to-leaf) conductivity decays exponentially with the xylem pressure, rather than the linear assumption in the previous stomatal models.
In this paper, we revised the model of stomatal conductance by Gao et al. [29,30] according to the prevailing literatures to reflect the tradeoff between the hydraulic efficiency and cavitation resistance. We then fit the revised model using the field diurnal measurements of leaf gas exchange of 24 shrub species in China to test the hypothesis that shrub functional groups differ significantly in their stomatal responses to changes in source (soil) and air moisture stresses. The results indicated that the three shrub functional groups (deciduous non-legume, deciduous legume, and Artemisia subshrubs) were significantly different in stomatal behavior. The sensitivities of stomatal conductance of the three shrub functional groups to relevant driving variables were calculated and discussed in the context of experimental evidences.

The Model
We revised the model of stomatal conductance by Gao et al. [30] so that the dependence of soil-to-leaf hydraulic conductance (K soil-to-leaf ) on xylem water potential (y x ) is hyperbolic rather than linear, and that the loss of xylem hydraulic conductance depends on the maximum soil-to-leaf conductivity (g p ). Specifically, where l is a parameter. Moreover, we assumed l is proportional to g p , i.e.l~C l g p , so that where C l is a constant across all species and all levels of the analysis. This formulation allowed us to balance the hydraulic conductance and the xylem safety, hence the greater the gp, the more vulnerable the xylem system. In addition to the above changes, the osmotic adjustment is now dependent on net photosynthesis A n rather than light intensity I p [41], so that the osmotic pressure of plant leaf and guard cells, p, is where p 0 is the baseline osmotic pressure (MPa), C i and C 0 are the intercellular and reference CO 2 concentrations, respectively, and j and K i are parameters. The model now takes similar treatment as Ball et al. [24,25]. The assumptions in the previous model [30] include: stomatal conductance, g s , is proportional to the leaf turgor pressure P turgor , where K y is the apparent compliance of guard cell structure; and transpiration is the product of stomatal conductance and the scaled vapor pressure deficit (D), i.e. the absolute vapor pressure deficit (VPD) divided by air pressure (P). The mass balance between leaf transpiration and soil-to-leaf flow gives the equation where y is the soil water potential. Solving the Equation (5) for y x and substituting the result in (4) give the following model where Equations (3, 6-9) constitute the complete model with 3 basic parameters of K y ,g p , and C l , and three osmotic adjustment parameters of p 0 , K i , and j.

Data Collection and Preparation
This analysis used the field diurnal gas exchange data of shrub leaves in northern China ( Table 1). Part of the data were collected by the authors and colleagues during 2002-2007, and others were reported in the literature from 1990 to 2003. All the field measurements were done in public lands with free access to researchers. No specific permissions were required since the experiments had no obvious impacts on the sites and did not involve any endangered or protected species. The data from literature were read from tables and charts. Measurement on one species at a site is called an observation, and an observation may be done in multiple days. Each observation produced a data table with a number of variables (columns) including stomatal conductance/resistance, net photosynthesis, air pressure, and vapor pressure deficit, at leaf surface. Each record (row) in the data table represented one replicate of measurements on one plant leaf surface at specific hour of the day. There were a total of 43 observations within 80 diurnal measurements involving 24 species in this analysis. Table 1 also listed the approximate geographic location, elevation, instrument used, and data source. Data measured by the authors and colleagues are provided in Table S1.
To facilitate the analysis, all datasets of the diurnal measurements were converted to the formats and units used by the portable photosynthesis system Licor-6400 (Licor, Nebraska, USA). To avoid the problem of repeated measurements on the same plant leaf, the records of all data tables were averaged for each leaf at a particular hour.

Procedures of Analysis
The observations were grouped into three functional groups: non-legume deciduous shrubs (DCDS), deciduous legume shrubs (LEGM, mostly in Caragana genus), and subshrubs (SUBS, Artemisia genus). The LEGM group, with their potential nitrogen fixation capability, may be advantageous over the DCDS when nitrogen is limiting. The subshrubs have characteristics of both woody and herbaceous plants (woody based but herbaceous primary growth).
The above model of stomatal conductance was fitted to the data at three hierarchical levels: observations (OBS, species cross sites), functional groups (PFT), and global (GLB, pooling all data into a single group). Parameter C l is treated as a super parameter for all the species and observations, so that the tradeoff between conductivity and xylem safety is realized. Model parameters were estimated by means of Bayesian Inference using Gibbs sampling with Markov Chain Monte Carlo (MCMC), implemented in the WinBUGS program [42,43] run within the free software of R [44].
Fitting the stomatal conductance model required the data for soil water potential which is not available. To get around the problem, we made an additional assumption that soil water potential is approximately constant for each day and the daily soil water potentials were estimated at the observation level. The assumption behind this treatment is that minimizing the difference between the model and the data should also allow us to estimate both model parameters and the unknown daily soil water potential [29,30]. The estimated daily soil water potentials from the observation level were used at the levels of GLB and PFT.
MCMC requires specification of prior probability distributions for the parameters. On the other hand, it is generally difficult for the iterations to converge to the physiologically meaningful results because of the nature of nonlinear models. We handled this difficulty by providing priors of uniform distributions within the specified physiologically meaningful ranges. This allowed us to place the following bounds to the parameters: K y [ 0:01,1: The parameter C l was only fitted at the observation level, so that the parameter at the other two levels (PFT and GLB) was fixed at the mean value from the observational fit. Since the previous studies have recognized the higher drought tolerance of shrubs than other functional types, we placed the bounds of daily soil water potential as y[ {0:03,{3:5 ½ MPa. These bounds were realized in the specification of uniform prior distribution. For example, a statement in the WINBUGS program, ,dunif (0.01, 1.5) specifies that the prior for the parameter K y is a uniform distribution between 0.01 and 1.5 mol m 22 s 21 MPa 21 . For each level of the hierarchical analysis, the MCMC was iterated for 5,000 times so that almost all parameters can converge.

The Statistics of the Fitted Models and Parameters
Summary statistics (Table 2) for the analysis at the three levels indicate that from global to shrub functional groups to individual observations, the deviances of the models decrease, while the number of parameters increases. The deviance for stomatal model is negative because most measured and predicted stomatal conductance are smaller than 1. Specifically, the deviance decreased from 21,423 to 21,804 to 23,981, and the standard error of the residual decreased from 0.170 to 0.158 to 0.095 mol m 22 s 21 , for GLB, PFT, and OBS levels, respectively. Smaller deviances indicate better fits of the models to the data, in the price of increased model complexity with more parameters. The difference in deviances between any two hierarchical levels has been shown to follow the x 2 distribution with the degrees of freedom equal to the difference in number of parameters between the two levels. Significant differences in model quality between the two hierarchical levels should be indicated in a significantly large x 2 value. The deviance tests (Table 2) among the three levels, viz. PFT vs. GLB, and OBS vs. PFT, indicate that the model for shrub functional groups (PFT) is a significant improvement over the global model (GLB), and that there is significant difference in behavior of stomatal conductance among the three shrub functional groups. Similarly, the model at the observation level is significantly better than the models at the level of functional groups. The trend of improved model fit from global to functional groups to observations is further demonstrated in the increased correlations between the model predictions and measurements ( Table 2). Specifically the Pearson's correlation coefficient between the measured and predicted stomatal conductance increased from 0.61 to 0.70 to 0.90 (Fig. 1) as the analyses went from global to functional groups to observation levels. The plot of the model predicted against the measured stomatal conductance (Fig. 1) showed the progressive improvement of the fitting from coarse to fine granularity.
The obtained parameters of the stomatal model at the observation level (Table 3) showed wide ranges of variations across species and sites. Differences in plant and leaf ages, water and nutrient status, instruments used, completeness and accuracy of the data, timing of measurements, and instability of nonlinear algorithms, might all have contributed to these large variations.
The super parameter C l is 3.1361.66 (m 2 s mol 21 ).
Despite that we used the independent uniform prior distributions for the model parameters, we detected significant negative correlations between the fitted parameters K y and g p (r = 20.44, p = 0.002), and between p 0 and m (r = 20.44, p = 0.0015), and positive correlations between parameters K y and p 0 (r = 0.33, p = 0.021), between g p and K i (r = 0.44, p = 0.073). The negative correlation between K y and g p indicated that the stiffness of guard cell structure and efficiency of hydraulic conductance exhibit some kind of coordination. The more efficient vertical xylem transport, the stiffer the guard cell structure to keep the stomata open. It also reflects the balance between stomatal sensitivity to turgor pressure and xylem sensitivity to cavitation.
The obtained soil water potentials at the observation level ( Fig. 2) were plotted as histograms of the three functional types. We found the mean and standard deviations of the soil water While p 0 is the baseline osmotic pressure, the combination of K i and j determines the extent and sensitivity of osmotic adjustment. The three functional groups have similar p 0 in the range of 1.87 to 1.94 MPa. The DCDS group is shown to have the smallest mean K i value (198.8 mmol CO 2 m 22 s 21 ), whereas the SUBS has K i of 319.3. The SUBS has the lowest j (9.9) followed by 23.8 for DCDS and 34.1 for LEGM. Referring to Equation (3), the efficiency of osmotic adjustment is largely determined by the ratio of j to K i since K i is much greater than the maximum net assimilation (,20 mmol m 22 s 21 ) in this case. The obtained model parameters at the global level (Table 4) mostly fall in the range of those at the level of functional groups, except for p 0 of 1.76, which is smaller than any of the three functional groups.

Behavior of the Stomatal Model
Based on the model parameters at the levels of PFT and GLB, calculated stomatal conductance was plotted as functions of soil water potential and dimensionless vapor pressure deficit (Fig. 3), with leaf net assimilation fixed at 6.5 mmol m 22 s 21 (approximately the mean value of net photosynthesis for these groups). Sharp comparison among the three shrub functional groups exists. The lower right corners of the panels depict the maximum stomatal conductance at favorable moisture conditions with y = 20.033 MPa (approximately field capacity) and D = 0.0035. At temperature of 30uC at the sea level, D of 0.0035 means 91% of relative humidity. At these conditions, the subshrub has the highest stomatal conductance of 1.7 mol m 22 s 21 followed by the 1.1 of the LEGM group. However, the DCDS group has the much smaller stomatal conductance of 0.44 mol m 22 s 21 . The GLB has the maximum stomatal conductance of 0.83 mol m 22 s 21 , slightly smaller than the average of the three functional groups. The ranking of these values is largely determined by the K y parameter (Table 4). When moisture conditions are favorable, greater guard cell compliance implies greater aperture of stomata, hence greater stomatal conductance.
The subshrub group has lower tolerance to soil water stress than the two deciduous shrub groups, as the stomatal conductance of SUBS closes at the soil water potential of 22.2 MPa. In contrast, the stomatal conductance of DCDS and LEGM groups decrease to zero at 23.2 and 23.5 MPa, respectively. However, the subshrub group has higher tolerance to vapor pressure deficit than the two deciduous shrub groups (Fig. 3). When the soil water  The tolerance to soil water stress is primarily determined by the osmotic pressure and the compliance of guard cell structure. The SUBS has low tolerance to soil water stress largely because of the greatest K y which makes the stomata sensitive to changes in soil water stress.
The spacing and slope of the contour lines (Fig. 3) represent the relative sensitivities of stomatal conductance with respect to soil water stress and vapor pressure deficit. The response of stomatal conductance to vapor pressure deficit is determined by the ratio of   K y to K soil{to{leaf , so that greater K y but smaller K soil{to{leaf makes the stomatal conductance more sensitive to VPD. A smaller K soil{to{leaf tends to cause insufficient water supply from roots to leaves under greater VPD, and the insufficient water supply decreases the leaf xylem water potential, turgor pressure, and stomatal conductance. A greater K y will make stomata more sensitive to the decrease in the turgor pressure, so that the stomatal conductance decreases faster with VPD. The LEGM and SUBS groups have high sensitivity to vapor pressure deficit, largely because the former has the smallest g p (maximum K soil{to{leaf ), and the latter has the greatest K y values. The DCDS group has a small g p (2.57 mmol m 22 s 21 MPa 21 ), however, it also has a small K y (0.16 mol m 22 s 21 MPa 21 ), so that the stomatal conductance is not as sensitive as LEGM or SUBS. We also used the fitted parameters at the observation level to calculate stomatal conductance of all the observations, and plotted the arithmetic means and one standard deviation below and above against soil water potential (Fig. 4) and scaled vapor pressure deficit (Fig. 5). The results confirmed the findings at the level of functional groups. The average stomatal conductance of DCDS, LEGM, and SUBS close approximately at 24, 23.6, and 23 MPa, respectively, so that the tolerance and sensitivity to soil water stress and vapor pressure deficit of the three functional groups follow the same ranks in Fig. 3.
Finally the comparison of mean K soil{to{leaf as a function of xylem water potential based on the parameters obtained at the observation level (Fig. 6) shows similar characteristics for the three functional groups. However, the calculated y x, 50 values,at which K soil{to{leaf decreased by a half of its maximum value, are 212.3,

Experimental Evidences Connected to the Results
The behavior of stomatal conductance, driven by soil moisture and vapor pressure deficit, depends largely on the stiffness of guard cell structure, and the hydraulic conductance of plant xylem which transports water from soil to leaves. A number of literatures reported the findings of the xylem hydraulic conductivity of various plant life forms.
Kocacinar and Sage [45] experimentally measured xylem hydraulic conductivities of 16 shrub species in west America and central Asia, and found the mean hydraulic conductivities are 3.24610 24 and 0.46610 24 kg m 21 s 21 MPa 21 for C 3 and C 4 species, respectively. If flow paths are about 2 m long, then the hydraulic conductance would be 9.0 and 1.3 mmol m 22 s 21 MPa 21 for C 3 and C 4 species, respectively. The synthesis of measurements on leaf-specific hydraulic conductivity of various plant growth forms [46] reveals that the leaf-specific whole plant conductance is in the range between 0.2 and 20 mmol m 22 s 21 MPa 21 , and that the desert subshrub tends to have higher hydraulic conductance than other growth forms. The measurement of vessel density, diameter, and leaf-specific hydraulic conductance of three Caragana species in northern China [47] showed that their soil-to-leaf hydraulic conductance varied between 0.5 to 20 mmol m 22 s 21 MPa 21 . The data of Iovi et al. [48] for the Mediterranean species showed that soil-to-leaf conductance varies within a wide range up to 40 mmol m 22 s 21 MPa 21 , and the herbaceous group has much higher hydraulic conductance than other groups. They also found a negative exponential decay of the hydraulic conductance with respect to decreased xylem water potential. Our analysis at the observation level yielded g p (maximum potential soil-to-leaf conductivity) in the range from 1.04 to 62.13 mmol m 22 s 21 MPa 21 , and the parameters of the three shrub functional groups are in the range from 2.05 to 5.66 mmol m 22 s 21 MPa 21 , comparable to the above experimental findings. Greater g p is an indicator of greater vessel lumen diameter, but not an indicator of greater vessel density, as the negative correlation between the two quantities found by Poorter et al. [49]. However, the greater diameter of the vessel members in general means higher potential loss of hydraulic conductivity due to cavitation and embolism of vessel members under xylem tension. This tradeoff is reflected in our obtained C l parameter which makes the loss of conductivity proportional to the maximum conductivity. For example, the SUBS group has much higher hydraulic conductance, so that the decrease in conductance with xylem water stress is faster than the other two groups. The result is similar to those found by Iovi et al. [48].
Using finite element analysis, Cooke et al. [50] indicated that a typical stomata aperture width increases from 7 to 15 mm when turgor pressure inside guard cells increases from 0 to 700 kPa. If we assume that stomatal conductance is directly proportional to stomatal aperture and that stomatal apertures of 7 and 15 mm approximately correspond to typical stomatal conductance of 0. The model parameter K y in our analysis is a product of this compliance and stomata density, so that the apparent compliance K y should have more variation range if the variation in stomatal density is considered. We found this parameter varies between 0.088 and 1.482 with a mean of 0.613 at the observation level, and between 0.16 to 1.42 mol m 22 s 21 MPa 21 at the functional group level. The result is comparable to Cook's finding.
From the findings by Leishman et al. [51,52], we know that the non-legume deciduous shrubs have lower stomatal conductance than the other two groups. The most deciduous shrubs in this study (both legume and non-legume) were found as xerophyte growing in relatively shallow soils with high water use efficiency. Their instantaneous photosynthesis can be high when soil water is abundant. However, due to the long-term drying conditions, they have to invest large amount of the assimilated products on constructing xylem for small-diameter vessels with lower g p , which allows them to sustain the excessive drought conditions. Most legume species are thin-leaved Caragana species, which might have something to do with their relatively greater K y than the nonlegume deciduous group.
Finally, the SUBS group (Artemisia), with their less lignified stems and leaves, is shown to have the greatest compliance, and the largest potential hydraulic conductance than the two deciduous shrub groups. Thus the subshrubs are shown here to be typical drought-avoiding in most occasions. In this sense, the subshrubs share some characteristics of herbaceous plants, with stomatal conductance more responsive to variation of moistures in soil and air [18,48]. In other words, the subshrubs behave as an opportunist, so that they have a large stomatal conductance to accommodate the great photosynthesis when moistures are abundant in soil, but choose to close the stomata under severe drought conditions. Implication for Macro-scale Ecosystem Dynamics The dynamics of leaf stomatal conductance is coupled with plant morphological structure and other variables such as leaf nutrient status and boundary layer thickness [28,30]. The present model should be applied and further tested within more complex ecosystem models to improve its connection with external variables.
Our hierarchical analysis showed that the functional group model is a significant improvement over the global level model, and the model at the observation level is a significant improvement over that at the functional group level. This result suggests the parameters at three levels can be used in ecosystems modeling at large, intermediate, and small scales. At regional and global scales in which we need to place all shrubs in one group, the global level parameters can be used. For small patch-scale ecosystems, parameters at the observation level are more appropriate. Modeling at intermediate mesoscales (watershed or landscape) may necessitate distinguishing the properties among shrubs types, thus the results at the level of functional groups may be applicable.
Our findings at the functional group level offered a clue to the equivocal issue in shrub-grass interaction. The process of shrub encroachment into grasslands has been understood as involving nonlinear processes in soil and plant communities, with strong positive feedbacks that consequently lead to the domination of shrubs in the previous grasslands [53]. The current understanding is that once shrubs get established and dominated in the grasslands, it is impossible for the process to reverse to restore the grasslands because of the altered competition at the community scale. However, recent studies gave a number of exceptional contradictory cases. In a 16-year exclosure (excluding from livestock grazing) established at the Erdos Sandland Ecosystem Station in northern China, it was found that a previously dominating shrub (Artemisia ordosica Krasch) decreased by 90% [54], but grasses and forbs increased substantially inside the exclosure. Another example is the Artemisia frigida Willd., a commonly observed shrub species in northern China grasslands. The proliferation of this species in the so-called typical steppes has been considered a result of overgrazing [55]. Enclosure studies in northern China showed that this shrub can be put down with appropriate management. Several studies showed that the enclosures with managed livestock grazing kept the dominance of the species less than 8% in the enclosures of 5, 14, and 25 years, with comparison to the 30% dominance in the heavily grazed sites [56,57,58]. These results contradict the current understanding of the irreversibility of shrub invasion. However, we noticed that the species in the above case studies are Artemisia subshrubs. Our results showed that there is a significant difference in stomatal behavior between the subshrubs and the deciduous shrubs. The deciduous shrubs have great advantages in resistance and tolerance to soil water stress over the grasses, which may contribute to their irreversible encroachment into grasslands. The subshrubs are less resistant and tolerant to soil moisture stress than the shrubs, and thus are less advantageous in the competition with the grasses. Our hypothesis is that when the grazing pressure is heavy, the soil cannot hold enough water because of the decreased grass roots biomass. The subshrubs start to gain dominance. However, when grazing pressure is reduced, root biomass start to increase, allowing fast infiltration of water. Consequently the soil layers tend to hold more water to allow better grass growth, and the advantage of the subshrubs over the grasses is weakened. Therefore, with appropriate management, the invasion of the subshrubs into grasslands might be reversed and the health of the grassland ecosystems might be recovered. This hypothesis has to be tested in more strictly designed experiments.

Supporting Information
Table S1 Data measured by the authors group with Licor 6400 portable gas analyzer. (XLSX)