Linking vegetation patterns to environmental gradients and human impacts in a mediterranean-type island ecosystem

Vegetation patterns at the landscape scale are shaped by myriad processes and historical events, and understanding the relative importance of these processes aids in predicting current and future plant distributions. To quantify the influence of different environmental and anthropogenic patterns on observed vegetation patterns, we used simultaneous autoregressive modeling to analyze data collected by the Carnegie Airborne Observatory over Santa Cruz Island (SCI; California, USA). SCI is a large continental island, and its limited suite of species and well documented land use history allowed us to consider many potential determinants of vegetation patterns, such as topography, substrate, and historical grazing intensity. As a metric of vegetation heterogeneity, we used the normalized difference vegetation index (NDVI) stratified into three vegetation height classes using LiDAR (short, medium, and tall). In the SAR models topography and substrate type were important controls, together explaining 8–15 % of the total variation in NDVI, but historical grazing and spatial autocorrelation were also key components of the models, together explaining 17–21 % of the variation in NDVI. Optimal spatial autocorrelation distances in the short and medium height vegetation models (600–700 m) were similar to the home range sizes of two crucial seed dispersers on the island– the island fox (Urocyon littoralis santacruzae) and the island scrub-jay (Aphelocoma insularis)—suggesting that these animals may be important drivers of the island’s vegetation patterns. This study highlights the importance of dynamic processes like dispersal limitation and disturbance history in determining present-day vegetation patterns.


