Introduction

The larval phase is the main dispersive stage of most demersal teleost marine fishes, controlling population dynamics and shaping connectivity patterns. As such, it plays a key role in large-scale ecological processes such as gene flow and biogeography1. The dispersal of larval fish is governed by two main mechanisms, ocean currents and larval behavior1. Although our understanding of ocean currents has greatly improved in the past decades allowing better reconstruction and prediction, a proper understanding of larval behavior has been more challenging due to a high degree of uncertainty and inter- and intra-specific variability in larval traits2.

In the past two decades, multiple studies demonstrated substantial swimming and orientation capabilities for larval fish of various species, which can affect their dispersal outcome3,4. However, most of these empirical studies focused on a single species at a single location, such that a generalized cross-species inference has never been attempted. This might be the reason why larval behavior is not implemented in most biophysical larval dispersal models5. Here, we apply a cross-species meta-analysis focusing on the fundamental behavioral trait of larval directional swimming.

Recent studies have repeatedly demonstrated that fish larvae influence their dispersal by swimming directionally3. Directional movement is a central component in animal movement ecology6,7, referring to the tendency of an individual to move along a straight path3,7. A movement is considered directional if it has a significant directional precision (or mean vector length), based on Rayleigh’s test of uniformity8. Directional movement is different from directed or oriented movement, in which by definition, there is an inherent use of external cues9. Similarly, unoriented movement is a situation in which there is no use of external cues for orientation.

Directional movement can be achieved using internal stimuli or with reference to external cues, with only the latter representing oriented movement6. It is currently unclear which type is used by fish larvae. This is critical since larval dispersal is a key process governing demographic connectivity, gene flow and biogeography of marine populations10. Incorrect representation of larval orientation in biophysical models can lead to inaccurate estimations of larval transport and connectivity11. To test whether fish larvae use external cues for directional movement, we compare observed movement patterns to those expected under a strict use of internal cues. Our analyses are based on the simplifying assumptions that a strict use of internal cues is expressed in a Correlated Random Walk (CRW) process, whereby consecutive movement directions or ‘bearings’ are auto-correlated11, and that oriented movement patterns, measured across multiple time steps (>20), are more directional (i.e., straighter paths characterized by higher mean vector lengths) compared to unoriented patterns12,13. Under these assumptions, we find a robust support for the use of external cues by fish larvae, both at the individual and at the species levels. We discuss cases that diverge from these assumptions, for which our methodology is inappropriate.

Directional movement has been demonstrated in fish larvae of more than 20 species3,14,15,16,17,18,19,20. The studied larvae were mostly wild larvae captured by light traps, although some were reared (Table 1). All were in the post-flexion stage of development, meaning the caudal fin was formed, and they had the ability to swim at speeds that meant they were moving in an inertial hydrodynamic environment21,22. Two formats of field trials were used to assess larval directional movement: (i) tracking larvae by scuba-divers (Scuba-Following23) and (ii) video recording of larvae inside the Drifting In Situ Chamber (DISC24). In both types of trials, the observed precision could be achieved by utilizing external (i.e., oriented) or internal (i.e., unoriented) cues. The difference relies on the orientation mechanism: larvae that utilize internal cues may use their proprioceptive system, similar to a gyroscope25, or simply keep moving toward the same general direction due to inertia26, giving rise to CRW11. In contrast, truly-orienting larvae use external cues such as the earth’s magnetic field27, resulting in patterns such as Biased Random Walk (BRW), whereby bearings correlate with a fixed external direction6,11. The distinction between the two is critical since external cues enable a more persistent directional movement over time, and it allows for corrections should the larvae be displaced6.

Table 1 Species used for the meta-analysis and relevant statistics comparing \({\hat{R}}_{\theta }\) quantile distribution with the null \({R}_{{\theta }_{0}^{{vm}}}\) and \({R}_{{\theta }_{0}^{r}}\)quantile distributions using chi-square (X2) goodness-of-fit test.

In the present work, we show that the mean vector lengths larval-fish field trials were significantly higher that those expected under CRW both at the individual and at the species levels, providing supporting evidence the use of external cues for directional movement by fish larvae.

Results

Out of 832 examined trials, 755 (91%) were indicated as directional (Rayleigh’s test, p < 0.05), and were thus considered for the quantitative analyses (Table 1). The results of the Correlated Random Walk-von Mises (CRW-vm) analysis (detailed in the Methods section) show that the observed mean vector length (\({\hat{R}}_{\theta }\)) of all tested species exceeded the mean expected under CRW (\({\bar{R}}_{{\theta }_{0}^{{vm}}}\), Fig. 1a, b, Table 1). DISC trials carried out in the Mediterranean Sea (e.g., Chromis chromis), exhibit generally lower \({\hat{R}}_{\varDelta \theta }\) compared to trials conducted in the Red Sea and Great Barrier Reef (e.g., Chromis atripectoralis; Fig. 1b, Table 1). Subsampled DISC trials exhibit similar quantile ranges to those of the scuba-following trials (dotted crosses in Fig. 1a), suggesting a solid comparability among the two sampling methods. In five species, confidence intervals overlap with \({\bar{R}}_{{\theta }_{0}^{{vm}}}\) (Fig. 1a, b, Table 1). Only two of these species, Chaetodon aureofasciatus and Amblyglyphidodon curacao do not display a clear indication for oriented movement, likely because of their distinctive depth-dependent behavior Supplementary note 1, Supplementary Figure S1. In the three other species (Platax teira, Epinephelus coioides and Caranx ignobilis), the small number of trials (2–15) per species may explain the confidence intervals overlap (Table 1).

