Reconstruction of environmental histories to investigate patterns of larval radiated shanny (Ulvaria subbifurcata) growth and selective survival in a large bay of Newfoundland

Abstract

We used otolith microstructure analysis to reconstruct the growth histories of larval radiated shanny (Ulvaria subbifurcata) collected over a 2-week period in Trinity Bay, Newfoundland. A dynamic 3-dimensional, eddy-resolving circulation model of the region provided larval drift patterns, which were combined with measurements of temperature and zooplankton abundance to assess the environmental history of the larvae. The abundance of juvenile and adult capelin (Mallotus villosus), the dominant planktivorous fish in this area, was monitored using five hydroacoustic surveys. The goal was to determine whether environmental histories are helpful in explaining spatial and temporal differences in larval shanny growth, measured as cumulative distribution functions (CDF) of growth rates. We found evidence for a selective loss of slower growing individuals and recognized considerable spatial differences in the CDF of larval growth rates. Consistent patterns in capelin abundance suggested that faster growing survivors, sampled at the end of the 2-week period, developed in areas of low predator densities. A dome-shaped relationship between temperature and larval growth was observed, explaining a significant but small amount of the overall variability (14%). Effects of experienced prey concentrations on larval growth rates could not be demonstrated.

Introduction

To develop more useful predictions of recruitment, it is critical to understand the processes that determine growth and survival of larval fish in the sea. To acknowledge the complexity of underlying mechanisms, it has been emphasized that process studies need to be based on individuals rather than on averaged populations (Rice et al., 1987; Miller et al., 1988). Each individual larva develops its own set of traits (e.g. growth rate) as a result of both intrinsic and external factors, which may confer a different probability of survival under variable environmental conditions. If we study a characteristic trait in surviving larvae and then try to reconstruct the environmental history as well as the prior distribution of these traits in the population, it may be possible to gain insight into processes influencing early survival (Letcher et al., 1996).

Daily growth rates may provide a trait appropriate to characterize surviving larvae, for they reflect the growth potential of an individual (intrinsic factors) as well as environmental influences, such as temperature or feeding conditions (Parma and Deriso, 1990). Growth rates are also thought to be indicative of a larva's vulnerability to predation (Rilling and Houde, 1999). In general, it has been hypothesized that individuals with the highest growth rates are more likely to survive, because larvae of a bigger size-at-age usually have better foraging and escape abilities (Miller et al., 1988; Folkvord et al., 1997). However, several studies suggested that faster growth may lead to a higher probability of encountering planktivorous fish (Litvak and Leggett, 1992; Pepin et al., 1992), which are seen by some as the most important source of mortality in older larvae (Hunter, 1981; Paradis and Pepin, 2001).

The reconstruction of a larva's growth history can be done most effectively by means of otolith microstructure analysis. Because increment widths are directly correlated with daily somatic growth rates (Moksness and Wespestad, 1989), a back-calculation of previous lengths-at-age is not necessary to recognize the effects of the environment on individual growth rates or their distribution within a population (Gallego et al., 1996). However, it is critical to acknowledge that the widths of subsequent increments are not independent but correlated with each other (Mosegaard et al., 1988). Pepin et al. (2001) estimated that, due to serial correlation, larval growth rates will not reflect the potential influence of the local environment for at least 3 days. Therefore, once we sample a larva, we need to reconstruct its environmental history, rather than study the processes based on the conditions on the day of capture. Furthermore, because advection greatly influences the environmental experience, a larva's drift history must be reconstructed as well, at least in a general sense. This can be achieved by coupling survey observations with a dynamic circulation model, which allows the schemes of larval transport or dispersal to be inferred by tracking computer-generated particles through the model domain (Gallego et al., 1999; Voss et al., 1999; Hinrichsen et al., 2001). Using these approaches, Pepin et al. (2003) studied larval growth patterns in Conception Bay, Newfoundland, and found indications of selective loss of faster growing individuals due to predation by capelin (Mallotus villosus), the dominant planktivorous fish in this area. However, the authors suspected that their representation of the predator environment was probably inadequate, and suggested that patterns in predator abundance may be more readily detected in larger study sites, where relative movements of pelagic fish would be of less importance.

Here, we study patterns of growth and survival of larval radiated shanny, Ulvaria subbifurcata (Storer, 1839), a small blennoid species spawning demersal eggs in spring and early summer (Green et al., 1987). Shanny larvae hatch after 35–40 days at average lengths of 6.5 mm TL, ascend to and stay in the surface layer until they settle at the end of the summer at average lengths of 18.4 mm TL (LeDrew and Green, 1975).

During three surveys carried out over a 2-week period in July 2000, larval fish and environmental data were collected in Trinity Bay, Newfoundland, a study site three times bigger than the adjacent Conception Bay. Juvenile and adult capelin abundance was monitored using five hydroacoustic surveys. Otolith microstructure analysis provided larval growth histories while larval drift was simulated using a fine-resolution, wind-driven, 3-D eddy-resolving circulation model.

