Next Article in Journal
Diversity of Olfactory Responses and Skills in Astyanax Mexicanus Cavefish Populations Inhabiting different Caves
Previous Article in Journal
Patterns of Rotifer Diversity in the Chihuahuan Desert
Previous Article in Special Issue
Biological Control of Salvinia molesta (D.S. Mitchell) Drives Aquatic Ecosystem Recovery
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Prioritizing Management of Non-Native Eurasian Watermilfoil Using Species Occurrence and Abundance Predictions

1
Bureau of Water Quality, Division of Environmental Management, Wisconsin Department of Natural Resources, 101 S Webster Street., Madison, WI 53707, USA
2
Bureau of Water Quality, Division of Environmental Management, Wisconsin Department of Natural Resources, 107 Sutliff Avenue, Rhinelander, WI 54501, USA
3
Bureau of Science Services, Wisconsin Department of Natural Resources, 2801 Progress Road, Madison, WI 53716, USA
4
Center for Limnology, University of Wisconsin-Madison, 680 N Park Street, Madison, WI 53706, USA
*
Author to whom correspondence should be addressed.
Diversity 2020, 12(10), 394; https://doi.org/10.3390/d12100394
Submission received: 11 August 2020 / Revised: 19 September 2020 / Accepted: 19 September 2020 / Published: 13 October 2020
(This article belongs to the Special Issue Ecology of Invasive Aquatic Plants)

Abstract

:
Prioritizing the prevention and control of non-native invasive species requires understanding where introductions are likely to occur and cause harm. We developed predictive models for Eurasian watermilfoil (EWM) (Myriophyllum spicatum L.) occurrence and abundance to produce a smart prioritization tool for EWM management. We used generalized linear models (GLMs) to predict species occurrence and extended beta regression models to predict abundance from data collected on 657 Wisconsin lakes. Species occurrence was positively related to the nearby density of vehicle roads, maximum air temperature, lake surface area, and maximum lake depth. Species occurrence was negatively related to near-surface lithological calcium oxide content, annual air temperature range, and average distance to all known source populations. EWM abundance was positively associated with conductivity, maximum air temperature, mean distance to source, and soil erodibility, and negatively related to % surface rock calcium oxide content and annual temperature range. We extended the models to generate occurrence and predictions for all lakes in Wisconsin greater than 1 ha (N = 9825), then prioritized prevention and management, placing highest priority on lakes likely to experience EWM introductions and support abundant populations. This modelling effort revealed that, although EWM has been present for several decades, many lakes are still vulnerable to introduction.

Graphical Abstract

1. Introduction

Non-native species are a leading driver of global environmental change. They can alter ecosystem structure and function and decrease global biodiversity [1,2,3,4]. They are economically costly and can pose hazards to human health [5,6]. Invasive species have been recorded from over half the extant phyla and divisions, and their modes of impact are as diverse as the invaders themselves [7,8,9,10]. Thus, the vulnerability of any given area to invasive species is a central concern for ecologists and natural resource managers.
Assessing lake-specific invasion vulnerability requires understanding three filters that mediate species invasions. The first filter relates to likelihood of species introduction, the second includes the probability that the introduced species establishes a self-sustaining population, and the third relates to its likely impact [11]. In other words, lakes in which a non-native species is likely to arrive, survive, and thrive are considered vulnerable. Predicting which lakes are vulnerable is helpful for management; it allows the informed and efficient application of limited resources for prevention and control, thereby minimizing adverse ecological and economic effects.
Predicting whether a lake is vulnerable first requires predicting where that species is likely to occur. Species occurrence is mediated by the first two invasion filters related to a species’ ability to arrive and survive. Important factors to consider include the potential for propagules to disperse, climate, and the availability of critical resources; these are reviewed elsewhere [12]. By mathematically relating species occurrences and spatial and environmental characteristics within a known sample of lakes, we may predict species distributions for the population [13]. Including dispersal-related factors is especially important when a species range is changing, as is the case for many non-native species, but dispersal is not routinely accounted for in species distribution models [14,15,16].
The second critical step in assessing lake-specific invasion vulnerability is to determine whether the species is likely to have adverse effects. The effects of non-native species are generally more difficult to forecast than their distributions and, as a result, they are not often predictively modelled [17,18,19]. A promising solution is to use abundance as a proxy for impact because it is more tractably modelled as a function of explanatory variables [14]. For most invasive species, abundance and impact are positively related, although the precise shape of that relationship can take different forms [19,20,21,22].
The central goal of this study was to describe lake-specific vulnerability to Eurasian watermilfoil (EWM), a non-native invasive macrophyte. We extended prior work to predict EWM occurrence by adding an empirical abundance model [23,24,25]. We predicted species occurrence and abundance as a function of water quality, land use, dispersal potential, geology, and climate variables. Predictors were selected for their link to EWM environmental suitability or aquatic invasive species (AIS) occurrence [24,26,27,28]. We then united predictions of occurrence and abundance in a prioritization framework to identify lakes at risk for developing abundant EWM populations. We separated vulnerable lakes into tiers of increasing management priority and thereby offer a simple tool to prioritize prevention efforts and management actions that reduce the impact of EWM.

2. Materials and Methods

2.1. Invasion History of Eurasian Watermilfoil

Eurasian watermilfoil is an invasive aquatic plant that can grow to high abundance in certain freshwater systems. Native to Europe, Asia, and North Africa, the precise date of introduction to the United States is unknown but probably occurred between 1880 and 1940 [29]. It has since spread throughout the continental US and Canada [30]. In Wisconsin, it is present in over 800 waterbodies and has steadily expanded its range northward following introduction in the south-central region during the 1960s [31]. When at nuisance levels, EWM forms thick mats at the surface that prevent navigation, reduce property values, and affect native species [32,33].
Eurasian watermilfoil is a costly and time-consuming management challenge [34]. From 2018 to 2020, over half of Wisconsin’s total funding for Aquatic Invasive Species control grants supported EWM control, amounting to over $500,000 USD annually. This figure excludes additional private expenditures for prevention, planning, and control. Eurasian watermilfoil management can also have significant non-target effects. Common management tools such as herbicides and mechanical harvesting have been associated with negative ecological effects on native plant communities [35,36,37]. Given the high economic and ecological cost of the species and its management, preventing the spread of EWM to uninvaded lakes is a priority. However, containing the species to the more than 1000 locations in which it is found in Wisconsin, or shielding the other 15,000 lakes from future introduction are both efforts that vastly exceed available funding. We used predictive modelling to identify vulnerable lakes where EWM is likely to be introduced and become abundant. We hope our findings will increase the efficiency of prevention and control funding, focusing it on the most vulnerable uninvaded waterbodies.

2.2. Occurrence and Abundance of Eurasian Watermilfoil

We used aquatic plant survey data obtained from the Wisconsin Department of Natural Resources (WDNR) to generate EWM occurrence and abundance models. WDNR staff and partners working under a long-term aquatic plant management research and monitoring program conducted standardized and repeatable plant surveys on a total of 657 waterbodies distributed across Wisconsin’s three lake-rich ecoregions. Lake surface area in the survey dataset ranged from 1.36 to 3958 ha, while watershed land use ranged from pristine to nearly entirely developed [38]. Surveyors performed their work during growing seasons from 1 May to 1 October 2005–2016, employing a grid-based point-intercept survey method to observe species’ presence/absence at many sites per lake. The sampling window was selected to minimize the variation in abundance over time during the growing season and thus reduce the influence of seasonality on abundance estimates (unpublished data). The number of points sampled per lake ranged from 32 to 4149, with a mean of 406 points [39]. We estimated species-specific abundance as littoral frequency of occurrence in the littoral zone using the following method. At each sampling point, observers used a double-sided bow rake attached to a 4.5 m pole to collect macrophytes from a ~0.3 m2 area. A similar rake head attached to a rope was used to collect macrophytes from sites deeper than 4.5 m [40]. All live macrophytes detached by the rake were identified to species [41,42]. We then defined littoral zones per lake based on sampled depths that were equal to or shallower than the 99th percentile of ordered depths at which macrophytes were observed. On average, 234 sample points fell within lake littoral zones. For each species, we calculated abundance as the number of occurrences divided by the total number of littoral points in each lake. To reduce the number of false absences resulting from sparse EWM populations that were below the plant survey’s limit of detection, we augmented the dataset’s presence/absence observations using a list of verified population records obtained from the WDNR. Verified records are confirmed by management staff who review observations originating from several sources, including citizen reports, routine monitoring by WDNR staff, and formal Aquatic Invasive Species detection surveys. This procedure added presence observations for 92 surveyed lakes with a 0 for EWM abundance, resulting in a total number of 385 lakes with verified populations out of a total of 657 surveyed.

