The Extent of Seasonally Suitable Habitats May Limit Forage Fish Production in a Temperate Estuary

The sustained production of sufficient forage is critical to advancing ecosystem-based management, yet factors that affect local abundances and habitat conditions necessary to support aggregate forage production remain largely unexplored. We quantified suitable habitat in the Chesapeake Bay and its tidal tributaries for four key forage fishes: juvenile spotted hake Urophycis regia, juvenile spot Leiostomus xanthurus, juvenile weakfish Cynoscion regalis, and bay anchovy Anchoa mitchilli. We used information from monthly fisheries surveys from 2000 to 2016 coupled with hindcasts from a spatially interpolated model of dissolved oxygen and a 3-D hydrodynamic model of the Chesapeake Bay to identify influential covariates and construct habitat suitability models for each species. Suitable habitat conditions resulted from a complex interplay between water quality and geophysical properties of the environment and varied among species. Habitat suitability indices ranging between 0 (poor) and 1 (superior) were used to estimate seasonal and annual extents of suitable habitats. Seasonal variations in suitable habitat extents in Chesapeake Bay, which were more pronounced than annual variations during 2000–2016, reflected the phenology of estuarine use by these species. Areas near shorelines served as suitable habitats in spring for juvenile spot and in summer for juvenile weakfish, indicating the importance of these shallow areas for production. Tributaries were more suitable for bay anchovy in spring than during other seasons. The relative baywide abundances of juvenile spot and bay anchovy were significantly related to the extent of suitable habitats in summer and winter, respectively, indicating that Chesapeake Bay habitats may be limiting for these species. In contrast, the relative baywide abundances of juvenile weakfish and juvenile spotted hake varied independently of the spatial extent of suitable habitats. In an ecosystem-based approach, areas that persistently provide suitable conditions for forage species such as shoreline and tributary habitats may be targeted for protection or restoration, thereby promoting sufficient production of forage for predators. Further, quantitative habitat targets or spatial thresholds may be developed for habitat-limited species using estimates of the minimum habitat area required to produce a desired abundance or biomass; such targets or thresholds may serve as spatial reference points for management.

The sustained production of sufficient forage is critical to advancing ecosystem-based management, yet factors that affect local abundances and habitat conditions necessary to support aggregate forage production remain largely unexplored. We quantified suitable habitat in the Chesapeake Bay and its tidal tributaries for four key forage fishes: juvenile spotted hake Urophycis regia, juvenile spot Leiostomus xanthurus, juvenile weakfish Cynoscion regalis, and bay anchovy Anchoa mitchilli. We used information from monthly fisheries surveys from 2000 to 2016 coupled with hindcasts from a spatially interpolated model of dissolved oxygen and a 3-D hydrodynamic model of the Chesapeake Bay to identify influential covariates and construct habitat suitability models for each species. Suitable habitat conditions resulted from a complex interplay between water quality and geophysical properties of the environment and varied among species. Habitat suitability indices ranging between 0 (poor) and 1 (superior) were used to estimate seasonal and annual extents of suitable habitats. Seasonal variations in suitable habitat extents in Chesapeake Bay, which were more pronounced than annual variations during 2000-2016, reflected the phenology of estuarine use by these species. Areas near shorelines served as suitable habitats in spring for juvenile spot and in summer for juvenile weakfish, indicating the importance of these shallow areas for production. Tributaries were more suitable for bay anchovy in spring than during other seasons. The relative baywide abundances of juvenile spot and bay anchovy were significantly related to the extent of suitable habitats in summer and winter, respectively, indicating that Chesapeake Bay habitats may be limiting for these species. In contrast, the relative baywide abundances of juvenile weakfish and juvenile spotted hake varied independently of the spatial extent of suitable habitats. In an ecosystem-based approach, areas that persistently provide suitable conditions for forage species such as shoreline and tributary habitats may be targeted for protection or restoration, thereby promoting sufficient production of forage for predators. Further, quantitative habitat targets or spatial thresholds may be developed for habitat-limited species using estimates of the minimum habitat area required to produce a desired abundance or biomass; such targets or thresholds may serve as spatial reference points for management.

INTRODUCTION
Trophic interactions among aquatic predators and prey are rarely incorporated in stock assessments (Skern-Mauritzen et al., 2016), yet quantification of such interactions is critical to advancing ecosystem-based management. In the mid-Atlantic, predators such as summer flounder Paralichthys dentatus, striped bass Morone saxatilis, and bluefish Pomatomus saltatrix use estuarine and offshore habitats to support their critical life functions, and their seasonal migrations and trophic interactions maintain connectivity between coastal and offshore ecosystems (Scharf et al., 2004;Latour et al., 2008;Overton et al., 2008). Although feeding habits of many predators are well studied, the distribution and abundance of prey species that comprise the forage base have received less attention (but see Arbeider et al., 2019;Woodland et al., 2021). Estuaries that support relatively high forage production offer spatially extensive habitat conditions that sustain recruitment, growth, and survival; conversely, forage production may be low in estuaries where habitat conditions are degraded. The relationship between population abundance and the extent of suitable habitats has been reported for many species (Holbrook et al., 2000;Parsons et al., 2014;Sundblad et al., 2014;Weber et al., 2017), but has not been widely explored for forage species. Such a relationship, however, may reveal conditions under which habitats limit forage production.
In this study, we focused on forage fishes in Chesapeake Bay because of the availability of temporally and spatially rich data for these taxa and because ecosystem-based approaches to management are currently pursued in the region (Atlantic States Marine Fisheries Commission [ASMFC], 2018;Freitag et al., 2018;Leslie, 2018). The health and sustainability of iconic fisheries in this system depend on sufficient production and availability of forage as well as effective management and protection from anthropogenic degradation of habitats. With the exception of Woodland et al. (2021), habitat conditions necessary to support forage production in this system are not well understood. Our objectives were to (1) quantify suitable habitats for several forage fishes in Chesapeake Bay, and (2) assess the relationship between the extent of suitable habitats and annual abundance of these species. We considered four numerically dominant forage fishes (Tuckey and Fabrizio, 2020): juvenile (age 0) spotted hake Urophycis regia, juvenile spot Leiostomus xanthurus, juvenile weakfish Cynoscion regalis and bay anchovy Anchoa mitchilli. Small-bodied fishes such as these are important components of the diets of resident and transient predators in Chesapeake Bay Latour, 2011, 2015). Because we selected taxonomically and ecologically disparate species, we expected that suitable habitats would be defined by habitat features that differed among species. If the extent of suitable habitats limits the production of forage fishes in Chesapeake Bay, then we would expect annual patterns in forage fish abundances to exhibit patterns similar to those for suitable habitats.
Static features such as substrate type are often used to characterize fish habitats because such features affect distribution and habitat use (Day et al., 1989;Fabrizio et al., 2013). Dynamic environmental conditions such as salinity, temperature, dissolved oxygen (DO), and depth also contribute to variations in the distribution and abundance of estuarine and coastal species. In river-dominated estuaries, river flow affects salinity and alters the extent of suitable habitats for juvenile fishes (Kostecki et al., 2010), many of which may serve as forage for predators. Temperature is a key determinant of habitat suitability for ectotherms because temperature governs critical processes such as metabolic rates, movement, and growth (Little et al., 2020). Low DO conditions are believed to limit suitable habitats for fishes, especially during summer when some estuarine and coastal systems exhibit prolonged seasonal hypoxia. In particular, abundance of fish is low in hypoxic (<2 mg O 2 /l) waters (Craig and Crowder, 2005;Tyler and Targett, 2007;Zhang et al., 2009;Buchheister et al., 2013;Glaspie et al., 2019), suggesting that individuals actively avoid hypoxic habitats. Other habitat features such as bottom-current velocities, water column stability, and salinity stratification contribute to hydrodynamic complexity of estuarine systems and as such, may shape variations in the spatial distribution and abundance of estuarine organisms (Manderson et al., 2011;Jenkins et al., 2015;Bever et al., 2016). Indeed, hydrodynamic models can be used to estimate habitat volume for estuarine species when coupled with information on physiological tolerances and bioenergetics requirements (e.g., Schlenger et al., 2013). Hydrodynamic models have also been used to assess the effect of sea-level rise on fishes that depend on marsh habitats for juvenile growth and survival (Fulford et al., 2014), to assess the extent of suitable habitats for fishes in coastal environments (e.g., Le Pape et al., 2003;Bever et al., 2016;MacWilliams et al., 2016), and to evaluate potential impacts of climate change on extent of habitats (Crear et al., 2020a,b).
Fish-habitat relationships are best derived from observations across broad spatial scales and long time periods (Gray et al., 2011;Lecours et al., 2015), thus, we quantified these relationships for each of the four species in Chesapeake Bay and its subestuaries during the 17-year period, 2000-2016. We developed an integrated modeling framework to couple information on the abundance of forage fishes with environmental conditions estimated from two models of Chesapeake Bay, rather than considering only those habitat features measured at the time of fish sampling. The primary data were monthly catches from fishery-independent surveys of forage fishes, hindcasts of dynamic environmental conditions (covariates describing salinity, temperature, current speed, depth, and DO conditions), and estimates of static habitat conditions (sediment composition and distance to shore). We applied a data-driven approach, boosted regression tree analysis (Elith et al., 2008), to select a subset of habitat covariates that were most influential in explaining fish relative abundance. Non-parametric suitability models using the histogram approach (Tanaka and Chen, 2015;Guan et al., 2016) were then constructed using the selected influential environmental covariates. Higher suitability is ascribed to conditions in which greater abundances of organisms are observed, and as such, habitat suitability models are processbased models. We used spatial distributions of the environmental covariates to examine seasonal habitat suitability throughout Chesapeake Bay and tributaries during 2000-2016, because most of the species we studied are seasonal migrants; all of the species studied use the Chesapeake Bay as a nursery, but the nursery function is temporally restricted. Finally, we assessed the role of habitat area in driving forage fish abundance across the 17 years of study.

