Study Areas

The Meramec, Gasconade, and Little Black river basins are within the Ozark Highlands ecoregion of Missouri (Fig. 1) and share similar physiographic and watershed features. In the interior of the region, dolomite and sandstone comprise the dominant bedrock, while the western outer regions are dominated by Mississippian limestone (Ozark Ecoregional Assessment Team 2003). All three watersheds have steep bluffs along streams, narrow valleys, and karst features, and many of their streams are spring-fed. Seasonal patterns of discharge are similar among all three streams (Fig. 2), but discharge in the Little Black River is much lower than the Gasconade and Meramec rivers because of its smaller watershed (990 km² and 7,268 km² for the Little Black and Gasconade rivers, respectively). We describe additional features of the Gasconade and Little Black rivers below; a description of the Meramec River basin can be found in Key et al. (2021).

Gasconade River.—The mainstem Gasconade River flows north for 436 river-km (rkm) before joining the Missouri River (Blanc 2001). Our habitat suitability models included about 800 rkm including the mainstem Gasconade River and three of its tributaries, Osage Fork, Big Piney River, and Roubidoux Creek (Fig. 1). These streams are not channelized or impounded, but in-channel gravel mining has altered and destabilized some segments (Blanc 2001), and decreased riparian vegetation has also contributed to channel instability and erosion (Jacobson and Primm 1997). Forty-six mussel species are reported from the Gasconade River basin (Blanc 2001).

Little Black River.—The mainstem Little Black River flows south 137 rkm into Arkansas before joining the Current River. Most of the Little Black River and its tributaries are within the Ozark Highlands, but the downstream portion of the mainstem flows through the Mississippi Alluvial Plains (Fig. 1; Wilkerson 2003). Because of the physiographic differences between the Ozark Highlands and Mississippi Alluvial Plains, we did not include that portion of the stream in our habitat suitability models. Our habitat suitability models included 120 rkm comprising the mainstem Little Black River and three of its tributaries, North and South prongs and Beaverdam Creek (Fig. 1). The Little Black River is highly altered with 13 impoundments and about 98 rkm of channelized streams (Wilkerson 2003). Thirty-nine mussel species are reported from the Little Black River basin (Wilkerson 2003).

Mussel Survey Dataset

We determined mussel bed locations in the Gasconade and Little Black rivers from the Missouri Department of Conservation mussel database (data available upon request to and subject to the approval of the Missouri Department of Conservation, 3500 East Gans Road, Columbia, MO 65201). This database includes mussel survey information for specific locations across Missouri, including GPS points, survey methods, lists of species found, and numbers of individuals found. We used mussel survey data from 1994 to 2013, following the MRHSM (Key et al. 2021). We filtered the data to include only timed-search samples; incidental collections, collections using a groping technique, or entries with missing sampling method were excluded. We considered sites within 180 m of each other to represent the same mussel bed (Lueckenhoff 2015; Schrum 2017; Key et al. 2021). Our resulting dataset included 130 unique mussel bed locations. We selected a subset of 85 mussel beds that had the highest species richness (hereafter species-rich mussel beds; Key et al. 2021) to develop our habitat suitability models. The remaining 45 mussel beds (hereafter validation mussel beds) were used to validate our models by determining how many of these beds fell within habitat deemed suitable by our model.

Figure 1.

Map of Missouri, USA, showing the Gasconade (pink), Little Black (dark pink), and Meramec (blue) river watersheds and major streams included in habitat suitability models. Other text and boundaries within the state boundary are Level IV ecoregions.

Generation of Hydrogeomorphic Variables

We derived 12 hydrogeomorphic variables for the Gasconade and Little Black rivers, including all 10 variables used in the MRHSM and two additional variables that we created (see below; Table 1). These variables represent habitat characteristics thought to correspond to suitable habitat for mussels at a riverscape scale, including bluff adjacency, presence of and proximity to gravel bars, lateral channel stability, low-flow surface water availability, and stream power index (Table 1; Key et al. 2021). Mussel beds in Ozark rivers are often found in the vicinity of bluffs, possibly because bluffs exert channel control and stabilization that is amenable to mussel establishment and persistence (Vannote and Minshall 1982; Key et al. 2021). Mussel beds often are associated with gravel bars, and the presence of persistent gravel bars after high-flow events can indicate channel stability (Bates 1962; Peck 2005; Zigler et al. 2008; Key et al. 2021). Lateral channel movement is indicative of bank erosion and sediment deposition, which can destabilize substrate and limit mussel occurrence (Strayer 1999; Strayer et al. 2004). Low-flow surface water availability is intended to represent a proxy for the existence of permanently watered areas that serve as refugia during drought periods (Table 1; Golladay et al. 2004; Key et al. 2021). Stream power is an index of potential energy in the channel and influences channel erosion and stability.

