The bay scallop Argopecten irradians, a short-lived bivalve mollusk, has three subspecies that occupy estuarine and coastal environments along the U.S. Atlantic and Gulf of Mexico coasts, usually in or near beds of submerged aquatic plants such as eelgrass Zostera marina (Brand 2006). The northern subspecies, the northern bay scallop A. irradians irradians, extends from Massachusetts to New Jersey, the southern subspecies, A. irradians concentricus, is found intermittently from North Carolina (under dispute) and again along the west coast of Florida to Louisiana, and the third subspecies, A. irradians amplicostatus, exists sparsely in the Gulf of Mexico along the coasts of Texas and Mexico (Clarke 1965). Partially because of the bay scallop's short lifespan, population size and commercial landings have had extreme yearly fluctuations, and their vulnerability to environmental and anthropomorphic change has led to a general decline throughout the species’ range (MacFarlane 1991; MacKenzie 2008; Tettelbach and Smith 2009). Nantucket, Massachusetts, is one of the few remaining locations where a wild fishery of bay scallops persists, but Nantucket landings from 1966 to 2010 showed similar decreasing trends.

In addition to an early spawn in late May through early July, a second seasonal spawn occurs frequently in the Nantucket Island (Massachusetts) northern bay scallop population. The early spawn results in a “classic” juvenile scallop that reaches a shell height of 40–55 mm before overwintering. It then forms a growth ring at that position in the spring, spawns twice during the next summer at year 1, and enters the fishery in the fall. Classic scallops will almost always die before entering another spawning season. The second spawn in late summer or fall produces a scallop known to islanders as a “nub,” with its growth ring close to the hinge. Juvenile scallops produced by the second spawn reach shell heights averaging 10 mm during their first autumn before growth ceases for the winter (Tettelbach et al. 2001). Growth then resumes in the spring and is marked by the formation of the first (nub) growth ring ≤ 10 mm from the umbo. This nub ring is the only growth ring found in late-spawned scallops until they form another one after surviving their second winter. While less than 40% of nub scallops spawn in their first summer, at least 50% are able to survive and spawn during their second summer before dying late in the season (Hall 2014). Nubs in Nantucket scallop population are subject to harvest in their second winter if they exceed 63 mm in shell height; however, they rarely survive into a third commercial season.

The relative frequency and existence of late-spawned northern bay scallops appear to vary greatly. Belding (1910) estimated that 10–20% of Massachusetts bay scallops survive to spawn a second year, suggesting that those living longer were born late in the first spawning season. Kelley (1981) reported adult Nantucket bay scallops spawning into the month of September. MacFarlane (1991) found a 9% occurrence of “ring-at-hinge” scallops in Pleasant Bay, Massachusetts, during fall 1979. Taylor and Capuzzo (1983) inferred fall spawning in bay scallops in Waquoit Bay, Cape Cod, Massachusetts, as did Tettelbach (1991) in Groton, Connecticut. Juvenile scallops < 20 mm averaged 2.59% of the populations in some parts of the Peconic Bays, Long Island, New York, during the winter of 1990–1991 after a nonbrown tide year (Tettelbach et al. 2001), but by October 1992, after a brown tide bloom, adults with growth rings 2–7 mm from the hinge approached 100% in one bay (Tettelbach et al. 1999). These may have been the result of the fall spawn of 1991. Bologna et al. (2001) observed settling juveniles of less than 15 mm shell height in October of 1998 and 1999, which strongly suggests there is a fall spawn in New Jersey scallops.

Leslie (1945) first introduced matrix population models to ecologists. While scalar population models rely only on counts or estimates of total population size, matrix models can describe populations in terms of ages, size-classes, or stages of development (Caswell 2001). They are now readily implementable to computer programs such as MATLAB. The matrix elements, or vital rates, of Leslie's model are derived from a life table of an organism and a table of age-specific fertility rates, and consist of fertility rates across the top row and survival probabilities from one age to another along the principle subdiagonal. The population growth rate (λ) can be determined from the principle eigenvalue of the Leslie matrix. The λ term is related to the more familiar r (intrinsic rate of natural increase) by the equations λ = e^r and r = lnλ. The right and left eigenvectors of the Leslie matrix represent the stable age distribution and the reproductive value of each age-class, respectively. An age-specific Leslie matrix can be multiplied by a column vector of the initial age distribution of the population to achieve a projection of future population structure. (Caswell 2001).

While Leslie (1945) constructed his matrices based upon age-classes, Leftovitch (1965) realized that some species could not easily be divided into age-groups. Rather, stages are often more recognizable. He added the vital rate of growth of one stage to another to the vital rates of fertility and survival. Two studies that used stages are Brault and Caswell (1993) and Smith et al. (2005). The former researchers designed a matrix model with four stages (yearlings, juveniles, reproductive adults, postreproductive adults) to study the demography of killer whales Orcinus orca. The latter evaluated variability in flood and precipitation as they affected the demography of the threatened floodplain plant, Boltonia decurrens, along the Illinois river, using seeds, seedlings, and various reproductive stages. Many others have applied Leftovitch's model to sizes rather than stages. Werner and Caswell (1977) found that size was a better predictor of population dynamics than age in the teasel Dipsacus sylvestris. Doak et al. (1994) established sizeclasses of the desert tortoise Gopherus agassizii in the Mojave Desert to evaluate various proposed management regimes for this threatened species. Although Nakaoka (1997) had previously used an age-based matrix model to study the demography of the infaunal clam, Yoldia notabilis, in Japan, he constructed a size-classified stochastic matrix model to assess the effect of fluctuating recruitment rates on the population. He did so because survivorship and reproduction in this species is more dependent upon size than age; it is easier to incorporate annual growth fluctuations into the model, and older age-classes with smaller sample sizes can be combined into one large size-class.

