Plant functional traits and climate influence drought intensification and land-atmosphere feedbacks.

The fluxes of energy, water, and carbon from terrestrial ecosystems influence the atmosphere. Land-atmosphere feedbacks can intensify extreme climate events like severe droughts and heatwaves because low soil moisture decreases both evaporation and plant transpiration and increases local temperature. Here, we combine data from a network of temperate and boreal eddy covariance towers, satellite data, plant trait datasets, and a mechanistic vegetation model to diagnose the controls of soil moisture feedbacks to drought. We find that climate and plant functional traits, particularly those related to maximum leaf gas exchange rate and water transport through the plant hydraulic continuum, jointly affect drought intensification. Our results reveal that plant physiological traits directly affect drought intensification and indicate that inclusion of plant hydraulic transport mechanisms in models may be critical for accurately simulating land-atmosphere feedbacks and climate extremes under climate change.

Earth’s land surface affects the atmosphere through exchanges of energy, water, and carbon (1, 2). A rapidly growing literature has documented extensive soil moisture feedbacks to subsequent precipitation and temperature in regions across the globe (2–7). These soil moisture feedbacks are critical land–atmosphere interactions in a changing climate because they likely play a strong role in drought intensification, heatwaves, and climate extremes (5, 8). For example, the number of hot days in many regions is strongly associated with preceding precipitation deficits (8) and droughts that occur with elevated temperatures, often called “hot droughts,” can be particularly damaging to ecosystems (9). Earth system models capture some of the central physical processes, such as latent and sensible heat fluxes and their drivers. Given this mechanistic strength, earth system models have been used extensively to quantify soil moisture feedbacks by running different scenarios with differing soil moisture boundary conditions and predicting how different boundary conditions affect water and energy fluxes (2, 10, 11). However, the representation of plant transpiration in earth system models during drought is notably coarse and could lead to large uncertainties in simulations of soil moisture feedbacks (12–14). Thus, to understand and predict the spatial patterns and impacts of future climate extremes, it will be crucial to represent the effects of vegetation processes on soil moisture feedbacks because vegetation feedbacks may mediate drought intensification (10, 15).

The influence of soil moisture on land–atmosphere interactions and drought intensification involves complex and multiscale dynamics of energy partitioning and circulation. A rich theoretical literature has predicted that soil moisture feedbacks should be more prominent in soil-moisture-limited (e.g., generally hotter and drier) rather than in energy-limited (e.g., solar radiation) regions (2). In addition, the influence of soil moisture on plant transpiration, which is the dominant component of latent heat (LH) fluxes in vegetated areas (e.g., forests) and plays a critical role in rewetting the convective boundary layer for the generation of condensation and precipitation, is critically important for land–atmosphere interactions (2). A prominent soil moisture feedback occurs when decreasing soil moisture suppresses plant transpiration, which drives declines in LH fluxes and increases in sensible heat (SH) fluxes, local air temperature (T), and vapor pressure deficit (VPD) (2). These soil moisture feedbacks can drive “drought intensification” where the declines in soil moisture lead to increasingly severe conditions for other components/drivers of drought, such as higher VPD, and can be quantified as the correlation between antecedent soil moisture and subsequent T/VPD/SH. The centrality of transpiration in surface energy exchange indicates that plant physiology and physiological traits, in addition to vegetation type and structure, may play a critical role in soil moisture feedbacks (13, 16, 17). The 2003 European heat wave provides preliminary evidence for trait-based effects on soil moisture feedbacks, as different land cover types, including differing forest composition, were associated with local temperature anomalies (18–20).

