A heuristic model of socially learned migration behaviour exhibits distinctive spatial and reproductive dynamics

Abstract

We explore a “Go With the Older Fish” (GWOF) mechanism of learned migration behaviour for exploited fish populations, where recruits learn a viable migration path by randomly joining a school of older fish. We develop a non-age-structured biomass model of spatially independent spawning sites with local density dependence, based on Pacific herring (Clupea pallasii). We compare a diffusion (DIFF) strategy, where recruits adopt spawning sites near their natal site without regard to older fish, with GWOF, where recruits adopt the same spawning sites, but in proportion to the abundance of adults using those sites. In both models, older individuals return to their previous spawning site. The GWOF model leads to higher spatial variance in biomass. As total mortality increases, the DIFF strategy results in an approximately proportional decrease in biomass among spawning sites, whereas the GWOF strategy results in abandonment of less productive sites and maintenance of high biomass at more productive sites. A DIFF strategy leads to dynamics comparable to non-spatially structured populations. While the aggregate response of the GWOF strategy is distorted, non-stationary and slow to equilibrate, with a production curve that is distinctly flattened and relatively unproductive. These results indicate that fishing will disproportionately affect populations with GWOF behaviour.

entrainment hypothesis, evolutionarily stable strategy, homing behaviour, non-stationary stock–recruitment relationship, socially learned migration, spatial population dynamics

Introduction

Movement is a critical yet poorly understood aspect of the ecology and management of animal populations (Kritzer and Sale, 2004; Nathan et al., 2008). Movement often occurs in response to short-term needs related to the availability of prey (Andrews et al., 2009), predation risk (Sheaves et al., 2015), reproductive opportunities (Warner, 1995), or environmental conditions (Johnson et al. 2002). The mechanisms underlying movement behaviour remain poorly understood and, in fact, are still being discovered (Righton et al., 2016). Moreover, because of the relationships between movement and the ecology of populations or communities, movement is fundamental to a number of contemporary management and conservation problems. As examples, failing to account for movement patterns can (1) impact fisheries management by biasing fisheries stock assessments (Cooke et al., 2016); (2) affect the efficacy of marine protected areas should fish move outside MPA boundaries (Tolimieri et al., 2009); and (3) underestimate the benefits associated with habitat restoration among spatially distinct habitats that are connected via movement (Levin and Grimes, 2002).

Seasonal migration is a notable form of movement behaviour that results in a relocation of fish on a much greater scale, and involves movement of much longer duration than those arising in its normal daily activities (Dingle and Drake 2007). Given the steep energetic costs and increased risk of predation associated with migration, numerous studies have explored the mechanisms driving migration in fishes, especially for species demonstrating site fidelity and homing (e.g. MacCleave et al., 1982; Berdahl et al., 2016).

There has been growing support for the hypothesis that migratory paths of fishes may be the result of learned behaviour. Although the importance of this mechanism is established for some whales (Valenzuela et al., 2009) and is well-accepted for migratory birds (e.g. Able and Able, 1998), its acceptance for fishes has been limited. As in the case of birds, young individuals may “learn by doing” and gain knowledge of the path and its way-points as the flock or school is guided by older, experienced migrators. Nearly 30 years ago, Dodson (1988) reviewed the role of social learning in fish migration, and observed that “the preoccupation of much fish orientation research with innate fixed patterns of behaviour on one hand and hydrodynamics on the other has led us to underestimate the possibility that orientation is a flexible process relying on developmental sequences, calibration of the motor-sensory interaction based on experience and the learning of environmental pattern.” A substantial body of evidence for learned migration behaviour of fishes has accumulated since then, including for cod (Rose, 1993), Nassau grouper (Bolden, 2000), and various reef-dwelling fishes (Colin, 1996). Petitgas et al. (2006) developed a logical framework for testing what they termed “the entrainment hypothesis.” An ICES Working Group applied this framework to a variety of fish stocks, and concluded that learned migration is a widespread phenomenon (ICES, 2007).

Given the challenges in observing migration patterns, there have been attempts at simulation modelling of migration patterns, challenged in part by lack of accepted mechanisms and algorithms. Several authors (e.g. Fagan et al., 2012; MacCall, 2012; Huse, 2016) have encouraged development of spatially explicit simulation models to explore the population dynamics associated with learned migrations. Unfortunately, detailed treatment of migration by actual fish stocks is subject to a daunting complexity of challenges (e.g. Fernö et al., 1998). Consequently, we propose that an understanding of the properties of possible mechanisms of learned behaviour can benefit from the use of demographically simplified, but spatially detailed models. To date, the only substantial attempt to develop such a model has been that of Secor et al. (2009), who presented a two-compartment spatial model incorporating learned migration behaviour, but with external designation of dominant and subordinate sites.