We derived estimates for the hydrogeomorphic variables from high-resolution, open-source datasets of aerial imagery and topography following Key et al. (2021) and summarized as follows. We began our workflow by defining the stream dimensions and location of the river channel and subsequently creating a stream centerline. We then generated points on the stream centerline at 10-m cross sections to create a spatially continuous dataset. After we defined our stream dimensions and stream centerline, we derived our 12 hydrogeomorphic variables and assigned the data to each point on the stream centerline. The 10-m points were then interpolated using natural neighbors into continuous grids representing our final hydrogeomorphic variables. All spatial analyses were performed in ArcGIS and projected to NAD 1983 UTM zone 15N (ESRI 2011; Key et al. 2021).

Figure 2.

(A) Mean monthly unit discharge (ft³/s/mile²) and (B) mean monthly discharge (ft³/s) for the Meramec, Gasconade, and Little Black rivers. Flow data are from  https://waterdata.usgs.gov/mo/nwis/rt (accessed January 23, 2023).

We used a combination of 1-m light detection and ranging (LiDAR) and 10-m digital elevation model (DEM) coverage of both study areas (MSDIS 2011) to generate six hydrogeomorphic variables related to bluff adjacency and stream power. Two bluff adjacency variables represented whether a bluff is present within one channel width of each bank from a mussel bed (binary variable) and, if so, the total bluff area (continuous variable). In addition to the bluff adjacency variables from the MRHSM, we generated two new variables representing the total bluff area within 500 m upstream or downstream of a mussel bed. We added these variables to explore whether bluffs located upstream or downstream, versus directly adjacent to the stream channel, are associated with mussel bed occurrence. Because of the limited availability of LiDAR for the Gasconade River, we used 10-m DEMs to extend the remotely sensed data across the entire drainage area. The 1-m horizontal resolution LiDAR tiles and 10-m DEMs were mosaicked into a single DEM and resampled to 10-m resolution for analysis (Key et al. 2021). We then used a slope and range criteria to define bluffs in the watersheds and performed a zonal search at each point on the stream centerline.

Table 1.

Justification and description of hydrogeomorphic variables evaluated in habitat suitability models for the Meramec, Gasconade, and Little Black rivers. Variables that were included in the final, best-fit models differed among streams and transferability levels (see text). “Type” refers to whether the variable was continuous or binary.

For the stream power variables, we estimated the watershed area at each stream centerline point and estimated the stream power index (SPI) using the bankfull elevation (see Key et al. 2021 for methods to determine bankfull elevation) as SPI = S × ln(A_d), where, for any location along the stream centerline, SPI is the stream power index, S is the channel slope, and A_d is the watershed area (Moore et al. 1991). Slope was averaged over a 500-m interval, 250 m upstream and 250 m downstream from each point on the stream centerline, and then smoothed using a 50-m moving average. The binary stream power variable was derived by classifying each pixel as either high or low stream power using the mean value as the break between the two classes.

We used National Agriculture Imagery Program leaf-off aerial imagery to generate the six remaining hydrogeomorphic variables related to gravel bars, lateral channel stability, and low-flow surface water availability. We derived three variables reflecting the presence of or proximity to gravel bars for each mussel bed (Table 1). We classified mussel bed locations as “gravel” if the stream reach was dominated by persistent gravel before and after a high-flow year, and as “pool” if the reach was dominated by water (binary variable). We then determined whether each mussel bed was located within 100 m of a gravel bar (binary variable), and we created a continuous variable representing the Euclidean distance from a mussel bed to the nearest gravel bar. We derived these variables with a differencing technique between two sets of aerial imagery. In the Little Black River, 2007 and 2015 were low-water years, while 2013 was a high-water year. In the Gasconade River, 2012 and 2014 were low-water years, and 2013 was a high-water year. Pixels that changed state (water or gravel) between the two images were classified as a pool, and pixels that did not change state were classified as either gravel or a pool. The gravel/pool class therefore does not represent the underlying sediment (gravel versus depositional sediments) but rather areas that had persistently exposed gravel bar versus areas that were predominantly water during low-flow conditions. Without ground-truthing, we cannot differentiate whether the areas classified as pools had gravel or depositional sediments.