2.3. Explanatory Variables

Our goal was to predict EWM occurrence and abundance on lakes statewide. We compiled information on factors thought to predict EWM occurrence and abundance if they were available for at least 7000 of the 9285 Wisconsin lakes greater than 1 ha surface area. Predictors represented a suite of factors related to water quality and lake morphometry as well as regional variables related to dispersal, land use, geology, and climate (Table S1, Supplementary Materials). We used watersheds delineated by Menuz et al. to calculate variables expressed at the watershed scale [43]. We extracted climate variables at lake centroids from WorldClim (worldclim.org). We described watershed geological characteristics using the whole-rock percentage of calcium oxide (CaO) in near-surface geology and soil erodibility (K-factor) [44,45]. High surface rock calcium oxide content relates positively to surface water conductivity, acid neutralizing capacity, and calcium content, whereas highly erodible soils are linked to increased runoff, sedimentation, and nutrient loading [44,46]. We calculated percent agriculture (crops and pasture) and percent urban land use in the watershed [47]. As a proxy for dispersal potential or propagule pressure, we computed the mean distance (km) between a lake’s centroid and all other Wisconsin lakes with EWM [48]. Factors related to vector pressure included nearby road density and lake surface area [26,49]. We calculated the density of vehicle roads (m/m2) in a 500 m buffer around each lake [50] and lake surface area extracted from the 24k hydrography dataset and maximum depth from the WDNR Register of Waterbodies [51,52]. Spatial analyses were conducted using the R packages ‘rgdal’ (v. 1.2–5), ‘raster’ (v. 2.5–8), ‘sp’ (v. 1.2–4), and ‘rgeos’ (v. 0.3–22) and ArcGIS (v. 10.2.2, Environmental Systems Research Institute, Redlands, CA, USA) [53,54,55,56]. Finally, a comprehensive database of limnological parameters provided water conductivity (μS/cm), alkalinity (mg CaCO3), pH, and satellite-estimated Secchi depth (m) for Wisconsin lakes [57].
Missing-at-random values comprised less than 3% of all observations. We imputed missing variables using predictive mean matching and the package ‘mice’ (v. 2.30) over 50 iterations [58]. We log-transformed highly skewed numeric variables and square-root-transformed skewed percentages (see Table 1). We computed variance inflation factors (VIFs) for each variable in the dataset using the function ‘vif’ in the package ‘car’ (v. 2.1–4) [59]. We sequentially excluded variables with the largest inflation factor until no inflation factor exceeded 10.

2.4. Predicting Eurasian Watermilfoil Occurrence

We built species distribution models using logistic regression in a generalized linear modelling framework applied to the EWM occurrence dataset. The procedure employs a maximum likelihood optimization algorithm to estimate intercept and slope parameters ( β 0   ,   β j ) for a set of j predictors ( X ) to determine the probability ( p ) that a given lake ( i ) has been invaded. The equation linearizes the response variable via a logit transformation.
y i =   logit   ( p i ) =   ln   p i 1 p i = β 0 + j = 1 n β i j X i j
The probability is subsequently calculated as follows:
P ( E W M   presence ) = e y / ( 1 + e y )
Model fitting was performed using Firth’s method of bias reduction in R using the function and package ‘brglm’ (v. 0.5–9) [60].
We used a 5-fold cross-validation procedure repeated 10 times to evaluate model performance. We randomly split the data into 5 approximately equal parts (termed a fold) and developed a model five times, once per each unique combination of N = 4 folds. For each run, we generated predicted values by applying the resulting model on the remaining fifth fold [61]. After each cross-validation, we evaluated model performance. First, we conducted a receiver operating characteristic analysis as follows [62]. The logistic regression equation produces an estimate of the lake-specific probability of EWM presence. One can use the probability of presence to determine predicted presence by selecting a probability threshold above which one would assume the species is present. The threshold may be set at any point along the range in probability from 0 to 1. Starting at 0, we incrementally increased the threshold value and compared the resulting list of predicted presences to the observed presences, plotting the percentage of true presences against true absences. A 1:1 relationship between true presences and true absences would describe a model that is no better than random chance; the area under the 1:1 line is 0.50. One may plot a curve relating true absences and presences (AUC) as the threshold increases from 0 to 1. Values for the area under the curve (AUC) above 0.5 reflect predictive power that is better than chance, and a value of 1 describes a model with perfect discrimination. We generated these receiver operating characteristic curves using the function ‘roc’ in package ‘pROC’ (v. 1.9.1) [63]. We also calculated overall model deviance, the percentage of deviance explained by the model (D2), and Tjur’s coefficient of discrimination, which can be interpreted much in the same way as an R2 value for linear regressions [64]. We averaged all model performance statistics across the 10 repeated cross-validations.
Next, we extended the model to predict probability of presence on all lakes larger than 1 ha in surface area (N = 9285). Using the threshold values generated to evaluate the model, we selected the value 0.246 to distinguish likely presences. We used this threshold to determine which lakes were most vulnerable to EWM invasion and mapped occurrence probability using ESRI software (v. 10.2.2, Environmental Systems Research Institute, Redlands, CA, USA) [55]. Our choice of threshold allowed no more than 10% of the predicted absences to be false. Emphasizing model specificity over sensitivity minimized the chance of classifying a waterbody as “safe” when it was, in fact, vulnerable; however, see [24]. In addition, erring on the side of overpredicting may reduce bias related to survey detection failures in the occurrence data. Our approach will result in overprediction of EWM occurrence, but a less cautious prediction may be made by selecting a different threshold for the occurrence data provided in Table S2, Supplementary Materials.

2.5. Predicting Eurasian Watermilfoil Abundance

Next, we developed a statistical model predicting EWM abundance using the same environmental predictors on the 657 surveyed lakes. Exploratory univariate plots revealed curvilinear and unimodal distributions, so we included quadratic transformations for all predictors. Littoral EWM abundance was highly heteroskedastic, overdispersed, and right-skewed; accordingly, we selected the beta distribution for the model, which has a flexible shape controlled by mean (μ) and precision (φ) parameters [65,66]. The model contains two submodels, i.e., one for the mean response and one for dispersion. We used extended beta regression models with bias correction to estimate mean and precision parameters as a function of predictors, thus modelling mean response, variable dispersion, and skewness [67]. Here, for any fixed μ, greater φ relates to decreased variability in the response variable [68]. The expected value and variance of the response variable y is determined by the following:
E ( y ) = μ
v a r ( y ) =   μ ( 1 μ ) 1 + φ
Using a logit transformation for the mean submodel and a log link for dispersion, the models are specified as follows:
logit ( μ i ) =   ln   μ i 1 μ i = β 0 + j = 1 n β i j X i j
ln φ i =   θ 0 + j = 1 n θ i j Z i j
The set of predictors ( X i , Z i ) may vary by submodel but need not be mutually exclusive. We allowed all predictors to contribute to both submodels. The range of the response variable included 0, so we transformed it using the Smithson and Verkuilen method cited in Cribari-Neto and Zeileis [69,70]. Bias-corrected model fitting was performed using the function and package ‘betareg’ [69,71].
We then applied a 5-fold cross-validation procedure as described in Section 2.4. After each cross-validation, we evaluated model performance using several metrics as follows: Pearson’s correlation coefficient (r) to reflect concordance between observed and predicted values, Spearman’s rank correlation (ρ) to test concordance among value ranks, and parameters m and b from a simple linear regression to further describe the relationship between observed and predicted values. Given perfect concordance among observed and predicted values, the intercept (b) and slope parameter (m) would be 0 and 1, indicating no bias and a comparable range of observed values at all points along the range of predicted values. Lower or higher values for b indicate model under- or overprediction, while a different m reflects a bias that may differ in magnitude along the range of predictions [72]. We calculated the root mean square error among observed and predicted values. In all cases, model performance measures were averaged across 10 repeated cross-validations. We used ESRI software to generate maps of observed and predicted abundance (v. 10.2.2, Environmental Systems Research Institute, Redlands, CA, USA) [55].

2.6. Defining and Prioritizing Management Targets