Estimation of Relative Abundance
Geo-referenced catches of forage fishes were obtained from two bottom-trawl surveys: the Virginia Institute of Marine Science Juvenile Fish Trawl Survey (hereafter, Virginia survey) and the Maryland Department of Natural Resources Blue Crab Summer Trawl Survey (hereafter, Maryland survey). The sampling domain of the Virginia survey includes waters greater than 1.2 m throughout Virginia tidal waters of the Chesapeake Bay and its major tributaries (Figure 1). Each month, from January to December, the Virginia survey sampled fishes from 111 stations selected from a random stratified survey design (Supplementary Table 1). A 30 semi-balloon bottom trawl was deployed for 5 min at each site; protocol details are available in Tuckey and Fabrizio (2016). The Maryland survey is primarily a shallow-water survey (mean depth = 2.1 m; 99.7% of sites < 5 m deep) that samples fishes from fixed sites in tributaries and sounds of the Maryland portion of the Chesapeake Bay (Figure 1). A 16 semi-balloon otter trawl is towed for 6 min at each site. In 2000 and 2001, sampling was conducted monthly from May through October at 37 sites in the Chester River, Choptank River, Eastern Bay, Patuxent River, Pocomoke Sound, and Tangier Sound. In 2002 and thereafter, 16 additional sites were sampled in Fishing Bay, the Little Choptank River, and the Nanticoke River (57 sites total; Figure 1 and Supplementary For each trawl tow, we expressed the relative index of abundance of juvenile (age-0) spotted hake, juvenile spot, juvenile weakfish, and all life stages of bay anchovy as the catch per unit effort (CPUE), where effort was estimated by the area swept by the net. Area swept was calculated by multiplying the geodetic distance of the tow by the effective net width, estimated as 55% of the headline spread based on gear tests performed in a flume tank. To relate forage fish abundance to suitable habitat areas and to ensure that CPUE best represented relative abundance, we restricted consideration of CPUE measures to those seasons in which individuals were available to the gear: spring (March, April, May) for juvenile spotted hake; summer (June, July, August) and fall (September, October, November) for juvenile spot and juvenile weakfish; and spring, summer, fall, and winter (December, January, February) for bay anchovy. Note that no sampling was completed in Maryland in winter, thus the bay anchovy CPUE index in winter was based on catches from Virginia waters only.
A Bayesian hierarchical method (Conn, 2010) was used to estimate baywide relative abundance for each species using the survey CPUEs. This method uses the coefficient of variation to weight the individual data sources and extracts a single annual index to represent the pattern exhibited by the multiple indices under the assumption that component indices are subject to process error (from variation in catchability, spatial distribution, etc.) and sampling error (i.e., within-survey variance; Conn, 2010). Annual baywide indices of relative abundance and their associated 95% credible intervals were estimated for the season(s) of interest for each forage species. We used WinBugs accessed through an R script (R Core Team, 2019) to perform these calculations.

Static Habitat Covariates
Two static variables, distance to shore (km) and percent fine sediment, were determined for locations at the midpoint of each trawl tow (Table 1). Fringing marsh and other shallow water areas may provide resources that enhance survival and growth of forage fish (e.g., refuge from predators and provisioning of food; Manderson et al., 2004;França et al., 2009;Boutin and Targett, 2019), and as such, distance to shore may influence fish habitat use. The shortest distance to the nearest shoreline was calculated, and in some cases, the nearest shoreline was an in-Bay island (Figure 2A). Seabed percent fine sediment, a key feature of fish habitat (Kritzer et al., 2016), was determined from a baywide surface grain-size distribution map ( Figure 2B) based on observed surface seabed grain size (Moncure and Nichols, 1968;Byrne et al., 1983;Kerhin et al., 1988;Velinsky et al., 1994;Maryland Geological Survey, 1996;Reid et al., 2005).

Bottom-Water Dissolved Oxygen
Dissolved oxygen concentrations (mg O 2 /l) were hindcast for bottom waters of the Chesapeake Bay and its tributaries using the methods of Du and Shen (2014)

Dynamic Habitat Covariates
Estimates of temperature, salinity, depth, and current speed were obtained from a three-dimensional model of the Chesapeake Bay developed using the UnTRIM hydrodynamic model Zanolli, 2002, 2005). This model takes advantage of the grid flexibility allowed in an unstructured mesh by gradually varying grid cell sizes, beginning with large grid cells in the Atlantic Ocean and transitioning to finer grid resolution in the smaller channels of the tributaries and the northern portion of the estuary (Figure 3). This approach allows for the model to accurately capture the bathymetry and shoreline at multiple spatial scales while maintaining suitable simulation times using a single high-end workstation computer. Further description of the inputs for the hydrodynamic model are provided in the Supplementary Material, along with a summary of the model validation to the data most relevant to this study. This validation demonstrated that the model accurately estimated temperature and salinity in the Bay and major tributaries under a wide range of environmental conditions (Supplementary Figures 1, 2); model accuracy was similar to that reported for a suite of Chesapeake Bay models evaluated by Irby et al. (2016). A large number of environmental variables were initially considered for use in developing fish habitat suitability models to eliminate a priori specification of environmental conditions that may be important for describing abundance and distribution of forage fishes. Twenty-two dynamic environmental variables were calculated, along with sediment composition and distance to shore ( Table 1). Environmental covariates were extracted from the hydrodynamic model and the DO model at the time and location of the individual tows (midpoint of the tow) to allow us to couple fisheries observations with hindcasts of environmental covariates. We considered hindcasts of environmental covariates at multiple temporal and spatial scales because such measures may provide more accurate predictions of habitat suitability (Lecours et al., 2015). Dynamic variables were extracted from the hydrodynamic model as instantaneous values at the time of each tow (i.e., the time of capture), and as the tidal-averaged (24.8 h) values in the tidal cycle preceding the capture event. Conditions preceding capture may influence the distribution of fishes (e.g., Jaonalison et al., 2020), and preliminary investigations with boosted regression trees suggested that models that considered tidal-averaged conditions performed better than those with comparable instantaneous values; henceforth, we considered tidal-averaged conditions. Depth-averaged conditions were obtained for salinity, temperature, and current speed ( Table 1). Tidal-averaged current speeds provide a metric of flow which may be used by some species to aid movements within the estuary (Brady and Targett, 2013). We also considered covariates describing near-bed conditions (1 m above the seabed), and maximum depth-averaged current speed. Because vertical or horizontal gradients in current speed may act to aggregate food near complex currents or fronts, we also considered these covariates. The vertical gradient in the current speed was calculated as the difference between the current speed one meter above the bottom and one meter below the surface; the horizontal gradient in the current speed was calculated as the maximum difference in current speed between adjacent model grid cells divided by the distance between the model grid cells.
Although habitat conditions near the seabed may be most relevant for understanding fish-habitat relationships for demersal species such as spotted hake, habitat use may also reflect overall conditions in the water-column because even demersal fishes are not confined to near-bed habitats. For example, salinity stratification may influence the supply of food or DO to the near-bed region sampled by the trawl, and may be an indicator of forage fish occurrence. Thus, we used covariates describing salinity and temperature stratification, that were calculated as the difference between instantaneous surface (top 1 m) and near-bed (1 m above seabed) conditions. In addition, because favorable habitats may be characterized by a range of conditions, we considered covariates based on the percent of time that nearbed conditions fell within a given range, for example, the percent of time that salinity exceeded 20 PSU. The percent of time within a given salinity or temperature range was calculated over the same time interval used for tidal averaging. We identified three salinity ranges (<10 PSU; 10-20 PSU; >20 PSU), and three temperature ranges (<10 • C; 10-20 • C; >20 • C) consistent with observed patterns in fish communities in Chesapeake Bay (Tuckey and Fabrizio, personal observation).

