Development and validation of a habitat suitability model for the non-indigenous seagrass Zostera japonica in North America

We developed a spatially-explicit, habitat suitability model that can be used to identify and predict areas at higher risk for non-native dwarf eelgrass (Zostera japonica) invasion. The niche-based model uses simple readily available environmental parameters (depth, near shore slope, and salinity) to quantitatively describe habitat suitable for Z. japonica invasion based on ecology and physiology from the primary literature. Habitat suitability is defined with values ranging from zero to one, where one denotes areas most conducive to Z. japonica and zero denotes areas not likely to support Z. japonica growth. Functional forms and equations for the ZJHSM were developed a priori, and the model was validated by comparison with multiple years of independent field-collected spatial Z. japonica maps from Yaquina Bay, Oregon, USA, an area that has well documented Z. japonica expansion over the last two decades. Sensitivity analysis performed to evaluate the contribution of each parameter to the model prediction revealed that depth was the most important factor. The highest suitability values for Z. japonica occurred in the mid to upper portions of the intertidal zone, with larger expanses occurring in the lower estuary. While the upper estuary did contain suitable habitat, most areas were not as large as in the lower estuary, due to inappropriate depth, a steeply sloping intertidal zone, and lower salinity. The lowest suitability values occurred below the lower intertidal zone, within the Yaquina River channel. Shallow sloping intertidal soft-bottom sediments appear to provide optimum habitat for this species. This model would allow resource managers to identify those areas at risk for future invasion, and develop proactive plans for limiting the opportunities for further introduction and spread in areas where it does not currently occur.