We extended the occurrence and abundance models developed on the 657 surveyed lakes to produce occurrence and abundance predictions for 9825 lakes over 1 ha in size. We then split lakes into categories of increasing invasion vulnerability by trisecting the ranges for predicted probability of occurrence and predicted abundance. This resulted in three groups of lakes with low, medium, and high vulnerability for EWM occurrence defined by occurrence probability thresholds at 0.496 and 0.746. For abundance, low, medium, and high vulnerability thresholds fell at 0.18 and 0.36. By cross-tabulating the priority categories for both presence and abundance, we constructed a 3-tier management priority matrix that may help direct work toward the highest-priority prevention and management targets.

3. Results

3.1. Occurrence Models

Logistic regression models predicting EWM occurrence in surveyed lakes performed well, with mean cross-validated AUC = 0.82. Variables positively related to occurrence probability included road density, surface area, maximum air temperature of the warmest month, and lake maximum depth, while factors that were negatively related included watershed % surface rock calcium oxide content, annual temperature range, and mean distance from all source populations (Table 1). Mean cross-validated deviance was 510, Tjur’s coefficient of determination was 0.33, and the amount of deviance accounted for by the model was 23%.
Using a threshold probability of 0.246 to distinguish likely presences from absences resulted in a model that was 72.6% accurate in its classifications (Figure 1, Table 2). We selected the threshold probability to minimize false negatives, resulting in high model specificity. Overall, the true positive rate (sensitivity) was 0.372, while the true negative rate was 0.97 (specificity). Modeled predictions per lake are provided in Table S2 to allow alternate classifications.
When the model was extended to all lakes over 1 ha in surface area, it revealed the potential for EWM to continue to expand in Wisconsin. Eurasian watermilfoil is predicted to occur in 2357 lakes where it is not currently known (Figure 1d). Vulnerable waterbodies were found statewide, but lakes with the greatest vulnerability for occurrence tend to occur in the south-central and southeast regions of the state.

3.2. Abundance Models

Eurasian watermilfoil abundance models explained a small but statistically significant portion of the observed variation (log-likelihood = 1856, pseudo R2 = 0.24; Table 3). Mean cross-validated deviance was −2894, root mean square error was 15.7%, while cross-validated correlation coefficients among observed and predicted abundance showed a moderate degree of concordance, with r = 0.43 and ρ = 0.41. Mean EWM abundance generally increased with conductivity, road density, and maximum air temperature, while abundance decreased with % lithological CaO and annual temperature range (Table 3). Quadratic terms indicating curvilinear relationships were important for some predictors. When all other variables were held at their mean values, abundance generally increased with increasing conductivity, then decreased at the very high end of the range. For mean distance from source, abundance first decreased, then increased, and for soil erodibility, abundance dropped at the low and high ends of the observed range, displaying a peak at intermediate values. All other predictor variables used in the model were not statistically significant. Significant precision terms were mostly negative, indicating variables associated with decreased variability in the mean response. As conductivity, alkalinity, and % watershed agriculture increased, so did the precision in the resulting estimates. Soil erodibility had a positive linear precision coefficient, while that for mean distance was negative; however, the degree to which precision changed varied across the measured range (i.e., the relationship was curvilinear).
Observed and predicted EWM abundances were highly correlated (r = 0.51, ρ = 0.52), although cross-validated performance indicated uncertainty in modelled predictions. With respect to observed abundance, predictions were relatively unbiased (b = −0.004) and reasonably consistent (m = 0.99). In the surveyed lakes, EWM abundance ranged from almost 0 to 1, with low abundances consistently observed in northern Wisconsin (Figure 2a). Predicted abundance mirrors that of observed abundance but tends to underpredict when EWM abundance is high (Figure 2b). Modelled predictions never exceeded an abundance of 0.46, but 33 out of the 657 lakes with EWM populations had observed abundance values ranging from 0.46 to 0.96.

3.3. Statewide Predictions

After developing the occurrence and abundance models on the N = 657 surveyed lakes, we extended model prediction to 9825 Wisconsin lakes over 1 ha in surface area (Figure 3). There was again a strong clustering of lakes at risk for abundant EWM populations in the southeast corner of the state, while northern lakes were generally less vulnerable to invasion, establishment, and abundant populations. Because occurrence is generally easier to predict on a large scale than abundance, we examined the relationship between predicted probability of occurrence and abundance, finding that predicted occurrence probability was positively related to predicted abundance (F = 0.00015, p < 0.001, Radj = 0.62), but that the relationship weakened when predicted occurrence probability was used to predict observed abundance (F = 99.07, p < 0.001, Radj = 0.13).

3.4. Prioritizing Management

We combined statewide predicted occurrence and abundance estimates to describe overall lake-specific vulnerability related to EWM’s ability to arrive, survive, and thrive. We examined occurrence and abundance predictions for the 2357 uninvaded lakes where the occurrence model predicted EWM was likely to occur and assigned lakes to one of three tiers of increasing management priority (Table 4, Figure 4). Tier 3 comprises the lowest-priority lakes, including waterbodies with a relatively low occurrence risk and which are also unlikely to support abundant populations. Over 78% of lakes likely to be invaded fall into Tier 3. Of slightly higher prevention and management priority are Tier 2 lakes with a moderate risk of occurrence and abundance, a high risk of abundance balanced by a low risk of occurrence, or a high risk of occurrence balanced by a low risk of abundance. Tier 2 comprises 13% of lakes likely to be invaded. The most vulnerable lakes make up Tier 1; only 9% of lakes likely to be invaded fall into this category. In these lakes, EWM is both likely to be introduced and grow to high abundance. Tier 1 lakes occur mostly in the southeast region of the state, with vulnerability decreasing as one moves in a northwesterly direction. It is our hope that this framework along with the waterbody list presented in Table S2 can be used to inform prevention and management actions to support the efficient use of limited prevention and control resources.

4. Discussion

Eurasian watermilfoil has been present in Wisconsin for at least 60 years and there are still many lakes that are vulnerable to invasion [73]. However, very few vulnerable lakes are likely to support abundant populations. While some predictors of occurrence and abundance were the same, factors related to environmental suitability, such as conductivity and soil erodibility, were uniquely predictive of abundance, whereas factors related to visitation and recreational value, such as surface area and maximum depth, emerged as unique predictors of occurrence [24,49,74]. Overall, the most vulnerable populations were found in the south and southeast regions of the state, with vulnerability decreasing northward. Lakes with a high likelihood of occurrence were not necessarily predicted to occur at high abundance, and some lakes predicted to attain high abundance did not have a high likelihood of occurrence. Our study highlights the importance of not conflating the likelihood of occurrence and abundance, which are different aspects of invasive species’ success and influenced by different factors.

4.1. Occurrence Models

Models were developed to maximize predictive power rather than test ecological relationships, but several interesting relationships indicated by model coefficients warrant discussion. Factors associated with species transport and arrival were important predictors of EWM occurrence. Higher occurrence probability was predicted for larger and deeper lakes with more nearby roads. Larger lakes generally attract more traffic, resulting in higher propagule pressure [49]. In addition, higher spatial heterogeneity in large lakes may reduce the number of stochastic extinction events and enhance population persistence [75]. In addition, roads make lakes more accessible to humans, who are implicated as an important vector of invasive species transport [76]. Finally, invaded lakes had a smaller average distance to other established EWM populations. Dispersal probability declines with distance, such that distance from established populations is often helpful in predicting invasive species occurrence [16,77]. Distance from established populations in the occurrence model may capture constraints related to dispersal, but it may also reflect spatially auto-correlated environmental conditions that are difficult to disentangle [48].
After an invasive species arrives at a location, its survival is in part determined by local environmental conditions. Early work to model the distribution for EWM found environmental factors to be more important than those related to human activity, while a later study identified the additional importance of human activity and dispersal [23,24]. New drivers revealed by this study further add to our understanding of what controls EWM occurrence and abundance; for example, calcium oxide surface rock content predicted low occurrence probability. This variable exhibited a strong spatial pattern associated with the presence of a large dolomite deposit in eastern Wisconsin. Marl lakes occur in abundance in this region, and they have unique biogeochemical qualities. Calcium carbonate in these high-alkalinity, high-pH lakes is plentiful, but rapidly co-precipitates with phosphorus and dissolved organic material. The low concentrations of free CO2, phosphorus, iron, and manganese in the water of marl lakes can limit macrophyte growth [78]. Maximum air temperature was positively associated with EWM occurrence, whereas annual temperature range was negatively related. Eurasian watermilfoil can grow in a wide range of temperatures with populations extending up to 68.8°N latitude [79]. Temperature is unlikely to be a limiting factor for this species in the spatial extent considered by this study. In addition, these climate factors may be co-linear with a set of uncaptured variables related to environmental conditions or invasion history.
Other environmental factors such as conductivity, alkalinity, and water clarity were not significant predictors of EWM occurrence. This may be due to the species’ broad environmental tolerances and the fact that environmental variables were mostly within published tolerance ranges [28]. That said, dissolved inorganic carbon is an important nutrient source for the species; alkaline lakes are generally considered more suitable. In our dataset, EWM never occurred when conductivity was below 16 μS/cm or when alkalinity was below 4 mg CaCO3/L. Around 6% of lakes statewide had conductivity below 16μS/cm and 7.5% had alkalinity below 4 mg CaCO3/L. Additional work is necessary to determine whether this threshold is biologically relevant for EWM or an artifact of the distribution of values within our sample.
Species occurrence models typically assume a population is at equilibrium, but this assumption is violated in the case of a range-expanding species like EWM. The inclusion of true absence data along with spatial variables related to dispersal allows us to reduce the bias that would otherwise be present [16]. Still, it is likely that our training set includes occurrences that have not been detected, leading to some degree of error in our predictions. While the extrapolated statewide predictions we present can be useful in planning prevention and management activities, it is important to treat them as indicators of a likely but uncertain future state.