Materials and methods

Study site, survey design, and data collection

Three surveys of Trinity Bay (48°N, 53.5°W) were conducted on board the research vessel CSS Wilfred Templeman between 18 and 30 July 2000. Trinity Bay is approximately 100 km in length and 30 km in width, and is one of the major bays along the Atlantic coast of Newfoundland (Figure 1). In its central area, a deep trench with a maximum depth of 630 m runs parallel to the coastline, while the sill at the bay's mouth has a maximum depth of 240 m. The bay is influenced by the inshore branch of the cold Labrador Current, which is strongest in spring (Davidson and de Young, 1995). Like other bays in Newfoundland, stratification and circulation are highly responsive to wind forcing (Templeman, 1966; Yao, 1986). In summer, southwesterly offshore winds prevail and induce upwelling. Conversely, periods of onshore winds are likely to suppress upwelling and intensify surface heating.

Figure 1

Map of the study area and the 29 sampling locations (crosses), laying on nine transects (numbers 1–9) of three to four stations (TB, Trinity Bay; CB, Conception Bay; BB, Bonavista Bay; Meteorological stations: ♢, St. John's airport; ♦, Bonavista).

Each survey covered the bay with a grid of 29 sampling stations along nine transects (three to four stations each) with an approximate distance between stations of 8 km (Figure 1). Each survey took 35–44 h to complete, and intervals between surveys were 5 and 4 days, respectively.

Larval fish were collected in single oblique tows of a 4 m² Tucker trawl deployed from the surface to 40 m depth. As previously shown (De Young et al., 1994; Pepin et al., 1995), this sampling layer was likely to contain more than 95% of the ichthyoplankton. The net had sections of 1000, 570, and 333 μm mesh nitex (Pepin and Shears, 1997) and was towed for 10–15 min at 1 m s⁻¹. The volume of water filtered during the tow was estimated with two General Oceanics flow meters, mounted at the mouth of the net. This Tucker trawl configuration has been shown to reduce sampling bias and avoidance for larvae up to 20–25 mm long (Pepin and Shears, 1997). On deck, samples were immediately preserved in 95% ethanol. Samples were drained after 1 week to renew the preservative. All fish larvae were subsequently separated from other planktonic organisms and identified to species or to the lowest taxonomic level possible. Up to 200 specimens per taxon and station were measured (standard length, nearest millimeter) using a dissecting microscope and a gridded background.

At each station, a 0.5 m diameter plankton net (70 μm mesh) was deployed to 40 m depth and hauled vertically at 1 m s⁻¹ to the surface. The zooplankton samples were preserved in 2% buffered formaldehyde. In the laboratory, samples were rinsed through 500 and 65 μm sieves to separate the microzooplankton fraction, which was then split six to eight times using a Motoda splitter and double-filtered seawater, until dilution was appropriate for recording particle numbers with a Coulter Multisizer II^® counter. Particles were discriminated at size intervals of 1.57 μm. Counts of all size classes >65 μm were subsequently summed to yield the overall abundance of microzooplankton per station (particles l⁻¹).

At each station, a Seabird-25 CTD (sampling rate of 8 Hz) was lowered at 1 m s⁻¹ to record temperature, salinity and fluorescence profiles of the water column down to a maximum depth of 500 or 5 m above the bottom. A value representative for the ichthyoplankton layer was obtained by averaging measurements of the top 30 m because variability in temperature and salinity in the cold intermediate water layer below this depth was negligible. Because we do not know the vertical position of larval radiated shanny accurately, the averaged temperature provides an index of the thermal environment rather than an accurate reflection of the true temperature experienced by the larvae. This may introduce error in underlying relationships as well as limit our ability to make inferences about some factors that are influencing larval growth rates.

Abundance of juvenile and adult capelin was estimated acoustically along transects while the ship traveled at speeds of 5–10 knots. The vessel was equipped with a calibrated hull-mounted split-beam 38 kHz transducer and an EK500 echosounder system. High-resolution (10 cm vertical bins, one ping per second) backscatter volume (Sv) measurements were acquired at a threshold of −100 dB to a depth of 200 m. Echogram files were subsequently edited and integrated using an Sv threshold of −85 dB. Echogram mark characteristics such as shape, density and target strength distribution were used to differentiate between fish and crustaceans (Simard and Lavoie, 1999). Capelin were distinguished from other fish backscatter using echogram characteristics combined with species composition and biological characteristics. Biological data originated from 26 oblique IYGPT trawls (Koeller et al., 1986) at four different sites in Trinity Bay conducted during the course of the study. From each IYGPT sample, up to 200 capelin were selected at random and measured for total length to the nearest millimeter.

Backscatter values (Sv) originating from capelin were integrated over the whole water column in 10 m depth intervals and 100 m horizontal bins to produce estimates of the back-scatter area (Sa) in m² m⁻² (MacLennan and Simmonds, 1992). Values originating from the top 10 m were not included, because capelin above or near the transducer are not accurately detected. The back-scatter area (Sa) was then scaled to numbers of individuals m⁻² by applying the length-frequency distribution of capelin, sampled during IYGPT trawls, and the target strength relationship 20×log L−73.1, where L is length in cm (Rose, 1998). We further estimated median capelin density per station and day, including only observations within 4 km of each ichthyoplankton sampling location. To ascertain the stability of capelin distribution patterns during the entire period, echo data were complemented by two additional dedicated hydroacoustic surveys along the same transects and at the same speeds during the interval in-between the sampling surveys.

Otolith microstructure analysis

Up to 15 radiated shanny larvae were chosen at random from each sample for otolith microstructure analysis. Each larva was assigned to a unique identification number and measured for total and standard length (to the nearest 0.1 mm) using an Optimas^® image analysis system connected to a dissecting microscope. Both sagittal otoliths of 359 larvae were extracted and each otolith was mounted on a drop of Crystal Bond^® thermoplastic cement on a microscope slide. Choosing randomly the left or right sagitta, unless one side proved impossible to read, the otolith was then ground to near mid-plane using a 0.3 μm lapping film. All otoliths were read under 500× magnification with an Olympus BH-2 compound microscope connected to the Optimas^® image analysis system. Individual increment widths were measured along the longest otolith axis.

Because of a greater potential for measurement error for the earliest narrow increments (Pepin et al., 2001), all otoliths from individuals younger than 12 days were read twice. For older larvae, measurements were repeated in every 10th individual. In both cases, otoliths were rejected if counts differed by more than 10%, which led to 12 exclusions. Care was taken that otolith readings were not ordered by station number, survey, or larval length. Each increment was considered to represent 1 day. Therefore, the last increment was excluded from analyses, because its formation might not have been completed on the day of sampling. Because mean widths and variances increase with age (Pepin et al., 2001), patterns in increment widths were analyzed using data standardized to zero mean and unit deviation

(1)

where x_ij is the increment width (μm) of the individual i at age j, and x_j and s_j are the mean and standard deviation of increment widths at age j, respectively, for all specimens observed. To avoid confusion with true somatic or otolith growth rates, we will hereafter refer to standardized increment widths as ‘relative growth rates’. Otolith growth in radiated shanny is closely linked to somatic growth (Pepin et al., 2001). Although there is some evidence that slower growing individuals have slightly larger otoliths at a given size (but older age) than individuals with high growth rates, daily increment widths provide a good approximation of daily growth rates. In the absence of evidence that relative growth rates are normally or log-normally distributed, there is no guarantee that the fundamental assumptions for using parametric analytical methods are satisfied. We therefore derived age-dependent, cumulative distribution functions (CDF) of standardized increment widths by using local non-parametric density estimations (Davison and Hinkley, 1997; Pepin et al., 1999; Evans, 2000). This method does not estimate a single distribution but a separate, locally weighted CDF for each possible value of x (age), assuming that observations nearest to the target x were most relevant for the distribution at x. The advantage of using CDF is that it allows to contrast changes in the overall distribution of relative growth rates, rather than dealing exclusively with the mean as is common for most statistical approaches. The 10, 50, and 90% probabilities were used to describe shape and changes of the CDF of increment widths. The difference between the 90 and 10% probabilities served as a measure of the variability of relative growth rates. To decide whether two medians differed significantly, the 95% confidence intervals around the median were estimated from 500 randomizations of the data. If confidence intervals did not overlap, the two medians were considered to be significantly different.

To study patterns of selective mortality, we contrasted CDF of relative growth rates on 18 and 24 July among larvae, which were sampled during the three subsequent surveys. Because of the degree of autocorrelation in adjacent increment widths (Pepin et al., 2001), we included not only the most recent but the last four increments formed prior to the day of comparison. For example, to analyze growth patterns of larvae present on 18 July, we compared increment widths 6–9 of all 10-day-old larvae caught during the first survey with those of all 16-day-old larvae caught during the second survey and of 21-day-old larvae caught during the third survey.

To determine to what extent relative growth rates show a functional response to the environmental conditions experienced by the larvae, a General Linear Model was used relating all increment widths formed during the study period (for which we have environmental information) to two independent variables: a 3-day average of both temperature and microzooplankton concentration experienced on the day of increment formation and the 2 preceding days as inferred from larval drift paths predicted by the circulation model.

Drift projections

Circulation patterns in Trinity Bay were simulated using Davidson et al.'s (2001) 3-D eddy resolving CANDIE model. The model solves 3-D non-linear Navier–Stockes equations (i.e. the equations of x, y, and z momentum as well as the density and continuity equations) on an f-plane using three standard oceanographic approximations: hydrostatic and rigid lid (Gill, 1982), and the Boussinesq approximation (Spiegel and Veronis, 1960). The equations are finite differenced on a 3-D, Arakawa C-grid with 1 km horizontal and 10 m vertical resolution, which allows the resolution of the dominant scale of eddies in the system under study. The model domain consisted of a realistic coastline geometry and bottom topography of Trinity Bay and the two adjacent bays (Bonavista Bay, Conception Bay). On the lateral land boundaries, a free-slip condition was applied to avoid a flux of momentum or density through the water–land interface, whereas a Neumann boundary condition (Greatbatch and Otterson, 1991) was used for all three open boundaries. Because no external influence on the model domain such as the inshore branch of the Labrador Current was applied (i.e. passive boundaries), the model represents wind-forced circulation in the presence of realistic vertical stratification and topography.

Winds were assumed to be spatially homogeneous over Trinity Bay (de Young et al., 1993). Because hourly measurements by Environment Canada at St. John's and Bonavista airports (Figure 1) were sufficiently similar, only the St. John's data were eventually used in the model (Pepin et al., 2002). All model runs were initialized at rest (i.e. no motion) using a horizontally uniform stratification based on the long-term mean (1957–1997) for July at Station 27 (5 km east of St. John's). Wind stress derived from velocities by the quadratic formulation of Large and Pond (1981) was smoothly introduced over 2 days using a hyperbolic tangent ramping function. The time measurement began on day 190 after winds reached 50% strength, and ended after 35 days of simulation on day 225.

Larval drift patterns were inferred from the paths of computer-generated particles tracked smoothly through a finite resolution velocity field of the model domain (Davidson and de Young, 1995). To allow for dispersion, the drifters were given a random-walk component to their displacement at every time step equivalent to a diffusion coefficient of 10 m² s⁻¹. This value is lower than estimated by de Young and Sanderson (1995) but is appropriate for a simulation with a circulation model with high spatial resolution as used here. Because we were interested primarily in the horizontal dynamics of particle movement and potential patterns in larval vertical migration were unknown, particles were constrained to one of two seeding layers (0–10 and 10–20 m), which were modeled separately, in contrast to water masses which were allowed to flow among layers (Davidson et al., 2001). Throughout this study, we will focus on the second layer (10–20 m), which is less responsive to wind forcing (Pepin et al., 2003) and may therefore provide a better approximation of the larval drift in the ichthyoplankton layer (0–40 m). The patterns for the two layers were similar but the variability in drift trajectories was slightly greater in the upper layer. However, the overall conclusions are not affected by our choice.

To assess the origin of larvae associated with a given sampling station at a given time, we used backward projections of particle tracks. In the first simulation over 6 days, particles were uniformly seeded on 19 July (500 m apart) and assigned to sampling stations on 25 July (the end of the second survey) only if they were located within 3 km of the station location within the model domain. A second simulation lasting 11 days started on 19 July to assign the source of larvae collected on 30 July (the end of the third survey).

Future drift of larvae sampled at a particular station was simulated by forward projections. A 12-day simulation was run by seeding 1000 particles within a radius of 1 km around every sampling station occupied by larvae on 19 July and a similar one on 24 July. In both cases, particle positions were recorded at 24 h intervals until 30 July.

Results

Environmental conditions

Between 15 and 30 July 2000, winds over Trinity Bay blew predominantly from southern or southwestern directions with speeds ranging from 5 to 25 km h⁻¹. During the first 7 days, strong offshore winds prevailed (average speed 19.8 km h⁻¹) and caused upwelling of cold, intermediate water masses along the bay's western shore (Figure 2a). This was reflected in the observations during the first survey (18 July), when average temperature in the 5–30 m layer at western shore stations was substantially lower (3.0°C) than in the rest of the bay (5.2°C). Offshore winds ceased abruptly on 22 July and remained moderate until the end of the study period (average speed=9.7 km h⁻¹). During the second survey (24 July), the upwelling pattern had vanished, but cold, intermediate water masses (3.0°C) were now observed in the surface layer of the inner half of Trinity Bay (to approximately 48.1°N, Figure 2b). The outer half of the bay was covered by warm water masses (5.2°C), which appeared to propagate south into the bay. During the third survey (July 29/30), an overall temperature increase (Figure 2c), along with a depression of all isotherms, indicated that the warm water masses had fully displaced the cold water body from the surface.

Figure 2

Environmental conditions in the surface layer of Trinity Bay during three subsequent surveys in July 2000 (spatial interpolation on a 1×1 km grid by kriging). (a–c) Average temperature (°C; 5–30 m); (d–f) average microzooplankton concentration (particles l⁻¹; 0–40 m).

Abundance of microzooplankton was lowest during the first survey with average concentrations of 41 particles l⁻¹ (Figure 2d). During the second survey, concentrations had increased on the western side of the bay and at the mouth. With 79 particles l⁻¹ on 29 July, average microzooplankton concentration was significantly higher (Welch's one-way ANOVA: F_2,49=23.3, p<0.001) than during both previous surveys (Figure 2f). The observed differences were associated with a weak tendency for decreasing particle size over the 2-week period.

Hydroacoustic estimates of capelin abundance (10–200 m) ranged from zero to 128 fish m⁻² with median densities of 0.04 fish m⁻², without major changes in the distribution between surveys (Figure 3). Data characteristically were highly variable between consecutive 100 m echo bins. On a broader spatial scale, densities appeared to be consistently higher along the western shore, lower along the eastern shore, and particularly low in central areas of the inner bay, where the deep trench is located. Median capelin abundance, calculated for a 4 km radius around stations, substantiated this impression and revealed significantly higher densities at stations along the western shore than at stations in central or eastern parts (Welch's one-way ANOVA: F_2,197=30.9, p<0.001).

Figure 3

Hydroacoustic estimates of capelin density during five surveys in Trinity Bay in July 2000 (spatial interpolation on a 1×1 km grid by kriging; contours, number of fish m⁻²; thin black lines, survey tracks.

Larval distribution

The ichthyoplankton community was dominated by capelin and radiated shanny larvae, which together comprised about 50% of all fish larvae. The spatial distribution of shanny larvae was very heterogeneous, with concentrations >10 larvae 1000 m⁻³ occurring in dense clusters limited to the outer areas of the bay (Figure 4). In inner areas (to approximately 47.9°N), abundance was consistently low during the entire period. The distribution during the second survey was suggestive of two separate larval cores, one with higher concentrations at the western shore and one with lower concentrations at the mouth of the bay. In the western core, the proportion of small larvae (<7.5 mm SL) was notably higher compared with the northern core, whereas large larvae (>11.5 mm SL) were more abundant in the northern than in the western core. Overall, the proportion of larvae in the larger size classes increased slightly over the period, yet relative abundance of larvae <8.5 mm SL remained high, indicating continuous production and emergence.

Figure 4

Spatial distribution (larvae per 1000 m³) and relative length frequencies of radiated shanny larvae encountered. The second survey distribution suggested a western (diamonds, black bars) and a northern core (triangles, white bars) of high abundance.

A population growth rate of 0.32 mm d⁻¹ (derived from age–length data on the day of capture) was used to estimate larval mortality between successive surveys. Larvae were assumed to grow 2 mm SL in the 5 or 6 days between surveys, given that measurements were to the nearest millimeter. Negative losses, indicating production, were found for the smallest larvae sampled (5 mm SL). Mortality rates were significantly greater than zero for larvae of 7–11 mm SL (F_1,26=15.1, P<0.01; Figure 5). However, differences among size classes were never significant (F_8,26=1.23, P>0.2).

Figure 5

Daily mortality of shanny larvae (error bars: 95% confidence intervals around the mean) in relation to their length during the first survey.

Circulation patterns and particle drift

Model outputs were indicative of two major periods in the circulation pattern coinciding with the main shifts in strength and direction of winds during the 2 weeks. Upwelling, induced by strong southwesterlies, characterized the first period with medium, northeasterly currents on the western shore and stronger, outflowing currents close to the eastern shore (Figure 6a–c). Currents at the mouth of the bay did not flow directly out of the region but instead flowed across the mouth, possibly allowing larvae to remain within the region. High-density water masses propagated eastward, covering about two-thirds of the bay on 21 July, the day prior to the sudden relaxation of offshore winds. During the second period, characterized by moderate and variable winds, high-density water masses gradually disappeared from the surface layer (Figure 6d–f). Also currents became weaker, with the exception of the western side of the bay, where a narrow near-shore jet (∼43 cm s⁻¹) developed on 22 July, rapidly dragging lighter and warmer water masses southward (Figure 6d). On the eastern shore, currents reversed and flowed into the bay for the rest of the period (Figure 6d–f). Overall, model predictions corresponded well to CTD measurements but were deviant for the inner part of the bay during the second survey. In this area, cold and high salinity waters still dominated the surface layer on 24 July (Figure 2b), while the model indicated a faster replacement by low density waters (Figure 6).

Figure 6

Selected images from the 3D, eddy-resolving CANDIE model, showing the major circulation patterns in 10–20 m depth (black arrows, speed and direction of the predicted currents; background shading, density anomaly with respect to the initial stratification).

Backward projections of particle trajectories suggested a different environmental history for the two cores where high shanny densities were encountered during the second survey. Figure 7a depicts the drifter positions on the seeding day, 19 July, and the locations (3 km around stations) where particles were found on 25 July. Red particles representing northern core stations originated from central areas of the bay's mouth and were associated with substantially higher average temperatures (6.2°C, Figure 7d) than blue particles representing western core stations (3.7°C, Figure 7d). The predicted origin of larvae in the latter was the inner half of the western shore: a region in which our sampling had indicated very low densities on 18 July (Figure 4). The backward projection of particles seeded on 19 July and collected on 30 July (Figure 7b) indicated that larvae comprising the only cluster of high abundance during the third survey probably originated from central outer areas of the bay and experienced intermediate temperatures ranging between 4.9 and 5.6°C (Figure 7d). The starting positions of particles indicated that only a few larvae would have drifted north from inner areas of the bay, and these would not have originated from the area of coastal upwelling. Instead, they would have drifted from inner central portion of the area.

Figure 7

Particle drift projections in model space to infer origin and fate of shanny larvae encountered. (a) Six-day backward run for particles corresponding to the two cores of high abundance stations during the second survey (circles); (b) backward run for particles corresponding to the region of high abundance observed during the third survey (circle; particle positions are plotted together for all 11 days of the simulation); (c) 6-day forward projection of particles seeded on 24 July, 1 km around western and northern core stations (circles); and (d) thermal history of larvae from the western (blue) and northern (red) cores as well as from the high abundance cluster observed during the third survey (brown) based on average particle drift tracks (black lines in a–c) and linearly interpolated (5–30 m) temperature fields. Solid lines, backward projections; dashed lines, forward projections; no environmental projections were made when average tracks left the area of environmental interpolation).

Forward simulations (24–29 July) revealed that larvae found in the western core during the second survey (blue particles, Figure 7c) were likely to drift north in a confined area along the western shore. Figure 7c illustrates that this drift pattern would not have transported notable portions of this group of larvae to the central part of the outer bay where high abundances were observed during the third survey (Figure 7b). In contrast, particles seeded in northern core stations (red particles, Figure 7c) moved slowly towards the central outer bay and into the area where high concentrations were observed on 29 July. These results suggest that shanny larvae sampled during the third survey should be more similar in their characteristics to those sampled during the second survey in the northern core than in the western core.

The drift simulations indicate that advective losses from the bay were small (<2%). Although particles associated with stations on the eastern side of Trinity Bay could have moved out of the system for part of the period, the circulation pattern would have returned them to the area by the end of the period (Figure 6). However, it is important to be cautious in interpreting advective losses because the model does not have a specific formulation for the variability in the inshore arm of the Labrador Current, which is likely to influence the advection of plankton into and out of the bay.

Because the forward drift projections suggest that the western and northern cores represent distinct aggregations of larvae, we used the average drift trajectory over the 2-week period to assess the cumulative distribution of juvenile and adult capelin densities that both cores would have encountered. We used the non-parametric local density estimator to generate a CDF of predator densities by combining drift trajectories with the hydroacoustic data with a bandwidth of 4 km. This analysis clearly indicates that larvae from the northern core would have encountered significantly lower capelin densities than larvae from the western core (medians lower by a factor of 10; Figure 8).

Figure 8

CDF of capelin densities encountered by larvae from the western and northern cores of abundance during the drift simulation from the second to the third survey using the non-parametric local density estimator based on a bandwidth of 4 km (see Figure 7 for details; drift trajectories were based on the mean daily position of particles released from each station within the core areas).

Growth patterns

Median relative growth rates of larvae present on 18 July, as derived from back-calculations, were lowest during the first survey, slightly higher during the second, and significantly higher in the CDF of survivors sampled during the third (Figs. 9a–c, 10a). Differences were most pronounced for larvae that were estimated to be 10–20 days old on 18 July. The difference in the distribution of relative growth rates between the second and third survey was greater for larvae present on 24 July (Figure 9d–e) than for those present on 18 July. On 29 July, larvae that had been 10–30 days old during the second survey had significantly higher growth rates than individuals in these age classes present during the second survey (Figure 10b). When contrasted with the first and second survey CDF, the lower percentile of the third survey CDF was notably shifted upwards between ages 10 and 20, indicating a selective loss of slower growing individuals throughout the study period (Figure 9a–c). The difference in the lower percentile between the second and third survey was substantially greater for 15–25 day old larvae present on 24 July (Figure 9d–e) than between the second and third survey on 18 July (Figure 9b,c). A broadening of CDF from the first to second survey was indicated by an upward shift in the upper percentiles of survivors present during the latter and continued to the third survey.

Figure 9

Survey comparison of the CDF (solid lines: 10th, 50th, and 90th percentiles, estimated using a local non-parametric density estimator with a bandwidth of 3 days) of standardized increment widths (SIW), formed on 14–17 July in relation to larval age on 18 July (a–c), and of those formed on 20–23 July in relation to larval age on 24 July (d–e).

Figure 10

Median of the CDF of standardized increment widths (SIW) and its 95% confidence intervals (estimated from 500 randomizations of the data) by sampling date in relation to larval age on (a) 18 July and (b) 24 July (based on the increments formed on the 4 days prior to these dates).

We re-analyzed differences in the distribution of relative growth rates for the two distinct aggregations separately, because they represented apparently separate drift and environmental histories (Figure 7a,d). Larvae from the western core had significantly lower median relative growth rates than those from the northern core (Figure 11b). The lower and upper percentiles of the CDF in the northern core were notably shifted up relative to those observed in the western core. The CDF of larvae from the northern core (Figure 11a) were more similar to that observed during the third survey for individuals present on 24 July (Figure 9e). This suggests that individuals present during the final survey represent a spatial subset of the overall population present during the second survey. The variability of relative growth rates was slightly higher in the northern versus western core for 5–22 day old larvae (Figure 11a).

Figure 11

(a) CDF of standardized increment widths (SIW) formed on 20–23 July in relation to larval age on 24 July for the western (gray circles) and northern (black triangles) core of high larval abundance observed during the second survey (see map; solid lines: 10th and 90th percentiles of the CDF, estimated using a local non-parametric density estimator with a bandwidth of 3 days), and (b) medians and their 95% confidence intervals (estimated from 500 randomizations of the data) of these in relation to larval age on 24 July (gray lines, western core; black lines, northern core).

For increments formed during the study period (n=1382), the General Linear Model showed a significant, dome-shaped relationship between relative growth rates and temperature (Figure 12). Adding microzooplankton abundance as a second independent variable did not result in a better fit of the model (Table 1). Although the variance explained by the model (14%) is statistically significant, it represents a small fraction of the overall variability in growth of larval radiated shanny.

Figure 12

Standardized increment width in relation to a 3-day average of temperature (°C) and microzooplankton abundance (particles l⁻¹) experienced on the day of formation and the 2 preceding days, based on larval drift paths inferred from average particles trajectories.

Table 1

Results of a General Linear Model analysis, relating standardized increment width (dependent variable) to a 3-day average of temperature (T) and microzooplankton abundance (MZP), experienced on the day of formation and the 2 preceding days, based on larval drift paths inferred from average particles trajectories.

Source	DF	SS	MS	F-value	P	Par estimate	SE
Model	4	230.6	57.7	57.6	<0.0001
Error	1378	1379.0	1.0
Type III
T	1	73.7	73.7	73.6	<0.0001	3.83	0.45
T2	1	115.6	115.6	115.5	<0.0001	−0.33	0.03
MZP	1	2.1	2.1	2.1	0.15	0.03	0.02
MZP×T	1	1.4	1.4	1.4	0.23	−0.005	0.004
Intercept						−10.26	1.47
Corr. total	1382	1609.6

R²=0.14.

Discussion

By combining otolith microstructure analysis with larval drift patterns and environmental measurements, we gained insight into processes that influenced growth and survival of radiated shanny larvae in eastern Newfoundland. The reconstructed thermal environments of the larvae were useful to interpret spatial differences in growth rates within the area, suggesting a dome-shaped relationship between temperature and growth. We also found evidence for a selective loss of slower growing individuals during the 2-week period.

Ambient temperature is recognized to be one of the major environmental factors influencing individual growth rates (Pepin, 1991; Heath, 1992), which is consistent with our findings showing a close correspondence between spatial patterns in surface temperature and shanny growth rates. This was most pronounced during the second survey, when the reconstructed thermal histories indicated considerable differences between the two cores of larval abundance. Given that growth rates can accelerate by 10% per degree centigrade (Pepin, 1991; Heath, 1992), we conclude that the magnitude of temperature difference had the potential to lead to the observed growth differences during the second survey. More importantly, the different thermal histories experienced by the two groups of larvae could not have been inferred from observations on the sampling day alone. Only the reconstructed temperature history provided a clear perspective of the differences in growth rates. This appears to confirm the need to assess environmental histories in order to understand the processes influencing larval populations (Gallego et al., 1996; Pepin et al., 2001).

The dome-shaped relationship between experienced temperature and relative growth rates is consistent with results presented by Campana and Hurley (1989) for cod and haddock larvae in the Gulf of Maine. However, in the absence of information on larval vertical distribution, the estimated temperatures should not be expected to reflect the true environment experienced by an individual larva, but to represent an index informative of the mesoscale temperature patterns conducive for a cohort's growth. The relationship suggests that the typical increase in summer surface temperature in coastal embayments can result in sub-optimal growing conditions for larvae of cold-adapted fish species. The addition of data on reconstructed prey abundance did not enhance our ability to explain variability in growth rates, a result similar to other field and theoretical studies unable to demonstrate effects of food levels on growth rates (Gallego et al., 1996, 1999; Letcher et al., 1996; Pepin et al., 2003). In this case, it should be noted that the automatic particle counter allowed only a gross estimate of microzooplankton abundance, whereas counting and staging of individual copepod nauplii might have been more effective to evaluate a potential relationship between larval growth and prey density. However, fluctuations in prey abundance followed the spatial scales of the physical environment, which is consistent with other studies in this area (Frank and Leggett, 1982). Therefore, overall patterns of zooplankton abundance should have been captured reasonably well by our sampling programme. Failure to detect a positive correlation between larval growth rates and prey abundance may therefore suggest that events at the scale of the individual are more important in determining its growth and survival potential (Pepin et al., 2003). Unfortunately, environmental variability on such small spatial and temporal scales has yet to be adequately described.

Despite the evidence that thermal history significantly influenced larval growth, the signal was surprisingly weak and explained only 14% of the overall variability. Intrinsic factors responsible for a larva's growth potential may therefore be more important in determining the variation in shanny growth rates. The large unexplained variability could also have resulted from partly unrepresentative environmental parameter values assigned to individual increment widths. Because spatial information was restricted to location of capture on the sampling day, temperature and zooplankton variables were estimated from the average position of all particles per station tracked backwards through the model domain. However, with every time-step the spatial and environmental uncertainty increased as a consequence of particle dispersion that was explicitly included in the model.

The selective loss of slower growing individuals is consistent with expectations based on size-dependent food web dynamics (Peterson and Wroblewski, 1984) and corroborates the findings of previous field and mesocosm studies (Meekan and Fortier, 1996; Folkvord et al., 1997). However, we were unable to detect significant size-dependent differences in estimated mortality rates. The selective removal of slower growth rates is in marked contrast to Pepin et al.'s (2003) study in Conception Bay where they found evidence for selective loss of faster growing individuals, which they assigned to larvae becoming more vulnerable to predation by capelin because their length was below the optimum in the predator size–prey size relationship (10%; Paradis et al., 1996). In Trinity Bay, the majority of shanny larvae also had not yet reached a length corresponding to 10% of the modal length of capelin sampled using the IYGPT trawl (Figure 13). What factors could explain our pattern of selective loss from the population?

Figure 13

Relative length-frequency distributions of capelin (n=3341) and larval radiated shanny (2811) sampled (larval shanny standard lengths SL were transformed to total lengths TL according to TL=1.115 SL, R²=0.98, n=270).

Dower et al. (1998, 2002) concluded that shanny larvae were not highly susceptible to starvation because the incidence of feeding and the amount of food in the gut indicate that most individuals achieve daily ingestion rates sufficient to sustain growth. Therefore, starvation is not likely to cause selective loss of slower growing individuals. An alternative hypothesis is that the two distinct larval concentrations (western and northern cores) not only had different growth potentials but also experienced different predation mortalities. Radiated shanny are produced from demersal eggs laid in nearshore waters after which the larvae emerge and disperse throughout coastal waters (Green et al., 1987). Data from the second survey indicate that smaller individuals predominated in the western area whereas larger individuals were more abundant in the northern core. In addition, integrated surface temperatures were generally lower along the western shore than in the offshore area of the northern core. In addition to the effects of temperature on growth rate, the median potential for encounter with capelin was more than 10 times higher for larvae from the western core than for those from the northern core. The high concentration of larvae on the western shore during the second survey had disappeared during the third survey, whereas our forward simulations of particle drift predicted that this distinct cluster should have remained in the area. Although patterns of larval dispersal unresolved by the sampling programme or circulation model might have caused some error in our inferences, the evidence on hand suggests the spatial distribution of growth rates coupled with the spatial overlap with predators has resulted in a greater overall loss of larvae from one region over another. Although size-dependent processes may have occurred, these spatial differences may explain the overall selective loss of slower growing individuals.

This field study demonstrates that past rather than recent environmental variability is an important consideration if we are to understand the growth characteristics of larval populations sampled in the open sea. It has also shown some substantial difficulties still besetting such field studies, particularly the inability to represent larval prey abundance on appropriate scales, and the complexity of including predators such as pelagic fish. Broad regions of different predation pressure could be described, but to further understand predation, we need to move beyond descriptive studies and try to estimate the spatial structure of encounter probabilities between fast moving predators and the larval fish.

We thank T. Shears for his efforts to coordinate the logistic and technical activities in the field and in the laboratory, as well as G. Maillet, G. Redmond, and H. Maclean for their assistance in the field. C. Mercer is thanked for her invaluable guidance and help in the analysis of otolith microstructure. The work could not have been carried out without the assistance and expertise of the officers and crew of the Wilfred Templeman. This study received funds from the International Bureau of the Bundesministerium für Bildung und Forschung (BMBF, Project Nr. CAN 01/007) to HB, the Natural Sciences and Engineering Research Council to PP, and the Department of Fisheries and Oceans. The constructive criticism provided by two anonymous referees and the editor was helpful in improving the presentation of our findings.