Effect of straying, reproductive strategies, and ocean distribution on the structure of American shad populations

Abstract

The use of species distribution models has proliferated, providing insights for sustainable management of migratory species in a globally changing environment. However, many of these models are based on statistical relationships developed from historical conditions that may not perform well under changing or even analogous conditions caused by climate change. In this paper, we used a mechanistic species distribution model called GR3D (Global Repositioning Dynamics for Diadromous Fish Distribution) to examine the integrated dynamics of American shad (Alosa sapidissima) populations across their native range along the Eastern U.S. coast, where the species demonstrates latitudinal variations in life histories and reproductive strategies. The initial design of the model was adapted to incorporate region-specific parameterization to fit the species ecology. Then, a sensitivity analysis was performed to test the influences of uncertain processes regarding American shad distribution at sea, straying and reproduction on key characteristics of the species distribution. The sensitivity analysis showed the influence of the Allee effect (i.e., “depensatory” process) and the homing rate (i.e., fidelity to the breeding sites) on the probability of presence and abundances among catchments and metapopulations estimated by the model. Contrary to the homing rate, the distance of straying did not change the estimated number of metapopulations or abundances. Homing strength, however, was quite influential. The integration of complex migration patterns during the marine phase (i.e., wintering and summering offshore areas) provided more likely estimates of the species' overall distribution. Overall, our study illustrated the utility of incorporating factors governing the large-scale distribution of migratory species to improve local management.

INTRODUCTION

Climate change and other anthropogenically driven stressors necessitate understanding current species ranges and prevailing future distributional shifts and the resulting interactions with human activities. Over the past decades, species distribution models have become an important tool for identifying risks to wildlife populations and to develop mitigation strategies and management approaches (Dawson et al., 2011; Hoegh-Guldberg et al., 2008). A suite of statistical approaches combining observed geographic distribution (i.e., occurrence/abundance) and climate data have been developed (e.g., Elith et al., 2011; Phillips et al., 2006) and used to predict future range and suitable habitats for many taxa including plants (Cole et al., 2011), amphibians and reptiles (Araújo et al., 2006), and freshwater and diadromous fish (Buisson et al., 2008; Lassalle et al., 2009). With advancement of mechanistic models, simulations can now be used to integrate species' functional traits and dispersal (Kearney & Porter, 2009; Pacifici et al., 2015). These efforts provide a better understanding of mechanisms driving biodiversity and ecosystems changes and particularly species distributional responses over large spatiotemporal scales.

Mechanistic species distribution models (M-SDMs) integrate population dynamics across a species range with spatially explicit environmental data triggering responses to environmental conditions (Kearney & Porter, 2009). Primary applications of M-SDMs have addressed species vulnerability or extinction rates (see Pacifici et al., 2015), though this approach also may aid in reducing ecological and biological uncertainties associated with the target species. Building and parameterizing a species-specific M-SDM requires a substantial amount of data and knowledge so that these initiatives often lead to the identification of critical knowledge gaps and opportunities for future assessments. Nonetheless, when conceptualized and fully parameterized, M-SDMs may provide a solid framework to test various ecological hypotheses, thereby narrowing the range of influential ecological processes (Gaucherel et al., 2014; Guisan & Zimmermann, 2000; Winsberg, 2003).

Diadromous species, migrating between freshwater and marine environments (McDowall, 1988) are particularly sensitive to environmental changes because they experience increasing human pressures acting on both aquatic environments. Many diadromous species have experienced severe population reductions and extirpations over their distribution range to the point that most populations are now considered of high conservation concern (Limburg & Waldman, 2009; Wilson & Veneranta, 2019). In this context, an M-SDM named GR3D (Global Repositioning Dynamics for Diadromous fish Distribution; Rougier et al., 2014) was built to assess the spatial population dynamics of diadromous species in a globally changing environment with estimations of population size through a species' distribution range. The model was first applied to allis shad (Alosa alosa Linnaeus 1758) which experienced a continuous decline across Europe over the 20th century until now. Through this first application, we enhanced our mechanistic understanding of drivers behind the species' distribution edges. Insights from this model application could be beneficial to other shad species with similar histories, such as the American shad (A. sapidissima Wilson 1811).

The American shad is an anadromous fish species native along the Eastern Atlantic Coast of North America where it is widely distributed from the St. Johns River (Florida) to the St. Lawrence River (Canada) (Limburg et al., 2003). In response to a combination of overfishing, inadequate fish passages at dams, water pollution, and climate change, the species has collapsed throughout its range (Limburg & Waldman, 2009). American shad spend most of their life at sea, feeding and growing before returning to rivers to spawn. The life-history traits and population dynamics of the American shad vary with latitude (Glebe & Leggett, 1981). Latitudinal variations include changes in spawning periodicity, growth rate, reproductive strategies, fecundity (Limburg et al., 2003), and genetic structuring (Hasselman et al., 2013). At the southern edge of the species range, in Florida, the spawning migration begins in January, occurring successively later (May or June) in northern populations (Limburg et al., 2003). Likewise, fish in the north reach older maximum ages and larger maximum sizes than populations in the south (Gilligan-Lunda et al., 2021). Populations spawning in the southern extent of their range are semelparous (i.e., reproduce once and then die) while those spawning north are iteroparous (i.e., capable of reproducing more than once in a life time) (Leggett & Carscadden, 1978). As a result, fish in southern populations have higher annual fecundity, but similar lifetime fitness is achieved by iteroparous fish in the north through repeat spawning (Leggett & Carscadden, 1978; Olney & McBride, 2003). By examining 33 populations of American shad, Hasselman et al. (2013) revealed that these patterns resulted in distinct genetic clusters of populations: a cluster of semelparous populations in the south and two clusters of iteroparous populations in the north. These same three clusters are used for managing the species in the absence of river-specific information (ASMFC, 2020). However, with as many as 165 potential populations in Canada and the United States (Zydlewski et al., 2021), some aspects of American shad ecology remain open to debate.