“Multistate” matrix models, either age- or sizeclassified, can be used to classify individuals and to follow their transitions either by geographic regions or by demographic group (Caswell 2001). They have often been applied to metapopulation models, where subpopulations are connected through larval dispersal, diffusion, or migration, or multiregional models when applied to human demography. For example, Wootton and Bell (1992) analyzed the response of the endangered California peregrine falcon Falco peregrines population to different management strategies by developing a model with two subpopulations linked only by migration, Barbeau and Caswell (1999) applied a multistate model to evaluate the effectiveness of seeding strategies of juvenile sea scallops Placopecten magellanicus, and Strasser (2008) applied the model to a theoretical two-patch metapopulation of the soft-shell clam Mya arenaria to evaluate the relative importance of “sources” and “sinks” in fishery management. When the models are used to combine multiple demographic classifications, they are specifically called multistate or multidimensional models, and the same mathematical principles apply to both (Caswell 2001).

Matrix population models were first applied to a marine mollusk, M. arenaria, in Ipswich Bay, Cape Ann, Massachusetts (Brousseau 1978a, 1978b; Brousseau et al. 1982; Brousseau and Baglivo 1984). After determining the clam's spawning cycle and gathering size-specific fecundity and mortality data, an age-based matrix model was formulated. Survival probability was divided into two factors: settlement rate (r_s) and survival from settlement to year 1 (b₁). Population growth rate in that study was most sensitive to changes in settlement rate, which is mainly composed of survival. The second most important influence on population growth rate was survival to age 1. Ripley (1998) investigated another population of M. arenaria in relatively pristine Barnstable Harbor, Cape Cod, Massachusetts, and compared it with one growing in the contaminated sediments of Boston Harbor. Ripley (1998) found that variability in settlement rate in this long-live species, due to the occurrence of occasional years with exceptionally high recruitment, was responsible for the persistence of the population. This conclusion was confirmed by Ripley and Caswell (2006) using stochasticity analysis of clams in Barnstable Harbor. Sporadic recruitment success in a species with broadcasting larvae becomes more important as environmental variability increases.

Other population studies of marine mollusks that have successfully employed matrix models are those of Malinowski and Whitlatch (1988), Weinberg (1989), Noda and Nakao (1996), and Barbeau and Caswell (1999). Malinowski and Whitlatch (1988) used matrix models at different larval settlement rates to analyze life history tactics of M. arenaria, the oyster, C. virginica, and the hard clam, Mercenaria mercenaria, then assessed the relative benefit of different management strategies on various stages of their life cycles. Since population growth rate was most sensitive to the survival of larval and juvenile stages (mainly due to predation), Malinowski and Whitlatch (1988) concluded that the best management practices were those that increased the survival of those stages. Weinberg (1989) used cohort-specific growth rates for a population of infaunal clams, Gemma gemma, in Little Narragansett Bay, Rhode Island, to make annual forecasts of population size from an age-based matrix model. Weinberg (1989) found that the model was realistic enough to describe the demography of this population, particularly because the species has little migration and no dispersing larval stage. Noda and Nakao (1996) reported the first documented case of recruitment limitation in a whole population of marine benthic animals using the Japanese snail, Umbonium costatum. Those investigators varied recruitment rate only in eight independent Leslie matrices that described the transition in age from age 0 to over 6 years, combining them to simulate long-term population dynamics. Barbeau and Caswell (1999) divided sea scallops into 10 size-classes, each with a range of sizes from juveniles to harvestable adults, to elucidate the best strategy by which to seed juveniles. They compared the contributions of predation and dispersal to the growth and survival of the juveniles, which had been seeded in four locations along the coast of Nova Scotia, and identified possible management strategies. Barbeau and Caswell (1999) concluded that reductions in predator densities would have the most dramatic effect on scallop survival, especially if the size of seeded scallops was increased.

In this study we addressed the effects of environmental stochasticity on the bay scallop, with particular interest in the role that its short life span plays in its population dynamics. Some studies have found that shorter life histories are most often affected by environmental stochasticity (Benton and Grant 1996; Tettelbach and Smith 2009). However, such longlived mollusk species as the gastropod, U. costatum, (Noda and Nakao 1996), the protobranch bivalve, Y. notabilis, (Nakaoka 1996), and the sea scallop (Barbeau and Caswell 1999) are also vulnerable to environmental variation, especially in certain life history stages. Random environmental fluctuations alter the vital rates in stochastic matrix models compared with the constant rates in deterministic models. Information attained from a deterministic matrix model (growth rate λ, a measure of fitness, reproductive values, stable age distribution, elasticity, and future projections) assumes a constant environment, but this assumption is often unrealistic for natural populations (Gourley and Lawrence 1977). Orzack and Tuljapurkar (1989) examined 25 life histories with differing patterns of variability and found a lognormal distribution of random environmental variation. They assumed that variability in vital rates depends upon the way in which organisms experience these changes.

Most mollusk species have life spans longer than that of the bay scallop's 2 years: Arctica islandica (up to 400 years); Panopea abrupta (up to 200 years); M. mercenaria (up to 40 years); Placopecten magellanicus, C. virginica, and Y. notabilis (up to 20 years); Mya arenaria (10–12 years); C. islandica and U. costatum (17 years); G. gemma (8 years) (Abele et al. 2009). Constructing and analyzing a matrix population model for such a short-lived species as the bay scallop is a challenging task that has not been undertaken, except for an initial attempt by the U.S. Environmental Protection Agency (Hinchey et al. 2004). The further addition of the demographic effects of two different bay scallop life history patterns, classic (spawned early in season) and nub (spawned late in season), has never been attempted.

We constructed a multistate, age-based, Leslie matrix population model to investigate the contribution of late summer and fall spawning to population dynamics of the bay scallop. Projection matrices were constructed with and without the contribution of nubs, and then using initial conditions for both good and poor nub years (based on population estimates from September surveys conducted from 2006 to 2010). Population growth rates, stable age distributions, and reproductive values were determined from each projection. Total population estimates for 10 years at 1-year intervals were compared between life histories with and without a second seasonal spawn. Elasticity analysis revealed that population growth rate was sensitive to changes in vital rates in the projections, and stochasticity analysis modeled the effects of a randomly changing estuarine environment on bay scallop population dynamics. If our analysis shows that spreading reproduction out to include a second seasonal spawn does enhance bay scallop population growth, our matrix population models can be used as powerful tools to assist in the management of this critically important shellfish resource.

