Environmental and anthropogenic influences on spatiotemporal dynamics of <fi>Alosa</fi> in Chesapeake Bay tributaries

American Shad (Alosa sapidissima), Hickory Shad (A. mediocris), and river herrings (Alewife A. pseudoharengus and Blueback Herring A. aestivalis) are anadromous pelagic fishes, which as adults spend most of the annual cycle at sea, but enter the coastal rivers in spring to spawn. Once as one of the most valuable fisheries along the Atlantic coast, Alosa populations have declined in recent decades and current populations are at historic lows. Various management actions have been conducted to restore the populations, and stocks in different river systems display different demographic trends. Demonstration of synthetic diagnostics on the factors impacting these populations is important to better conserve this species group. We developed a Bayesian hierarchical spatiotemporal model to identify the population trends of these species among rivers in the Chesapeake Bay based on results of surveys conducted by the Virginia Department of Game and Inland Fisheries and Maryland Department of Natural Resources and to identify environmental and anthropogenic factors influencing their distribution and abundance. The hierarchical model structure helped to diagnose river-specific population trends and impacts of surrounding factors, and decrease uncertainties in rivers with less samples available. The results demonstrate river-specific heterogeneity of spatiotemporal dynamics of these species and indicate river-specific impacts of multiple factors, including water temperature, river flow, chlorophyll a concentration, and total phosphorus concentration, on their population dynamics. Atlantic Multidecadal Oscillation and Gulf Stream meanders displayed significant influence on the inter-annual trends of Alosa species in rivers with more data available. The results would help to develop riverand species-specific management strategies to recover these spe-