Fig. 1: Results of Correlated Random Walk-von Mises (CRW-vm) and Correlated Random Walk resampling (CRW-r).
figure 1

CRW-vm analysis at the species level based on the diagram in Fig. 2c for the scuba-following trials (a), and for the DISC trials (b) with various number of samples-per-trial: Nobs = 21 (a), 300, 90, and 180 (b). Crosses represent means ± 95% CI of the observed (RΔθ, Rθ) pooled by species. Crosses with dashed lines in b represent species with Nobs different than 180, see Table 1 for more details. The 5th, 10th, 20th,…., 90th, and 95th quantiles of the \({R}_{{\theta }_{0}^{{vm}}}\) distribution is represented in grey contours, with thicker contours for the 5th and 95th quantiles. Thick blue line in a and b represents \({\bar{R}}_{{\theta }_{0}^{{vm}}}\). Species names are ordered top-down according to their Rθ means and correspond in color to their respective crosses in a and b. Crosses with dotted lines in a represent DISC trials, which were subsampled to Nobs = 21; for these species, colors match the names in b. The species’ full names are provided in Table 1. Chi-square goodness of fit plots comparing the \({\hat{R}}_{\theta }\) quantiles distributions (black bars) to the null quantile distributions (grey bars) for CRW-vm (c) and CRW-r (d). e A heatmap of all trial counts binned at 5% increments according to their \({\hat{R}}_{\theta }\) quantiles in the CRW-vm and CRW-r analyses.

Fourteen species have enough trials (Ntrials ≥ 20) for a within-species chi-square test (proportions: 0–50%, 50–75%, and 75–100%). Of these, eleven species exhibit significant indication for oriented movement for both CRW-vm and CRW-resampling (CRW-r) analyses (chi-square test, P < 0.05, Cohen’s W > 0.5; Table 1). The chi-square test results support a significant indication for straighter movement than expected under CRW as the mean of effect sizes (Cohen’s W) is significantly larger than 0.5 for CRW-vm and CRW-r (One-tailed one-sample t-test, t > 1.93, P < 0.05).

CRW-vm and CRW-r quantile analyses provide a significant indication for straighter movement than expected under CRW because the \({\hat{R}}_{\theta }\) quantiles are significantly skewed towards higher values compared to the null (\({R}_{{\theta }_{0}^{{vm}}}\) and \({R}_{{\theta }_{0}^{r}}\), Fig. 1c, d; chi-square test, P < 0.0001, Cohen’s W > 0.5). If larvae are using only internal cues for directional swimming, exhibiting CRW (i.e., unoriented movement), their density distribution should be centered around the medians (Fig. 1e, white circle). The high larval concentration on the top-right quarter of Fig. 1e provides a strong indication for straighter movement than expected under CRW, suggesting the use of external cues for directional movement by fish larvae (i.e., oriented movement).

In addition, CRW- wrapped Cauchy (CRW-wc) analysis–which includes a wrapped Cauchy fundamental distribution (instead of von Mises)–produces similar results to that of the CRW-vm analysis. Specifically, the means of all species fall above the \({\bar{R}}_{{\theta }_{0}^{{wc}}}\) curve, with mean quantiles higher than 50 (Supplementary note 2 and Supplementary Figure S2).

Discussion

The set of quantitative analyses used in the study indicate a significantly straighter movement of larval fish than expected under CRW. These analyses are based on comparing the observed movement patterns to null distributions expected under CRW, produced by two complementary approaches: theoretical fundamental distribution (CRW-vm and CRW-wc) and resampling (CRW-r). This combined approach provides support for oriented movement by fish larvae.