Selection of Influential Habitat Covariates
Boosted regression trees (BRTs) were used to select a subset of influential habitat covariates that explained variations in the nominal catch rates of fish (Breiman et al., 1984). The regression tree algorithm uses recursive partitioning to explain variation in the response (nominal catch rates), that is, observations are repeatedly split into increasingly homogeneous groups based on threshold values of the predictors (habitat covariates; Breiman et al., 1984). Cross-validation was used to assess model fit and to ensure that the resultant trees were applicable to out-of-sample observations; cross-validation was achieved by fitting the tree to a subset of the data (the training set) and fit was assessed using the remaining data (the test set). Furthermore, the performance of regression tree algorithms may be improved with ensemble methods such as boosting, which aggregates multiple trees to enhance the stability of the resultant model (Knudby et al., 2010). We used a Poisson response to model the number of fish captured per tow with the gbm.step procedure (R package 'dismo, ' Elith et al., 2008;Elith and Leathwick, 2017;R Core Team, 2019). Catches from the Virginia survey were used without modification (numbers per 5-min tow), but catches from the Maryland survey were expressed in 5-min-tow equivalencies rounded to the nearest integer. All habitat covariates were standardized to permit direct comparison of covariate importance (Schielzeth, 2010).

Optimization of Boosted Regression Tree Models
Prior to fitting the BRTs to the 2000-2016 observations, we optimized the model-fitting parameters of the BRTs by exploring the combination of parameter values that produced the lowest deviance of the cross-validated data sets (Elith et al., 2008;Cameron et al., 2014). To determine optimal parameters values for the BRTs and because optimization is computationally intensive, we used a subset of observations (2010-2012; n = 4,604 tows) that represented notably different environmental conditions (2011 was a wet year compared with 2010) as well as large differences in the observed relative abundance of forage species. BRT parameters were optimized separately for spotted hake, spot, weakfish, and bay anchovy using the gbm.step procedure in R (dismo package, R Core Team, 2019). Model fitting failed to produce at least 1,000 trees for bay anchovy using all of the 2010-2012 observations, so we optimized BRT parameters using data from each season separately. For bay anchovy in fall, BRT modeling yielded less than 1,000 trees and was thus unreliable (Elith et al., 2008) and not considered further.
Optimization focused on selection of the learning rate and tree complexity, two model-fitting parameters assigned by the analyst. The learning rate determines how quickly the model approximates the observed data (Miller et al., 2016), and the tree complexity represents the level of interaction possible among the predictors. The third parameter selected by the analyst is the bag fraction, or the proportion of observations used to train the model. Observations in the training subset are selected randomly without replacement for each model run and the remaining observations are used for cross-validation. Preliminary investigations suggested that a bag fraction of 0.75 was reasonable. Using this bag fraction, we fitted a series of trees to a range of learning rates (lr = 0.0005, 0.0050, 0.0075, 0.0100, 0.0200, 0.0300, 0.0400, 0.0500, 0.0750) and tree complexities (tc = 1, 3, 5, 10), similar to Cameron et al. (2014). For each species, we considered only those BRTs for which at least 1,000 trees were fit and selected parameter values that reduced the cross-validated deviance (Elith et al., 2008). We identified the optimal tree complexity for each species by graphically examining the change in cross-validated deviance across learning rates. Next, using the selected tree complexity, we identified the learning rate that produced the minimum cross-validated deviance.
To select influential habitat covariates, we fitted BRT models for each species for the period 2000-2016 (N total = 25,333 tows; N Virginia = 20,326 tows, N Maryland = 5,007 tows) using a bag fraction of 0.75 and values of the optimized species-specific learning rates and tree complexities determined by optimization. Optimized model-fitting parameters (lr and tc) varied among species: juvenile spotted hake lr = 0.01 and tc = 10; juvenile spot lr = 0.02 and tc = 10; juvenile weakfish lr = 0.005 and tc = 10; and bay anchovy (spring, summer, winter) lr = 0.02 and tc = 3. In addition, optimization runs indicated that the six covariates describing percent time were least informative, so these were not considered further. Therefore, a suite of 18 covariates (16 dynamic, 2 static; Table 1) were considered for the BRT models. Estimates of variable influence and scree plots produced by the gbm.step procedure (R Core Team, 2019) allowed us to identify and select a subset of influential covariates for each species, and for bay anchovy for each season (spring, summer, winter). Care was taken to consider only those covariates that did not exhibit high correlations with other influential covariates, that is, only those covariates with r 2 < 0.8 were considered in the subsequent calculation of habitat suitability indices (HSIs). This approach avoids overweighting of the HSI for a particular habitat condition.

Habitat Suitability Models
Habitat suitability models were used to assign habitat suitability scores and to quantify the extent of suitable habitat for forage fishes throughout the Chesapeake Bay and its tributaries from 2000 to 2016 for each season. Habitat suitability models were estimated with the non-parametric histogram approach because this approach makes no assumption about the nature of the relationship between environmental features and fish abundance (Guan et al., 2016). Briefly, thresholds of environmental conditions that resulted in a gradient of suitability from least suitable (0) to most suitable (1) were identified for each influential habitat covariate for each species. The HSI, calculated as the mean of two or more covariate-specific suitability indices (SIs), also ranged between 0 and 1 to ease interpretation.

Suitability Indices
We estimated SIs for the range of observed values for each of the influential habitat covariates identified by the species-specific BRTs. We used the Tanaka and Chen (2015) approach to estimate SIs but applied a disjoint clustering method to identify 'natural clusters' of the habitat covariates for the histogram approach; we implemented this method with the FastClus procedure in SAS/STAT R . Tanaka and Chen (2015) fixed the number of individual bins to 10 for each habitat covariate, but we found that this resulted in bins with few observations (<5) or narrowly defined environmental limits (e.g., bottom temperature between 16.2 and 16.3 • C). Thus, we allowed the number of bins to vary (but not exceed 10), and restricted cluster sizes to a minimum of 40 observations; in all cases, the smallest cluster included at least 54 observations, allowing a reasonable description of average abundance (and hence, relative suitability) in each cluster. SIs were estimated for each cluster and habitat covariate using where SI ij is the suitability index for cluster j of habitat covariate i, CPUE ij is the average catch (fish/km 2 ) observed in cluster j of habitat covariate i, and CPUE i,min and CPUE i,max are the minimum and maximum average catches observed across all clusters of habitat covariate i (Tian et al., 2009;Chang et al., 2012). This formulation allows the SIs to range between 0 and 1.0, with 1.0 indicating the most suitable condition and 0 the least. More explicitly, each covariate cluster, defined by a range of values, was associated with an SI score.

