Near-global spawning strategies of large pelagic fish

Domain: earth_science
Published: 2025-09-01
Editorial outcome: accept

Abstract
Understanding the spawning strategies of large pelagic fish could provide insights into their underlying evolutionary drivers, but large-scale information on spawning remains limited. Here we leverage a near-global larval dataset of 15 large pelagic fish taxa to develop habitat suitability models and use these as a proxy for spawning grounds. Our analysis reveals considerable consistency in spawning in time and space, with 10 taxa spawning in spring/summer and 9 taxa spawning off Northwest Australia. Considering the vast ocean expanse available for spawning, these results suggest that the evolutionary benefits of co-locating spawning in terms of advantageous larval conditions outweigh the benefits of segregated spawning in terms of reduced competition and lower larval predation. Further, tropical species spawn over broad areas throughout the year, whereas more subtropical and temperate species spawn in more restricted areas and seasons. These insights into the spawning strategies of large pelagic fish could inform marine management, including through fisheries measures to protect spawners and through the placement of marine protected areas.

Introduction
Large pelagic fish, such as tuna and billfish, are key species in marine food webs1. They are important socioeconomically, supporting valuable fisheries in subtropical and tropical national waters and in the high seas2. These fisheries provide income and food for many countries, particularly island nations of the Global South2,3. Many large pelagic fish also play important ecosystem roles in carbon cycling through predation, diel vertical migration, and sinking of feces and carcasses4. Their high mobility, migrating across vast ocean expanses and traversing both national boundaries and the high seas, presents challenges for their effective fisheries management and conservation. A deeper understanding of their life history—including how these species reproduce in time and space—could ultimately aid in their management.

Life history strategies of pelagic fishes are highly variable and can be inferred from co-variation in traits, such as body size, longevity, growth rates, and fecundity. In scombrid fishes (tunas, bonitos, mackerel), most life history variation is explained by traits describing body size, swimming speed, and reproductive schedule, with body size being the most influential5. At one end, opportunistic strategists, such as tropical tunas (e.g., skipjack tuna - Katsuwonus pelamis), have a fast life history strategy. These species are often small to medium in size, grow rapidly, mature early, and have short life spans. At the other end, periodic strategists (e.g., Atlantic bluefin tuna - Thunnus thynnus) predominantly inhabiting temperate waters have a slow life history. These species are often large in size, grow slowly, mature late, and have long life spans, leading to slower population turnover5–7. Similarly, billfishes such as marlin, spearfish (Istiophoridae) and swordfish (Xiphiidae), also have a slow life history strategy because they are large bodied, late-maturing and slow-growing8.

Fundamental to the life history strategy of large pelagic fish is their spawning strategy in time and space. Tropical scombrid species with a fast life history generally have a longer spawning season than subtropical and temperate scombrids with a slow life history5. Although, it is unknown whether this trend is followed more generally by other large pelagic fish. By contrast, there is much more limited information on the extent of spawning grounds of large pelagic fish in relation to life history variation. Muhling et al.9 interestingly noted the tendency for three bluefin tuna species with slow life histories to spawn in spatially restricted spawning grounds. A synthesis of the spawning strategy in time and space of large pelagic fish could provide valuable insights into evolutionary processes governing spawning.

Many species tend to return to one or several natal grounds to spawn at similar times each year. For example, some spawning grounds host aggregations of multiple large pelagic fish species, such as off Northwest Australia for skipjack and Southern bluefin tuna, and off Japan for skipjack and Pacific bluefin tuna10–12. This co-location of spawning is despite these species having contrasting life history strategies. The extent to which this strategy of co-location of spawning in space and time is prevalent among a broader suite of large pelagic fish taxa remains unknown.

The best data on near-global spawning dynamics is restricted to key commercial tuna and billfish species, with limited information for other non-targeted species. Multiple studies have found restricted temperature preferences for large pelagic fish, suggesting co-location of spawning. Reglero et al.13 found that seven tuna species preferred spawning in warm waters and Schaefer14 found that tuna species spawn in waters >24 °C. Muhling et al.9 suggested that tuna larvae tend to be found within narrower temperature windows than adults, suggesting that there could be considerable overlap in larval distributions. Reglero et al.13 also found that most tuna spawned in waters of intermediate values of mesoscale activity from eddies, supporting the triad hypothesis of Bakun (2006)15 that links favorable spawning areas to the physical environment, and suggests that co-location might be a valuable strategy. This evidence for common environmental drivers of many species and the co-location of their spawning in time and space indicates that large pelagic fish could exploit advantageous environmental conditions for adults and larvae9,16.

There could also be substantial evolutionary benefits for large pelagic fish species if their spawning is segregated in time and space. This could reduce density-dependent food limitation and thus competition for food by adults, potentially increasing spawning and enhancing recruitment17. Such segregated spawning could also enhance the food availability for developing larvae, reduce egg and larval cannibalism by conspecifics, and decrease predation on eggs and larvae by other large pelagic fish species and invertebrates, all of which could lead to an increase in survivorship, a major factor limiting recruitment18. Therefore, both strategies—species either segregating or co-locating spawning in time and space—could have evolutionary benefits for large pelagic fish species.

