Predicting suitable habitat for dreissenid mussel invasion in Texas based on climatic and lake physical characteristics

Eurasian zebra and quagga mussels were likely introduced to the Laurentian Great Lakes via ballast water release in the 1980s, and their range has since expanded across the US, including some of their southernmost occurrences in Texas. Their spread into the state has resulted in a need to revise previous delimitations of suitable dreissenid habitat. We therefore assessed invasion risk in Texas by 1) predicting distribution of suitable habitat of zebra and quagga mussels using Maxent species distribution models based upon global occurrence and climate data; and 2) refining lake-specific predictions via collection and analysis of physicochemical data. Maxent models predicted a lack of suitable habitat for quagga mussels within Texas. However, models did predict the presence of suitable zebra mussel habitat, with hotspots of suitable habitat occurring along the Red and Sabine Rivers of north and east Texas, as well as patches of suitable habitat in central Texas between the Colorado and Brazos Rivers and extending inland along the Gulf Coast. Although predicted suitable habitat extended further west than in previous models, most of the Texas panhandle, west Texas extending toward El Paso, and the Rio Grande valley were predicted to provide poor zebra mussel habitat suitability. Collection of physicochemical data (i.e., dissolved oxygen, pH, specific conductance, and temperature on-site as well as laboratory analysis for Ca, N, and P) from zebra mussel-invaded lakes and a subset of uninvaded but high-risk lakes of North and Central Texas, did not refine model predictions because there was no apparent distinction between invaded and uninvaded lakes. Overall, we demonstrated that while quagga mussels do not appear to represent an invasive threat in Texas, abundant suitable habitat for continuing zebra mussel invasion exists within the state. The threat of continued expansion of this poster-child for negative invasive species impacts warrants further prevention efforts, management, and research.


Introduction
Zebra mussels (Dreissena polymorpha Pallas, 1771) were likely introduced from Ponto-Caspian Eurasia to the Laurentian Great Lakes inadvertently via ballast water release around 1986 (Hebert et al. 1989), and closely related quagga mussels (Dreissena bugensis Andrusov, 1897) followed within the next several years (Mills et al. 1996). Both dreissenid mussel species spread rapidly throughout the Great Lakes basin due to high fecundity (Keller et al. 2007) as well as considerable natural and anthropogenic dispersal potential of their free-floating larval life stages (Bossenbroek et al. 2001;Sieracki et al. 2014). The Great Lakes have since served as a beachhead for dreissenid mussel invasion throughout much of North America, including the Hudson and Mississippi River basins (Strayer et al. 1996;Cope et al. 1997). Today, both species have extended their ranges throughout the United States. The initial incursion of zebra mussels into Texas occurred in Lake Texoma on the Texas-Oklahoma border around 2009 (TPWD 2017). Recent in situ analyses in Lake Texoma suggested that environmental conditions (especially relatively high water temperatures) encountered as dreissenid mussels expand their range southward could increase growth rates, decrease time to maturity, and promote invasion (Churchill et al. 2017). Given recent spread of zebra mussels further into the state in reservoirs such as Belton Lake and Lake Waco on the Brazos River, delimiting suitable dreissenid habitat within Texas represents a critical need for managing the southward expansion of these invasive mussels.
Management concern about dreissenid mussels stems from their ability to act as ecosystem engineers, manipulating physical habitat by forming dense colonies attached to hard substrates (including native mussels) as well as transforming turbid, eutrophic systems into clear waters through voracious filter feeding behavior, which can contribute to increases in aquatic vegetation and profound shifts in native communities (reviewed in Nakano and Strayer 2014). Additional consequences of zebra mussel invasion include expense of physical removal of colonies from recreational equipment, power plants, municipal water facilities, dams, and other human infrastructure, as well as the application of pesticides to prevent or slow their reintroduction and recolonization of such structures (Aldridge et al. 2006). Overall, the economic burden of zebra mussel invasion approaches hundreds of millions of dollars annually (Strayer et al. 1996;Caraco et al. 1997). Surveillance and rapid response efforts such as the 100 th Meridian Initiative (http://www.100thmeridian.org) have emphasized the ongoing need for large scale management of the dreissenid invasion.
High dissolved calcium requirements for dreissenid shell development represent a potentially critical constraining factor for dreissenid mussel colonization (Koutnik and Padilla 1994). Indeed, a risk assessment conducted by Whittier et al. (2008) defined risk based solely on calcium concentrations, distinguishing risk categories of "very low" (< 12 mg/L), "low" (12-20 mg/L), "moderate" (20-28 mg/L), and "high" (> 28 mg/L). Within Texas, the Brazos River basin and its neighboring basins to the southwest (Colorado River and Rio Grande) have naturally high levels of dissolved calcium, especially over the Permian plateau (VanLandeghem et al. 2012;Israël et al. 2014;Sharma et al. 2014). More recent assessment of quagga mussel survival, growth, and reproductive potential in waters of the western United States by Davis et al. (2015) noted that a single indicator such as calcium may oversimplify risk considerations and advocated considering additional factors and scales (e.g., whole lake vs. microhabitat). Therefore, we conducted research to refine understanding of potential dreissenid mussel distribution in Texas at multiple scales through both a state-wide distribution modeling approach using globally available climatic predictors as well as physicochemical data collection at single-lake scales.
Overall, the goal of our research was to assess potential distribution of the dreissenid mussel invasion in Texas. Specific objectives of our research included: 1) predicting the general distribution of suitable habitat for zebra and quagga mussels in Texas using Maxent models; and 2) refining lakespecific zebra mussel (i.e., the only dreissenid known to occur in Texas) predictions via collection of physicochemical data.