Previous work that quantitatively distinguished between oriented and unoriented movement, e.g., Correlated versus Biased Random Walk (CRW vs. BRW), have used: (i) step-length dynamics and their interactions with bearings28, (ii) displacement-related metrics for gauging oriented movement, e.g., net squared displacement12,29, or (iii) net displacement6,28. In contrast, we focus on bearings, constructing CRW-expected pattern. Specifically, we examine the movement dynamics based on the relationship between the Rθ and RΔθ. In our analyses, every larva is gauged against its own potential of autocorrelated directional swimming (CRW) and can be readily compared against other larvae from either the same or different species (Fig. 2). Note that our analyses do not distinguish between the different types of directed movement (i.e., BRW vs. BCRW). In addition, while previous methods determined indication for oriented movement if a given sequence exhibited orientation (e.g., displacement) above the 95% CI range expected for CRW, we employ a different quantitative approach in which we compare the observed distribution of trials’ quantiles to the distribution expected under CRW. The consideration of continuous quantiles (1–100%) rather than the dichotomous determination of within or outside the 95% CI range, increases the sensitivity of the analysis. Moreover, we provide a quantitative method for a meta-analysis of multiple individuals from multiple species that combines both theoretical and resampling approaches6,12,28,30,31.

Fig. 2: Protocol of the two quantitative methods applied in the study: Correlated Random Walk-von Mises (CRW-vm; left section), Sequence Resampling (CRW-r; right section).
figure 2

a Each trial includes a bearings sequence (θ; green), from which the turning angles sequence (Δθ; magenta) is computed. For visual clarity, a includes only seven bearings (Nobs = 7), while in the actual trials Nobs is larger. b (i) To generate the CRW-vm null distribution, we sample (with replacement) from a von Mises distribution, generating for each kappa value (K:0,1,2,3…, 399) 1000 Δθ sequences and corresponding θ sequences at the same length as Nobs per trial. (ii) For each of these simulated sequences, mean vector lengths (R) of Δθ and θ are computed (RΔθ, Rθ), such that a distribution of Rθ (grey vertical distributions in c) is generated for each kappa (yellow dots). The 5th, 10th, 20th,…., 90th, and 95th quantiles of each of these Rθ distributions is computed, and contours connect these quantiles across the different Rθ distributions (black contours in c), with the mean represented as a thick blue line. Samples with high \({\hat{R}}_{\theta }\) compared to \({R}_{{\theta }_{0}^{{vm}}}\), above and further than the blue line in the cyan colored area, represent a tendency for oriented movement. Complex movement patterns are not suitable for analysis using our methodology and can occur in the area below \({R}_{{\theta }_{0}^{{vm}}}\), e.g., one-sided bias, as well as in the area above \({R}_{{\theta }_{0}^{{vm}}}\), e.g., composite correlated random walks or bi-model turning angle distribution (zig zag) (for more details see Supplementary note 4). Note that the schematics in c represent the fact that the mean vector length (straightness) of the trajectory decreases from directed through CRW, to completely random (Simple Random Walk) or complex displacement-reducing movement patterns. (iii) (\({\hat{R}}_{\varDelta \theta },\) \({\hat{R}}_{\theta }\)) is computed per experimental trial and plotted on the phase diagram both individually and by species. (iv) Quantiles of the \({\hat{R}}_{\theta }\) within the \({R}_{{\theta }_{0}^{{vm}}}\) are computed using a 2-D interpolation. (v) Chi-square tests are applied to compare the \({\hat{R}}_{\theta }\) and the \({R}_{{\theta }_{0}^{{vm}}}\) distributions; and the effect size of the chi-square test (Cohen’s W) is computed. (vi) A one sample t-test is applied to examine the significance of the effect size (Cohen’s W > 0.5) across all species. Correlated Random Walk- resampling (CRW-r; right section). For each trial, Δθ is computed (a), as shown in d, which represents a histogram of Δθ from the example sequence Supplementary. Then, Δθ sequence is sampled without replacement 100 times (vii), and 100 θ sequences are generated (viii). Next, (e) R is computed for each of these 100 θ sequences (\({R}_{{\theta }_{0}^{r}}\)) (ix), and the quantile of the \({\hat{R}}_{\theta }\) within the \({R}_{{\theta }_{0}^{r}}\) distribution is calculated (iv), in this example quantile = 78%. Then, stages (iv–vi) are applied the same as for the CRW-vm method, ultimately computing the strength of the differences between the \({\hat{R}}_{\theta }\) and \({R}_{{\theta }_{0}^{r}}\) quantile distributions.

Movement simulations often assume a von Mises or wrapped Cauchy distributions of Δθ, which are widely used to simulate CRW processes6,11,13,22,32,33. Yet, these distributions may not perfectly represent the true underlying distribution of Δθ. The strength of the sequence resampling (CRW-r) analysis is that it makes no assumptions regarding the distribution of Δθ, but instead reconstructs θ sequences by resampling the empirical Δθ6,12,34. Since these sequences are often short (Nobs = 21), they may not properly represent the true underlying distribution. Yet, the combination of resampling and theoretical distribution methods provides a complementary view and both support oriented movement in fish larvae. A subsampling of trials with high number of observations (Nobs = 180) produce comparable results to those of trials with a low number of observations (Nobs = 21) in terms of the range of quantiles for a given species (Supplementary note S3, Supplementary Figure S3).