For the lateral channel stability variable, we created two polygons representing the stream banks based on visual cues such as shadow, vegetation, and scour lines in the leaf-off imagery from 1990–95 and 2015. We defined each point on the stream centerline as unstable if the channel moved > 10 m between the two time periods or stable if the channel moved ≤ 10 m (binary variable; Key et al. 2021).

We derived two variables for low-flow surface water availability index using imagery taken at the time of lowest discharge available (2007 for the Little Black River and 2012 for the Gasconade River). We performed a focal search to estimate the number of pixels classified as water that were adjacent to the focal pixel (continuous variable). We then used the median value to categorize high- and low-surface water availability in a low-flow period to create the binary variable. We acknowledge that these variables do not directly represent vulnerability to drying because water depth is not accounted for (see Key et al. 2021). However, we used these variables as proxies for drought refugia because bathymetric data were not available from our imagery.

Habitat Suitability Models

We used maximum entropy modeling (MaxEnt; Phillips and Dudik 2008) to generate habitat suitability models for mussels in the Meramec, Gasconade, and Little Black rivers. This method uses presence-only occurrence data in combination with environmental data layers to produce a model of habitat suitability spanning a specified geographic area (Phillips and Dudík 2008). For models that included the Meramec River, we used the presence-only occurrence data and environmental layers created by Key et al. (2021). We spatially constricted our habitat suitability models to each drainage and used the location of species-rich mussel beds in each river system in combination with the hydrogeomorphic variables determined for those locations. We used the same settings in MaxEnt as used for the MRHSM (Key et al. 2021). Specifically, we set the run type to bootstrap to generate training and test occurrence data (80% and 20% of the species-rich mussel bed locations, respectively), and we ran models with 10,000 background points and 5,000 iterations. None of the hydrogeomorphic variables included in each model were correlated with each other (correlation coefficient < 0.40).

We converted the raw model results to a binary map of suitable and unsuitable reaches based on the equal test sensitivity and specificity logistic threshold of each model (Key et al. 2021). The equal test sensitivity and specificity logistic threshold is a commonly used threshold that sets the sensitivity equal to the specificity (Cao et al. 2013; Phillips 2017; Key et al. 2021). After suitable and unsuitable reaches were delineated, we used a buffer of 40 m to separate the suitable and unsuitable habitats to account for areas of transition (following Key et al. 2021). We then used jackknife analysis and the test gain values to assess the relative contribution of each hydrogeomorphic variable and to determine which variables were most important for model fit (Phillips 2017; Key et al. 2021). We used a stepwise model selection approach of our hydrogeomorphic variables to select the best-fit model. The area under the receiver operating curve (AUC) values from MaxEnt provided relative values for comparing the performance of models that were built with the same data (Phillips et al. 2006). Therefore, we selected variables that led to higher AUC values and contained hydrogeomorphic variables with sizeable individual effects on model results when others were removed (following Elith 2002 and Key et al. 2021). We also created response curves to investigate the relationships between suitable and unsuitable habitats (y-axis) and our hydrogeomorphic variables (x-axis). The results for our continuous hydrogeomorphic variables were presented as curves spanning the range of values for that layer, whereas the binary hydrogeomorphic variables were presented as two bars representing the binary. The range in values of the continuous—or bars of the binary—hydrogeomorphic variables were classified as suitable if they were equal to or above the equal test sensitivity and specificity logistic threshold on the response curves. For low-flow surface water availability and stream power, higher values in the response curves represented more contiguous surface water availability during low-flow conditions and higher stream power, respectively. Although AUC values provided comparisons of model performance, they did not provide a measure of the accuracy of habitat suitability (Jiménez-Valverde and Lobo 2006). Therefore, we used the location of the validation mussel beds (not used in model development) to assess the accuracy of our best-fit models. We calculated the percentage of the validation mussel beds that fell within a reach predicted to be suitable by the best-fit models across the entire spatial extent. For models that included the Meramec and Gasconade rivers or the Meramec and Little Black rivers, we also calculated validation per drainage as the proportion of validation mussel bed locations that fell within a reach predicted to be suitable for each drainage, separately.

Transferability

We separated our methods of transferability into three categories representing different levels of dependence on the original MRHSM: Level 1, transferring the original model; Level 2, transferring the modeling framework; and Level 3, adapting the modeling framework (Fig. 3).