A deeper understanding of the variation in life history strategies of large pelagic fish could provide the foundation for their improved management. For example, if many species co-locate their spawning in time and space, then this could be used in fisheries management to help avoid harvesting spawning females and to elucidate stock structure19,20. Detecting co-located spawning grounds, could also inform the location of marine protected areas (MPAs) or other effective area-based conservation measures (OECMs). Although debates persist on the potential benefits of MPAs for highly migratory fish species21,22 and challenges exist in protecting such species throughout their long migrations23, species such as tunas and billfish are likely to receive some benefit from the protection of their spawning areas. As countries seek to meet the Global Biodiversity Framework target of conserving 30% of land, waters and seas by 203024, knowledge of the spawning areas of commercial and non-commercial species are likely to be valuable, especially considering the lack of data available for pelagic ecosystems25.

Assessing the spawning strategies in time and space is complicated by the paucity of available data. Collecting spawning data during field surveys is labor-intensive and expensive. It usually involves sampling individual females to assess gonadal state or eggs and larvae from the water column and is sometimes complemented by local knowledge and observations26. Realized spawning areas—defined by the presence of eggs—are difficult to identify because they are patchy and ephemeral, particularly for highly migratory species27. For these species, the presence of larvae is often used as an indicator of spawning events9. Larval data and habitat suitability models can then be used to identify potential spawning times and locations28. Most data on the larvae of large pelagic fish species are collected from one-off surveys, involve individual species, cover local-to-regional scales, and use different methods29. A notable exception is the extensive dataset of Nishikawa et al.30, which provides an unprecedented window into the spawning strategies of large pelagic fish species. This dataset spans 40°N to 40°S at 1°×1° resolution, includes seasonal data collected during surveys from 1956 to 1981, and has 63,017 samples collected with similar methods. The seminal work of Reglero et al.13 on environmental requirements for spawning habitats also used the Nishikawa et al.30 observations and other datasets at a coarser resolution (5°×5°) and investigated seven tuna species. However, they did not assess spawning seasonality, nor did they address the question of the co-location of spawning by different species.

Here we leverage previously digitized30,31 larval data as a proxy for spawning to answer questions about the generality of spawning strategies—over time and space—for 15 taxa of large pelagic fish of both commercial and non-commercial value. We match these larval data with environmental variables from Earth System Model (ESM) outputs to build larval habitat suitability models and thus identify potential seasonal larval hotspots in the Indian and Pacific Oceans. To provide insights into spawning strategies, we assess the degree to which larvae for multiple species are segregated versus those that are co-located in hotspots across time and space. We investigate how consistent these spawning patterns are across hemispheres and across a slow-fast continuum of life history strategies5. Finally, we test the influence of latitude—specifically, whether the spawning season of large pelagic fish is longer5 and their spawning grounds are more restricted spatially9 at higher latitudes than those at lower latitudes. This work could help inform conservation and sustainable use of these large pelagic fish, through informed fisheries management measures and protected area placement.

Results
Larval habitat suitability
Using larval habitat suitability models for the pelagic fish taxa based on Boosted Regression Trees (BRTs, see Methods), the optimal models (Area Under the receiver operator Curve, AUC > 0.80; Supplementary Table 1) showed that high-probability areas of larval occurrence across the 15 taxa were mostly located in the tropics and subtropics of the Western Pacific and the Eastern Indian Oceans (Figs. 1–4). However, the BRTs showed notable similarities and differences in the spatio-temporal distributions of different species of larvae (Figs. 1–4). For example, habitats conducive to skipjack (Fig. 1a–d) and yellowfin tuna larvae (Fig. 1e–h) were concentrated across the Indo-Pacific, particularly in the Northwestern Pacific in boreal spring/summer, and east of Papua New Guinea and northwest of Australia in austral spring/summer. Suitable habitats for Southern bluefin tuna larvae (Fig. 2m–p) were concentrated off Northwest Australia in austral spring/summer, whereas suitable habitats for larval Pacific bluefin tuna were in the East China Sea in boreal spring/summer. Shortbill spearfish (Fig. 3e–h), albacore (Fig. 1i–l), and blue marlin (Fig. 4e–h) showed a similarly wide distribution, with larvae mostly around 20°N and 20°S in the Pacific, but also in the Indian Ocean off Northwest Australia. Suitable habitats for sailfish (Fig. 4i–l) and frigate tuna (Fig. 2a–d) were concentrated in the Indo-Pacific region and off Japan, as well as the equatorial Eastern Pacific. The sauries (Fig. 2i–l), however, were confined to the extremities of the modeled region—25–40°N in boreal winter/spring and 25–40°S west of South America in austral winter/spring.