Background
Two of the biggest alterations of estuarine and marine biodiversity and ecosystem function are likely to be driven by climate change (e.g.temperature and pH) and anthropogenic species introductions.It is likely that the impacts of these separate processes may work in a non-linear manner, thus climate change may exacerbate the spread of non-native species resulting in alterations of ecosystem services (Bierwagen et al. 2008;Hershner and Havens 2008).Consequently, there is intense interest in modeling the distribution of non-native species with respect to their ecological and physiological tolerances.Descriptive and mechanistic modeling approaches have been used to predict distributions.The descriptive approach uses habitat suitability models which relate species occurrences to a suite of environmental variables.Habitat suitability models are based on ecological niche theory and describe the relationships between species occurrence or abundance and a set of environmental data (Guisan and Thuiller 2005).This type of model can be used as a screening or risk assessment tool by resource management agencies in order to efficiently utilize limited resources.The mechanistic model approach focuses on predicting invasion dynamics or ecosystem response (Gallien et al. 2010) via empirically derived relationships.Mechanistic models may be inherently more satisfying from an ecological perspective but they tend to be site specific and data intensive.Often the physical and biological data required to drive these models are not readily available and expensive to collect.Consequently, the large scale application of numerical mechanistic ecosystem response models for invasion ecology is generally more limited.Previous efforts to develop a spatially explicit numerical population model (Almasi and Eldridge 2008) of Zostera japonica Ascherson and Graebner colonization was plagued by limited data availability for model development and validation.Although the Almasi and Eldridge (2008) model output suggests that disturbance, competition and survival were important variables, most of these parameters were not well constrained by empirical or physiological data, and the model was subjected to limited validation.Consequently, the predictive capacity of the Almasi and Eldridge model has not been verified over interannual time scales and as a result has not been adopted for resource management.
Niche-based habitat suitability models have been widely used as management screening and risk assessment tools to predict the potential distributions of invasive species over large spatial scales (Thuiller et al. 2005;Inglis et al. 2006;Gallien et al. 2010;Jimenez-Valverde et al. 2011).Here we present a niche-based habitat suitability model of the non-native seagrass Zostera japonica that is expanding its range along the Pacific coast of North America (Shafer et al. 2014).
Zostera japonica is believed to have been unintentionally introduced with aquaculture products to North America, and is currently distributed in bays and estuaries in British Columbia (Canada), and Washington, Oregon and northern California (USA).A recent review on the science and management of Z. japonica in North America (Shafer et al. 2014) noted that invasion of new areas is continuing at a rapid pace.Although several studies have documented Z. japonica expansion (Gaeckle et al. 2001;Young et al. 2008;Young et al. 2015), the underlying factors responsible for the rapid increase in areal extent have not been identified.Moreover, several studies suggest that this species has only colonized a small fraction of the available suitable habitat throughout its potential range in North America (Harrison and Bigley 1982;Shafer et al. 2008;Shafer et al. 2011).Recent research has focused on the physiological ecology of Z. japonica.For example, it has been determined that Z. japonica is eurythermal with optimal growth at 20 °C (Shafer et al. 2008) and a lethal thermal threshold at 35 °C (Kaldy and Shafer 2012).North American Z. japonica populations have optimal growth at salinity 20 but are euryhaline, and can tolerate daily freshwater immersion (Kaldy 2006;Shafer et al. 2011).Recent work suggests that vertical zonation patterns can be controlled by temperature, not light limitation (Shafer and Kaldy 2014;Kaldy et al. 2015).The potential for interaction with the native seagrass Z. marina has also been explored.Some work suggests negative interactions (e.g.Bando 2006), but more recent work suggests that the scale of competitive interactions is limited to specific microtopographic features (Hannam and Wyllie-Echeverria 2015).Although some aspects of plant ecophysiology are better understood, the ecological consequences (e.g.ecosystem metabolism) and economic impact of Z. japonica colonization are poorly understood (Shafer et al. 2014).The transition from bare mudflat to a nonnative seagrass habitat alters the physical and biogeochemical process (e.g.net primary productivity, nutrient cycling, etc.) associated with the habitat; additionally negative economic impacts to intertidal Manila Clam aquaculture have been documented (Patten 2014).Consequently, understanding the potential for continued colonization of unvegetated mudflat in North America is important from both an ecosystem conservation and economic perspective.
Management of Z. japonica in North America varies between States and Countries and has changed over time (see Shafer et al. 2014 for review).Recent management actions in the USA have been predominantly geared toward control or eradication but have generally been unsuccessful at the estuarine scale (Shafer et al. 2014).One aspect that has hampered effective management of this non-native species has been a lack of predictive tools to identify potential areas at risk for invasion which has suppressed the ability to develop early detection and monitoring protocols.Development of a model to identify areas at risk of Z. japonica colonization would allow resource agencies to focus limited monitoring funds on areas susceptible to colonization, thereby utilizing those resources in a more efficient manner.Rapid discovery of new populations of nonnative species greatly increases the probability of eradication or elimination (Myers et al. 2000;Simberloff 2003).
Our objective was to develop and test a simple niche-based habitat suitability model for Zostera japonica as a risk assessment tool to predict areas expected to be susceptible to colonization.The model was specifically developed to utilize three relatively simple environmental variables (salinity, depth and bathymetric slope) that are widely available or can be readily calculated for many estuaries.The simple data requirements allow application of the model to areas where there are insufficient data available to develop more complex mechanistic models.

Model Overview
The Zostera japonica Habitat Suitability Model (ZJHSM) model is a spatially-explicit, grid-based model that was developed to predict suitable locations for the invasion and growth of dwarf eelgrass (Z.japonica).The terminology and model evaluation techniques were adapted from Pollack et al. (2012).The model assumes that the three driving variables (depth, bathymetric slope and salinity) quantitatively describe suitable habitat for Z. japonica invasion.Each variables is assigned a dimensionless suitability index value that represents the relationship between the environ-mental variable and the habitat requirements of Z. japonica (noted as ZSI n ).Each ZSI is represented quantitatively as a series of linear suitability curves, with a minimum value of 0 for unsuitable to 1.0 for optimal habitats.The shape of each suitability curve was based on the expert opinions of the co-authors, reviews of the relevant peer-reviewed literature, and more than a decade of research on the physiological tolerances of Z. japonica, collected from multiple locations within its introduced range in North America (Kaldy 2006;Shafer 2007;Kaldy and Shafer 2012;Shafer and Kaldy 2014;Shafer et al. 2008Shafer et al. , 2011Shafer et al. , 2014;;Kaldy et al. 2015aKaldy et al. , 2015b)).Suitability curves are formulated as stepfunctions with linear approxi-mations between each step.An overall suitability index is calculated as the geometric mean of the ZSI values and represents the overall suitability for invasion (noted as HSI).Data and equations are imported into a GIS and applied to specific georeferenced locations within Yaquina Bay, Oregon, where the distribution of Z. japonica has been well studied.
A portion of the Yaquina River, located in Oregon between its confluence with the Pacific Ocean and river mile 14 south of Toledo, was selected as the model domain area (Figure 1).This area was chosen due to the documented expansion of Z. japonica (Young et al. 2008) and the large number of habitat surveys and related environmental data available to build and evaluate the model hypothesis (Young et al. 2015).

