Seasonal patterns of habitat suitability and spatiotemporal overlap within an assemblage of estuarine predators and prey

Subtropical estuaries often experience seasonal climatic shifts, which influence environmental conditions such as water temperature and salinity. The distribution of estuarine organisms depends on the availability of suitable environmental conditions, the extent of which may differ across seasons. Furthermore, ecological interactions among organisms, such as predation and competition, often occur when multiple species overlap spatiotemporally and can influence their distributions and abundances. Determining the habitat associations and overlap patterns among predator and prey species across seasons can thus enhance our understanding of ecosystem function. This study used a long-term monitoring data set to examine the habitat suitability and overlap patterns between 4 predator and 4 prey species in a northwestern Gulf of Mexico estuary. Seasonal (fall and spring) habitat suitability models and maps were created, which allowed for de lineation of spatiotemporal overlap between predators and prey. Environmental drivers of habitat suitability differed among species, but temperature and salinity were frequently the most important variables. Seasonal differences in habitat suitability were observed for all species, but these differences were more notable for prey than for predators. Seasonality in habitat suitability was likely driven primarily by seasonal changes in environmental conditions and movements re lated to reproductive activities. Spatiotemporal overlap between predators and prey was observed primarily in whichever season each prey species was more abundant and in regions where highly suitable habitat was most ubiquitous for each predator. These results suggest that overlap be tween these predators and prey is highly seasonal, implying seasonality in the diets of predators and predation pressure on prey.


