We compiled data on vegetation response to Tamarix removal and associated environmental variables from three published studies with sites that were distributed across the Upper Colorado, Lower Colorado, and Rio Grande river basins in the southwestern US (González et al. 2017a, b; Sher et al. 2018; Fig. 1). These sites covered a range of abiotic environmental conditions as well as biotic starting points (González et al. 2017a). In total, there were 260 sites where Tamarix had been subject to active removal and/or biocontrol and 132 reference sites. Removal treatments included all that are typical for this species, including biological control, prescribed burning, extraction or mulching with heavy machinery, chainsaw cutting (usually with herbicide applied to stumps afterwards to prevent resprouting; referred to as cut-stump), and biological control (Gonzalez et al. 2017b).
Independent variables in our analysis included both aspects of management and abiotic features measured at each site (González et al. 2017a, b; Sher et al. 2018; Table 1). “Site type” was our categorical variable used to explore the differences in response among the four removal methods and two categories of non-removal, reference sites. When more than one removal method was used, the one considered highest disturbance was used (González et al. 2017a). Of the 132 non-removal, reference sites, 53 were “desirable” (i.e., dominated by native vegetation; mean Tamarix cover = 0.84%, maximum cover = 8.6%) while 79 were “undesirable” (as deemed by managers due to Tamarix abundance, mean = 46.6%; González et al. 2017a). Other data on the management of sites included whether herbicides (of any type for any purpose) were applied, and whether active introduction of desirable vegetation by planting or seeding was used. Details regarding herbicides or planting techniques/approaches were usually unavailable and so could not be included in the models. Environmental data available for these sites included information about climate, geography, hydrology, and soils (Table 1). Some variables in the source data had missing data for some rows; if more than 20% of data was missing, we excluded the variable from analysis, and if less than 20% was missing we filled in missing values with the row median. The variable “time since Tamarix removal” was not included both because it was not available for all sites, and it was highly non-significant for explaining Salix cover (F1,198=0.011, p = 0.92).
Table 1: List of all explanatory variables included in analysis.
Management & vegetation variables
|
Change in Tamarix cover
|
Tamarix cover before removal
|
Site type (positive and negative references, restoration methods)
|
Herbicide used at any point (including cut-stump treatment)
|
Active introduction of native plants
|
All-Salix/S. exigua cover before removal
|
Climate variables
|
Minimum average temperature (average of 30 years prior to removal date)
|
Maximum average temperature
|
Average precipitation
|
Precipitation year sampling
|
Minimum temperature year sampling
|
Maximum temperature year sampling
|
Geography variables
|
Longitudinal slope
|
Elevation above sea level
|
Distance to nearest paved road or railway
|
Distance to nearest unpaved or paved road/railway
|
Hydrology variables
|
Distance to river’s edge
|
River width
|
Elevation above water level
|
Soil variables
|
Electrical conductivity
|
pH
|
Percent sand
|
For our dependent variables, we used two measures of Salix response to Tamarix removal treatments: change in cover over time and final year cover (i.e., final year of monitoring for a given site). By testing the response of both change and “final” measures of Salix cover, we were able to investigate the response of Salix both by 1) considering Salix response as a function of the efficacy of Tamarix control and 2) by using presence of Tamarix at a single point in time as a predictor. The latter is the more common approach used, but it only tells us whether there is a correlation between cover of the two tree species, which could be due to abiotic environmental factors, not competitive exclusion. Whereas the former analysis of change is a better test of the impact of removal (Sher et al. 2018).
Our first variable, “final” cover, represented the last time data was collected for each site (n = 392). Cover was a standardized measure of percent abundance, calculated for each plant species found at each site (González et al. 2017a); because SWFL can nest in any Salix species, we used the sum of cover for all Salix species (henceforth “all-Salix”; Table 1). We also separately examined Salix exigua (aka coyote willow; sandbar willow). S. exigua was the single most abundant Salix species across sites (Fig. 2) and also has several functional differences from other Salix species, including the ability to reproduce asexually via resprouting (Rood et al. 2011).
Change in cover (for both all-Salix and S. exigua) was calculated as the difference between the final year measurement and a control site measurement. Of the 392 total sites, 180 had “before-after” controls, with data for pre- and post-restoration conditions, and could thus be analyzed in terms of vegetation cover change over time (for more detail on this variable see Sher et al. 2020).
It is important to note that no SWFL populations were present at any of the sites used in this study; since the SWFL is a federally endangered bird, it would be illegal to remove Tamarix at a site known to support breeding populations. Thus, this study is assessing potential (not current) habitat, as well as Salix cover in general as an intrinsic restoration goal. The range of the SWFL overlaps some of our more southern study sites (specifically in southern Nevada, Arizona, and New Mexico; Hatten 2016), but its range is likely to eventually encompass all sites as the SWFL adapts to climate change. Many bird species have already shifted their ranges northward due to increasing temperatures (Spence and Tingley 2020). Thus, determining habitat for SWFL beyond the current breeding locations and range is important for understanding the ongoing effects of Tamarix removal.
Statistical Analysis
To answer our first question about whether there was evidence that Salix cover increased with decreasing Tamarix cover, we investigated Salix change over time as a function of Tamarix change over time with a linear mixed-effect model (LME; Zuur et al. 2009) using river reach as a random effect to account for spatial auto-correlation; this was done for each treatment type (biological control alone, cut-stump method, heavy machinery and burning, and both undesirable and desirable controls) separately (respectively n = 33, 67, 33, 14, 11, 22) and combined (n = 180). To explore whether the Salix response to change in Tamarix cover was moisture-dependent, we then used we used linear mixed models with backward stepwise selection of the independent variables that included both main effects of Tamarix change over time, hydrological variables, and interactions between them, with river reach as random effect (n = 180). The backwards selection for all linear mixed models was done using Akaike’s Information Criterion (AIC) to identify significant variables to be used in the final model.
To answer our second question regarding which variables best predicted change over time of Salix more generally, we used linear mixed models with backward stepwise selection to predict response of both all-Salix and S. exigua cover using all independent variables, with river reach as a random effect (n = 180). Finally, we identified the conditions that predicted Salix and S. exigua cover regardless of whether Tamarix was removed or not, using all independent variables and all sites (n = 392). Our sample size was too small relative to the variance among sites to consider all possible interaction terms of our independent variables, so only interactions that tested our specific hypothesis about the role of water were analyzed, separately from these general explorations. In addition to exploring the interactions between change in Tamarix and hydrological variables, we also tested interactions between herbicide use and climate variables, because herbicide was selected in all change-over-time models but showed a high degree of variability.
Several of the environmental variables exhibited multicollinearity (Appendix A); elevation above sea level, average precipitation, and average minimum and maximum temperatures were all correlated predictors of climate. Likewise, distance of each site to a paved road was correlated with distance to any road and to the site’s longitudinal slope. To address this issue, each model was first created using all independent variables with a maximum likelihood calculation, then stepwise selection using AIC. For change-over-time analyses we used the final-year values for all abiotic environmental variables, as they did not change significantly over the duration of the study (results not shown). The same methods were used for both the final-year and change over time models, but with some different starting variables due to differences in available data between final-year-only and change-over-time datasets (e.g., Tamarix change over time vs. Tamarix final cover as explanatory variables; Appendix C). See Appendix B for starting models: 21 for final cover models and 42 for change over time. There were more of the latter since we could build models with Tamarix cover before removal and alternatively with Tamarix cover change over time as explanatory variables, whereas the final cover models only had Tamarix cover before removal. Each explanatory variable was scaled (constrained to (-1, 1)). Final mixed models were constructed using restricted maximum likelihood calculations. For each model we calculated conditional and marginal adjusted R2 and tested model significance using a log-likelihood test to compare with a null model with only the random factor (river reach) as an explanatory variable. Finally, in order to determine whether any variables significantly explained Salix cover variation on their own, we conducted univariate tests of each variable on final-year cover, using river reach as a random variable.
All analyses were performed using R version 3.3.6 with RStudio version 1.2.5042 using the functions "lme” and “stepAIC” in the nlme package (Pinheiro et al. 2020).