Suitability Indices
Depth (ZSI Depth ) Depth is one of the most important factors affecting seagrass distribution worldwide, due to their dependence on light to drive photosynthetic processes (Dennison 1987;Duarte 1991;Koch 2001).In its native range, Z. japonica has been reported to grow as deep as 3-7 m, although it typically grows at depths < 1 m (Hayashida 2000;Nakaoka and Aioi 2001;Abe et al. 2010).Unfortunately, very few (if any) of the publications on Z. japonica from Asia report the tidal datum; consequently the actual depth distribution in Asia is somewhat unclear.Throughout its established range in North America, Z. japonica is found primarily in mid-to upper-intertidal zones (Shafer et al. 2014;Young et al. 2015).Although there are currently no documented cases of Z. japonica growing below mean sea level in North America (Shafer et al. 2014), Z. japonica should be able to grow to a similar depth as Z. marina based on photosynthetic characteristics (Kaldy et al. 2015).The lower depth margin for Z. marina in Yaquina Bay is ~ -2 Mean Lower Low Water (MLLW), with most biomass centered around 0 m MLLW (Young et al. 2015).Consequently, we chose the depth interval between 0 m MLLW and the High tide mark (+ 2.3 m MLLW) to represent the potential Z. japonica colonization zone.
The relationship between Depth and its ZSI is formulated as a linear step function (Figure 2A).Breakpoints in the step-functions were determined according to the depth distributions described in Kaldy (2006).Values between the step-functions were linearly interpolated (that is, we calculated a linear equation between neighboring steps).Equations one through eight (Table 1) calculate ZSI Depth values for the model domain.
Percent Slope (ZSI % slope ) Field observations suggest that the zonation patterns of Z. japonica are related to the slope of the intertidal zone; steeply sloping intertidal zones will provide less potential habitat than broad, relatively flat areas (Shafer et al. 2014).In the ZJHSM, we represented this phenomenon as the percent slope calculated between the upper limit of the bathymetry data (+2.3 m MLLW) and the mean lower low water lines (% Slope).Table 1.Equations for used for calculating ZSI values for each of the variables in the ZJHSM.

Variable
Equation Eq# Depth (m) Depth≤ 0.6 or Depth> 2.4 0 (1) 0.6 <Depth ≤ 0.7  The relationship between % Slope and its ZSI is formulated as a linear step function (Figure 2B).Breakpoints in the step-functions were determined by expert opinion based on field surveys and observations (Shafer 2007).Values between the step-functions were linearly interpolated (that is, we calculated a linear equation between neighboring steps).Equations nine through twelve (Table 1) calculate ZSI % slope values for the model domain.