INTRODUCTION
Predators are integral to shaping the structure and function of ecosystems through direct predation and indirect, behaviorally mediated effects on prey populations, conspecifics, and other predator species (Paine 1980, Wootton 1993, Heithaus et al. 2008. Habitat degradation has been widespread in coastal and estuarine ecosystems, which can negatively influ-ence predator populations and disrupt ecosystem function (Lotze et al. 2006). The development and enactment of ecosystem-based fisheries management (EBFM) has become increasingly common as many fish populations have been degraded globally (Marshall et al. 2019). The central concept of EBFM is to consider various aspects of the surrounding ecosystem when determining best practices for managing fish populations, including anthropogenic impacts, species interactions, habitat requirements, and climate scenarios (Pikitch et al. 2004). One of the most critical considerations when creating these holistic strategies is the delineation of essential fish habitat, defined as the waters required by a particular species to feed, reproduce, and reach maturity. By delineating essential habitats for multiple species in a given area and across multiple seasons, managers can protect regions that efficiently benefit as many species as possible (Cook & Auster 2005, Moore et al. 2016.
Estuarine ecosystems are particularly dynamic, with vast spatiotemporal heterogeneity in environmental conditions due to river inflow, wind, and tides. The diversity of habitats provided by estuaries is reflected by the diversity of species assemblages that inhabit them (Day et al. 1989). Estuarine communities are generally structured by each species' physiological tolerances and habitat associations, which can range from freshwater-dwelling to fully marine (Jassby et al. 1995, Marshall & Elliott 1998. In temperate and subtropical estuaries, changes in temperature, salinity (via river flow), and water depths (via wind and tides) occur seasonally, and these patterns exhibit interannual variability. Therefore, the extent of essential habitat for a given species can change dramatically based on climatic conditions. The effective delineation of essential fish habitat in this region thus requires long-term monitoring data and modeling that considers spatiotemporal factors in addition to physicochemical conditions. While it is important to determine the environmental and habitat conditions that are ideal for an estuarine assemblage, biotic factors such as species interactions are also often important in shaping community structure (Connell 1961). Predation is a particularly important driving force in most biological communities, with interactions between predators and prey influencing species distributions and abundances (Clark et al. 2003, Trainor et al. 2014. The individual effects of abiotic and biotic factors on species distributions can be challenging to disentangle, especially when data regarding the diets and behaviors of predators are limited. Inferring the potential for species interactions by comparing the extent of essential habitat between predators and prey can allow for some understanding of the biotic forces shaping ecosystem function in the absence of other data. This study used a long-term monitoring data set to examine and compare the habitat associations of an assemblage of estuarine predators and prey in the Galveston Bay Complex (GBC), Texas. Four predators were chosen that exhibit a variety of life-history strategies and exist along a gradient of freshwater habitat use. In approximate order from most to least freshwater-associated, the predators of interest are alli gator gar Atractosteus spatula, bull shark Carcharhinus leucas, black drum Pogonias cromis, and spotted seatrout Cynoscion nebulosus. Adult spotted seatrout exhibit active movement away from intense freshwater inflow events (Callihan et al. 2015) and reductions in metabolic scope and swimming velocity at salinities below 10 ppt (Wohlschlag & Wakeman 1978). Adult black drum are found ubiquitously throughout estuarine salinity gradients and can tolerate extreme salinity conditions . Juvenile bull sharks are thought to associate strongly with freshwater regions as nursery grounds during early life but have been observed moving freely throughout estuaries and exhibiting preferences for a wide range of salinities (Froeschke et al. 2010. Alligator gar are considered freshwater-dwelling fish, residing primarily in rivers and streams. They are known to exist in estuaries, most often near sources of freshwater inflow, but the duration and extent of estuarine use of this species is poorly understood (Buckmeier et al. 2013).
These model estuarine predator species also differ in their prey preferences and feeding strategies. Spotted seatrout exhibit an ontogenetic shift in diet, with juveniles and sub-adults feeding on crustaceans, polychaetes, and small fishes, and adults becoming more piscivorous with age (Darnell 1958, Wenner & Archambault 1996. Adult black drum feed primarily on oysters and the mollusks, bivalves, crustaceans, and polychaetes that are associated with oyster reefs (Overstreet & Heard 1982, Rubio et al. 2018. Juvenile bull sharks are thought to be primarily piscivorous, mainly consuming teleost fishes in the families Mugilidae, Clupeidae, Sciaenidae, and Ariidae (TinHan 2020), but will also consume penaeid shrimps and crabs (Snelson et al. 1984). The diet of alligator gar residing in estuarine systems is un known, but studies of populations in freshwater suggest a primarily piscivorous diet consisting of freshwater and estuarine fishes (García de León et al. 2001, Robertson et al. 2008.
Estuaries host diverse assemblages of potential prey, but some are particularly important due to their high densities and broad distributions throughout the ecosystem. Four of the most common species of potential prey that co-occur with the focal predators in this study are Gulf menhaden Brevoortia patronus, At lan tic croaker Micropogonias undulatus, brown shrimp Farfantepenaeus aztecus, and white shrimp Lito penaeus setiferus. These species are highly abun -dant in Gulf of Mexico (GOM) estuaries and are commonly found in the diets of estuarine predators (Akin & Winemiller 2006), so they are the most likely prey items to be shared by the predators in this study. These prey species have annual life cycles, with adults of each species moving to coastal waters to spawn in the fall (Atlantic croaker, Gulf menhaden, and brown shrimp) or spring (white shrimp) (Muncy 1984, Cowan & Shaw 1988, Rogers et al. 1993, Brown-Peterson et al. 2017. Using these 8 species as model predators and prey in a northern GOM estuary, the objectives of this study were to (1) create habitat suitability models and maps for each predator and prey species using environmental and spatiotemporal predictor variables, (2) delineate the extent of highly suitable habitat (HSH) for each species on a seasonal basis, and (3) determine the degree of overlap in HSH between predators and prey. A random forest modeling ap proach was applied to species presence/absence and environmental data from a combination of gear types (gillnet, seine, and trawl) to effectively address each of the aforementioned objectives. The results of this study will provide insights into the extent and degree of spatio temporal overlap between predators and prey across 2 seasons (spring and fall) and 3 decades (1986− 2018), which elucidates when and where ecological interactions between them are most likely to occur.

Study area
This study was conducted in the GBC, Texas, in the northwestern GOM from 1986−2018 (Fig. 1). The GBC consists of 5 major sub-bays, each with unique physicochemical conditions. The major source of freshwater into the GBC is the Trinity River, which enters the estuary in the northeast region, Trinity Bay. Barrier islands separate the GBC from the GOM, and 2 tidal inlets allow for transfer of seawater at the Houston Ship Channel (Central Bay) and San Luis Pass (West Bay). The majority of the GBC is shallow, with an average depth of 2 m (Orlando et al. 1993), with the notable exception of the Houston Ship Channel (12 m depth; USACE 2019) which bisects the GBC along its longitudinal axis.

Field data collection
To examine the habitat associations and spatiotemporal overlap patterns of predators and prey in the GBC, habitat suitability models were created using data from the Texas Parks and Wildlife Department Coastal Fisheries Division's long-term Marine Resource Monitoring Program (Martinez-Andrade et al. 2009). A stratified cluster sampling design was implemented, with each major bay system along the Texas coast treated as a stratum. In each stratum, 20 bag seines and 20 trawls were conducted per month (sampling period: 1 mo), and 45 gillnets were conducted per season (sampling period: 1 season; spring or fall). Sampling locations were chosen from grids (one minute latitude by one minute longitude) and were drawn independently and without replacement for each sampling period. The spring season for gill net sampling began the second full week of April and continued for 10 wk, while the fall season began the second full week of September and continued for 10 wk. Gillnet sets occurred overnight, starting within 1 h before sunset and ending within 4 h after sunrise. water quality variables, including temperature, salinity, dissolved oxygen, and turbidity were collected at the deep end of the gillnet immediately following deployment. Water depth was recorded at the shallow end of the net (near shore) and the deep end (away from shore). Bag seine and trawl sampling occurred yearround, with half of the monthly samples collected during the first 2 wk of each month and the other half collected during the last 2 wk. No grid was sampled more than once per month. Bag seine and trawl sampling occurred during the day, between 30 min before sunrise and 30 min after sunset. Bag seines were 18.3 m long, 12.2 m wide (ensured using a limit line), and 1.8 m in height, with 19 mm stretched nylon mesh in the wings and 13 mm stretched nylon mesh in the bag. Bag seine tows were 15.2 m long and were conducted parallel to the shoreline. Surface abiotic water quality variables (temperature, salinity, dissolved oxygen, and turbidity) were collected 3.1 m from shore in the sampling area. Water depth was recorded at the shoreline end of the bag seine (shallowest) and the offshore end (deepest).
While gillnet and seine sampling were restricted to grids that include the shoreline, trawl sampling occurred throughout the middle regions of the bay. Catch data from trawl samples were not included in this analysis, but environmental variables collected at the time of sampling were used. Abiotic water quality variables, including temperature, salinity, dissolved oxygen, and turbidity, were collected from 0.3 m off the bottom before trawling began. The shallowest and deepest water depths encountered during the tow were recorded.

Data analysis
The number of individuals of each species collected per sample (per gillnet set or per seine tow) was calculated for gillnet and seine catch data and then transformed to presence/absence. Gillnet catch data were used for analysis of predators, while seine catch data were used for prey. Gillnets were generally selective for larger, more mobile species such as the predators included in this study, and the seines captured all prey species effectively in similar locations as gillnets (along shorelines). Catch data from trawl sampling were not used for creating habitat suitability models of prey because of the lack of spatial overlap between gillnet and trawl sampling, which would reduce comparability between predator and prey models. Furthermore, seines captured a similar size range of prey compared to trawls (i.e. representative size ranges of all prey species), so seine data was the optimal choice for prey habitat suitability modeling. Abiotic water data from trawl sampling were used for interpolating predictions across the GBC. For all data sets, a season factor was created, with April, May, and June being designated as spring and September, October, and November as fall. A strata factor was also created to be used as a spatial variable in future habitat suitability modeling. Five strata of approximately equal size (250− 300 km 2 ) were created in ArcGIS (ESRI) to represent the 5 major sub-bays within the GBC and were linked with all samples conducted within each stratum (Fig. 1). Decade (1980s, 1990s, 2000s, 2010s) was also calculated as a factor to examine long-term temporal effects. Finally, data were censored to include only samples collected after 1985, since collection methods were not ubiquitously standardized until 1986.
Seasonal (fall and spring) habitat suitability models were created for the 4 predator species (spotted seatrout, black drum, bull sharks, and alligator gar) and 4 prey species (Atlantic croaker, Gulf menhaden, brown shrimp, and white shrimp). Random forest modeling was chosen for this analysis, as it has been shown to perform well when using large data sets with many predictor variables to model species distribution and habitat suitability (Lawler et al. 2006). The basis of random forest modeling is the classification tree, which is a nonparametric model that recursively splits a binary response variable into groups, as homogeneously as possible, based on predictor variables. The most commonly used threshold value for determining classification is 0.5, meaning the class with the majority of votes 'wins' for a given observation (Liaw & Wiener 2002). Random forest models build an ensemble of classification trees, with each tree using a randomly selected subset of data points and predictor variables, and calculate predictions using a model-averaging method (see Breiman 2001 for additional detail). This method decouples the effects of predictor variables on the response, which eliminates the need to conduct variable selection and is effective for data sets with many potentially interactive (collinear) predictors.
Two random forest models were created for each of the 8 species: one using data from the spring and one using data from the fall (total of 16 individual models). The response variable for each model was the presence/ absence of that species, with gillnet data used for predators and seine data used for prey. The predictor variables used in each model were water temperature (°C), salinity (psu), dissolved oxygen (mg l −1 ), turbidity (NTU), shallow water depth (m), deep water depth (m), decade (categorical), and stratum (categorical). No transformations were necessary for predictor variables, since random forest models make no distributional assumptions (Cutler et al. 2007). Additionally, no time-lagged predictor variables were included, so the models assume that the conditions at the time of sampling predict habitat use sufficiently. While time-lagged predictor variables are useful in some marine species distribution models (e.g. Olden & Neff 2001, Wang et al. 2018, it is likely that abiotic conditions in shallow subtropical estuaries like Galveston Bay change too rapidly and frequently for this method to apply. For each model, 500 unpruned classification trees were built, each using a random bootstrapped selection of data points with replacement. In a typical bootstrap sample, approximately 63% of the original data points occur at least once in the selection (Cutler et al. 2007). Each tree used 3 of the 9 predictor variables (e.g. the square root of the number of predictors), which is a common method in random forest modeling (Liaw & Wiener 2002, Cutler et al. 2007). The default threshold value of 0.5 was used for all models (Liaw & Wiener 2002).
Data points that are not used in the building of each tree, the out-of-bag (OOB) data, are used as a testing data set for that tree iteration (approximately 37% of the original data points each iteration). This process allows for the calculation of an error estimate by running the OOB data through the created tree and comparing the proportion of predicted classifications to the true proportions. This calculation provides an estimate of the model prediction error, similar to the use of cross-validation. Five accuracy metrics were calculated for each model: the overall average error rate, sensitivity (proportion of presences correctly classified), specificity (proportion of absences correctly classified), the true skill statistic (TSS), and the area under the receiver operating characteristic curve (AUC). The TSS is calculated as follows: (sensitivity + specificity) − 1. It is an assessment of a model's predictive capacity which incorporates commission (false positive) and omission (false negative) errors and is insensitive to species prevalence. TSS values near 1 indicate nearly perfect model performance, while values near 0 indicate the model performs similarly to random chance (Allouche et al. 2006). Similarly, AUC provides an independent measure of each model's accuracy, ac counting for both omission and commission errors by plotting the relationship between the true positive rate (sensitivity) and the false positive rate (1 − specificity) and calculating the area under the curve. An AUC value of 1 indicates nearly perfect model accuracy, while an AUC value of 0.5 represents random chance model performance (Fielding & Bell 1997). To examine the importance of each variable in predicting the response, a variable importance metric is calculated using the Gini method. When node splitting occurs in a classification tree, there is a Gini impurity criterion for each node which is always smaller than that of the parent node. The sum of the decrease in the Gini index for each variable across all 500 trees is the variable importance (Breiman et al. 1984). To visualize the variable importance (see Fig. 2), the Gini metric was extracted for each predictor variable in each model. That value was then divided by the maximum value in that model, such that the most important variable was given a value of 1 and all others were relative to the most important. All random forest modeling and calculation of accuracy metrics were conducted in R v.3.5.1 (R Core Team 2018) using the packages 'randomForest' (Liaw & Wiener 2002) and 'pROC' (Robin et al. 2011).
Once all habitat suitability models were constructed, the probability of presence of each species in each season was predicted based on previously unseen abiotic data that had not been used during model building. For the predators, this included abiotic data from trawl and seine sampling, while for the prey it included abiotic data from trawl and gillnet sampling. Using these new abiotic data points allowed for prediction of animal presence across the entire bay system, since trawl sampling occurred in the central open-water regions of the GBC and seine and gillnet sampling occurred exclusively along the shorelines. Interpolation of habitat use in open-water regions using models built with data exclusive to shorelines may introduce some uncertainty to these maps, primarily concerning differences in depth and access to structured habitats. However, these models apply to open-water reasonably well since they do not include structured habitat variables and are thus based primarily on abiotic water quality variables. To create seasonal habitat suitability maps for each species, these predicted values were plotted in ArcGIS (ArcMap v.10.7.1) and projected to a standard coordinate system for spatial consistency (NAD 1983 UTM Zone 15N). Inverse distance weighted interpolation was then used to create smoothed maps using a cell size of 100 × 100 m, a power of 1, and a fixed search radius with a minimum of 15 points. The final habitat suitability maps (2 for each species, total of 16) were then visualized using 10 bins of increasing probability of presence, each greater by 0.1 (see Figs. 3 & 4).
To delineate the extent of HSH for each species in each season, all cells in each habitat suitability map were extracted that contained values greater than the average probability of presence of that species across both seasons ([average probability of presence in fall + average probability of presence in spring] / 2). This created a species-specific cutoff value for HSH, which was necessary considering the wide range in probability of presence across species. The extent of HSH was mapped and calculated for each species in each season. Finally, to determine the amount of overlap in HSH between predators and prey, overlay analysis was used to extract any areas where HSH of both species occurred. The area of that overlap was calculated for each iteration. This overlap analysis was conducted for all predator and prey combinations in each season, resulting in 8 overlap maps per predator (see Figs. 5, 6, 7 & 8, Table 3).

Summary of data
The data set used to model predator presence/ absence consisted of 2882 individual gillnet sets over 33 yr (1986−2018) The average values of each abiotic variable used in the modeling and mapping procedures are detailed in Table 1. On average, temperatures were higher in the spring than the fall in all strata, and there was a northeast to southwest gradient of temperature from Trinity Bay (lowest) to West Bay (highest). Salinity was, on average, higher in the fall than the spring in all strata and exhibited a similar gradient from Trinity Bay (lowest) to West Bay (highest). Dissolved oxygen was slightly higher in the fall than spring and was highest in Trinity Bay and lowest in West Bay in both seasons. Turbidity was higher in the spring than fall and was highest in East Bay and Trinity Bay in both seasons. Water depths (shallow and deep) were generally shallowest in West Bay, followed by East Bay, Central Bay, Trinity Bay, and Northwest Bay (deepest) in both seasons.

Seasonal habitat suitability
Accuracy metrics for each habitat suitability model are detailed in Table 2 However, the fall model had far lower specificity (0.02) and lower TSS (0.02) and AUC (0.59) values than the spring model (specificity = 0.31, TSS = 0.27, AUC = 0.79), meaning absences were predicted less accurately in the fall than the spring. This likely arose because spotted seatrout were absent in fewer gill-nets in the fall (absent in 12% of gillnet sets) than in the spring (15%). Overall, spotted seatrout were very common, resulting in high sensitivity of both models (0.97 and 0.99 in the spring and fall, respectively). The most important variable influencing the probability of spotted seatrout presence differed between the spring (salinity) and the fall (temperature; Fig. 2). Salinity was far more important than all other variables in the spring (i.e. dominant), but temperature was only slightly more important than other variables in the fall (i.e. non-dominant). The seasonal habitat suitability maps for spotted seatrout reflect their abundance and ubiquity throughout the GBC, with a high probability of presence in all strata in both seasons (Fig. 3A,B). However, the probability of presence in the spring was reduced in Trinity Bay along the shoreline closest to the mouth of the Trinity River (Fig. 3A). Similar to spotted seatrout, black drum had a lower rate of absences in the fall (absent in 16% of gillnets) than in the spring (21%), resulting in lower specificity in the fall (0.09) than the spring (0.18). The average error rate for black drum was slightly higher in the spring model (21%) than the fall model (16%), and the AUC was lower in the fall (0.66) than the spring (0.70). The TSS was also lower in the fall (0.07) than the spring (0.13). Since black drum were common in the gillnet catch data, both models had high sensitivity (0.96 and 0.98 in the spring and fall, respectively). The most important variable predicting the presence of black drum was salinity in the spring and temperature in the fall (Fig. 2). However, differences in variable importance scores among most variables were minor in both the spring and in the fall, suggesting the absence of a single dominant variable in either season. The seasonal habitat suitability maps show a reduction in black drum presence in the spring compared to the fall, especially in the Trinity Bay and Northwest Bay strata (Fig. 3C,D).
The spring model for bull sharks had a higher average error rate (18%), a lower AUC (0.81), but a higher TSS (0.34) compared to the fall model (error = 14%, AUC = 0.83, TSS = 0.31). Due to the high rate of absences of bull sharks in this data set (absent in 78 and 82% of gillnets in the spring and fall, respectively), the specificity (spring = 0.94, fall = 0.96) was higher than the sensitivity (spring = 0.40, fall = 0.35) for both seasons. Temperature was the most important variable in both seasons and had a far higher importance score compared to the next-highest, salinity (Fig. 2). The probability of bull shark presence was generally low throughout the GBC but was highest in the Trinity Bay and East Bay strata  Table 2. Accuracy metrics from seasonal random forest models for predators and prey. The average error rate of the model is based on out-of-bag (testing) data, sensitivity is the proportion of presences correctly classified, and specificity is the proportion of absences correctly classified. True skill statistic (TSS = sensitivity + specificity − 1) and area under the receiver operating characteristic curve (AUC = relationship between true positive rate and false positive rate) are measures of model predictive accuracy (Fig. 3E,F). Minor seasonal differences in bull shark presence were observed, with higher probability of presence in the Northwest stratum in the spring compared to the fall, and higher probability of presence in the West Bay stratum in the fall compared to the spring (Fig. 3E,F). The spring model for alligator gar had a higher average error rate (26%) and a lower AUC (0.78) compared to the fall model (error = 21%, AUC = 0.81). The TSS value was higher in the spring (0.40) than the fall (0.36). There was a moderately high rate of alligator gar absences (absent in 64 and 73% of gillnets in the spring and fall, respectively), resulting in higher specificity (spring = 0.87, fall = 0.91) than sensitivity (spring = 0.53, fall = 0.45) for both seasons. Salinity was the most important variable predicting alligator gar presence in both seasons, and the second most important variable, temperature, had a far lower importance score (Fig. 2). Similar to bull sharks, the probability of alligator gar presence was highest in the Trinity Bay and East Bay strata across both seasons but was higher in the Northwest and West Bay strata during the spring compared to the fall (Fig. 3G,H).
Atlantic croaker were far more abundant in the spring (present in 70% of seine pulls) than the fall (24%), resulting in higher sensitivity (0.89) than specificity (0.39) in the spring and a lower sensitivity (0.16) than specificity (0.91) in the fall. The average error rate was similar across seasons (26% in the spring and 25% in the fall), but the AUC and TSS values were slightly lower in the fall (AUC = 0.67, TSS = 0.10) than the spring (AUC = 0.73, TSS = 0.28). Temperature was the most important variable influencing Atlantic croaker presence in both seasons, but it was not dominant (i.e. other variables had similar importance scores; Fig. 2). The highest probability of Atlantic croaker presence occurred in the Trinity Bay and East Bay strata in both seasons, but this probability was much higher throughout the GBC in the spring compared to the fall (Fig. 4A,B).
Similar to Atlantic croaker, Gulf menhaden were more abundant in the spring (present in 51% of seine pulls) than the fall (11%) and had slightly higher sensitivity (0.70) than specificity (0.65) in the spring and far lower sensitivity (0.08) than specificity (0.98) in the fall. The average error rate was higher in the spring (32%) than the fall (11%), but the AUC was similar across seasons (0.75 in the spring and 0.76 in the fall). However, the TSS was much higher in the spring (0.35) than the fall (0.07). With almost equal importance scores, temperature and deepest depth were the most important variables predicting Gulf menhaden presence in the spring. Salinity was the most important variable in the fall, but it was not dominant (Fig. 2). Gulf Menhaden probability of presence was highest in open-water regions throughout the GBC and was generally much higher in the spring than the fall (Fig. 4C,D).
Brown shrimp were more abundant in the spring (present in 79% of seine pulls) than the fall (43%), resulting in higher sensitivity (0.94) than specificity (0.42) in the spring and lower sensitivity (0.45) than specificity (0.70) in the fall. The spring model had a Temperature was the most important variable for both seasons, but it was not dominant (Fig. 2). Brown shrimp probability of presence was highest in the West Bay and East Bay strata in both seasons, and this probability was generally higher in the spring than the fall (Fig. 4E,F). In contrast to the other prey species, white shrimp were more abundant in the fall (present in 79% of seine pulls) than the spring (19%), and the sensitivity (0.09) was lower than the specificity (0.98) in the spring while the sensitivity (0.95) was higher than the specificity (0.25) in the fall. The average error rate was equal between the seasons at 19%, but the AUC and TSS values were lower in the spring (AUC = 0.70, TSS = 0.07) than the fall (AUC = 0.75, TSS = 0.20). The most important variable influencing white shrimp presence was temperature in both seasons, but it was not dominant (Fig. 2). The probability of presence of this species was generally low throughout the GBC in the spring, but slightly higher in the Trinity Bay stratum. In the fall, high probability of white shrimp presence was observed throughout the GBC, but it was lowest in the Central Bay stratum (Fig. 4G,H).

Spatiotemporal overlap between predators and prey
Areal extent of spatiotemporal overlap for all predator− prey combinations are detailed in Table 3. Spotted seatrout exhibited extensive overlap (de fined as > 300 km 2 ) in HSH with Atlantic croaker, Gulf menhaden, and brown shrimp in all strata but Trinity Bay in the spring (Fig. 5A,C,E), but had very little overlap (defined as < 50 km 2 ) with those prey species in the fall (Fig. 5B,D,F). In contrast, extensive overlap occurred between spotted seatrout and white shrimp HSH in the fall throughout the GBC (Fig. 5H), but no overlap between them occurred in the spring (Fig. 5G).
HSH for black drum overlapped extensively with that of Atlantic croaker, Gulf menhaden, and brown shrimp in the spring, primarily along the coastlines and in some open water regions of the East, West, and Central Bay strata (Fig. 6A,C,E). No overlap be tween black drum and white shrimp HSH was ob served in the spring (Fig. 6G) Table 3. Area (km 2 ) of overlap in highly suitable habitat between predators and prey by season overlap between black drum and Atlantic croaker, Gulf menhaden, or brown shrimp HSH was observed ( Fig. 6B,D,F). Extensive overlap in HSH between black drum and white shrimp was ob served throughout the GBC in the fall, primarily in East and West Bay and the shorelines of the other strata (Fig. 6H). Extensive overlap in HSH between bull sharks and Atlantic croaker, Gulf menhaden, and brown shrimp occurred in all strata in the spring, but most ubiquitously in Trinity and East Bay (Fig. 7A,C,E). Little overlap in HSH was observed between bull sharks and those prey species in the fall, with the most occurring between bull sharks and brown shrimp in the East and West Bay strata (Fig. 7B,D,F). No overlap in HSH occurred between bull sharks and white shrimp in the spring (Fig. 7G), but extensive overlap occurred between them in the fall in all strata except the Northwest Bay (Fig. 7H).
Similar to bull sharks, extensive overlap in HSH between alligator gar and Atlantic croaker, Gulf menhaden, and brown shrimp occurred in the spring in Trinity and East Bay, but also in the Northwest and West Bay strata (Fig. 8A,C,E). Very little overlap was observed between Alligator Gar and those species in the fall, with most overlap occurring with brown shrimp in the open-water areas of East Bay (Fig. 8B,D,F). No overlap in HSH occurred between Alli gator gar and white shrimp in the spring (Fig. 8G), but extensive overlap was observed between them in the Trinity and East Bay strata in the fall (Fig. 8H).

DISCUSSION
The seasonal habitat suitability models and maps produced in this study exemplify the variability in habitat use patterns among predators and prey in the GBC and delineate the extent of spatiotemporal overlap between them. The influence of environmental variables on the distribution of each species differed, but temperature or salinity were almost always the most important. Some predator species were highly abundant and had moderate to highly suitable habitat throughout the GBC, such as spotted seatrout and black drum, while others were less abundant and more restricted in the extent of their suitable habitat, such as bull sharks and alligator gar. Seasonal shifts in habitat suitability occurred for all predator and prey species, but this pattern was more notable for prey than for predators. Overlap in HSH between predators and prey occurred primarily during whichever season the prey species was most abundant and in regions with extensive HSH for the predator. Determining the ideal environmental conditions and degree of overlap between these predators and prey allows for inference regarding the biotic and ecological interactions that shape this estuarine ecosystem.
When sympatric predators overlap spatiotemporally, interactions between them can affect their ecological roles through behavioral alterations and changes to local predation pressure (Sih et al. 1998. Determining the habitats, regions, or environmental conditions that support multiple sympatric predators can thus have implications for ecosystem function. Habitat suitability maps of predators in the GBC revealed major differences in the abundance and distribution of the 4 study species. Spotted seatrout were the most abundant and ubiquitous, with high habitat suitability throughout the system in both seasons, followed by black drum, which were slightly less abundant and exhibited clear regions where habitat suitability was reduced. Alligator gar and bull sharks were the least likely to be present throughout the system, and their extent of HSH appeared to be restricted primarily to the Trinity Bay and East Bay strata across both seasons. This spatial separation of bull sharks and alligator gar to specific regions of the estuary suggests that the distribution of ecological dynamics, such as predation pressure on shared prey, is not even across the ecosystem. In areas where all 4 predators experience high habitat suitability, such as East Bay and Trinity Bay, interactions between species and conspecifics may play a larger role in shaping ecosystem function compared to areas where only 1 or 2 predators are abundant, such as West Bay. In river-dominated estuaries like the GBC, salinity is often an important factor driving community structure since it can vary dramatically based primarily on river inflow and tidal mixing (Jones et al. 1990, Ley et al. 1999, Jenkins et al. 2015. The Trinity Bay and East Bay strata had the lowest average salinity levels (see Table 1) due to inflow of freshwater from the Trinity River in the northeast corner of the GBC. While spotted seatrout were relatively abundant throughout the system, they experienced lower habitat suitability in the Trinity Bay stratum in the spring. Salinity was the dominant variable with the highest importance score for spotted seatrout in the spring but was less important in the fall. This slight seasonal difference in habitat use likely arose due to the reduced salinity in that region in the spring compared to the fall. With higher amounts of freshwater entering the system through the Trinity River in the spring, the probability of presence of spotted seatrout was reduced in that area, suggesting movement away from low-salinity waters. Similar behavior occurs in Louisiana estuaries, with spotted seatrout actively retreating from low-salinity regions after the release of fresh water into the system (Callihan et al. 2015). On an oyster reef in a Texas estuary, multiple adult spotted seatrout were observed migrating synchronously away from the reef during a period of extremely low salinity (TinHan et al. 2018). Reproductive activities may also contribute to the spatial distribution of spotted seatrout within estuaries, with adults generally spawning in higher salinity regions in the spring and summer (Helser et al. 1993, Saucier & Baltz 1993. The reduction in habitat suitability of adult spotted seatrout in the low-salinity regions of the GBC during the spawning season suggests environmental and biological mechanisms are likely both driving the species' distribution. Black drum were less ubiquitous in the GBC compared to spotted seatrout and exhibited more restricted areas of high habitat suitability. While water depth was not the most important variable in habitat suitability models, shallow areas close to structure and shorelines appeared to harbor high probabilities of black drum presence, which may be indicative of the feeding strategy of this species. Black drum primarily consume benthic invertebrates such as oysters, shrimps, crabs, polychaetes, and bivalves (Overstreet & Heard 1982, Rubio et al. 2018, which are generally abundant and accessible in shallow waters. Other variables, including temperature and salinity, influenced the distribution and abundance of black drum, but the effect of these variables was not obvious and differed seasonally. Black drum were caught less frequently in the spring than the fall, resulting in lower probability of presence, and thus less suitable habitat in the GBC in the spring. Adult black drum, once sexually mature, move from within estuaries to the coastal zones and continental shelf to reproduce during the spring (Nieland & Wilson 1993). The movement of adults out of the estuary to spawn was most likely the major factor driving reduced probability of presence in the spring, not necessarily a reduction in habitat suitability due to changing environmental variables.
Ontogenetic changes in habitat-use patterns, such as movements related to reproductive status, are common in coastal and estuarine fishes (e.g. Mellin et al. 2007, Knip et al. 2011. While bull sharks are euryhaline and can be found in estuaries, coasts, and open-ocean systems, they use freshwater regions of estuaries as juvenile nursery grounds while they grow and develop . The average total length of bull sharks collected by gillnets in this data set was 1013 mm in the fall and 1055 mm in the spring, which is well below recorded lengths at maturity for both males and females (Branstetter & Stiles 1987), meaning most bull sharks in this study were juveniles or subadults. The observed restriction of their highest habitat suitability to the Trinity and East Bay strata suggests that the bull sharks in this study primarily remain in low-salinity waters in the GBC, most likely to reduce predation by larger predators that cannot tolerate freshwater exposure. However, the most important variable influencing the distribution and abundance of bull sharks was temperature, which exhibited a similar pattern to salinity and was lowest in the Trinity and East Bay strata. This combined temperature− salinity influence was likely responsible for seasonal differences in habitat use. In the fall, bull shark probability of presence was slightly increased in West Bay, despite salinity being the highest observed across all seasons and strata. The average temperature (25.05°C) was similar to Trinity and East Bay temperatures in the spring (25.83 and 25.91°C, respectively), suggesting there may be a temperature optimum for juvenile bull sharks in the GBC around 25°C, which supersedes the effect of salinity. A previous study found evidence for a strong, negative impact on juvenile bull sharks during a short period of very cold temperatures (Matich & Heithaus 2012). It is possible that the episodic freshwater inflow at the base of the Trinity River in the fall creates prohibitively cold conditions, and juvenile bull sharks may find refuge in the warmer West Bay region.
The most freshwater-associated of these study species, alligator gar, is generally found within rivers and streams, and the results of this study suggest that the species primarily inhabits regions within the GBC estuary that are generally dominated by low salinities. Furthermore, the most important variable in habitat suitability models was salinity, so the distribution of alligator gar in river-dominated estuaries like the GBC is dominated by the effect of freshwater inflow from the major river. Seasonal differences were observed, whereby alligator gar probability of presence was higher throughout the GBC in the spring, when salinities were generally lower in all strata compared to the fall. The expansion of alligator gar into regions of the GBC that were farther from the freshwater source suggests that more expansive low-salinity waters allowed for an increase in ideal habitat throughout the ecosystem. Previous work has suggested an estuarine subpopulation of alligator gar that remains near the base of the Trinity River (Buckmeier et al. 2013), but this study only included individuals captured within the river. While our results provide evidence for a large number of alligator gar within the GBC and seasonal expansion of suitable habitat throughout the entire ecosystem, further re search is required to determine the residency of estuarine individuals in comparison to their riverine conspecifics.
While all predator species exhibited some degree of seasonal change in habitat suitability, this seasonality was more evident for the prey species as dramatic differences in their probability of presence were observed between the spring and fall. Seasonal shifts in prey presence were likely due to the timing and nature of their reproductive strategies, with adults of each species moving to coastal waters to spawn in the fall (for Atlantic croaker, Gulf menhaden, brown shrimp) or spring (for white shrimp) (Muncy 1984, Cowan & Shaw 1988, Rogers et al. 1993, Brown-Peterson et al. 2017). Once they have reproduced, adults move back into the estuary for the remainder of the year, and the juveniles they produced will recruit into the population the following season. This annual inshore−offshore life cycle re sults in very high abundance during periods of growth and development and sharply reduced abundance during periods of reproduction. Because these prey species were chosen specifically due to their high abundance throughout the GBC, the probability of presence for each species was generally high in all strata, especially in whichever season they were more abundant.
Spatiotemporal overlap of HSH between predators and prey occurred most extensively during the spring for all prey but white shrimp. Because each prey species was highly abundant and widely distributed during a single season in the GBC, the spatial extent of overlap in HSH was dictated primarily by the habitat suitability of the predators. For example, overlap between Atlantic croaker and black drum occurred almost exclusively in the spring (when Atlantic croaker were abundant) and along the shorelines of East, Central, and West Bay (where black drum probability of presence was highest). Therefore, the temporal distribution of overlap was determined by the seasonality of the prey species while the spatial distribution of overlap was determined by the habitat requirements of the predator. If we assume that overlap in HSH between predators and prey results in a high rate of physical contact be -tween them, these results provide the spatial and temporal conditions where ecological interactions, such as predation, are most likely to occur between each combination of these predator and prey species.
This analysis of spatiotemporal overlap in HSH between predators and prey suggests that the availability of potential prey differs significantly between spring and fall. While there are many other potential prey species in the GBC (other small finfish, benthic invertebrates, etc.), the 4 prey species chosen for this study are the most abundant and widespread based on the catch data used in this analysis. Seasonal changes in prey availability may influence the diets of these predators, which may alter vital rates such as growth and mortality on a seasonal basis and influence population dynamics (Zhou et al. 2013). Seasonality in resource availability has been linked to changes in the diets and health of estuarine fishes, including black drum and bull sharks. In an estuary south of the GBC along the Texas coast, Baffin Bay, emaciated Black drum were observed during a period of intense drought and hypersaline conditions during the fall of 2012 (Olsen 2016). The prolonged hypersalinity in the estuary was linked to reductions in food availability, as well as limited movement of black drum related to foraging activities , Breaux et al. 2019). In the Florida Everglades, bull sharks can adjust their diets during periods of drought (i.e. dry seasons), when prey from freshwater marshes are pulsed into the upstream regions of the estuary (Matich & Heithaus 2014). Investigating the seasonality of predator diets in the GBC would complement the results of this study by determining whether a seasonal difference in spatiotemporal overlap between predators and prey translates to seasonal changes in diet and condition.
The results of this study enhance our understanding of the factors shaping the distribution and habitat suitability of predators and prey in the GBC and suggest when and where interactions between them are most likely to occur. However, the extent and distribution of suitable habitat for these species is likely to be affected by future climate scenarios (Fujiwara et al. 2019). Rainfall events are predicted to become less frequent but more intense in many regions (Pachauri et al. 2014), which would alter environmental conditions such as salinity and temperature in the GBC and other river-dominated estuaries. Since these environmental variables, and their seasonality, are critical to determining the extent of HSH and overlap between predators and prey, any alterations via anthropogenic disturbances are likely to disrupt ecosystem functioning. Future management strate-gies should incorporate variability in climatic and environmental conditions when considering best practices for protecting commercially and recreationally exploited species. Further research investigating the spatiotemporal overlap between predators and prey in coastal ecosystems, and the environmental drivers behind those habitat use patterns, should enhance our understanding of how changing climatic conditions will alter ecological dynamics and the roles of predators and prey in their shared ecosystem.