Biophysical drivers of larval suitability
Despite spatial differences in larval habitats for each taxon, we found prominent similarities in their ranges of suitable biophysical environments, and this varied little across seasons. Most of the larvae were found in warm surface waters (Supplementary Fig. 1). The tropical tunas (e.g., skipjack and yellowfin tunas) preferred warmer waters than the temperate tunas (e.g., Southern and Pacific bluefin tunas). Most larval taxa generally preferred temperatures >25 °C, especially yellowfin and skipjack tunas, albacore, shortbill spearfish and blue marlin. However, slender tuna was generally found between 20 and 25 °C, the sauries observed in waters <25 °C, and frigate tuna had a broad temperature range, but with the greatest presence in waters >25 °C.

Probabilities of larval occurrence were higher for most taxa in waters of moderate large-scale flow activity (generally <5 m2 s−2) measured using mean kinetic energy (see Methods). However, some taxa, like frigate tuna and sailfish were found across broad ranges of mean kinetic energy (Supplementary Fig. 2). Probabilities of larval occurrence were generally higher in areas with: pH levels between 8.1 and 8.2 (Supplementary Fig. 3, except for sauries and frigate tuna found across broader ranges of pH levels); lower concentrations of chlorophyll-a (<0.2 g m−3; Supplementary Fig. 4), nitrate (<5 mmol m−3; Supplementary Fig. 5), phosphate (<0.5 mmol m−3; Supplementary Fig. 6), and ammonium (<0.4 mmol m−3; Supplementary Fig. 7); higher salinities (>30 ppt; Supplementary Fig. 8); weak thermal (<0.01 °C km−1; Supplementary Fig. 9) and salinity (<0.001 ppt km−1; Supplementary Fig. 10) gradients; and dissolved oxygen levels from 6.0 to 8.0 mg L−1 (Supplementary Fig. 11, except for frigate tuna found in broad ranges of dissolved oxygen concentrations). Probabilities were also higher in regions with surface mixed layer depths of ~50 m (Supplementary Fig. 12), in areas ~500 km from the nearest coastline (Supplementary Fig. 13), and in mean water depths of 0–6 km (Supplementary Fig. 14).

Seasonality
Mean probabilities of larval occurrence varied seasonally, but often with similar patterns in both hemispheres (Fig. 5). Three seasonal patterns were clear: spring/summer peaks, winter/spring peaks, and little seasonality. The most common seasonality was a spring/summer peak in the mean larval probability (Fig. 6) for albacore, striped marlin, and sailfish in both hemispheres. Peaks in warmer months for skipjack tuna, yellowfin tuna, bigeye tuna, and blue marlin were more evident in the northern hemisphere. Although Pacific bluefin tuna and swordfish larvae were mostly observed in the northern hemisphere and Southern bluefin tuna larvae in the southern hemisphere, all three species also followed a spring/summer pattern. Others species exhibited a winter/spring peak, including larvae of sauries, which had higher mean larval probabilities during winter/spring in both hemispheres. Frigate and slender tuna larvae exhibited the same pattern but only for the northern and southern hemispheres, respectively. Some species displayed little seasonality, such as longfin escolar and shortbill spearfish.

Spatiotemporal life-history and interhemispheric patterns
Spatial and temporal spawning patterns across hemispheres were revealed by evaluating the relationship between the Spatial Aggregation Index and Seasonal Index (see Methods) for each taxon (Fig. 6; F-statistic(25 degrees of freedom) = 11.14, p = 0.002649, r2 = 0.28, 95% confidence intervals = [0.20, 1.19]). Some taxa spawned seasonally but were spatially dispersed, such as striped marlin and swordfish. Frigate tuna exhibited spatially aggregated spawning patterns with little seasonality. Finally, spawning patterns of Pacific and Southern bluefin tunas were both spatially aggregated and seasonal. The Spatial Aggregation and Seasonality Indices were strongly related to the latitudinal preference of taxa. There was a significant positive relationship between the Spatial Aggregation Index and the mean latitude across taxa (Fig. 7a; F-statistic(25 degrees of freedom) = 11.8, p = 0.002014, r2 = 0.32, 95% confidence intervals = [0.71, 4.26]; considering each hemisphere separately for each taxon). Thus, tropical taxa (e.g., yellowfin tuna, skipjack tuna and bigeye tuna) generally had dispersed spawning areas, whereas subtropical and temperate taxa (e.g., sauries, slender tuna, Pacific bluefin tuna, and Southern bluefin tuna) had more tightly aggregated areas. There was an even stronger relationship between the Seasonality Index and mean latitude across taxa (Fig. 7b; F-statistic(25 degrees of freedom) = 32.42, p = 0.000006266, r2 = 0.56, 95% confidence intervals = [0.11, 1.10]). Thus, tropical taxa (e.g., yellowfin tuna, skipjack tuna and bigeye tuna) spawned broadly throughout the year, while subtropical and temperate taxa (e.g., Pacific bluefin tuna, slender tuna, sauries) exhibited more concentrated spawning periods. Generally, taxa found across both hemispheres had similar values for both Spatial Aggregation and Seasonality Indices in each hemisphere.

We found clear relationships between the Spatial Aggregation and Seasonality Indices and major taxonomic groups as well as with different life history strategies (Table 1). Tuna spawning areas are more spatially aggregated than billfish (Fig. 8a), whereas the degree of seasonality of tuna spawning was more variable than billfish (Fig. 8b). Further, we observed similar degrees of spatial aggregation between slower- and faster-growing taxa (Fig. 8b), but slower-growing taxa generally had greater seasonality in their spawning than faster-growing ones (Fig. 8d).

Larval hotspots
We did a Principal Component Analysis (PCA) based on the habitat suitability model outputs (i.e., probabilities of occurrences of each taxon) to identify the primary seasonal larval hotspots (see Methods). We used the first two PC axes and found similarities and differences in spatiotemporal patterns among taxa and between hemispheres (Fig. 9). Areas with high PC scores reflect areas with high probabilities of occurrences across different taxa, signifying potential larval hotspots. Also shown are the correlations of the PC axes across the taxa, where high, positive correlations correspond to taxa likely spawning in areas with positive PC scores and high, negative correlations corresponding to taxa likely spawning in areas with negative PC scores. In the northern hemisphere, there were high PC1 (Principal Component 1) scores in boreal spring/summer, particularly in the East China Sea south of Japan (Fig. 9c, e). In the southern hemisphere, high PC1 scores were found off Northwestern Australia in austral spring/summer (Fig. 9a, g). We observed generally high, positive correlations between PC1 scores and many taxa (red grids on the left in Fig. 9), indicating that these taxa could have co-located spawning in areas of high PC1 values. Thus, spawning of skipjack tuna, yellowfin tuna, blue marlin, and swordfish were consistently co-located across seasons, although swordfish did not follow the general pattern in boreal summer.

PC2 revealed spatial patterns unobserved in PC1. Positive PC2 scores revealed hotspots in the eastern tropical Pacific in boreal autumn and winter (Fig. 9b, h), north and east of Papua New Guinea in austral spring (Fig. 9h), and south of Japan in boreal winter to summer (Fig. 9b, d, f). Frigate, yellowfin, and bigeye tuna were most strongly positively correlated with PC2 across all seasons indicating that these taxa co-occur in the PC2 hotspots. Negative PC2 scores were observed in latitudinal bands at ~20°N and ~20°S across seasons in the western to central Pacific, but this extended to the Southern Indian Ocean in austral summer, which appear to be larval hotspots for taxa like albacore and shortbill spearfish across seasons, striped marlin in austral summer, and blue marlin in boreal spring.

Methods
Preparing seasonal historical larval and environmental data
A 1°×1° grid from 40°N to 40°S was created spanning the Indian and Pacific Oceans and all data were resampled into this grid. The digitized catch per unit effort (CPUE) categories from the Nishikawa data31 were converted to presence-absence data (i.e., CPUE categories >0 were considered presences) because there were few sampling points with categories >1 (Supplementary Fig. 15). Following Nishikawa et al.30, we defined seasons as one month delayed from standard calendar seasons (i.e., boreal summer is July–September and austral summer is January–March). Historical climate data for the following surface predictors were downloaded from the Ocean Model Intercomparison Project Phase 2 (OMIP2)81: (1) temperature (°C); (2) oxygen (mol m−3); (3) pH; (4) chlorophyll-a (mmol m−3); (5) salinity (ppt); (6) mixed layer thickness (m); (7) nitrate (mmol m−3); (8) phosphate (mmol m−3); (9) ammonium (mmol m−3); (10) zonal (east-west) velocity (m s−1); and (11) meridional (north–south) velocity (m s−1) (Supplementary Fig. 16–26).

OMIP2 outputs are the result of six simulation cycles of a 61-year forcing dataset (1958–2018), resulting in 366-year simulation outputs. While the Nishikawa dataset runs from 1956 to 1981, we prepared the mean, seasonal layers of each surface predictor through 1963–1981 from the last simulation cycle of the OMIP2 ensemble members (Supplementary Table 2), removing the first five years of the simulation cycle to account for the possible overshoot from the previous simulation cycle82. Using the temperature and salinity data, we calculated broad-scale thermal (°C km−1; Supplementary Fig. 27) and salinity (ppt km−1; Supplementary Fig. 28) gradients using VoCC R package83. We calculated mean kinetic energy (m2 s−2) as a proxy for large-scale surface flow dynamics84, using zonal (u) and meridional (v) geostrophic velocities (Supplementary Fig. 29):

MKE = 1/2 (u2 + v2) (1)

The mean depth layer was prepared using the current gridded, sub-ice bathymetric dataset of General Bathymetric Chart of the Oceans (GEBCO; Supplementary Fig. 30)85. The distance of each grid cell to the nearest coastline was calculated using the natural earth86 and sf87 R packages (Supplementary Fig. 31). Finally, we combined larval and environmental data to prepare seasonal datasets following the Nishikawa data (i.e., January–March, April–June, July–September, and October–December).

We intersected all layers with the grid to provide response and predictor values in each of the 15,551 cells of 1°×1°. Values of predictors for each cell were assigned based on geographic location and time period (for seasonal environmental predictors). Each grid cell inherited the nearest known value for the predictors. Larval and environmental data described were prepared with the same resolution. All analyses were done in the R statistical computing environment (version 4.4.3)88 and the Climate Data Operators (CDO) software89. Packages used to analyze the data are available in Supplementary Table 3. All output from habitat suitability models is available at Buenafe et al.90.

Building historical habitat suitability models
Many different techniques can be used to build habitat suitability models91, including Generalized Additive Models for larval habitat suitability models13,92. However, Boosted Regression Trees (BRTs)—a machine learning technique—are generally more powerful at detecting interactions between predictors79. Since our aim was to describe seasonal, historical larval hotspots, we used BRTs to build the habitat suitability models rather than other techniques. BRT models were built using the dismo R package93. Predictors in the BRTs were oceanographic and biophysical drivers that have been identified as characterizing potential spawning location of fish larvae in previous studies: (1) longitude94,95; (2) latitude94,95; (3) season13; (4) mean surface temperature28,95; (5) mean surface dissolved oxygen96; (6) mean surface pH97; (7) mean surface chlorophyll-a28,45,95; (8) mean surface salinity28,94; (9) mean mixed layer thickness95; (10) mean surface nitrate concentration98; (11) mean surface phosphate concentration98; (12) mean surface ammonium concentration98; (13) broad-scale thermal gradients28,45; (14) broad-scale salinity gradients28; (15) mean kinetic energy48; (16) mean depth95; and (17) distance to the nearest coastline95. While many of the predictors co-vary, BRTs, as do any other modeling approaches, are robust to a moderate amount of covariation79.

We used the presence-absence data (i.e., 1 – presence; 0 – absence) of the following fish taxa as the response: (1) skipjack tuna (Katsuwonus pelamis); (2) yellowfin tuna (Thunnus albacares); (3) albacore (Thunnus alalunga); (4) bigeye tuna (Thunnus obesus); (5) frigate tuna (Auxis rochei and A. thazard); (6) Southern bluefin tuna (Thunnus maccoyii); (7) Pacific bluefin tuna (Thunnus orientalis); (8) slender tuna (Allothunnus fallai); (9) blue marlin (Makaira nigricans); (10) shortbill spearfish (Tetrapturus angustirostris); (11) swordfish (Xiphias gladius); (12) striped marlin (Kajikia audax and K. albida); (13) sailfish (Istiophorus platypterus); (14) longfin escolar (Scombrolabrax heterolepis); and (15) sauries (Cololabis saira, C. adocetus, and Scomberesox saurus). We categorized these taxa according to their oceanic environment99 and life-history strategy5 (Table 1). We also reported the number of presences of each taxon from 40°N-40°S as reported in Nishikawa et al.31. All taxa had the same number of samples (n = 13,275) in the Nishikawa dataset (i.e., CPUE for each taxon was measured in all the samples). As each taxon was modeled using individual BRTs, models for taxa that were more common (i.e., larvae were present in more samples) were more robust than for species that were relatively rare (i.e., taxa that had more zeros) (Table 1). Although the larval dataset has information on three more taxa (little tuna, bonitos, and black marlin), the BRTs did not converge for these taxa because of too few presence data.

The BRT method combines regression trees with boosting, an ensemble learning method, iteratively aggregating individual weak learners to a final strong learner79. We built the models using cells with presence/absence data in our domain (n = 6,166), dividing them into an 80–20 train-test split. Using the training data, we executed a grid search to determine the configuration of hyperparameters that would yield the optimal model for each taxon (i.e., model with the highest Area Under the receiver operating Curve, or AUC, of the testing data). We ran the grid search across this range of different hyperparameters: (i) tree complexity (i.e., interaction depth): 1 or 2; (ii) learning rate: 0.005 to 0.01 in increments of 0.001; and (iii) bag fraction: 0.5 or 0.75. We reported the optimal model for each taxon by employing a five-fold cross-validation to search for the optimal number of trees while ensuring that each BRT was built with at least 1000 trees (Supplementary Table 1)79. Using the optimal models, we presented the seasonal habitat suitability maps of each taxon (i.e., 15 taxa × 4 seasons = 60 maps).

To increase the confidence in the models, we overlaid a 10°×10° grid onto the Indian and Pacific Oceans. We then only retained 1°×1° grid cells found within the larger 10°×10° cells that had at least 5% of its area covered by 1°×1° sampling cells. By doing so, we restricted the extrapolations to areas relatively close to where larvae were sampled. We also restricted the resulting habitat suitability maps of each taxon to their respective adult distributions. The range of the adult distributions are defined using predicted occurrences from AquaMaps (probabilities of species occurrence ≥0.01 are considered as presences100). Because not all taxa were identified to the species-level in Nishikawa et al.30, we followed Buenafe et al.31 and calculated the mean of the adult probabilities of species occurrence of the species listed within each taxon. Spatial dependence among grid cells was not considered in the models.

Describing interhemispheric seasonality and spatiotemporal aggregation
To visualize interhemispheric seasonality across taxa, we calculated the mean probabilities for each season in each hemisphere. As a measure of spawning seasonality, we defined the “Seasonality Index” for each taxon in each hemisphere as the standard deviation (across seasons) of the mean probabilities. Higher Seasonality Index values indicate greater seasonality in spawning. As a measure of spatial aggregation in spawning, we defined the “Spatial Aggregation Index” for each taxon in each hemisphere as the mean (across seasons) of the coefficient of variation (i.e., standard deviation over the mean) of the probabilities. The coefficient of variation scales the standard deviation with the different magnitudes of probabilities for each taxon. Higher Spatial Aggregation Index values indicate greater aggregation in spawning (i.e., less uniformly distributed). The Seasonal Index and the Spatial Aggregation Index are compared across taxonomic groups and life-history strategies (Table 1). We assessed the relationship between the Seasonality and Spatial Aggregation Indices across taxa. Further, we quantified the relationships of the indices across different predictors, such as latitude (specifically, the absolute value of latitude), taxonomic groups, and life history strategies (Table 1). These relationships were quantified by fitting linear models for the different responses (typically, the indices) and predictors (latitude, taxonomic groups, life history strategies) (Supplementary Table 4).