Habitat Suitability Indices
Habitat suitability indices were calculated by expressing the HSI as an average of multiple covariate-specific SIs; we restricted the number of covariates in the HSI to those that were most influential as determined by BRTs ( Table 2). The HSI for a given set of covariates may be expressed as an arithmetic or geometric mean of the individual SIs (e.g., Brown et al., 2000;Tanaka and Chen, 2015). A single averaging approach to estimate the HSIs, however, may not be appropriate for all species (e.g., Yu et al., 2019), so we explored both models of the mean. The arithmetic mean model for the HSI is where SI 1 is the suitability index for habitat covariate 1, SI 2 is the suitability index for covariate 2, and so forth; and p is the number of covariates considered (e.g., Hess and Bay, 2000). The geometric mean model for the HSI is HSI gm = p SI 1 × SI 2 × SI 3 × · · · SI p (e.g., Layher and Maughan, 1985;Lauver et al., 2002;Tian et al., 2009). The geometric mean index applies the concept of a 'limiting factor' whereby a low SI for a single covariate results in a low HSI gm (Zajac et al., 2015). HSI calculations were performed in SAS R or Matlab (MathWorks Inc.).
Habitat suitability index models were calibrated by using fish catches at each tow location and graphically examining the relationship between the HSIs and the average relative abundance for each species-season combination (e.g., Tanaka and Chen, 2015). We used the 5% trimmed mean as a measure of the average because trimmed means are insensitive to the occasional extreme catches observed for some species. For a properly calibrated HSI, the mean relative abundance of forage fish is expected to increase as habitat conditions approach optimal for the species, that is, as the HSI increases from 0 to 1.0.

Verification of Modeling Approach
We verified the use of BRTs for selection of covariates and evaluated the reliability of the two averaging formulations of HSI for forage fishes using bootstrap resampling (Efron and Tibshirani, 1986). About 70% of the fisheries observations (N = 18,121) comprised the training data set, and the remaining ∼30% (N = 7,212) was used as the test (or verification) data set. The SurveySelect procedure in SAS/STAT R was implemented to randomly select samples without replacement using a stratified design to ensure representation across years, seasons, and geographic areas (Maryland and Virginia). Training and test data sets were constructed separately for each species, and bay anchovy sets were constructed separately for spring, summer, and winter. For each training data set, we fitted BRTs using the same bag fraction and species-specific lr and tc as before; we then selected influential covariates, and modeled the HSI am and HSI gm . Note that the BRTs for each data set may have indicated a different number of influential covariates, as well as a different suite of influential covariates, than what was identified by the original BRT model fitted to observations from 25,333 tows. The resulting HSI models were applied to each of the test data sets to estimate the predicted HSIs. Due to computational intensity, 10 cross-validation data sets were generated (consistent with Pennino et al., 2020). The expected performance of the HSI am and HSI gm for each species across all seasons was evaluated with the root mean square error (RMSE), calculated as the standard deviation of the residuals (where residuals represent the difference between the predicted HSI and observed HSI for each location where fish were sampled). For bay anchovy, we estimated RMSEs for spring, summer, and winter individually. We used a paired t-test implemented in the glm procedure in SAS R to assess differences in mean RMSEs, and retained the HSI formulation that exhibited the lower mean RMSE for further analyses.

No change Yes
Model fitting failed for bay anchovy when data were considered across seasons, so we fit boosted regression trees (BRTs) for each season separately; for bay anchovy in fall, BRT modeling failed to produce at least 1,000 trees and was thus not considered further. DO is dissolved oxygen. Extent of suitable habitat was calculated for each season as the sum of the areas throughout Chesapeake Bay with habitat suitability indices ≥ 0.5; SE is the standard error of the mean. The change in extent of suitable habitats evaluates the monotonic change (increase, decrease, no change) from 2000 to 2016; the relationship between abundance and extent of suitable habitat is based on results from a nonparametric regression (see Table 3).

Estimation of the Extent of Suitable Habitat
The extent of suitable habitat for each species was quantified (objective 1) by calculating HSIs from the environmental covariates at each hydrodynamic model grid cell for each season from 2000 through 2016. To facilitate GIS visualization of seasonal habitat conditions and calculation of the area of suitable habitat for each species, we used the median of the daily values of the covariates to represent the seasonal average for a given model grid cell, season, and year. Median values for each habitat covariate at each model grid cell were then used to estimate the HSIs for each species. In this manner, we mapped the species-specific seasonal HSIs at the spatial resolution of the hydrodynamic model because processes operating at small The model fitted to the data was Y i = β 0 + β 1 X 1i + ε i where Y i is the ranktransformed response, i = 1 to N, N is the sample size, β 0 is the overall average response (intercept), β 1 is the regression coefficient (slope), X 1i is the value of the predictor for observation i, and ε i is the unexplained random error. Extent of suitable habitat was calculated for each season as the sum of the areas throughout Chesapeake Bay with HSI values ≥ 0.5. spatial scales may be masked when environmental conditions are averaged over large spatial scales (Windle et al., 2012). The areas of the individual grid cells where HSI exceeded a given threshold of suitability (i.e., HSI ≥ 0.5, HSI ≥ 0.6, HSI ≥ 0.7, HSI ≥ 0.8) were summed to obtain an estimate of the extent of 'suitable' habitat throughout the Chesapeake Bay and its tributaries for each species, season, and year. Preliminary graphical investigations revealed that annual patterns in suitable habitat extents were similar among the 0.5, 0.6, and 0.7 thresholds; the 0.8 threshold yielded values that were too low to be useful. The 0.5 threshold, used to describe habitat suitability for the eastern oyster Crassostrea virginica in Chesapeake Bay (Theuerkauf and Lipcius, 2016) and for pelagic sharks in Australia (Birkmanis et al., 2020), was used for subsequent analyses and presentation.

Relationship Between Suitable Habitat Extent and Relative Abundance of Forage Species
We hypothesized that annual and seasonal changes in the area of suitable habitat affect the abundance of forage species in this temperate ecosystem. To address this hypothesis, we related the annual time series of suitable habitat with annual estimates of baywide relative abundance for each species (objective 2). We limited the exploration of these relationships to the season during which a particular species was most vulnerable to the trawl gear: juvenile spotted hake in spring, juvenile spot in summer and fall, juvenile weakfish in summer and fall, and bay anchovy in spring, summer, and winter. We rank transformed abundance indices because of the small number of observations (n = 17 years except for bay anchovy in winter when n = 16) and fit non-parametric regressions to model the relationship between rank abundance and extent of suitable habitat. For this analysis, we assumed the relationship was stationary, that is, the effect of suitable habitat extent on the abundance of forage fishes remained stable from 2000 to 2016 (e.g., Zeng et al., 2018). In addition, we tested the null hypothesis that the seasonal extent of suitable habitat was constant in Chesapeake Bay between 2000 and 2016; as before, we rank transformed the extent of suitable habitat (defined as areas with HSI ≥ 0.5). Computations for non-parametric regression analyses were performed with the rank and glm procedures in SAS R . Residual plots indicated a reasonable fit of the rank regression model to the data.

