Pelagic Habitat Partitioning of Late-Larval and Juvenile Tunas in the Oceanic Gulf of Mexico

Tunas are ecologically important in pelagic ecosystems, but due to their high economic value, large-bodied species are overfished. Declines in fishery landings of large-bodied tuna species in the Gulf of Mexico (GoM) are expected to increase fishing pressures on unmanaged, small-bodied tuna species, whose life history traits are less known. While predicting spawning stocks and recruitment success typically focuses on estimates of larval abundances, juveniles may provide a better estimate of future adult stock sizes, as they are more likely to survive to adulthood because mortality rates scale inversely with body size. However, distributional studies on juveniles are rare, leading to a gap in our understanding of tuna ecology. In the present study, tuna early life stages were collected across the GoM from January–September 2011. The size class examined in this study, representing large larvae and small juveniles, is larger than that of previous larval tuna studies in the GoM. Faunal composition, abundance, frequency of occurrence, and diel catchability were investigated. Generalized additive models (GAMs) were used to examine spatiotemporal distributions of the family Scombridae and the three most-abundant tuna species in the GoM’s epipelagic waters with respect to location, oceanographic features, and temporal change. In total, 11 of the 16 scombrid species inhabiting the GoM were collected, with small-bodied tuna species (Euthynnus alletteratus, Thunnus atlanticus, Auxis thazard) dominating the assemblage. Overall, scombrids were caught at higher abundances and frequencies at night than during the day, demonstrating that nighttime sampling generates a more accurate representation of faunal abundance and distribution. Abundance and presence–absence GAMs identified a coastal group (E. alletteratus and A. thazard) associated with productive continental shelf/slope environments (low salinity, higher chlorophyll a concentrations, nearer to shelf break) and an oceanic group (represented by T. atlanticus) associated with offshore, oligotrophic habitats (high salinity, lower chlorophyll a concentrations, further from shelf break). These results demonstrate that over a broad spatiotemporal domain, large larvae and juvenile tunas partition pelagic habitat on the mesoscale in addition to the temporal partitioning of adult spawning. These factors are important for spatially and temporally explicit modeling aimed at predicting tuna stock sizes.

INTRODUCTION Scombridae (i.e., the tunas, mackerels, and bonitos) are of high ecological and economic importance in pelagic ecosystems. While they are top-level predators that contribute to the pelagic food web and ecosystem structure, function, and stability (Matthews et al., 1977;Collette and Graves, 2019), they also support valuable commercial and recreational fisheries worldwide. Fisheries management and conservation efforts require information on the population dynamics of tuna early life stages in addition to spawning adults. Larval abundance indices and distribution data are used to predict spawning stock biomass, spawning patterns (location and time), spawning habitat quality, and recruitment success (Scott et al., 1993;Hsieh et al., 2006;Richardson et al., 2010). Moreover, information on larval spatiotemporal distribution and abundance provides insight into the factors that influence survival, growth, and recruitment.
Mortality rates are inversely related to body size in bony fishes; thus, the mortality rate is much higher for larvae than for juveniles (and adults;Hendriks, 1999;ICCAT, 2016a). Therefore, it is important to survey the abundance and distribution of juvenile fishes, as they represent the members of the surviving year class. High taxonomic uncertainty and limited knowledge regarding the distributional patterns of late-larval and juvenile tunas have led to an "operational taxonomic unit" gap in our understanding of tuna ecology. Thus, understanding the biological, ecological, and spatiotemporal distribution information of juvenile tunas can provide new data for fisheries management efforts and will increase our understanding of critical juvenile habitat. The lack of data on these important life history stages (larger larvae and smaller juveniles) limits adult population prediction and management.
Despite their prevalence in the GoM, small-bodied tuna species (e.g., Euthynnus alletteratus and Thunnus atlanticus) are relatively understudied , and as a result, there are currently no federal management plans or stock assessments for these fishes (ICCAT, 2016b). Limited knowledge regarding their basic ecology, biology, and distribution and abundance patterns has hindered our ability to manage small-bodied tuna species that may be heavily fished in the future and/or subjected to future large-scale Abbreviations: CW, Gulf Common Water classification; LCOW, Loop Current Origin Water classification; MIX, Intermediate Water classification; ONSAP, Offshore Nekton Sampling and Analysis Program; S EPI, minimum salinity in the epipelagic zone; T EPI , maximum temperature ( • C) recorded in the epipelagic zone.
anthropogenic disturbances, such as oil spills from increasingly deeper oil and gas extraction activities. Thus, it is essential to study late-larval and juvenile size classes in order to enhance our knowledge on these small-bodied species and their future populations in the GoM.
The GoM's highly dynamic and complex pelagic ecosystem contains hydrographic features (e.g., Loop Current, fronts, and mesoscale eddies) that can influence the development and survival of early life stages. The GoM is a semi-enclosed oceanic system that connects the Caribbean Sea to the Atlantic Ocean by the Loop Current, which transports warm water into the Gulf through the Yucatan Channel and makes an anticyclonic turn before exiting through the Straits of Florida to become the Florida Current and then the Gulf Stream (McEachran and Fechhelm, 1998). The extensions of the Loop Current have strong seasonal and annual variability, which alter the current's location, flow patterns, temperature, and hydrographic features (Molinari, 1980;Nakata et al., 2000), and in turn, affects organismal behaviors and distributions. The boundary of the Loop Current is a highly dynamic region with meanders and strong convergence and divergence zones that can generate cyclonic and anticyclonic eddies (Olson and Backus, 1985). In the northern GoM, the Mississippi River empties large quantities of nutrients into the GoM, creating a zone of high primary productivity near the river's mouth (Le Fevre, 1986;Grimes and Kingsford, 1996) that is sometimes transported offshore by interacting eddies.
Previous studies showed that specific oceanographic features provide favorable conditions for larval T. thynnus survival and recruitment success, such as moderately warm, offshore oligotrophic waters that are outside the Loop Current and corresponding eddies (Muhling et al., 2010). While most studies have focused on larval T. thynnus, mesoscale features and the freshwater inflow from the Mississippi River have also been associated with the distributions of Auxis spp., Thunnus spp., and Katsuwonus pelamis (Lang et al., 1994;Richardson et al., 2010;Lindo-Atichati et al., 2012;Habtes et al., 2014;. However, species-specific analyses remain incomplete for these taxa, as most studies only describe distributions and abundances based on genus level (e.g., Auxis spp. and Thunnus spp.). Therefore, understanding the influence of habitat parameters on the spatiotemporal dynamics of smallbodied tuna early life stages is essential for assessing their population status within the GoM.
The objectives of this study were to determine the faunal composition and assemblage structure of scombrids throughout the oceanic domain of the northern GoM and to characterize the spatiotemporal distributions of the most-abundant larval and juvenile scombrids in the GoM's epipelagic waters with respect to location, oceanographic features, and temporal change using generalized additive models (GAMs).

Sample Collection and Processing
Late-larval and juvenile scombrids were collected across the GoM during three research cruises from January to September 2011, as part of the NOAA-supported Offshore Nekton Sampling and Analysis Program (ONSAP). The ONSAP was created to assess the composition, abundance and distribution of deepwater invertebrates and fishes in the oceanic GoM that could have been impacted by the Deepwater Horizon oil spill (April-September 2010). Scombrids were collected using a 10-m 2 mouth area, 3-mm mesh Multiple Opening/Closing Net and Environmental Sensing System (MOCNESS; Wiebe et al., 1985) at a subset of established Southeast Area Monitoring and Assessment Program (SEAMAP) stations (Sutton et al., 2016;Cook et al., Unpublished).
Full details of the sampling methodology are provided in Cook et al. (Unpublished), but a brief description is as follows: the MOCNESS, a six-net, discrete-depth sampling system, surveyed specific depth strata in the water column from the surface down to 1500 m depth, with deployments centered around solar noon (day sampling) and midnight (night sampling). The depth strata were: 0-200 m (epipelagic), 200-600 m (upper mesopelagic), 600-1000 m (lower mesopelagic), 1000-1200 m (upper bathypelagic), and 1200-1500 m (upper bathypelagic). A Tsurumi-Seiki-Kosakusho (T.S.K.) magnetically sensed flowmeter was used to calculate the water volume filtered by each net; this value was then used to standardize abundances per unit effort (presented as no. individuals 10 −5 m −3 ). Samples were fixed in 10% buffered formalin:seawater at sea and later transferred to 70% ethanol:water in the laboratory.
Larval and juvenile scombrids were identified to the lowest taxonomic level possible using morphological characteristics, body shape, myomere counts, and pigmentation patterns (Richards, 2005;Pruzinsky, 2018). Late-larval and juvenile Auxis thazard were identified by the presence of a distinct lateral midline of pigmentation along the tail. In cases of taxonomic uncertainty due to the lack of larval pigmentation or the juvenile stage being morphologically undescribed (e.g., Thunnus albacares), specimens were identified to genus level only. Quality assurance/quality control were conducted with leading scombrid taxonomic experts John Lamkin (NOAA NMFS, Miami) and Aki Shiroza, M.S. (NOAA NMFS, Miami), in order to ensure the accuracy of larval identifications. Standard length (SL) measurements to the nearest 0.01 mm were taken for all specimens.

Catch Data
Although the water column was sampled from the surface to 1500 m depth, scombrid early life stages primarily inhabit epipelagic depths (Richards, 2005); therefore, statistical analyses were conducted with quantitative samples collected in the upper 200 m of the water column. Standardized abundances and percent frequency of occurrence (F o ) were calculated for each species. Standardized abundances were derived by dividing the sum of the raw count of individuals by the sum of the volume of water filtered, and F o was determined by dividing the total number of trawls in which a taxon occurred by the total number of trawls in the epipelagic. Scombrid size-frequency plots were examined to investigate variation in size classes.

Spatiotemporal Distributions in the Epipelagic: GAMs
Scombrid abundance and presence-absence GAMs were fitted using the gamlss package (Rigby and Stasinopoulos, 2005) R software (R Core Team, 2019) to examine the distributions of the family Scombridae and the three most-abundant species (Euthynnus alletteratus, Thunnus atlanticus, Auxis thazard) in relation to a suite of oceanographic, spatial, and temporal variables. GAMs allow for non-linear relationships between response and multiple explanatory variables using additive smoothing functions (Zuur et al., 2009). The variables considered for the full models were: water mass type (following Johnston et al., 2019), sea surface height anomaly (SSHA), minimum salinity in the epipelagic zone (S EPI ), sea surface chlorophyll a concentrations (Chl a), distance to the nearest 200-m isobath, maximum temperature recorded in the epipelagic (T EPI ), Julian date (since January 1, 2011), and diel cycle (day or night sampling). Water masses were identified as Gulf Common Water (CW), Loop Current Origin Water (LCOW) or an intermediate type (MIX) based on the mean recorded temperature between 200 and 600 m depth collected by in situ MOCNESS sensors (Johnston et al., 2019). S EPI and T EPI were also collected from in situ MOCNESS sensors. SSHA was derived from E.U. Copernicus Marine Service Information (CMEMS 1 ), and Chl a data were downloaded from NASA Ocean Color Group's Moderate Resolution Imaging Spectroradiometer (MODIS Aqua 2 ; Nasa Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group, 2018). Distance to the nearest 200-m isobath was calculated using the marmap package in R (Pante and Simon-Bouhet, 2013) and were derived from the General Bathymetric Chart of the Oceans (GEBCO 3 ).
Prior to each analysis, collinearity of the explanatory variables was examined using a pair-plot or pairwise scatterplot. The inclusion of variables that appeared to covary in the pair-plots were verified using variance inflation factors (VIF > 5 reflected highly correlated variables) (Zuur et al., 2010). Because the survey period covered January to September only, T EPI and Julian date were linearly correlated (Supplementary Figure 1; tau = 0.73); thus, one of the collinear variables (i.e., T EPI ) was dropped in order to create the "full model." The model did not converge when including Chl a, so it was removed from the full model as well. SSHA was also dropped from the full model, as the water mass classifications were used for simplicity.
The full model contained five explanatory variables: water mass, S EPI , distance to the nearest 200-m isobath, Julian date (2011), and diel cycle. Water masses denoted where an individual scombrid was collected based on water type classification (LCOW, CW, MIX water). S EPI was indicative of coastal runoff and riverine input. Distance to the nearest 200-m isobath (km) was considered indicative of coastal influence and geographic location. Julian date denoted intra-annual temporal change by indicating when the specimen was collected in 2011, and diel cycle was used to investigate differential catch patterns exhibited during the day and at night. Scombrid abundance data were modeled using the negative binomial distribution (NBI; Rigby and Stasinopoulos, 2005) and presence-absence data were modeled using the binomial distribution. Smoothers were fitted for distance to the nearest 200-m isobath, Julian date, and S EPI using a penalized Betasplines (pb) smoother. In each case ≤50 iterations of the Rigby and Stasinopoulos algorithm (RS method) and ≤200 iterations of the Cole and Green algorithm (CG method; Rigby and Stasinopoulos, 2005) were used to fit the models. Log e volume filtered was included in each model as an offset term to allow for differences in catch effort.
Term selection for each model was conducted by backwardselection using AICc scores. If the difference between the full model and reduced models AICc scores (dAICc) was <2, the models were considered to be equivalent and the removed variable did not affect scombrid abundance or occurrence. If the dAICc was from 2 to 4, the explanatory variable was considered to have marginally affected scombrid abundance or occurrence, and if the dAICc was >4, the explanatory variable was considered an important determinant of scombrid abundance or occurrence (Burnham and Anderson, 2002). The resulting fitted models were validated by visual examination of the quantile residuals and plotted against the observed data and against each explanatory term included in the full model.
Abundance data can be affected by encountering random, large aggregations or patches of fishes, which can skew the results, while presence-absence data comes with the cost of losing information about the actual, observed abundances. Thus, modeling both abundance and presence-absence data provided validation for the observed patterns in both models and highlighted where the models agree or disagree.

Catch Summary
A total of 326 larval and juvenile scombrids were collected from 890 quantitative tows (net fished correct depth strata with valid flow data) from the surface to 1500 m ( Table 1). The majority of scombrids (75%, n = 245) were collected in the epipelagic zone during both day (n = 100) and night (n = 107) quantitative tows ( Table 2). Scombrids were collected in 35.8% of all epipelagic trawl samples, which was more frequent than the pelagic trawl samples from the surface to 1500 m (15%).
Overall, 85.0% of individuals were identified to species, 13.5% to genus only, and 1.5% to family only. The larval and juvenile scombrid assemblage in this survey comprised 11 of the 16 scombrid species previously reported from the GoM (Richards, 2005). Only four specimens of the endangered T. thynnus were collected, all of which were caught below the epipelagic zone (>200 m depth). One T. thynnus larva (6.4 mm SL) was collected on August 25, 2011, outside of their known spawning time period.
The most-abundant species were also collected at higher F o at night than during the day in the epipelagic zone (   day and night samples, with 2.0% of the day trawls and 10.3% of the night trawls. Katsuwonus pelamis exhibited a F o in 2.0% of day trawls and 8.4% of night trawls. The remaining, rare-event species and taxa did not exhibit higher nighttime abundances or occurrences, though less than five individuals were collected from each of these species in the epipelagic zone ( Table 2).

Spatiotemporal Distributions: Generalized Additive Models
The fitted GAM modeling the total abundance of Scombridae included Julian date, water mass, diel cycle, and distance to the nearest 200-m isobath as important explanatory variables (Table 3). Overall, abundances began to increase in April and May and peaked in August (Figures 2A,C). More individuals were caught in CW (Figures 2D,F), further from the shelf break ( Figures 2G,H), and at nighttime (Figures 2I,J). The fitted GAM modeling Scombridae occurrences included Julian date as an important variable, and water mass as a marginally important variable (Table 3). Results aligned with the abundance GAMs, in which there was a higher probability of catching scombrids later in the year (Figures 2B,C) and in CW (Figures 2E,F).
Euthynnus alletteratus was only captured in CW; therefore, E. alletteratus GAMs were only fitted to samples collected Important explanatory variables (dAICc > 4) are bolded and marginally important variables (2 < dAICc < 4) are italicized and bolded. "CW only" and "LCOW/CW only" indicate when water mass types were excluded from analyses. from CW, and water mass was removed as a variable from the full models. The fitted GAM modeling the abundance of E. alletteratus contained Julian date, S EPI , and distance to the nearest 200-m isobath as important explanatory variables (Table 3). Euthynnus alletteratus abundances increased throughout the year, with the highest peak occurring around August (Figures 3A,C). Higher abundances were associated with lower S EPI (Figures 3D,F), in which catches occurred in waters as fresh as 29.47 and most specimens (n = 59) were collected in S EPI < 34. More individuals were also collected nearer to the shelf break FIGURE 3 | Term plots for the best-fitting models for Euthynnus alletteratus data fitted to the abundance data (A,D,G) and presence-absence data (B,E,H). Standardized abundance plots (C,F,I) are also presented for all important and marginally important variables. Data that had to be excluded from the model due to missing MOCNESS sensor data are indicated by crosses.
(19.79 to 227.31 km from the shelf break), with the majority of specimens within 180 km of the isobath (n = 62, Figures 3G,I).
The fitted GAM modeling E. alletteratus occurrences identified Julian date and S EPI as important, and distance to the nearest 200-m isobath as marginally important variables ( Table 3). The presence-absence model results were similar to the abundance GAMs, in which individuals were more likely to be caught later in the year (Figures 3B,C), in low S EPI (Figures 3E,F), and nearer to the shelf break (Figures 3H,I).
Only one Thunnus atlanticus specimen was collected in MIX water; thus, the GAMs fitted to T. atlanticus only included samples from LCOW and CW. The fitted GAM modeling the total abundance of T. atlanticus included Julian date and S EPI as important and diel cycle as marginally important variables ( Table 3). Thunnus atlanticus abundances showed higher abundances beginning in June and continuing through September (Figures 4A,C). High abundances were positively correlated with higher S EPI (Figures 4D,F), with catches occurring only in S EPI between 33.82 and 36.13. The majority of specimens (n = 37, 88%) were caught in water with S EPI > 35. More individuals were collected at nighttime (Figures 4G,H). The fitted GAM modeling T. atlanticus occurrences included Julian date as important and S EPI as marginally important variables ( Table 3). The presence-absence model results were similar to the abundance GAMs, in which there was a higher probability of catching T. atlanticus from June to September (Figures 4B,C) and in higher S EPI (Figures 4E,F).
Auxis thazard was only collected in CW; therefore, only samples from CW were included, and the water mass variable was removed from the full models. The fitted GAM modeling the total A. thazard abundance included diel cycle as an important variable (Table 3), with more individuals collected at night (Figures 5A,C). The fitted GAM modeling A. thazard occurrences included diel cycle as important ( Table 3) and distance to the nearest 200-m isobath as marginally important variables ( Table 3). The presence-absence model results were similar to the abundance GAMs, in which there was a higher probability of catching A. thazard at night (Figures 5B,C) and mid-shelf as an outer neritic species (Figures 5D,E).

Assemblage Structure and Spatiotemporal Distributions
The present study identified the faunal composition and assemblage structure of late-larval and juvenile scombrids throughout the northern GoM and characterized the spatiotemporal distributions of the most-abundant species in epipelagic waters with respect to geographic location, oceanographic features, and temporal change.
Of the 16 scombrid species that inhabit the GoM, 11 species were collected in the present study. The collection of a high proportion of the species that occur in the GoM FIGURE 4 | Term plots for the best-fitting models for Thunnus atlanticus data fitted to the abundance data (A,D,G) and presence-absence data (B,E). Standardized abundance plots (C,F,H) are also presented for all important and marginally important variables. Data that had to be excluded from the model due to missing MOCNESS sensor data are indicated by crosses.
FIGURE 5 | Term plots for the best-fitting models for Auxis thazard data fitted to the abundance data (A) and presence-absence data (B,D). Standardized abundance plots (C,E) are also presented for all important and marginally important variables. Data that had to be excluded from the model due to missing MOCNESS sensor data are indicated by crosses.
is likely a function of the continuous surveying method, the longevity of sampling (9-month uninterrupted sampling period, January-September 2011), and the day/night sampling during these cruise series. The advantage of continuously sampling a large spatial area and environmental conditions over a relatively long time period (9 months) increases the likelihood of collecting a variety of species. These large-scale faunal surveys remain extremely rare in pelagic research; thus, the data presented here are both robust and unique.
The dominant tunas of the Tribe Thunnini that occur in the GoM, primarily Thunnus spp., but also Auxis spp., Euthynnus alletteratus, and Katsuwonus pelamis Richardson et al., 2010;Habtes et al., 2014;, were collected in the present study. Thunnus spp. have consistently been reported as the dominant tuna taxon in ichthyoplankton surveys during the GoM's spring and summer months Lindo-Atichati et al., 2012;, although, the most-abundant taxa captured in this study were small-bodied tuna species: E. alletteratus, T. atlanticus, and Auxis thazard. It is possible that the preponderance of smallbodied taxa in our study relative to previous studies reflects higher survivorship of their early life history stages than largerbodied tunas, the latter showing a more r-selected reproductive strategy (sensu Pianka, 1970).
Scombrids collectively preferred CW conditions over the warmer and high SSHA waters within the MIX water and LCOW. Of the three most-abundant species, T. atlanticus was the only species collected in LCOW in addition to CW. Lindo-Atichati et al. (2012) found that larval E. alletteratus had similar abundances between CW and LCOW in the GoM, larval Auxis spp. (i.e., A. thazard and A. rochei) were found along the boundaries of anticyclonic features and within the GoM's CW, and Thunnus spp. and T. thynnus larvae were more abundant in the boundaries of anticyclonic features in the GoM. It has been proposed that year-round inhabitant species (e.g., T. atlanticus) have broader habitat preferences in the GoM than T. thynnus and are better able to tolerate warm features, such as the Loop Current and warm eddies (Muhling et al., 2010;Teo and Block, 2010;Lindo-Atichati et al., 2012). Variability in the abundances within CW can also relate to smaller-scale oceanographic features (e.g., cyclonic eddies) in addition to changes in environmental variables that were investigated in this study (i.e., S EPI , Chl a).
Abundance and presence-absence GAMs additionally identified a coastal group (E. alletteratus and A. thazard) and an oceanic group, represented by T. atlanticus. The coastal group was associated with more productive continental shelf and slope environments (low S EPI , high Chl a, nearer to shelf break), while the oceanic group was associated with offshore, oligotrophic habitats (high S EPI , low Chl a, further from shelf break).
Nearshore environments and areas with high terrestrial runoff are associated with low S EPI and high Chl a, which relates to the coastal life history traits of E. alletteratus and A. thazard. Euthynnus alletteratus preferred nearshore environments with lower S EPI (<34), and A. thazard was characterized as an outer neritic species, as it preferred areas along the outer shelf break. Maximum freshwater discharge into the GoM from the Mississippi River occurs in the spring. The plume is characterized by lower S EPI and higher Chl a. Their preference of nearshore environments with lower S EPI suggests that areas with high freshwater inflow are suitable habitats for the larvae and smallbodied juvenile size classes of coastal tunas. As nutrients increase from runoff, primary and secondary production increase, and in turn provide food for larvae along the continental shelf regions (Le Fevre, 1986;Grimes and Finucane, 1991;Grimes and Kingsford, 1996). Thus, riverine discharge likely maximizes growth and survival of E. alletteratus and A. thazard early life stages.
The GoM experienced a highly productive year in 2011 (Muller-Karger et al., 2015), with high Chl a waters and low S EPI , especially along the coast due to increased runoff from the Mississippi River. These favorable and highly productive conditions in the nearshore environment could contribute to the numerical dominance of E. alletteratus and A. thazard observed in this study.
Offshore waters typically exhibit higher S EPI and low to moderate Chl a. Larval and juvenile T. atlanticus were collected in offshore waters, with higher S EPI . Spawning in these open-ocean environments increases the initial survival of their eggs and larvae due to the reduction in ichthyoplankton predators compared to the coastal waters, though this may be offset at exogenous feeding by lower food supply. Similar to adult T. atlanticus, adult T. thynnus also spawn in waters with lower surface Chl a (0.10-0.16 mg m −3 ) and higher S EPI (35.5-37.0; Teo et al., 2007). Other pelagic fishes (e.g., swordfish, Xiphias gladius) also utilize warm, oligotrophic waters for spawning (Teo et al., 2007). Some tuna species have larvae adapted to living in these nutrientpoor environments, utilizing appendicularians for food at the beginning of their piscivorous early life stages ). It appears that T. atlanticus early life stages prefer to remain in, and have adapted to, living in areas with increased S EPI and decreased Chl a.  also found that T. atlanticus larvae preferred intermediate to high salinities, ranging from 31 to 36, and identified T. atlanticus as the most common true tuna from surveys conducted in the northern GoM (June and July, 2007 to 2010), accounting for 81% of the Thunnus larvae. High abundances of T. atlanticus also observed in the present study suggest that T. atlanticus is the most-abundant true tuna species inhabiting the GoM.
Temporal changes in scombrid abundances were associated with species-specific spawning preferences, in which high peaks in abundances in April, June, and August were influenced by A. thazard, E. alletteratus, and T. atlanticus, respectively.
An increase in abundance in April (around day 115) was dominated by A. thazard specimens. Auxis thazard spawns at sea surface temperatures (SSTs) of 21.6 to 30.5 • C, with mass spawning between 25.0 and 26.0 • C (Rudomiotkina, 1984;ICCAT, 2016a). The SSTs in April reached the massive spawning temperature range for A. thazard, which provides an explanation for the high observed abundances of this species.
Euthynnus alletteratus produces several spawning batches per reproductive season (Chur, 1973;Rudomiotkina, 1986;ICCAT, 2016a), which explains the numerous peaks in abundances throughout the year (increasing from March to September and a large peak in August). Additionally, spawning occurs when waters are the warmest in the GoM (preferably greater than 25 • C), from April to November. Surface temperatures reached about 25 • C in April, when spawning begins to occur for this species (Chur, 1973;Rudomiotkina, 1986;ICCAT, 2016a). Thus, the temporal changes observed through this 9-month survey influenced the high abundance of E. alletteratus in the GoM.
For T. atlanticus, spawning in the GoM typically occurs between June and September (Collette, 2010), particularly when SSTs reach 27 • C (Juarez and Frías, 1986). Water temperatures reached 27 • C in June when abundances strongly increased (with highest catches in August), aligning with the previously reported spawning preferences. A few specimens were collected in September, as temperatures began to drop. These patterns relate to the high abundances noted in the ichthyoplankton surveys conducted by  in June and July in 2011.
In addition to the three most-abundant taxa, Katsuwonus pelamis, Auxis rochei, Sarda sarda, Acanthocybium solandri, Thunnus thynnus, Scomber colias, and Thunnus albacares were collected in the present study. Some species, such as the nonresident T. thynnus, were rare in these collections, most likely due to their restricted and shorter spawning seasons in the GoM compared to other resident spawners (ICCAT, 2016a). However, it is interesting to note the collection of a 6.4 mm SL T. thynnus specimen in late August. A comparison with estimated growth rates suggests that the specimen was approximately 9 days old (Malca et al., 2017), indicating that this specimen was spawned in mid-August. While this is an observation of only a single specimen, our findings suggest that the spawning period for T. thynnus, which is currently estimated as occurring between April and June in the GoM, may extend through to August in a limited fashion. Moreover, all T. thynnus larvae were collected below 200 m depth, which indicates a level of connectivity between 'classic' epipelagic fishes and deeper pelagic waters over the course of development. Conducting additional larval surveys outside of the "typical" spawning period would help elucidate the extent to which spawning occurs at other times of year in this highly valuable species.
This study highlights the value of large-spatiotemporal scale surveys of oceanic ecosystems such as the GoM, given that scombrid species spawn at different times throughout the year and under different environmental conditions. These large-scale and long-term surveys provide information on the variance in the ecology of species that inhabit an area, can identify commonalities among the faunal assemblage overtime (among years and/or months), and can highlight typical and rare species occurrences in a region.

Size Classes
The 10-m 2 MOCNESS (3-mm mesh) primarily caught large larval and small juvenile scombrids. The average size specimen collected in this study (11.8 mm SL) was greater than other ichthyoplankton surveys in the GoM, such as Muhling et al. (2012) that collected larvae between 2.5 and 5.0 mm SL and DEEPEND's complementary ichthyoplankton cruises that collected an average-sized specimen of 4.6 mm SL .
The MOCNESS is a larger gear type compared to bongo and neuston nets (with mesh sizes ranging from 335 to 1200 µm) that are typically used in ichthyoplankton surveys to catch small larvae (Richards et al., 1984;Lindo-Atichati et al., 2012;Habtes et al., 2014;. All else being equal, avoidance is a function of mouth size and extrusion is a function of mesh size (MacLennan, 1992). The larger mouth area, mesh sizes, and faster tow speeds of the MOCNESS ostensibly reduce the ability of these larger larvae and small juveniles to avoid the gear, and in turn, these larger nets are more effective at capturing larger fishes (Kashkin and Parin, 1983). Thus, this study collected fewer planktonic larvae compared to other ichthyoplankton surveys, as smaller individuals were extruded through the MOCNESS's larger mesh, and in turn, the assemblages of larger larvae and small juveniles were represented.
The larger individuals collected in this study, existing in a higher Reynold's number environment than small larvae, have increased swimming abilities and mobility due to the development of locomotive features and the need to sustain their high energetic needs and high metabolic costs. More mobile individuals can actively locate and capture prey more easily, adjust their distributions within the water column (horizontally and vertically), and, in turn, increase their growth and survival rates (Werner and Gilliam, 1984). Understanding the assemblage structure and distributional patterns of these larger larvae and small juvenile size classes are critical for fisheries management and conservation efforts as these size classes represent the cohort that has survived the high-mortality "gauntlet" experienced by small larvae (Anonymous, 1984;ICCAT, 2016a).
Comparison of the relative dominance between sampling strategies in the same place and at similar times of year can also potentially inform taxon-specific mortality rates at different size classes. While larval assemblages are used as a proxy for predicting spawning stock biomass, late-larval and juvenile assemblages should also be used to predict future stocks, as these individuals provide a more accurate representation of the surviving year classes of scombrids in the GoM, and subsequently, those individuals (and genetic lineages) that have a higher chance of persisting to adulthood. It is an advantage to use multiple sampling methods (e.g., neuston nets, bongos, MOCNESS) in order to gain a more complete analysis of scombrid ecology and size class assemblages.

Diel Catchability in the Epipelagial
It is important to understand diel differences in day and night catch rates during surveying, as quantitative data are used in scombrid stock assessments. In this study, scombrid early life stages were collected at higher abundances and higher frequencies at night than during the daytime in the epipelagic zone. Increased catches at night are most likely a result of net detection and avoidance during the daytime . Diel differences in feeding activity also influence catch rates. Most larval and juvenile scombrids feed during the day Tanabe, 2001;Morote et al., 2008), when they are more active. Sensing the nets, and swimming away from them, is more likely to occur during the day when early life stages are more active. Thus, higher catches at night may be related to lower activity levels, decreased swimming activity, and reduced ability to visually detect and avoid nets (Takashi et al., 2006).
While daytime sampling may be appropriate for smaller larvae, it is evident that for the size classes collected in this study (larger larvae and smaller juveniles) it is more beneficial to sample at night in order to collect a more accurate representation of abundance in the epipelagic. These larger individuals are able to actively move across larger spatial scales, and their distributions are more likely to be behavioral compared to a planktonic existence exhibited by larvae.
There is additional supporting evidence of higher capture rates of larval tuna at night.  noted an increase in abundances of T. atlanticus prior to sunset and after dawn in the northern GoM, though their sampling protocol did not include night (midnight) sampling, obfuscating comparison to the present study. In the vicinity of Puerto Rico and the Virgin Islands, Hare et al. (2001) collected more K. pelamis larvae at night, but larval T. atlanticus and E. alletteratus abundances were similar between day and night tows. Previous studies involving larval Auxis spp. also observed higher catches at night in the Philippines, Atlantic Ocean, and Pacific Ocean (Wade and Bravo, 1951;Matsumoto, 1959;Strasburg, 1960;Klawe, 1963).
Recognizing catch differences between day and night sampling is important for fisheries management. Results from this study indicate that it is more appropriate to sample the epipelagic zone at night in order to collect quantitative abundance data that more accurately reflect the true abundance of large larvae and smaller juvenile scombrids in an area.

CONCLUSION
Bongo nets have proven effective at catching smaller larvae, neuston nets effectively catch slightly larger larvae, and hookand-line sampling catches adults. However, an effective sampling method does not currently exist for sampling these large larvae and small juveniles, and traditional gear types undersample these size classes. Utilizing multiple sampling methods to target tuna early life stages can improve long-term assessments of recruitment, spawning, and stock biomass of tunas. Such largescale surveys taken over seasonal cycles provide invaluable information regarding spawning, recruitment, and survivorship rates throughout a year. This study provided a cumulative quantitative analysis and more accurate representation of the scombrid cohort that survived the high mortality that is typically experienced by small larvae. It is important to continue exploring additional modifications of sampling methods for these size classes in order to gather ecological data on these poorly studied, small-bodied tuna species.
Scombrids have a wide variety of life history strategies and spatiotemporal distributions that are often dictated by adult spawning and migratory behaviors. Through spawning, adults establish the initial broad distribution of eggs and small larvae, and the larger larvae and small juveniles modify these distributions through their own behavior. Different seasonal and horizontal distributional patterns existed among the species examined in this study. Horizontal distributions were closely linked with physical characteristics of the water column and mesoscale oceanographic features. Oceanic species (Thunnus atlanticus) preferred more oligotrophic habitats (high S EPI , low Chl a, further distance from shelf break), while coastal species (Euthynnus alletteratus and Auxis thazard) preferred more productive continental shelf and slope environments (low S EPI , high Chl a, nearer to shelf break).
Overall, this study quantified the habitat preferences of latelarval and juvenile scombrids in the northern GoM. Results from this study demonstrate the partitioning of pelagic habitat by tunas, from late larvae to adults, particularly for small-bodied tuna species (e.g., E. alletteratus and T. atlanticus) that do not have any current stock assessments or management plans in place. By understanding habitat preferences of tuna early life stages, we can protect critical spawning grounds and nursery habitats and aim to improve management and conservation efforts regarding scombrid populations in the GoM.

AUTHOR CONTRIBUTIONS
NP and TS designed the study; NP analyzed the samples and organized the database; NP and RM performed the statistical analysis and prepared the figures. All authors contributed to manuscript preparation and revision, and read and approved the submitted version.