Previous implementation of oriented (e.g., BRW)22,35,36,37,38 and unoriented (e.g., CRW)22,32,39 movement patterns in biophysical models of larval dispersal, demonstrated a significant effect on dispersal trajectories, settlement success and connectivity patterns22,32,35,36,37. The mathematical methodology behind such implementation has been extensively described11. It is largely based on a sampling of swimming directions from a von Mises or wrapped Cauchy distributions centered around the direction of the cue source for oriented movement, or around the swimming direction in the previous time step for autocorrelated movement22,32,33. In both cases, kappa or rho, the concentration parameters of the von Mises and the wrapped Cauchy distributions, respectively, govern the simulated precision of directional movement. Several studies used empirical species-specific data from in situ orientation trials to parametrize orientation behavior22,36,40. Yet, the combination and inter-dependence between oriented and autocorrelated directional movements were not implemented. The methodology presented here provides the basis for such implementation given the consideration of both oriented and unoriented movement under a single framework.

Animal orientation patterns are variable between species, between conspecifics, and within a single individual, with uncertainty and variability associated with environmental factors, availability of orientation cues, and the internal state of the organism41. In this work, we do not distinguish between higher order behaviors such as Biased Random Walk (BRW), Biased-Correlated Random Walk (BCRW), complex movement patterns such as chemotaxis42,43 or infotaxis44, and combinations of movement strategies30,45,46. We also disregard the spatial context of the movement trajectories and the possibility that 3-D movement may influence resultant 2-D orientation patterns47. Such complex behaviors are inappropriate for testing using our methodology. For example, animals that exhibit zig zag movement, with turning angles drawn from two distributions (centered around negative and positive values) or CRW movement sequences that are based of variable kappa could exhibit higher mean vector lengths compared to those expected under CRW (see Supplementary note S4, Supplementary Figure S4, and Supplementary Table S1 for more details). Other types of behaviors could be manifested in a relative reduction of \({\hat{R}}_{\theta }\) compared to what is expected given \({\hat{R}}_{\varDelta \theta }\), as observed in C. aureofasciatus and A.curacao. Hence, the interpretations provided in our study are based on the simplifying assumptions that a strict use of internal cues is expressed in a CRW process, and that oriented movement patterns result in more directional (or straighter) paths compared to unoriented patterns. We tested the possible effect of irregular behaviors on our results by removing trials that contained such irregular patterns and rerunning the analyses on this subset. The results still show a significantly straighter movement than expected under CRW, suggesting oriented movement by fish larvae (Supplementary note 5, Supplementary Figure S5).

In addition, spatiotemporal variability of the external cues should be considered and further examined in future studies. In our analyses we assumed that directional cues would come from a fixed direction throughout each trial. However, this assumption may not necessarily represent the conditions in the field. Part of the limitation lies in the fact that we cannot know for sure what are the directional cues utilized by fish larvae. For example, this assumption is correct if a larva was following a celestial cue1,2. In contrast, if a larva was following an auditory cue coming from a coral reef patch, and the larva was swimming or drifting, the direction of the cue source would change across the trial. However, considering the short durations of our trials (<15 min), currents and larval swimming speeds (<25 cm/s), the distance from shore in which the trials normally took place (>50 m), and the fact that potential sources of directional cues are often distributed parallel to the shore (e.g., fringing reef), we do not expect this to have a strong influence on the resultant movement properties of our trials.

Moreover, it is important to note that larvae in the field are often found in groups22, and larval swimming and orientation in groups are more efficient compared to single larvae16,22. In our results it seems that Chromis atripectoralis groups exhibited straighter paths and higher tendency for directed swimming compared to individual larvae based on both experimental methods. Yet, it is unknown if this effect is present in other species as well.

Previous studies indicated similarity in orientation patterns for scuba-following and DISC trials, with scuba-following exhibiting generally higher \({R}_{\theta }\) values16,48. The results summary in Table 1 suggests that DISC trials exhibit a higher tendency for oriented movement compared to the scuba-following trials. However, this is likely a result of the difference in the Nobs. For DISC trials with low Nobs, the \({\bar{R}}_{{\theta }_{0}^{{vm}}}\)values are high, resulting in a more linear relationship, compared to trials with higher Nobs (Fig. 1a, b). Indeed, a subsampling of the DISC data to obtain Nobs = 21 results in a general reduction of RΔθ, Rθ, and the quantile of the trial. Notably, the scuba-following and subsampled DISC trials exhibit similar quantile ranges (Fig. 1a). In addition, it seems that DISC trials conducted in the Mediterranean Sea are characterized by lower RΔθ compared to DISC trials from the Red Sea and the Great Barrier Reef. This might be related to differences in cue availability, that is affected by environmental conditions such as water turbidity and overcast sky conditions19. The fact that both the original and subsampled DISC datasets exhibited quantile distribution significantly higher than what is expected under CRW further supports the indication for directed movement by fish larvae.