First principles and physiological modeling studies predict, and observational studies have confirmed, that plant functional traits mediate the magnitude of water fluxes and how water fluxes change in response to falling soil moisture and increasing atmospheric vapor pressure deficit. Multiple studies have highlighted that plant-level water fluxes are influenced by key plant traits such as maximum photosynthetic capacity and the vulnerability of plant water transport/hydraulic tissues to increasing water stress (21–23). Variation in water fluxes due to species composition and a diversity in trait strategies is apparent at the whole forest ecosystem scale as well (24). Thus, plant traits that govern leaf gas exchange and plant hydraulic transport could potentially influence drought intensification and mediate future climate extremes. For example, communities of plants with traits that cause rapid depletion of soil moisture, such as high photosynthetic and hydraulic transport rates, could lead to a stronger drought intensification feedback. However, the extent to which plant functional traits influence soil moisture feedbacks, and which functional traits are most important, remains largely unknown.

Here, we examine the soil moisture feedbacks to drought intensification—diagnosed as the correlations between antecedent soil moisture and subsequent T, VPD, or SH (). We use a network of 40 eddy covariance towers in temperate and boreal forests, multiple datasets of plant functional traits, a mechanistic vegetation model, and a recent satellite-derived quantification of soil moisture feedbacks to ask: 1) How does the strength of feedbacks between soil moisture and T, VPD, or SH vary across climate gradients in temperate and boreal forests? 2) To what degree do differences in climate and/or plant functional traits mediate drought intensification? 3) Which plant functional traits emerge as potentially important for land–surface feedbacks affecting drought intensification? 4) Does plant diversity, as a proxy for trait diversity, explain spatial patterns of continental-scale soil moisture–precipitation feedbacks?

Results and Discussion

We found substantial variability in soil moisture feedbacks across flux sites in the temperate and boreal forests, with strong variation in the temporal correlations between previous soil moisture and subsequent temperature, VPD, and sensible heat flux during the growing season (Fig. 1 and ). However, the signs and magnitudes of feedbacks were generally consistent within a given site across different time scales and meteorological response variables. This site-specific consistency of feedbacks is largely expected given the lack of large changes in vegetation composition or structure at individual sites. Positive feedbacks of drought intensification—negative correlations between soil moisture and subsequent T/VPD/SH—were observed in all biomes, but were more prevalent at warmer sites (Fig. 1). In this scenario, low soil moisture could feed back to increase T/VPD/SH, which further increases the evaporative demand and further decreases soil moisture. Negative feedbacks, where decreasing soil moisture preceded decreasing T/VPD/SH, were predominantly observed at cooler sites. These patterns were robust across a diversity of time-integration windows ranging from 1 to 14 d and to the method used to define the plant-active/growing season across sites (Fig. 1 and ).

Fig. 1.

Substantial differences in land–atmosphere feedbacks are present across 40 temperate and boreal forest flux-tower sites. Pearson correlation coefficient between soil moisture in a previous period (1-, 7-, or 14-d average, indicated in green at top of each column; see ) and temperature, VPD, or SH in a subsequent period (1-, 7-, or 14-d average, indicated in purple at top of 3-column groups). All factorial combinations of averaging periods yield 9 columns per variable. Each flux site is a row (for full names, see ) and are organized from coldest MAT (Top) to warmest (Bottom). Colors indicate the sign (positive/negative) and color/sizes of dots indicate the strength of the correlation. Averaging windows without enough data are shaded in gray.

Climate and plant functional traits explained between 54% and 67% of the observed cross-site patterns in drought intensification feedbacks (Fig. 2). Using a suite of model selection algorithms and variance decomposition, we determined that climate and traits jointly explained about half of this variation; climate alone explained about 33% and traits alone about 17% (Fig. 2). Mean annual temperature (MAT) of the site was the most important climate variable. For functional traits, we examined the community-weighted mean and SD in maximum photosynthetic rate (Amax), specific leaf area (SLA), wood density, water potential at 50% loss of branch hydraulic conductivity (P50), minimum water potential (PsiM), and hydraulic safety margin, which capture a range of key resource acquisition and water transport strategies (). Of these traits, the community-weighted mean Amax and SD of P50 were among the most critical (Fig. 2 and ). While traits were less important in the combined Akaike Information Criterion (AIC) weight across all models (), a result that is consistent with the variance decomposition analysis (Fig. 2), they nevertheless were important predictors and displayed similar magnitudes in their standardized effect as mean annual temperature (Fig. 2).