4.2. Abundance Models

Once probability of introduction and survival is known, the remaining question is one of impact: if the species is introduced and survives, is it also likely to cause problems? While impact is arguably the most important filter to consider, it is often the most difficult to predict [20]. Eurasian watermilfoil has been associated with several effects on native flora, macroinvertebrates, habitat, and water quality [35,80,81,82,83]. Eurasian watermilfoil has also been associated with recreational impairment and decreased lake property value [37,84,85]. That said, the magnitude of its socioeconomic and ecological effects are most likely directly related to abundance [19].
In contrast to the EWM occurrence model, several local environmental variables were uniquely predictive of abundance (Table 3). Conductivity was positively associated with M. spicatum abundance. Centered and standardized predictors allow model coefficients to be interpreted as effect sizes, and conductivity had one of the largest effects observed. The influence of conductivity is strong; it is often considered a “master factor” driving landscape-scale species distributions [74]. In particular, the strong effect makes sense for EWM, which can use bicarbonate as a source of carbon dioxide, thus attaining a competitive advantage in high-conductivity and high-alkalinity lakes [86].
Soil erodibility displayed a concave-down curvilinear relationship such that lakes in watersheds with moderate soil erodibility supported more abundant populations. EWM is tolerant of nutrient enrichment and prefers moderately enriched lakes [28]. Phosphorus in surface water is often derived from rock and soils and is generally higher when watersheds are comprised of highly erodible soil [87]. Erodible soil is also often favored for agriculture, where exogenous additions of fertilizer further enrich surface waters [88]. As nutrient content in water increases, so does primary productivity; abundant plant populations are often found in enriched waters. However, beyond a certain level of enrichment, filamentous and planktonic algae begin to dominate, resulting in water with low light transmission. Macrophytes then decrease in abundance, unable to survive in these turbid, highly enriched waters [89,90].
Temperature, calcium oxide surface rock content, and distance from established populations were significant predictors of both occurrence and abundance. Higher air temperature and lower annual temperature range were associated with higher abundance. Eurasian watermilfoil has a relatively high optimum temperature for photosynthesis, and warmer temperatures likely enhance productivity and expansion [91]. Surprisingly, as distance from invaded systems increased, so did abundance. This may be explained by a relationship between distance from established populations and time since introduction—reflecting the “boom” typically exhibited by EWM populations preceding a decline that often occurs around 10–15 years following invasion [31,92,93].

4.3. Management Prioritization: Uniting Occurrence and Abundance

Many lakes are vulnerable to EWM introduction, but of the 2357 lakes in which EWM is likely to occur, only 245 have the highest priority for prevention and management given their combined risk of occurrence and abundance. Moreover, of the 245 Tier 1 lakes singled out for action, only 29 (12%) are in the highest abundance risk category. This makes sense, given that we know species abundance distributions for non-native and native species are similar—they are “commonly rare and rarely common” [94]. This also aligns with our understanding of EWM abundance distributions—of 388 surveyed lakes with EWM populations, only 50 (13%) are in the highest abundance category (>0.36). In light of our predictions, we recommend enhancing monitoring and prevention programs on the 29 highest-priority lakes to help offset large future costs and improve management efficiency [34]. The remaining Tier 1 lakes that are vulnerable to EWM introduction but which have a lower risk for abundant populations may present opportunities for public engagement through community-based science and prevention.
A key implication of our study is that even though EWM invaded Wisconsin decades ago, there remain many uninvaded lakes that are predicted to be suitable for establishment, but significantly fewer lakes where the species is expected to reach high abundance. Given this finding, these models may be an important tool for prioritizing resources. Considering that the species has been present and spreading for decades, one might expect EWM to be largely saturated on the landscape, having already established itself in most suitable lakes. The lesson here applies to other invasive species—just because an invader has been long present in a region does not mean that it has established itself in all suitable habitat. There is value in efforts to prevent or slow further spread, even for long-established non-native species.
Management of EWM can be costly. Placing equal priority on all lakes fails to take advantage of the fact that, for most lakes, EWM abundance is low. Yet managers are justifiably risk-averse, and proactive management to keep populations small in case the population can attain high abundance is understandable. The approach we present here, to combine occurrence and abundance predictions, allows researchers and lake managers to empirically determine a management strategy appropriate for a recent introduction based on evidence; a wait-and-see approach may be appropriate when the likelihood of achieving abundant populations is low, while a more aggressive or proactive approach would make sense when the predicted abundance is high. Better justification for planning a proactive versus reactive strategy will hopefully result in more efficient allocation of limited management funding.
Understanding multiple levels of the invasion process helps generate realistic predictions of system-specific risk [11,95]. Few efforts to date have produced risk assessments that integrate occurrence and abundance, but see [19]. We could improve on this assessment of vulnerability with a more nuanced understanding of species’ impacts across systems. We selected risk thresholds to determine management priority by simply trisecting the observed range in predicted occurrence and abundance, assuming adverse effects have a simple positive linear relationship to abundance. An improvement to this work would be to select prioritization thresholds in consideration of the actual shape of the impact–abundance relationship(s) [19,21,22].
A final caveat in using these predictions to direct management effort is the following: the EWM occurrence model explained 32% of observed variation and correctly classified 75.6% of surveyed lakes. The abundance model explained 25% of the observed variation. At a statewide scale, occurrence and abundance predictions allow strategic use of limited resources despite the uncertainty in the model. However, when evaluated at the individual lake scale, error in the model should be recognized to appropriately characterize risk. The most conservative option for a less vulnerable Tier 3 lake with energetic and willing volunteers may still be to support active prevention.
In conclusion, understanding lake-specific vulnerability in light of predicted occurrence and abundance is a promising approach to effective management of AIS across a landscape. Prevention and control funding and volunteer effort are precious resources that should be directed to the most vulnerable waterbodies. Smart prevention and prudent control guided by a clear prioritization framework will save money, minimize non-target effects, and hopefully result in better management outcomes.

Supplementary Materials

The following are available online at https://www.mdpi.com/1424-2818/12/10/394/s1, Table S1 summarizes the variables used to predict EWM occurrence and abundance. Table S2 displays observations and modelled predictions for EWM occurrence and abundance for all 9287 study lakes. Also presented are the final management priority tier assigned by the prioritization framework and the relative lake-specific risk for EWM occurrence and high abundance. Lakes are identified by the WDNR waterbody identification code, lake name, and county.

Author Contributions

A.M. and E.R.K. conceptualized the work, A.M. and M.J.V.Z. designed the study, A.M. assembled and analyzed the data and designed visualizations, A.M. wrote the original draft. A.M., C.L.H., S.V.E., E.R.K. and M.J.V.Z. reviewed and edited subsequent drafts. A.M. and M.J.V.Z. secured funding for the work. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Science Foundation through the Graduate Research Fellowship Program under Grant #DGE-1256259 and the North Temperate Lakes Long-Term Ecological Research Program under cooperative agreement #DEB-1440297, and the Wisconsin Department of Natural Resources.

Acknowledgments