METHODS

Model formation and analysis.—We used a Leslie agebased matrix model with the addition of multistate and stochastic analyses. The six stages used in the bay scallop model were: (1) classic (early spawned) newly settled juvenile, (2) year-1 classic adult, (3) year-2 classic adult, (4) nub (late spawned) newly settled juvenile, (5) year-1 nub adult, and (6) year-2 nub adult.

In our multistate matrix model, reproductive products did not move physically from one region to another, but from one life history strategy to another. Our multistate matrix (M) consisted of two “regions,” which represented classic (early spawning only) and nub life histories (late spawning only), with four submatrices: (A₁) simple classic life history giving rise only to classics (thus remaining in that region), (A₂) simple nub life history giving rise only to nubs (thus remaining in that region), (M_1→2) year-1 classic late spawning that gives rise to nubs, and (M_2→1) year-2 nub early spawning that gives rise to classics (matrix, Figure 1a; life cycle graph, Figure 2). We compared a simple age-classified Leslie matrix model for classic bay scallops with a multistate one, basing its construction on that of Barbeau and Caswell (1999) for sea scallops in Lunenburg Bay, Nova Scotia. Initial deterministic Leslie matrices for early spawned classics only (submatrix A₁), latespawned nubs only (submatrix A₂), and the multistate twospawn model (M), not taking into account environmental variation or uncertainty in assumed parameters, were used to estimate eigenvalues. The principal eigenvalue, i.e. the largest real positive root, represents the population growth rate (λ) (Caswell 2001). The right eigenvector (w_i), indicates stable age distribution at time t ± 1, and the left eigenviector (v*_i), indicates reproductive value of each age-class.

FIGURE 1.

(a) Generalized multistate population matrix (M) of the Nantucket bay scallop population. A₁ is the basic classic life history matrix; F_2E and F_3E are fertilities of stages 2 and 3 during the early spawn only. A₂ is basic nub life history matrix; F_5L and F_6L are fertilities of stages 5 and 6 during the late spawn only. M_1→2 is the contribution of classics to the nub life history; F_2L is the fertility of stage 2 during the late spawn only. M_2→1 is the contribution of nubs to the classic life history; F_6E is the fertility of stage 6 during the early spawn only, (b) Multistate matrix with proposed values for vital rates substituted into the matrix described in (a).

FIGURE 2.

Fife cycle graph for the Nantucket Harbor bay scallop population. A₁ is the classic life cycle; stages 1, 2, and 3 are year-0, year-1, and year-2 classics, respectively. A₂ is the nub life cycle; stages 4, 5, and 6 are year-0, year-1, and year-2 nubs, respectively. Crossovers between life cycles represent contributions of one life cycle to the other: F_2L is the second (late) spawn of year-1 classics; F_6E is the early spawn of year-2 nubs.

TABLE 1.

Vital rates chosen for the matrix population model for the northern bay scallop. Fecundity is number of eggs per individual per spawning event. Fertilization rate, larval survival, and settlement rate, expressed as percent, are multiplied with fecundity to yield Fertility (F) for the top row of the matrix. Probability (P) of juvenile and adult survival is expressed as percent for the major subdiagonal of the matrix.

Vital rates.—Values chosen for the vital rates of fertility (F) and probability of survival (P) in our model are described below, summarized in Table 1, and diagrammed in Figure 1b.

To calculate F for the model, we multiplied our estimated values for fecundity, fertilization rate, larval survival, and settlement rate and applied them to each age-class of early spawned classics and late-spawned nubs (juvenile, year-1, and year-2) (Table 1).

The first stage in our bay scallop model (age-0) began with recently settled spat (juveniles), because this was the first stage when scallops could be quantified. Fertility thus involved several component processes spanning from initial spawning through larval settlement. Because of the difficulty in assessing bivalves at an earlier time, matrix models often begin with this stage (Ripley 1998; Ripley and Caswell 2006).

Fecundity is defined as the number of eggs produced by one individual per spawning event (Llodra 2002; Barber and Blake 2006). Values in the literature for the bay scallop range from 2 × 10⁶ (Belding 1910) to 23.7 × 10⁶ eggs (Bricelj et al. 1987) per spawning. We chose the more conservative rate of 2 × 10⁶ (Belding 1910) for year-1 classics because it is still considered a valid indicator of egg production in the wild (Tettelbach et al. 2011). There was no difference in fecundity during the first and second spawning periods of year-1 classics during 2010, the only year tested (ANOVA: P = 0.64, F = 0.2167, n = 40) (Hall 2014), so we used the same value for both. We assumed that 100% of year-1 classics released eggs during each period. Only 39% of year-1 nubs released eggs and then only in the late spawning period (Hall 2014), so we calculated their fecundity for the model as zero for the early spawn and 7.8 × 10⁵ (39% of classic fecundity) for the late spawn. Bricelj et al. (1987) calculated the average fecundity of a year-2 scallop as 8.6 × 10⁶. However, we observed that only 9% (SD, ±4.3%) of ovaries of year-2 classics examined histologically showed any spawning and only in the early spawn (Hall 2014), so we calculated their fecundity as 1.8 × 10⁵ (9% of classic fecundity) in the early spawn and zero in the late. Year-2 nubs were at least twice as fecund as year-1 classics during the early spawn, but only 4% spawned in the early period (thus, a fecundity of 4% of 4 × 10⁶, or 2.4 × 10⁵), but equal in fecundity to classics during their second spawn, though only 44% of them spawned (thus, 44% of 2 × 10⁶, or 8.8 × 10⁵).

The mean percentage of eggs fertilized for our model was based upon the method of Tettelbach et al. (2013), who determined the regression equation from data provided by Lundquist and Botsford (2004) Figure 7 :