The marine phase of the American shad life cycle remains enigmatic. After a short freshwater residency, most juveniles move into marine waters to form schools along the coast (Greene et al., 2009; Neves & Depres, 1979), while some individuals remain in the lower estuarine reaches of their natal rivers to spend their first summer (Hoffman et al., 2008). Juvenile emigration toward the sea is commonly associated with fall (Limburg, 2001; Zydlewski & McCormick, 1997; Zydlewski et al., 2003) although there is a high degree of inter- and intrapopulation variability, with juveniles from northern rivers emigrating seaward earlier in the year (Greene et al., 2009; Walther & Thorrold, 2010). Previous studies suggested that American shad then travel to offshore areas to find suitable temperature conditions for feeding and growth, displaying some seasonality in their life at sea (Dadswell et al., 1987). Most American shad are found off the Florida coast, in the Mid-Atlantic Bight (MAB) and as far as the Scotian Shelf during winter. Then, fish move northward to offshore summering areas in the Bay of Fundy, St. Lawrence estuary, or coastal waters of Newfoundland and Labrador (Dadswell et al., 1987). Populations from many river systems are hypothesized to travel and mix in these regions, although little direct evidence exists (Dadswell et al., 1987; Gahagan & Bailey, 2020). The current lack of basic knowledge relating to the marine distribution and movements at sea makes it an important direction of investigation.

After 3–5 years of marine residency, mature American shad may return to spawn in the river where they originated (homing) or they may disperse to a new river to attempt reproduction at non natal sites (i.e., straying; Quinn, 1993). American shad enter their natal rivers and tributaries to spawn with high fidelity (Melvin et al., 1986), though some population mixing through straying is necessary to have resulted in the life-history variability and population genetic structuring observed in American shad, particularly in systems that were colonized after glaciation (Hasselman et al., 2013). Straying allows individuals to move to new suitable areas and potentially serves as a buffer against changes in environmental conditions or habitat quality (Hendry et al., 2004). A straying rate of 3% is commonly accepted for American shad, but this is based on a single study of one river in Canada (Melvin et al., 1986).

The potential benefit of straying is also linked to how a species is distributed in space and time. If adults stray to new or unoccupied rivers, there may be no net benefit from this straying because individuals are less likely to find a suitable breeding habitat and mate (Hendry et al., 2004). Reproduction benefits from the numbers or densities of conspecifics so that average individual reproductive success is expected to increase at high densities (“Allee effect”; Stephens et al., 1999).

The Allee effect, a concept that broadly links species fitness to population size, describes a nonlinear relationship between population density and the per capita growth rate of a population (Stephens et al., 1999). At low population density, the Allee effect leads to reduced reproduction or survival and can result in critical population size below which the population cannot persist (Courchamp et al., 2008). This depensatory process has not been mechanistically explored in American shad populations. It was, however, recognized as a possible explanation for population collapses and extirpations of other fish species, including allis shad (Rougier et al., 2012) and Atlantic cod (Kuparinen et al., 2014; Kuparinen & Hutchings, 2014).

In this study, our goal was to provide a deeper understanding of the main biological factors governing the large-scale distribution of American shad to tackle the knowledge gaps related to American shad ecology and population dynamics, especially at sea. To achieve this, we adapted the existing GR3D model to fit American shad ecology in its native range. Then, we ran a sensitivity analysis to test likely combinations of uncertain biological processes on model outputs related to American shad distribution and metapopulation structuring. Main knowledge gaps and data limitations to test in the sensitivity analysis were informed by consulting species experts. In particular, we tested the most likely scenarios of ocean migrations while assessing the influence of assumed reproductive strategies and straying/homing on American shad population dynamics across its native range.

MATERIALS AND METHODS Presentation and adaptation of the M-SDM design GR3D main features

GR3D, as described by Rougier et al. (2014, 2015), is an individual-based model developed for allis shad to evaluate the species range-shift into a climate change context. The GR3D model integrates three types of entities: one for fish, and two for the environment in which fish evolved that is, the “sea” and “continental” compartments. The continental compartment is composed by a set of “river basins” while the sea compartment is split into “sea basins” located in front of each river basin (Rougier et al., 2014).

Fish progress with a seasonal time step and move from one to the other compartment according to the stage at which process occurs meaning that the model is spatially and temporally explicit (Rougier et al., 2014, 2015) (Figure 1).

FIGURE 1. Conceptual diagram of the timing and scheduling of fish movements along the four types of entities included in the model physical environment from birth (Y0) to reproduction. Migrations between each “box” are indicated by arrows: dotted arrows symbolize seasonal migrations between offshore basins and full arrows migrations of adults and juveniles from rivers to the shores. Adults spawn simultaneously in each river basin in spring and produce a stock of juveniles that start their migration to sea and enter the inshore basins.

GR3D covers the entire life cycle of an anadromous species by integrating six key fundamental biological processes along with their thermal requirements: (1) reproduction, (2) downstream migration, (3) growth, (4) survival, (5) maturation, and (6) upstream migration (Rougier et al., 2014). Though variability of life-history traits may have a genetic dimension, GR3D only models variability in these processes as plastic responses to changing environmental conditions (Rougier et al., 2014). Plastic responses to temperature influenced growth of individuals at sea, survival of spawners (both pre- and post-spawning), and early-life stage survival in rivers. Growth, survival, spawning dispersal in river, and reproduction are modeled as stochastic processes.

When applying GR3D to the new case study, only minor modifications to the overall model framework were required to fit American shad ecology and life history across its native range. Major updates relative to this case application included a temperature-dependent post-spawning survival and the addition of new spatial entities corresponding to offshore basins where shad seasonally occurs. The physical environment is now split into four entities instead of two: the “river basins,” the “inshore basins” (formerly the “sea basins”), “the wintering offshore basins” where shad spend winter and spring, and the “summering offshore basin” to which shad migrate in spring to spend summer and fall. Following this new marine environment, a seventh biological process named (7) migrations between summering and wintering habitats at sea was added. A conceptual model depicting the timing and scheduling of the model applied to the new species is provided in Figure 1. Consecutive changes in population dynamics are described in detail in the next section while detailed descriptions of the initial model design and modeling choices were provided in Rougier et al. (2014, 2015).

Overview of the population dynamics within the model Reproduction