Previous studies demonstrated that increase in sampling length (or reduction in sampling frequency) leads to an obstruction of the short-term persistence, making it appear more BRW, and leading to a decrease in the animal’s apparent speed49,50. Similarly, the computation of movement properties (path length, or the difference between observed and CRW-expected path length) given an incremental change in step-length shows that a directed movement is manifested especially across large scales3,4. In our data, due to the limited number of observations (Nobs), such systematic analyses could not be applied. Future research should explore the effect of trial method, trial duration, sample frequency, and geographical location on the movement properties.

Overall, oriented movement allows a more persistent movement along a mean bearing over time compared to unoriented movement6,13,26,28,32. Such persistent movement results in an increased behavioral displacement relative to the water in which the larva swims, which in turn, modifies larval dispersal distance, recruitment success, and connectivity patterns5,22,32,35,36. Modifications are site- and species-specific, and are based on the fact that oriented movement allows the larvae to depart from entraining currents, and reach their settlement habitat more efficiently22,32,35,37. In our analyses, the high Cohen W effect size value (>0.5) suggests a high strength and consistency of oriented movement across individual larvae and pooled by species. It is important to consider intra-specific variation in orientation behavior, which may be a manifestation of bet hedging strategy51. Such variation should be further studied and implemented in biophysical models of larval dispersal.

Previous studies have demonstrated the need to consider larval orientation to understand the observed settlement success and connectivity40. It is therefore important to implement experimentally obtained orientation patterns into biophysical models, and the analytical approach presented here can increase the biological realism of larval dispersal models. Yet available orientation trials are of limited duration and longer trials are needed to study how orientation and its related cues change across time, space, and ontogeny. Future work should thus focus on isolating specific orientation cues, as well as studying cool-water, non-perciform taxa of commercial importance (e.g., Pleuronectiformes and Gadiformes) to promote a sustainable management of marine fish populations. Our meta-analysis suggests that fish larvae make substantial use of external cues for oriented movement, enabling them to find their way in the seemingly featureless, open ocean.

Methods

General methodological approach

To examine if larvae utilize external cues (i.e., oriented movement) to swim in a directional manner (i.e., significant mean vector length), we develop two complementary analyses that compare the empirically observed directional precision (i.e., mean vector length) with the null distribution expected under a strict use of internal cues (i.e., unoriented movement). The empirically observed directional precision is quantified as the mean vector length (R) of larval bearings (θ) (Fig. 2a), herein \({\hat{R}}_{\theta }\). The angular differences between consecutive bearings, herein turning angles (Fig. 2a; Δθt = θtt-1), are used to generate two null distributions of Rθ expected under the unoriented movement of Correlated Random Walk (CRW; \({R}_{{\theta }_{0}}\)), based on the two analyses: Correlated Random Walk-von Mises (CRW-vm) and Correlated Random Walk- resampling (CRW-r), described below. The first is theoretical and is based on a von Mises distribution of simulated Δθ (Fig. 2b, c); the second is empirical, and is based on resampling the Δθ within each trial (Fig. 2d, e). These two analyses are complementary because the first can generate an unlimited number of trajectories but is based on a theoretical distribution rather than on observations, whereas the second is based on a finite number of observations. In addition to these two main analyses, we apply a third analysis, the Correlated Random Walk-wrapped Cauchy, herein CRW-wc, which is similar to CRW-vm, with the only difference of using wrapped Cauchy distribution instead of von Mises. The reason for applying CRW-wc is that it was shown to represent well animal movement in some cases33. Notably, we consider the simple cases of undirected movement pattern with a turning angle distribution centered at 0 (CRW), testing if the mean vector length of the trial’s sequence is higher than that expected under CRW. If true, that would be an indication for a directed movement pattern (i.e., BRW or BCRW), or an indication for more complex behaviors (discussed in Supplementary note 4).

Statistics and reproducibility

Quantitative analyses are applied to directional trials, i.e., larval bearing sequences (\(\hat{\theta }\)) that are significantly different from a uniform distribution based on the Rayleigh’s test8 (p < 0.05). We compute the quantiles in which the observed precision (\({\hat{R}}_{\theta }\)) of each trial falls within the null distributions (\({R}_{{\theta }_{0}^{{vm}}}\)and \({R}_{{\theta }_{0}^{r}}\), computation explained below), and compare these quantile distributions with the null quantile distributions using a chi-square test, gauging the observed directional precision \({\hat{R}}_{\theta }\) against the potential of autocorrelated precision (\({R}_{{\theta }_{0}}\)). We employ the simplifying assumptions that a strict use of internal cues is expressed in a CRW process, and that oriented movement patterns result in more directional (or straighter) paths compared to unoriented patterns. Under these assumptions we expect that the empirical \({\hat{R}}_{\theta }\) will exceed the autocorrelated pattern \({R}_{{\theta }_{0}}\) for individuals that apply oriented movement, whereas for an unorienting individual, \({\hat{R}}_{\theta }\) is expected to be closer to \({\bar{R}}_{{\theta }_{0}}\)6 (Fig. 2c). Note however, that it is often difficult to distinguish between oriented and unoriented movement over a short duration (Fig. 2c); the differences between these two types of movements are much more apparent over long time duration, with oriented movement achieving greater displacement compared to unoriented movement45. \({\hat{R}}_{\theta }\) values less than \({R}_{{\theta }_{0}}\)may result from complex behaviors such as one-sided bias (left or right), representing the utilization of internal cues (Fig. 2c, Supplementary Information section S1). In addition, our methodology is not appropriate when movement patterns are complex, e.g., CRW composite, in which there is more than a single CRW pattern per a given sequence, and zig zag patterns, in which consecutive turning angles are drawn from von mises distributions centered around positive and negative values, respectively (Supplementary note 4).