Materials and methods
Objective 1: Predict the general distribution of suitable dreissenid mussel habitat in Texas using Maxent models Species distribution models for zebra and quagga mussels were created with Maxent version 3.4.1, a machine learning tool that compares the probability distributions of species presence and local environmental data to create a model that can be projected in geographic space . In other words, Maxent combines species occurrence data and environmental covariates to produce a heat map that can be interpreted as a visual representation of habitat suitability. Among its strengths, Maxent has been praised for its strong performance compared to other species distribution modeling methods (Elith et al. 2006), and it is amenable to datasets consisting exclusively of presence-only data (Elith et al. 2011;Wittmann et al. 2016).
Separate models were developed for zebra and quagga mussels. To collect data for Maxent modeling, we accessed global zebra and quagga mussel occurrence data through the Global Biodiversity Information Facility online database (http://www.gbif.org), and we supplemented this information with known zebra mussel occurrences in Texas (N = 11 from Figure 1). No records for quagga mussel occur within Texas. We used global occurrence data (rather than limiting our data to Texas only) to ensure that the most complete representation of the zebra and quagga mussel niches were included in model development. Overall, we compiled 13,297 total global occurrences of zebra mussels and 1,069 global quagga mussel occurrences. To reduce bias that may be generated by uneven sampling effort, we rarified occurrence data before model implementation by converting the occurrence points to a raster file with the same cell size as our environmental data (10 arcminute, approximately 340 km 2 ), then Figure 1. Physicochemical data survey lakes. Sites categorized by TPWD (at the time of this study in 2016) as "infested" (the water body has an established, reproducing population) or "positive" (zebra mussels or their larvae have been detected on more than one occasion despite lack of evidence of a fully established, reproducing population) are indicated by red triangles and included: Lakes Austin, Belton, Bridgeport, Dean Gilbert, Lavon, Lewisville, Ray Roberts, Stillhouse Hollow, Texoma, Travis, and Waco. Sites categorized by TPWD as zebra mussel "negative" are indicated by green circles and included: Lakes Aquilla, Buchanan, Georgetown, Granbury, Granger, Hubbard Creek, Inks, Lady Bird, LBJ, Limestone, Marble Falls, Palo Pinto, Pflugerville, Possum Kingdom, Proctor, and Whitney. back to a points file, resulting in a maximum of one point per cell (McDowell et al. 2014). Data conversion and all further mentioned visualizations were performed in ArcGIS 10.2.2 (Environmental Systems Research Institute, Redlands, California, USA). As a result of rarefication, the working zebra mussel occurrence dataset included 2,080 points, and the working quagga mussel occurrence dataset included 318 points.
Environmental data used in the model included the 19 Bioclim layers (http://www.worldclim.org/bioclim), biologically meaningful and globally continuous layers generated from annual trends in temperature and precipitation (Hijmans et al. 2005), with 10-acrminute (approximately 340 km 2 ) resolution. Other predictor variables (e.g., calcium) have obvious appeal for modeling dreissenid mussel habitat suitability, but the unavailability of such variables in the format of a continuous global data layer precluded their incorporation into our modeling effort. For the environmental variables that were available, we assessed collinearity (Dormann et al. 2013) for zebra mussels and quagga mussels separately, retaining only those predictor variables with Pearson's correlation | r | ≤ 0.7 in all pairwise comparisons (calculated using R version 3.5.1). Thus, the predictor variables used in zebra mussel models included the Bioclim layers representing annual mean temperature, isothermality, temperature annual range, mean temperature of wettest quarter, mean temperature of warmest quarter, precipitation of wettest month, precipitation of driest month, and precipitation seasonality. Quagga mussel models included mean diurnal range, temperature seasonality, maximum temperature of warmest month, mean temperature of wettest quarter, mean temperature of coldest quarter, precipitation seasonality, and precipitation of wettest quarter.
Because the purpose of our model development was to generate predictions rather than evaluate overall Maxent performance, we generally opted for default Maxent software settings (Phillips and Dudík 2008). However, we did increase the maximum allowable model iterations to 5,000 based on pilot model runs in which models didn't appear to converge on optimal solutions within the default 500 iterations. Overall, we produced 100 replicate models for each mussel species, each trained with a randomly-selected 80% of rarefied occurrence data and evaluated with the remaining 20%. The reported results represent the average of the 100 models produced for each mussel species. We interpreted the Maxent logistic output as the probability of mussel habitat suitability found around the globe. In addition to visual inspection and description of the Maxent output, we assessed model performance using area under the receiver operating characteristic curve (AUC), where AUC = 0.5 indicates the model predicts outcomes no better than random, and AUC ≥ 0.7 indicates strong predictive power (Hosmer and Lemeshow 2000). Because our models were global in scale, multivariate environmental similarity surfaces (MESS; Elith et al. 2010) and mobility-oriented parity (MOP; Owens et al. 2013) outputs generated automatically during Maxent implementation were expected to indicate no environmental extrapolation. Model uncertainty resulting from our 80% subsampling routine was assessed by examining standard deviations of 100 replicate models for each species.
Objective 2: Refine lake-specific predictions based on collection of physicochemical data.
Quagga mussels are not known to occur in Texas, but to build upon the predictions of suitable zebra mussel habitat produced through Maxent modeling associated with Objective 1, we collected physicochemical data from invaded lakes and a subset of identified high-risk lakes of North and Central Texas (N = 27; Figure 1). Because movement of recreational boats and other anthropogenic vectors represents a primary means of zebra mussel dispersal (Bossenbroek et al. 2001), we prioritized lakes based on number of public boat launches as an index of human use, as well as proximity to lakes already known to be invaded. We further refined our list based on consultations with management partners (Monica McGarrity, Texas Parks and Wildlife Department, personal communication), and data availability from previous physicochemical data compilations (VanLandeghem et al. 2012;Dawson et al. 2015).
Surveyed lakes included eight lakes classified at the time of this study as "infested" with zebra mussels (i.e., water body has an established, reproducing population) by Texas Parks and Wildlife Department (TPWD). These lakes were: Belton, Bridgeport, Dean Gilbert, Lewisville, Ray Roberts, Stillhouse Hollow, Texoma, and Travis. We also surveyed Lakes Austin, Lavon, and Waco, where zebra mussels or their larvae had been detected on more than one occasion despite lack of evidence of a fully established, reproducing population (termed "positive" by TPWD). Finally, we surveyed a suite of sixteen negative sites (where zebra mussels had not previously been reported at the time of this study) across the Brazos and Colorado River basins, including Lakes Aquilla, Buchanan, Georgetown, Granbury, Granger, Hubbard Creek, Inks, Lady Bird, Lyndon B. Johnson, Limestone, Marble Falls, Palo Pinto, Pflugerville, Possum Kingdom, Proctor, and Whitney.
Surveys occurred October 12-16 and 20-21, 2016. At one to three public access points (depending upon availability) at each lake, we recorded water chemistry conditions including dissolved oxygen (DO), pH, specific conductance, and temperature using a Hach HQd water chemistry probe or YSI 556 multiparameter probe. The Hach instrument reported actual conductivity, which was converted to specific conductance according to the formula, where SC is specific conductance at 25 °C, AC is reported conductivity, T is the sample temperature, and r (= 0.0191) is a temperature correction coefficient (Miller et al. 1988).
Additionally, at each sampled location within a lake, we collected two 500-mL water samples for laboratory determination of total nitrogen (total N), total Kjeldahl nitrogen (TKN), and total phosphorus (total P) using a Hach DR3900 spectrophotometer, and calcium-hardness (measured as CaCO 3 ) using a Hach Digital Titrator, all according to manufacturer instructions. We estimated inorganic N in each sample by subtracting TKN from total N (and assuming that ammonia-N was negligible).
To manage the threat of transporting invasive species between invaded and uninvaded inland waters, all research equipment was soaked in 10% bleach solution for a minimum of 10 minutes between sites. Due to the sensitivity of water chemistry meters, they could not be soaked in bleach, but were instead sprayed with 10% bleach solution, then immediately rinsed with clean water.
Statistical analyses were conducted with Dell Statistica, version 13 (Dell Inc., Tulsa, Oklahoma, USA). Replicate samples within each lake were averaged for each water quality variable, and the average value was used for analysis. We applied Kolmogorov-Smirnov tests to examine normality of distributions of each measured factor, and log transformations were applied when necessary. Pearson correlation analysis was then applied to identify and eliminate redundant variables (| r | ≥ 0.7). After a priori elimination of correlation between predictor variables, principal component analysis (PCA) was used to determine whether any of the remaining variables were related to zebra mussel presence/absence. Specifically, biplots of the appropriate principal components were visually inspected to assess patterns of water quality distributions across our survey lakes and identify differences between lakes with and without zebra mussels. In addition, MANOVA was applied to case scores for the appropriate principal components using the presence or absence of zebra mussels as an independent (grouping) factor. The pH sensor in one of the probes malfunctioned during sampling at two sites (Granbury, Proctor) so these data were not collected. To be able to include these sites in the multivariate analysis, their pH values were imputed using k-nearest neighbors analysis (k = 3) based on regression against temperature, which had the highest correlation coefficient with pH across all lakes (Pearson r = −0.50).

Additional methodological considerations
We summarize our approach as a combination of predictions at a broad scale (i.e., species distribution modeling) and fine scale (i.e., lake-specific chemistry measurement), and we don't presume that these two scales of analysis would relate to one another. Therefore, we did not perform analysis to test for a statistical relationship between the two scales.
We also note that all analyses completed in this study reflect our understanding of dreissenid mussel distributions in Texas at the time of data collection in 2016, but the status of some lakes changed (i.e., several sites we considered zebra mussel negative in 2016 have since become positive) over the time period when this manuscript was written and peer reviewed. For a current map of zebra mussel distribution in Texas, see https://tpwd.texas.gov/huntwild/wild/species/exotic/zebramusselmap.phtml. While we consider changes since 2016 and their relevance in the Discussion, we did not update our analysis in response to these changes because, particularly for lake-specific physicochemical measurements, we could not rule out the possibility that changes in lake conditions influenced the change in zebra mussel status.

Objective 1: Predict the general distribution of suitable dreissenid mussel habitat in Texas using Maxent models
We produced a Maxent model using the Bioclim environmental layers and global zebra mussel occurrence data to predict the extent of suitable zebra mussel habitat in Texas. The result of this analysis was a heat map in which shading indicates the logistic output of the Maxent model; warmer colors were interpreted as relatively suitable habitat, and cooler colors were interpreted as less suitable habitat (Figure 2). The average area under the receiver operating characteristic curve of 0.918 indicated that the model accurately distinguished global zebra mussel occurrences and provides confidence in our ability to make predictions in uninvaded sites in Texas. The most important environmental layers (i.e., > 10% contribution to Maxent model predictions) included precipitation of the driest month, annual mean temperature, and isothermality. Predicted climatic "hotspots" of suitable zebra mussel habitat were concentrated along the Red and Sabine Rivers of the northern and eastern Texas borders and extended into central Texas between the Colorado and Brazos Rivers as well as along the Gulf Coast. Most of the Texas panhandle, west Texas extending toward El Paso, and the Rio Grande valley were predicted to provide poor zebra mussel habitat suitability.
We also produced a Maxent model to predict the extent of suitable quagga mussel habitat in Texas. The average area under the receiver operating characteristic curve of 0.970 indicated that the model accurately distinguished global quagga mussel occurrences. The most important environmental layers included precipitation seasonality, mean temperature of the coldest quarter, precipitation of the wettest quarter, and temperature seasonality. No suitable quagga mussel habitat was identified in Texas (i.e., all Maxent scores < 0.5; Figure 3).
Because we produced Maxent models using global zebra and quagga mussel occurrence data, multivariate environmental similarity surfaces (MESS) and mobility-oriented parity (MOP) outputs confirmed that no environmental extrapolation occurred, as expected, indicating little intrinsic model uncertainty due to transferability. Furthermore, the 80% subsampling routine employed during the production of 100 replicate models for each species resulted in limited variation between models; low standard deviations among quagga mussel models (maximum standard deviation < 0.09 on the logistic scale) and zebra mussel models (maximum standard deviation < 0.03) suggested low levels of uncertainty in overall model predictions for both species.

Objective 2: Refine lake-specific predictions based on collection of physicochemical data from identified high-risk lakes of North and Central Texas
We collected physicochemical data from zebra mussel invaded and uninvaded lakes in North and Central Texas to analyze habitat suitability on a fine scale and build upon the distribution model produced in Objective 1 (Table 1). Of the 9 water quality variables considered, 3 were normally distributed (temperature, pH, DO). Other variables (specific conductance, Ca, total N, inorganic N, TKN, and total P) were logtransformed to improve normality. Correlation analysis showed that total N was highly correlated with TKN (r = 0.98) and Ca was highly correlated with specific conductance (r = 0.91); thus, TKN and specific conductance were not used in further analyses. Results of PCA indicated that the first three principal components had eigenvalues > 1 and together accounted for 67% of data variability. Visual inspection of biplots of these components did not reveal any clear separation between zebra mussel positive and negative lakes (Figure 4). The multivariate test of significance on the first three components using presence or absence of zebra mussels as grouping factor also failed to show differences between positive and negative lakes (Wilks λ = 0.936; F 3,23 = 0.524, p = 0.67). Analysis of water quality variables in 27 study lakes. Variables that predominated in each component (|factor loading| ≥ 0.50) are shown on the appropriate axes. Individual lake data are represented by symbols, with open circles representing lakes without previously reported incidences of zebra mussels (absent, 16 lakes), and solid circles those known to harbor the invasive species (present, 11 lakes) at the time of sampling (October 2016). No separation between the two lake groups is evident in either of the biplots. Ca, calcium, N, nitrogen; P, phosphorous.

Discussion
We characterized the potential habitat of nonindigenous invasive dreissenid mussels in Texas. Using global occurrence data and continuous environmental layers available at a global scale (i.e., temperature-and precipitation-based variables), we produced a Maxent species distribution model that predicts that much of north and central Texas contains suitable habitat for zebra mussels (Figure 2). Results of water quality analysis of lakes in this region failed to discriminate between zebra mussel positive and negative lakes at the time of sampling. Together, these results suggest that much of Texas is at risk for invasion by zebra mussels, with potential invasion hotspots occurring around the Red, Sabine, Neches, Trinity, Brazos, and Colorado Rivers and their surrounding reservoirs. The model suggested that west Texas (i.e., Rio Grande, Pecos, upper Brazos and upper Colorado Rivers) has a low risk of zebra mussel establishment based on our estimate of habitat suitability. The model predicted that annual mean temperature and isothermality are important drivers of zebra mussel habitat suitability; therefore, west Texas may experience temperature fluctuations over the course of the year that are too intense to support zebra mussel populations. Additionally, the model indicated that precipitation of the driest month represents an important driver of zebra mussel habitat, which may simply reflect the fact that zebra mussels are obligately aquatic organisms; increased precipitation likely corresponds to more availability of aquatic habitat, and west Texas may appear climatically to have less available habitat. Nevertheless, the aquatic habitats that are available in west Texas are characterized by relatively high calcium concentrations (VanLandeghem et al. 2012;Israël et al. 2014;Sharma et al. 2014), a known predictor of zebra mussel success (Cohen 2005). Indeed, all but two of the calcium measurements (derived from calcium carbonate equivalents) recorded in this study exceeded the "high risk" categorization (i.e., > 28 mg/L) of Whittier et al. (2008). Furthermore, because zebra mussels can withstand salinities up to 10 ppt (Ludyanskiy et al. 1993), the relatively high salinities of west Texas watersheds are not likely to serve as barriers against expansion. On the contrary, these conditions seem favorable for further zebra mussel expansion, ironically as increasing salinity appears to be reducing suitability for native mussels (e.g., Hart et al. 2019). The susceptibility of the western regions of the state to zebra mussel invasion deserves further study.
Our distribution model predictions generally align with previous efforts to predict the potential range of zebra mussels. Drake and Bossenbroek (2004) developed predictions of zebra mussel habitat using a different modeling approach, a genetic algorithm for rule-set production (GARP; Stockwell and Peters 1999), which could contribute to differences in predictions. Drake and Bossenbroek (2004) notably developed separate models using climate data alone or combined with underlying geologic data as a proxy for water chemistry. In comparison with both models, our prediction of suitable zebra mussel habitat extends further northwest, nearing the Texas panhandle. Gallardo et al. (2013) presented maps of predictions of zebra mussel habitat based on Maxent models developed with native-range occurrence data only, European range data only, or North American data only, but not a model developed using all available data. Of their presented maps, the North American model predicted the largest amount of suitable habitat in Texas, although thresholding applied in their application led Gallardo et al. (2013) to predict a lack of suitable habitat in central and most of coastal Texas. The Maxent-based predictions of Quinn et al. (2014) were developed with global zebra mussel occurrence data but also appeared to fail to predict most of central Texas as suitable zebra mussel habitat due to applied thresholding. Thus, overall, our predictions of suitable zebra mussel habitat in Texas extend further than previously developed models, which is likely explained by the accumulation of more occurrence data as the invasion has continued (e.g., Václavík and Meentemeyer 2012) but may also be the result of rapid evolution and niche expansion during the invasion process (e.g., Cox 2004).
Texas does not appear to contain climatically suitable habitat for quagga mussel (Figure 3). Our predictions notably differ with the Maxent models of Quinn et al. (2014), which predicted limited suitable habitat near the Texas-Louisiana border and southern Texas along the Rio Grande. The most important environmental layers in our Maxent model included temperature seasonality and mean temperature of the coldest quarter, which could suggest that temperature fluctuations prevent Texas from providing suitable habitat for quagga mussels. As with zebra mussels, the quagga mussel model emphasized the importance of precipitation factors, including precipitation of the wettest month and precipitation seasonality. This could again suggest that based on climate data alone, the Texas landscape is not expected to have surface water habitat capable of supporting populations of the obligately aquatic quagga mussel. However, the use of natural indicators of water availability as predictors within species distribution modeling efforts focusing on aquatic organisms has been questioned previously (Barnes et al. 2014) because this may result in misleading model outputs; although the climatic variables of Texas do not naturally promote the maintenance of large, lentic ecosystems, the presence of anthropogenic reservoirs, canals, and even small recreational and aesthetic ponds could provide suitable habitats that the model cannot predict. Furthermore, the fact that all the large surface water systems in Texas are anthropogenic reservoirs rather than natural lakes may result in hydrologic dynamics (e.g., responses to precipitation events) that differ from the patterns and processes typical of natural lakes, including the Great Lakes ecosystem where zebra and quagga mussels were first introduced in North America. The effect of these differences on invasion dynamics and species distribution model performance warrant future consideration.
It is critical to note the caveat that distribution models assume that modeled organisms are at equilibrium with their environment (i.e., not demonstrating range expansion or contraction). Given recent, rapid spread of zebra mussels in Texas, this assumption has likely been violated. Previous modeling experiments by Václavík and Meentemeyer (2012) demonstrated that the predictions of distribution models may change as an invasion progresses and new sites are colonized. Therefore, while our model represents an accurate (AUC = 0.918) prediction of the distribution of zebra mussel habitat based on current knowledge of zebra mussel occurrences, predictions could change as zebra mussel invasion continues.
Similarly, quagga mussel habitat predictions may yet expand into the state. The lack of equilibrium in both zebra mussel and quagga mussel distributions may also help explain why habitat suitability maximums below 1 were predicted (i.e., 0.71 for zebra mussels; 0.87 for quagga mussels). Non-equilibrium conditions could inhibit the ability of the Maxent models to confidently identify habitat of maximum suitability. This phenomenon could also indicate that the two mussel species act as habitat generalists with wide fundamental niches. Distribution models may provide some insight into the relative population success that zebra and quagga mussels would experience in various parts of the state if introduced (e.g., Wittmann et al. 2016), but further research into the relationship between dreissenid mussel performance and distribution model results would strengthen this claim.
Overall, we have demonstrated that suitable habitat for continuing zebra mussel invasion exists within Texas. Our predictions are further supported by recent zebra mussel expansions within the state. We developed our distribution models based on zebra mussel distributions in 2016; since then, eight sites which we classified as negative have become classified as "infested" by TPWD: Canyon, Eagle Mountain, Georgetown, Lady Bird, Livingston, Lyndon B. Johnson, Pflugerville, and Randell. Additional "positive" designations have occurred for multiple lakes, including, Dunlap, Fishing Hole, Fork, Granger, Grapevine, McQueeney, Placid, Ray Hubbard, Richland Chambers, Walter E. Long, and Worth. Notably, several lakes considered negative at the time of our physicochemical data collection are now infested or positive: Georgetown, Granger, Lady Bird, Lyndon B. Johnson, and Pflugerville. We did not update our analysis in response to changes in zebra mussel status because, particularly for lake-specific physicochemical measurements, we could not rule out the possibility that changes in lake conditions influenced the change in zebra mussel status. Future studies may consider repeating our analysis to see if our findings remain consistent. For now, we conclude that the spread of zebra mussels to these sites supports our conclusions that much of Texas is at risk.
Continued expansion of zebra mussel, a poster child for negative invasive species impacts, represents a considerable management challenge in Texas. As previously mentioned, most water bodies in Texas are anthropogenic reservoirs rather than natural lakes, and the tendency of reservoirs to increase invasion rates relative to natural waters (Johnson et al. 2008) exacerbates this challenge. Furthermore, zebra mussels have promoted shifts in phytoplankton communities through selective filter feeding in the Laurentian Great Lakes (Vanderploeg et al. 2001) and inland reservoirs (Yu and Culver 2000), and a similar fate may await the waters of Texas. Zebra mussel invasion also poses an acute threat to unionid mussels through competition for resources and direct interference through colonization on unionid shells (Schloesser et al. 1996), and these native mussels represent a threatened group in Texas like many other parts of 7 7 North America (e.g., Ford et al. 2014). Combined with potential impacts to other flora and fauna (Nakano and Strayer 2014), the continued spread of zebra mussels in Texas represents a critical biodiversity concern. Finally, zebra mussel infestation can challenge human water use, including drinking water treatment and power generation (Connelly et al. 2007), representing a potentially large economic burden.
Managers interested in applying our results for dreissenid invasion risk assessment could consider sites where both climate and in-lake chemistry appear suitable for dreissenid establishment to be the highest-risk sites, while sites that appear suitable at one scale of analysis or another (e.g., climatically suitable but not demonstrating suitable water chemistry) could be considered lower-risk for dreissenid mussel invasion. Overall, our model predictions of extensive zebra mussel habitat in Texas paired with the potential negative consequences of an expanding invasion warrant further prevention efforts, management, and research.