Fig. 2.

Climate and plant functional traits both influence land–atmosphere feedbacks and affect drought intensification. Groupings show the feedbacks between soil moisture (SM) and Temp, VPD or SH. (A) Variance explained by the full model (white), variance jointly explained by climate and traits in variance decomposition (tan), explained by climate alone (yellow) and by traits alone (green). (B) Standardized variable coefficients across all models for MAT and community-weighted traits of: Amax.m, P50.sd, Amax.sd, and PsiM.m. See for trait names.

Extensive previous research has documented the impacts of aridity or vegetation presence versus absence, such as deforestation, on surface energy fluxes (2, 7, 25–29). Here, we demonstrate that plant physiological characteristics are important for land–atmosphere interactions beyond simply vegetation presence/absence or structure. We further show that plant functional traits, particularly those related to maximum carbon uptake and water transport, significantly influence soil moisture feedbacks across a diversity of temperate and boreal biomes. Using a hydraulic trait-based forest model (21), we performed several sensitivity tests to explore the potential mechanisms that underlie plant trait-mediated soil moisture feedbacks (). We found that higher community-weighted mean Amax increased the maximum transpiration/LH fluxes during wet periods, led to faster soil moisture drawdown and sharper declines in LH fluxes during dry periods, and yielded a stronger drought intensification feedback (), consistent with the cross-site regression coefficient for Amax in Fig. 2. Higher diversity of species’ stem hydraulic vulnerability traits (P50 values) in an ecosystem, a hydraulic trait thought to capture drought tolerance through the water potentials that can be tolerated in the xylem (30), tended to decrease LH fluxes in wetter periods but sustain LH fluxes much longer during a dry-down, leading to weaker drought intensification feedbacks. Results from the mechanistic vegetation model analysis are consistent with the cross-site trait-based patterns in Fig. 2 (). Overall, the combination of community-weighted average Amax and diversity in P50 is consistent with the soil moisture–drought intensification feedback framework, although further research is needed to fully understand the detailed processes and mechanisms.

The statistical and model selection procedures performed here identify important and explanatory variables, but have several limitations and uncertainties. The multivariate models demonstrate statistical linkages, which in complex situations such as the analyses here are often not straightforward to interpret. The statistical linkages between site climate and trait metrics can reveal patterns and associations, but are subject to uncertainties in the input variables (e.g., scaling individual species’ traits to whole communities, appropriate climate variables and uncertainties in using gridded climate datasets), output variables (e.g., quantification of feedbacks), and their statistical connections (e.g., unaccounted-for confounding variables). In particular, more work is needed using ecosystem-level experiments and process-based land surface model simulations to examine the processes and mechanisms through which community-level plant traits might directly affect feedbacks.

We next examined if plant diversity mediated regional-scale soil moisture feedbacks and drought intensifications. A recently published rigorous quantification of the soil moisture–precipitation feedback over the continental United States observed stark differences between the western and eastern United States (Fig. 3 ) (6). As a preliminary exploration, we tested if the interaction between tree species diversity and local climate might explain spatial patterns in these feedbacks. We consider only woody-plant-dominated grid cells and included grid cell leaf area index (LAI) as a covariate to account for vegetation structure differences across sites. We used tree species richness as a predictor variable because community-weighted trait maps are not available for our primary functional traits at these scales. Tree species richness and LAI explained 36% and 57% of the spatial variation in soil moisture feedbacks (P < 0.00001 for both feedback estimates; Fig. 3). All else equal (e.g., accounting for mean soil moisture), regions with lower tree diversity tended to show stronger signals of drought intensification; i.e., that below average soil moisture anomalies were associated with decreased probabilities of subsequent precipitation. This finding was robust to accounting for the effect of spatial variations in average soil moisture across the United States (). Furthermore, we observed a significant interaction between richness and LAI that indicated a saturation of the effect of species richness on soil moisture feedbacks at high LAI values (Fig. 3 and ).

Fig. 3.

Tree species richness patterns are associated with soil moisture–precipitation feedbacks. A and B are the quantification of the “dry” soil moisture–precipitation feedback from ref. 6 using 2 different satellite soil moisture datasets where values below 1 denote low soil moisture associated with lower probability of subsequent rainfall and above 1 low soil moisture associated with higher probability of subsequent rainfall. (C) Tree species richness (no. species). (D) The relationship between the soil moisture feedback map in A and tree species richness from (C) with the blue line showing the best fit generalized additive model (deviance explained = 51%) and 95% confidence intervals.

Our results demonstrate that plant functional diversity can influence drought intensification and climate extremes. Current land surface models likely capture many of the key climate effects on soil moisture feedbacks, but do not represent the effects of decreasing soil moisture on transpiration with mechanistic fidelity, potentially adding substantial uncertainty to carbon cycle feedbacks (12). Our findings emphasize that simulating diversity in plant hydraulic transport and its mechanistic effects on stomatal conductance (e.g., refs. 13 and 21) in land surface models may be crucial for capturing drought intensification and climate extremes because this would provide a mechanistic linkage between soil moisture and plant transpiration. Fortunately, large-scale hydraulic-enabled models are actively being developed (31, 32). This improved simulation of soil moisture feedbacks, including diversity in plant hydraulic strategies, will be critical for providing rigorous predictions of climate extremes and climate impacts.

Methods

Eddy Covariance Site Data.

We used the tier 1 (publicly available) eddy covariance data from the recently released FLUXNET2015 dataset. This contains data from >200 eddy covariance sites from around the world with energy, water, and carbon fluxes and meteorological data and has undergone a standardized set of quality control and gap filling (33). We used the daily data for all analyses here, as our goal was to quantify multiday to multiweek soil moisture feedbacks. Sites that met the following criteria were included in our analyses: 1) forest sites, 2) no disturbance within 10 y before the onset of eddy covariance measurements, 3) on-site soil moisture measurements, and 4) available plant functional trait data for at least 2 functional traits for 80% of the tree community (see below). This led to a final list of 40 sites and 352 site years across 4 forest biomes—deciduous broadleaf forest, evergreen needleaf forest, mixed forest, and evergreen broadleaf forest—which were almost entirely in temperate and boreal regions ().

Trait and Climate Data.

Based on an extensive review of the literature and publications from individual eddy covariance sites, we compiled lists of the dominant tree species present in the footprint of each tower and, where possible, the relative abundance of these species. Abundance data were available for around 60% of sites. We then compiled the functional traits for species present in these sites, drawing on the Global Wood Density Database (34) for wood density data, the dataset presented in ref. 35 for light-saturated maximum photosynthetic rate and SLA, and the Xylem Functional Traits dataset (36) for data on the water potential at 50% loss of stem hydraulic conductivity and the hydraulic safety margin. These traits influence a diversity of important plant functions, including maximum gas exchange rate and the vulnerability of water transport to drought conditions. Furthermore, these traits are actively being incorporated into next-generation land surface models (13, 31, 32, 37) and thus their importance for drought intensification is useful to test. At each site, we calculated the community-weighted (i.e., composition weighted) trait mean and SD for each trait with values for more than 80% of the plant community. The mean provides an estimate of the average plant strategy in the ecosystem and the SD provides insight into the diversity of strategies/traits present.

In addition, we compiled the reported age of the canopy trees and the mean annual temperature and mean annual precipitation reported as metadata in the FLUXNET2015 dataset for each site as important climate variables that might influence cross-site patterns. We used the site latitude and longitude to extract the growing season (defined here as June–August for northern hemisphere and December–February for southern hemisphere) temperature and precipitation from Climatic Research Unit data (38) for the 1980–2008 period. We further extracted annual and growing season potential evapotranspiration (PET) for each site (averaged over the 1980–2008 period) from the Sheffield dataset (39) and included site mean annual PET, mean growing season PET, mean annual precipitation minus mean annual PET, and mean growing season precipitation minus mean growing season PET as potential climate predictor variables. Growing season temperature was also predictive of feedbacks, but was less strongly correlated with feedbacks than MAT () and thus was not included in the final models. See for all predictor variables.

Soil Moisture Feedback Methods.

We restricted our data to the plant-active season using 2 different methods. In both, we removed all data where a quality flag indicated that a measurement of soil moisture, temperature, vapor pressure deficit, or sensible heat may have been flawed. We further verified this by plotting annual time courses of all variables for all sites to ensure that all data were within reasonable bounds. In the first method, we implemented a daily average temperature threshold of 15 °C for temperate and 10 °C for boreal sites and further removed 1 wk on either end of each plant-active season. In the second method, we implemented the cross-site method proposed by Gu et al. (40) and considered a growing season to be between the “stabilization day” and “downturn day” of annual GPP (via the nighttime partitioning method) time series. These methods yielded largely similar results (Fig. 1 and ), and thus we present the temperature threshold (first) method in the main text.

We followed a number of previous studies in estimating the strength of soil moisture feedbacks and drought intensification as the correlation between previous soil moisture and subsequent T, VPD, or SH (3, 8, 41, 42). Quantifying land surface feedbacks is inherently challenging because 1) observational data are typically at much smaller spatial scales (e.g., an eddy covariance footprint) than the relevant atmospheric processes, 2) effects may arise from autocorrelation of dependent variables (e.g., synoptic patterns in precipitation), and 3) lagged correlations among variables may not mean causation (2, 42). Challenge 1 is largely unavoidable for site-level analyses where plant trait data are available, but we believe our approach is still reasonable for several reasons. First, this limitation has not precluded the detection of likely soil moisture feedbacks in previous studies (3, 43). Second, inferring feedbacks from site-level data are reasonable if the sites are representative of broader regional vegetation types (thereby capturing the main mechanisms for regional-scale feedbacks), which is largely the case of our focal flux towers. Finally, we analyze broader continental-scale datasets as well, which show a consistent pattern with flux-tower data (Fig. 3). In addition, soil depth or properties such as soil texture, rooting depth, and belowground plant or mycorrhizal traits could all influence land–atmosphere interactions, but were not considered here due to a lack of available data. Challenge 2 is a much larger problem for the detection of soil moisture–precipitation feedbacks, and analyses focusing on soil moisture–temperature are likely to be more robust than soil moisture–precipitation feedbacks (2, 42). This is because soil moisture–temperature feedbacks are directly influenced by surface radiation budgets (e.g., Bowen ratio), are less likely to be driven by an external variable (e.g., sea surface temperatures), and are less influenced by statistical challenges of precipitation persistence (2). We further aimed to minimize the problems of challenge 2 by using an array of time windows and removing seasonal patterns in variables (see below). Challenge 3 is an important caveat in all soil moisture feedback analyses, and causality can only be fully prescribed with mechanistic model simulations that vary soil moisture boundary conditions, something that cannot be done with observational datasets. However, some statistical techniques assess “predictive causality” (e.g., Granger causality), which is more rigorous than simple correlation, and is used for the satellite-based feedbacks assessed here (6).

presents an illustration of how we quantified feedbacks by showing the calculation from a given day i at a given site. Within each growing season at each eddy covariance site, we first calculated all nonoverlapping 1-, 7-, and 14-d averages for T, VPD, and SH. From each period, we then calculated the preceding 1-, 7-, and 14-d period averages for soil moisture. Where preceding soil moisture data were not within the plant-active season, those dates were dropped from the analysis (for example, day 2 of a growing season would have a value for the preceding 1-d period, but not the preceding 7- or 14-d period, and thus would not be included in those analyses). In the case of data gaps, a majority of days with data were required to be included. The 3 time periods of dependent variables (T, VPD, SH) and 3 time periods of the independent variable (preceding soil moisture) led to 9 combinations over which we calculated Pearson correlation coefficients, where negative correlations indicate “drought intensification” because declines in soil moisture are correlated with subsequent increases in T/VPD/SH. We used this suite of different time periods to 1) ensure that the patterns in drought intensification that we observed were robust to different methods, and 2) to capture effects at longer timescales where autocorrelation of T, VPD, and SH are much lower (44). In addition, we calculated these correlations after removing a seasonal pattern from all variables via fitting a quadratic regression for each year and performing the correlations on the residuals for a random subset of 10 flux sites. The feedback correlations on the deseasonalized data were quite similar to the correlations used in Fig. 1 (R2 = 0.65, P < 0.001), indicating that seasonal or synoptic patterns are not likely behind site differences. Because the correlations were quite consistent across different time period combinations (e.g., Fig. 1), we averaged all 9 time period combinations to yield an average feedback strength for each site for the model selection analyses.

Model Selection.

We performed model selection to determine the most important predictors in cross-site patterns of soil moisture feedbacks. See for all predictor variables. Because multivariate analyses can be confounded by collinear predictor variables, we first removed collinear predictors following previously published methods (45) by calculating the pairwise Pearson correlation coefficients among all predictor variables. Whenever 2 predictor variables exhibited a correlation of >0.5, each was then compared against the dependent variable and the predictor with the lowest correlation with the dependent variable was dropped. The final model set included 6 predictor variables and was consistent across all 3 dependent variables: mean annual temperature, community age, and community-weighted traits of: mean maximum photosynthetic rates (Amax.m), SD of the water potential at 50% loss of stem hydraulic conductivity (P50.sd), SD of maximum photosynthetic rates (Amax.sd), and mean minimum water potential experienced (PsiM.m).

We then used the model selection technique of “all possible models” that calculates all potential combinations of predictor variables and ranks them by Akaike Information Criterion (46). This approach is considered superior to and more robust than stepwise techniques and has been used in soil moisture feedback analyses previously (6). From this extensive family of models, we then calculated the standardized variable coefficients, Akaike variable weights, and the prevalence of each variable in the top models (defined as all models within 3 AIC of the best model). No single model emerged as the most parsimonious model (e.g., <−3 AIC from another model), and thus presenting the Akaike variable weights and prevalence of variables in the top models provides a holistic picture of the importance of each predictor variable.

Finally, we performed variance decomposition on the most complex model within the top models for each model selection analysis (i.e., soil moisture-T, soil moisture-VPD, or soil moisture-SH). This model had the form:Variance decomposition allowed us to directly quantify the amount of variance explained by the full model (all predictor variables), climate (i.e., mean annual temperature) alone, traits alone, and climate and traits jointly. This was done by calculating the semipartial correlation—the correlation of 2 variables (e.g., Feedback∼MAT) with variation from all other variables removed only from the dependent variable (e.g., removing the variation of Amax.m and P50.sd from Feedback above). This was done separately for climate (MAT alone) and traits (Amax.m and P50.sd) to calculate the R2 from each semipartial correlation and the remaining R2 is considered the variance explained by the joint variables (i.e., shared between the 2 sets of variables). See for mechanistic simulations of trait differences in a hydraulic model.

Continental US Soil Moisture Feedback Analysis.

The analysis of Tuttle and Salvucci (ref. 6; results shown in their figure 1A and B) identified the causal influence of soil moisture on the occurrence of next-day precipitation. Precipitation is often autocorrelated, with precipitation events occurring in clusters over timescales of weeks, months, or years. This autocorrelation may be due to local (e.g., soil moisture feedbacks) or far-field influences (e.g., synoptic weather patterns, sea surface temperatures). No matter the cause of the autocorrelation in precipitation, this signal will be passed on to soil moisture, due to the direct coupling between the two (i.e., when it rains, soil moisture increases). Therefore, any lagged correlation between soil moisture and subsequent precipitation will be contaminated by autocorrelation in precipitation, regardless of the cause (47).

To isolate the effect of soil moisture on the occurrence of next day precipitation, Tuttle and Salvucci (6) used generalized linear modeling (specifically, probit regression) within a Granger causality framework. Precipitation occurrence (a binary variable: rain vs. no rain) was modeled as a function of sinusoids that represented various modes of interannual variability and seasonality, and up to 4 d of past precipitation (see ref. 6 for details). For each given location, models with all possible combinations of the independent variables were separately evaluated against observed precipitation (i.e., “all possible models” or “all possible regressions”), and the model with the combination of independent variables that yielded the lowest AIC was chosen as the best “restricted” model for that given location. Then, a similar probit regression model was evaluated, with 1-d lagged, seasonal soil moisture anomaly as the independent variable (along with 4 d of lagged atmospheric pressure) and the residual precipitation occurrence from the “restricted” model as the dependent variable. This second, “unrestricted” model tested whether or not soil moisture exerted a statistically significant influence on next-day precipitation occurrence, after accounting for seasonality, interannual variability, and precipitation persistence. The ratio of the predicted precipitation probability from the models with and without soil moisture as an independent variable provided a quantification of the impact of soil moisture on the occurrence of next-day precipitation (i.e., results shown in ref. 17, Fig. 1 ). The precipitation data used in Tuttle & Salvucci (2016) are from the North American Land Data Assimilation System, Phase 2 (48). 2 soil moisture datasets derived from Advanced Microwave Scanning Radiometer for the Earth Observing System satellite observations ref. 49 in Fig. 1; ref. 50 in Fig. 1) were evaluated against the precipitation data. The analysis was restricted to the contiguous United States from June 2002 to June 2011 and to months without freezing air temperatures.

Because community-weighted trait maps are not currently available for our key traits at continental scales, we used a tree species richness map for the continental United States (51). In addition, we drew upon a biome map (52), gridded monthly soil moisture data from the Global Land Data Assimilation System (53), and gridded data of leaf area index from Moderate Resolution Imaging Spectroradiometer from July 2005 (54), which was chosen to represent the peak growing season in the continental United States, and central in time for our analysis in the soil moisture feedback datasets. We included these ancillary datasets either as masks or covariates in our analyses (see below).

For each of the 2 estimates of soil moisture–precipitation feedback, we restricted our analyses to only grid cells dominated with woody vegetation (cover types 1–9) and removed grid cells with values of 1 (no statistically significant feedback, which constituted 38% of woody vegetation grid cells). We then performed a multivariate ordinary least squares regression of the feedback as a function of species richness, average leaf area index, and their interaction term, which was the model that had the lowest AIC. We tested for spatial autocorrelation in our data with a Moran’s I test (55), which was not significant. Given the saturation of the feedback at higher values of richness, we further fit a generalized additive model with feedback as a function of richness.

To ensure that our results were not driven by a correlation of average soil moisture with the soil moisture feedback estimates, we also ran analyses where we first regressed the soil moisture feedback against average soil moisture from 1989 to 2008 and then used the residuals of this regression against tree species richness, which was still significant (R2 = 0.19, P < 0.00001).

Statistics.

All statistics and analyses were conducted in the R computing environment. Specifically, we used the dredge function in the MuMIn (56) package for model selection, the ppcor (57) package for variance decomposition, the nlme (58) package for generalized additive model analyses, and the RnetCDF (59), raster (60), and rworldmap (61) packages for spatial analyses and map plots. We examined all ordinary least squares linear regressions for to ensure that they met key assumptions using quantile and other diagnostic plots.