where y = mean percent of eggs fertilized and x = density of spawners per square meter.

We calculated the density of spawners by using the population density of the survey year with the largest percent of adult scallops (0.89 scallops/m² in 2007: Hall 2014) and the mean survival of adult classic and nub scallops from two unpublished caged studies (54%). The resulting spawning density was 1.65 scallops/m², determined by dividing the mean density of adults from the 2007 September survey (0.89 scallops/m²) by the mean survival through the summer spawning period (0.54).

In broadcast spawners such as bay scallops, fertilization success (percent of eggs released that actually become fertilized) also depends upon tidal current flow. Pennington (1985) showed that percent fertilization in green sea urchins Strongylocentrotus droebachiensis was lower in fast currents (>0.2 m/s) than in slow ones (<0.2 m/s). Lundquist and Botsford (2004) used that study to develop two simulated negative exponential sperm distributions for their model: a broad exponential distribution for slow current speeds and a narrow exponential distribution for high current speeds (see their Figure 7). Using Lundquist and Botsford's (2004) model, Tettelbach et al. (2013) employed the broad exponential sperm distribution based on current speeds of <0.2 m/s recorded in the Peconic Bays to determine percent of eggs fertilized (i.e., fertilization success). We used the narrow exponential sperm distribution after determining that mean current speeds in Nantucket Harbor (0.43 m/s) exceeded 0.2 m/s (P. Boyce, Maria Mitchell Association, unpublished data). We also assumed a slightly aggregated spawner distribution as did Tettelbach et al. (2013). Using a spawner density of 1.65 scallops/m² for Nantucket Harbor, we calculated a mean fertilization rate of 10.29% for Nantucket Harbor (equation 1).

Larval survival has not been directly estimated for bay scallops, but is thought to be very low. Vance (1973) predicted that any invertebrate with a pelagic larval stage will experience survival rates below 1% during that time. Previous studies found that the average for four scallop species was only 0.54% (Vahl 1981; LePennec et al. 1998; Martinez et al. 2007; Soria et al. 2010) and that for seven other bivalve species was 0.26% (Strathmann 1985; Hines 1986; Malinowski and Whitlatch 1988; Rumrill 1990). We chose a larval survival rate of 0.35% for the model because that was the rate calculated by Soria et al. (2010) for the rock scallop Spondylus calcifer, whose 14-day larval period is similar to that of the bay scallop. We chose a settlement rate of 1%, which is an estimated maximum for marine invertebrates with planktotrophic larvae to survive metamorphosis and settlement (Thorson 1966).

Probability of survival (P) was calculated for two postsettlement stages: (1) juvenile to age 1, and (2) age 1 to age 2. Survival of nubs and classics was based on both the literature and on results of earlier unpublished studies.

Juvenile survival varies widely depending upon whether or not the scallops are caged. Caged juvenile classics and nubs in Nantucket Harbor had a similar survival of ∼ 75% from late fall (early postsettlement period) through late spring (prespawning period) (Hall 2014). In order to obtain an estimate of survival to be used in the population model, we compared those results with ones found in the literature. Wild bay scallop juvenile survival ranges from 25% (Irlandi et al. 1995) to 70% (Tettelbach 1990) for classics and 16.5% for nubs (Tettelbach et al. 2001). The latter figure represents an average of two values reported from two locations in the Peconic Bays after the winter of 1990–1991. Survival values of other juvenile bivalves range from 1% in My a arenaria (Ayers 1956) to 50% in Mercenaria mercenaria (Zarnoch and Schreibman 2008). Brousseau et al. (1982) surveyed the literature for 11 species and found that the average juvenile survival in the wild was 17%. We assessed the wild, unfished population of juvenile scallops, so we assumed 40% survival for classics, an average of the wild bivalve data mentioned above. Although Tettelbach et al. (2001) observed 16.5% mortality for nubs, we assumed, based on our caged studies, that their survival was equal to that of classic juveniles (40%). We also assumed no fishing mortality in juveniles since no fishing occurs until after spawning at year 1.

Adult survival from year 1 to year 2 was assessed in both caged and wild scallops in Nantucket (Hall 2014), with caged nub adults exceeding survivorship of classic adults (39% versus 6%) over 2 years. Based upon our survey results, mean survival of wild scallops from age 1 to age 2 in Nantucket Harbor was 2% for classics and 13% for nubs (Hall 2014), the difference presumably was due to fishing pressure on classics. Extremely low salinities and burial in sediments in winter can further reduce wild scallop survival at times (Tettelbach et al. 1985; Tettelbach 1990). For the population model, the probability of year-1 classics surviving to year 2 was assumed to be 2%, while the same probability for nubs was assumed to be 13%, both based only on our field data. It was not possible to separate fishing and natural mortality in making these assumptions.

Population projections were first made from a simple 3 × 3 Leslie matrix of early spawned classic scallops only (submatrix A₁) then from the 6 × 6 multistate matrix (M), adding the effect of late spawning. For comparison, matrix A₂ (late spawning only) was also analyzed. Initial abundance vectors shown in Table 2 for classics only, nubs only, and multistate with both contributions were multiplied by their corresponding matrices to produce population projections 10 years into the future at 1-year intervals, using the equation:

where A is the projection matrix and n_t is the initial population vector.

Initial vectors were based on estimations for each age-class in the harbor, determined from population surveys (authors’ unpublished data). The year 2008 was chosen as the “low-classic-high-nub year,” when 7,921,177 (39.3%) of the estimated scallop population in Nantucket Harbor were year-1 nubs, while year-1 classics made up only 3.6% of the population. The “high-classic-low-nub year” was 2009 when an estimated 1,495,072 (10.3%) were year-1 nubs, while year-1 classics made up 73.9% of the population.

TABLE 2.

Values for initial vectors (n_t) which were multiplied by the Leslie matrices (A) to produce population projections using the equation n(_{t+ 1}) = An_t for 2008, a “low-classic—high-nub year,” and 2009, a “high-classic—low-nub year.”