Salinity (ZSI salinity )
Experimental research on Z. japonica physiology has shown that this species is extremely euryhaline, although it does have optimal growth and photosynthesis at salinity 20 (Shafer 2007;Shafer et al. 2011;Kaldy and Shafer 2012).Although Z. japonica is capable of tolerating short-term fluctuations in salinity over the course of daily tidal cycles, sustained periods of salinities under 10 may stress Z. japonica (Shafer et al. 2011;Kaldy and Shafer 2012).We hypothesized that sustained exposure to low salinity conditions over three to four months would limit Z. japonica growth.Because salinity values change seasonally, we categorized the salinity data into two seasons: October through April (the rainy season) and May through September (dry season), and calculated an average salinity value during each season from daily modeled values.The relationship between Salinity and its ZSI is formulated as a linear step function (Figure 2C).
Values between the step-functions were linearly interpolated (that is, we calculated a linear equation between neighboring steps).Equations 13 and 14 (Table 1) calculate ZSI salinity values for the model domain.
Habitat Suitability (HSI) The suitability of a particular location for Z. japonica invasion (Habitat Suitability Index: HSI) is determined as the geometric mean of the ZSI values for the three component variables (similar to the restoration suitability index calculated by Pollack et al. (2012).If any component ZSI is 0 (unsuitable), then HSI is 0 and the habitat is considered unsuitable for invasion.HSI is calculated as:

Data Collection and Preparation
Several data sets developed by the US EPA, Western Ecology Division, Newport, OR were used in preparation and standardization of the three environmental variables (depth, percent slope and salinity) to be utilized in the creation of habitat suitability index raster files.The Yaquina Bay Topobathy Digital Elevation Model (DEM) was compiled using the Topogrid command in Arcinfo workstation 7.2.1 based on the ANUDEM 4.6.3method for hydrologically correcting digital elevation models (Young et al. 2008).Model data were interpolated from soundings and independent verification found intertidal the DEM bathymetry was accurate to ± 0.25m (P.The percent slope values were determined by calculating the degree of maximum elevation change from the topobathy DEM using the Slope tool in ArcGIS®, and the resultant values were stored as a percent for each cell within the raster file (ESRI 2011 andUSEPA 2003).The percent slope raster was output to a 2 m horizontal resolution for model data consistency.
Salinity variable values were derived by modeling the tidal range within Yaquina River study area by (1) the collection of South Beach sonde tidal data, (2) the division of the project area into five zones, and (3) the application of salinity information from three separate sonde locations throughout the Yaquina Bay study area (Critesers, Oregon oyster, and Oregon State University (OSU) Dock) (Figure 3)) to those zones.Tidal data from the South Beach sonde were originally obtained in 1hr increments from NOAA (2013) and documented the tidal state at the mouth of the Yaquina Bay estuary.The tidal data from a series of diurnal cycles were then graphed together with the 2004 temperature data from a series of temperature stations located in Yaquina Bay.Details of the sediment temperature dataset are described in Specht (2013).Continuous records of sediment temperature at 15 min intervals were made with in-situ solid-state temperature loggers (Onset Stowaway® TidbiT® Temperature Logger) placed at 5 cm depth.The temperature peaks and troughs were compared to the tidal peaks and troughs to estimate the tides occurring throughout the bay.This was done to increase the fitness of a bathtub-type tide modeling approach which uses only two variables, elevation and tidal data, to map inundation levels (NOAA 2010; Schmid et al. 2013).Next, five zones were derived for the model area by calculating river distances from the mouth of the bay to the project area extent, and dividing the area into equidistant sections (Figure 3).
Finally, continuous records of water column salinity collected every 15 min during 2007 at OSU dock, Oregon Oyster and Critesers (Figure 3) using YSI 6000 and YSI 6600 sondes were used as model inputs.Pre-and post-deployment calibrations ensured that sonde data met quality assurance requirements.Since only three sonde locations were available in the model area, the salinity data for each sonde were grouped and averaged to 1 hour temporal resolution and mapped to the appropriate zone based on following calculations.The salinity values were then divided into two separate datasets based on the temporal range for the dry season (a) Season 1 = May-Sept 2007, and the rainy season (b) Season 2 = Jan-April 2007 and Oct-Dec 2007 (Kaldy 2006).Then, the salinity data in each zone was averaged for the season based on weekly salinities.A raster dataset was then generated from the polygons for each seasonal salinity range using the same spatial extent, horizontal resolution (2 meters), and projection as the depth and percent slope variable raster files.

Habitat Suitability Index Model Variables
Raster files were created for each environmental variable (Depth, Salinity, %Slope) by applying the respective ZSI curves to each grid cell.The spatial domain was limited to study sites where Z. japonica is known to occur.Subsequently, all individual environmental ZSI raster were clipped to the model area as shown in Figure 1 and output to the geographic coordinate system; UTM, Zone 10, North NAD 1983.The HSI output values were calculated by applying equation 1 to the HSI rasters via ArcGIS® Raster Calculator.HIS results were divided into four categories: 0≤HSI≤0.25,0.25<HSI≤0.50,0.50<HSI≤0.75,and 0.75<HSI≤1.0,representing low, low-medium, medium-high, and highly suitable habitat, respectively.

Model Evaluation and Statistical Analyses
Observed Z. japonica habitat distribution was mapped by EPA for each of the following seven years : 1997: , 1998: , 1999: , 2001: , 2004: , 2007: , and 2009: (Clinton et al. 2007) and was used to determine if the ZJSHM habitat quality estimates were associated with Z. japonica presence (our hypothesis was that cells with higher HSI scores would be associate with Z. japonica presence).The spatial extent of the mapped distributions for 1997, 1998, 1999, and 2007 comprise the entire model area, while the spatial extent in 2001, 2004, and 2009 was smaller, and did not cover the central to eastern portion of the model area (Figure 1).
We evaluated our model using known presence/ absence data for Z. japonica collected in Yaquina Bay by developing several contingency tables (=confusion matrices) that explored associations between actual Z. japonica occupancy, and the HSI values our model generated.We used field data collected in 2007 for our analysis since it was the most extensive dataset.First, we developed a 2 × 2 contingency table that tested the conditional independence of habitat quality and Z. japonica occupancy.The modeled HSI values were aggregated and summarized into two categories: high quality habitat (HSI ≥ 0.5) or low quality habitat (HSI < 0.5).The presence or absence of Z. japonica was recorded for each cell and aggregated across the spatial domain.A Pearson's χ 2 goodness of fit test was used to determine if Z. japonica occupancy was conditionally independent of habitat quality.Next, using a 2 × 4 contingency table (Agresti 2012), we determined if each HSI bin was conditionally independent of seagrass presence or absence.Pearson's χ 2 goodness of fit test was used to determine conditional independence of Z. japonica and habitat quality.This method provided a quantitative method to further evaluate whether our model could determine a relationship between model-generated HSI values, and field-collected Z. japonica distributions.Methods followed Agresti (2012) and Simonoff (2003), statistical analyses used α = 0.10, and all analyses were conducted in JMP Pro 11.0.0 (SAS Institute).

Model Sensitivity Analyses
In order to determine how robust the overall HSI equation was to the variables and parameter values included in the model, three levels of sensitivity analyses were conducted.First, overall model sensitivity was determined by applying Equation (1) to two scenarios: 1) increasing all the variables by 10%, and 2) then again by decreasing all the variables by 10%.The amount of area (in ha) for each of HSI bins was recorded for the altered models, and statistical differences between each altered model and the baseline model were determined using a χ 2 goodness of fit test.
Second, we systematically removed a parameter from the model in order to determine how the predicted model index values responded to the exclusion of specific model parameters.The sensitivity analysis evaluates how the total area occupied by specific HSI bins changed value from a 3-parameter model to a 2-parameter model scenario.Specifically, we determined the total area of HSI values affected by removing one parameter (e.g., XX hectares of suitable habitat (HSI-values from 0.75-1.0)changed to unsuitable with the removal of one parameter), which reflects the relative importance of each parameter to the model structure.We increased the spatial domain of the model to include the channel and surrounding areas so that we could have a larger area over which to test the model.This method quantifies the importance of each parameter to overall model performance and allows for a more rigorous evaluation of model outputs (Pollack et al. 2012;Swannack et al. 2014).We chose a representative year (2007) for our sensitivity analysis to explore general trends in model sensitivity.
Finally, in order to determine how sensitive overall HSI values were to variable inclusion, each possible combination of 2-parameter model was run (e.g., salinity and slope, depth and salinity, etc.), and results were compared to the baseline model using a χ 2 goodness of fit test a specific variable to model performance.All statistical analyses used α = 0.10 and were conducted using JMP Pro 11.0.0 (SAS Institute).

