## Abstract

The lack of a clear relationship between spawning output and recruitment success continues to confound attempts to understand and manage temporally variable fish populations. This relationship for a common reef fish is shown to be obscured by nonlinear processes in operation during the larval phase. Nonlinear responses of larval fish to their noisy physical environment may offer a general explanation for the erratic, often episodic, replenishment of open marine populations.

Understanding the nature of feedback control between the size of a population and its successful replenishment is fundamental to the understanding of the mechanics of population change. Fisheries assessors in particular use various models (for example, Beverton-Holt) to try to predict recruitment levels from reproductive stock size (1). Such fitted curves, however, have been plagued by extreme variance in recruitment, and explaining this bivariate scatter has remained a largely unrealized goal for the last 100 years. Some researchers have attributed the unexplained variability to sampling error and the misuse of analytical techniques (2), while others have suggested that unpredictable and high levels of mortality during an intervening larval phase will necessarily obscure the stock-recruitment relationship (3). Because an understanding of the fate of larval fishes has remained elusive, it has proven difficult to evaluate the relative contributions of these sources of error. Here, we offer a new approach to investigating processes at work during the larval phase.

For many benthic marine species, observations of rates of larval supply have been recorded for considerable lengths of time (4). Such records are typically erratic, dominated by a small number of high-magnitude events (5). Attempts to relate larval numbers to concurrently measured environmental variables using linear techniques—that is, regression analysis—have been largely unsuccessful. Episodes of strong larval supply have often been termed random and thought to be unpredictable with available oceanographic and meteorological data. However, we provide evidence that pulses in larval supply may arise from a nonlinear response of larval fishes to a relatively small number of key physical forcing variables.

We focused on a population of *Pomacentrus amboinensis*, a common damselfish species found on the Great Barrier Reef, Australia. Daily observations of spawning output (visual census of eggs released from nests) and recruitment (counts of recently settled juveniles) have been recorded for 2 years, and larval supply (mature larvae collected with light traps) for 3 years, at Lizard Island. As is characteristic of fisheries problems, density-dependent feedback between stock size and subsequent recruitment accounts for only ∼10% of the variation in the relationship between these two variables at the scale of a single reef (6).

We applied techniques from nonlinear time series analysis to determine whether the fluctuations in the larval supply time series are indicative of the action of nonlinear processes or are random features of the data. The approach was to construct a series of algorithms, ranging from global/linear to local/nonlinear, to forecast the data based on lagged-coordinate embeddings (7). The underlying philosophy is that, if nonlinear algorithms are found to outperform their linear counterparts, this may be taken as evidence that the processes that generated the time series were themselves nonlinear (8). Alternatively, if the underlying processes were linear or stochastic, linear models should provide better forecasts, as they use global information to generate predictions. For forecasting algorithms, we used the method of S-maps (9). S-maps generate forecasts for each point in the embedded time series by first excluding that point (the predictee), and then fitting a linear surface to the remaining data (the predictors). This technique is similar to linear kernel regression (10) in that the degree of nonlinearity in these mappings is given by the extent to which the contribution of predictor points to a fitted surface is a function of their distance from the predictee. This weighting procedure is controlled by a tuning parameter theta, where theta equal to zero gives a global linear map, and increasing values of theta give increasingly local and nonlinear algorithms (11).

Results of the S-map analysis are given in Fig. 1. Forecastability in both egg production and recruitment remained essentially constant across the spectrum from linear to nonlinear models; daily dynamics of spawning and recruitment were well described by a simple linear autoregressive model (12). However, for the larval supply series, prediction improved significantly (*P* < 0.05, Z test) as the model was tuned toward nonlinear solutions. This suggested that, in addition to purely linear or stochastic events, there was a nonlinear component of larval phase processes that was obscuring the spawning-recruitment relationship.

To test this, we examined three distinct physical-biological interactions: a lunar or semilunar entrainment of egg release from the nest, the transport environment experienced by mature larvae returning to the reef, and the effect of turbulence levels on young, first-feeding larvae (13). We chose, respectively, the following physical variables as representative of these processes: nighttime irradiance lagged 19 days (corresponding to average larval life), cross-shelf wind speed lagged one day, and average daily wind speed at a range of time lags. Model generation involved fitting linear surfaces; model nonlinearity was controlled by adjusting the neighborhood size to which the surfaces were fit (14). To increase the number of nonzero points in the time series, we considered all pomacentrid species in this analysis (15). Table 1 gives summary statistics for all the multivariate models constructed, and representative model forecasts are given in Fig. 2.

We began by modeling lunar entrainment, a deterministic component of larval supply. There was little improvement in model performance when nonlinear models were considered in comparison with linear regression. We then added cross-shelf wind speed, a stochastic component, into the model. In this case, nonlinear maps offered some increase in forecasting performance (Table 1) (16).

Choosing a physical variable to reflect events occurring early in larval life was more difficult. Such effects are plausible, given the suggestion that retention mechanisms may work to keep larvae local to their native reef (17), and otolith studies describing the age of these fish have been reported (6). However, because not all larvae are precisely the same age, and further because both the time required for yolk-sac resorption and the time to inevitable mortality in the absence of successful feeding are unknown, we proceeded conservatively by considering a wide range of lags (from 1 to 30 days) for daily averaged wind speed. Viewed linearly, there was no significant difference in model performance between any of these lags.

We then tuned all 30 lags toward nonlinear model solutions, creating a rho-lag-neighborhood surface (Fig. 3A). The slight improvement in correlation coefficient at intermediate lags, which was evident in linear models, increased significantly (Z test;*P* < 0.05) as the models were made more nonlinear. Notably, only models coinciding with the age of the larvae offered improvement over the bivariate nonlinear model. Otolith analysis of these fishes has indicated that they are from 16 to 23 days old, with a mean age of 19 days; the strongest effects we identified ranged from 15 to 17 days and centered at 16 days, 3 days after release from the nest.

To extract the relationship between larval supply and each of the three forcing variables in this model, we employed a neighborhood averaging procedure (18). For nighttime irradiance (Fig. 3B), maximum larval abundance was associated with a full moon 19 days prior. For cross-shelf wind speeds (Fig. 3C), larval abundance was positively associated with weak onshore winds. Recent findings have demonstrated that mature reef fish larvae are strong swimmers (19), and we hypothesize that lower wind speeds may promote easier orientation and increase directional swimming in larvae as they approach benthic habitats for the first time. The relationship between larval supply and 16-day lagged wind speed is shown in Fig. 3D. Recent hypotheses have proposed a dome-shaped relationship between larval feeding and turbulence levels (20). Our analysis supports this, and we suggest that wind effects are indeed important to early larval survival. In this case, the optimal window for early larval survival appears to have been approximately 4 to 8 meters per second.

It seems clear that, in this system, variability in larval supply originates from both deterministic (lunar entrainment) and stochastic processes (wind stress). The key elements of nonlinearity appear to be the response of the larvae to these stochastic variables and the sequential nature of processes operating at different times in larval life. Such sequential action serves to couple physical variables, causing their total impact on larval supply to be multiplicative. In other systems, one may expect the most relevant physical variables and lags, and the responses to those variables, to vary. However, it is likely that these key elements of nonlinearity will remain. We therefore offer a dynamical explanation for the episodic nature of larval supply: the nonlinear amplification by biological processes of stochastic physical forcing.