We made an approximation of the carrying capacity of Nantucket Harbor by comparing the commercial scallop landings (in bushels) during a survey year where the vast majority of scallops were adults (2007) with those of a previous year (1980), which had the maximum recorded landings for the harbor. We estimated the total population of bay scallops in the harbor in 1980 was 110 million individuals by setting up an equation comparing the ratio of the 2007 estimated scallop population to that year's commercial landings with the landings in 1980 (assuming an average of 400 scallops/bushel). The two adjacent years, 2008 and 2009, were compared with or without a second spawn to see whether its addition made a significant difference to the population growth rate and to the time during which each would reach an approximation of carrying capacity of the harbor.

The matrices were then subjected to sensitivity analysis to test the effects on population growth rate of changing the vital rates. The sensitivity matrix is calculated by multiplying the right and left eigenvectors together and reveals the sensitivity of the population growth rate (λ) to changes in rates of fertility and survival. Sensitivity of λ to both vital rates tends to decline with age (Caswell 2001). It is more straightforward to express sensitivities in an elasticity matrix rather than by using absolute numbers, because its elements are proportional (ranging from zero to one). Elasticity can be calculated by multiplying the corresponding elements of the sensitivity matrix by those in each original projection matrix and dividing their product by λ (Caswell 2001). We prepared elasticity matrices for classics only (submatrix A₁), nubs only (submatrix A₂), and the full multistate model with both classics and nubs (matrix M).

Stochasticity analysis tested the effects of a varying random environment on the population growth rate. The vital rates P₁ (juvenile survival) and F₂ (classic adult fertility) were first each varied alone in a stochastic model and then varied together. The stochastic growth rate (λ_s) is the dominant eigenvalue calculated for each stochastic matrix. The mean of the stochastic growth rates () resulting from 2,000 iterations (and its CI) was calculated for each case (A₁ [classics only] and M [multistate with both classics and nubs]). The stochastic growth rate is a function of the intrinsic rate of increase (λ = e^r) and the mean generation time, and represents fitness in a randomly varying environment (Orzack and Tuljapurkar 1989; Caswell 2001). One should use a gamma (γ) or lognormal distribution for fecundity and a beta (β) distribution for survival in order to randomize the vital rates for stochastic simulations (Benton and Grant 1996; Caswell 2001). We drew 2,000 vital rates for stochasticity analysis randomly from a (1) β distribution of possible values of the vital rate to which population growth rate was most sensitive (P_1; early juvenile survival) and (2) lognormal distribution of the next highest rate (F₂, year-1 classic fertility) (Benton and Grant 1996; Caswell 2001).

Values of P₁ following the β distribution ranged from a minimum of zero to a maximum of 1.0. This early juvenile survival covered the period from early postsettlement to early summer of the following year, when the scallops would enter the reproductive population. It was thus equivalent to the bay scallop's recruitment rate. Although the observed survival rate of juvenile bay scallops maintained in published field studies and subject to predation, but not to fishing, are 44% to 93% (Pohle et al. 1991; Irlandi et al. 1995), there have been reports of nearly unsuccessful recruitment of early spawned (classic) juveniles in situations such as during the 1995 brown tides caused by Aureococcus anophagefferens in the Peconic Bays (Tettelbach et al. 1999) and the 2009 rust tide outbreak caused by Cochlodium polykrikoides in Nantucket Harbor, Massachusetts (Hall et al. 2011).

Values for F₂ following the lognormal distribution ranged from a minimum of 0.777 to a maximum of 2,150 eggs/individual per spawn. These values were determined by multiplying together the lowest figures for fecundity, fertilization rate, larval survival, and settlement rate obtained from the literature (Belding 1910; Brousseau 1978b; Soria et al. 2010; Tettelbach et al. 2011) and then doing the same with the highest figures (Thorson 1966; Bricelj et al. 1987; LePennec et al. 1998; Tettelbach et al. 2011).

The following equations were used to determine the parameters of the two distributions:

(1) P₁ β distribution

where µ = 0.4 is the mean of P₁ (Table 1) and σ = 0.25 is the SD, taken as one-quarter of the range of P₁.

(2) F₂ lognormal distribution

where m = 7.2 is the mean of F₂ (Table 1), v is the variance of F₂, and σ = 1.98 is the SD of ln(F₂) taken as one-quarter of the range of ln(F₂). By rearranging equation (4), parameter µ was calculated using the formula

RESULTS

Deterministic Model

The upper left quadrant (A₁) of the matrix (Figure 1a) signifies the contribution of the classic life history alone to the bay scallop population. Its principal eigenvalue (λ) indicated that the population growth rate was 1.6979, implying that the population was expected to increase even with only one seasonal spawn. When the entire multistate matrix M was used (with all four quadrants of Figure 1a), the principal eigenvalue (λ = 1.7114) increased by only 0.81% over that of the classic life history alone. Thus, the population was not expected to grow significantly more in a constant environment with the addition of the second seasonal spawn. The lower right quadrant (A₂) of the multistate matrix (Figure 1a) signifies the contribution of the nub life history alone to the bay scallop population. In this case λ = 1.1260, a population growth rate 34% less than either a classic life history alone or one with two seasonal spawns.

TABLE 3.

Left and right eigenvectors (v₁, W₁) for matrices analyzed in the Nantucket bay scallop population indicating reproductive values and stable age distributions, respectively, for classic life history with only a single early spawn (matrix A₁), nub life history with only a late spawn (matrix A₂), and full life history (multistate matrix M) with two spawns per season.