Influential Habitat Covariates
Environmental conditions and habitat features that comprised suitable habitats varied among species and ranged between four and seven, depending on species ( Table 2). Water depth and one of the current speed metrics were consistently identified as influential covariates for all species; one of the temperature covariates was influential in describing suitable habitats for forage fishes in spring, summer, and fall, but was not selected for describing suitable habitats in winter ( Table 2). At least one salinity metric defined suitable habitats for juvenile spotted hake and bay anchovy, and distance to shore explained suitable habitats for juvenile spot, juvenile weakfish, and bay anchovy ( Table 2). Bottom DO conditions delineated suitable habitats for juvenile spot in winter and bay anchovy in winter and spring ( Table 2). Conditions at sites sampled in Maryland waters differed from those in Virginia waters: in Maryland, sampled habitats tended to be shallower, closer to shore, warmer in summer, and cooler in fall than habitats sampled in Virginia (Supplementary Figure 3). Most notably, Maryland sites exhibited lower bottom DO concentrations in summer than Virginia sites (Supplementary Figure 3). In addition, relative to sites in Virginia, Maryland sites exhibited lower surface salinities, less stratification in terms of salinity and temperature, lower current speeds and less stratification in current speeds (Supplementary Figure 3).

Verification and Calibration of the Habitat Suitability Index Modeling Approach
Bootstrap analyses verified that BRTs were useful for selection of influential covariates; in general, the same or similar covariates were identified as most influential among the 10 bootstrap realizations. For juvenile spotted hake, the HSI formulation based on the geometric mean (HSI gm ) provided the best approximation to the original HSI estimated for each sample as indicated by the significantly lower RMSE (t = 4.56, P < 0.01). Unlike results for juvenile spotted hake, the HSI based on the arithmetic mean (HSI am ) performed better for bay anchovy, regardless of season (t spring = -3.08, P spring < 0.01; t summer = -4.01, P summer < 0.01; t winter = -3.27, P winter < 0.01). Although we found no evidence for a difference in the mean RMSEs of the HSI am and the HSI gm for juvenile spot and juvenile weakfish (t spot = -1.13, P = 0.27; t weakfish = -1.65, P = 0.12), we used the HSI am for these species because the mean RMSE of the HSI am was consistently less than the mean RMSE of the HSI gm .
FIGURE 4 | Relationship between HSI and trimmed mean catches of juvenile spotted hake, juvenile spot, juvenile weakfish, and bay anchovy in Chesapeake Bay, 2000-2016. For juvenile spotted hake, the HSI gm is shown, whereas the HSI am is shown for other species. HSI values were binned for ease of plotting such that bin 0.0 includes all observations where 0 ≤ HSI < 0.2, bin 0.2 includes all observations where 0.2 ≤ HSI < 0.4, and so forth. Note differences in y-axes.
Relative abundance as measured by the trimmed mean of the four species increased with increasing values of HSIs (Figure 4), indicating proper calibration of habitat suitability models. The ranges of observed HSI values across years were 0-0.95 for juvenile spotted hake in spring; 0.12-0.98 for juvenile spot in summer; 0.04-0.86 for juvenile spot in fall; 0.09-0.99 for juvenile weakfish in summer; 0-0.89 for juvenile weakfish in fall; 0-0.92 for bay anchovy in spring; 0-0.92 for bay anchovy in summer; and 0.04-0.98 for bay anchovy in winter.

Juvenile Spotted Hake
We detected a strong seasonal pattern in the extent of suitable habitats for juvenile spotted hake in Chesapeake Bay: little to no suitable habitat was available in summer and fall, increased in winter (mean winter = 467.0 km 2 or 4.3% of the total area modeled), and was greatest in spring (mean spring = 2,046.0 km 2 or 18.8% of the total area; Table 2, Figure 5A, and Supplementary Figure 4A). In spring, the extent of suitable habitat was greater in the Virginia portion of the Chesapeake Bay system than in the Maryland portion (Supplementary Figure 5A). Suitable habitats for spotted hake in spring were relatively deep, away from shorelines, and where salinity stratification was greater than 4.9 psu ( Figure 5A); these habitats were characterized by tidal-averaged bottom temperatures ranging between 5.3 and 14.2 • C. Annual changes in the extent of suitable habitats varied by as much as 90.3% in winter ( Figure 6A) and as much as 38.8% in spring (Figure 6B) across the time period of study. The extent of suitable winter habitat was higher during years when waters in the region began to warm earlier in the year (e.g., 2012) and was lower when warming was delayed (e.g., 2011; Figure 6A).

Juvenile Spot
The extent of suitable habitat for juvenile spot displayed a persistent seasonal pattern, with relatively greater extents of  suitable habitats in spring (mean spring = 3,169.4 km 2 or 29.1% of the total area) and summer (mean summer = 4,650.7 km 2 or 42.6% of the total area) and relatively little suitable habitat in fall and winter (mean fall = 775.4 km 2 or 7.1% of the total area; mean winter = 864.0 km 2 or 7.9% of the total area; Table 2, Figure 5B, and Supplementary Figure 4B). Suitable habitats for juvenile spot were primarily found in shallow areas near shorelines in spring, and in the deeper portions of the Chesapeake Bay and its tributaries in summer ( Figure 5B). In tributaries and embayments such as Mobjack Bay, the extent of suitable habitat for juvenile spot in spring was greater than that in summer by as much as 47.2%. Suitable habitat extents in spring and summer were greater in the Maryland portion of the Chesapeake Bay system compared with the Virginia portion (Supplementary Figure 5B). The extent of suitable habitats for spot was driven by distance to shore, depth, temperature stratification, bottom DO, and horizontal gradients in current speed, and was not well described by a single covariate (Table 2). Nevertheless, suitable habitats exhibited bottom DO levels less than 4.0 mg O 2 /l in summer and less than 5.3 mg O 2 /l in fall, suggesting that juvenile spot used habitats that may be considered marginal or unsuitable for other species. Low bottom DO conditions are associated with stratified water column conditions and indeed, juvenile spot were more likely to occur in habitats where tidal-averaged temperature stratification exceeded 2.2 • C in summer and 2.7 • C in fall.

Juvenile Weakfish
A notable seasonal pattern was evident in the amount of suitable habitat available for juvenile weakfish: little suitable habitat was available in winter (mean winter = 667.3 km 2 or 6.1% of the total area), increased in spring (mean spring = 2,446.3 or 22.4% of the total area) and was greatest in summer (mean summer = 4,394.8 or 40.3% of the total area; Table 2, Figure 5C, and Supplementary Figure 4C). Throughout the Chesapeake Bay system, the extent of suitable habitat in fall was, on average, 19.1% greater than in spring, but this pattern was not observed in 2000, 2010, and 2012 (Supplementary Figure 4C). In summer, suitable habitats for juvenile weakfish were found close to shorelines; in fall, suitable habitats were away from shore, near the mouth of the Potomac River and the lower Chesapeake Bay (Figure 5C). In Virginia waters, we observed similar extents of suitable habitats in summer (mean summer = 2,056.7 km 2 or 18.9% of the total area) and fall (mean fall = 2,007.7 km 2 or 18.4% of the total area), but in Maryland waters, suitable habitat extents declined by an average of 61.2% from summer (mean summer = 2,338.0 km 2 or 21.4% of the total area) to fall (mean fall = 906.3 km 2 or 8.3% of the total area; Supplementary Figure 5C). Regional differences in fall were clearly depicted by seasonal maps: HSI values in the mainstem of the bay were greater in waters south of the Rappahannock River than north of the Rappahannock River ( Figure 5C). Suitable habitats for juvenile weakfish were characterized by distance to shore, depth, bottom temperature, and current speed stratification ( Table 2).