We apply our methods to a database of 835 in situ orientation trials gathered on larvae of 21 species from eight families of perciform fishes at various tropical and warm-water locations in East Asia14,15, Australia16,17,18, Mediterranean52 and Red Sea19, synthesized from eight previously published studies (Table 1). These studies examined the orientation behavior of settlement-stage17 and pre-settlement-stage18 larvae15 of reef-14, non-reef20, and pelagic15 fish species. In addition, some of these studies examine whether larval fish use directional information from the sun for oriented movement19,52, as well as the difference in orientation patterns between individuals versus groups of larvae16.

The methodology used for these in situ studies can be divided into two main categories. First, with direct observations through Scuba-Following, where a larva is released in the pelagic environment and tracked by scuba divers for 10 min, during which swimming direction is recorded every 30 s, resulting in 21 observations (Nobs = 21); for the full protocol, see23,53. Second, with observations using the Drifting In Situ Chamber (DISC54). For each DISC trial, a larva is placed into a circular chamber, and its position is recorded for 15–20 min with an upward-looking camera fixed circa 50 cm below the chamber. The first 3–5 minutes of each DISC trial are considered as acclimation time and are excluded from the analysis, whereas the residual 10–15 min are the actual observations used for the analysis.

The number of observations per DISC deployment (Nobs = 90, 180 or 300; see Table 1) varies with the recording frequency of larval positions, ranging from 2 to 10 seconds (Table 1). Some of the DISC trials had missing data due to the fact that the position of the larva is not always identified during the manual digitization of the DISC trials, due to the small size of the larva and due to unfavorable visibility conditions54. Trials that had Nobs smaller than the required number were not used for the analyses. For the DISC trials, Nobs that had to be larger than 90% of the maximal Nobs designated per group (i.e., Nobs > 81, 162, 270). Trials with Nobs higher than the maximal Nobs were trimmed to contain the maximal Nobs per species, retaining the later-in-time data. For the scuba-following trials, the number of observations had to be Nobs > 20 due to the sensitivity of the analysis to a low number of observations. In other words, a low number of observations limits the capacity of the quantitative analyses to distinguish between oriented and unoriented movement patterns (see Supplementary note 3, Supplementary Figure S3). Importantly, both methods were shown to be robust in terms of artifacts and biases55,56, and have been tested together demonstrating high consistency in larval orientation results16,48.

Each orientation trial includes a sequence of larval swimming directions, termed bearings (θ) (Fig. 2a). For the DISC trials, θ are the cardinal directions of larval positions within the DISC’s chamber55. The angular differences between θ of consecutive time steps (t) are defined as Δθ (Δθt = θtt-1), such that for every θ sequence of a given length (N), there is a respective Δθ sequence of length N-1 (Fig. 2a). Directional precision with respect to external and internal cues is computed as the mean vector length of bearings (Rθ) and of turning angles (RΔθ), respectively54. Values of mean vector length (R) range from 0 to 1, with 0 indicating a uniform distribution of angles and 1 indicating that all angles are the same.

We used two quantitative approaches to examine if larvae exhibit oriented movement: the Correlated Random Walk- von Mises and Correlated Random Walk- wrapped Cauchy (CRW-vm and CRW-wc) analyses and the CRW resampling (CRW-r) analysis. Both types of analyses are based on the assumption that trajectories of animals that strictly use internal cues for directional movement are characterized by a CRW pattern. Hence, their capacity for directional movement is exclusively dependent on the distribution of their turning angles (Δθ)57. In contrast, for an external-cues orienting animal, for which movement directions are correlated with an external fixed direction, the mean vector length of the observed bearings, \({\hat{R}}_{\theta }\), is expected to exceed that of a CRW, \({R}_{{\theta }_{0}}\)6. Both analyses compare \({\hat{R}}_{\theta }\) against the expected \({R}_{{\theta }_{0}}\), but the first type computes \({R}_{{\theta }_{0}^{{vm}}}\)and \({R}_{{\theta }_{0}^{{wc}}}\)using theoretical von Mises and wrapped Cauchy distributions of Δθ, and the second type computes \({R}_{{\theta }_{0}^{r}}\) by producing 100 new θ sequences per individual trial (larva) by multiple resampling-without-replacement of the Δθ.

A key principle for both analyses types stems from the fact that the mean vector length of bearings (Rθ) is inherently dependent on the mean vector length of turning angles (RΔθ)28. In other words, an animal with a high capacity for unoriented directional movement, i.e., a narrow distribution of Δθ, is likely to yield a high Rθ, even if it makes absolutely no use of external cues for oriented movement. Hence, in both analyses \({\hat{R}}_{\theta }\) is gauged against a distribution of \({R}_{{\theta }_{0}}\), given its respective mean vector length of turning angles \({\hat{R}}_{\triangle \theta }\). The open-source software R58 with the package circular59 is used for all analyses in this study.

Correlated Random Walk-von Mises (CRW-vm)

In this analysis, we first generate the directional precision (R), expected for unoriented CRW movement using the theoretical von Mises distribution (\({R}_{{\theta }_{0}^{{vm}}}\)). The CRW bearings sequences (\({\theta }_{0}^{{vm}}\)) are generated by choosing a random initial bearing, followed by a series of Nobs-1 turning angles (\({\triangle \theta }_{0}^{{vm}}\)) in bearing direction; drawn at random (with replacement) from a von Mises distribution (Nrep = 1000). The length of \({\theta }_{0}^{{vm}}\) sequence is according to the number of observations in our four types of experimental trials: Nobs = 21 for the scuba-following, and 90, 180 and 300 for the DISC (Table 1). The directional precision of the von Mises distribution is dependent on the concentration parameter, kappa. Kappa values ranging from 0 to 399 are applied at 1-unit increments to cover the entire range of directional precision from completely random (kappa = 0), to highly directional (kappa = 399). Next, the directional precision of the bearings (Rθ) and the turning angles (RΔθ) are computed for each simulated sequence of θ (Fig. 2a–c).

These respective pairs of values (RΔθ, Rθ) provide the basis for generating the expected relationship between \({R}_{{\theta }_{0}^{{vm}}}\) and \({R}_{{\triangle \theta }_{0}^{{vm}}}\). Then, for any given kappa value, the following quantiles are computed: 5th, 10th, 20th,….,90th, and 95th (grey vertical distributions in Fig. 2c). Next, smooth spline functions are fitted through all respective quantiles, generating the \({R}_{{\theta }_{0}^{{vm}}}\)quantile contours, which represent the null expectation under CRW. This expected (RΔθ, Rθ) correspondence creates a phase diagram (Fig. 2c), based on which the observed θ patterns are gauged. The procedure is repeated four times to match the among-study differences in the number of θ observations per trial (i.e., Nobs = 21, 90, 180, and 300; see Table 1).

To examine if the observed larval movement patterns differ from those expected for unoriented movement (CRW-vm), we compute RΔθ and Rθ for each individual trial (\({\hat{R}}_{\triangle \theta }\) and \({\hat{R}}_{\theta }\)). We then place these values in the phase diagram and examine their positions with respect to \({R}_{{\theta }_{0}^{{vm}}}\) (Fig. 2c). Larvae with \({\hat{R}}_{\theta }\) substantially higher than \({\bar{R}}_{{\theta }_{0}^{{vm}}}\), are considered to have a higher tendency for a straighter movement than expected under CRW, suggesting oriented movement such as BRW and BCRW (Fig. 2b, c)6,28. Larvae with \({\hat{R}}_{\theta }\) values substantially below \({\bar{R}}_{{\theta }_{0}^{{vm}}}\)indicate irregular patterns such as a one-sided drift (right or left). A larva is considered directional if the bearing sequence (\(\hat{\theta }\)) is significantly different from a uniform distribution based on the Rayleigh’s test (p < 0.05)8. Non-directional larvae are characterized by low \({\hat{R}}_{\triangle \theta }\) and \({\hat{R}}_{\theta }\), and thus will be situated at the bottom left area in the phase diagram (Fig. 2c). 95% confidence interval (CI) was computed for each species’ trials (\({\hat{R}}_{\triangle \theta }\), \({\hat{R}}_{\theta }\)). The difference (∆R) between \({\hat{R}}_{\theta }\) and \({\bar{R}}_{{\theta }_{0}^{{vm}}}\) was computed per each trial and pooled by species to assess the tendency for a straight movement compared to that expected under CRW, as an indication for oriented movement.

The quantile (Q) of each trial is then computed based on the location of (\({\hat{R}}_{\triangle \theta }\), \({\hat{R}}_{\theta }\)) within the null quantiles’ contours in the phase diagram (Fig. 2c), using a 2-D interpolation such that X = RΔθ, Y = Rθ, and Z = Q (Akima R package60; Fig. 1b). 2-D interpolation is used once more to overlay the (\({\hat{R}}_{\triangle \theta }\), \({\hat{R}}_{\theta }\)) of two species with a different number of observations (Nobs = 300: Premnas biaculeatus, Nobs = 90 Chromis atripectoralis) on the DISC’s phase diagram (Nobs = 180), which represents most of the DISC trials.