In the “classics only” matrix (A₁), reproductive values (left eigenvector) were lowest for year-2 classic scallops (0.1814) and highest for year-1 adults (2.1794), whereas year-0 juveniles had an intermediate value of 0.5134. In the “nubs only” matrix (A₂), reproductive values were lowest for year-0 juveniles (0.5018) and almost the same for year-1 and year-2 nubs (1.4124 and 1.4260, respectively). When the entire multistate matrix model (M) was used, allowing the two life histories to interact, the reproductive values of year-0 classics and nubs were very different (0.9706 and 0.0155, respectively) and those of year-1 classics and nubs were 4.1525 and 0.0663, respectively. Within year-classes of the multistate model, the closest reproductive values were between those of year-2 nubs (0.5394) and year-2 classics (0.34030) (Table 3).

In the “classics only” matrix (A₁), year-0 juveniles made up by far the largest portion of the population (81%), whereas year-1 classics comprised 19% and year-2 classics only 0.22% of the population at the stable age distribution (right eigenvector). In the “nubs only” matrix (A₂), year-0 juveniles made up the largest portion (72%), whereas year-1 nubs comprised 25% and year-2 nubs contributed 3% to the population at the stable age distribution. When the entire multistate matrix (M) was used, showing the contribution of each life history to the other, year-0 juvenile classics and nubs were found in the greatest quantities (29% and 50%, respectively). Year-1 classics made up 7% of the population, while year-1 nubs were 12%. Year-2 classics and nubs made up 0.08% and 0.9%, respectively, of the population at the stable age distribution (Table 3).

FIGURE 3.

Future projections of total Nantucket bay scallop population based on initial conditions in two survey years (2008 and 2009) and comparing three models (full multistate matrix M, classics only submatrix A₁. nubs only submatrix A₂) for (a) 2008 with low-percent classics—high-percent nubs, and (b) 2009 with high-percent classics—low-percent nubs. Initial conditions were based on September population survey data. Population size is expressed in log₁₀ values. The dashed line is a rough approximation of Nantucket Harbor's carrying capacity (110,000,000 scallops; log₁₀ = 8.05) based on 1980 commercial landings.

The multistate matrix, allowing contribution of nubs to the classic life history and of classics to the nub life history, allowed the population to increase by nearly three orders of magnitude in 10 years, regardless of which of the two initial conditions was used (Figure 3). The classic life history alone, with no second spawn, allowed the population to increase by two orders of magnitude in 10 years, while the nub life history, with only a late spawn, increased less than one order of magnitude regardless of which of the two initial conditions was used. Population growth with two seasonal spawns was greater beginning with the initial conditions for 2008 (low classic-low nub) than with those for 2009 (high classic-low nub). These population projections, however, represent the Nantucket bay scallop's biotic potential for exponential growth, are density independent, and do not take into account the carrying capacity of its estuarine environment. When a rough approximation of carrying capacity was considered (dashed line in Figure 3), the population reached carrying capacity with two seasonal spawns (multistate M model) in 1 to 2 years and in the classics alone model (A₁) with one early spawn in 3 to 4 years (depending on initial conditions), but in the nubs alone model (A₂) with one late spawn carrying capacity was not reached within 10 years.

FIGURE 4.

Elasticity matrices for (a) A₁ classics only submatrix, (b) A₂ nubs only submatrix, and (c) M full multi state matrix. Elasticity values range from zero to one; the larger the value, the more sensitive λ is to changes in that vital rate. Values underlined with double line are highest; those underlined with a single line are next highest in elasticity.

The elasticity matrices indicated that, for the “classics only” (A₁) life history, λ was essentially equally sensitive to juvenile survival (0.4998) (P₁ in Figure 4a) and to the fertility of year-1 adults (0.4993) (F₂ in Figure 4a). In the nubs only” (A₂) life history, λ was more sensitive to the survival of juveniles (0.4607) than to the fertility of year-1 adults (0.3820). In the two matrices, the population growth rate was not especially sensitive to changes in either fertility or survival of year-2 scallops. When both life histories were combined in the multistate matrix, M (Figure 4c), the population growth rate was most sensitive to changes in P₁ (classic juvenile survival, 0.4866) and F₂ (year-1 classic fertility, 0.4785) was a close second. The population growth rate was relatively insensitive to juvenile and year-1 nub survival (0.0131 and 0.0081, respectively) and particularly insensitive to changes in nub fertility (0.0005 and 0.0004 for year-1 and year-2, respectively).

Stochastic Model

When survival rates (P₁) for 2,000 juveniles were randomly selected from a β distribution, the mean of the stochastic population growth rates () of the resulting “classics only” matrices (A₁) was 1.6232. Using the same P₁ values, the mean of the stochastic growth rates of the full multistate matrices (M) increased to 1.7154. The mean of the stochastic growth rates of the full matrix was 5.7% higher than that of classics alone matrix, but nearly the same as the deterministic growth rate (1.7114). When fertility rates (F₂) for 2,000 classic adults were randomly selected from a lognormal distribution, the mean of the stochastic growth rates of the “classics only” matrices was 1.0053 and the mean of the stochastic growth rates () of the multistate matrices was 1.4652, an increase of 45.8%. When juvenile survival (P₁) was allowed to fluctuate together with classic adult fertility (F₂), increased by 58.3% from the classic matrices (0.8762) to the multistate matrices with two spawns (1.3866). All values with their 95% CIs are shown in Table 4. Box plots in Figure 5 compare population growth rates in all matrix models with and without the second spawn.

TABLE 4.

Population growth rates (λ) of classic (A₁) versus full multistate (M) bay scallop life histories using the deterministic model and then stochastic models varying vital rates P₁ (classic juvenile survival) and F₂ (classic year-1 fertility). First P₁ and F₂ varied alone and subsequently both vital rates varied together. Stochastic growth rates calculated from 2,000 iterations of vital rate(s), thus expressed as mean and lower and upper limits to 95% CI (in parentheses). The percent increase from A₁ classic (one spawn) to M multistate (two spawns) is shown in the right column.

FIGURE 5.

