Biological processes and physical oceanography are often integrated in numerical modelling of marine fish larvae, but rarely in statistical analyses of spatio-temporal observation data. Here, we examine the relative contribution of inter-annual variability in spawner distribution, advection by ocean currents, hydrography and climate in modifying observed distribution patterns of cod larvae in the Lofoten–Barents Sea. By integrating predictions from a particle-tracking model into a spatially explicit statistical analysis, the effects of advection and the timing and locations of spawning are accounted for. The analysis also includes other environmental factors: temperature, salinity, a convergence index and a climate threshold determined by the North Atlantic Oscillation (NAO). We found that the spatial pattern of larvae changed over the two climate periods, being more upstream in low NAO years. We also demonstrate that spawning distribution and ocean circulation are the main factors shaping this distribution, while temperature effects are different between climate periods, probably due to a different spatial overlap of the fish larvae and their prey, and the consequent effect on the spatial pattern of larval survival. Our new methodological approach combines numerical and statistical modelling to draw robust inferences from observed distributions and will be of general interest for studies of many marine fish species.
Spatio-temporal variability in distributions of early life stages (ELS; i.e. eggs and larvae) of marine fish is shaped by complex interactions between the effects of abundance, fecundity and behaviour of spawners, ocean circulation patterns, and spatial and temporal patterns of mortality (which are in turn dependent on availability of food and predator aggregations). These interacting processes operate across various scales to establish the seemingly stochastic patterns observed . Further, changing climate regimes have the potential to disrupt these interactions. Climate influences distributions, for instance, by affecting larval dispersal pathways (through speed and direction of currents) , through temperature-dependent processes (growth, food availability, metabolic rates and energetic cost of larvae decisions) , and by affecting spawner distribution and spawning phenology .
Great progress has been made within the field of biophysical modelling of ELS of marine fish since its humble beginnings some decades ago. The development of computational tools allows linking high-spatial-resolution numerical (hydrodynamic) models with individual-based models (IBMs). This coupled approach allows for simulating oceanographic transport and ELS dynamics (i.e. growth, mortality and behaviour) to investigate the processes shaping fish distributions (see reviews in [5–8]). Independently, recent studies have demonstrated the importance of including spatial aspects in statistical analyses of inter-annual dynamics of plankton , ELS [10,11] and fish juveniles . Particle-tracking/IBMs of ELS often use observational data to calibrate the biological parameters of the models in order to be able to qualitatively reproduce the observed patterns (e.g. [13–15]). In some cases, particle-tracking estimates have been used to throw light on the causes of inter-annual variability in recruitment . However, despite the valuable advances made within particle-tracking/IBMs of ELS, the output is rarely analysed in unison with observational data to draw quantitative inferences by means of statistical methods. Here, we apply spatial statistics to link modelled distributions with measured observations, and evaluate which factors influence spatial and inter-annual variation in larvae abundance.
Northeast Arctic (NEA) cod (Gadus morhua) in the Lofoten–Barents Sea (figure 1a) is out used as our study system. NEA cod population dynamics and the spatio-temporal pattern of juvenile survival have been shown to change with climate [18,19]. Biophysical modelling studies have also confirmed that NEA cod larvae drift routes and individual growth can vary inter-annually with hydrographical regimes [2,20,21]. NEA cod is thus a well-suited model organism to test methods for determining the importance of biotic factors and environmental conditions in modifying the spatial distribution at the larval stage, and whether the effects differ with climate.
We investigate the large-scale spatial patterns of larvae advected from their spawning grounds by currents, and how these patterns are affected by spawner distribution, larval mortality and different environmental covariates of the larval habitat. Our approach is to link the simulated distributions of virtual particles (hereafter named drifters) with the observed distributions of larvae in July from 1986 to 1991. To do that, we incorporate drifter abundance for each year and grid location as covariates in a statistical spatial analysis, and include temperature, salinity and an index of convergence/divergence (CD) activity as additional covariates. Additionally, the study period comprises contrasting conditions of the North Atlantic Oscillation (NAO ). The NAO refers to inter-annual and inter-decadal variability in atmospheric forcing, which affects sea temperature, the strength of inflow to the Barents Sea and the population dynamics of NEA cod [17,19,22] (figure 1b). Consequently, the larvae experience both spatially and inter-annually variable ocean climate. To investigate possible non-stationary effects of the investigated predictor variables, we hypothesize that their functional forms may change under different climate periods.
2. Material and methods
(a) Biological data
Abundance data for cod larvae were obtained from yearly cruises carried out off northern Norway and in the western Barents Sea during the months of June–July from 1977 to 1991 using a pelagic trawl, sampling a depth from 50 m to the surface. The trawl net codend had a 4 m liner with 5 mm (stretched) hexagonal mesh . Note that, for simplicity, we refer to all sampled young-of-the-year fish as larvae, although an undetermined proportion of them had metamorphosed to a post-larval stage [24,25]. Their lengths varied between 12 and 54 mm. Acoustic spawning aggregation data were obtained during annual spring surveys carried out from 1986 and onwards off the Lofoten archipelago (in the centre of the blue box in figure 1a; figure 2a), where spawning lasts from early March to mid-May (see electronic supplementary material for details). The Lofoten area was divided into three sub-areas (Vestfjorden, Yttersida and Nord; figure 2a) based on Ellertsen et al. . Centres of gravity (CsG; abundance-weighted mean latitude and longitude) of the spatial distribution of spawners were calculated for each sub-area and year. These CsG were used as starting points for the release of drifters (i.e. eggs) in the particle-tracking model (see details below).
(b) The ocean circulation model
Data on ocean currents, sea temperature and salinity were obtained from a coupled model system consisting of an ocean circulation model (MI-POM) and an ice model (MI-IM). MI-POM is the Norwegian Meteorological Institutes version of the Princeton Ocean Model . The MI-IM model is described by Røed & Debernard . The coupled model system has a nested set-up as illustrated by Melsom & Fossum . The coarse grid domain covers the Arctic Sea, Nordic Seas and the North Atlantic Ocean, with a resolution of 20 km. The finer grid domain within this has a mesh size of 4 km, and covers the Barents Sea and most of the Norwegian Sea (the inner black rectangle in figure 1a). Owing to storage limitations, model results were interpolated to an 8 km grid. Note that the 8 km grid is the geographical reference used in this study and the basis for the statistical spatial analysis. More information about the model initialization and data assimilation is provided in the electronic supplementary material.
Since ocean circulation has a non-deterministic response to atmospheric forcing, a 10-member ocean circulation ensemble was run, where each member only differed in the atmospheric forcing. The ensemble was constructed by using forcing fields from the 10-day forecast produced twice a day at European Centre of Medium-Range Weather Forecasts (http://ecmwf.int). Detailed information about the construction and the validation of the ocean model ensemble is available in .
For all environmental variables used as covariates in the statistical analyses, one field per year was calculated for fixed positions and averaged over the 10 ensemble members over the period 23 June to 7 July at 10 m depth. For salinity and temperature, an averaged field per day was first calculated for the 10 ensemble members, and an average over the two weeks was later computed. A CD index was calculated as a proxy of potential retention directly from the velocity field with the function: 2.1where , and u and v are the x- and y-components of the current field. Because CD values at small spatial scales can change very rapidly over a short time period (i.e. hours), we used the standard deviation of CD over the aforementioned two weeks as the best proxy of vertical activity in the study period.
(c) The particle-tracking model and initialization of virtual drifters
Egg/larvae drift was computed based on a semi-Lagrangian advection model forced by currents from the ensemble members. Particle trajectories were first computed based on velocity fields from each of the ensemble members separately and then merged together to represent the average particle concentration from the 10 realizations. These averaged results were used in the statistical analysis. Owing to the fine mesh size of 4 km for the model, high-resolution current fields were simulated, which resolved ocean eddies. However, there are processes at still finer scales that were not resolved. These processes were parametrized by adding Gaussian noise to the model results (see electronic supplementary material for details). Owing to a potential overestimation of on-shore transport in the near-shore areas, a padding zone of 0.4 km was included in the particle-tracking model when an advection time step of 30 min was applied. Particles that were within this distance from the coastline, and not stranded, were moved outside the padding zone.
In the particle-tracking model, drifters were released with a Gaussian distribution centred at the CsG with a standard deviation of 4 km (0.5 grid units). For each sub-area and for each year, 15 000 drifters were released around the sub-area's CsG for that year (figure 2a). To account for variation in the timing of spawning, drifters were released on 15 March, 1 April and 15 April. The spatial transport of ELS of cod is also depth-dependent. At 10–20 m depth, the Norwegian Coastal Current (NCC) retains almost all the pelagic and passive-drifting stages close to the coast. Above 10 m, ELS are probably advected northwards by the Norwegian Atlantic Current (NAC; figure 1a) [2,20]. Thus, to account for depth variation, drifters were released at depths of 10 and 20 m. As there is little evidence of diel migration owing to 24 h daylight during summer , drifters were passively transported at these fixed depths by the NCC and the NAC based on results from the 10-member ocean model ensemble. The resulting larvae distribution was recorded at the time of the larval survey (around 1 July; i.e. after 90 days of drifting). In total, 270 000 drifters were released each year and standardized afterwards based on findings obtained from Ellertsen et al. . The standardization took into account (see the electronic supplementary material): phenology of spawning, depth distribution (assumed constant across years), and inter-annual variation in the distribution of spawners between sub-areas and in total spawner abundance (figure 2b).
(d) Statistical analysis of cod larvae distribution
Prior to analysis, the natural logarithm of observed larvae abundances (plus a constant, k = 0.5) was interpolated by means of ordinary kriging over a regular grid. To allow direct comparison of the different types of data, we used the same standardized grid (8 × 8 km; figure 1a) for the kriging as for the ocean circulation model. Kriged data for cod larvae were regressed against co-located covariates using generalized additive models (GAMs) . We used two types of GAM formulations: (i) purely additive, assuming that the effect of each covariate was stationary (that is, that the functional form of the effect did not change over the years), and (ii) threshold (non-additive), to test the hypothesis that the effect of a covariate changed in relation to an external threshold variable (e.g. different climate periods). The additive formulation used was 2.2where Lt,(x,y) is kriged cod larvae abundance (log-scale) in year t at location x,y in the grid, at is a year-specific constant, s is a two-dimensional non-parametric smoothing function describing the geographical effect (thin plate regression spline with maximally 27 degrees of freedom (28 knots)), f, g, h and j are one-dimensional non-parametric smoothing functions (cubic splines with maximally three degrees of freedom (four knots)), D is the natural logarithm of the density of drifters (plus a constant, k = 0.5) from the particle-tracking model, T is temperature, S is salinity, CD is the CD index calculated for the given grid location over two weeks, and ɛ is an error term.
To test the hypothesis of different effects of covariates under different climatic periods based on the NAO index, we fitted non-additive threshold GAMs (TGAMs)  to the data. TGAM is a semi-parametric regression where the shape of the exploratory function may change according to whether an external covariate is below or above a threshold value. The change in functional form may occur in one or several of the terms of the TGAM (in our model, from s1 to s2, from f1 to f2, from g1 to g2, from h1 to h2 and/or from j1 to j2): 2.3
We initially applied an NAO index threshold, n, on each function separately (i.e. on effects of location, drifter abundance, temperature, salinity and CD index), while the other functions were identical in low and high NAO regimes. We fitted series of GAMs and TGAMs with different levels of complexity (i.e. from one to five covariates and with all possible combinations). We considered all threshold values of NAO (=0, 0.5, 1, 1.5 and 4) that could uniquely delimit the potential temporal regimes from 1986 to 1991. Each threshold value was individually applied to each of the covariates for each model, resulting in a total of 114 models considered. The best model, GAM or TGAM, was determined by computing the (genuine) cross-validation (CV; a random sample of 10% of the data were used for these calculations, but threshold values were kept fixed; see  for details). Low CV values indicate the best model compromise between model complexity and fit to the observed data. After identifying the best models with only one threshold variable, we considered selected models with up to two threshold variables, as described in §3.
The residuals did not show temporal dependency but were spatially correlated. This is unlikely to bias estimates of the additive effects, but results in an overestimation of the significance level of the covariates, thus invalidating the standard p-values . We used a non-parametric wild bootstrap approach to get an accurate estimate of p-values and confidence intervals of covariates . For the best-fitted model based on the CV minimization, the residuals were extracted, rescaled to the same variance as the estimated scale parameter of the model and randomly replaced among years (sampling with replacement), before being used as response variable in the model (having the same set of covariates as the original model; see details in ). As an additional test of the significance of the threshold effects, we used an alternative CV approach, in which we left out data for whole years at a time (instead of a random 10% of the data). We then confirmed that the final selected threshold model (but without a year-specific intercept) had higher predictive power for data cases (years) not used when fitting the model compared with alternative formulations without thresholds.
(a) Spawning aggregations
Both the spatial distribution (CsG) and abundance within sub-areas varied among years (figure 2a,b; see spawner distribution and CsG for all the years in the electronic supplementary material, figure S1). The CsG were generally more offshore from 1987 to 1989 compared with 1986, 1990 and 1991 in Yttersida and Nord. Additionally, CsG were in the outer area of Vestfjorden from 1987 to 1989 and 1991. Pairwise correlations among CsG showed inter-annual consistency in the spatial variation: x-axis positions (note spatial coordinates refer to the oceanographic model domain; see figure 1a) were highly and positively correlated between Vestfjorden and Yttersida (Spearman, r = 0.92, p < 0.05), while y-axis positions were negatively correlated between Nord and Vestfjorden (r = −0.74, p < 0.05). With the exception of 1987 and 1991, the relative contribution of each sub-area to total annual cod abundance in Lofoten was quite constant (figure 2b). In 1987, 51 per cent of the fish spawned in the Nord sub-area, compared with 15 to 25 per cent for most years. In 1991, the year with highest total spawner abundance, only 6 per cent spawned in this sub-area.
(b) Particle-tracking model
Maps of simulated larvae distributions after 90 days of drifting from the mean date of spawning (1 April) and from each sub-area show only the predicted effects of inter-annual variation in the current pattern (see electronic supplementary material, figure S2). By taking into account variation in the relative contribution of each sub-area, the model predicted extensive spreading of drifters in the Barents Sea for 1987 (see electronic supplementary material, figure S3a). After taking into account variation in total spawner abundance, extensive spreading was also predicted for 1991 (electronic supplementary material, figure S3b). Higher abundance of drifters in the near-shore areas is a general pattern in the distributions obtained.
(c) Spatial analyses
We analysed the distribution of larvae estimated from survey data (figure 3a; see all the years in the electronic supplementary material, figure S4, and its inter-annual variation in electronic supplementary material, figure S5) using additive or non-additive models with one to five covariates. Besides drifter density estimates described above (figure 3b; all years are shown in the electronic supplementary material, figure S3b), temperature, salinity and CD index were used as covariates (figures 3c–e; all years are shown in the electronic supplementary material, figure S6). Results show that models with non-additive threshold effects (TGAM) were more parsimonious (as given by CV values) compared with purely additive models (see electronic supplementary material, table S1). The most parsimonious TGAMs generally identified a threshold of NAO = 1, corresponding to two short-term climatic periods: 1986–1988 (NAO < 1, ‘cold’ period) and 1989–1991 (NAO > 1, ‘warm’ period). The preliminary best models included all five covariates and an NAO = 1 threshold on either: (i) spatial location, the best model (see electronic supplementary material, table S1, model 80); (ii) temperature, the second best model (see electronic supplementary material, table S1, model 76); or (iii) drifter abundance, the third best model (see electronic supplementary material, table S1, model 78). These models were sequentially combined, applying an NAO = 1 threshold on two or three covariates, and resulted in three models (models 81–83). Two of these models included salinity, although salinity was not significant (p = 0.18) when applying a non-parametric bootstrap on model 81. Spatial heterogeneity in the correlation between salinity and temperature as shown in our study (figure 3f) justifies the inclusion of both temperature and salinity in our models. The final model (see electronic supplementary material, table S1, model 83; figure 4) explained 84 per cent of the observed variance in larvae abundance.
The final model is shown in figure 4. The partial effect of year (p < 0.05; figure 4a) shows an increasing trend, partly coinciding with the trend in spawner abundance (figure 1b). The relationship between drifter abundance and observed larval abundance is asymptotic (p < 0.05; figure 4b). Temperature effects differed between the climatic periods (1986–1988 and 1989–1991; the two effects p < 0.05; figure 4c). Specifically, during the warm period, high abundances of larvae were found in waters with temperatures around 9.5°C, while in the cold period the effect of temperature was stronger and the highest abundances were found in temperatures above 10°C. The effect of the geographical location showed a clear and significant difference between periods, suggesting positive anomalies in the southwestern parts of the study area in the period 1986–1988 (p < 0.05; figure 4d) and in the northeastern parts in the period 1989–1991 (p < 0.05; figure 4e). Note that these effects remain while accounting for local environmental covariates. Finally, there was some evidence for an effect of CD activity (p = 0.073; figure 4f). To assess the unique contribution of each covariate, we calculated the difference in percentage of variance explained between the full model and a model with the given term removed. The spatial term and the partial effect of year contributed the most (25 and 19 per cent, respectively), while temperature, drifter abundance and CD activity explained 11.8, 7.2 and 1.8 per cent, respectively.
Our study reveals how the abundance and distribution of spawners, ocean circulation and hydrographical environmental conditions in combination modify the spatial distribution of ELS, and how the effects of these factors differ between two periods with contrasting states of the NAO. The best model obtained explains a major fraction of the observed variance and includes drifter abundance (which integrates effects of spawner abundance and distribution, and larval advection) and the CD index as additive covariates, and temperature and geographical location as non-additive covariates—their effects changing between a high and low NAO index. We show that by combining ocean circulation models and statistical techniques, we can reproduce the spatio-temporal dynamics of NEA cod larvae, and reveal the environmental factors shaping their distribution. Thus, this study provides a novel and promising methodological framework for integrating quantitative information on spawner aggregations with ocean circulation and hydrography in statistical spatial analyses of pelagic ELS.
(a) Implications of spawner distribution for advection and survival
The starting position of drifters fundamentally controls the directions and distance of their drift paths. Therefore, the spatio-temporal distribution of spawners (when, where and how many) is crucial for proper initialization of particle-tracking models . Although cod spawning phenology and geographical location is a fine-tuned evolutionary strategy developed through generations to maximize survival , large-scale inter-decadal displacements of spawning aggregations have been described and attributed to both climate  and fishing effects on size structure . From the mid-1980s, a warm period during which cod populations were demographically eroded, the Lofoten grounds were the most important area for the spawner aggregations. Our results add to this picture by emphasizing the effects of small-scale variation within the Lofoten archipelago. The inter-annual variation in the CsG was consistent among sub-areas, indicating a certain spatial synchrony in spawning aggregation patterns. However, we cannot establish whether the temporal variation of the fine-scale structure in the spawning areas was related to climate, mating behaviour or changes in the population density over time.
Our results confirm that spatio-temporal variation in spawner aggregations within and among sub-areas around Lofoten combined with varying circulation patterns impact the resultant spatial distribution of eggs and larvae. In 1987, for example, the higher contribution of spawners at Nord (the more offshore sub-area) facilitated a more extensive spreading of ELS. Such spreading may ensure that a sufficient number of eggs and larvae encounter environmentally favourable conditions ( and references therein), and thus may have contributed to the higher larval survival found that year (figure 4a). High survival in 1991, when high spawner abundance caused extensive spreading of drifters into the Barents Sea, is consistent with this explanation.
(b) Spatial distribution of early life stages of cod
The asymptotic relationship observed between drifter abundance and observed larval abundance could be related to lower survival of ELS very close to the coast, which is where the highest concentrations of drifters occurred. Stranding of eggs onto the coast, higher predation during the drifting process and higher mortality at high larval densities (e.g. because of food limitation) are all factors that may have contributed to low survival. However, it is worth noting that the unresolved currents in the coastal areas owing to the confluence of multiple drivers interacting at fine to meso scales near shore could potentially overestimate drifters along the coast . Although a padding zone was included in the particle-tracking model to mitigate this problem, overestimation of drifter density along the coast is certainly possible. Nevertheless, the spatial pattern of our results is consistent with other studies applying different ocean models [2,20].
Having accounted for drifter abundance (i.e. the combined effect of spawner distribution and larval advection), the statistical model captures the potential effect of the other covariates on larval distribution. The most novel and striking result is the different contribution of temperature and spatial location under two short-term periods with contrasting climate conditions (1986–1988 and 1989–1991). Previous studies have shown that both temperature and recruitment success of cod in the study area are higher during NAO-positive years . Our analysis adds to this picture by showing that the generally positive temperature–larval abundance relation was weaker in years with high NAO values. During the generally warm (NAO-positive) 1989–1991 period, high abundances of larvae were found in water masses with a relatively wide range of temperatures. During the colder 1986–1988 period, high larval abundances were only found in the warmest waters. The causal basis for this difference is not clear, but we note that temperature is positively associated with food availability and that there is a better temporal match between the development of cod larvae and their zooplankton prey in warm years [24,37]. We speculate that in cold periods, zooplankton of appropriate size is only available in sufficient quantity in the warmest water masses, whereas in warm periods, zooplankton conditions are favourable over larger areas (and a wider temperature range). This is supported by recent evidence showing that zooplankton biomass predicts survival of ELS of NEA cod , and how prolonged temporal overlap between fish larvae and their prey in warm years enhances cumulative growth and survival .
Lastly, the partial effect of geographical location of larvae shows a clear and significant difference between the two periods. Based on earlier work, the southwestern distribution of larvae could potentially be explained by a reduction in transport during cold periods (NAO-negative years) [2,17]. However, as far as the drifters capture the current transport, the spatial pattern summarizes the variability left over by the selected environmental covariates, and therefore could indicate that the spatial pattern in mortality differs between the two periods. Several factors could cause such a shift in mortality, such as changes in the distribution of predators or prey. Higher zooplankton biomass in NAO-positive years has indeed been shown for the ‘downstream’  but not the ‘upstream’ (L. C. Stige 2011, unpublished data) part of our study area, consistent with such a role of prey. A potential role of predators is supported by the results of Ciannelli et al.  suggesting that cannibalism may modify the distribution pattern of 1-year-old cod in the Barents Sea. On the other hand, the climatic influence on the transport cannot be disregarded. For example, during windy conditions (more frequent in NAO-positive years) eggs are mixed deeper in the water column, where NCC favours a more downstream distribution . Further, our results also showed that CD activity could affect the observed larval distribution. The highest CD is found in the southern part of the shallow Barents Sea, where the NAC is deflected northwards and the NCC spreads out from the coast (figure 1a). This process contributes to a front developing between the water masses, creating a convergence zone with increased concentration of phytoplankton and zooplankton . Finally, CD activity in this area could slow down the northward transport of larvae and reinforce the southern distribution observed in cold periods.
(c) Implications and future challenges for an integrative approach
Disentangling the drivers of larval dispersal and the consequences for the spatio-temporal patterns of ELS distribution is an inherently biophysical problem. This requires an interdisciplinary effort, integrating tools from different scientific disciplines. This study provides a new biophysical framework to link observational information from spawner distributions to larval distributions, taking into account variability in ocean circulation and the hydrographical environmental conditions of the larval habitat. Although we provide some methodological advances, such as the use of a model ensemble to include uncertainty into the numerically modelled larval distributions and covariates, there is also room for improvement. Key challenges for future studies include improvement in the near-shore oceanographic modelling, accounting for the size structure of the mature stock, inclusion of larval behaviour linked to observations (e.g. vertical migration linked to observations of food availability) and accounting for the lack of synopticity in the ELS surveys.
Although the circulation pattern is the main physical factor shaping the spatial distribution of larvae, climate indirectly modifies this distribution by affecting larval survival (e.g. [25,37]). It is worth noticing that inter-annual variation in spatial distribution of spawners at scales of less than 100 km can affect the spread of larvae into the larval habitat. In essence, our study highlights that larvae were distributed over a wider temperature range, over a larger area and more downstream in warm periods owing to a combination of (i) higher abundance of spawners, (ii) an increase of ocean transport, (iii) higher overall survival (possibly because of better food availability and/or escape from predation), and (iv) a likely change in the spatial pattern in natural survival, hypothetically caused by changes in spatial distribution of predators or prey. ELS of many marine fish species experience different environmental conditions throughout their distribution ranges. The combination of hydrodynamical and statistical modelling is a promising tool for integrating information about oceanographic processes, small-scale ecological mechanisms and large-scale observations to draw statistically robust inferences about the processes shaping ELS distribution and survival. Our study shows that such a tool is likely to provide novel knowledge, even for highly well-studied species and systems.
M.H., Y.G., D.H., G.O. and L.C.S. received support from the Research Council of Norway (projects no. 173487 and 178434). M.H. thanks M. Llope for the advice and discussions on an early version of this manuscript and for help in the data analysis, Lee Hsiang Liow and Jason D. Whittington for reviewing an advanced version, and two anonymous reviewers for their helpful comments.
- Received April 10, 2011.
- Accepted May 27, 2011.
- This Journal is © 2011 The Royal Society