Reproduction was modeled to occur simultaneously in all river basins each spring. Variation in American shad spawning phenology occurs along the Eastern U.S. coast (Limburg et al., 2003) but was not included since it would have required changing the model time steps, from seasonal to monthly. The number of recruits produced by the spawning stock in a river basin was modeled as a density-dependent process using a modified Beverton and Holt stock-recruit curve (BH-SR; Beverton & Holt, 1957; Appendix S1: Equations S1.a–f). The relationship between recruitment mortality and temperature followed a dome-shaped curve (Rosso et al., 1995) with recruit survival decreasing outside the range of optimal temperatures defined for the species (Appendix S1: Equation S1.f and Table S1).

Growth, maturation, and spawners' dispersal into the rivers

Growth of individuals was modeled using a von Bertalanffy growth function (hereafter VBGF; Von Bertalanffy, 1938) with stochasticity of seasonal growth increments described by a lognormal distribution (Appendix S1: Equation S2.a and Table S1). Assuming that growth depends on water temperature (Brown et al., 2004; Gilligan-Lunda et al., 2021), we provided an extension of the common VBGF function to account for seasonal variations in growth as described by Bal et al. (2011). The seasonal growth increment was modulated by water temperature (T) using a modified Lehman function (Svirezhev et al., 1984; Appendix S1: Equation S2.a).

Fish remained in the ocean until they reached the size of maturity and then entered a river basin and spawned. For “homers” (i.e., fish returning to natal rivers to spawn), the odds of migrating into the natal river was derived from a constant probability (Appendix S1: Table S1). For “strayers” (i.e., individuals spawning in rivers other than natal origin), we assumed that the distance between the natal river and the new destination river was the only driver. This is a simplifying assumption though attractiveness could result from multiple factors such as discharge or presence of conspecifics (e.g., Vrieze & Sorensen, 2001; Appendix S1: Equation S3.a). The probability of straying was described by a logit kernel function on the standardized distance between the natal and the new rivers (Appendix S1: Equation S3.b and Table S1).

Survival

At each time step, individual survival was estimated using its location and employing annual mortality at sea and in river. The natural mortality coefficient in river was assumed to be a function of water temperature. In river, mortality occurred for adults before or after they spawned depending on prior spawning history (i.e., iteroparity). In comparison with the original version of the model (Rougier et al., 2014), the probability of repeat spawning for American shad was imposed through post-spawn survival. For juvenile recruits, the probabilities of survival were computed both from hatch to 14 days post-hatching (dph) and from 14 dph to downstream migration. Survival of juvenile recruits was a function of temperature and was included in the stock-recruitment relationship through a dome-shape curve (Rosso et al., 1995) (Appendix S1: Table S1). Once they migrated to sea, juvenile survival to age 1 adult was calculated based on the annual mortality coefficients at sea as previously described.

Model parameterization for the American shad Distributional and environmental data Species distribution