INTRODUCTION adults, these fishes spend most of the annual cycle at sea, but enter freshwater in spring to spawn (Turner and Limburg 2016). Emigration from natal rivers to oceanic waters occurs during the first year of life (Turner and Limburg 2016). These Alosa species serve as important forage fishes for anadromous and marine predators, such as Striped Bass (Morone saxatilis), Spiny Dogfish (Squalus acanthias), Atlantic Cod (Gadus morhua), and Pollock (Pollachius virens), and seaward-emigrating young-of-year fish encounter a gauntlet of marine predators (Smith and Link 2010). Invasive predators such as Blue Catfish (Ictalurus furcatus) and Flathead Catfish (Pylodictis olivaris) also may impact Alosa populations (Schmitt et al. 2017). Hence, these species interactions present a clear trophic link between inland and marine production (MacAvoy et al. 2000, Walters et al. 2009).
American Shad, known as America's "founding fish," and river herring once supported the largest and most important commercial and recreational fisheries along the Atlantic coast (McPhee 2003). Since colonial times, the blockage of spawning rivers by dams and other structures (Perley 1852, Kosa and Mather 2001, Hall et al. 2011, as well as habitat degradation, climate change, overfishing, and bycatch (Köster et al. 2007, Limburg and Waldman 2009, Hasselman et al. 2015, have severely depleted their populations (ASMFC 2007(ASMFC , 2012. The coastwide decline of Alosa stocks culminated in the Atlantic States Marine Fisheries Commission (ASMFC) requiring moratoria on alosid fisheries unless a state or jurisdiction developed a sustainable fishery management plan (ASMFC 2009(ASMFC , 2010. States and jurisdictions are also required to implement fisheries-dependent and independent monitoring programs to outline spawning and nursery habitat, identify local threats to habitat, and develop plans for mitigation and restoration in each river system (ASMFC 2009(ASMFC , 2010. Two most serious threats identified include barriers to migration and climate change (Capossela 2013, Hilton et al. 2014. Fish passage efforts, including dam removal and fishway construction, have continued on the east coast as just one action to try to reverse the decline in Alosa populations (Chesapeake Bay Program 2019, Martin 2019. Between 1989 and 2017, a total of 3746 stream miles have been reopened to the migration of fish in the Chesapeake Bay watershed (Chesapeake Bay Program 2020).
A deep understanding of the critical environmental and anthropogenic factors influencing populations through space and time in river systems helps to identify river-specific threats and aid in future management. Complex life history, combined with geographic differences in population dynamics, complicates assessment of these species on a coastwide scale (ASMFC 2007(ASMFC , 2012. Ideally, they should be assessed and managed by incorporating river-specific heterogeneity (ASMFC 2012), which is likely due to riverspecific environmental factors, including anthropogenic structures and activities. Different recovery trends were observed for river-specific stocks with data available (ASMFC 2012(ASMFC , 2017. Ogburn et al. (2017a) also argued for including heterogeneity in Alosa management in the Chesapeake Bay based on differences in genetic and juvenile abundance data. Given the wide distribution of shads and river herrings along the Atlantic coast, a bay-scale study to diagnose the factors that influence spatiotemporal dynamics of these populations, and a study to develop models that consider spatial complexity is critical for management of these populations. Such a study should help to better understand population status at the local scale, provide a river-specific, comprehensive illustration of major threats to Alosa, and develop potential strategies to enhance stock-specific recovery.
The ongoing fisheries-dependent and independent monitoring programs conducted by states and jurisdictions provide opportunities to explore river-specific population trends and identify threats to the local habitat. Analyses on such spatial structured ecological data often require complex model structures (Zuur et al. 2009). Hierarchical models have been increasingly used in analyzing complex nested data (Zuur et al. 2009, Harrison et al. 2018. Hierarchical models can handle complex interactions by allowing parameters to vary at more than one level via an introduction of random effects and perform partial pooling by borrowing information from other groups (Zuur et al. 2009, Harrison et al. 2018. Bayesian approaches are ideally suited for constructing complex hierarchical models, and parameter estimation is straightforward (Gelman et al. 2014a). The advantages of Bayesian hierarchical models emerge as complexity increases, when, for example, spatiotemporal variability needs to be modeled explicitly (Banerjee et al. 2014).
In this study, we developed quantitative models to identify population trends of alosines among the Chesapeake Bay tributaries based on data from surveys conducted by the Virginia Department of Game and Inland Fisheries (VDGIF) and Maryland Department of Natural Resources (MDNR). We identified environmental and anthropogenic factors influencing alosine populations in each river system. We developed models including river as the upper level of hierarchy and species within river as the lower level. The hierarchical model structure helped to diagnose river-specific population trends and impacts of surrounding factors and decrease uncertainties in rivers with less samples available (Zuur et al. 2009, Harrison et al. 2018. Environmental factors considered included both global-scale factors, such as global climate oscillations (Nye et al. 2012, Lynch et al. 2014, and local-scale factors, such as water temperature and river flow. Anthropogenic factors represented by water quality data also were considered (Bilkovic et al. 2002, Madin 2017. The ASMFC organized data workshops to solicit survey data that are credible to be used among states. Here, we followed the list of data available in the Chesapeake Bay region that the ASMFC deemed meaningful in both spatial coverage and survey design. The river herrings spawning in the Chesapeake Bay tributaries are part of the mid-Atlantic stocks of both Alewife and Blueback Herring, which have suffered the most severe declines (Palkovacs et al. 2014). However, due to insufficient data, the Chesapeake Bay region was poorly represented in coastwide assessments (ASMFC 2012), and fisheries were placed under moratoria. Understanding population dynamics of alosines among the Chesapeake Bay tributaries could improve future coastwide stock assessments and management. This study provided a framework for developing hierarchical models to analyze complex ecological problems, and the methods could be applied to other species across different river systems.

Data
The VDGIF has conducted boat electrofishing surveys for alosines in the James, York, and Rappahannock rivers and their tributaries during February to June since 1994 (Fig. 1, ASMFC 2016). Electrofishing surveys use a boat-mounted generator and voltage regulator to put an electric shock into the water to collect fish. The shock affects fish near the boat's electric field that is an area within a few feet around the boat. Electroshocked fish swim toward booms dangling off the front of the boat, where they are temporarily stunned and netted by biologists. The primary objectives of these surveys are to establish a long-term time series of relative abundance indices, to monitor the demographic structure of the spawning runs of American Shad, Hickory Shad, and river herrings, and to relate recruitment indices to relative year-class strength and age structure of spawning adults. Data recorded include catch, effort, environmental factors such as water temperature, dissolved oxygen and salinity, and biological information such as sex, spawning condition, length, weight, and age.
The MDNR has collected fishery-dependent samples of shads and river herrings from commercial pound and fyke nets in the Nanticoke River from 1989 to present (Fig. 1), providing relative indices of abundance for these species. Pound nets are fixed finfish entrapment net devices with a fence leader that interrupts fish movements and a heart that funnels fish into the pound; the nets consist of an arrangement of netting or wire supported upon stakes or piles and have the head ropes or lines above the water or in a frame that is supported by floats and anchors (MSBE 2020). Fyke nets are hoop nets with one or more wings or a leader attached and held in place with anchors or stakes to help guide finfish into the hoop net (MSBE 2020). These surveys typically begin in March and end in April, although there are a few years where they began in February and/or ended in June. Prior to 1997, American Shad and Hickory Shad data were not entered into a digital format, so the data from 1989 to 1996 include only Alewife and Blueback Herring. Data recorded include catch, effort, and gear type, environmental factors such as water temperature, dissolved oxygen, and salinity, and biological information such as sex, length, weight, and age. Due to data format and quality issues, MDNR data only from 1997 to present were used to investigate the population dynamics of the four-species group.
Water quality data from the Chesapeake Bay Program were incorporated into our models to indicate selected anthropogenic influences (Chesapeake Bay Program 2017). The chlorophyll a (chl a), total nitrogen (TN), and total phosphorus (TP) concentrations were obtained from the nearest sampling station on the same sampling date as the VDGIF and MDNR collections.
Daily mean river flow records were obtained from the U.S. Geological Survey (USGS) National Water Information System (NWIS) also at the nearest sampling station on the same sampling date. River flow was not considered for the Nanticoke River in Maryland, because the pound and fyke nets were set for multiple days and the flow data at the time of fish swimming into the net were not available.

Explanatory variable selection
Potential explanatory variables for species abundance in each river that were considered in the respective models are listed in Table 1. The distributions of variables and the correlations between each pair of variables in each river are shown in Fig. 2. There were strong linear correlations (|r| > 0.6, Grewal et al. 2004) between water temperature and dissolved oxygen in all rivers; between water temperature and chl a, between water temperature and TP, between chl a and salinity, between chl a and river flow in the Rappahannock River, and between TN and TP in the There were 46 stations on the James River, 15 stations on the Rappahannock River, 13 stations on the York River, and 44 stations on the Nanticoke River. Some stations were in the tributaries of the respective rivers. In some tributaries, the sample size was too small to distinguish any pattern, so tributaries are included with their main river. Shapefile data are from OpenStreetMap. Plot is made using QGIS.
v www.esajournals.org James and York rivers. The strong correlations between water quality variables in the Rappahannock River were uncertain due to small sample size for chl a, TN, and TP concentrations. Only one variable among water temperature and dissolved oxygen and one among TN and TP concentrations were incorporated into analyses, and otherwise, collinearity could severely distort model estimation (Dormann et al. 2013). Water temperature and TP concentration were recorded for more sampling trips, so these two variables were incorporated into models. Although month and water temperature were correlated with each other, preliminary analyses demonstrated that models with month and water temperature performed better than models with only one of them, so both were incorporated to indicate the effect of temperature and fish behavior. About half of the sampling trips in all rivers lacked salinity records, and 90% of salinity records were within 0 to 0.12, so salinity was not incorporated into the model. Year was included in the analysis to explore inter-annual variation. The range of sampling duration in each river was limited (Table 1), so we used catch-per-trip instead of catch-per-hour as the response variable. Considering the possible effect of fishing effort, sampling duration and sampling method were regarded as explanatory variables. The distributions of sampling duration, chl a, TP, and river flow were skewed positively, so they were normalized through log-transformation. Distance of each sampling location from the corresponding river mouth along the river network was calculated using the R package RIVERDIST (Tyers 2017). RIVERDIST reads river network shapefiles and provides tools for distance calculation along river networks. The river shapefiles were created within QGIS3.10 using OpenStreetMap dataset (https://www.openstreetmap.org).
The four rivers displayed different ranges of water quality variables (Fig. 2). For example, the Nanticoke River had a larger range of chl a (2.85-32.04 μg/L) than other rivers. The Nanticoke River is characterized as a "Low Urban, High Agriculture" river, where nitrogen and phosphorous compounds are added to the river through drainage from surrounding agricultural lands, so that nitrogen and sediment levels are very high in comparison with other rivers in Maryland (MDNR 2014). More ingestible nitrogen and phosphorous may trigger phytoplankton bloom, indicated by a high chl a concentration. Records for chl a and TP were available for only 41% of the sampling trips; therefore, when these water quality factors were taken into analyses, nearly 59% of sampling trips would be filtered out. To obtain more accurate results, we developed a set of models for the full dataset (2642 sampling trips) to detect the spatiotemporal patterns and the effects of environmental factors without regard to the water quality variables. Another set of models was developed and fit to the filtered data (1090 sampling trips) with chl a and TP records Year 1994-2017 1994-2017 2000-2002, 2004-2012, 2014 1997-2001, 2003-2009, 2011-2014, 2016 v www.esajournals.org available to explore the effects of these water quality factors.

Model structure
Details on data used in the models are listed in Table 2. There were high percentages of sampling trips with zero catch for each species in each river. Hurdle (Cragg 1971), and zero-inflated models (Lambert 1992) were developed to analyze data with excessive zeros. Both models can be considered as mixture models, but they interpret zeros differently. Zero-inflated models divide extra zeros into two types: structural and sample zeros. Structural zeros come from reality, while sample zeros come from chance or mistake. Hurdle models consider all zeros as structural zeros. We conducted a preliminary analysis to compare the performances of hurdle and zero-inflated models with different distributions when fitted to the full data and found zero-truncated Negative Binomial hurdle model performed better (see Appendix S1 for details), so we focused on this model in the following analyses.
The basic hurdle model has two sub-models: first, a binomial probability sub-model governs zeros/non-zero observations, and second, a zerotruncated sub-model for the non-zero values. The hurdle model has been widely used to deal with the high percentage of zero observations in fishery data (Li andJiao 2015, Bi et al. 2018). The probability sub-model assumed that the events of capturing or not capturing fish belonging to a specific species (s) in a sampling trip in a specific river (r) followed a binomial distribution with a logit link: where p s,r is the probability of catching a fish belonging to species s in river r; β p s,r is the intercept for each species in each river; x i is the ith explanatory variable. Sampling method, year, and month were modeled as categorical variables, and their effects followed normal distributions with a mean of zero. The effects of other continuous variables were modeled using smoothing functions defined through a first-order random walk (RW1) process (Rue and Held 2005), which assumes that f(X t ) − f (X t−1 )~N(0, τ −1 ); where X is a covariate, t is the tth value of X, τ is a precision parameter. The RW1 process was constrained to sum to zero and rescaled to have typical variance 1 when τ = 1. The intercept, year, month, and other continuous variables were modeled to follow two-level hierarchies with river as the upper level and species within river as the lower level; except for sampling method, which followed one hierarchy, species, because it was a categorical variable and was unique for each state. For details in hierarchical structure, see Appendix S2.
The positive catch sub-model assumed that the positive number of fish caught for a specific species (s) in a sampling trip in a specific river (r) followed a zero-truncated negative binomial distribution of mean = μ s,r with a log link: For models fitted to the full data, potential explanatory variables in the probability submodel (x i ) and in the positive catch sub-model (x j ) included sampling method, year, month, sampling duration, water temperature, river flow, distance to river mouth along river or longitude and latitude. For models fitted to the filtered data, potential explanatory variables included chl a and TP in addition to the above variables. Null models with species as the only hierarchical level were developed to evaluate the significance of incorporating hierarchical spatial effects.

Model fitting and comparison
A relatively new Bayesian approach, the integrated nested Laplace approximations (INLA) methodology, was used to fit the models. The INLA methodology can decreases the computational costs of the traditional Markov Chain Monte Carlo techniques (Besag et al. 1991, Rue and Held 2005, Rue et al. 2009a, Held et al. 2010). Its R interface allows us to develop sophisticated models more easily (Rue et al. 2009b). The default and recommended settings for priors were adopted (Held et al. 2010, Lindgren and Rue 2015, Simpson et al. 2017, Fuglstad et al. 2018). These priors were vague priors or approximations of "noninformative" priors that had little influence on the posterior distributions, and results were mostly derived from the data. For details on the priors, see Appendix S2.
We run a forwards stepwise variable selection process, in which a step involved testing the addition of each covariate separately. We started with the model only incorporating the variable year, because of our interest in the temporal trends of populations. Support for models with different explanatory variables was compared based on the deviance information criterion (DIC, Spiegelhalter et al. 2002) and Watanabe-Akaike information criterion (WAIC, Watanabe 2010).
The DIC is defined as: where D is the posterior mean of the deviance of the model, and p D is the effective number of parameters in the model (Spiegelhalter et al. 2002). The WAIC is defined as: where LPPD is the log posterior predictive density (Watanabe 2010). The DIC has been recognized for its tendency to under-penalize and select over-parameterized models (Plummer 2008). It can produce negative estimates of the effective number of parameters in a model. The WAIC is fully Bayesian and uses the entire posterior distribution, so it is recommended over the DIC criterion (Watanabe 2010, Gelman et al. 2014b. A smaller value of DIC or WAIC indicates a better explanatory quality of the model (Gelman et al. 2014b).
Model fit was measured using the conditional predictive ordinate (CPO, Geisser and Eddy 1979) that is an approximation of leave-one-out cross-validation: CPO i ¼ πðy obs i jy Ài Þ where y −i denotes the observations y with the ith component removed. It expresses the posterior probability of observing the value of y i when the model is fitted to all data except y i . Cross-validation avoids the problem of overfitting but remains tied to the data (Gelman et al. 2014b). The logarithmic score (LS) is defined by and a smaller value of the LS indicates a better predictive fit. If LS, DIC, and WAIC all show the same preference for a model, we have more evidence that the preference is correct.

Influence of global climate oscillations
Global climate oscillations-including the North Atlantic Oscillation (NAO), the Atlantic Multidecadal Oscillation (AMO), and the Gulf Stream index-were examined to see whether they were related to the population trends of shads and river herrings. These climate indices have been associated with changes in sea surface temperature, precipitation, wind fields, sea ice formation, phytoplankton, and zooplankton abundance, and consequently with changes in the abundance of marine fishes (Taylor 1995, O'Brien et al. 2000, Drinkwater et al. 2003. The NAO index was taken from the National Center for Atmospheric Research (NCAR 2019). The winter (i.e., December-March mean) NAO index was incorporated into our analysis because the signal:noise ratio of the NAO was strongest in winter (Hurrell et al. 2003). The annual AMO index based on the Kaplan Sea Surface Temperature dataset was from the Earth Science Research Laboratory (ESRL 2019). The annual Gulf Stream North Wall (GSNW) index, a metric of the latitudinal position of the Gulf Stream as it leaves the northeast coast of the United States as meanders, was obtained from personal communication with Dr. Arnold H. Taylor at Plymouth Marine Laboratory in the UK. The relationships between these climate indices and the standardized relative abundance indices of different Alosa species in different rivers were analyzed through their cross-correlation function (CCF), a statistical approach commonly used to investigate the possible time-lagged dependence between two variables (Probst et al. 2012). A P value ≤ 0.05 of the CCF was considered significant.

Model comparison
The DIC, WAIC, and LS values for models with different combinations of explanatory variables are presented in Tables 3, 4. When fitted to the full data, models with two-level hierarchies performed better than null models with species as the only level of hierarchy. The best-supported model (M7) retained year, longitude, month, water temperature, log-transformed river flow, log-transformed sampling duration, and sampling method as explanatory variables of the probability sub-model; year, latitude, month, Notes: L-Dur, log-transformed sampling duration; Temp, water temperature; L-Flow, log-transformed river flow; Lat, latitude; Lon, longitude; DIC, deviance information criterion; WAIC, Watanabe-Akaike information criterion; LS, logarithmic score. Models M follow two-level hierarchies with river as the upper level and species within river as the lower level, and models NM with species as the only level of hierarchy. DIC z , WAIC z, and LS z are for probability sub-model, DIC y , WAIC y, and LS y are for positive catch sub-model. Top performing models in each step are listed for brevity. Notes: Temp, water temperature; L-Flow, log-transformed river flow; L-Chl-a, log-transformed chl a concentration; L-TP, logtransformed total phosphorus concentration; Lat, latitude; Lon, longitude; DIC, deviance information criterion; WAIC, Watanabe-Akaike information criterion; LS, logarithmic score. Models F follow two-level hierarchies with river as the upper level and species within river as the lower level, and models NF with species as the only level of hierarchy. DIC z , WAIC z, and LS z are for probability sub-model, DIC y , WAIC y, and LS y are for positive catch sub-model. Top performing models in each step are listed for brevity. water temperature, longitude, log-transformed sampling duration, and log-transformed river flow as explanatory variables of the positive catch sub-model (Table 3).
When fitted to the filtered data, models with two-level hierarchies also performed better than null models with species as the only level of hierarchy. Model F7 was chosen as the bestsupported model for the filtered data, and retained year, longitude, water temperature, logtransformed river flow, month, log-transformed chl a concentration and sampling method as explanatory variables of the probability submodel, year, latitude, sampling method, longitude, water temperature, log-transformed TP concentration, and month as explanatory variables of the positive catch sub-model (Table 4).

Impacts of explanatory variables
Model results were plotted using the R package ggplot2 (Wickham 2016). The impacts of sampling method, sampling duration, water temperature, river flow, and month on catch probability and positive catch derived from model M7 fitted to the full dataset were combined to estimate their impacts on fish relative abundance (see Appendix S2 for details), which are shown in Fig. 3. The pound net was associated with greater relative abundance for Hickory Shad, Alewife, and Blueback Herring (Fig. 3a), indicating a greater catchability. The MDNR survey is a fishery-dependent survey, and commercial watermen determine the type, number, and location of the nets based on their fishing needs, so the pound and fyke nets tend to be set in locations with more fish. Among the two nets used, pound nets were used twice as often as fyke nets and were observed to catch more fish. Sampling duration showed a positive effect on relative abundance of American Shad in the James, Rappahannock, and Nanticoke rivers; Hickory Shad in the Nanticoke River; Alewife in the Rappahannock River; and Blueback Herring in the James, Rappahannock, and Nanticoke rivers (Fig. 3b). Water temperature between 10°C and 20°C was associated with greater relative abundance of American Shad and Hickory Shad in all rivers, except for the large uncertainties at low temperature in the Rappahannock and Nanticoke rivers; lower water temperature was associated with greater relative abundance of Alewife, while higher water temperature was associated with greater relative abundance of Blueback Herring (Fig. 3c). River flow displayed diverse effects on relative abundance for different species in different rivers; river flow effects in the York River had large uncertainties due to small sample size (Fig. 3d). The effect of month on relative abundance for different species showed river-specific heterogeneity: for American Shad, the month effects peaked in May in the Rappahannock and York rivers and peaked in April in the James and Nanticoke rivers; for Hickory Shad, month effects peaked in April in all rivers; for Alewife, month effects peaked in March in the Rappahannock, York, and Nanticoke rivers and peaked in April in the James River; for Blueback Herring, month effects peaked in May in the James, Rappahannock, and Nanticoke rivers, and peaked in April in the York River (Fig. 3e).
The impacts of chl a and TP concentrations on catch probability and positive catch derived from model F7 fitted to the filtered data were combined to obtain their impact on fish relative abundance, which are shown in Fig. 4. The impacts of these water quality factors displayed large uncertainties due to the small sample sizes and high variance of the metrics over time. There were only 11 sampling trips in the Rappahannock River and 46 sampling trips in the York River with water quality data ( Table 2). The relatively larger sample sizes in the James River (669 sampling trips) and Nanticoke River (364 sampling trips) helped to derive more accurate estimates of the effects of these variables and what they index. In the James River, relative abundance of American Shad was greater under lower chl a concentration (Fig. 4a); relative abundance of Alewife was greater under lower TP concentration (Fig. 4b). In the Nanticoke River, relative abundance of American Shad and Alewife was greater under lower chl a concentration (Fig. 4a); relative abundance of Alewife was greater under lower TP concentration (Fig. 4b).
Large uncertainties blurred the impacts in other cases.
The effects of latitude and longitude on catch probability and positive catch derived from model M7 fitted to the full data were computed together to obtain the changes of relative abundance indices over latitude and longitude (Fig. 5). In the James River, relative abundance of American Shad was greater in the northern inland area; relative abundance of Hickory Shad was greater in the Appomattox River; relative abundance of Alewife and Blueback Herring was greater near the coastal area. In the Rappahannock River, relative abundance of Alewife was greater in the inland area; relative abundance of Blueback Herring was greater near the coastal area. In the York River, relative abundance of Alewife was greater in the coastal area. Two main sampling sites in the James River (i.e., Belle Isle and Bosher's Dam) that were located in the fall zone and observed to have zero probability of capturing Hickory Shad and Alewife and slight probability of capturing Blueback Herring showed low spatial effects on the three species.

Standardized relative abundance indices and correlations with long-term climate oscillations
The estimates of intercepts and year effects from two sub-models of model M7 fitted to the full dataset gave the standardized relative abundance indices for the four species in all rivers ( Fig. 6; see Appendix S2 for details). Each species displayed different population trends in the respective rivers. The standardized indices and geometric mean survey indices showed similar patterns, but with some differences. Hickory Shad was usually caught in small numbers by fyke and pound nets, since they tended to actively avoid these gears as observed by sampling crews. The data captured by the MDNR survey for Hickory Shad might not represent v www.esajournals.org their abundance effectively, potentially reflected by the low relative abundance indices of Hickory Shad in the Nanticoke River. We correlated the standardized relative abundance indices of each species in each river with long-term climate indicators. The annual trends of the respective Alosa species were significantly cross-correlated with the climatic variables only in a few cases: In the Rappahannock River, the GSNW index of the previous year displayed a significant positive effect on the standardized relative abundances of Alewife, the GSNW index of the current year displayed a significant positive effect on the standardized relative abundances of Blueback Herring, the AMO index of the previous year displayed a significant negative effect on the Fig. 6. Z scores of standardized relative abundance indices and geometric mean survey indices of the four species in the four rivers derived from model M7 fitted to full data. Lines represent posterior mean values; ribbons represent 95% credible intervals. standardized relative abundances of Blueback Herring; and in the Nanticoke River, the AMO index of the previous year displayed a significant negative effect on the standardized relative abundances of Alewife and Blueback Herring (P < 0.05).

Impacts of explanatory variables
Our analyses demonstrate the river-specific effects of environmental factors on relative abundance indices of Alosa species in selected tributaries of the Chesapeake Bay. Previous studies have documented that the optimal spawning temperature for Alewife is 4-19°C and for Blueback Herring is 14-22°C (O'Connell and Angermeier 1999, Ogburn et al. 2017b). Our results indicate that relative abundance of Alewife adult was greater under lower water temperature; in comparison, relative abundance of Blueback Herring adult was greater under relative higher water temperature, which is consistent with these previous findings. It also was found that relative abundance of Alewife peaks in earlier months, when water temperature is lower, while that of Blueback Herring peaks in later months, which is consistent with previous observations (Plough et al. 2018).
Variations in river flow have shown effects on the energy expenditure of migrating Alosa adults and spawning habitat availability (Haro et al. 2004, Walsh et al. 2005. High river flow could decrease feeding efficiency and migration distance of anadromous clupeids (Limburg 1996, Haro et al. 2004. Low river flow could also influence movements of Alosa populations, and it was observed that adult river herring moved downstream out of spawning areas on the Choptank River during low and declining flows (Ogburn et al. 2017b). Large river systems rarely become dewatered, and thus, fish might be less susceptible to stranding than in small systems (Kosa and Mather 2001). Previous studies also hypothesized the effect of river flow to be site-and speciesspecific (Tommasi et al. 2015). Alewife tend to use low-flow habitats, while Blueback Herring could use both high-and low-flow habitats (Walsh et al. 2005), which is consistent with the effects of river flow on relative abundance indices of these species, especially in the Rappahannock River.
More primary production tends to benefit fish larvae and juveniles (Grimes and Finucane 1991). However, extremely high chl a concentration might cause hypoxia and eutrophication and decrease trophic transfer efficiency, which would negatively influence shad and river herring populations (Havens et al. 2000, Cai et al. 2011. Although with large uncertainties due to small sample size, our results indicate that Alosa species, American Shad in the James and Nanticoke rivers and Alewife in the Nanticoke River in particular, are negatively influenced by extreme high chl a and TP concentrations. More information on water quality is necessary to obtain more accurate results.

Management implications
The sampling surveys conducted by VDGIF and MDNR provide data on the spatiotemporal dynamics of Alosa populations and allow us to assess environmental and anthropogenic effects on their populations in different river systems. A key benefit of this study is that each species in each river system is found with its threats characterized, and a corresponding habitat restoration is able to be developed for each species within a river. Stocks in different rivers display different trends: Blueback Herring in the James River, American Shad, Hickory Shad, and Blueback Herring in the Rappahannock River have been increasing in recent years; American Shad and Alewife in the James River, American Shad in the Nanticoke River are still at lows; time-series trends of stocks in the York River are difficult to distinguish due to missing observations in many years. For stocks not showing a recovery trend, extra habitat restorations need to be considered.
The impacts of environmental and anthropogenic factors on river-specific stocks provide useful information to develop habitat restoration strategies. American Shad in the James River are negatively influenced by high water temperature, low river flow, and high chl a concentration. It was found that within Virginia land use pressure is highest along the James River at Richmond, followed by the confluence of the James and Chickahominy rivers and the peninsula separating the James River from the York River (Hilton et al. 2014), where with higher relative abundance of American Shad (Fig. 5). Land use around rivers is likely associated with water pollution including nitrogen and phosphorous yields, sediment load, and thermal effluent (Hilton et al. 2014), which would negatively affect American Shad. Land use regulations associated with water quality are expected to help to recover this shad stock (Hilton et al. 2014). Alewife in the James River are negatively influenced by high water temperature and high TP concentration; American Shad in the Nanticoke River are negatively influenced by high water temperature and high chl a concentration. Similarly, land use regulations associated with water quality would benefit the two stocks.
The estimated spatial effects indicate habitat preference of each species within rivers. Spatial locations with greater effects on relative abundance of Alosa species could be considered by managers as they consider how to focus on habitat restoration efforts. For example, restorations in inland areas of the James and Nanticoke rivers could help to recover American Shad stocks; restorations in coastal areas of the James River could help to recover Alewife stock. State and federal agencies deal on a regular basis with time-of-year restrictions for in-stream work to protect anadromous fish spawning runs. The month effects estimated from our model could be factored into such decision making. To be more specific, restrictions for in-stream work in April would benefit American Shad stocks in the James and Nanticoke rivers; restriction in March and April would benefit Alewife stock in the James River.
Inter-annual demographic trends of some Alosa stocks are related to large-scale climate oscillation-AMO and Gulf Stream meanders. These types of relationships are observed for river herrings in the Rappahannock and Nanticoke rivers. The discontinuity of samplings in the York River makes it difficult to observe any potential relationships. A positive phase of the AMO index may bring about warming water, increased drought severity, loss of flood plain spawning habitat, and negatively affect Alosa populations (NCDMF 2014). A larger GSNW index indicates the Gulf Stream meanders follow a more northerly track and would attract fish to the river in a higher latitude. Under global warming scenarios, Alosa populations may move northward, or the spawning timing may shift to an earlier month. The situation will be more severe for Alewife, because they tend to stay in habitats with lower temperature. These shifts may lead to other factors affecting their abundance and cause changes in managements, such as changes in Time of Year restrictions for instream work. The specific effects of climate change should be noticed and further investigated.
The ASMFC completed a coastwide benchmark stock assessment in 2012, concluding the overall coastwide population of river herring stocks on the U.S. Atlantic coast is depleted to near historic lows (ASMFC 2012). As indicated in the stock assessment, river herring should be assessed and managed by individual river systems; however, for the majority of rivers, data were not available to conduct a model-based stock assessment (ASMFC 2012). Instead, temporal trend analysis was used to identify recovery patterns in the available data (ASMFC 2012). River herring were assessed at a river scale for three stocks only, including the Monument River in Massachusetts, the Nanticoke River in Maryland, and the Chowan River in North Carolina (ASMFC 2012). With more data being collected and powerful statistical approach being developed (such as the hierarchical models) in the future, it will be possible to conduct modelbased stock assessment in more river systems, and incorporate river-specific heterogeneities in population dynamics to facilitate development of coastal stock assessment.

Recommendations for data collection
The fisheries-dependent and independent monitoring programs conducted by states and jurisdictions provide a database for river-specific trend analyses. However, there are some limitations in the datasets used in this study. First, data in the York and Nanticoke rivers have gaps in some years, which hinder us from getting yearly continuous trends of stocks in these two rivers. More consistent samples help to derive timeseries trends. Second, these surveys take weekly sampling strategies. However, anadromous spawning runs are generally highly episodic (Ogburn et al. 2017b), so the current sampling strategies often yield high variabilities. To capture more accurate trends, higher-frequency sampling strategies are needed. Some techniques such as imaging sonar can record run counts over shorter periods and provide more informative data (Ogburn et al. 2017b). Third, the current spatial coverage is insufficient. The numbers of sampling stations in the Rappahannock and York rivers are limited; the sampling stations in the Nanticoke River are close in space. Additional separate stations are needed to cover a larger spatial field and improve the picture of population trends. Fourth, our results indicate that Alosa species are likely negatively influenced by high chl a and TP concentrations, but the results have large uncertainties due to small sample size on the water quality data. It is suggested that water quality metrics be taken in association with fish samples. In addition, Alewife have been observed to move on and off of spawning grounds during the spawning season (McCartin et al. 2019). It is likely that fish move back and forth between tidal and non-tidal areas during the spawning season, particularly as flow increases or decreases. However, tidal status was not recorded in most sampling trips, and based on the trips with tidal status, its impact was not significant. More records in tidal status in the future sampling may help to diagnose its impact. Fifth, there are other datasets including long-term data in the Potomac River and shorter time-series data in several other rivers (Northeast, Choptank, Patapsco, etc.). There are also juvenile abundance indices data available in the Chesapeake Bay region. We consider to incorporate these data into analyses in the future. Future management could be better informed by models derived from more consistently collected and more informative data.
Although it was expected that the relationship between the number of fish caught and sampling duration is non-negative, our analyses show that the number of fish caught in a sampling trip is lower under longer sampling duration for Alewife in the James River and Hickory Shad in the Rappahannock River. In the James River, the low catch number under greater sampling duration in the data occurred in May 1994 and April 2014 in inland areas. In the Rappahannock River, the low catch number under greater sampling duration in the data occurred in March 2000 near coastal areas. Sampling crews sampled established stations consistently based on a predetermined time duration. Variation in the sampling duration could occur if anglers were in the area or the current moved the boat through the station at faster rates, either resulting in shorter samples time wise. Given the variation in current velocity, sampling duration might not be a perfect index for sampling effort, and a more appropriate index such as a formula combining pedal time (number of minutes electric current is applied during electrofishing surveys) and area covered is needed for future surveys.

Model improvements
The present study developed a computationally efficient and statistically powerful approach to analyze river-specific heterogeneity of the spatiotemporal dynamics of Alosa and identify riverspecific impacts of environmental and anthropogenic factors on these stocks in some Chesapeake Bay rivers. Spatial studies of species within river systems are complex and difficult to implement, because spatial correlation can exist between rivers and within rivers. Under the current model structure, the within river effects of latitude and longitude could be conflated with the between river effects of location. A one-dimensional process linked at the points where streams join, and a twodimensional process linked to more general spatial features can be combined to address spatial effects at different scales in the future. More sampling locations across the region help to extend the model to space.
In summary, Alosa stocks display different trends in selected river systems of the Chesapeake Bay. Environmental and anthropogenic factors, including water temperature, river flow, chl a concentration, and TP concentration, play riverspecific impacts on these populations. The riverspecific impacts of these factors should be considered when we detect population status in specific rivers and develop potential recovery strategies. River-specific stock trends could help to facilitate development of coastal population models.