Comparison of mean population growth rates (λ) of bay scallop populations in Nantucket Harbor with (M multistate) and without (A₁ classic) the second spawn of (1) the deterministic model, (2) a stochastic model (Stoch) using a beta distribution of juvenile classic survival (P₁), (3) a stochastic model using a lognormal distribution of year-1 classic adult fertility (F₂), and (4) a stochastic model varying P₁ and F₂ together. All vital rates for which growth rates were calculated were drawn randomly from a distribution of 2,000 possible values. Box plots show median values (solid horizontal line), 25% and 75% percentile values (box outline), and 2.5% and 97.5% percentile values (whiskers). The dashed line represents the growth rate at which the population is neither increasing nor declining (λ = 1).

DISCUSSION

When a deterministic matrix population model was employed, the presence of a second seasonal spawn in the Nantucket bay scallop increased the population growth by less than 1% over that of the classic early spawning life history. When stochasticity was added to the model, varying the two vital rates to which population growth rate was most sensitive, the bay scallop population growth increased by 58.3% when the second spawn was added. The lower limits of stochastic estimates (Figure 5) were much greater when a second spawn was included, suggesting there is a much lower frequency of population decrease caused by environmental variability.

Elasticity analysis revealed that the population growth rate was most sensitive to the survival of classic juveniles (P₁), and the fertility of year-1 classic adults (F₂) was a close second. We were surprised that adding a second spawn increased λ when fertility was varied by so much greater a percentage than when juvenile survival was varied alone (45.8% versus 5.7%, respectively) (Table 4). However, since the first age-class in our model was newly settled spat, fertility had to be calculated with a combination of fecundity, fertilization and settlement rates, and the probability of larval survival, almost all of which are subject to wide fluctuations. Malinowski and Whitlatch (1988) found the highest sensitivity to be survival of the youngest age-classes in three commercially important bivalves (Mya arenaria, Mercenaria mercenaria, C. virginica). In general, sensitivity of survivorship values in those species is two orders of magnitude greater than those for fecundity. Doak et al. (1994) found a similar difference in elasticity values between survival and fecundity, but in the case of the desert tortoise, elasticity of survival rates was highest in the largest size-class. Nakaoka (1997) also found the highest sensitivity for survivorship of the largest size-class of the bivalve, Y. notabilis. Greater sensitivities in older age- or size-classes seem to be the norm in such especially long-lived species as the desert tortoise and Y. notabilis.

Because classic juvenile bay scallop survival rate (P₁) and classic age-1 adult fertility (F₂) had the highest calculated elasticities in our study, we chose to vary only those vital rates in performing the stochastic analyses: first P₁ and F₂ each alone and then P₁ and F₂ together. Adding a second spawn increased the population growth rate by 5.7% under conditions of stochastic juvenile survival. There was an increase of 58.3% when the means of the stochastic growth rates were compared between classic and full multistate life histories, using random distributions of both juvenile survival and adult fertility. Although the variability within the estimates of for the two life histories was greater than the difference between the two estimates, there are some meaningful implications. The distribution of population growth rates for each revealed that, while the majority of values were below one (indicating population decline) in the classics only (A₁) stochastic model, no values below one were found in the multistate (M) stochastic model with two seasonal spawns. In all life histories (classic alone, nub alone, and multistate), variability in stochastic growth rates, as shown by CIs, decreased substantially when a second spawn was added (Figure 5).

Juvenile bay scallop survival in this study represented survival from settlement to recruitment into the reproductive population at age 1, which is therefore analogous to recruitment rate in other studies. Marine bivalves, most of which have planktotrophic larval stages, experience considerable recruitment fluctuations (Brousseau 1978a; Strathmann 1985; Mann 1988; Ambrose and Lin 1991; LePennec et al. 2003). Bay scallop populations also have large variations in recruitment; thus, the stochastic estimate of population growth seems to be more relevant. Variable recruitment rates were modeled in stochastic matrix analysis of the subtidal snail, Umbonium costatus, (Noda and Nakao 1996), the bivalve, Y. notabilis, (Nakaoka 1997), and the soft-shell clam, My a arenaria (Ripley 1998; Ripley and Caswell 2006). In the first study (Noda and Nakao 1996), an unstable age distribution indicated sporadic recruitment, and recruitment densities showed yearly fluctuations. The second study (Nakaoka 1997) also noted large interannual recruitment fluctuations and rates that were highly skewed towards smaller values (thus, lognormally distributed). In the last study, Ripley and Caswell (2006) found that the great variability in recruitment rate actually increased the stochastic growth rate of My a arenaria in Barnstable Harbor, Massachusetts, when only those rates are varied. The great persistence of the population may be due to the ability of this soft-shell clam to live long enough to experience at least one good recruitment event in a lifetime.

The deterministic population growth rate (λ) for bay scallops in Nantucket Harbor was 1.7114 with a second seasonal spawn and 1.6979 without. However, it is not uncommon for species with variable reproduction or recruitment to have deterministic estimates of growth rate that are less than or just slightly over one. Population growth rates in the teasel, a biennial plant found at later stages of old-field succession, ranged from 0.275 to 2.605 depending upon location (Werner and Caswell 1977). The population growth rate of G. gemma, a small clam living in Rhode Island sandflats, fluctuated on either side of one because this species does not undergo reproduction each year (Weinberg 1989). Noda and Nakao (1996) calculated an average λ of 0.9268 for the subtidal snail, U. costatum, studied for 8 years when the year with extremely high recruitment (1982) was removed. With that year included, the average population growth rate rose to 1.078. In Nakaoka's (1996) review of population growth rates of 10 studies using matrix population models, four species had λ values less than one: Striped Bass Morone saxatilis, desert tortoise, and the marine mollusks, Y. notabilis and U. costatum (cited above). The other species (jack-in-the pulpit Arisaema triphylllum, savanna grasses Andropogon semiberbis and A. brevifolis, seaweed Ascophyllum nodosum, gorgonian coral Leptogorgia virgulata, and red deer Cervus elephas) had λ values equaling or slightly above one. None of the latter populations are increasing substantially; rather, they are either stable or just holding on.