Bay Anchovy
The estimated extent of suitable habitat for bay anchovy varied seasonally, with the greatest extent of suitable habitat area in the Chesapeake Bay system in spring (mean spring = 4,976.9 km 2 or 45.6% of the total area), followed by summer (mean summer = 3,809.0 km 2 or 34.9% of the total area) and winter (mean winter = 2,906.3 km 2 or 26.6% of the total area; Table 2, Figure 5D, and Supplementary Figure 4D). Although this seasonal pattern was also observed for Virginia waters (mean spring = 3,117.2 km 2 or 28.6%; mean summer = 2,774.4 km 2 or 25.4%; mean winter = 1,839.8 km 2 or 17.0%), suitable habitat extents in Maryland waters were greatest in spring (mean spring = 1,859.7 km 2 or 17.0%) and lower but similar in summer and winter (mean summer = 1,034.6 km 2 or 9.5%; mean winter = 1,066.5 km 2 or 9.8%; Supplementary Figure 5D). The greatest extents of suitable habitats in tributaries occurred in spring; in summer and winter, the upper reaches of major tributaries were particularly unsuitable for bay anchovy ( Figure 5D). Suitable habitat conditions occurred in areas that reflected complex relationships between multiple environmental and geophysical covariates and consistently included sediment composition, depth, and distance to shore.

Relative Baywide Abundance and Relationship to Extent of Suitable Habitat
Seasonal indices of relative abundance for forage fishes in Chesapeake Bay varied across years, and interannual patterns were generally similar for the Virginia and Maryland surveys, particularly for juvenile spotted hake in spring, juvenile spot in summer and fall, and juvenile weakfish in summer (Figure 7). Inconsistent patterns in interannual changes in relative abundance between Virginia and Maryland waters were observed for juvenile weakfish in fall and bay anchovy in spring and summer (Figure 7), suggesting that seasonal processes affecting the distribution of these species varied across regions of the Chesapeake Bay. We note that the mean index of relative abundance for juvenile weakfish in fall was several orders of magnitude lower in Maryland waters than in Virginia. Similarly, the mean indices of relative abundance for bay anchovy in spring and summer were an order of magnitude lower in Maryland waters than in Virginia. Spatial differences in relative abundance reflected the geographic differences in the extents of suitable habitats which were notably lower in Maryland than in Virginia for these species-season combinations (Supplementary Figures 5C,D).
Two contrasting relationships were detected between the ranked baywide relative abundance index and the extent of suitable habitat for forage fishes. We observed a significant positive relationship between seasonal ranked baywide relative abundance and extent of suitable habitat for juvenile spot in summer (Tables 2, 3 and Figure 8A). The baywide relative abundance index for juvenile spot in summer was highly variable with contrasting periods of low and high abundance ( Figure 9A). Extents of suitable summer habitat for juvenile spot exhibited no significant linear pattern across time (Tables 2, 3 and Figure 9A). We observed a similar relationship for bay anchovy in winter: the extent of suitable habitat in winter was a significant determinant of the ranked relative abundance of bay anchovy ( Table 3 and Figure 8B). The baywide relative abundance of bay anchovy varied widely in winter and exhibited no obvious temporal pattern ( Figure 9B); extents of suitable winter habitat for bay anchovy showed no systematic change through time (Tables 2, 3 and Figure 9B).
More commonly, we were unable to detect a significant relationship between the area of suitable habitat and the ranktransformed estimate of baywide relative abundance of forage fishes (Tables 2, 3). For juvenile spotted hake in spring, we found no indication that the extent of suitable habitat was limiting, except perhaps in 2002 when the area of suitable habitat declined below 1,600 km 2 and the ranked abundance index was among the lowest observed in the time series. This, however, may be coincidental. The extent of suitable spring habitat for spotted hake varied without trend since 2000 (Tables 2, 3). We found no evidence of an effect of the extent of suitable habitat on the ranked relative abundance of juvenile spot in fall (Tables 2, 3). The extent of suitable habitat for juvenile spot in fall exhibited no trend through time (Tables 2, 3) and was markedly less than in summer. Between 2000 and 2016, the extent of suitable habitat for juvenile weakfish increased significantly in summer and fall (Tables 2, 3 and Figures 10A,B). The relative abundance of juvenile weakfish, however, exhibited no detectable response to increases in the extent of suitable habitats in either summer or fall during the study period (Tables 2, 3). Relative abundances of bay anchovy in spring and summer were highly variable (Figure 7) and annual estimates were imprecise. Although the extent of suitable spring and summer habitats for bay anchovy increased significantly since 2000 (Tables 2, 3 and Figures 10C,D), we did not detect a response in the ranked relative abundance of bay anchovy to changes in the extent of suitable habitats in spring or summer (Tables 2, 3).

DISCUSSION
Our modeling framework combined the power of machine learning to identify influential habitat covariates with the flexibility of non-parametric approaches to characterize habitat suitability and the capabilities of GIS to quantify and depict suitable (and unsuitable) habitats for forage fishes in Chesapeake Bay from 2000 to 2016. We coupled catch information from fishery surveys with static features of the environment and outputs from models of dynamic conditions to depict suitable habitats in Chesapeake Bay. In an ecosystem-based approach, these habitats may be targeted for protection (e.g., by limiting habitat alterations or other human activities that may incidentally kill or injure forage fishes) or restoration (e.g., by improving water quality conditions), thereby promoting production of sufficient forage for predators. Importantly, our modeling approach for building forage-fish habitat suitability models for the Chesapeake Bay was calibrated and verified, thereby allowing estimation of habitat suitability for tributaries and embayments in Chesapeake Bay that are not routinely sampled by fishery surveys (e.g., Mobjack Bay, Potomac River). Furthermore, our results allow resource managers to focus protective measures on critical habitat areas in Chesapeake Bay, for example, shallow areas near coastlines in spring and summer which persistently support suitable habitats for multiple forage species. Our finding that shallow areas near coastlines are important habitats for some forage species is consistent with the observed ontogenetic habitat shift of juvenile weakfish, which move from salt marsh tributaries to shallow habitats near coastlines in summer (Boutin and Targett, 2019). We found annual patterns in suitable habitat extent that mirrored those of baywide relative abundance for two of the species examined, juvenile spot in summer and bay anchovy in winter; as such, estimates of the minimum habitat area required to produce a desired abundance (or biomass) of FIGURE 7 | Seasonal standardized abundance indices for forage fishes in Maryland (blue) and Virginia (red) waters of Chesapeake Bay, 2000-2016: juvenile spotted hake in spring, juvenile spot in summer and fall, juvenile weakfish in summer and fall, and bay anchovy in spring, summer, and winter. Seasons were March-May for spring, June-August for summer, September-November for fall, and December-February for winter. Seasonal relative abundance indices (mean catch per unit effort expressed as number of fish/km 2 ) were standardized to a mean of 1.0 across the 17 years, thus, patterns of abundance can be compared between states, but these standardized abundance indices do not reflect differences in estimated mean catch rates within a given year. Note that bay anchovy were not sampled in Maryland in winter, and thus, only the standardized index for Virginia is depicted.
forage fish can be used to establish quantitative habitat targets (Kritzer et al., 2016) or spatial thresholds that may serve as spatial reference points for management (Reuchlin-Hugenholtz et al., 2016). Quantitative habitat targets and spatial reference points for bay anchovy and juvenile spot in the Chesapeake Bay warrant further consideration.
Suitable seasonal habitat extents for forage species exhibited annual changes reflecting annual-scale heterogeneity in habitat conditions in Chesapeake Bay, with persistent seasonal signals that varied among species. Current speed, water depth, and either temperature or DO were important covariates for the four species we examined, and distance to shore was important for three species; thus, suitable habitat conditions resulted from a complex interplay between water quality and the geophysical properties of the habitat. Variations in seasonal extents were more pronounced than annual variations in suitable habitat extent, supporting the notion that the Chesapeake Bay serves as a nursery area for juvenile fishes and that the nursery function is temporally restricted and varies among species (e.g., spring for juvenile spotted hake, and summer for juvenile spot and weakfish). For some species, the extent of suitable seasonal habitat increased since 2000 (summer and fall habitats for juvenile weakfish, and spring and summer habitats for bay anchovy), whereas for other species, extents varied annually with no clear trend. None of the species examined were at the southern limit of their geographic range, and as waters of the Chesapeake Bay continue to warm (Hinson et al., 2021), we expect that suitable habitat extent may increase for species with broad thermal tolerances such as spot, weakfish, and bay anchovy. Other climate-related effects predicted for the region include increased precipitation and sea-level rise, both of which may alter salinity and salinity stratification. Salinity mediates the thermal tolerances of some fishes (e.g., Lankford and Targett, 1994;Nepal and Fabrizio, 2020) and could serve to limit suitable habitats in the future. When coupled with field observations, laboratorybased investigations of the interactive effects of salinity on the thermal tolerances for these forage species could be informative (e.g., Lankford and Targett, 1994).
The relationship between the extent of suitable habitat and relative abundance of forage species was species-dependent and when present, varied seasonally. Such relationships have not been widely explored for aquatic species. Although manipulative field experiments permit exploration of these relationships on small spatial scales (e.g., Parsons et al., 2014), only a few  Table 3). Values of HSI am ≥ 0.5 were considered suitable habitat. studies attempted to relate extent of suitable habitat and relative abundance of aquatic species using large-scale, field-based observations (e.g., Holbrook et al., 2000;Le Pape et al., 2003;Sundblad et al., 2014;French et al., 2018;Yu et al., 2019). For example, Yu et al. (2019) used a graphical assessment to note the consistency between declines in the relative abundance of neon flying squid Ommastrephes bartramii and the spatial shrinkage of suitable habitats in the northwest Pacific Ocean. Sundblad et al. (2014) used linear regression to examine relationships between relative adult abundance and the percent of the study area that was suitable nursery habitat for 12 populations of pikeperch Sander lucioperca and European perch Perca fluviatilis. Here, we used non-parametric regression to statistically evaluate the strength of such relationships for forage fishes in Chesapeake Bay during 17 years and found a positive relationship between suitable habitat extent and baywide relative abundance of juvenile spot in summer and bay anchovy in winter, indicating that environmental and geophysical conditions affect the carrying capacity of the Chesapeake Bay for these species during a portion of the year. As such, our results are consistent with the findings of a meta-analysis of the correlation between abundance and indicators of habitat suitability (Weber et al., 2017).
We found no evidence that suitable habitat extents were limiting in Chesapeake Bay for juvenile spotted hake in spring, juvenile spot in fall, juvenile weakfish in summer and fall, and bay anchovy in spring and summer. These results suggest that seasonal suitable habitat extent exceeded that necessary to support these populations, and that other factors such as predation mortality (e.g., Minello et al., 1989), food availability (e.g., Tableau et al., 2016), or degradation of egg and larval habitats (e.g., Sundblad et al., 2014) may have contributed to changes in relative abundances. For example, suitable habitat extent for juvenile weakfish in summer and fall increased significantly since 2000, but was not significantly related to changes in relative abundance of juvenile weakfish, suggesting that factors other than those considered in the suitability model exerted a greater role in driving annual fluctuations in abundance. Indeed, abundance and growth of juvenile weakfish are greatest in habitats that contain abundant mysid resources and that are found adjacent to large expanses of salt marsh (Boutin and Targett, 2019). Consideration of the spatial distribution and annual variation in mysid abundance may be required to further delineate suitable habitats for juvenile weakfish in Chesapeake Bay. Additionally, the lack of a significant relationship between annual relative abundance of juvenile weakfish and suitable habitat extent may have resulted from changes in predation mortality of juvenile weakfish. This hypothesis is consistent with the observed increase in natural mortality rates for this species in the 2000s (Atlantic States Marine Fisheries Commission [ASMFC], 2019; Krause et al., 2020); although the sources of increased natural mortality in weakfish remain unclear, increased predation is believed to have played a role (Atlantic States Marine Fisheries Commission [ASMFC], 2019). Abundances of other forage fishes may be limited by similar trophic interactions with their predators or prey, or by unsuitable egg and larval habitats in Chesapeake Bay.
Although habitat extents limited the relative abundance of juvenile spot in summer, no such limitation occurred in fall. We hypothesize that the summer-to-fall decoupling of the relationship between habitat extent and relative abundance of juvenile spot may be partly explained by temperature. Mean water temperatures in Chesapeake Bay are greatest in late August-early September and the increased metabolic rates and energy demands of predators during this time may increase predation mortality on juvenile spot. Rising water temperatures during this time may also affect the phenology of spot emigration resulting in earlier emigration and fewer juvenile spot remaining in the estuary during fall in warmer years. We were unable to detect a change in the extent of suitable fall habitats for juvenile spot between 2000 and 2016, but warming water temperatures in fall associated with directional climate change may alter future habitat conditions in fall for this species. Continued monitoring of fall abundances and a better understanding of the cues that trigger spot emigration in fall will be necessary to address this hypothesis.
Another possible reason for our inability to detect a relationship between suitable habitat area and baywide relative abundance of some forage fishes is that the suite of static and dynamic habitat covariates we considered failed to adequately describe suitable habitat conditions. For example, the abundance of predators, availability of prey, and proximity to biogenic habitats such as seagrass or oyster reefs may help to shape the distribution and abundance of forage fishes. In addition, seasonal relative abundance indices for many of the forage species were imprecise and hence were statistically invariable across time (based on 95% credible intervals of the baywide hierarchical index). Thus, a 'good' year with a relatively high mean index of abundance was not statistically discernible from a 'poor' year with a relatively low index; this lack of contrast may have hampered our ability to uncover a relationship between relative abundance indices and suitable habitat extents. These results suggest that if the abundance of forage fishes in Chesapeake Bay is changing, the temporal and spatial intensity of sampling by current fisheries surveys is insufficient to detect changes. Alternatively, and more likely, abundance may be stable but highly variable across years, as is typically exhibited by forage fishes (Hilborn et al., 2017). Continued monitoring will be critical to detect directional changes in abundance.
Seasonal variation in the geographic location of suitable estuarine habitats may be linked to variation in freshwater input (e.g., Rubec et al., 2019). In Chesapeake Bay, freshwater input influences salinity and salinity stratification; however, we identified additional hydrodynamic covariates such as temperature stratification and current speed that contributed to variation in suitable habitats. For example, for bay anchovy, the suitability of winter habitats was partly determined by bottom DO and the horizontal gradient in tidal-averaged current speeds. Water temperature was not considered in the HSI model for bay anchovy in winter, however, we cannot rule out temperature as an important covariate describing habitat conditions for bay anchovy in winter because bottom temperature and DO were correlated. Interestingly, the HSI model for juvenile spot also included bottom DO and a measure of temperature stratification, suggesting that multiple hydrodynamic factors are required to describe suitable habitats.
In this study, we considered information from two fisheryindependent surveys, one of which yielded fewer annual observations but sampled shallow-water habitats across a large portion of the estuarine salinity gradient. Overall, our fisheries observations were collected at a relatively fine spatial resolution (>100 sites sampled/month) and high temporal intensity (monthly), thereby minimizing biases due to seasonal or shortterm habitat use (e.g., by sampling only one or two months each year). By using fisheries observations from multiple surveys, we were able to detect the effect of bottom DO on the relative baywide abundance of juvenile spot in summer. Low bottom DO conditions in summer are more prevalent and persistent in Maryland waters of the Chesapeake Bay; in particular, the Chester River, Eastern Bay, Choptank River, Little Choptank River, and Patuxent River exhibited bottom DO levels in summer that were markedly lower than what was observed in Virginia waters. These observations extended the range of bottom DO conditions over which these effects could be investigated. In addition, the Maryland survey sampled sites in shallow habitats close to shore, and these conditions were not well sampled in Virginia waters. Together, the two surveys provided observations from a greater range of environmental and bathymetric conditions commonly encountered in the Chesapeake Bay region.
Hydrodynamic models and other numerical models of environmental conditions can provide information on dynamic habitat features that are not measured at the time of sampling and represent a significant advancement toward refining spatial relationships between fish and their environment (e.g., Crear et al., 2020a). Consideration of such information may yield habitat models with greater predictive accuracy (Scales et al., 2017). In our study, we used tidal-averaged conditions to develop habitat suitability models that reasonably reflected the relationship between (daily) environmental conditions and relative abundance (at the tow level) of forage fishes in Chesapeake Bay. This fine-scale approach is preferable to one that uses seasonal averages of habitat conditions to build habitat suitability models (Scales et al., 2017). Indeed, the temporal resolution of environmental covariates used to build the suitability model affects the scale of inference. For example, Woodland et al. (2021) described annual patterns in the distribution and abundance of forage fishes and invertebrates relative to patterns of predation and environmental conditions in the Chesapeake Bay and its major tributaries. Seasonal changes, however, could not be addressed because several habitat conditions were represented by annual means (e.g., average discharge from tributaries, and the Atlantic Multidecadal Oscillation [AMO] index). Large-scale climatic changes as indexed by the AMO affect mean abundances of bay anchovy and juvenile spot (Woodland et al., 2021), and our findings for bay anchovy in winter and juvenile spot in summer are consistent with findings in Woodland et al. (2021). Specifically, our results suggested that variations in environmental conditions contributed to the observed variation in relative abundance of these forage species. Although we did not consider large-scale climate indices per se, we did examine smallscale environmental indicators of climate change (temperature, salinity) and demonstrated how these changes affect habitat suitability and relative baywide indices of abundance for juvenile spot in summer and bay anchovy in winter. Unlike Woodland et al. (2021) who found greater relative abundance of bay anchovy in the upper bay than in the lower bay, our indices of relative abundance for bay anchovy in summer were an order of magnitude lower in Maryland than those observed in Virginia waters, likely because we lacked samples from deep sites (>2.7 m) in Maryland which were considered in Woodland et al. (2021). Also, the time scale of our studies differed. Similar to Woodland et al. (2021), we found that bay anchovy and juvenile spot exhibited higher relative abundances in southern tributaries than in northern tributaries of the Chesapeake Bay, but our spatial depiction of suitable habitat conditions throughout the system allowed us to examine the fine-scale spatial distribution of suitable habitats. Such depictions are helpful for identifying the geographic scale at which habitat protection or restoration must be implemented.
As the number of habitat descriptors available to researchers increases, variable reduction techniques are critical to the selection of covariates that help explain the variation in observed abundance and distribution of aquatic organisms. For example, satellite imagery, ocean observing systems, and hydrodynamic models yield a multitude of environmental descriptors of habitat and these data are commonly used to inform fish habitat models. Similar to Georgian et al. (2019), we used BRTs to identify influential covariates from a large suite of possible covariates. Because tree complexity can play a large role in improving the outcome of cross-validation and BRT model fitting, BRT parameters should be optimized using the data under consideration. Many researchers either fail to optimize regression trees or when optimization is implemented, only a single parameter is optimized (typically learning rate; e.g., Georgian et al., 2019;Yu et al., 2020) after using the default bag fraction (0.75) and selecting an arbitrary value for tree complexity (typically between 2 and 5; e.g., Georgian et al., 2019;Pennino et al., 2020;Yu et al., 2020). We recommend optimization of BRTs using the approach we applied here or the R package gbm.auto (Dedman et al., 2017).
Across the four species we examined, the geometric mean formulation of the HSI was best for only a single speciesjuvenile spotted hake. Although the HSI gm is widely used, this formulation may penalize the index too harshly for mobile species that can tolerate broad variations in environmental conditions, including sub-optimal conditions for limited periods of time. For instance, in areas where bottom temperature exceeds a given species' thermal tolerance, the SI for bottom temperature may be close to 0; in this case, the value of the HSI gm will also be close to 0, but other environmental conditions in these areas may be suitable, even optimal, and thus, the overall suitability may not be well indexed by the HSI gm . Because the appropriate HSI formulation depends on species, we recommend use of data-driven analyses and assessment of model performance to inform selection (Chang et al., 2012;Tanaka and Chen, 2015;Yu et al., 2019; this study).
Our fisheries observations from two surveys that sampled across a broad geographic area of the largest estuary in the United States represented 17 years of monthly sampling and reflected the breadth of habitat conditions that fishes were likely to encounter in Chesapeake Bay. Although survey designs differed (stratified random design in Virginia, and fixed-site design in Maryland), the integration of information from such surveys can provide models with good predictive performance as long as observations are spatially extensive (Soranno et al., 2020). Furthermore, when the number of observations used to fit the model is sufficiently large, then projections for unsampled areas within the same time frame are considered valid interpolations (Elith and Leathwick, 2009;Soranno et al., 2020). Our projections of HSIs in areas not sampled by the Maryland or Virginia surveys (e.g., Mobjack Bay, Potomac River) were based on 25,333 observations, and as such, were valid interpolations for assessment of habitat conditions in non-sampled areas during the timeframe of the study (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). Our interpolations did not 'extend beyond the conditions represented by the data used to fit the model, ' and thus we avoided extrapolations to areas where novel combinations of predictors occur (Conn et al., 2015).
Finally, we note that the uncertainties associated with habitat suitability modeling and resulting projections are not typically assessed (Elith and Leathwick, 2009), although such uncertainties are useful for understanding the limitations of model-based results for conservation and fisheries management. Uncertainties may arise from model specification (e.g., the type of model used, or covariates omitted from the model) and from the observations used to fit the model (e.g., samples may not represent the population of interest, or sample size may be inadequate). Ensemble approaches have been used to partially mitigate model uncertainty, but ensemble models do not fully overcome the limitations of the individual component models (Elith et al., 2010). Model and observational uncertainties may affect habitat suitability projections in different ways, and uncertainty analysis for HSI models warrants further research (Zajac et al., 2015) and engagement with environmental statisticians.

DATA AVAILABILITY STATEMENT
The data analyzed in this study are subject to the following licenses/restrictions: William & Mary ScholarWorks (https: //scholarworks.wm.edu/) provides open public access to the monthly abiotic data collected as part of the VIMS trawl survey (https://doi.org/10.25773/97cg-7w42), the seasonal habitat suitability maps (https://doi.org/10.25773/d37g-1h89), and the estimated areas of suitable habitat (https://doi.org/10.25773/ d37g-1h89) for 2000 to 2016. Maps and data are under a CC-NC-SA 4.0 license. Maryland trawl survey data should be requested from the original data source.

AUTHOR CONTRIBUTIONS
All the authors conceptualized the research, performed the analyses, interpreted the results, read and approved the published version of the manuscript. MF prepared the original draft of the manuscript. TT, AB, and MM contributed to writing and editing of the manuscript. MF and AB contributed to the project administration and visualization. MF and TT contributed to the funding acquisition.

FUNDING
This study was funded by the NOAA Chesapeake Bay Office, grant NA17NMF4570156 to the authors.

ACKNOWLEDGMENTS
The UnTRIM code used in the hydrodynamic modeling was developed and provided by Professor Vincenzo Casulli (University of Trento, Italy). We thank J. Shen (Virginia Institute of Marine Science, VIMS) and B. Marcek (US Fish & Wildlife Service) for extending the model of bottom DO conditions. We are grateful to J. Buchanan (VIMS) and V. Nepal (VIMS) for application and verification of the DO model, to W. Lowery (VIMS) for preparing Figure 1, and to R. Dixon (VIMS) and S. Smith (VIMS) for assistance with verification of the HSI modeling approach. We thank C. Walstrum (MD DNR) for providing fisheries data from the Maryland trawl survey. We appreciate the steadfast sampling in Virginia tidal waters by Captains H. Brooks, V. Hogge, and W. Lowery, and the scientific crew of the VIMS Juvenile Fish Trawl Survey that participated in trawl sampling between 2000 and 2016. Finally, we thank the reviewers for comments that helped to improve the manuscript. This manuscript is Contribution No. 4051 of the Virginia Institute of Marine Science, William & Mary.