Introduction
Vegetation patterns are determined by a suite of factors including environmental filtering, disturbance history, dispersal limitation, and biotic interactions (Clements 1916;Pickett et al. 2009). Understanding the relative importance of these processes underpins habitat restoration (Hobbs and Norton 1996), prediction of biotic responses to climate change (Sala et al. 2000), and assisted colonization (Ricciardi and Simberloff 2009). Slow rates of change and unknown site histories make it difficult, however, to attribute pattern to process in many ecosystems. Local-scale studies have shown that environmental filtering (Keddy 1992) is a significant driver of vegetation patterns (e.g. Cornwell andAckerly 2009, Kraft et al. 2008), but relationships between environmental gradients and vegetation in these studies are surprisingly weak given the global-scale correspondence between temperature, moisture, and vegetation type (Whittaker 1975).
A recent focus in ecology has been on the distinction between environmentally-driven ('niche') patterns and stochastic ('neutral') patterns in vegetation composition and biophysical traits (Clark 2009, Hubbell 2001, Rosindell et al. 2012. One implication of this juxtaposition is that many of the processes driving ecosystem assembly are simply unknowable (e.g. undocumented historical events), and therefore are best represented as random, or near-random (Clark 2009). In a rarely disturbed ecosystem made up of long-lived plant species, many pattern-generating processes may indeed be unknowable, yet not all ecosystems meet these criteria. In a more simple system with well-understood biotic and abiotic components, it may be possible to attribute patterns to processes beyond simple correlations with environmental conditions.
Island ecosystems have been used in ecology for decades to address questions of ecosystem assembly (e.g. Simberloff 1976), in part because of their relative simplicity and clear boundaries. Mediterranean-type island ecosystems offer the additional benefit of supporting a mosaic of dominant plant functional types (grasses, shrubs, and trees; deciduous and evergreen), making vegetation patterns readily apparent. Unfortunately, islands in the mediterranean-type climates of the world tend to be located near dense human populations, and so much of the scholarship surrounding mediterranean-type islands has focused on invasion ecology (e.g. Lloret et al. 2005) or the conservation of rare or endangered island endemics (e.g. Morrison et al. 2011). Recent restoration efforts are changing this paradigm, and now globally there are several protected, recovering mediterranean-type islands or parts of islands, including California's Channel Islands (Channel Islands National Park, www.nps.gov/chis) and South Australia's Investigator Group archipelago (Investigator Group Conservation Park, www.environment.sa.gov.au).
Coincident with the recent protection of these island ecosystems, advances in airborne remote sensing now allow for the collection of spatially-explicit, high-resolution information on vegetation structure from light detection and ranging (LiDAR) and function from imaging spectroscopy (Asner et al. 2007). When combined with spatial statistical analyses, these observations can reveal patterns and trends that would not be feasible to document in the field (Dahlin et al. 2013). Tools from spatial statistics are relatively new to ecological studies, and their potential uses are still being explored (Levin 1992;McIntire and Fajardo 2009). Although community assembly studies typically focus on individual species (e.g. Tilman 1985) at relatively small scales (Leibold et al. 2004), the same suite of processes that affect biodiversity can influence patterns of vegetation structure, quality, or chemical composition (McGill et al. 2006), which can be estimated spatially via remote sensing.
Here we address three questions critical to understanding ecosystem processes on Santa Cruz Island (SCI), one of the Channel Islands in California: (1) Which environmental conditions and geographic trends correspond most closely to vegetation patterns on SCI? (2) Is historical grazing intensity a strong predictor of vegetation patterns? And (3) can the variation and spatial patterning that remain once known gradients and impacts have been accounted for be attributed to other processes (e.g. fire history, dispersal limitation, other anthropogenic impacts)?

Study area
SCI is 38 km off the coast of Santa Barbara, California (Fig. 1a) and is the largest and highest of the Channel Islands (250 km 2 , maximum elevation = 753 m). Today, vegetation on SCI consists of grasslands, chaparral, coastal wetlands, sage scrub, oak woodlands, riparian woodlands, Bishop pine (Pinus muricata) forests, and isolated stands of ironwood (Lyonothamnus floribundus ssp. asplenifolius). The island is home to only four terrestrial mammals (excluding humans and bats): the island fox (Urocyon Landscape Ecol littoralis santacruzae), the island harvest mouse (Reithrodontomys megalotis santacruzae), the SCI deer mouse (Peromyscus maniculatus santacruzae), and the island spotted skunk (Spilogale gracilis amphiala). Many birds occupy the island, but most relevant to this study is the island scrub-jay (Aphelocoma insularis). The island fox and island scrub-jay have recently been the focus of major conservation efforts (Coonan et al. 2010;Morrison et al. 2011).
Human occupation of the Channel Islands began with Native Americans, with the oldest artifacts dated [11,000 YBP (Rick et al. 2001). The introduction of livestock did not take place until 1830 (Ord 1956), but by the mid-1800s sheep, cattle, horses, goats, and pigs had all been introduced. Many of the livestock populations, especially sheep, expanded beyond what could be managed and many parts of the island were overgrazed. In the 1960s and 1970s an estimated 180,000 sheep were dispatched from the island (Junak et al. 1995). Since 1978 The Nature Conservancy (TNC) has managed the western portion of SCI; the US National Park Service (NPS) acquired the remainder of the island in 1997. Conservation management of the island has focused on the removal of non-native species, such as sheep, cattle, and pigs. Today all introduced mainland vertebrate species have been eradicated and vegetation has been recovering for the last few decades (Morrison et al. 2011). While much of the native vegetation has survived in refugia or returned with the removal of livestock, parts of the island have converted to dense stands of fennel (Foeniculum vulgare), a common California invasive Landscape Ecol plant, many of the grasslands are dominated by European annuals (e.g. Avena fatua), and some areas remain unvegetated. Since a large fire in 1931 (4850 ha;Junak et al. 1995), only a few small burns have occurred.

Airborne remote sensing
We collected high-resolution data over the entirety of SCI in August 2007 using the Carnegie Airborne Observatory Beta system (CAO Beta; decommissioned in 2009; http://cao.ciw.edu). CAO Beta combined three instrument sub-systems into a single airborne package: (1) the Jet Propulsion Laboratory's airborne visible/infrared imaging spectrometer (AV-IRIS; Green et al. 1998), (2) a small-footprint waveform LiDAR scanner, and (3) a global positioning system-inertial measurement unit (GPS-IMU) (Asner et al. 2007). For this study, the system was operated at an altitude averaging 3.4 km, providing spectroscopic measurements at 2.2 m spatial resolution. The LiDAR sub-system was configured for discrete-return laser point recording with a laser spot spacing of 1.5 m. In this semi-arid environment, the August timing of the flight meant that most annual grasses and forbs had senesced, while woody plants, most of which are evergreen, still had their leaves.
The AVIRIS data were converted to at-sensor radiances by applying radiometric corrections developed during sensor calibration in the laboratory. Apparent surface reflectance was derived from the radiance data using an automated atmospheric correction model, ACORN 5LiBatch (Imspec, Palmdale, California, USA). Atmospheric correction algorithm inputs included ground elevation (from LiDAR), aircraft altitude (from the GPS-IMU), viewing and solar geometry, atmosphere type (temperate), and estimated visibility in kilometers. The code uses a MODTRAN look-up table to correct for Rayleigh scattering and aerosols. Water vapor was estimated from the 940 and 1,140 nm water vapor features in the radiance data (Asner and Heidebrecht 2003).
To quantify variation in vegetation on a continuous scale, we calculated the normalized difference vegetation index (NDVI; Tucker 1979) across the island. NDVI was calculated using AVIRIS bands 49 (830 nm) and 33 (674 nm). NDVI is known to correspond to species differences in structure and chemical content (e.g. chlorophyll and foliar nitrogen content) in California vegetation (Gamon et al. 1995). Due to the late summer flight time, here high NDVI values ([0.6) primarily correspond to contiguous areas of large evergreen trees and shrubs with high leaf area indices, while low NDVI values (0.12-0.2) likely correspond to senesced grasses, and intermediate values correspond to subshrubs, shrubs, and pixels that may contain a mixture of shrubs and grasses. Based on our knowledge of the island, we determined that areas with NDVI less than 0.12 were not vegetated (roads, river beds, and unvegetated slopes, Supplemental Fig. S1).
The LiDAR sensor allows for high-resolution measurements of vegetation height (Asner et al. 2007), but due to the small stature of much of the vegetation on SCI, LiDAR point-cloud optimization was adjusted to classify relatively small changes in 'first returns' as vegetation, not ground. As a result, some steep cliffs were misclassified as tall vegetation. We used a two-step process to remove these areas from further analysis. First, large areas with NDVI less than 0.12 were removed from further analysis. Then, remaining pixels with NDVI less than 0.3 but with a vegetation height measured as greater than 5 m (unlikely to be living trees) were masked so that NDVI and vegetation height equaled zero. Pixels with a local ground slope (3 pixel kernel) greater than 50°w ere also removed from further analysis, as were areas containing large human-built structures or large stands of non-native trees. In order to reduce the size of the data set and to reduce very fine-scale variation, the remaining NDVI, vegetation height, and ground digital elevation model (DEM) pixels were then aggregated to 25 m resolution, resulting in a data set of 300,772 pixels.

Data distribution
The overall distribution of NDVI values appeared to be multi-modal (Fig. 2). In lieu of transforming these data, we stratified the data by LiDAR-derived height (\0.25 m, 0.25-1.5 m, and [1.5 m), producing nearnormal distributions in NDVI (Fig. 2), and subsequent analyses were performed separately on three heightstratified NDVI distributions. Using a detailed vegetation classification for SCI (The Nature Conservancy 2007), we found qualitatively that the three height classes correspond to different plant communities (Fig. 1b, c). The short height class (\0.25 m) includes Landscape Ecol annual and perennial grasses and patches subshrubs like of buckwheat (Eriogonum spp.) and California sagebrush (Artemisia californica). The medium height class (0.25-1.5 m) is comprised of most of the shrub species, including lemonadeberry (Rhus integrifolia), several manzanita species (Arctostaphylos spp.), and the island scrub oak (Quercus pacifica). The tall height class ([1.5 m) includes some large Q. pacifica, as well as several tree species such as Quercus tomentella, L. floribundus, and P. muricata. For the rest of this study we will discuss variation in NDVI within these three height classes. Variation in NDVI within a height class therefore does not correspond to changes in plant structural types (grass versus shrubs versus trees), but within a height class variation in NDVI may correspond to species differences as well as to intraspecific heterogeneity caused by environmental conditions, stressors, or other factors. While NDVI is a coarse metric of vegetation change, it has been used for decades to assess vegetation heterogeneity at a variety of scales (e.g. Townshend and Justice 1986;Gamon et al. 1995;Muratet et al. 2013).

Environmental variables
We mapped 43 potential environmental predictor variables, which can be divided into topographic and substrate variables (for a full list of variables see electronic supplementary Table S1, a subset of maps is shown in Fig. 3, frequency distributions for all nonbinary variables are in Fig. S2). Topographic variables are those calculated from the LiDAR-derived DEM and include slope, aspect, soil wetness index, curvature metrics, and incident solar radiation (insolation). All variables were generated using IDL ENVI v. 4.8 (ITT, Boulder, Colorado, USA), ArcGIS v. 9.3 (ESRI, Redlands, California, USA) or a combination of the two. Slope, aspect, and the curvature metrics were each calculated using a local kernel (3 9 3 pixel area, *0.56 ha) and a general kernel (13 9 13 pixel area, *10.56 ha). Aspect was converted to a linear variable, ''northness' ' [cosine(aspect)]. Soil wetness index (SWI) was calculated as where A S is the specific catchment area (area upslope of the pixel), and b is the slope at the pixel (k = 3) (Moore et al. 1991). The six curvature measuresprofile convexity, plan convexity, longitudinal convexity, cross-sectional convexity, minimum curvature, and maximum curvature-together orthogonally describe the geomorphology of an area (Wood 1996). The root mean squared error (RMSE) variable measures the variation of the data within a given kernel, an approximation of the surface roughness within a pixel. Insolation was calculated for the summer and winter solstices and the equinoxes in ArcGIS (Rich et al. 1994). Substrate variables were derived from a scanned and orthorectified geology map of SCI provided by the Pacific Section of the Society for Sedimentary Geology (Weaver and Nolf 1965). Substrate type and subtype polygons were hand digitized at 1:15,000 scale, then aggregated into 20 geologic types and converted into binary categorical variables (Fig. 3e).
Land use history/grazing intensity A map of historical grazing intensity was digitized by Perroy (2009), modified from Schuyler (1993) and originally developed from Santa Cruz Island Company records (Fig. 3f). Based on these historical records, the island was divided into 24 polygons with variable grazing intensity (low, medium, or high). As with substrate, these grazing polygons were converted into categorical variables. While some information was available about the relative grazing intensity of the Landscape Ecol zones, we treated each zone as a separate binary variable. Many other agricultural practices occurred on the island, including viticulture and growing hay for livestock (Junak et al. 1995), however spatially explicit information about the distributions of these activities was not available. Fig. 3 A selection of the variables included in the predictive models. a elevation (ranges from sea level to 753 m), b northness (kernel (k) = 3 9 3 pixel area; ranges from -1 to 1), c minimum curvature (k = 3; -39.42 to 12.78), d maximum curvature (k = 13; -1.37 to 5.55), e substrate (see Appendix  Table S1 for full list of substrate types), f grazing intensity, g (scaled X) 2 (0 to 5.55), h (scaled X) 2 *(scaled Y) (-2.93 to 7.85) Landscape Ecol Geographic trends Santa Cruz Island's size, orientation, and coastal location allow for different weather and climate patterns to occur across the island. To test the effects of broad geographic trends, we included the east-west (x) direction and north-south (y) direction, using the scaled central coordinates (eastings and northings) for each pixel (mean = 0, s.d. = 0.5; Gelman 2008), as well as third-degree polynomial terms (Legendre and Legendre 1998). The resulting geographic trend portion of the model took the form: where the b 1-9 values are the estimated model coefficients. These polynomial terms modify the centered and scaled eastings and northings so that, for example, X 2 is a surface where the eastern and western ends of the island are numerically similar. Fig. 3g, h are examples of these geographic trends, which could correspond to variation in weather patterns across the island, temperature gradients, or other unmapped environmental gradients.

Statistical modeling
To understand the separate and combined effects of different environmental conditions and geographic trends, along with the influence of spatial autocorrelation, we took a multi-step approach. All statistical analyses were performed in R (R Development Core Team 2011). First, we combined all data sets, resulting in a stack of 76 possible explanatory variables, plus vegetation height and NDVI. To reduce processing time in statistical analyses we took an even 5 % subsample of the 300,772 pixels, leaving 15,039 pixels across the island. Because the island is not a rectangle, and since many pixels were masked due to clouds, lack of vegetation, or human impacts, this did not result in a perfectly uniform sampling across the island, but sampling was still very dense, with an average minimum distance between subsampled pixels of 68 m (s.d. = 31 m). The subsampled data were then divided into the three height classes: short (\0.25 m, n = 8,783), medium height (0.25-1.5 m, n = 4,241) and tall ([1.5 m, n = 2,015). Average distances between points in different height classes were slightly longer than the overall average distance, with an average minimum distance between points of 81 m for the short height class (s.d. = 45 m), 108 m (s.d. = 74 m) in the medium height class, and 113 m (s.d. = 115 m) in the tall height class. The large standard deviation in the tall height class is due to a right-skewed distribution of minimum distances, with 29 points (1.4 % of the data in the tall height class) located more than 500 m from a neighbor in the same height class. After subsetting and splitting the data set, all non-binary data were scaled (mean = 0, s.d. = 0.5; Gelman 2008). For each of the three NDVI data sets (short, medium, and tall), to minimize multicollinearity, if two of the independent variables were found to have an absolute correlation coefficient greater than 0.5, the independent variable that was less correlated to NDVI was removed from further analysis. With the reduced set of independent variables, ordinary least squares regression (OLS) was used to develop an initial model, and models were further reduced using a reverse stepwise algorithm based on the Akaike Information Criterion (AIC: Akaike 1974; stepAIC: Venables and Ripley 2002). This algorithm calculates the AIC for a model with all possible variables, then re-calculates the model with each variable removed, choosing the model with the largest drop in AIC. This continues until the model AIC no longer decreases. Since this is a somewhat permissive algorithm (Venables and Ripley 2002), in some cases the OLS model still contained some non-significant predictors (p [ 0.05) and which were removed in the final OLS model.
We used Moran's I computed over all directions (isotropic) to test for spatial autocorrelation (Haining 2003). In all three height classes we found significant spatial autocorrelation in the OLS model residuals (p ( 0.001), and so we computed isotropic correlograms for each dataset. To account for this spatial autocorrelation in the regression models we used simultaneous autoregression (SAR; Haining 2003). Recent reviews have shown this approach to be successful when used on simulated ecological data, in comparison to other methods (Dormann 2007;Kissling and Carl 2008;Beale et al. 2010).
From the correlograms of the OLS models' residuals we estimated the distance where Moran's I was near zero and used these as the maximum neighborhoods in the SAR models with a spatial error model (spdep: Bivand et al. 2012). Our original neighborhood estimates still produced models with spatially autocorrelated residuals and so we iteratively tested distances and distance functions in a range around the initial estimates until the global Moran's I was apparently minimized. Ideally this process would be automated to test every possible distance function and distance within some threshold distance of the initial estimate, but since each model took several hours to run, this process was performed manually. We then removed any nonsignificant predictor variables from the model (p \ 0.05), recalculated, tested for spatial autocorrelation, and produced a model that included all significant environmental, geographic, and land-use variables and the final neighborhood distance that minimized spatial autocorrelation of the residuals.
To compare the explanatory power of different types of factors, we used partial regression analysis (Legendre and Legendre 1998;Lichstein et al. 2002;Dahlin et al. 2012). We separated out sets of the remaining significant variables to calculate the relative importance of environmental conditions, grazing intensity, geographic trends, and spatial autocorrelation in explaining variation in NDVI in the three height classes.

Environmental condition correlations
The NDVI models all began with a total of 76 possible explanatory variables. The strongest individual correlations were with variables relating to incident solar radiation (insolation and northness) and substrate type. Table 1 shows the 20 strongest relationships, and all of the correlations are reported in electronic supplementary Table S1. Northness has a positive relationship with NDVI (north-facing slopes are greener) and insolation correlations are negative (sunnier places have less green vegetation). Several substrate types also had strong correlations-for example, the SCI Volcanics had a negative relationship with NDVI in all three height classes, while the Monterey Formation had a positive relationship with all three. In nearly all cases grazing polygons classified as having high grazing intensity had a negative correlation with NDVI.

OLS-SAR modeling
From the original set of 76 explanatory variables, removal based on multicollinearity resulted in a total of 47 variables being included in the short NDVI model, 51 in the medium height NDVI model, and 46 in the tall NDVI model. Following the AIC-based model selection procedure, the models contained 34, 30, and 19 variables, respectively. The short NDVI OLS model explained the most variation, with a coefficient of determination of 0.36. The medium NDVI model explained less (r 2 = 0.30), while the tall NDVI model explained the least variation (r 2 = 0.22). All three models had significant spatial autocorrelation in the The final SAR models improved explanatory power over the OLS models by an average of 8 %. Including spatial autocorrelation also decreased the number of significant predictor variables in all three models-the final SAR models contained 28, 18, and 14 predictor variables, respectively. The short NDVI SAR model explained the most variation (r 2 = 0.48, an increase of 12 % from the OLS model), the medium height NDVI model explained 38 % of the variation (an increase of 8 %), and the tall NDVI model explained 27 % of the variation (an increase of 5 %). Figure 4 shows the coefficient values for all of the variables that appeared in two or more of the final models (all values reported in electronic supplementary Table S1). High elevation was associated with lower NDVI values in the short and medium NDVI classes, but was not a significant predictor of tall NDVI. Maximum curvature (ridges) and surface RMSE were associated with lower NDVI values across all three height classes. Sunnier areas were associated with less green vegetation in the short height class, while NDVI in the tall class had a slight positive relationship with summer insolation. The strongest substrate predictors that appeared in multiple models were the Quaternary alluvium (GEO1; positive), Blanca Volcaniclastics group (GEO9; negative), and the Monterey formation (GEO14; positive). Broad geographic trends were minor components of all of the models, mostly in the north-south direction (Y 2 and Y 3 ), but east-west in the short height class (X 2 Y). Many of the grazing intensity zones were important predictors, with low and mid-level grazing intensity always associated with higher NDVI values, while high intensity grazing led to lower NDVI values in the medium and tall height classes.
We used correlogram range distance estimates to approximate neighborhood distances for the SAR models, then adjusted them as necessary to minimize Moran's I in the residuals. We weighted the distance matrices using a diminishing function of 1/distance. In addition, we tested an even distance function (all distance weights = 1) but found that these models had reduced performance compared to the diminishing distance function (data not shown). The tall NDVI model had the longest optimal spatial autocorrelation distance, with a 1,200 m radius (*452 ha), while the short and medium height models had similar distances of 600 m (113 ha) and 700 m (154 ha), respectively (Fig. 5). In contrast, the vegetation classification map (TNC 2007) delineated much smaller vegetation patches, averaging 2.3 ha in size (s.d. = 8.8), with only five of the 10,649 classified vegetation patches having an area greater than 100 ha.
To compare the different model components we calculated OLS coefficients of determination with different sets of variables removed (Fig. 6). When we compared explained variation between the models, 25 % of the variation in the short vegetation model Fig. 4 Simultaneous autoregressive model coefficient values for the 14 variables that were significant in at least two of the models. LUH land use history, low, mid, and high refer to relative grazing intensity, GEO Substrate. See Appendix Table S1 for all coefficient values was linked to environmental conditions, whereas measured environmental conditions explained 18 % of the variation in the medium height model, and only 8 % of the variation in the tall model. In all three models spatial autocorrelation explained 7.5-13 % of the total variation, while historical grazing intensity explained 7.7-10 % of the variation across the three models.

Discussion
We used airborne remotely sensed data and spatial statistics to show how patterns in vegetation greenness (NDVI) vary with environmental conditions, known land use history (historical grazing), geographic trends, and spatial autocorrelation. This analysis allowed us to demonstrate that environmental conditions like topography and substrate type do play a role in plant patterns on SCI, but that land use history and spatial patterning were also important drivers of current vegetation patterns.

The importance of environmental conditions
Of the three NDVI models considered in this study, patterns in the short vegetation were most closely linked to environmental conditions. Given the August flight time (near the end of the dry season in California), this is likely a reflection of environmental B neighborhood size for medium height NDVI model. C average home range of island fox (Roemer et al. 2001). D maximum documented foraging distance of island scrub-jay (Atwood 1980). E neighborhood distance in tall NDVI model.  Table S2 Landscape Ecol controls on the variability of grass and forb senescence, where some herbaceous plants had senesced completely, while others, perhaps in more shaded areas or with more access to water, were still green. As expected, high insolation was associated with lower NDVI values and alluvial areas (Quaternary alluvium, GEO1) were associated with higher NDVI values. Variation in greenness of the shrubs and subshrubs of the medium height class was most closely tied to local (kernel (k) = 3) slope and northness, with steep, south-facing slopes having lower NDVI values in this height class than other areas. Much of the dense shrub vegetation on SCI is concentrated on the interior north-facing slope, which is also a different substrate, SCI schist (GEO20), a comparatively high-nutrient substrate on the island. A photo from 1993 in Junak et al. (1995, p. 20) shows dense woody vegetation in this area which suggests that this plant community persisted through times of heavy grazing pressure. This persistence could be due to a positive feedback whereby the substrate type and aspect allowed dense woody vegetation to thrive, and then the combination of steep slopes and thick vegetation prevented overgrazing, allowing vegetation to persist, thus further stabilizing the soil, allowing plants to grow larger, and so on.
NDVI patterns in the tallest vegetation on SCI were the least well explained by environmental conditions, with substrate types being the strongest environmental drivers. The Monterey formation (GEO14) had higher than otherwise expected NDVI (e.g. a positive coefficient value in the SAR model), while areas of Blanca volcaniclastics (GEO9) had lower than expected NDVI values in the tall height class. While the Blanca substrate may inherently support more sparse vegetation and is known to be landslide-prone, there are some patches of oaks and pines growing within it, suggesting another explanation: it is located on steep, south-facing slopes in one of the most heavily grazed areas of the island (Junak et al. 1995;Perroy 2009), so this association is likely another positive feedback where soil erosion due to intense grazing is exacerbated by substrate instability and aspect.

Land use history
Grazing zones were important predictors of NDVI-in all cases areas with intense grazing had lower than otherwise expected NDVI values, while low and medium intensity grazing led to higher NDVI values. In particular grazing zones 23 and 24 (Fig. 3f) were significant predictors of low NDVI for shrub and tree vegetation, where sheep were removed approximately 10 years after they were dispatched from the rest of the island. NDVI patterns in the tallest vegetation on SCI had the strongest link to grazing intensity. This result was largely driven by a few low-and mid-intensity grazing areas (LUH3, 5, and 12), where tall vegetation appears to be more dense and green than in other areas. This result corresponds well with the conclusions of Perroy (2009) who showed that mapped grazing intensity was the strongest determinant of post-grazing recovery across SCI. Actual grazing intensity was undoubtedly more heterogeneous than the map suggests (Fig. 3f), but this gives us a coarse approximation of the major human-driven impact on the island. While the finding that high grazing intensity led to less green vegetation is not surprising, it is noteworthy that our modeling approach was able to identify which areas had high grazing impacts without this information being provided a priori.

Geographic trends
We expected to find a significant east-west trend due to the movement of fog to and from the California coast, but in the short NDVI model the importance of the scaled X 2 Y geographic gradient (Fig. 3h) suggests that the two ends of the island are actually similar in a way that is not captured in the other environmental conditions. This unexpected trend is possibly a signal from past land-use: both of these areas were historically hayed (Symes andAssociates 1922 cited in Junak et al. 1995) and are now dominated by European annual grasses (The Nature Conservancy 2007). This similarity could also be due to persistent weather patterns acting on the extreme ends of the island, or due to some other unmapped environmental condition or historical event.
The medium and tall NDVI models had negative relationships with trends in the north-south direction (Y 2 and Y 3 , respectively). In the case of the medium height NDVI model this suggests that the center of the island's shrubs are greener than otherwise expected (i.e. if this variable were not included in the model), while among tall vegetation the southern part of the island is greener than otherwise expected.

Spatial autocorrelation
All of the models included a substantial spatial autocorrelation component, although interpreting these patterns is challenging. Neighborhood distances in spatial autocorrelation analysis are not necessarily equivalent to patches, nor do they have clear causal relationships to processes (Legendre and Legendre 1998). We can, however, look for possible explanations for these patterns through reviews of the literature and other sources of information.
The medium height class is largely comprised of fruit and acorn bearing shrubs, and the short height class may be in the process of being colonized by these woody plants, making these two height classes likely candidates for dispersal limitation by animals. The relatively short neighborhood distances for the short and medium height NDVI models (600 and 700 m, respectively, Fig. 5a, b) suggest that dispersal limitation may be a key driver of vegetation patterns on the island. Two of the main dispersers of woody plant seeds are the island fox and the island scrub-jay. The average home range of the island fox has been reported as 55 ha (or a circle with radius 418 m; Fig. 5c), but it varies from 15 to 87 ha (219-526 m radius; Roemer et al. 2001). Berries are a common part of its diet (Crooks and Van Vuren 1995), and like other fox species, it marks its territory with scat, making the island fox an important dispersal agent for fruiting plants. Island scrub-jays have been observed travelling up to 650 m (Fig. 5d) to obtain acorns (Atwood 1980). These acorns are then scatter-cached, often buried in bare patches of soil. While many birds visit SCI, few are permanent residents or large enough to consume acorns, though many likely consume berries (Collins 2011). Given the large size of the island and therefore range of potential scales of spatial autocorrelation, the correspondence between the neighborhood distances that emerged in the short and medium height models and these two animals' home ranges suggest that 13 and 11 % of the variation in the short and medium height NDVI models, respectively, could be related to dispersal limitation. This correspondence is notable because the link between spatial autocorrelation and dispersal limitation is often suggested (e.g. Beale et al. 2010), but it is very difficult to quantify at the landscape scale. Other possible explanations for these scales of spatial autocorrelation could relate to past land use or fire history, though we were unable to find records of other possible impacts with a similar spatial scale.
In contrast, the shrubs and trees that comprise the tall height class appear to be survivors. The majority of these plants likely established during the era before the 1990s when the island was heavily occupied by sheep, cattle, and pigs. While there are no clear links between the neighborhood distance (1200 m, *450 ha; Fig. 5e) and known processes, this scale is in the same realm as the approximate home ranges of the island's now eradicated feral pigs (290 ± 125 ha; unpublished data, cited in Parkes et al. 2010) and at the lower end of the range of published feral pig home range sizes globally (Saunders and McLeod 1999). Feral pigs can have major impacts on ecosystem assembly, both through seed predation and through rooting (Parkes et al. 2010), which could potentially have an impact on the quality and structure of the vegetation in this height class. Other possible explanations for this scale of spatial autocorrelation could include patchy fog moving across the island, historical fires, or human impacts not documented in the grazing intensity map.

Sources of uncertainty
Overall the SAR models were able to explain less than half of the variation in greenness within the three height classes. We expect that the remaining variation is due to a combination of factors. Though our substrate map was highly detailed, it still undoubtedly left out subtle variations in geology and outcrops. We elected not to include a soils map in this analysis as soils maps are often derived from aerial photographs and are therefore based on vegetation cover, however differences in soil type likely do influence vegetation patterns. In addition, we have no information about the spatial patterns of past fires across SCI or human impacts other than grazing intensity. Given that SCI appears to be in a state of post-disturbance recovery, it is likely that much of the variation is due to heterogeneous recovery-small differences in environmental gradients, where seeds germinate, and how fast they grow will add to the uncertainty in this analysis.
It is also possible that our methods, including subsampling, added to the uncertainty in our models. As computing power becomes more readily available we hope to test these methods on larger, denser data sets. While NDVI appeared to effectively capture the heterogeneity across the system, and has been shown to be a good measure of diversity in California ecosystems (Gamon et al. 1995), it is also possible that, especially in the tall height class, NDVI was saturated, inhibiting our ability to map variation. However, as Fig. 2 shows, the vast majority of the vegetation on SCI falls below an NDVI of 0.6, where relationships between NDVI and vegetation properties tend to break down (Gamon et al. 1995).

Conservation implications
Despite the remaining uncertainty in the models, understanding the relative importance of the different factors driving ecosystem assembly can have direct impacts on conservation management decisions. In California's Channel Islands, one proposed strategy to reduce the extinction risk of the island scrub-jay is to establish a population of jays on neighboring Santa Rosa Island (SRI; Morrison et al. 2011). In 2011 the last of the introduced herbivores were dispatched from SRI, and it has been suggested that the introduction of scrub-jays could hasten re-establishment of native oaks and pines across that island (Morrison et al. 2011). Our study suggests that scrub-jays likely do act as 'ecosystem engineers' (Byers et al. 2006) in these systems by transporting acorns long distances. We also show, however, that other processes are at work as the island's vegetation recovers, many of which are poorly understood. In planning restoration projects, particularly in highly eroded areas, other more intensive approaches may be necessary, including starting seeds in a favorable environment, watering outplantings, inoculating soil with appropriate mycorrhizal fungi, or preventing herbivory. These actions may be critical to plant survivorship, especially in severely disturbed, semi-arid ecosystems. This study highlights areas that could benefit from more intensive interventions, and shows that some relatively intact areas of vegetation do still remain on SCI.

Conclusions
Our study showed that environmental conditions like topography and substrate type do play a role in determining the distribution of vegetation on SCI. We also found that heterogeneity in grazing intensity was a significant predictor of vegetation patterns, and interestingly, we found that our relatively simple approach allowed us to distinguish between high and low intensity grazing areas. We also found that geographic trends were significant predictors, and in particular trends that showed the east and west ends of the island to be similar, despite different substrates, possibly indicating similar outcomes from past land use. Finally, we found that the scale of spatial autocorrelation emerging from the models corresponded closely to the home range sizes of animals living on the island. While this correspondence does not demonstrate a direct link, it does suggest that seed dispersal by animals may be a significant factor determining the current vegetation patterns across the island.