Herring (Clupea spp.) present a unique opportunity to extend the current understanding about learned migration in fish. Herring spawn annually in mostly known spawning sites, then migrate away from shore as juveniles to join mixed-age schools in open waters. Despite a lack of quantitative exploration of their mixing, dispersion, and migration back to historical spawning grounds, there are some leading hypotheses about these mechanisms. Learned migration is strongly evident in populations of North Atlantic herring (Clupea harengus); a convincing case was made by Corten (2002). Corten hypothesized that recruiting year classes learn the migration pattern from older fish and, furthermore, that a change in migration pattern was usually associated with a relatively low ratio of experienced older fish to first-time spawners. Huse et al. (2002, 2010) quantified and confirmed this tendency in Norwegian spring spawning herring by examining time series of age compositions. Space-time occurrence models of Icelandic herring suggest a “wisdom of the crowd” mechanism whereby larger, more effective schools lead to improved navigation to spawning grounds (Macdonald et al., 2017). Additionally, the migratory paths of North Atlantic herring stocks following recovery differed from those preceding their fishery collapses (Óskarsson et al., 2009), suggesting learned, rather than innate, behaviour.

An independent line of evidence for learned migration behaviour in Pacific herring (Clupea pallasii) is provided by traditional knowledge (TK). TK is typically transmitted inter-generationally in situ and in vivo among practitioners, and may be based on fine-scale, repeated observations. TK has proven valuable in documenting continuity and change in local species’ behaviour and ecological conditions over significant time scales (Berkes et al., 2000, St. Martin et al., 2007). TK of Pacific herring, which have long-established aboriginal fisheries, also supports the hypothesis of learned migratory behaviour. Thornton et al. (2010) explicitly noted the convergence of local and traditional knowledge with the entrainment hypothesis of Petitgas et al. (2006), and documented the traditional “importance of mature fish leading first-time spawners to various spawning grounds.” Recently, Chief Gidansta (Guujaaw) of the Haida Nation, stated that the impact of intense fishing on age structure was associated with loss of migratory knowledge: “Once herring lost the elders they lost their way to their spawning grounds” (communicated to author PL on 18 Jan 2017). Similar observations exist in indigenous oral traditions in Southeast Alaska (author HK) and British Columbia (author RJ). While it can be difficult to distinguish between spatial behaviours associated with reduced abundance and those associated with fishery impacts on age composition, the congruence of knowledge types across herring species is striking.

Understanding the mechanisms associated with migration to spawning grounds is of primary importance for Pacific herring, as many of the social, ecological, and economic services associated with Pacific herring are tied to their spatially distinct spawning grounds (Levin et al., 2016). Behavioural mechanisms driving the appearance of spawn (eggs) along coastlines should be identified, particularly in light of potential interactions with fisheries, which for some herring stocks are substantial. Migration mechanisms and behaviours have consequences for spawning patterns in species that undertake reproductive migrations, such as herring (Jørgensen et al., 2006), and for such species, identifying plausible migration and dispersion dynamics can inform ecosystem-based fisheries management.

Here we extend the work of Petitgas et al. (2006) and Secor et al. (2009) to address the behaviour and dynamics of a simple, but more spatially extensive, metapopulation by comparing two heuristic models of possible herring-like recruitment and migratory homing strategies. Both strategies employ acquired homing behaviour with some straying, and differ only in the mechanism governing how young fish are initially recruited to their spawning sites (sites that subsequently become homing targets). Importantly, neither model requires philopatry (homing to the natal site), such as is the well-known case for salmon. Habitat-quality of candidate sites is variable in space, and the role of sites as dominant or subordinate is plastic, and will be seen to be an emergent property of the model. We structure the model loosely on the spatial ecology of Pacific herring, but the model structure and results are generalizable to any age- and spatially structured species.

Methods

Model of spatial distribution of recruits

Two models of the spatial distribution of recruits are explored. Under the diffusion reproductive strategy (DIFF), pre-spawning recruits are distributed in fixed proportions among neighbouring habitats and subsequently tend to home to those spawning locations without regard to the presence of older fish. Under the “Go With the Older Fish” (GWOF) strategy, recruits join a school of older migrating fish and adopt that school’s migration behaviour, with subsequent homing behaviour oriented to the migration path and spawning location of those older fish. Thus, first-spawning GWOF recruits are allowed to occupy the same set of nearby spawning habitats as DIFF recruits, but adopt them in proportion to the relative abundances of the older fish currently using those habitats. Many recruits will return to their natal site under DIFF and GWOF, but this happens probabilistically. GWOF is different from entrainment in that GWOF describes a specific, quantifiable mechanism of entrainment that is related to recruitment dynamics, whereas entrainment is a generic pattern of migrational population behaviour.

Although we use herring as an example, we do not attempt to portray a specific fish stock. Our primary purpose is to introduce the GWOF strategy by means of a simple illustration, to describe some of its unusual properties, and to contrast it with a roughly equivalent, but behaviourally neutral diffusive model, DIFF. Greater detail with additional reproductive strategies and a wider set of model parameter values could be explored, and additional specifications could be made to apply this model to more specific populations.

The simulation model