Level 1: Transferring the original model.—To transfer the original model from the Meramec River to the Gasconade and Little Black rivers, we utilized the species-rich mussel bed locations and hydrogeomorphic variables used in the MRHSM (Key et al. 2021) and species-rich mussel bed locations and hydrogeomorphic variables that we derived for the Gasconade and Little Black rivers. Specifically, we combined the species-rich mussel bed locations and hydrogeomorphic variables to include the spatial extent of both the Meramec and Gasconade rivers or the Meramec and Little Black rivers. We started this level of transferability with all 10 hydrogeomorphic variables from the MRHSM (Fig. 3). We then used the jackknife analysis and stepwise model selection approach as described previously to find the best-fit model. The results from the best-fit model were converted to a binary map and validated using the validation mussel bed locations. We considered model transfer to be successful if the best-fit model had a test AUC ≥ 0.70 and total validation ≥ 0.70 (Fig. 3). If either the test AUC or model validation was < 0.70, we progressed to Level 2 (Fig. 3).

Figure 3.

General framework used to test the transferability of the Meramec River habitat suitability model to the Gasconade and Little Black rivers, including three levels of transferability and the spatial extent, required data, and evaluation criteria for each level.

Table 2.

Results for the Meramec River habitat suitability model (Key et al. 2021) and evaluation of transferability of that model to the Gasconade and Little Black rivers. All models are best-fit models for each river and level of transferability. We considered model transfer successful if the best-fit model had a test AUC ≥ 0.70 and total validation ≥ 0.70. Equal test sensitivity and specificity logistic threshold is the value used to delineate suitable and unsuitable habitats for all variables in each model. Total validation is the proportion of validation mussel bed locations that were identified by the model as suitable habitat across all rivers in the model. Validation/drainage is the proportion of validation mussel-bed locations identified as suitable habitat within each river.

Level 2: Transferring the modeling framework.—To reduce our dependence on the MRHSM, we transferred only the modeling framework by building our MaxEnt models without data from the Meramec River. Specifically, we used only species-rich mussel bed locations and 10 hydrogeomorphic variables in either the Gasconade River or Little Black River (Fig. 3). Following the same methodology as Level 1, we determined the best-fit model and validated those results with the validation mussel bed locations. Again, if the best-fit model had a test AUC ≥ 0.70 and total validation ≥ 0.70, it was considered a successful transfer of the modeling framework. If either the test AUC or model validation was < 0.70, we progressed to Level 3 (Fig. 3).

Level 3: Adapting the model.—If neither of the previous transferability methods produced an adequate model, we adapted the modeling framework used in the MRHSM by including two additional hydrogeomorphic variables, longitudinal bluff adjacency upstream and downstream. Similar to Level 2, models were built only with the species-rich mussel bed locations and hydrogeomorphic variables from the Gasconade or Little Black rivers. We started this level with the 10 original hydrogeomorphic variables plus the two additional bluff adjacency variables that we created. We followed the same stepwise model selection approach as Levels 1 and 2 to find the best-fit model and then created the binary suitability map and validated the results with the locations of the validation mussel beds. If the best-fit model had a test AUC ≥ 0.70 and total validation ≥ 0.70, we considered this a successful adaptation to the modeling framework (Fig. 3). If either the test AUC or total validation was < 0.70, we considered model transfer unsuccessful (Fig. 3).

Gasconade River

Level 1: Transferring the original model.—The best-fit habitat suitability model for Level 1 had a test AUC of 0.69 (Table 2). The best-fit model included six hydrogeomorphic variables: lateral channel stability, distance to gravel bar, gravel/pool class, stream power index, bluff adjacency area, and low-flow surface water availability index. Jackknife analysis indicated that bluff adjacency area, distance to gravel bar, gravel/pool class, and lateral channel stability contributed significantly to the final model (Table 3).

An equal test sensitivity and specificity logistic threshold of 0.41 separated habitats into suitable and unsuitable. Response curves indicated that suitable habitat was represented by reaches with small bluffs, 0–700 m or > 1,500 m from gravel bars, areas with a low-flow surface water availability index > 0 but < 10, and intermediate stream power. While 68% of our validation mussel bed locations were found in areas identified as suitable in both the Meramec and Gasconade rivers, only 24% of the validation mussel beds were found in areas identified as suitable in the Gasconade River alone. Because the test AUC and total validation were < 0.70, we considered model transfer to the Gasconade River unsuccessful at this level and continued to Level 2 (Table 2).