Identifying historical larval hotspots
To summarize the major spatial patterns in larval distribution, including the high-probability areas common to multiple taxa (i.e., larval hotspots), we performed a Principal Component Analysis (PCA) for each season on the probability of occurrence for all larval taxa. The input matrix for each season comprised the 15 larval taxa (as columns) and the resulting probabilities for each of the grid points from the BRTs (as rows). The PCA was performed on the correlation matrix of this input matrix. The first two PC axes explained 35% of the total variance (Supplementary Fig. 32), suggesting that the identified patterns explained a substantial proportion of the larval distribution patterns across taxa, although there remains considerable unexplained variation. Grid cells with high PC scores reflect areas with high probabilities of occurrences across different taxa and are thus potential larval hotspots. Information on how to perform and interpret PCA can be found in Legendre and Legendre101. We also calculated the Spearman correlation coefficient between the individual seasonal larval maps and each of the first two PCA axes to determine the relationship between each taxon and the major seasonal patterns of larval distribution.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Thank you very much for your work. All the points have been addressed with clarity and thoroughness. Congratulations. I suggest publication with no further changes.", "reviewer_3": "The authors made considerable effort to address previous reviewer comments. They provided solid and pertinent results about large pelagic fish spawning habitats, that will bring a significant contribution to the fields of marine ecology, population modelling, conservation, and fisheries management.\nI recomment this version to be accepted for publication in Nature Communications, although the authors could address the following concern:\nThe previous reviewer 2 commented on the pertinence of using Mean Kinetic Energy, computed from 1° resolution data, to describe mesoscale activity. The authors aknowledged this but did not modify the interpretation of the relationship between MKE and larvae probability of presence, referencing relationships between spawning habitats and mesoscale flow dynamics such as Elman-type upwelling and surface fronts and mesoscale eddies.\nI am not an expert on the subject, but the L378-386 paragraph would deserve more justification about the ability of MKE, calculated at this spatial scale, to describe such mesoscale dynamics.\nThe rest are minor comments:\nL81-83: Maybe reformulate to “The extent to which this strategy of co-location of spawning in space and time is prevalent among a broader suite of large pelagic fish taxa remains unknown.”\nL238: Maybe reformulate 1st sentence to:\n“The variation of spatial and temporal spawning patterns between major taxonomic groups and life history strategies was investigated using taxa-specific Spatial Aggregation Index and Seasonality Index (see Methods).”\nSection “Spatiotemporal life-history and interhemispheric patterns” L237-270:\nI would suggest first presenting the global results: taxa/hemisphere-specific SAI x SI results (Fig 5), then presenting more detailed results: relationship of indices with major taxonomic group, life history, latitudinal preference.\nL253-254:”There was a significant positive relationship between the Spatial Aggregation Index and the mean (absolute) latitude for the hemispheres combined (Fig. 4a; r2 = 0.32, p = 0.002)”\nThe method associated to this result would deserve better description (what is being compared: taxa/hemisphere-specific SAI and mean absolute latitude values, what test (Pearson’s R ?), how mean absolute latitude is computed).\nSection “Larval hotspots”\nReading this section, we don’t understand how the PCA is performed and how this approach enables to identify major seasonal hotspots.\nBefore presenting results on PC1 and PC2 scores, explain how these scores should be interpreted, that they reflect the 2 main different patterns of species co-occurence, etc.\nL335-336: “species with slower/faster life history”: change to “species with slower/faster growth”\nL371: change to “while sauries and slender tuna prefer cooler temperature”"}, "author_response": "Response to reviewers’ comments\nWe thank the Reviewers and the Editor for your comments and feedback on the manuscript. We really appreciate the time and effort you’ve put into this work. Addressing the comments have greatly improved the results and the narrative of the paper.\nReviewer #1 (Remarks to the Author):\nThank you very much for your work. All the points have been addressed with clarity and thoroughness. Congratulations. I suggest publication with no further changes.\nThank you for your positive comments! Your feedback has been invaluable in improving the manuscript.\nReviewer #3 (Remarks to the Author):\nThe authors made considerable effort to address previous reviewer comments. They provided solid and pertinent results about large pelagic fish spawning habitats, that will bring a significant contribution to the fields of marine ecology, population modelling, conservation, and fisheries management.\nI recomment this version to be accepted for publication in Nature Communications, although the authors could address the following concern:\nThank you! They are very kind words. We have addressed your comments below and have made the necessary changes in the manuscript.\nThe previous reviewer 2 commented on the pertinence of using Mean Kinetic Energy, computed from 1° resolution data, to describe mesoscale activity. The authors aknowledged this but did not modify the interpretation of the relationship between MKE and larvae probability of presence, referencing relationships between spawning habitats and mesoscale flow dynamics such as Elman-type upwelling and surface fronts and mesoscale eddies.\nI am not an expert on the subject, but the L378-386 paragraph would deserve more justification about the ability of MKE, calculated at this spatial scale, to describe such mesoscale dynamics.\nThanks for picking this up. We’ve added a few more sentences in this paragraph to provide more context and justification on using MKE:\n“AddiTonally, the probability of larval occurrence was highest in areas with moderate mean kineTc energy values, characterized by intermediate large-scale dynamic flow. This is consistent with the opTmal environmental window hypothesis, which suggests spawning habitats are found in areas of moderate ocean mixing15,46. In areas of low kineTc energy, and thus low turbulence, limited mixing leads to low producTvity and encounter rates between larvae and their prey, leading to poor feeding condiTons and recruitment13,47. Conversely, in areas of high kineTc energy and thus turbulent condiTons, there is a greater mixing of nutrients into the euphoTc zone, but the lack of stability causes larvae to struggle to capture prey in these turbulent environments13,47,48. However, at intermediate kineTc energies, there is moderate mixing and sufficient nutrients to promote primary producTvity, and larval encounter rates are sufficiently high and successful for favorable feeding condiTons and strong recruitment. Further, strong advecTon associated with high kineTc energy can reduce recruitment through larval advecTon into poor feeding and recruitment areas48,49. Thus, larval survival is generally higher in areas characterized by moderate flow15, but can also be reduced by predaTon, as predators can also aggregate in these areas50.”\nWe have also added this in the Caveats:\n“Fi[h, the spaTal resoluTon of the calculated mean kineTc energy is 1°, which might not adequately reflect finer-scale local mixing processes. Thus, relaTonships between species occurrence and mean kineTc energy should be interpreted with cauTon.”\nThe rest are minor comments:\nL81-83: Maybe reformulate to “The extent to which this strategy of co-location of spawning in space and time is prevalent among a broader suite of large pelagic fish taxa remains unknown.”\nChanged to the suggested text.\nL238: Maybe reformulate 1st sentence to:\n“The variation of spatial and temporal spawning patterns between major taxonomic groups and life history strategies was investigated using taxa-specific Spatial Aggregation Index and Seasonality Index (see Methods).”\nWe have adopted the suggested text. However, given your following comment below, we have moved this to the second paragraph.\nSection “Spatiotemporal life-history and interhemispheric patterns” L237-270:\nI would suggest first presenting the global results: taxa/hemisphere-specific SAI x SI results (Fig 5), then presenting more detailed results: relationship of indices with major taxonomic group, life history, latitudinal preference.\nThank you. We have now restructured the Spatiotemporal interhemispheric and life-history patterns. We have also re-ordered the Figures:\n“SpaTal and temporal spawning paderns across hemispheres were revealed by evaluaTng the relaTonship between the SpaTal AggregaTon Index and Seasonal Index (see Methods) for each taxon (Fig. 6; F-staTsTc (25 degrees of freedom) = 11.14, p = 0.002649, r2 = 0.28, 95% confidence intervals = [0.20, 1.19]). Some taxa spawned seasonally but were spaTally dispersed, such as striped marlin and swordfish. Frigate tuna exhibited spaTally aggregated spawning paderns with lidle seasonality. Finally, spawning paderns of Pacific and Southern bluefin tunas were both spaTally aggregated and seasonal. The SpaTal AggregaTon and Seasonality Indices were strongly related to the laTtudinal preference of taxa. There was a significant posiTve relaTonship between the SpaTal AggregaTon Index and the mean laTtude across taxa (Fig. 7a; F-staTsTc (25 degrees of freedom) = 11.8, p = 0.002014, r2 = 0.32, 95% confidence intervals = [0.71, 4.26]; considering each hemisphere separately for each taxon). Thus, tropical taxa (e.g., yellowfin tuna, skipjack tuna and bigeye tuna) generally had dispersed spawning areas, whereas subtropical and temperate taxa (e.g., sauries, slender tuna, Pacific bluefin tuna, and Southern bluefin tuna) had more Tghtly aggregated areas. There was an even stronger relaTonship between the Seasonality Index and mean laTtude across taxa (Fig. 7b; F-staTsTc (25 degrees of freedom) = 32.42, p = 0.000006266, r2 = 0.56, 95% confidence intervals = [0.11, 1.10]). Thus, tropical taxa (e.g., yellowfin tuna, skipjack tuna and bigeye tuna) spawned broadly throughout the year, while subtropical and temperate taxa (e.g., Pacific bluefin tuna, slender tuna, sauries) exhibited more concentrated spawning periods. Generally, taxa found across both hemispheres had similar values for both SpaTal AggregaTon and Seasonality Indices in each hemisphere.\nWe found clear relaTonships between the SpaTal AggregaTon and Seasonality Indices and major taxonomic groups as well as with different life history strategies (Table 1). Tuna spawning areas are more spaTally aggregated than billfish (Fig. 8a), whereas the degree of seasonality of tuna spawning was more variable than billfish (Fig. 8b). Further, we observed similar degrees of spaTal aggregaTon between slower- and faster-growing taxa (Fig. 8b), but slower-growing taxa generally had greater seasonality in their spawning than faster-growing ones (Fig. 8d).”\nL253-254:”There was a significant positive relationship between the Spatial Aggregation Index and the mean (absolute) latitude for the hemispheres combined (Fig. 4a; r2 = 0.32, p = 0.002)”\nThe method associated to this result would deserve better description (what is being compared: taxa/hemisphere-specific SAI and mean absolute latitude values, what test (Pearson’s R ?), how mean absolute latitude is computed).\nWe have now rephrased this sentence, removing “(absolute)” to be clearer:\n“There was a significant posiTve relaTonship between the SpaTal AggregaTon Index and the mean laTtude across taxa (Fig. 7a; F-staTsTc (25 degrees of freedom) = 11.8, p = 0.002014, r2 = 0.32, 95% confidence intervals = [0.71, 4.26]; considering each hemisphere separately for each taxon).”\nWe have also added more information on these comparisons in the Methods:\n“We assessed the relaTonship between the Seasonality and SpaTal AggregaTon Indices across taxa. Further, we quanTfied the relaTonships of the indices across different predictors, such as laTtude (specifically, the absolute value of laTtude), taxonomic groups, and life history strategies (Table 1). These relaTonships were quanTfied by fipng linear models for the different responses (typically, the indices) and predictors (laTtude, taxonomic groups, life history strategies) (Table S3).”\nSection “Larval hotspots”\nReading this section, we don’t understand how the PCA is performed and how this approach enables to identify major seasonal hotspots.\nBefore presenting results on PC1 and PC2 scores, explain how these scores should be interpreted, that they reflect the 2 main different patterns of species co-occurence, etc.\nWe have added some more information on how we used PCA to identify larval hotspots in the methods:\n“To summarize the major spaTal paderns in larval distribuTon, including the high-probability areas common to mulTple taxa (i.e., larval hotspots), we performed a Principal Component Analysis (PCA) for each season on the probability of occurrence for all larval taxa. The input matrix for each season comprised the 15 larval taxa (as columns) and the resulTng probabiliTes for each of the grid points from the BRTs (as rows). The PCA was performed on the correlaTon matrix of this input matrix. The first two PC axes explained 35% of the total variance (Fig. S32), suggesTng that the idenTfied paderns explained a substanTal proporTon of the larval distribuTon paderns across taxa, although there remains a considerable unexplained variaTon. Grid cells with high PC scores reflect areas with high probabiliTes of occurrences across different taxa and are thus potenTal larval hotspots. InformaTon on how to perform and interpret PCA can be found in Legendre and Legendre (2012)100. We also calculated the Spearman correlaTon coefficient between the individual seasonal larval maps and each of the first two PCA axes to determine the relaTonship between each taxon and the major seasonal paderns of larval distribuTon.”\nWe have added more information in the “Larval hotspots” section in the Results:\n“We did a Principal Component Analysis (PCA) based on the habitat suitability model outputs (i.e."}]

Decision letter