Results
In the Yaquina Bay model area, depth ranged from -12.5 to +2.3 m MLLW (Figure 4), percent slope ranged from 0% to 64.2% (Figure 5), and the average weekly salinity conditions for the 2007 rainy season (October-April) ranged from 7.0 to 27.4 (Figure 6).Salinities were only limiting for Z. japonica from October through April (the rainy season).During the dry season (May -September), ZSI salinity values were all 1, so we limited our analyses to include only the time period from October -April.The general trend shows lower and more variable salinities in the upper part of the Yaquina Bay, as well as steeper slopes and narrower intertidal zones.Whereas, the central and lower portions of the estuary transition into wider, more shallow intertidal zones with higher salinities due to the hydrologic connection between the Yaquina River and the Pacific Ocean.
The higher HSI values generally occurred throughout the model area within the mid to upper portion of the intertidal zone with larger expanses occurring mostly in the lower estuary within Sally's Bend and Idaho Flat (Figure 7).While the upper estuary did contain some suitable habitat, most areas were not as large as those in the lower estuary due to inappropriate depth, a steeply sloping intertidal zone, and lower salinity.The lowest HIS values occurred below the lower intertidal zone, specifically within the Yaquina River channel.Although salinity conditions were suitable, it ranked low as a result of increased depth.
Z. japonica is consistently more prevalent in the lower estuary model area than in the upper estuary.The largest distribution of mapped Z. japonica was 22.25 ha during 2007 (Figure 1), while the smallest areal extent occurred in 1999 which contained 1.07 ha (Table 2).Even though some years cover a smaller extent, mapped Z. japonica in 2009, 2004, and 2001 represents the second (9.06 ha), third (7.26 ha) and fourth (3.17 ha) largest distributions respectively.
For each of the seven years, we superimposed the ZJHSM to distinguish the habitat suitability among sites (Brooks 1997).Comparison between model results and the field-collected distributions of Z. japonica indicated that the habitat suitability rankings of the ZJHSM were relatively accurate at predicting spatial distribution of Z. japonica in Yaquina Bay (Table 2).For example, in 2007 there was 22.25 ha of Z. japonica mapped, and the ZJHSM projected that 89.8% was in medium to highly-suitable habitat (i.e., HSI values from 0.5-1), and 10.17% was present in the low-medium category.Most of the mapped Z. japonica distributions (45.3 ha or 94.5%) was found within the med-high and high HIS value zones, which supports the model's ability to identify and predict areas at risk for Z. japonica invasion.Very little of the mapped Z. japonica distributions (2.6 ha or 5.5%) was found within the low-medium and low HIS value zones.For all years except 2007, less than 0.3 ha of Z. japonica was mapped in areas predicted to be low quality habitat (0<HSI<0.25)(Table 2).2).The more representative mapping years 1997, 1998, 1999 and 2007 which cover the entire model area range from a mean HSI value of 0.77 to 0.79 with the 2007 mapped year containing the largest mapped dataset (Table 2).The contingency analysis of the 2 × 2 table (performed on the 2007 data) indicated that there was an association of Z. japonica with habitat quality (χ = 10.916,p = 0.001) (Figure 8A).Likewise, a 2-tailed Fisher's Exact Test indicated that the probability of absence was greater for low-medium quality habitat than for the higher quality habitat (p = 0.0214).Model results were not perfectly predictive.Specifically, ten percent of Z. japonica occupied the low-medium habitat category (Figure 8A).Further, there were over 97 ha of unoccupied quality habitat in Yaquina Bay.When the model results were subdivided into four bins based on HSI score (low: < 0.25, low-med: 0.25-0.5,med.-high: 0.5-0.75, and high > 0.75), the test of independence was rejected (χ = 12.669, p = 0.0054), indicating that high quality habitat was more associated with Z. japonica presence than absence (Figure 8B).The collective sensitivity analyses indicated that the low-med.category was the most influenced by changes in parameter values, and that Depth had the most influence on model results, which was expected given the correlation between depth and aquatic vegetation (Table 3).When Depth was increased by 10%, the low-med.category increased by 29.9% (Figure 9A), but other categories did not change much (< 10%), when Depth was decreased by 10%.Likewise when Depth removed from the model, 1,056 ha of low quality habitat was reallocated to better habitat and overall values of HIS increased (Table 3).Specifically, there was a 3,266% increase in low-medium habitat (an increase from 19 ha to 662 ha), a 260% increase in medium-high (168 ha to 260 ha), and a 132% increase in highly suitable habitat (242 ha to 564).In contrast, HIS values were less sensitive to the removal of Salinity or %Slope (Table 3, Figures 9B and C).The low-med category increased by 14% when Salinity was decreased by 10%, and increased by 15% when Salinity was increased.When Salinity was removed, a total of 115 ha changed, and overall HIS decreased.There was a negligible increase in low quality habitat (0.52%, 6 ha).Highly suitable habitat decreased by 14% (35.8 ha), and medium-high habitat decreased by 47% (79.5 ha).These reductions of higher quality habitats were all repartitioned into the low-medium category (Table 3).When %Slope was altered by  ±10%, there was less than 6% change in overall area across all categories, with the least change in low quality habitat (Figure 9C).When %Slope was removed, 123 ha of habitat changed and overall suitability increased; eight and 115 ha of low-medium and medium-high habitat, respectively, were converted to 121 ha of highly suitable habitat, resulting in a 50% increase of highly suitable habitat.When the model was applied after each variable was increased by 10%, the Low-Med category increased by 18%, whereas there was a negligible shift (< 1%) in the other categories (Figure 9D).Decreasing the value of each variable by 10% decreased habitat in each of the four habitat categories, except Low, which increased by 1%.

Discussion
For natural resource agencies faced with combatting the spread of invasive species, understanding where and how these species will disperse is crucial.This issue is further complicated when species habitat is dynamic and influenced by large scale processes such as climate change.
These rapid and stochastic drivers can make it difficult to effectively target areas for monitoring of population expansion.The ZJHSM is a generalized model for determining potential locations suitable for colonization by Z. japonica, including those areas where the actual distribution of the species is unknown.Our approach focused on identifying the spatial continuum of habitat quality for Z. japonica¸ and we were able to elucidate an association between Z. japonica presence and higher quality habitat, doing so with a three parameter HSI model, representing critical important factors controlling its distribution of this species.

Model Variables and Sensitivity
Although many environmental variables influence seagrass growth and distribution, we were able to utilize a minimal subset and accurately predict Z. japonica distribution.However it is critically important to quantify the impact of variables by thoroughly exploring the parameter space through rigorous evaluative techniques, such as sensitivity analysis (Brooks 1997;Pollack et al. 2012).Our sensitivity analyses concluded that Depth was the most important model factor controlling the predicted distribution, and that marginal habitat is the most likely to change if parameter values change.Depth is recognized as a critical factor for the growth and distribution of submerged aquatic vegetation due to limited light availability at greater depths (Dennison 1987;Duarte 1991;Koch 2001).The ZJHSM was most sensitive to Depth being removed or altered.When Depth was removed, the overall amount of suitable habitat increased because of the increased number of cells that had suitable slopes and salinity values.Likewise, increasing Depth by as little at 10% can increase lower quality habitat by 30%.Since Z. japonica is confined to the mid-and upper intertidal zones (Kaldy 2006;Shafer et al. 2014;Young et al. 2015), and since the habitat quality projections are sensitive to changes in depth, it is essential to include water depth in order to accurately project the spatial distribution of the species.Z. japonica populations are physiologically tolerant of a wide salinity range (Shafer et al. 2011;Kaldy and Shafer 2012) and can be seasonally exposed to and tolerate low salinity from direct precipitation and from flood events (Kaldy et al. 2015b).Removing Salinity did not greatly affect overall HSI values since salinity was in the highly suitable range throughout most of the study area, except in the upper reaches of the Yaquina River.Although salinity gradients may not strongly contribute to Z. japonica suitability in Yaquina Bay, salinity is known to be an important factor affecting plant distributions in estuarine environments, so it is likely an important driver for the spatial dynamics of the species at larger scales.Experimental data indicate that seasonal exposure to pulses of low salinity or fresh water could stimulate germination of the Z. japonica seed bank (Kaldy et al. 2015b), suggesting that periodic exposure to low salinity conditions could be important in initiating new populations from seed.
Although the model accurately predicts the associations between habitat quality and Z. japonica distribution in Yaquina Bay, there are numerous caveats.Several reviews highlight the impact of omission and commission errors that can influence ecological niche modeling and geographic prediction efforts (Fielding and Bell 1997;Peterson et al. 1999).Omission of areas actually inhabited represents a failure of the modeling effort to extend to all ecological conditions which support populations.Commission errors are the erroneous inclusion of areas actually uninhabited.Commission error includes two components: real commission error, in which combinations of ecological conditions not actually within the species' niche are included, and apparent commission error, resulting from species' absences owing to interspecific interactions and historical factors, including limited colonization ability, speciation patterns, and local extinction (Peterson et al. 1999).Apparent commission error represents a real feature of species' distributional ecology: not all habitable areas are inhabited (Peterson et al. 1999).More recent work suggests that for models involving invasive species, model testing and evaluation should focus on predicting presence with a high level of confidence and place less emphasis on efficient predictions of species absence (Jimenez-Valverde et al. 2011).Consequently, our focus has been centered on developing a model that provides robust and accurate distribution predictions.We attempted to reduce error by using a multi-year dataset, and also by including the entirety of Yaquina Bay in our model domain.Our results did show that Z. japonica inhabited poorer quality habitats (HSI < 0.5), but it is not clear if those areas are actually high quality habitat misclassified as poor quality, or if Z. japonica can survive across a gradient of habitat quality.Further, ZJHSM returned high quality habitat where no Z.japonica was present -these results could be a misclassification of habitat quality, or could represent areas of high quality habitat that have not been colonized.If the latter is accurate, our model can be used as a screening tool to identify areas that are at higher risk for invasion.
Although this model was developed to try to accurately quantify Z. japonica habitat quality as parsimoniously as possible, it is flexible and other variables could easily be integrated into the model framework as data become available.For example, physiological studies have shown that temperature is an important factor affecting Z. japonica distribution (Shafer et al. 2008;Kaldy and Shafer 2012), and a previous model used temperature to predict occurrence of Z. japonica in Humboldt Bay, CA (Weltz 2012).Temperature has recently been suggested to provide a control on Z. japonica growth rates (Kaldy et al. 2015a).However, due to the limited spatial coverage of existing temperature data in Yaquina Bay (Specht 2013), temperature data were not directly used in this model, but were used to validate the close correlation between temperature and tidal elevation (e.g.depth).
Incorporating a temperature duration dependent growth factor (Kaldy et al. 2015a), or salinity dependent seed germination (Kaldy et al. 2015b) may be future refinements that increase the model predictive capacity.Shoreline orientation may also influence Z. japonica distribution through an indirect effect on temperature; southfacing shorelines would likely be warmer than shorelines oriented in other directions.Substrate type could also be considered as an additional variable in future versions of the model.Most seagrasses, including Z. japonica, inhabit soft, unconsolidated sediments.Within the Yaquina estuary, there are relatively few areas that exhibit a rocky substrate that would preclude Z. japonica colonization, therefore the inclusion of this variable would have added little discriminatory capability.However, in other areas with more heterogeneous substrates, including this parameter in the model could be a useful addition.
Future research could also explore the functional forms of the component HSI relationships as well as considering weighting different components.
At this point, the model does not include a recruitment/colonization component.Although we acknowledge the importance of propagule dispersal and recruitment in the establishment of new populations and potential range expansions, we know very little about the processes of seed dispersal and recruitment in seagrasses (Orth et al. 2006).Wind, currents and migratory waterfowl have been suggested as possible mechanisms of Z. japonica seed dispersal (Shafer et al. 2014).Ongoing research is being conducted to investigate the factors affecting seed germination in Z. japonica, and it might be possible to add a seed recruitment component in the model in the future.

Figure 1 .
Figure 1.Mapped distribution of Zostera japonica in Yaquina Bay, Oregon in 2007.

Figure 2 .
Figure 2. Relationship between Zostera japonica habitat suitability indices and (A) depth within Yaquina Bay (Depth), (B) degree of elevational change from mean high water to mean low water (% Slope), and (C) mean salinity during non-summer months.

Figure 9 .
Figure 9. Results of sensitivity analyses quantifying differences in area affected by altering combinations of parameter values.Results are reported as % change from the baseline model.(A) (B), and (C) represent results from changing one variable (Depth, Salinity, and %Slope, respectively) by +10% (gray bars), or -10% (black bars) while leaving the other two unaltered.(D) represents a paramterization where all variables were altered by +10% (gray bars) or -10% (black bars)Sensitivity analyses was performed for 2007 and applied throughout the entire modeled region.

Table 2 .
Equations (A) Comparison of seven years of the total number of hectares of HSI values divided into four suitability categories compared to total area of field-mapping efforts and (B) mean ± 1 S.D. of overall HSI values across years.

Table 3 .
Results from sensitivity analysis comparing the number of hectares occupied in each of four ISI bin of the full 3parameter model (Full Model) and respective 2-parameter models where one variable was removed (e.g., -Depth is the model with Depth removed).The numbers in parentheses are the percent change in area from the full model to the 2-parameter.