Correlated Random Walk-wrapped Cauchy (CRW-wc)

Although von-Mises distribution is the most commonly used circular distribution for simulating CRW11, the wrapped Cauchy distribution well represents the underlying distributions in some cases33. To examine the sensitivity to the underlying distribution of our method, we repeated the exact same protocol of CRW-vm, only with a wrapped Cauchy distribution instead of von-Mises, and respectively, using the rho concentration parameter instead of kappa, with values (n = 400) ranging between 0 and 0.999, representing rho’s minimum and maximum values.

Correlated Random Walk- resampling (CRW-r)

In this analysis, we generate \({R}_{{\theta }_{0}}\) expected under a strict use of internal cues of CRW pattern using resampling of the turning angles (Δθ) per trial sequence (i.e., \({R}_{{\theta }_{0}^{r}}\)). Specifically, for every trial sequence, \({R}_{{\theta }_{0}^{r}}\) is computed by generating 100 θ sequences from 100 resampled Δθ sequences (Nrep = 100, without replacement) from the empirical Δθ (Fig. 2d). The \({R}_{{\theta }_{0}^{r}}\) sequence length is equal to the number of observations in each trial (Nobs). Next, Rθ for each of the resampled sequences (\({R}_{{\theta }_{0}^{r}}\)) and the quantile in which \({\hat{R}}_{\theta }\) falls within \({R}_{{\theta }_{0}^{r}}\), are computed (Fig. 2e). The quantile represents the proportion of \({R}_{{\theta }_{0}^{r}}\) which is smaller than \({\hat{R}}_{\theta }\) for each trial.

Meta-analysis chi-square goodness-of-fit tests

To test if \({\hat{R}}_{\theta }\) was significantly higher than what is expected under the null, we used chi-square tests to compare the \({\hat{R}}_{\theta }\) quantile distributions (observed trial counts) with the null (\({R}_{{\theta }_{0}^{{vm}},}\,{R}_{{\theta }_{0}^{{wc}}}\)and \({R}_{{\theta }_{0}^{r}}\)) quantile distributions (simulated sequences counts). For applying chi-square tests on larvae of different species pooled together, we used the following quantile partitioning: 0–50%, 51–70%, 71–90%, and 91–100%. For applying chi-square tests on larvae pooled by species, we used the following quantile partitioning: 0–50%, 51–75%, 76–100%. The reason for the differences is the minimal number of samples limitation of the chi-square goodness-of-fit test, which is a minimum of 5 samples per expected bin61. This limitation allows a minimum of 20 samples in the species chi-square test, and 50 samples in the chi-square test for all larvae pooled together. Importantly, the analysis is done on counts of the individual trials’ values, thus there is no information loss due to pre-analysis pooling of data.

To test whether \({\hat{R}}_{\theta }\) are significantly higher than what is expected under the null across species, we computed the effect size (Cohen’s W) of each chi-square test per species, and examined if the effect sizes distribution is significantly higher than 0.5 using a one-sided one-sample t-test, after ensuring normality of effect sizes distribution using Shapiro-Wilk test62. Effect size of Cohen’s W ≥ 0.5 represents a strong effect size for the chi-square goodness-of-fit test63. This analysis was applied to trials that contained single larva rather than groups, as grouped larvae were shown to orient differently than single larvae16.

To examine the correspondence between the \({\hat{R}}_{\theta }\) quantiles of the two methods: CRW-vm and CRW-r, the quantiles (Q) data were binned at increments of 5% for the two analyses, creating a 20 × 20 cell matrix (M). Then, the matrix was filled based on the \({\hat{R}}_{\theta }\) quantiles of the larvae, such that for example, a given larva with QCRW-vm = 99% and QCRW-r = 96%, will be counted in M20,20. Whereas, a given larva with QCRW-vm = 52% and QCRW-r = 46%, will be counted in M11,10. Based on this matrix, the corresponding heatmap contoured plot in Fig. 1e was produced, using the R package plot3D64.

If the two methods of analysis (CRW-vm and CRW-r) provide significant test results for a given species, this can be regarded as evidence for oriented movement under our simplifying assumptions. If both methods fail to reject the unoriented movement null hypothesis, it seems likely that external cues are not used for directional movement. However, if the two methods provide differing test results, no definitive conclusion about how directional movement is maintained can be reached.

Subsampling of the DISC trials was carried out, obtaining subsampled sequences with Nobs = 21. This was done to examine the effect of variation in Nobs on the analyses results (i.e., quantiles and Rθ), and to compare between the DISC and the scuba-following trials under a similar Nobs. Subsampled trials underwent the same filtering scheme as the regular trials in terms of the mean vector significance (Rayleigh’s test) and the available number of observations.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.