Historical records of American shad in rivers used to run simulations were gathered from Walburg and Nichols (1967). Presence/absence data were reported in 98 river basins ranging from the St. Johns River (30.40° N, Florida) to the George River (58.61° N, Northern Canada; Figure 2). Due to a lack of information about the distribution of shads in mixing zones at sea, positions (latitudes and longitudes) of the wintering and summering habitats (“offshore basins”) were based on Dadswell et al. (1987) and approximated to be areas of 10,000 km² presumed to be at preferred depths for the species (between 50 m for shallowest areas [e.g., Bay of Fundy] to 200 m for deeper areas [e.g., MAB for Mid-Atlantic Bight]; Neves & Depres, 1979; Figure 2b). Bathymetric contours used to position the wintering and summering areas along the Atlantic Coast, from Florida extending to the Gulf of Maine, were derived from the United States Geological Survey (USGS; https://www.usgs.gov/). For the Bay of Fundy, bathymetric data were provided by the Commonwealth of Massachusetts Bureau of Geographic Information (MassGIS; https://www.mass.gov/orgs/massgis-bureau-of-geographic-information).

FIGURE 2. Physical environments of the GR3D (Global Repositioning Dynamics for Diadromous Fish Distribution) model tested in the sensitivity analysis showing the two hypotheses tested into the experimental design. (a) The “River-Inshore” hypothesis holds that shad do not mix into the ocean and stay in the vicinity of their natal river to grow. The physical inshore environment was constructed as a set of 98 “inshore basins” associated with each of the 98 rivers. (b) The “River-Inshore-Offshore” hypothesis assumed fish formed mixed schools in the ocean overwintering in suitable habitats (named wintering and summering “offshore basins”). Summering and wintering “offshore basins” were described from south to north and broadly corresponded to summer and winter aggregations defined by Dadswell et al. (1987). Therefore, the physical environment included “offshore basins” where juveniles and subadults grew and “inshore basins” from which spawners strayed.

Temperature data

Observed temperature time series were used to link life-history traits of the American shad to water temperatures. In-river temperature was estimated by averaging monthly atmospheric data and sea surface temperature (SST) at the outlet (hereafter called “SST_outlet”). Observed atmospheric data were extracted from the Climate Research Unit time series dataset (CRU TS v.4.03; crudata.uea.ac.uk; Harris et al., 2020) which consists of a monthly grid of land-based climate observations gathered from 1901 to 2018 with a resolution of 0.5° (~62 km). SST for the “inshore” and “offshore” basins (“SST_outlet” and “SST_offshore,” respectively) came from the CHORE_AS (CMCC Historical Ocean Reanalysis) dataset. This is a reanalysis performed by the CMCC (Centro Euro-Mediterraneo sui Cambiamenti Climatici; http://c-glors.cmcc.it/; Yang et al., 2017) reconstructing ocean variables from 1900 to 2010 at a resolution of 0.5°. Monthly time series of air temperature and SST were averaged over 3-month periods to provide seasonal temperatures.

Biological input parameters

Model parameters were calibrated using available field or laboratory data from the literature or were derived from expert knowledge (Appendix S2).

Each fish in the model was assigned to a river corresponding to its natal origin and sorted into one of the three regions based on the reproductive strategies and locations (sensu Hasselman, 2010). These regions were termed “semelparous” (from Florida to Cape Fear), “southern iteroparous” (north of Cape Fear to the Hudson River, NY), and “northern iteroparous” (north of the Hudson River) following the ASMFC (2020) classification (Appendix S3).

Sensitivity analysis

The influence of two key biological processes for which uncertainty remains, straying and reproduction, was examined on the overall American shad population distribution by changing: (1) the presence of an Allee effect, (2) the post-spawn survival, (3) the distance at which 50% of strayers find a destination catchment, and (4) the spawning fidelity (Figure 3). The range of input parameters was guided by knowledge gaps identified by the questionnaire addressed to experts and workshop discussions.

FIGURE 3. Conceptual diagram of the experimental design used to investigate the effect of reproduction and straying on the species distribution and metapopulation structure. Gray circles represent the biological processes targeted by the sensitivity analysis with the boxes corresponding to the uncertain processes and their modalities. The three model outputs examined in the sensitivity analysis appear in the box on the right. Cst, constant Allee effect; GR3D, Global Repositioning Dynamics for Diadromous Fish Distribution; H, high; I, intermediate; L, low; RI, River-Inshore; RIO, River-Inshore-Offshore; SP, size proportional; TD, temperature-dependent.

Fifty-year simulations (from 1900 to 1950) were performed to test the influence of straying and reproduction on population distribution and structure. An ANOVA-based global sensitivity analysis was conducted to evaluate model sensitivity to changes in input parameters (Figure 3). More specifically, ANOVA assessed the average contribution (main effect) of each input parameter to the overall outcome variance by also considering the interactions (Ginot et al., 2006; Saltelli et al., 2000).

Ecological assumptions behind the sensitivity analysis Simulation of shad migration routes and growth at sea

The sensitivity analysis consisted of two distinct physical environments at sea. Considering the currently limited knowledge regarding the marine stage of shad, we included two scenarios representing hypotheses about shad migration routes at sea. The first hypothesis, the “River-Inshore” hypothesis (“RI”), assumed that shad remained in the vicinity of their natal river to grow, as previously suggested for allis shad in several European rivers (Nachón et al., 2019; Taverny, 1991; Taverny & Elie, 2001). Accordingly, we defined the “inshore basins” located along the coast at the main entrance of each of the 98 rivers (Figure 2a). Under this assumption, juveniles emigrated from their natal river in summer to enter the “inshore basin” and remained there until returning to rivers at maturity (Figure 2a).

An alternative hypothesis, the “River-Inshore-Offshore” hypothesis (“RIO”), hypothesized that shad followed seasonal trophic migration routes at sea. This would be congruent with tagging studies showing that shad from different origins mixed at sea in spatially aggregated areas (Dadswell et al., 1987). Since there is a lack of data suggesting otherwise, we considered the juveniles to remain in the vicinity of their natal river (i.e., in the “inshore basins”) to spend summer and fall and then move to overwintering sites (called “wintering offshore basins”) to spend winter. The “offshore basins” were spatially structured based on the three regions where shad seasonally occur. “Wintering offshore basins” included (1) a region south of Florida, (2) the MAB, and (3) the Scotian shelf. During spring and summer, shad moved to areas that were called “summering offshore basins” including the (4) Inner Bay of Fundy, (5) St. Lawrence estuary, and (6) Newfoundland and Labrador (Figure 2b). Shad therefore experienced different thermal conditions depending on migration routes and locations at which fish aggregate.

The choice of the wintering offshore basin when leaving a river basin was driven by the closest distance separating the river basin from the wintering areas. Fish then migrated to one of the three “summering offshore basins” according to migration routes defined by Dadswell et al. (1987). Fish native to the southern, middle, and to a lesser degree the northern extent of their range and overwintering off Florida, in the MAB or in the Scotian Shelf mixed into the Bay of Fundy, while northernmost populations which spend winter in the Scotian Shelf and moved as far north as the St. Lawrence River estuary or coastal waters of Newfoundland and Labrador (Figure 2b; Appendix S4).

Growth parameters were calibrated to fit the observed growth curves for the three regions regardless of fish sex (ASMFC, 2020; Gilligan-Lunda et al., 2021). Observed growth curves were based on a set of 11 river systems for which length-at-age data were available for each of the three regions (Appendix S3). For our analysis, 10 rivers were kept to fit with the rivers included into the model environment (i.e., exclusion of the North River). Accounting for variability in growth and age at maturity between sexes, length at maturity was calibrated using maturity ogives defined in ASMFC (2020) (Appendix S1: Equations S4–S7b).

Spawner fidelity

For spawner fidelity, we considered either high or low rates, meaning that a large or lesser proportion of mature fish returned to their natal river to spawn. The homing rate was set to a high value of 0.95 (“H” in Figure 3) according to the straying rate (Stray = 1 − p_hom) commonly provided in the literature for the American shad (Melvin et al., 1986; Walther & Thorrold, 2010) while the low value of 0.75 (“L” in Figure 3) was derived from the European species with a lower fidelity to their natal systems (Rougier et al., 2014).

The relative distance to which half of the strayers were likely to stray was set to high, intermediate, and low values (“H,” “I,” and “L,” respectively, in Figure 3), of 88.3, 53.6, and 19 km, respectively, based on the answers provided by experts (group averages) (Appendix S2).

Reproduction dynamics

Since no evidence for or against an Allee effect in American shad populations currently exists, three scenarios were tested. Options included absence of an Allee effect, a constant Allee effect, or an Allee effect proportional to the basin size. A constant Allee effect assumed that the minimum number of spawners to successfully spawn was the same among all the 98 rivers (i.e., a constant Allee effect, simplified as “cst” in Figure 3). The second assumed the minimum number of spawners to be proportional to the river basin size (i.e., “size proportional”; simplified as “SP” in Figure 3). This assumption reflected the likely scenario that encounter rates decrease in large rivers when spawner numbers are low.

Post-spawn survival

Latitudinally correlated variations in the percent of repeat spawners (Leggett & Carscadden, 1978) are hypothesized to be the consequence of environmental conditions during early-life stages of American shad. It is hypothesized that southern rivers with warmer temperatures and more stable environments foster semelparity, while unpredictable changes in temperature and flow conditions in northern rivers favor individuals that may repeat spawning (Leggett & Carscadden, 1978). Therefore, we assumed that iteroparity was a consequence of survival of spawners after reproduction (Sp_sp) through a linkage with the water temperature (i.e., survival was temperature-dependent, simplified as “TD” in Figure 3). The probability of survival was then expressed as a logistic function as follows:[Image Omitted. See PDF]where log(19)/(T_minObs − T) corresponded to the slope of the logistic function; $T_{\text{minObs}}$ was the temperature above which there was no repeat spawning; T_j was the temperature in the river in summer and T_r the averaged summer temperature for all the river basins considered; and β the theoretical maximal survival. $T_{\text{minObs}}$ was set to 25°C according to the mean temperature observed in the southernmost St. Johns River for which absence of repeat spawners was reported (Walburg & Nichols, 1967). Post-spawn survival was thereby set to low values for highest temperatures and high values for coldest temperatures. Low post-spawn survival decreased the chance of fish to return seaward successfully and reproduce the following year.

To assess the influence of post-spawn survival on the distribution and structure of the populations, we tested an alternative hypothesis that an average of 10% of individuals survived spawning across the distribution range, as has been proposed in the allis shad GR3D application (Rougier et al., 2014).

Model outputs in the sensitivity analysis

Simulations were examined through (1) the log-likelihood statistic (hereafter “LogL”) (2) spawner abundances and (3) the number of emerging metapopulations (Figure 3). The log-likelihood is an indicator of how well the model, with its specific parameter set, fits the observations. Here, the log-likelihood was defined as the probability of presence in a river basin j and was calculated as the proportion of spawning leading to a mean recruitment of more than 50 recruits over the last 10 years of simulations, in accordance with Rougier et al. (2015). The resulting log-likelihood values indicated which hypothesis better fits the observed historical distribution of American shad. The higher values, the better the model fits the observed distribution.

The spawner abundance estimates were a function of the number of spawners entering the river in spring across all rivers and provided an estimate of the size of the overall spawning runs. We posited that each river contained a discrete local population with its own internal dynamics and that a metapopulation resulted from the interactions between local populations through dispersal of individuals. This perspective does not specify any extinction–recolonization dynamics that must arise (Hanski & Simberloff, 1997; Kritzer & Sale, 2004).

We considered local populations in a same metapopulation when they exchanged more than 5% of individuals. This threshold was derived from a general agreement in genetic studies, stating that 5% of the genetic variability was due to differences among populations, although formal evidence to support this value is still lacking. Exchanges of fish among the 98 rivers were computed based on the natal origin of individuals (i.e., rivers where individuals were born) and the destination river (i.e., where spawners migrated to reproduce). When straying occurred, fish entered a different river from their natal river to reproduce and were labeled as “allochthonous” while fish returning to their natal river to spawn were labeled as “autochthonous.” The percent of “autochthonous” fish was computed in each river for each seasonal time step. Rivers with the lowest percent of “autochthonous” were then clustered with the rivers with which they shared the most individuals. Computation ended when all groups shared less than 5% of individuals.

Experimental design

We used a complete factorial design to assess the average contribution of each input factor, to the overall outcome variance and first-order interaction using an ANOVA (e.g., Ginot et al., 2006). The complete design required 360 simulations (Figure 3; three parameters with three levels, two parameters with two levels, and 10 replications). Over the 360 ((3 × 3 × 2 × 2) × 10 = 360) individual simulations, one was retained to build the reference parameter set that best described the population dynamics of the American shad. This parameter set was used to provide simulations of the shad historical distribution and abundances in the Results.

Simulations were run over 1900–1950 as a reference period during which ecosystem shifts due to human activities were less pronounced compared with the past 50 years (Steffen et al., 2015). To initialize the model, a virtual shad population of 500,000 juveniles was introduced in each river basin, regardless of the basin size. To ensure each river basin was populated and to limit the influence of initial conditions on simulations, the evolution of the 500,000 juveniles was simulated during 100 years with constant temperature conditions (i.e., average temperatures of the 1901–1910 decade). Then, the model was run from 1900 to 1950 with updated seasonal temperature time series described in Overview of the population dynamics within the model. At each seasonal time step, the model provided estimates of shad abundances in each river basin as well as information on population status, dynamics, and spatial distribution (e.g., colonization range and number of river basins colonized).

The ANOVAs were performed using the open-source R software (R v.4.1.0; R Core Team, 2021). For each model output y, first-order indices (SI) and total sensitivity indices (TSI) (Faivre et al., 2013; Saltelli, 2008) were computed to measure the amount of variance explained by each parameter $x_i$ alone and assess their contribution with all the other inputs respectively. The first-order sensitivity index for a parameter $x_i$ measures the deviance of the factor $x_i$ and was expressed as follows:[Image Omitted. See PDF]The TSI of parameter $x_i$ described the sum of all sensitivity indices associated with the main effect of the factor and the interactions involving it.[Image Omitted. See PDF]Both SI and TSI vary between 0 and 1 with TSI being higher than SI as it includes interaction effects. In this case, SI and TSI were expressed in percentage, suggesting that the highest indices had the greatest impact on the model outputs y.

RESULTS Simulated historical distributions and metapopulation structuring

Simulations of American shad distribution were run using the set that included an Allee effect proportional to the catchment size, high homing rate associated with low straying distance, and a survival of spawners linked to water temperature. All the river basins with historical presences of American shad were accurately “populated” by the model up to the Restigouche River (48° N) under both hypotheses regarding trophic migration patterns at sea (i.e., fish stay in the vicinity of their natal river “RI” vs. fish travel to mixing zone at sea “RIO” hypothesis; Figures 2 and 4).

FIGURE 4. Comparison between (a) historical presences and absences of American shad across its native range (i.e., 1900–1950) and the simulated distributions under (b) “River-Inshore” and (c) “River-Inshore-Offshore” hypotheses, which assumed shad stayed in the vicinity of their natal river to grow or moved offshore into wintering and summering areas, respectively. Simulated distributions in (b, c) were taken from summer 1915 considering the reference parameter set. The abundances of shad in each river are represented by a continuous color scale, with each shade being an increase of 0.1 × 106 spawners. Light red indicates abundances of 0 to 0.9 × 106 and dark red indicates abundances greater than 1 × 106 spawners. Hatched river basins represent historic absences and were non-populated basins during simulations. In (a), gray basins represent historic presences.

The log-likelihood and the overall spawner abundances were higher under the “RIO” simulations than under the “RI” hypothesis (Table 1). The model predicted $9.29\pm 0.57\times {10}^6$ spawners to enter rivers across the overall range in the “RIO” hypothesis compared with $7.38\pm 0.56\times {10}^6$ under the “RI” hypothesis (Table 1).

TABLE 1 Simulated values of the three model outputs considered under both hypotheses of American shad distribution at sea, that is, shad stayed in the vicinity of their natal river (“River-Inshore” hypothesis; RI) versus shad migrated to more offshore areas (“River-Inshore-Offshore” hypothesis; RIO) after the 360 simulations, with LogL being the log-likelihood, N_m the number of emerging metapopulations, and Sp_A the number of spawners (mean ± SD).

Model outputs	Physical environment for shad at sea
	RI	RIO
LogL	$-18.2\pm 17.4$	$-12.4\pm 8.88$
Nm	$19.9\pm 7.43$	$21.4\pm 7.84$
SpA	$7.38\times {10}^6\pm 0.56$	$9.29\times {10}^6\pm 0.57$

Sensitivity analysis main effects

ANOVAs between the output variables LogL, N_m, and ${\mathrm{Sp}}_{\mathrm{A}}$ and the input parameters explained 86%, 99%, and 66% of the variability respectively under the “RI” hypothesis; and 89%, 99%, and 58% of the model deviance under the alternative hypothesis “RIO” (Appendix S5: Table S1). When considering first-order interactions, the deviance explained by the model was always higher, ranging between 59% and roughly 99% of the model variability. In both cases, the lowest deviances were observed for spawner abundance ${\mathrm{Sp}}_{\mathrm{A}}$ while the highest scores were for the number of emerging metapopulations N_m (Appendix S5: Table S1).

The ANOVAs indicated that both the straying rate and the Allee effect significantly impacted the three output variables under both hypotheses “RI” and “RIO” (p < 0.05, see Appendix S5 for detailed statistics) while the survival after reproduction only had a significant effect on LogL and the spawner abundances when considering the “RIO” hypothesis. The straying distance did not explain a significant amount of variability in any model outputs (Appendix S5: Table S2).

Influence of reproduction and straying on American shad spatial distribution and metapopulation structure

Allee effect predicted lower probability of presence in a river basin compared with simulations that included no Allee effect. The lowest LogL was obtained when the Allee effect was set to the same value across the river basins. Similarly, lower homing fidelity resulted in lower log-likelihood (Figure 5). Allee effect explained 85% of the variability of the log-likelihood (LogL) and was the most influential parameter (Figure 6 and Appendix S6: Table S1), though homing explained 12.4% as a main effect and 23.3% when considering interactions and limited movement of fish at sea (“RI” hypothesis) (Figure 6).

FIGURE 5. Variability in (a) the number of metapopulations emerging (Nm); (b) the spawner abundances (SpA) (×106); and (c) probability of presence in historically occupied rivers (LogL) depending on the four processes tested and their related hypotheses: Allee effect (constant, no effect, or size proportional), homing (low vs. high), post-spawn survival (constant or temperature-dependent), and straying distance (high, medium, and low). The 12 panels indicate the input factors used to test the sensitivity of model outputs with each modality listed along the x-axis for the “River-Inshore” (“RI”) hypothesis and the “River-Inshore-Offshore” (“RIO”) hypothesis. Results from the 360 simulations are represented by individual gray datapoints. The boxes and midline represent the first and third quartiles and the median, respectively. Whiskers symbolize the range of values, with the outliers symbolized by black datapoints at their extremities.

FIGURE 6. Total sensitivity index (TSI) values for the log-likelihood (LogL), the number of emerging metapopulations (Nm), and the overall spawner abundances (SpA) under the “River-Inshore” (“RI”) and “River-Inshore-Offshore” (“RIO”) hypotheses calculated from the ANOVAs with the Allee effect (“allee”), homing rate (“homing”), post-spawn survival (“postSpawnSurvival”), and the straying distance (“strayingDistance”). Higher percentages mean the model outputs are more sensitive to the input factor considered.

Variations in homing rate had the greatest influence (96.5% and 98.8% under the “RI” and “RIO” hypotheses, respectively) on the number of emerging metapopulations N_m (Figure 6). A high homing rate of 0.95 significantly increased the number of metapopulations $\left(\overline{N_{\mathrm{m}}}=28\right)$ compared with a low homing rate of 0.75 $\left(\overline{N_{\mathrm{m}}}=13\right)$ (Figure 5). While this trend is unsurprising since homing rate was directly linked to colonization of new habitat, the 50% reduction in the number of emerging metapopulations is disproportionally larger than the decrease in homing rate. Similarly, greater homing fidelity resulted in higher spawner abundances under both hypotheses (Figure 5). The abundances were higher when considering migration to offshore areas at sea (“RIO”). In contrast, spawner abundances decreased when Allee effect was proportional to the catchment size, or when post-spawn survival was constant rather than temperature-dependent (Figure 5). The influence of Allee effect and spawner survival on spawner abundances was stronger when fish were assumed to stay in the vicinity of their natal river (“RI” hypothesis). Together these factors explained around 23% of the variability in spawner abundances (Figure 6).

DISCUSSION

Range-wide studies integrating the complex life history of a species have the potential to improve information available to scientists and managers for decision-making (Guisan et al., 2013). By exploring likely assumptions regarding uncertain processes in American shad population dynamics, our study brought modeling support for possible causes and consequences of variability in American shad life history. Our results indicated that further investigations on the species dispersal dynamics, particularly the migration routes at sea and straying rate, may be beneficial for understanding American shad population functioning, as both factors strongly influenced the predicted distribution of spawners and the metapopulation structure and may affect management of stocks locally and coast-wide.

New perspectives regarding shad distribution at sea

In this study, we either assumed that individuals from multiple populations might group in offshore areas for growth (“RIO” hypothesis) in accordance with the assumption of Dadswell et al. (1987), or simply remain in the vicinity of their natal rivers as it was suggested for other shad species across Europe (“RI” hypothesis) (Rougier et al., 2014; Taverny & Elie, 2001). Although drastically contrasted in their ecological interpretations, both assumptions seemed plausible regarding the low level of existing knowledge on American shad distribution at sea. A third hypothesis close to the “RIO” hypothesis might also make sense with one study suggesting that fish from all the river basins travel to the Gulf of Maine together before returning to their natal rivers to spawn (Neves & Depres, 1979). Despite the evident limitations, the present study supported the most complex assumption regarding shad migrations at sea, as the “RIO” hypothesis provided better estimates of the historical species distribution. These findings suggested that a detailed characterization of oceanic migrations patterns for this species might be worth pursuing in future observational studies. More precisely, further research is needed to better characterize both the habitat use by juvenile and adult shad at sea and their abundances, identify the natal origin of individuals composing the mixed-stock fisheries occurring in the Inner Bay of Fundy (Waldman et al., 2014) and predict future changes in marine migrations with climate change (Gilligan-Lunda et al., 2021). New telemetry studies involving surgical implementation of tags into alosine have already provided encouraging results to follow continuous movement of shad histories from spawning ground to the sea and back (e.g., Davies et al., 2020; Gahagan & Bailey, 2020). Conducting more targeted studies using stable isotope ratios in tissue (e.g., Trueman et al., 2012) with different turnover rates or stomach contents coupling with ongoing tracking methods of individuals at sea such as eDNA (Boustany et al., 2008; Tillotson et al., 2018; Yamamoto et al., 2016) or otoliths as well as implementing archiving programs are actions that would be valuable for future stock assessments for which data are still limited to incidental catches only (ASMFC, 2007, 2020).

New insights on American shad metapopulation structure

From an evolutionary perspective, metapopulation structure is thought to confer advantages to populations by increasing the persistence of the species as a whole (Harrison, 1991). A population with high productivity would thus be able to “rescue” a declining population (“rescue effects”; see Carlson et al., 2014; Gotelli, 1991 for review) or even recolonize habitats in which the population would have disappeared, if prioritized for conservation (Bowlby et al., 2020).

Our study suggested a strong metapopulation structure among American shad populations (i.e., 11–30 metapopulations simulated) which was supported by moderate fidelity to breeding sites. High fidelity decreased the exchange of individuals among populations, thus increasing the number of metapopulations, whereas lower fidelity suggested higher straying and strengthened the metapopulation structure.

However, our model found more structuring in populations than the three eco-regions previously identified by genetics studies (Hasselman et al., 2013). This implicitly suggests that the current straying rate in nature may be higher than suspected, possibly driven by human activities through spawning habitat losses and fragmentation or influenced by stock transfers as it was shown for the species (Hasselman et al., 2013) and other alosines species (see McBride et al., 2015). On that, straying would be favorable in human-dominated landscape and might become a better adaptive strategy over time. Indeed, geographic proximity between rivers may facilitate straying, as it is commonly observed that more exchanges occur between relatively close rivers for several species including salmonids (Jonsson et al., 2003), estuarine-dependent weakfish (Thorrold et al., 2001), and allis shad (Martin et al., 2015) and some degree of mixing is observed within the Northeast United States for American shad. This pattern is thought to be the consequence of more similar environmental conditions between adjacent rivers (Jolly et al., 2012; Keefer & Caudill, 2014). This implies that straying is promoted by shorter distances between occupied rivers though there is now evidence of longer distance dispersal in some shad populations (Davies et al., 2020; Nachón et al., 2019). Thus, by increasing the distance to which spawners might stray in our simulations, we expected to reduce the number of metapopulations because we enabled exchanges of individuals over longer distances. In contrast, our results predicted that metapopulations were insensitive to straying distance. We noted, however, that our model incorporated 98 river basins, corresponding to a simplified topology of rivers. This simplification might bias the statistical distribution of distances between catchments for the Eastern U.S. coast toward high values and limit the influence of straying distance.

Importance of the Allee forces in population dynamics

The sensitivity analysis affirmed the influence of the Allee effect and its interactions with the homing rate on the species population dynamics. Integrating an Allee effect further constrained the simulated species distribution, leading to a lower fit to observed historical presences. This conclusion underscored the interest in considering the risk of a “depensatory” effect in population extirpations as a major research focus in the future, especially because the “depensatory” mechanism has not been acknowledged as crucial in American shad spatial distribution. Such an effect might be exacerbated considering the current low abundances of many stocks, especially in large rivers (Kramer et al., 2009). Recent studies conducted on a French allis shad population suggested that the stock-recruitment relationship including an Allee effect gave rise to a demographic collapse and, when combined with overexploitation, could explain the population collapse in the Gironde–Garonne–Dordogne (Rougier et al., 2012). Evidence for population extirpations related to an Allee effect was observed for the Northwest stocks of the Atlantic cod Gadus morhua where increasing mortality was sufficient to alter density-dependent population regulation and impede stock recovery (Kuparinen & Hutchings, 2014).

Allee forces, even if recognized as powerful drivers, do not always result in population extirpations in the short term. In some case, they are buffered through a combination of truncation of the migration season or geographic range resulting in the maintenance of populations at low levels of abundance (Stephens & Sutherland, 1999). For schooling species, such as shads, forming aggregations of individuals or rebuilding schools through mixing with other schools is thought to mitigate the expected Allee effect by maintaining higher population densities (e.g., Jorge & Martinez-Garcia, 2023).

Role of iteroparity

Greater post-spawn survival and resulting iteroparity led to higher spawner abundance throughout the distribution range. This conclusion is straightforward since we considered that survival probability of spawners increased when temperature decreased based on higher mortality in the species southernmost distribution range. For ectotherms, like fish, temperature is one of the proximate factors influencing major life-cycle processes (Angilletta et al., 2002; Huey & Berrigan, 2001). In response to environmental changes, fish might exhibit plasticity by adjusting their behavior and therefore be better adapted to a wide variety of ecological situations. An example of this is found in salmonids, for which phenotypic variability in fecundity, eggs size (e.g., Fleming & Gross, 1990), and embryo development rates (e.g., Sparks et al., 2017) were exhibited in response to thermal variability. Hence, geographic variations in iteroparity may result from either phenotypic plasticity (Pigliucci, 2006) or be direct consequences of genetic modifications (i.e., adaptation; Fox et al., 2019).

In our model, the latitudinal gradient of repeat spawning observed for the species was assessed through a water temperature-dependent survival probability, implying that spawner survival after reproduction was a plastic response to temperature variations. This assumption allowed us to partially “mimic” the latitudinal variability of repeat spawning and therefore supported the postulation that iteroparity is a plastic response rather than a specific life-history trait for American shad.

Missing drivers related to human activities

One limitation in the present study was the absence of explicit consideration of human influences on natural patterns in GR3D. Dam construction or removal, stocking, pollution or cleanup, and fishing are current examples of human influences on fish population dynamics and ecosystem processes (Brown et al., 2013; Dias et al., 2019; Mattocks et al., 2017). For instance, failure of fish passage or low passage efficiency over dams may increase the mortality of fish through migratory delay, entrapment (Keefer et al., 2012), or increased predation. Also, previous research suggested a strong relationship between downstream passages, spawner abundances, and rates of iteroparity (Castro-Santos & Letcher, 2010; Kahnle & Hattala, 2012; Leggett et al., 2004; Stich et al., 2019). Delays accrued at fishways extend the upstream migration and increase the related energy expenditures, while also leading to more passage-induced mortality (Castro-Santos & Letcher, 2010) and forced semelparity (Baktoft et al., 2020; Rubenstein, 2021; Zydlewski et al., 2021). In contrast, efforts to restore fish migration may positively influence in-river survival and iteroparity if upstream and downstream passage are efficient (Brown et al., 2013) or at least sufficiently high (Stich et al., 2019; Zydlewski et al., 2021). Thus, we expected our estimates to be conservative because factors associated with human influence were not included in our model. Nonetheless, spawner abundances simulated by our model were 13% of coast-wide abundances estimated in the literature under an historic scenario (i.e., no dam) (Zydlewski et al., 2021). Further investigation into the stock-recruitment relationship that involved a scale factor $\mathrm{\lambda }$ (defined as a model parameter; see Appendix S1) limiting the maximal number of recruits produced each year, would allow reaching the historical abundances levels reported across the rivers in previous studies. Overall, integrating such human-induced factors into modeling attempts (e.g., Stich et al., 2019) in addition to natural trends observed in pristine populations would improve understanding of the effect of human activities on population dynamics.

Possible management implications

American shad has historically been managed at the watershed scale by states and federal agencies due to the presumed high fidelity of the species to natal rivers (ASMFC, 2007, 2020). However, absence of management above single rivers can compromise population viability, local conservation investments, and effectiveness of management actions in some cases. Thus, understanding metapopulation structure is fundamental to address adequate management and conservation measures (Akcakaya et al., 2007). For instance, Regan et al. (2003) and Regan and Auld (2004) studied the effects of two major threats, seed predation and fire, on the persistence of an endangered shrub of Australian eucalypt forest. Focusing on single small population, Regan et al. (2003) concluded that predation reduction would improve the long-term persistence of shrub populations while the metapopulation approach suggested that fire management would be addressed in priority. Overall, incorporating information about the species dynamics from range-wide studies, as presented in this study, may provide new insights for decision-making.

In Europe, first applications of the GR3D model to allis shad predicted a species range-shift poleward by 2100, suggesting that allis shad would be able to cope successfully with raising temperatures imposed by climate change. These results were in line with simulations from empirical species distribution models which found newly favorable basins at the northern edge of the distribution in 2100 under the same climatic scenario (Lassalle et al., 2009). The authors also identified that the minimal temperature at which eggs and larvae survive was a key factor for the maintenance of populations at northern latitudes (Rougier et al., 2015). Some northern rivers, such as the Elbe, the Weser, or the Rhine were not populated by the model, although these rivers hosted historical shad populations (Lassalle, 2008). If underestimated, this process linked with young stages survival may affect the success of future conservation schemes involving stocking activities across the species range (Rougier et al., 2015). Thus, as the model is refined, conservation actions can be implemented with a greater understanding of the potential species responses. Comparatively, for American shad, our approach pointed toward new conclusions regarding the metapopulation aspects as suggested for river herrings (Hare et al., 2021). In particular, in simulations, the straying rate was positively correlated with the number of metapopulations, which are known to contribute to the resilience and persistence of local populations facing environmental stochasticity (Bowlby et al., 2020). Though this study provided some modeling evidence to support the causal relationship between life-history variability and population dynamics, more studies may benefit interstate management programs and conservation efforts at the scale of the species' range (Hare et al., 2021; Kritzer et al., 2022).

In conclusion, our study suggested that global approaches such as the regional and coast-wide initiative conducted by the Atlantic States Marine Fisheries Commission in the last benchmark stock assessment (ASMFC, 2020) may be widely beneficial to management. While most management actions are still directed locally, actions at national and global scales are required to mitigate the climate change effect on species stocks (Kritzer et al., 2022; Ouellet et al., 2022). GR3D explicitly integrates range-limit processes such as dispersal and enables us to forecast how changing conditions may affect population dynamics since multiple processes are influenced by temperature. As a next step, by simulating American shad distribution under climate change, we expect a northward expansion of populations based on temperature preferences and dispersal capacities that would strengthen our conclusions on the role of dispersal as a response to climate change by colonizing new suitable watersheds farther north outside state or country borders. As such, promoting large-scale and long-term management initiatives may help anticipate the impacts of climate change and increase the species viability in a globally changing environment.

AUTHOR CONTRIBUTIONS

Camille Poulet, Patrick Lambert, and Géraldine Lassalle contributed to the study conception. All authors contributed to the model conceptualization, design, and manuscript writing. A first draft was written by Camille Poulet, and then all authors commented on a previous version of the manuscript. Camille Poulet and Patrick Lambert contributed to the model development. Statistical analysis was performed by Camille Poulet.

ACKNOWLEDGMENTS

This study was co-funded by Agence de l'Eau Adour Garonne and Nouvelle Aquitaine Region within the SHAD'EAU/FAUNA projects (coordination: Françoise Daverat [INRAE]). We thank the Fundación para la Investigación del Clima (FIC; https://www.ficlima.org/) for providing the CRU database temperature extraction. We would like to thank Ryan Woodland (UMCES), Kayla Smith, and Cara Hodkin (SUNY-ESF) for their early contributions to this project. Any use of trade, firm, or product names is for descriptive purposes and does not imply endorsement by the U.S. Government.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT

Data, script, and GR3D US code version (Poulet, 2023) produced for publication are available from the Recherche Data Gouv Repository: https://doi.org/10.57745/BHLNQA.