Level 2: Transferring the modeling framework.—The best-fit habitat suitability model for Level 2 had a test AUC of 0.82 (Table 2). The best-fit model included six hydrogeomorphic variables: lateral channel stability, distance to gravel bar, gravel/pool class, stream power index, bluff adjacency area, and low-flow surface water availability index. Jackknife analysis indicated that bluff adjacency area, distance to gravel bars, gravel/pool class, and low-flow surface water availability contributed significantly to the model (Table 3).

An equal test sensitivity and specificity logistic threshold of 0.34 separated habitats into suitable and unsuitable. Response curves indicated that suitable habitat was represented by reaches classified as pools, near small bluffs, 100–300 m or > 1,250 m from persistent gravel bars, with low-flow surface water availability and stream power indices > 0, and laterally stable channels (Fig. 4 and Table 3). Because the test AUC was > 0.70 and total validation was 0.82, we concluded that transfer of the modeling framework to the Gasconade River was successful, and we did not evaluate Level 3 transferability (Table 2).

Table 3.

Results of the jackknife analyses for the final, best-fit models of the original Meramec River habitat suitability model and transfer of that model to the Gasconade River and Little Black rivers. An asterisk (*) indicates hydrogeomorphic variables that contributed significantly to the model based on the jackknife analyses.

Figure 4.

Response curves for hydrogeomorphic variables that contributed significantly to transferring the Meramec River modeling framework to the Gasconade River at Level 2. The dashed line represents the equal sensitivity and specificity logistic threshold used to delineate suitable and unsuitable habitats.

Little Black River

Level 1: Transferring the original model.—The best-fit habitat suitability model for Level 1 had a test AUC of 0.64 (Table 2). The best-fit model included the same six hydrogeomorphic variables as for the Gasconade River: lateral channel stability, distance to gravel bar, gravel/pool class, stream power index, bluff adjacency area, and low-flow surface water availability index. Jackknife analysis indicated that distance to gravel bar, gravel/pool class, low-flow surface water availability index, stream power index, and lateral channel stability contributed significantly to the final model.

An equal test sensitivity and specificity logistic threshold of 0.42 separated habitats into suitable and unsuitable. Response curves indicated that suitable habitat was represented by reaches near small bluffs, 0–500 m or > 2,000 m from gravel bars, with higher low-flow surface water availability, low-intermediate stream power indices, and laterally unstable channels. Sixty-four percent of validation mussel bed locations were found in areas identified as suitable by the model, and validation in the Little Black River was 60% (Table 2). Because the test AUC and total validation were < 0.70, we considered model transfer to the Little Black River unsuccessful at this level and continued to Level 2.

Level 2: Transferring the modeling framework.—The best-fit habitat suitability model for Level 2 had a test AUC of 0.74 (Table 2). The best-fit model included the same six hydrogeomorphic variables as Level 1. Jackknife analysis indicated that distance to gravel bar, gravel/pool class, stream power index, and lateral channel stability contributed significantly to the final model.

An equal test sensitivity and specificity logistic threshold of 0.48 separated habitats into suitable and unsuitable. Response curves indicated that suitable habitat was represented by reaches near small bluffs, with high water availability, lower stream power indices, and areas classified as pools. The test AUC was > 0.70, but because total validation was only 60%, we concluded that transfer of the modeling framework to the Little Black River was unsuccessful and continued to Level 3 (Table 2).

Level 3: Adapting the model.—The best-fit model for Level 3 had a test AUC of 0.72 (Table 2). The best-fit model included downstream bluff adjacency area, distance to gravel bar, gravel/pool class, low-flow surface water availability, and stream power index. Jackknife analysis indicated that downstream bluff adjacency area, gravel/pool class, low-flow surface water availability, and stream power index contributed significantly to the final model (Table 3).

An equal test sensitivity and specificity logistic threshold of 0.44 separated habitats into suitable and unsuitable. Response curves indicated that suitable habitat was represented by reaches classified as pools, with higher surface water availability and lower stream power indices (Fig. 5 and Table 3). Suitable habitat also was represented by reaches with any amount of downstream bluff area and persistent gravel bars at any distance. The best-fit model at this level had a total validation of 0.80 (Table 2). Because the test AUC was > 0.70 and total validation was 0.80, we concluded that transfer of the model to the Little Black River at this level was successful.