We would like to thank Paul Frater for his review of our mathematical notation and the staff at the Wisconsin DNR for their help collecting data, including Therese Ashkenase, Martha Barton, Dana Bigham, Shaunna Chase, Michael Fell, Paul Frater, Elizabeth Haber, Raffica La Rosa, Michelle Nault, Meghan Porzky, Erin Ridley, Chris Repking, Katie Roth, Jesse Schwingle, Nicholas Shefte, Kari Soltau, and Kelly Wagner. We thank DNR central and regional staff for supporting and participating in this work, including Tim Asplund, Heidi Bunk, Mary Gansberg, Kevin Gauthier, Susan Graham, Ted Johnson, Jim Kreitlow, Brenda Nordin, Scott Provost, Carroll Schaal, Alex Smith, and Pamela Toshner.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Mack, R.N.; Simberloff, D.; Lonsdale, W.M.; Evans, H.; Clout, M.; Bazzaz, F.A. Biotic invasions: Causes, epidemiology, global consequences, and control. Ecol. Appl. 2000, 10, 689–710. [Google Scholar] [CrossRef]
  2. Simberloff, D. How common are invasion-induced ecosystem impacts? Biol. Invasions 2011, 13, 1255–1268. [Google Scholar] [CrossRef]
  3. Vitousek, P.; DAntonio, C.; Loope, L.; Rejmanek, M.; Westbrooks, R. Introduced species: A significant component of human-caused global change. N. Zea. J. Ecol. 1997, 21, 1–16. [Google Scholar]
  4. Wilcove, D.S.; Rothstein, D.; Dubow, J.; Phillips, A.; Losos, E. Quantifying threats to imperiled species in the United States. Bioscience 1998, 48, 607–615. [Google Scholar] [CrossRef] [Green Version]
  5. Neill, P.; Arim, M. Human Health Link to Invasive Species. In Encyclopedia of Environmental Health; Elsevier BV: Amsterdam, The Netherlands, 2011; pp. 116–123. [Google Scholar]
  6. Pimentel, D.; Zuniga, R.; Morrison, D. Update on the environmental and economic costs associated with alien-invasive species in the United States. Ecol. Econ. 2005, 52, 273–288. [Google Scholar] [CrossRef]
  7. AquaNIS Editorial Board. Information System on Aquatic Non-Indigenous and Cryptogenic Species. World Wide Web Electronic Publication. 2015. Available online: http://www.corpi.ku.lt/databases/aquanis (accessed on 23 April 2020).
  8. Integrated Taxonomic Information System (ITIS). 2020. Available online: http://www.itis.gov (accessed on 23 April 2020).
  9. Invasive Species Specialist Group (ISSG). The Global Invasive Species Database. 2015. Available online: http://www.iucngisd.org/gisd/ (accessed on 23 April 2020).
  10. Simpson, A.; Eyler, M.C.; Sikes, D.; Bowser, M.; Sellers, E.; Guala, G.F.; Cannister, M.; Libby, R.; Kozlowski, N. A Comprehensive List of Non-Native Species Established in Three Major Regions of the United States: Version 2.0; U.S. Geological Survey: Reston, VA, USA, 2019. [CrossRef]
  11. Vander Zanden, M.J.; Olden, J.D.; Thorne, J.H.; Mandrak, N.E. Predicting occurrences and impacts of smallmouth bass introductions in north temperate lakes. Ecol. Appl. 2004, 14, 132–148. [Google Scholar] [CrossRef] [Green Version]
  12. Guisan, A.; Thuiller, W. Predicting species distribution: Offering more than simple habitat models. Ecol. Lett. 2005, 8, 993–1009. [Google Scholar] [CrossRef]
  13. Rushton, S.; Ormerod, S.J.; Kerby, G. New paradigms for modelling species distributions? J. Appl. Ecol. 2004, 41, 193–200. [Google Scholar] [CrossRef]
  14. Elith, J.; Graham, C.; Anderson, R.P.; Dudík, M.; Ferrier, S.; Guisan, A.; Hijmans, R.; Huettmann, F.; Leathwick, J.R.; Lehmann, A.; et al. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 2006, 29, 129–151. [Google Scholar] [CrossRef] [Green Version]
  15. Guisan, A.; Lehmann, A.; Ferrier, S.; Austin, M.; Overton, J.M.C.; Aspinall, R.; Hastie, T. Making better biogeographical predictions of species’ distributions. J. Appl. Ecol. 2006, 43, 386–392. [Google Scholar] [CrossRef]
  16. Václavík, T.; Meentemeyer, R.K. Invasive species distribution modeling (iSDM): Are absence data and dispersal constraints needed to predict actual distributions? Ecol. Model. 2009, 220, 3248–3258. [Google Scholar] [CrossRef]
  17. Dick, J.T.A.; Alexander, M.E.; Jeschke, J.M.; Ricciardi, A.; MacIsaac, H.J.; Robinson, T.B.; Kumschick, S.; Weyl, O.; Dunn, A.M.; Hatcher, M.J.; et al. Advancing impact prediction and hypothesis testing in invasion ecology using a comparative functional response approach. Biol. Invasions 2013, 16, 735–753. [Google Scholar] [CrossRef] [Green Version]
  18. Hall, J.L.W.; Anderson, R.D.; Killen, W.D.; Hosmer, A.J.; A Brain, R. Assessment of periphyton, aquatic macrophytes, benthic communities, and physical habitat in midwestern United States streams coinciding with varying historical concentrations of atrazine. J. Environ. Sci. Health Part A 2014, 49, 1091–1099. [Google Scholar] [CrossRef] [PubMed]
  19. Latzka, A.W.; Hansen, G.J.A.; Kornis, M.; Zanden, M.J.V. Spatial heterogeneity in invasive species impacts at the landscape scale. Ecosphere 2016, 7, e01311. [Google Scholar] [CrossRef] [Green Version]
  20. Parker, I.; Simberloff, D.; Lonsdale, W.; Goodell, K.; Wonham, M.; Kareiva, P.; Williamson, M.; Von Holle, B.; Moyle, P.; Byers, J.; et al. Impact: Toward a framework for understanding the ecological effects of invaders. Biol. Invasions 1999, 1, 3–19. [Google Scholar] [CrossRef]
  21. Strayer, D.L. Non-native species have multiple abundance-impact curves. Ecol. Evol. 2020, 10, 6833–6843. [Google Scholar] [CrossRef] [PubMed]
  22. Yokomizo, H.; Possingham, H.P.; Thomas, M.B.; Buckley, Y.M. Managing the impact of invasive species: The value of knowing the density-impact curve. Ecol. Appl. 2009, 19, 376–386. [Google Scholar] [CrossRef] [Green Version]
  23. Buchan, L.A.J.; Padilla, D.K. Predicting the likelihood of Eurasian watermilfoil presence in lakes, a macrophyte monitoring tool. Ecol. Appl. 2000, 10, 1442–1455. [Google Scholar] [CrossRef]
  24. Roley, S.S.; Newman, R. Predicting Eurasian watermilfoil invasions in Minnesota. Lake Reserv. Manag. 2008, 24, 361–369. [Google Scholar] [CrossRef] [Green Version]
  25. Tamayo, M.; Olden, J.D. Forecasting the vulnerability of lakes to aquatic plant invasions. Invasive Plant Sci. Manag. 2014, 7, 32–45. [Google Scholar] [CrossRef]
  26. Decker, K.L.; Allen, C.R.; Acosta, L.; Hellman, M.L.; Jorgensen, C.F.; Stutzman, R.J.; Unstad, K.M.; Williams, A.; Yans, M. Land Use, Landscapes, and Biological Invasions. Invasive Plant Sci. Manag. 2017, 5, 108–116. [Google Scholar] [CrossRef] [Green Version]
  27. Mikulyuk, A.; Sharma, S.; Van Egeren, S.; Erdmann, E.; Nault, M.E.; Hauxwell, J. The relative role of environmental, spatial, and land-use patterns in explaining aquatic macrophyte community composition. Can. J. Fish. Aquat. Sci. 2011, 68, 1778–1789. [Google Scholar] [CrossRef]
  28. Smith, C.S.; Barko, J.W. Ecology of Eurasian watermilfoil. J. Aquat. Plant Manag. 1990, 28, 55–64. [Google Scholar]
  29. Eiswerth, M.E.; Donaldson, S.G.; Johnson, W.S. Potential environmental impacts and economic damages of Eurasian watermilfoil (Myriophyllum spicatum) in Western Nevada and Northeastern California. Weed Technol. 2000, 14, 511–518. [Google Scholar] [CrossRef]
  30. EDDMapS. Early Detection & Distribution Mapping System; The University of Georgia-Center for Invasive Species and Ecosystem Health: Athens, GA, USA, 2017. [Google Scholar]
  31. Trebitz, A.S.; Nichols, S.A.; Carpenter, S.R.; Lathrop, R.C. Patterns of vegetation change in Lake Wingra following a Myriophyllum spicatum decline. Aquat. Bot. 1993, 46, 325–340. [Google Scholar] [CrossRef]
  32. Boylen, C.W.; Eichler, L.W.; Madsen, J.D. Loss of native aquatic plant species in a community dominated by Eurasian watermilfoil. Hydrobiologia 1999, 415, 207–211. [Google Scholar] [CrossRef]
  33. Provencher, B.; Lewis, D.J.; Anderson, K. Disentangling preferences and expectations in stated preference analysis with respondent uncertainty: The case of invasive species prevention. J. Environ. Econ. Manag. 2012, 64, 169–182. [Google Scholar] [CrossRef]
  34. Zipp, K.Y.; Lewis, D.J.; Provencher, B.; Zanden, M.J.V. The spatial dynamics of the economic impacts of an aquatic invasive species: An empirical analysis. Land. Econ. 2019, 95, 1–18. [Google Scholar] [CrossRef]
  35. Mikulyuk, A.; Kujawa, E.; Nault, M.E.; Van Egeren, S.; Wagner, K.I.; Barton, M.; Hauxwell, J.; Zanden, M.J.V. Is the cure worse than the disease? Comparing the ecological effects of an invasive aquatic plant and the herbicide treatments used to control it. FACET 2020, 5, 353–366. [Google Scholar] [CrossRef]
  36. Di Nino, F.; Thiébaut, G.; Muller, S. Response of Elodea nuttallii (Planch.) H. St. John to Manual Harvesting in the North-East of France. Hydrobiologia 2005, 551, 147–157. [Google Scholar] [CrossRef]
  37. Torn, K.; Martin, G.; Kotta, J.; Kupp, M. Effects of different types of mechanical disturbances on a charophyte dominated macrophyte community. Estuar. Coast. Shelf Sci. 2010, 87, 27–32. [Google Scholar] [CrossRef]
  38. Omernik, J.M.; Chapman, S.S.; Lillie, R.A.; Dumke, R.T. Transactions of the Wisconsin Academy of Sciences, Arts and Letters. Ecoregions Wis. 2000, 88, 77–103. [Google Scholar]
  39. Mikulyuk, A.; Hauxwell, J.; Rasmussen, P.; Knight, S.; Wagner, K.I.; Nault, M.E.; Ridgely, D. Testing a methodology for assessing plant communities in temperate inland lakes. Lake Reserv. Manag. 2010, 26, 54–62. [Google Scholar] [CrossRef] [Green Version]
  40. Hauxwell, J.; Knight, S.; Mikulyuk, A.; Nault, M.E.; Porzky, M.; Chase, S. Recommended Baseline Monitoring of Aquatic Plant in Wisconsin: Sampling Design, Field and Laboratory Procedures, Data Entry and Analysis, and Applications; Wisconsin Department of Natural Resources: Madison, WI, USA, 2010.
  41. Crow, G.E.; Hellquist, C.B. Aquatic and Wetland Plants of Northeastern North America. Vol 1. Pteridophytes, Gymnosperms and Angiosperms: Dicotyledons; University of Wisconsin Press: Madison, WI, USA, 2000. [Google Scholar]
  42. Crow, G.E.; Hellquist, C.B. Aquatic and Wetland Plants of Northeastern North America. Vol. 2. Angiosperms: Monocotyledons; University of Wisconsin Press: Madison, WI, USA, 2000. [Google Scholar]
  43. Menuz, D.R.; Ruesch, A.S.; Diebel, M.W. 1:24K Hydrography Attribution Metadata; Wisconsin Department of Natural Resources: Madison, WI, USA, 2013.
  44. Olson, J.R.; Hawkins, C.P. Predicting natural base-flow stream water chemistry in the western United States. Water Resour. Res. 2012, 48, 1–19. [Google Scholar] [CrossRef] [Green Version]
  45. Moebius-Clune, B. Soil Health Initiatives of the USDA Natural Resources Conservation Service (NRCS); RePEc; USDA: Washington, WA, USA, 2017.
  46. Wang, G.; Wu, B.; Zhang, L.; Jiang, H.; Xu, Z. Role of soil erodibility in affecting available nitrogen and phosphorus losses under simulated rainfall. J. Hydrol. 2014, 514, 180–191. [Google Scholar] [CrossRef]
  47. Jin, S.; Yang, L.; Danielson, P.; Homer, C.; Fry, J.; Xian, G. A comprehensive change detection method for updating the National Land Cover Database to circa 2011. Remote Sens. Environ. 2013, 132, 159–175. [Google Scholar] [CrossRef] [Green Version]
  48. Allouche, O.; Steinitz, O.; Rotem, D.; Rosenfeld, A.; Kadmon, R. Incorporating distance constraints into species distribution models. J. Appl. Ecol. 2008, 45, 599–609. [Google Scholar] [CrossRef]
  49. Reed-Andersen, T.; Bennett, E.M.; Jorgensen, B.S.; Lauster, G.; Lewis, D.B.; Nowacek, D.; Riera, J.L.; Sanderson, B.L.; Stedman, R. Distribution of recreational boating across lakes: Do landscape variables affect recreational use? Freshw. Biol. 2000, 43, 439–448. [Google Scholar] [CrossRef]
  50. Open Street Map. Roads Data for Wisconsin. Downloaded from GeoFabrik Website in June, 2014. 2014. Available online: http://www.geofabrik.de/data/download.html (accessed on 2 October 2020).
  51. Ruesch, A.S.; Menuz, D.R.; Diebel, M.W. 1:24K Wisconsin Hydrography Dataset Creation Toolset; Wisconsin Department of Natural Resources: Madison, WI, USA, 2013.
  52. The Register of Waterbodies; Wisconsin Department of Natural Resources: Madison, WI, USA, 2008.
  53. Bivand, R.; Keitt, T.; Rowlingson, B. rgdal: Bindings for the ’Geospatial’ Data Abstraction Library. R Package Version 0.8–13. 2019. Available online: https://cran.r-project.org/web/packages/rgdal/index.html (accessed on 13 October 2020).
  54. Bivand, R.; Rundel, C. rgeos: Interface to Geometry Engine—Open Source (’GEOS’). R Package Version 0.2–2. 2018. Available online: https://cran.r-project.org/web/packages/rgdal/index.html (accessed on 13 October 2020).
  55. Ryden, K. Environmental Systems Research Institute Mapping. Am. Cartogr. 1987, 14, 261–263. [Google Scholar] [CrossRef]
  56. Hijmans, R.J. Raster: Geographic Data Analysis and Modeling. 2017. Available online: https://rspatial.org/raster/ (accessed on 2 October 2020).
  57. Papes, M.; Vander Zanden, J. Wisconsin Lake Historical Limnological Parameters 2010, 1925–2009. Available online: https://doi.org/10.6073/pasta/66320ff8063706f6b3ee83a0ef3ef439 (accessed on 7 October 2020).
  58. Van Buuren, S.; Groothuis-Oudshoorn, K. mice: Multivariate Imputation by Chained Equations inR. J. Stat. Softw. 2011, 45, 1–67. [Google Scholar] [CrossRef] [Green Version]
  59. Fox, J.; Weisberg, S. An R Companion to Applied Regression, 2nd ed.; SAGE Publications: Thousand Oaks, CA, USA, 2011. [Google Scholar]
  60. Firth, D. Bias reduction of maximum likelihood estimates. Biometrika 1993, 80, 27–38. [Google Scholar] [CrossRef]
  61. Fielding, A.H.; Bell, J.F. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 1997, 24, 38–49. [Google Scholar] [CrossRef]
  62. Hosmer, D.W., Jr.; Lemeshow, S.; Sturdivant, R.X. Applied Logistic Regression, 3rd ed.; Wiley: Hoboken, NJ, USA, 2013. [Google Scholar]
  63. Robin, X.; Turck, N.; Hainard, A.; Tiberti, N.; Lisacek, F.; Sánchez, J.C.; Mueller, M. pROC: An open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinform. 2011, 12, 77. [Google Scholar] [CrossRef] [PubMed]
  64. Tjur, T. Coefficients of Determination in Logistic Regression Models—A New Proposal: The Coefficient of Discrimination. Am. Stat. 2009, 63, 366–372. [Google Scholar] [CrossRef]
  65. Chen, J.; Shiyomi, M.; Hori, Y.; Yamamura, Y. Frequency distribution models for spatial patterns of vegetation abundance. Ecol. Model. 2008, 211, 403–410. [Google Scholar] [CrossRef]
  66. Ferrari, S.L.P.; Cribari-Neto, F. Beta Regression for Modelling Rates and Proportions. J. Appl. Stat. 2004, 31, 799–815. [Google Scholar] [CrossRef]
  67. Simas, A.B.; Barreto-Souza, W.; Rocha, A.V. Improved estimators for a general class of beta regression models. Comput. Stat. Data Anal. 2010, 54, 348–366. [Google Scholar] [CrossRef] [Green Version]
  68. Hunger, M.; Baumert, J.; Holle, R. Analysis of SF-6D Index Data: Is Beta Regression Appropriate? Value Health 2011, 14, 759–767. [Google Scholar] [CrossRef] [Green Version]
  69. Cribari-Neto, F.; Zeileis, A. Beta Regression inR. J. Stat. Softw. 2010, 34, 1–24. [Google Scholar] [CrossRef] [Green Version]
  70. Smithson, M.; Verkuilen, J. A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. Psychol. Methods 2006, 11, 54–71. [Google Scholar] [CrossRef] [Green Version]
  71. Grün, B.; Kosmidis, I.; Zeileis, A. Extended Beta Regression inR: Shaken, Stirred, Mixed, and Partitioned. J. Stat. Softw. 2012, 48, 1–25. [Google Scholar] [CrossRef] [Green Version]
  72. Potts, J.; Elith, J. Comparing species abundance models. Ecol. Model. 2006, 199, 153–163. [Google Scholar] [CrossRef]
  73. Couch, R.; Nelson, E. Myriophyllum spicatum in North America. In Proceedings of the First International Symposium on Watermilfoil (Myriophyllum spicatum) and Related Haloragaciae Species, Vicksburg, MS, USA, 23–24 July 1985; Aquatic Plant Management Society: Vicksburg, MS, USA, 1985. [Google Scholar]
  74. Vestergaard, O.; Sand-Jensen, K. Alkalinity and trophic state regulate aquatic plant distribution in Danish lakes. Aquat. Bot. 2000, 67, 85–107. [Google Scholar] [CrossRef]
  75. Santamaría, L. Why are most aquatic plants widely distributed? Dispersal, clonal growth and small-scale heterogeneity in a stressful environment. Acta Oecol. 2002, 23, 137–154. [Google Scholar] [CrossRef]
  76. Johnstone, I.; Coffey, B.; Howard-Williams, C. The role of recreational boat traffic in interlake dispersal of macrophytes: A New Zealand case study. J. Environ. Manag. 1985, 20, 263–279. [Google Scholar]
  77. Havel, J.E.; Shurin, J.B.; Jones, J.R. Estimating dispersal from patterns of spread: Spatial and local control of lake invasions. Ecology 2002, 83, 3306–3318. [Google Scholar] [CrossRef]
  78. Rech, P.H.; Wetzel, R.G.; Thuy, N. Distribution, production and role of aquatic macrophytes in a southern Michigan marl lake. Freshw. Biol. 1971, 1, 3–21. [Google Scholar] [CrossRef]
  79. GBIF.org. 7 Aug 2020. GBIF Occurrence Download. Available online: https://www.gbif.org/occurrence/search?taxon_key=5361760 (accessed on 7 October 2020).
  80. Carpenter, S.R. Enrichment of Lake Wingra, Wisconsin, by submersed macrophyte decay. Ecology 1980, 61, 1145–1155. [Google Scholar] [CrossRef]
  81. Kovalenko, K.E.; Dibble, E.D.; Slade, J.G. Community effects of invasive macrophyte control: Role of invasive plant abundance and habitat complexity. J. Appl. Ecol. 2010, 47, 318–328. [Google Scholar] [CrossRef]
  82. Madsen, J.D.; Sutherland, J.W.; Bloomfield, J.A.; Eichler, L.W.; Boylen, C.W. The decline of native vegetation under dense Eurasian watermilfoil canopies. J. Aquat. Plant Manag. 1991, 29, 94–99. [Google Scholar]
  83. Wilson, S.J.; Ricciardi, A. Epiphytic macroinvertebrate communities on Eurasian watermilfoil (Myriophyllum spicatum) and native milfoils Myriophyllum sibiricum and Myriophyllum alterniflorum in eastern North America. Can. J. Fish. Aquat. Sci. 2009, 66, 18–30. [Google Scholar] [CrossRef]
  84. Horsch, E.J.; Lewis, D.J. The Effects of Aquatic Invasive Species on Property Values: Evidence from a Quasi-Experiment. Land Econ. 2009, 85, 391–409. [Google Scholar] [CrossRef] [Green Version]
  85. Olden, J.D.; Tamayo, M. Incentivizing the Public to Support Invasive Species Management: Eurasian Milfoil Reduces Lakefront Property Values. PLoS ONE 2014, 9, e110458. [Google Scholar] [CrossRef] [Green Version]
  86. Hutchinson, G.E. The chemical ecology of three species of Myriophyllum (Angionspermae, Haloragaceae). Limnol. Oceanogr. 1970, 15, 1–5. [Google Scholar] [CrossRef]
  87. Verheyen, D.; Van Gaelen, N.; Ronchi, B.; Batelaan, O.; Struyf, E.; Govers, G.; Merckx, R.; Diels, J. Dissolved phosphorus transport from soil to surface water in catchments with different land use. Ambio 2015, 44, 228–240. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  88. Parry, R. Agricultural Phosphorus and Water Quality: A U.S. Environmental Protection Agency Perspective. J. Environ. Qual. 1998, 27, 258–261. [Google Scholar] [CrossRef]
  89. Egertson, C.J.; Kopaska, J.A.; Downing, J.A. A Century of change in macrophyte abundance and composition in response to agricultural eutrophication. Hydrobiology 2004, 524, 145–156. [Google Scholar] [CrossRef]
  90. Smith, V.H. Eutrophication of freshwater and coastal marine ecosystems a global problem. Environ. Sci. Pollut. Res. 2003, 10, 126–139. [Google Scholar] [CrossRef]
  91. Grace, J.B.; Wetzel, R.G. The production biology of Eurasian watermilfoil (Myriophyllum spicatum L.): A Review. J. Aquat. Plant Manag. 1978, 16, 1–11. [Google Scholar]
  92. Carpenter, S.R. The decline of Myriophyllum spicatum in a eutrophic Wisconsin lake. Can. J. Bot. 1980, 58, 527–535. [Google Scholar] [CrossRef]
  93. Newman, R.M.; Biesboer, D.D. A decline of Eurasian watermilfoil in Minnesota associated with the milfoil weevil, Euhrychiopsis lecontei. J. Aquat. Plant Manag. 2000, 38, 105–111. [Google Scholar]
  94. Hansen, G.J.A.; Zanden, M.J.V.; Blum, M.J.; Clayton, M.K.; Hain, E.F.; Hauxwell, J.; Izzo, M.; Kornis, M.S.; McIntyre, P.B.; Mikulyuk, A.; et al. Commonly rare and rarely common: Comparing population abundance of invasive and native aquatic species. PLoS ONE 2013, 8, e77415. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  95. Vander Zanden, M.J.; Olden, J.D. A management framework for preventing the secondary spread of aquatic invasive species. Can. J. Fish. Aquat. Sci. 2008, 65, 1512–1522. [Google Scholar] [CrossRef]