The model is intended to be primarily heuristic and is implemented deterministically, which allows relatively clear interpretation and comparison. Extensions to allow for stochastic recruitment when this model is included in management strategy evaluation-type analysis are straightforward. The model is biomass-based, so there is no explicit consideration of fish growth or age structure, except through the time-lag between being spawned and entering the spawning population (recruits enter the reproductive stock at age 3). The biomass is distributed among a linear array of habitat sites consisting of 50 sites of varying quality (N = 50; all variables are defined in Table 1), arranged in a circle such that site 50 is adjacent to site 1. The spatial set of habitat qualities (H) for three of the 500 simulations are shown in the top row of Figure 1. Each site i in year t is occupied by local biomass B_i,t. For grammatical convenience, biomass may be referred to as being in “tons” and time is in “years,” although the scaling is arbitrary.

Table 1.

Variables used in the equations underlying the simulation models.

Symbol	Definition	Unit	Value/range of values used in simulations
N	Number of spawning sites	Spawning site	50
Bi,t	Spawner biomass for site i and year ta,b	Tons	Calculated [Equation (1)]
Si	Function of biomass that determines adult biomass after natural mortality and adult movement	Tons	Calculated [Equation (2)]
Ri,t	Biomass of recruits to spawning site i during year t	Tons	Calculated [Equations (4) and (5)]
Pi,t	Production of recruits (in biomass) at a spawning site	Tons	Calculated [Equation (3)]
β	Population density dependence	Unitless	0.005, 0.01, 0.015
Α	Population productivity	Unitless	0.25, 0.5, 0.75
Et	Exploitation ratec	Unitless	0–0.45
Hi	Habitat quality for spawning site i	Unitless	−1–1

a

${\underline{B}}_t$ denotes the vector of spawning biomass during year t (over spawning sites).

b

B_{GWOF/DIFF,I,t}denotes the GWOF and DIFF components when the population includes two components.

c

E_MSY is the exploitation rate at which expected catch is maximized.

Figure 1.

Three of the 500 simulations for alpha = 0.5 and beta = 0.01. The top row shows habitat quality (the dotted horizontal line is replacement level), the middle row shows the equilibrium spatial distribution of unfished biomass, and the bottom row shows equilibrium distribution of biomass when total biomass has been reduced to 50% of unfished biomass. The sites are arranged circularly so that S51 = S1; here the vertical line at S51 shows where the sites start again (S52 = S2, S53 = S3,…, Sn + 50 = Sn).

At each site, spawning biomass consists of surviving adult biomass, S, a small proportion of which arrives from adjacent sites; plus the number of recruits, R; less the fishery catch, which is determined by an exploitation rate $E\mathrm{ }$in the year prior. The exploitation rate is uniform across sites, but may vary in time. Thus, we calculate spawning biomass in each year at each site as:

$$B_{i,t}=S_i({\underline{B}}_{t-1})+R_i({\underline{B}}_{t-1})-E_{t-1}\hspace{0.17em}{B}_{i,t-1}$$

(1)

Applying E uniformly across all sites mimics a scenario where the fishery targets a fish aggregation away from the spawning sites, during a pre-spawning period, as is observed for Pacific herring. This also effectively simulates a life history stage where the fish migrate away from the spawning sites, and return annually to spawn. We simulate scenarios across a range of time-invariant values for E from 0 to 0.45 to understand the response of populations with distinct migration tactics to fishing mortality. The natural annual survival rate is specified arbitrarily as 70%, and for the base model 20% of adults diffuse to adjacent spawning sites. Site-specific survival and adult diffusion are combined to determine the net surviving biomass at each site i as:

$$S_i({\underline{B}}_t)=0.7(0.1B_{i-1,t}+0.8B_{i,t}+0.1B_{i+1,t})$$

(2)

For purposes of spatial accounting, the recruitment function is separated into two steps: first is production of young fish, and second is the subsequent arrival of those young fish as recruits to a set of nearby sites. All density dependence in recruitment occurs at the time of production and is site-specific in scope. We use a Ricker stock-recruitment relationship (SRR) with a 3-year time-lag between spawning and recruitment to the modelled population, a time-invariant rate of density-dependence per unit biomass that is constant among sites (β, here given a value of 0.01; see Supplementary Material for sensitivity of the primary summary statistics to values for β of 0.005 and 0.015, Supplementary Figures S1 and S2, respectively). We use a site-varying habitat quality parameter, H, for the density independent term, so the production of recruits P from spawning site i at time t is:

$$P_{i,t}=\mathrm{ }\alpha B_{i,t-3}\mathrm{e}\mathrm{x}\mathrm{p}(H_i-\beta B_{i,t-3})$$

(3)

The parameter α determines population productivity; here it was set at 0.5 so maximum yield occurs at a harvest rate somewhat near an exploitation rate of 0.2 in the DIFF model. Results are consistent using alternate values of α (see Supplementary Figures S3 and S4 showing results for values of α = 0.25 and 0.75, corresponding to alternative values for the exploitation rate corresponding to maximum yield). Results for a version of the model utilizing an alternative SRR, the Beverton-Holt SRR, are also given in the Supplementary Material (Supplementary Figure S5). Habitat quality H is time invariant, having originally been drawn from a Uniform (−1, 1) distribution without spatial correlation, with the 50 sampled values normalized to have a mean of 0 and a variance of 1/3.

Recruits arrive at a set of sites in the vicinity of and centred on their natal site i. In the diffusive DIFF model, the recruits produced at time t and site i [Equation (3)] arrive at nearby sites in a 5-site triangular distribution centred on site i, i.e. recruits to site i may be spawned at the closest 5 sites, as:

$$R_i({\underline{B}}_t)= \frac{1}{9}{P}_{i-2,t}+\frac{2}{9}{P}_{i-1,t}+\frac{3}{9}{P}_{i,t}+\frac{2}{9}{P}_{i+1,t}+\frac{1}{9}{P}_{i+2,t}$$

(4)

In the GWOF model, recruits are distributed in proportion to the adult biomasses (B_i,t) occupying those five sites:

$$R_i(t)=\frac{B_{i-2,t}{P}_{i-2,t}}{{\sum }_{j=i-4}^i{B}_{j,t}}+\frac{B_{i-1,t}{P}_{i-1,t}}{{\sum }_{j=i-3}^{i+1}{B}_{j,t}}+\frac{B_{i,t}{P}_{i,t}}{{\sum }_{j=i-2}^{i+2}{B}_{j,t}}+\frac{B_{i+1,t}{P}_{i+1,t}}{{\sum }_{j=i-1}^{i+3}{B}_{j,t}}+\frac{B_{i+2,t}\hspace{0.17em}{P}_{i+2,t}}{{\sum }_{j=i}^{i+4}{B}_{j,t}}$$

(5)

The specification of five sites in Equations (4) and (5) is intended to preserve approximate similarity in the spatial scope of the two mechanisms. Importantly, use of Equation (4) for DIFF or Equation (5) for GWOF is the only difference between the two models. For both models, recruits have a loose tendency to return to the vicinity where they last spawned or were recruited. Thus, the models accommodate an assumption of imprecise natal homing, but that assumption is not required. We also simulated diffusion in two dimensions within an array, where adults and recruits can diffuse in any of four cardinal directions from their natal site. Those simulations return similar results (see Supplementary Figure S6). Additional simulations with a broader diffusion kernel width for the recruit and adult life stages also returned qualitatively similar results (recruits disperse to 7 sites rather than 5—Supplementary Figure S7; adults disperse to 5 sites rather than 3—Supplementary Figure S8).

To simulate the model, we randomly seeded each site with a starting adult biomass and then ran the model for 700 years. Although equilibrium can be approached very slowly in these models, this number of years produced convergence in nearly all cases. No oscillatory solutions were encountered for the parameter values used here. For GWOF simulations, the final equilibrium distribution may vary with initial abundance by site, so multiple initial distributions (five) were considered for each draw of values for $H_i$⁠. These replicates captured the consequences of the initial model seed on the final distributions of biomass.

Effective spawning sites

Because of the role of habitat diversity in the ecology and management of herring, for example in portfolio effects (Siple and Francis, 2016), we are interested in understanding the impact of fishing on the diversity of habitats used by herring for spawning, and how migration strategies might mediate that relationship. One way to quantify habitat diversity is to estimate the number of effective spawning sites. The Shannon entropy index (I) is frequently used in ecology to describe the uncertainty in the identity of a species randomly sampled from a community. Here, we use it to compare spatial occupancy patterns and to calculate the effective number of sites occupied among two migration scenarios. For the occupancy of a local spawning site, Shannon entropy is

$$I=\mathrm{ }\sum p_i\mathrm{l}\mathrm{n}\left(p_i\right)$$

(6)

where p_i is the proportion of the total biomass in site i, and the summation is over all sites (Jost, 2006). This equation has the mathematical form of an expected value, in this case the expected value of ln(⁠$p_i$⁠). We use I to calculate the effective number of spawning sites, N, by exponentiating

$$N=e^{-I}$$

(7)

where N is the number of sites leading to the value of I if all sites were of equal proportion.

Testing for evolutionary stability

It is important in behavioural ecology and evolutionary game theory to establish whether a behavioural trait can be established within a population and persist evolutionarily in the presence of an alternative trait, i.e. whether it is an evolutionarily stable strategy (ESS, Maynard Smith, 1982). There are only two candidate behaviours in this set of models, but it is not clear a priori whether the GWOF strategy is a viable alternative to the DIFF strategy. For ESS testing, we constructed a combined version of the model with separate DIFF and GWOF components simultaneously inhabiting the same set of locations, and each subject to competition from the other. The dynamics of these two components are identical in all respects except for their spatial recruitment strategy; importantly, they contribute equally to each other’s density-dependent decline in local Ricker SRR reproductive success, i.e.:

$$P_{k,i,t}=\alpha B_{k,i,t-3}\exp \left(H_i-\beta {(B}_{GWOF,i,t-3}+B_{DIFF,i,t-3}\right))$$

(8)

where $P_{k,i,t}$ is the recruits from spawning site i during year t for component k (k = GWOF or DIFF). The recruits from each component are then allocated to site based on their respective model [Equations (4) and (5)] such that the only interaction between the components is in their Equation (8). The biomasses used to allocate GWOF recruits to spawning sites [Equation (5)] are the sum of the biomass of the two components. A behavioural trait is considered to be evolutionarily stable if it does not allow invasion by an alternative strategy. The ability to invade in this deterministic model can be tested by using simulations in which the initial abundance of one component is 0.1% of that of the other (with the less abundance component being the “invading” component), as well as by initially setting the two components to be of equal size. The invasion is successful if the abundance of the invading component increases. In the case that each strategy is able to invade the other, the ESS may be a Nash Equilibrium forming an evolutionarily stable set (Thomas, 1985) where the two strategies coexist in proportion such that there would be no benefit from switching strategies.

Results

The overall spatial distributions of the unfished populations are qualitatively similar for the DIFF and GWOF strategies (Figure 1), reflecting the same underlying set of H_i. For reference, a value of H = −0.51 produces replacement at near-zero biomass (Figure 1). In general, between 9 and 15 sites have values of H_i that fall below −0.51 and thus are “sinks” (range for central 90% of 500 simulations, Figure 1). Such sink sites are populated only because of immigration and recruits from other sites. Both strategies have highest densities at the most favourable sites (high H_i) and lower densities at less favourable sites, but the GWOF strategy leads to greater variation in biomass across the 50 sites (Figure 2). The DIFF strategy leads to a slightly larger aggregate unfished equilibrium biomass than the GWOF strategy for the same set of values of H_i. At unfished equilibrium, the GWOF strategy leads to higher biomass in favourable sites relative to the DIFF model.

Figure 2.

A single model simulation of herring biomass (Bi) across a range of habitat quality (Hi) for DIFF (left column) and GWOF (right column), under scenarios where the stock is unfished (top row) and fished (bottom row). Darker shading of some symbols results from overlapping symbols.

The two behavioural strategies showed distinct responses to fishing mortality. Under a constant fishing rate that reduced total biomass to half of equilibrium unfished biomass (0.5B₀), the GWOF strategy results in biomass being concentrated at a smaller number of productive sites than the DIFF strategy (Figure 1; bottom panels). The exploitation rates that maintain biomass at 0.5B₀ differ among strategies: 0.133 (median; 90% intervals 0.127–0.142) for DIFF, and 0.089 (median; 90% intervals 0.066–0.125) for GWOF, indicating a difference in stock productivity at a population size of 0.5B₀. This difference occurs despite identical habitat properties. The relative spatial distribution for the DIFF strategy differs little between the fished and unfished cases but, under exploitation, the relative distributions differ substantially for the GWOF strategy, with complete collapse at some sites and high, near-unfished biomass levels at other sites (Figure 2, lower right panel). In other words, fishing-induced declines in biomass are distributed relatively evenly under the DIFF strategy, but are highly uneven under GWOF.

Figure 3 shows the relationships between highest per-site biomass, effective number of occupied spatial sites, and catch at equilibrium as a function of fishing mortality, and biomass relative to unfished biomass. The results in Figure 3 are based on 500 sets of values for H_i (and five replicate initial distributions for GWOF; results for DIFF were insensitive to the initial distributions). Almost all of the variation in Figure 3 is attributable to the effects of the local variation in productivity rather than the implications of different initial distributions of biomass.

Figure 3.

Comparison of relative maximum biomass (maximum biomass across sites divided by the maximum local site biomass in unfished state), expressed as a percentage (a, d), effective number of occupied spawning sites (b, e), and catch (c, f), as a function of fishing mortality (a to c), and relative biomass, i.e. total spawning biomass divided unfished total spawning biomass, expressed as a percentage (d to f). The solid lines are medians from 500 replicate simulations, the dark shading covers the central 50% of the distributions, and light shading covers the central 90% of the distributions. Results are shown for GWOF in red and DIFF in blue, and are based on 500 sets of Hi values (and five initial distributions of abundance for GWOF).

The rate at which the biomass at highest biomass site declines with fishing is faster for DIFF than GWOF (Figure 3a and d). The overall population response to exploitation under the GWOF strategy is an immediate decline in the number of effective spawning sites (Figure 3b); conversely, increasing population abundance involves adding spawning sites more or less linearly (Figure 3e). Initial exploitation under the DIFF strategy shows a less immediate reduction in the number of effective spawning sites (Figure 3b) and, when recovering from severe depletion, DIFF shows a more rapid addition of spawning sites as biomass increases (Figure 3e). Conversely, under the GWOF strategy, spawning sites are added at a higher rate at high abundance and at a lower rate at low abundance (Figure 3e). The DIFF strategy leads to a production curve with a clearly defined MSY (Figure 3c and f), at a median exploitation rate of E_MSY = 0.216. The GWOF strategy results in much lower equilibrium productivity, with a median E_MSY that is 65% of the E_MSY for the DIFF strategy (0.14). The production curve for GWOF is flatter than for DIFF, with 90% of MSY achieved for exploitation rates of 0.140–0.376 for GWOF compared with 0.140–0.300 for DIFF. Maximum production occurs at a much lower proportion of unfished biomass for GWOF (27.0%) compared with DIFF (42.5%). The flattening of the production curve and the lower relative biomass at which maximum production is achieved under the GWOF strategy can be explained as follows. As overall abundance is reduced, some sites of lower habitat quality are abandoned and no longer contribute to production. At the same time, densities at the sites of highest quality (and high biomass) remain at relatively constant levels. Thus, surplus production from high quality sites remains nearly constant over a large range of exploitation rates yielding relatively small difference in catch (Figure 3c) and biomass (Figure 3f) over a large range of fishing mortalities.

Results are consistent using alternate productivity rates (α) (see Supplementary Figures S3 and S4 showing results for values of α = 0.25 and 0.75), and are even less sensitive to the density dependence parameter (β) (Supplementary Figures S2–S5). However, and as expected, the quantitative results, particular for GWOF are markedly impacted by the value for α, because lower values for α lead to a higher proportion of sites that are sinks. Results for a version of the model utilizing an alternative SRR, the Beverton-Holt SRR, are also described in Supplementary Figure S5.

The ESS analysis of the Ricker case shows that when the GWOF and DIFF strategies occur concurrently, they are each able to invade the other strategy (i.e. their populations grow when their abundance is very low relative to that of the alternative strategy), and a Nash Equilibrium occurs at a mixed population of about 30% GWOF and 70% DIFF for the set of habitats and parameters used here (Supplementary Figure S9). The relative abundances at Nash Equilibrium depend on details such as the spacing of favourable and unfavourable habitats, but this analysis indicates that both behaviours can and may even be likely to co-exist in a population. The yield function is relatively robust to the initial proportion of the GWOF and DIFF animals, but the relative proportions of each component after 700 time-steps depends, as expected from Supplementary Figure S9, on the initial proportion of GWOF vs DIFF animals (Supplementary Figure S10). A more thorough ESS analysis would require a stochastic model.

Discussion

The models incorporating the DIFF and GWOF strategies are heuristic, simplified representations of fish populations. The GWOF model is motivated by the concept of “entrainment” described by Petitgas et al. (2006) and others. However, GWOF is a specialized case of entrainment, wherein the mechanism behind the behaviour is associated with age structure, i.e. the proportion of adults in the population. The model we develop here is also motivated by scientific and traditional knowledge of Pacific herring. The information transfer in GWOF may occur passively if the young fish join and follow the older fish, learning the migration details in the process—details of that process remain poorly understood, and are omitted from our simple model. In contrast, the behaviourally neutral DIFF model is motivated by diffusion models found in physics. These contrasting models illustrate the potential for distinctive population behaviours under what we consider to be plausible alternative assumptions.

The GWOF strategy exhibits social learning of migration behaviour, which is a particular case of collective navigation (Berdahl et al., 2018). The GWOF and DIFF strategies incorporate homing to a preferred spawning site, but in neither case is that site necessarily the natal site. Thus, neither of these strategies requires an explicit assumption of philopatry to maintain a spatially coherent and productive population. In the base version of the model, DIFF and GWOF movement is limited to a nominal local window of five sites, which helps to maintain local populations. The GWOF strategy results in stronger feedback at the site level, until density dependence limits recruitment. In contrast, the DIFF strategy has weak site-level feedback, that is, the SRR is weakly density dependent, and thus emigration is higher from local areas of high biomass. The spatial distances over which either mechanism would operate in actual fish populations are presently unknown, but our explorations using various dispersal distances indicate that the observed patterns persist with reasonable alternatives (Supplementary Figures S5 and S6).

An important implication of our results is the tendency toward local site extinction when fishing interacts with the GWOF strategy. This local loss of sites is potentially overlooked by standard stock assessment procedures. Stocks with a GWOF strategy are likely to exhibit a “retrospective pattern” (Mohn, 1999) in stock assessments that assume a stationary SRR because of the slow, usually multi-decadal approach to reaching equilibrium. Productivity may be over-estimated, leading to inadvertent overfishing and unexpected further declines (Figure 1) if much of the input data for an assessment comes from a period of increasing fishing pressure. In the case of Pacific herring, where a GWOF strategy is supported by traditional knowledge, sustainable harvest rates for Pacific herring stocks should be evaluated over periods of high and low stock productivity, where the evaluation is stock specific. This is important because estimates of the exploitation rate corresponding to MSY for Pacific herring populations differ among populations and over time.

For the GWOF strategy, the typical response to increased fishing pressure is initial high productivity for few years followed by a decades-long transitional decline associated with progressive local disappearances. Numerous such declines and local disappearances have occurred in the last three decades, necessitating implementation of a minimum stock size threshold in Pacific herring management.

In recent decades, only one of the five major Pacific herring management areas in British Columbia shows evidence of being highly productive, while three have experienced persistent low production and low biomass states for a period of one to two generations (∼6–11 years). During this time, average depletion levels, defined as the ratio of estimated spawning biomass to unfished equilibrium biomass for a given year, were below 0.30 (DFO, 2017), and commercial fishery closures were implemented (DFO, 2016). Likewise, spawning no longer occurs in southerly spawning sites in the Strait of Georgia, though the overall status of the herring in the Strait of Georgia management area has not declined accordingly (DFO, 2016). As another example, herring fisheries in southeast Alaska have experienced a mix of some local increases and multiple fishery closures due to low abundance (Thornton and Hebert, 2015). Overall, the combined biomass of southeast Alaska herring stocks has increased steadily over multiple decades, with a recent decline during the last 5 to 7 years (Hebert, 2017). Several Alaskan herring stocks have shown opposite trends in biomass to those in British Columbia, despite similar harvest rate strategies (10–20% harvest rates and associated biomass thresholds; Dressel and Cleary, 2017). For instance, the biomass of Sitka Sound, Craig, Seymour Canal, and Hoonah Sound herring in southeast Alaska increased from 1980 to approximately 2010 (Hebert, 2017), whereas the biomass of herring stocks off Haida Gwaii, Central Coast, West Coast Vancouver Island, and Prince Rupert declined (DFO, 2016). The behaviour of the GWOF model indicates that it may be quite difficult to define the appropriate clusters of spawning sites that constitute a stock for purposes of evaluating overall productivity. Local sites or clusters of sites can be much more or less productive than the overall stock, and under the GWOF strategy their abundance and productivity is conditional due to processes of capture or loss of recruits from surrounding locations. In the presence of fishing, lowering of stock productivity is predicted by the progressive local disappearances (Figure 3b and d) and flattening of the equilibrium yield curves (Figure 3c and e) under the GWOF strategy.

As they appear in this model, the DIFF and the GWOF strategies both allow extensive mixing, and neither strategy would be likely to lead to sufficient genetic isolation to allow development of genetically distinct stocks even though local spawning populations appear to be geographically distinct and persistent. It may be that spawning groups under the GWOF strategy are defined by a shared migrational “meme” (defined as a behaviour or behavioural pattern passed on by non-genetic means, sensuDawkins, 2006). In this case, a stock would have no genetic basis, but would consist of a collection of mutually compatible memes. Consequently, the GWOF strategy would be able to produce many of the peculiar properties of herring populations described by Dickey-Collas et al. (2009) (complex population structure without genetic differentiation, variable migration patterns, etc.) as well as the archaeological spatial patterns of Pacific herring distribution described by McKechnie et al. (2014). The GWOF strategy also is consistent with the phenomenon of multiple stock components simultaneously using the same spawning site, as reported by Johannessen et al. (2009) for the Lindås herring (a fall spawning group in western Norway), which appear to have both resident coastal and migratory oceanic components or sub-stocks. Collective migration can emerge as the result of several mechanisms (Berdahl et al., 2018), and produce several of the spatial patterns observed in herring populations. Traditional knowledge and oral histories of the Pacific herring fishery at Sitka, Alaska, also indicate that local spawning may consist of mixed or sequential waves of separate sub-stocks (author HK).

The shape of the production curve for the DIFF strategy (Figure 3f) suggests that a conventional aggregate biomass stock-recruitment model with stationary parameter values (perhaps with an extra “shape” parameter—the local Ricker SRR shape does not necessarily apply at the population level) may be able to approximate the SRR for the entire population. The GWOF strategy presents a more difficult modelling problem in that the SRR for the overall population may not resemble the local site SRR. Locally, the proportional dispersal of the GWOF strategy leads to more recruits following adults back to high abundance sites (relative to DIFF), increasing density dependence according to the SRR. However, the overall SRR has no consistent long-term form or parameterization, i.e. it is non-stationary because the spatial distribution undergoes prolonged rearrangement in response to fishing. Although a modern age-structured stock assessment utilizing statistical deviations from an arbitrary SRR may be able to estimate recent recruitments and biomasses fairly accurately, the form and parameters of that arbitrary SRR are not reliable predictors of future recruitment or spawning distribution in space. For example, unfished biomass (B₀) is often used as a stationary reference point for calculating relative depletion. However, if the stock follows a GWOF strategy, the estimated value of B₀ is likely to be lower than the true unfished biomass because many original spawning sites may have long been abandoned, and are no longer accounted for by the fitted SRR (depending on the time series used). Moreover, under the GWOF strategy, a relaxation of fishing pressure does not lead to immediate recolonization of abandoned sites, and reestablishment at former spawning sites may require decades. Tlingit indigenous knowledge confirms that heavy reduction fishing (1925–1955) resulted in abandoned herring spawning sites in Southeast Alaska, which have taken 30–50 years to recover or have yet to recover (authors HK and TT). There are many other examples where movement ecology has important implications for conservation and management of fisheries (Westley et al., 2018).

On a millennial time scale, an “entrainment” mechanism such as GWOF has clear advantages if a fish stock is faced with low-frequency changes in habitat suitability or structure, such as the post-glaciation changes in the northern oceans: new habitats can be colonized and old habitats abandoned flexibly as changing conditions dictate. At a centennial time scale, this strategy may allow stocks to take advantage of relatively brief reproductive opportunities such as those of the episodic Bohuslän herring where very large fish populations have appeared and as quickly disappeared episodically for many centuries (Alheit and Hagen, 1997). The extremely rapid population growth and temporarily high productivity of Bohuslän herring may have been partly due to capture of recruits from other Northeast Atlantic herring stocks via a GWOF mechanism. Although herring spawning locations may vary substantially on an annual scale, on longer time scales a GWOF mechanism may provide an indirect way to track approximate relative suitability of spawning habitats, approximating density-dependent habitat selection and approaching an ideal free distribution with associated “basin model” population behaviour (MacCall, 1990), even if the individual fish have no direct means of evaluating habitat suitability. At interdecadal time scales, this population behaviour should lead to abundance-dependent expansion and contraction of geographic range similar to that predicted by MacCall’s (1990) model, as an essential feature of the GWOF mechanism is that there is positive feedback (reinforcement) of population abundance in proportion to the abundance of adult fish. For example, this is consistent with the spatial pattern of recovery of the North Sea herring observed by Schmidt et al. (2009), where abundance at traditional core sites was maintained, but with a growing contribution from minor or newly recolonized peripheral sites.

The GWOF strategy may have a further benefit that is not captured in the simple model of this article. If some individuals are subject to adverse transport during early life, which seems inevitable, a portion of the reproductive output may inhabit locations that are non-viable for reproduction. The DIFF model is relatively mechanical, and has no mechanism to recover such straying individuals. The resulting mortality would be incorporated in the density independent term of the local SRR. In contrast, the GWOF strategy may allow the adults to behaviourally “sweep up” a portion of those straying recruits, possibly with increasing efficiency at higher adult densities. Logically, the GWOF strategy would be likely to have a somewhat lower density independent mortality rate in the SRR, but the magnitude of any such advantage would be case specific. The GWOF strategy also potentially carries with it benefits often associated with reproductive homing, i.e. increased foraging and reproductive success, and predator avoidance (Brown and Dreier, 2002). These benefits may have led to the persistence of GWOF behaviour even in populations where many individuals were following a DIFF strategy, a pattern confirmed by our ESS analysis.

Previous work has demonstrated that socially learned migration is sensitive to shifts in population structure (Fagan et al., 2012). Our simple biomass-based model does not contain sufficient age structure to explore important demographic details of learning and migration. These details could, however, be incorporated in an operating model used for management strategy evaluation (Punt et al., 2001). A full model would include case-specific treatment of fundamental life history parameters, including stochastic properties and spatial structure. Additionally, a full GWOF model would require tracking the ratio of inexperienced to experienced migrators as done by Huse et al. (2002, 2010), or even the number of migrations previously undertaken by a cohort if the strength of migratory navigation increases with multiple migratory experiences. Because exploitation leads to shorter expected lifespans and risk of degraded transmittal of migratory knowledge, inclusion of these demographic details may be especially important in attempts to determine optimal management policies.

Our models suggest that the GWOF strategy, and to a lesser extent the DIFF strategy, offers additional dimensions for management consideration. Long-term management performance will vary not only with target fishing intensity, but with its timing in the annual cycle (e.g. fishing on mixed feeding schools vs. site-oriented harvest of reproducing schools) and presumably also with alternative spatial distributions of harvest. Our GWOF model indicates that productivity of fished populations may decline over decades due to spatial reorganization of the population. Our analysis indicates that fishing reduces the effective number of spawning sites used if the GWOF migration strategy is prevalent in a spatially structured population, shown here when biomass is reduced to half of the unfished level (Figure 3). Management emphasizing maximum yield may need to balance preferential harvest of high-productivity sites against risk of site abandonment, or identification of target population biomass levels that maintain the effective number of spawning sites at some desired number, while recognizing that even with no harvest, individuals may abandon spawning sites for unknown reasons. If the GWOF mechanism is operative, preservation of spawning sites that are valued socially or culturally (or more generally, the establishment of marine protected areas) may merit consideration not only of the possible source or sink properties of that site, but also of the exploitation patterns at other interacting sites. Another topic of current relevance is the response of fisheries and resources to climate change. The GWOF strategy should enable a population to adapt to changing patterns of productivity at potential spawning sites, provided that fishery exploitation or increases in natural mortality rates do not prevent the adoption of new reproductive sites or interfere with inter-generational transmission of migratory knowledge.

The GWOF hypothesis provides a simple and transparent mechanism for learned migration behaviour or “entrainment.” In addition to GWOF providing an extension to the body of published scientific literature on learned migration behaviour, it also is based on an independent body of traditional knowledge that often remains unrecognized or under-used in modern quantitative population modelling. We have shown that including GWOF in spatial models of population dynamics can generate spatial patterns and dynamics that are distinct from those produced by behaviourally neutral diffusion, and, moreover, that these distinct patterns resemble those exhibited by Pacific herring and other clupeids. The GWOF hypothesis provides insights into population and fishery behaviour that should be useful to the assessment and management of Pacific herring. The GWOF model generates a wealth of hypotheses for future testing, analysis, and exploration.

Acknowledgements

The Packard Foundation, through its Ocean Modelling Forum, and the Pew Charitable Trusts (Ocean Science Division and the Pew Forage Fish Conservation Initiative) provided funding for TF, AP, and PL. The comments of five reviewers and the editor improved the manuscript and are much appreciated. We especially thank the Haida Nation and the Sitka Tribe of Alaska for their collaboration and engagement with their traditional knowledge.