When the bay scallop population was projected into the future at the growth rate calculated by our deterministic matrix model, adding a second spawn made a difference in the time it took for the population to increase to the carrying capacity. In our stochastic models, population growth rates without a second spawn often fell below one, indicating a declining population (Figure 5). However, adding a second spawn while varying the vital rates of juvenile survival and year-1 adult fertility allowed the population to increase by over 50% from a single-spawn to a two-spawn life history, a much greater increase than that seen in the deterministic model. Although the stochastic growth rate with a two-spawn life history was not as great as the deterministic estimate, we believe that it was more realistic, considering the scallop's variable environment. In the bay scallop, a second spawn could thus be a mechanism to prevent the population from declining to extinction.

Reproductive values of classic early spawned bay scallops in Nantucket followed the typical pattern of being lowest in juveniles, reaching their highest values in adults at first age of reproduction, and then declining with increasing age (Caswell 2001). In contrast, late-spawned nubs, which can spawn twice in their second summer, had their highest reproductive value as age-2 adults, slightly surpassing year-1 nubs, less than 40% of which could spawn once but only late in their first summer (Hall 2014), and both had greater reproductive value than did juveniles. It is possible that juveniles in both cases were underrepresented in our model because we were not able to assess stages earlier than newly settled spat. Malinowski and Whitlatch (1988) found that all three bivalves (Mya arenaria, Mercenaria mercenaria, and C. virginica) showed increasing reproductive values between age 0 and age 2 (age of first reproduction), followed by a slow decline until death.

The stable age distribution represents the constant age structure achieved when the population reaches an asymptote at the given growth rate (λ) (Caswell 2001). This distribution is ergodic (independent of initial conditions). The stable age distribution of Nantucket bay scallops for both classic and nub life histories revealed that juveniles less than 1 year old were found in the highest abundance (81% and 72% for classics and nubs, respectively). Year-1 adults were next in abundance (19% and 25%, respectively), while year-2 adults were lowest (0.2% and 2.9%, respectively). In the multistate matrix, representing the combined contributions of classic and nub life histories, juvenile classics and nubs also made up the largest percentage of the population (30% and 50%, respectively) (Table 3). Relative numbers of scallops in each cohort differed widely from year to year in the September survey. For example, year-1 classics made up 4% and year-1 nubs 39% of surveyed scallops in 2008, while year-1 classics made up 74% and year-1 nubs 10% of surveyed scallops in 2009 (Table 2). Thus the modeled and observed age-class distributions were not the same.

Nantucket Harbor is an unpredictable environment and experiences sudden changes in temperature and salinity as well as episodic blooms of harmful algae and population explosions of predators. Bay scallops are especially sensitive to these changes during larval and juvenile stages when they experience wide fluctuations in recruitment. When only one vital rate (P₁) was varied, representing classic juvenile survival, stochastic population growth rate increased 5.7% from a classics-only life history to one incorporating a second seasonal spawn. When two vital rates were varied together (P₁ and classic adult fertility F₂), the difference in the mean of the population growth rates () between the two life histories was 58.3% (Table 4). The substantial larval survival component of F₂ must have interacted with juvenile survival (P₁) to produce the great increase in population growth rate between life histories with and without a second spawn.

The second spawn of the bay scallop population in Nantucket Harbor appears to illustrate the life history strategy of bet-hedging (Slatkin 1974). In bet-hedging, there is a tradeoff between current reproductive effort (and fecundity) and future survivorship (Caswell 1980; Goodman 1982; Strathmann 1985; MacDonald et al. 1987). Some examples of bet-hedging are spreading reproduction over multiple years in biennial plants (Klinkhamer and DeJong 1983), prolonged dormancy in annual desert plants (Philippi 1993), delayed hatching in anostracans living in vernal pools (Simovich and Hathaway 1997; Philippi et al. 2001), hatching asynchrony in nesting birds (Laaksonen 2004), and extended hatching of larval king crabs (Stevens 2002). Bet-hedging appears to be particularly advantageous in stochastic environments, where adaptations to unpredictable environmental conditions ensure that at least some offspring will find suitable conditions in which to survive (Menu et al. 2000; Krug 2009). Bet-hedging has been observed in several scallop species. “Dribble spawning” (partial spawning over a long period of time) in the sea scallop ensures that at least some larvae will survive in an uncertain environment (Langton et al. 1987). The Australian saucer scallop Amusium balloti also employs partial spawning to provide multiple opportunities for larval development and recruitment success (Joll and Caputi 1995). Tettelbach et al. (2001) hypothesized that the ability of late-spawned New York bay scallops to survive and spawn a second year suggests that this is important to ensure the population persists when early recruitment has failed. Bishop et al. (2005) stated that the adoption of an iteroparous (more than one spawning event) life style by North Carolina bay scallops is an example of bethedging. In that case, it is the early spawn that persists in spite of a large fall spawn because of the latter's risk of recruitment failure from autumn events such as harmful algal blooms and hurricanes. Multiple spawnings in Florida bay scallops throughout the year could allow for a greater reproductive output over their lifetimes, another example of bet-hedging (Geiger et al. 2010).

Bimodal spawning in northern bay scallop populations thus appears to be a strategy for coping in a stochastic environment. Fluctuations in recruitment, combined with a short lifespan, have caused the bay scallop to be particularly vulnerable to environmental change. The progeny of the late spawn may be exposed to a different set of environmental conditions than that of the early spawn and may experience greater growth and survival. By spreading out reproductive effort over more than one spawning event in a season, the northern bay scallop may be capable of increasing its chances for successful recruitment. This life cycle diversity also gives the Nantucket bay scallop population its resilience, or the ability to recover quickly from environmental perturbations. The importance of the second spawn to the persistence of northern bay scallop populations has significant management implications. The greater the survival of the second spawn's progeny (nub scallops), the more they can contribute to the population in their second spawning season. If recreational and commercial harvest of first-year nubs is again prohibited, a major source of their mortality could be eliminated. With that accomplished, a successful second spawn would be able to buffer the effects of the bay scallop's highly variable recruitment rate and prevent possible population collapse.