Figure 1. Eurasian watermilfoil (EWM) occurrence predictions for N = 657 surveyed lakes. Maps depict true presences (a), true absences (b), false absences (c), and false presences (d). Binomial response generalized linear models (GLMs) fit with Firth’s bias reduction method [24].
Figure 1. Eurasian watermilfoil (EWM) occurrence predictions for N = 657 surveyed lakes. Maps depict true presences (a), true absences (b), false absences (c), and false presences (d). Binomial response generalized linear models (GLMs) fit with Firth’s bias reduction method [24].
Diversity 12 00394 g001
Figure 2. Eurasian watermilfoil abundance on N = 657 surveyed lakes. Maps depict abundance observed during aquatic plant field surveys (a) and predicted using beta regression (b).
Figure 2. Eurasian watermilfoil abundance on N = 657 surveyed lakes. Maps depict abundance observed during aquatic plant field surveys (a) and predicted using beta regression (b).
Diversity 12 00394 g002
Figure 3. Occurrence and abundance models developed on N = 657 and extended to 9825 lakes over 1 ha in surface area. Eurasian watermilfoil occurrence predictions (a) generated using a bias-reduced GLM. The probability threshold distinguishing likely presences from absences was set at 0.246 to minimize the chance of false absences. Eurasian watermilfoil abundance predictions (b) generated using a beta regression model.
Figure 3. Occurrence and abundance models developed on N = 657 and extended to 9825 lakes over 1 ha in surface area. Eurasian watermilfoil occurrence predictions (a) generated using a bias-reduced GLM. The probability threshold distinguishing likely presences from absences was set at 0.246 to minimize the chance of false absences. Eurasian watermilfoil abundance predictions (b) generated using a beta regression model.
Diversity 12 00394 g003
Figure 4. Management priority tiers based on modelled risk of EWM occurrence and abundance following the prioritization matrix in Table 4. Lakes shown are those not known to contain EWM currently, but where occurrence models predicted it is likely to occur. Tier 1 (orange) lakes are the most vulnerable. They are likely to support EWM at high abundance and have top priority for prevention and management. Tier 2 lakes include those with moderate risk for EWM invasion and abundance. Tier 3 lakes are less likely to have EWM populations and are also unlikely to support abundant populations.
Figure 4. Management priority tiers based on modelled risk of EWM occurrence and abundance following the prioritization matrix in Table 4. Lakes shown are those not known to contain EWM currently, but where occurrence models predicted it is likely to occur. Tier 1 (orange) lakes are the most vulnerable. They are likely to support EWM at high abundance and have top priority for prevention and management. Tier 2 lakes include those with moderate risk for EWM invasion and abundance. Tier 3 lakes are less likely to have EWM populations and are also unlikely to support abundant populations.
Diversity 12 00394 g004
Table 1. Estimated coefficients for Eurasian watermilfoil (EWM) occurrence model developed using data from 657 plant survey lakes. Coefficients expressed as odds ratios calculated for centered and scaled predictors, profile confidence intervals in parentheses. Negative relationships are indicated by odds ratios less than 1.
Table 1. Estimated coefficients for Eurasian watermilfoil (EWM) occurrence model developed using data from 657 plant survey lakes. Coefficients expressed as odds ratios calculated for centered and scaled predictors, profile confidence intervals in parentheses. Negative relationships are indicated by odds ratios less than 1.
PredictorCoefficient
Intercept0.13 *** (0.07−0.21)
Road density (log (m/ha +1))1.93 *** (1.42−2.74)
Surface area (log ha)1.72 *** (1.33−2.31)
Maximum air temp. (°C × 10)1.69 ** (1.23−2.40)
Maximum depth (log m +1)1.55 ** (1.20−2.06)
Conductivity (log μS/cm)1.47 (0.72−3.20)
Alkalinity (log mg CaCO3 +1)1.44 (0.67−3.03)
Soil erodibility (kwfact)1.17 (0.93−1.49)
Watershed urban ( % )1.07 (0.77−1.53)
pH1.06 (0.80−1.41)
Secchi depth (log m +1)0.85 (0.62−1.16)
Watershed agriculture ( % )0.81 (0.58−1.10)
CaO ( % )0.74 ** (0.57−0.92)
Annual temp. range (°C × 10)0.64 * (0.42−0.92)
Mean distance source (log m)0.61 *** (0.45−0.82)
* p < 0.05, ** p < 0.01, *** p < 0.001.
Table 2. Cross-validated confusion matrix relating observed to predicted EWM occurrence in the 657 surveyed lakes. We set the threshold for classifying the predicted probability as presence at 0.246, a value that minimizes false absences to less than 10% of all predicted absences.
Table 2. Cross-validated confusion matrix relating observed to predicted EWM occurrence in the 657 surveyed lakes. We set the threshold for classifying the predicted probability as presence at 0.246, a value that minimizes false absences to less than 10% of all predicted absences.
Observed
AbsentPresent
PredictedAbsent10011
Present169377
Table 3. Coefficients estimated using a beta regression model for abundance data on 657 surveyed lakes. Coefficients and standard errors for mean (logit link) and precision (log-link) submodels describe patterns and variability in the EWM abundance response. Negative linear coefficients for the mean submodel indicate factors that predict lower EWM abundance, negative quadratic coefficients describe concave-down relationships, and positive quadratic coefficients are concave-up. For a fixed mean estimate in the precision submodel, a larger precision coefficient indicates lower variance in the response.
Table 3. Coefficients estimated using a beta regression model for abundance data on 657 surveyed lakes. Coefficients and standard errors for mean (logit link) and precision (log-link) submodels describe patterns and variability in the EWM abundance response. Negative linear coefficients for the mean submodel indicate factors that predict lower EWM abundance, negative quadratic coefficients describe concave-down relationships, and positive quadratic coefficients are concave-up. For a fixed mean estimate in the precision submodel, a larger precision coefficient indicates lower variance in the response.
Mean SubmodelPrecision Submodel
LinearQuadraticLinearQuadratic
PredictorsEstimateSEEstimateSEEstimateSEEstimateSE
Intercept−3.71 ***0.28 2.11 ***0.31
Conductivity (log μS/cm)1.27 ***0.24−0.52 ***0.11−1.74 ***0.310.61 ***0.14
Road density (log (m/ha +1))0.34 *0.14−0.060.06−0.240.170.120.07
Alkalinity (log mg CaCO3 +1)0.310.23−0.010.11−0.73 *0.300.41 **0.14
Maximum air temp. (°C x 10)0.28 *0.11−0.110.100.000.120.150.11
Mean distance source (log m)0.220.160.19 *0.09−0.300.18−0.23 *0.10
Watershed agriculture ( % )0.100.130.030.08−0.31 *0.150.020.09
Watershed urban ( % )0.100.140.020.04−0.230.160.000.04
Maximum depth (log m +1)0.070.11−0.100.060.090.130.110.07
Soil erodibility (kwfact)0.040.08−0.23 **0.09−0.030.090.37 ***0.10
Surface area (log ha)0.020.190.020.050.090.22−0.020.06
Secchi depth (log m +1)−0.080.10−0.090.060.000.110.030.06
pH−0.110.110.010.050.270.15−0.060.06
CaO ( % )−0.23 *0.100.060.050.220.12−0.060.05
Annual temp. range (°C x 10)−0.54 *0.21−0.050.070.420.230.050.07
Log-likelihood1856
Df58
Pseudo R20.24
* p < 0.05, ** p < 0.01, *** p < 0.001.
Table 4. Three-tiered management priority matrix for uninvaded lakes that are vulnerable to EWM invasion (N = 2357 lakes where probability of occurrence is above the threshold of 0.246). Tier 1 lakes are most likely to have EWM and support abundant populations; they are assigned the top management priority. Tier 2 lakes have moderate risk for occurrence and abundance or a high risk for either occurrence or abundance. Tier 3 lakes comprise the lowest-priority group, i.e., they are still vulnerable to invasion but have a lower probability of occurrence and are less likely to support abundant populations.
Table 4. Three-tiered management priority matrix for uninvaded lakes that are vulnerable to EWM invasion (N = 2357 lakes where probability of occurrence is above the threshold of 0.246). Tier 1 lakes are most likely to have EWM and support abundant populations; they are assigned the top management priority. Tier 2 lakes have moderate risk for occurrence and abundance or a high risk for either occurrence or abundance. Tier 3 lakes comprise the lowest-priority group, i.e., they are still vulnerable to invasion but have a lower probability of occurrence and are less likely to support abundant populations.
Abundance
High
(0.36–1)
Med.
(0.18–0.36)
Low
(0–0.18)
PresenceHigh (0.75–1)28188227
Med. (0.50–0.75)170644
Low (0.25–0.50)0281171
Prevention and Control PriorityTier 1Tier 2Tier 3

Share and Cite

MDPI and ACS Style

Mikulyuk, A.; Hein, C.L.; Van Egeren, S.; Kujawa, E.R.; Vander Zanden, M.J. Prioritizing Management of Non-Native Eurasian Watermilfoil Using Species Occurrence and Abundance Predictions. Diversity 2020, 12, 394. https://doi.org/10.3390/d12100394

AMA Style

Mikulyuk A, Hein CL, Van Egeren S, Kujawa ER, Vander Zanden MJ. Prioritizing Management of Non-Native Eurasian Watermilfoil Using Species Occurrence and Abundance Predictions. Diversity. 2020; 12(10):394. https://doi.org/10.3390/d12100394

Chicago/Turabian Style

Mikulyuk, Alison, Catherine L. Hein, Scott Van Egeren, Ellen Ruth Kujawa, and M. Jake Vander Zanden. 2020. "Prioritizing Management of Non-Native Eurasian Watermilfoil Using Species Occurrence and Abundance Predictions" Diversity 12, no. 10: 394. https://doi.org/10.3390/d12100394

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop