Plot selection and manipulations
Our experimental design consisted of 15 - 1 ha blocks spatially
distributed across the study area selected to have similar total
percentage vegetation cover within a block, but to vary in the
proportion of perennial grass and mesquite cover among blocks.
Potential experimental block locations were identified using two
approaches that used either percentage bare ground or change in
vegetation through time; we then combined the two approaches for
final plot selection. First, percentage bare ground across the
grassland-shrubland gradient was estimated for contiguous 20 m x
20 m grid cells using an NDVI analysis of a 2003 QuickBird
imagery. Each grid cell was classified into one of five classes of
bare ground cover (0 to <20%, 20 to <40%, 40 to <60%, 60
to <80%, 80 to 100%) used to delineate contiguous 1-ha regions
for each class. Thirty grid cells from each class (30 cells x 5
classes = 150 cells) were randomly selected as the preliminary
pool of block locations. These grid cells were then reclassified
as either low (0 to 30%), medium (33 to 67%), or high (67 to 100%)
bare ground. Second, the change in NDVI from 2003 to 2010 from
another Quickbird image was used to capture vegetation change
following a series of years with above-average precipitation (2004
to 2010; (52)). The 2010 Quickbird image was separated into the
same 20 m x 20 m grid cells as the 2003 image used for bare ground
classes, and NDVI was calculated for each grid cell. The
difference in NDVI from 2003 to 2010 was used to classify each
cell into either low (0-33%), medium (33-67%), or high percentage
vegetation cover (67-100%). Ten potential locations were randomly
identified for each of the nine combinations of percentage bare
ground and change in NDVI classes. Four grid cells were selected
in each of the nine classes for a total of 36 locations to ensure
that candidate block locations were arrayed across the study area.
Ground-truthing of bare gap size distribution using a standard gap
intercept method was conducted for the final pool of 30 candidate
block locations (44). The final 15 block locations were chosen to
characterize a gradient in grass and mesquite shrub cover. At each
block location, four 15m x 25m plots were identified with similar
percentage grass and shrub cover with a focal patch in the center.
Plots were arrayed perpendicular lengthwise to the dominant
southwest to northeast wind direction, and were >15 m apart to
minimize between-plot treatment effects.
Experimental design
Each plot within a block was randomly assigned to one of four
connectivity treatments: (1) plant scale where all mesquite plants
within and surrounding the focal patch were killed in place to
modify competitive interactions between woody plants and recovery
of perennial grasses and other herbaceous plants with no direct
effects on horizontal transport by wind and water, (2) patch scale
where ConMods were located in bare soil interspaces between plants
in the focal patch to reduce gap size and to modify transport of
water, soil, nutrients, litter, and herbaceous seeds, (3) both
patch- and plant-scale manipulations were conducted in each focal
patch, and (4) no manipulations [controls].
Block and plot selection were completed in June 2012 followed by
the characterization of initial vegetation cover in all plots in
June 2013 when treatments were initiated. For plant-scale
manipulations, mesquite plants were killed using a direct
application of a clopyralid herbicide (Reclaim; Dow AgroSciences)
in early summer 2013 ca. one month after leaf out (53). This
herbicide was chosen because it does not affect herbaceous plants,
is relatively nontoxic to wildlife (see MSDS D03-106-606), and
readily biodegradable with a half-life of 1.4 days. Effects of
dead mesquite roots on soil moisture and soil organic matter were
temporary because decomposition rates of roots in arid systems are
equal to or faster than mesic systems (54). Herbicide
reapplication was conducted in August 2013 to ensure the mortality
of mesquite plants surrounding the focal patch as well as woody
plants within a larger area were treated (ca. 15m x 15m) because
mesquite plants have an extensive horizontal rooting system (55).
Herbicide spot treatment of mesquite was done annually where
needed; there were few plants with regrowth and all basal; none
after 2017.
For patch-scale manipulations, each focal patch was divided into
1m2 grid cells, and ConMods (Connectivity Modifiers; (23)) were
spatially distributed in grid cells to reduce gap sizes between
live perennial plants to < 0.5 m. The minimum number of ConMods
was added to a patch to meet this criterion. This size is
sufficient to minimize erosional losses across a patch (43).
Reducing bare gap size in a patch below this threshold allows the
accumulation of litter, seeds, and soil to stabilize the soil
surface and improves soil water conditions sufficiently to allow
herbaceous plant establishment within ConMod structures (23).
These fine-mesh (1 cm2), 50 cm wide x 20 cm tall obstructions
reduce horizontal transport of material by wind and water, yet
they do not influence cover of living plants outside of a patch
(22, 51). For plots examining the combined effects of plant- and
patch-scale treatments, mesquite plants were killed in place by
the herbicide and ConMods were located inside the focal patch. For
Control plots, mesquite plants were not treated by herbicide and
ConMods were not placed inside the focal patch. A similar array of
1 m2 contiguous grid cells was located across each focal patch.
Sampling locations were selected to reduce gap sizes between live
plants to <0.5 m.
Plant-scale manipulations (+WPI)
At the plant scale, repeat photos of ten randomly selected ConMods
and ten randomly selected control points were taken annually at
peak growing season (September) from 2013 to 2021. Standardized
overhead photos were used to estimate perennial grass cover as our
response variable, and cover of forbs (annual and perennial),
annual grasses, and litter as fine-scale factors that could
influence grass cover by providing favorable microhabitats for
grass seedlings (litter) or by competing with perennial grass
plants (forbs, annual grasses). Overhead photos were collected
from 2013 to 2021, and image analysis software was used to
classify 100 random points by species, soil or litter (USDA
SamplePoint software; (56)). Species cover in each year was summed
to obtain functional group cover: annual grasses or forbs. For
perennial grasses, we analyzed each species separately and the
total. Rain use efficiency was calculated as the amount of
perennial grass cover produced per unit of rainfall (26).
Patch-scale manipulations (+WPa)
Lateral photos obtained from the same locations as overhead photos
were used to show the redistribution of soil and litter by wind
and water from bare interspaces to ConMods or to herbaceous plant
canopies. Lateral imagery was collected from 2013 to 2017 at
ground level from each cardinal direction to estimate the vertical
accumulation of soil and litter (hereafter "litter").
The area between the tops and sides of each ConMod’s outer rods
and the top contour of the combined litter and soil surface was
determined (methods described in (23) (SimgaScan pro 5.0: Systat
Software, Inc. San Jose, CA USA) using Trace Measurement Mode with
Area and Distance measurement options. The area difference between
years for the same ConMod was a measure of the change in vertical
accumulation of litter or soil. Lateral photos were discontinued
after 2017 when density and cover of herbaceous plants made the
analysis of litter and soil accumulation difficult and prone to
error.
Landscape-scale factors
Factors at the landscape scale were initial shrub cover based on
the 2011 Quickbird image (described above), and annual
precipitation. Daily precipitation data from 13 weather stations
spatially distributed throughout the study area were summed for
the water year (1 October to 31 September) preceding vegetation
and material accumulation measurements, and used as the
precipitation variable. For each block, precipitation from the
nearest weather station was used in the analysis.
Statistical analyses
To examine both the within- and across-block (site) effects of
treatment on perennial grass, RUE, Bouteloua eriopoda, and
Sporobolus flexuosus cover, we constructed two mixed-effects ANOVA
models (nlme R package; (57)). Across-site differences in
treatment effects were tested by including site and treatment as
fixed effects, and time and block as random effects. Within-site
differences in treatment effects were tested using only treatment
as a fixed effect. To determine which treatment differences are
statistically significant, a post-hoc comparison of significant
mixed-effects ANOVA was conducted using the lsmeans R package
(58).
Models of perennial grass cover through time were compared using
simple linear models, logistic models, and split linear models.
Logistic models were fit using the drc package in R (59) and the
pseudo r-squared was calculated as the squared correlation of
fitted values to observed values. For split line regressions, the
breakpoint was determined for each treatment set using the
segmented function of base R. Psuedo r-squared was used to assess
how well each model explains the observed data and, following the
approach described in (27), AIC was used to assess each model’s
predictive ability and characterize perennial grass response to
treatment. AIC support for linear models (lowest value) was
interpreted as evidence for linear internal dynamics (20), support
for logistic models as evidence for logistic growth associated
with a rapidly changing driver (60), and support for split line
models as evidence of alternative attractors (61).
Covariates assessed in this study were examined across three
distinct temporal (initial, lagged 1-year, and current) and
spatial (plant, patch, and landscape) scales. Initial conditions
included: the plant-scale effects of the 2017 perennial grass
cover (%) measured from overhead imagery of ConMods, patch-scale
influences of vertical accumulation of soil and litter (cm2)
estimated from lateral photos, and landscape-scale processes
represented by estimating the percentage of shrub cover in each
block during 2011. Lagged relationships at the plant scale were
represented using the percent cover of annual grasses, forbs, and
litter estimated from overhead imagery in the previous year.
Current conditions at the plant scale were reflected through
annual grass cover, forb cover, and litter cover in the current
year. For landscape-scale factors, we considered total water year
precipitation (cm). Univariate relationships between perennial
grass cover and covariates were assessed by fitting separate
linear regression models for each year and treatment. The
model-based slopes and standard errors of the slopes were then
plotted to illustrate changes in univariate relationships through
time.
To examine which factors were important in predicting perennial
grass cover through time, we compared repeated-measures regression
models that estimated cover for each treatment using the MuMIn
package (62). Competing models were constructed by assembling
unique combinations of 8 explanatory variables with different
characteristic spatial and temporal scales (initial conditions:
plant scale (perennial grass cover in 2017), patch scale (soil and
litter accumulation in 2017), landscape-scale (woody plant cover
2011); current conditions: plant scale (litter) and landscape
scale (precipitation), and time lag variables (t-1): plant scale
(annual grass cover, forb cover, litter cover). For each model,
the additive effects of variables were fixed effects, and
time|block were treated as random effects. Additionally, a null
model with zero slope for each treatment was included in each
suite of competing model to establish a baseline for comparisons.
We used Akiake’s information criterion (AIC) to compare models
estimating grass cover through time, and because this approach
balances parsimony against model fit and complexity, models with
lower AIC can typically be regarded as more reflective of the
underlying internal dynamics (27). For each treatment, the
best-performing models were selected (ΔAIC < 2 indicating a
high degree of empirical support (63)). From those, models with
infrequently occurring variables (occurring in < 30% of subset)
were removed and the model with the lowest AIC was selected. For
each selected model, the estimated relative importance of each
variable was calculated by summing the AIC weights across all
models containing variable x (64). All analyses were conducted in
R (version 4.2.2: R Core Team 2022). We interpreted the strength
of evidence using p ≤ 0.05 (65).
