Taxonomic Richness and Diversity of Larval Fish Assemblages in the Oceanic Gulf of Mexico: Links to Oceanographic Conditions

Biodiversity enhances the productivity and stability of marine ecosystems and provides important ecosystem services. The aim of this study was to characterize larval fish assemblages in pelagic waters of the northern Gulf of Mexico (NGoM) and identify oceanographic conditions associated with areas of increased taxonomic richness (TF) and Shannon diversity (H’). Summer ichthyoplankton surveys were conducted in the NGoM in 2015 and 2016 using neuston net (surface layer; upper 1 m) and oblique bongo net (mixed layer; 0–100 m) tows. Over 17,000 fish larvae were collected over the two-year study, and 99 families of fish larvae were present. Catch composition in the surface layer was relatively similar to the mixed-layer catch, with carangids (jacks), scombrids (mackerels, tunas), and exocoetids (flyingfishes) being numerically dominant, while deep-pelagic species, including myctophids (lanternfishes), gonostomatids (bristlemouths), and sternoptychids (marine hatchetfishes), were present almost exclusively in the mixed layer samples. Generalized additive models were used to evaluate the effect of oceanographic conditions on ichthyoplankton abundance and biodiversity. Salinity and sea surface height (SSH) were the most influential oceanographic conditions, with higher occurrence, higher TF, and higher H’ all present in areas of lower salinity, and lower SSH. This study highlights the ecological importance of cyclonic mesoscale features and areas of upwelling as areas of increased biodiversity for larval fishes, and also indicates that the mixed layer in the NGoM is essential habitat for deep-pelagic fishes during the early life interval.


INTRODUCTION
Pelagic fishes play an important role in open ocean ecosystems, and changes in their abundances can impact community structure and ecosystem stability (Cury, 2000;Myers, 2003;Myers and Worm, 2003). Declines in the abundances of pelagic fishes are often attributed to overfishing (Ward and Myers, 2005) but other types of anthropogenic disturbance (e.g., habitat loss or degradation) and climate change also influence their distribution and abundance (Lehodey et al., 2006;Rijnsdorp et al., 2009). New management approaches that focus on ecosystem-level processes rather than single stocks or species are necessary to effectively mitigate past overexploitation and better understand the drivers of community change in pelagic ecosystems (Pikitch et al., 2004).
Research on the early life stages of pelagic fishes is important because it can provide information on spawning locations, spawning stock biomass, and population-level processes (Houde, 2002). Unfortunately, studies on larvae and juvenile fishes during the first few months of life are limited or non-existent for many open ocean species despite the fact that biological data on these stages are needed to better assess and monitor recruitment variability. Temporal and spatial trends in the distribution and abundance of fish larvae can be used to identify environmental factors that affect early life survival (Nonaka et al., 2000). Moreover, changes in the distribution, abundance, and assemblage composition can also be indicative of changing oceanographic conditions (Hernandez et al., 2010;Carassou, 2012), including anthropogenic disturbances such as the Deepwater Horizon oil spill (Rooker et al., 2013;Kitchens and Rooker, 2014). To date, information on the early life ecology and the environmental drivers of abundance of pelagic fishes in the Northern Gulf of Mexico (NGoM) is incomplete for most taxa. This is particularly true when considering the numerically dominant deep-pelagic taxa, and such information is needed to fill in data gaps regarding factors that influence the distribution, abundance, and population dynamics of pelagic species (Richardson, 2008).
The pelagic environment in particular provides unique challenges for locating areas of high biodiversity because the geographic location of mesoscale features and associated conditions are dynamic in time and space (Marchese, 2015). As a result, management of pelagic ecosystems requires multifaceted approaches that couple ecology and oceanography (Game et al., 2009;Lewison et al., 2015). Despite increased awareness of the importance of biodiversity, our understanding of biological communities in pelagic ecosystems is incomplete (Mittermeier et al., 2011). Identifying areas of high taxonomic richness (T F ) and diversity and the oceanographic conditions that create or maintain areas of elevated biodiversity are critical because species-rich ecosystems are considered more stable and less likely to collapse compared to species-poor ecosystems (Bakun, 2006;Worm et al., 2006). Increased biodiversity also has a positive impact on ecosystem services and functions, such as resource use efficiency, nutrient cycling, and higher fisheries yields, and can stabilize ecosystems against regime shifts (Gamfeldt et al., 2014;Rocha et al., 2015). As a model system, the NGoM offers many advantages for evaluating the biodiversity and structure of larval fish assemblages. Most notably, the oceanic component of this region is generally considered oligotrophic, but gets occasional injections of nutrient discharges from the Mississippi-Atchafalaya River System (MARS) that lead to higher productivity (Dagg and Breed, 2003). This supports primary and secondary production and high fishery yields of "coastal pelagic" taxa (Browder, 1993). Surrounding the MARS plume, densities of fish larvae may reach up to 20 times higher than reported for other areas of the GoM (Grimes and Finucane, 1991;Richards, 1993). In addition, the Loop Current and associated mesoscale features in the NGoM can concentrate fish eggs and larvae, particularly along fronts associated with divergent (cyclonic), and convergent (anticyclonic) eddies (Richards, 1993;Shulzitski et al., 2015). These mesoscale features play an important role in regulating the spatial distribution of ichthyoplankton (Karnauskas et al., 2013), and a higher northern intrusion of the Loop Current has been shown to increase the abundance of fish larvae in the NGoM (Lindo-Atichati, 2012). The Loop Current is generally associated with higher salinity and warmer waters. In particular, cyclonic features (cold core) often enhance production through upwelling, leading to increased foraging opportunities for fish larvae, and are thus assumed to serve as critical nursery habitat for several taxa of pelagic fishes (Richardson et al., 2010).
Here, we assess the attributes of the NGoM as early life habitat of pelagic fishes, including deep-pelagic taxa, with a special emphasis on identifying areas, and oceanographic conditions that support larval fish assemblages with high T F and Shannon diversity (H'). When determining biodiversity of the pelagic environment, it is well recognized that mesopelagic fauna (depth range: 200 to 1000 m) of both invertebrates and fishes, frequent the upper 200 m of the water column, or epipelagic zone during all life stages through diel vertical migration (Richards, 1993). In response, deep-pelagic fish taxa may be important determinants of T F and H' in the epipelagic zone, and an objective of this study was to quantify linkages between fishes typically associated with these two different zones of the water column. We also coupled T F and H' with physicochemical and biological factors using generalized additive models (GAMs) to evaluate the relative importance of oceanographic conditions on biodiversity, which provides a means for identifying regions and conditions that support species-rich assemblages of fish larvae in the NGoM. We hypothesize that biodiversity hotspots for larval fishes (high T F and H') in the NGoM will occur primarily in convergence zones (frontal features).

Sample Design
Ichthyoplankton surveys were conducted in June and July over two consecutive years (2015,2016) in a sampling corridor that ranged from 27.0-28.0 • N and 88.0-91.0 • W. The sampling corridor contained 48 stations located on transects at both 27.0 • N and 28.0 • N, with stations along each transect approximately 15 km apart (Figure 1), which represents an area sampled continuously for the past decade to assess larval recruitment variability of pelagic fishes, including billfishes, tunas, dolphinfishes, and swordfish (e.g., Rooker et al., 2012Rooker et al., , 2013Kitchens and Rooker, 2014;Cornic et al., 2018). Nearsurface sampling was conducted with a 1 × 2 m neuston net rigged with a 1200 µm mesh. Neuston net tows were conducted in the upper 1 m of the water column (surface layer) at each station, and each tow was approximately 10 min in duration. In addition, oblique bongo net tows were conducted from between 0-100 + m of the water column (mixed layer) at each station; paired bongo nets were rigged with 333-µm mesh and 500µm mesh nets. Although different mesh sizes were used for surface and mixed layer sampling, catch composition is known to be similar between the mesh sizes and gears with similar tow profiles (Richards, 1993;Randall et al., 2015), which allows for general comparisons of assemblage structure and diversity between the two distinct regions of the water column. All tows were performed at a vessel speed of approximately 2.5 knots, and the volume of water sampled during each tow was determined by equipping nets with General Oceanics flowmeters (Model 2030R, Miami, FL, United States).
Sargassum (kg/m 3 ) collected in the neuston nets and invertebrates in the neuston and combined bongo nets (kg/m 3 ) were separated, weighed, and recorded. Samples from neuston and combined 333-µm mesh and 500-µm bongo tows were preserved in a 100% ethanol solution for transport back to the lab. Sea surface temperature (SST, • C), salinity, and dissolved oxygen (mg/L) were measured at the surface of each station using a Sonde 6920 Environmental Monitoring System (YSI Inc., Yellow Springs, Ohio, United States). Other oceanographic parameters at each station were determined using remotely sensed data accessed through Copernicus Marine Environmental Monitoring Service 1 and the marine geospatial ecology toolbox (version 0.8a44) in ArcGIS (version 10.0). Sea surface height (SSH, m) data were calculated weekly at a resolution of 1/4 degree using satellite altimetry measurements (GLOBAL_ANALYSIS_PHYS_001_020) from Copernicus 2 . Distance to the Loop Current was estimated by measuring the linear distance from the edge of the feature (km), based on the 20-cm SSH contour following Randall et al. (2015) using the Spatial Analyst toolbox in ArcGIS. Water depth at each station was estimated from NOAA National Geophysical Data Center using the GEODAS US Coastal Relief Model Grid with a grid cell size of 3 arc-s 3 .
Samples from each station were sorted under Leica MZ stereomicroscope in the laboratory and fish larvae were isolated and preserved in 70% ethanol solution. All fish larvae were identified to family through visual identification following keys in Richards (2006), with family used as the taxonomic level for assessing biodiversity (Hernandez et al., 2013). Although T F and H' were estimated from identification to family level, genetic approaches such as High resolution Melting Analysis and Polymerase chain reaction were often used to determine species identification for certain taxa, which provided confirmation of assignments to the family level for all individuals assayed (Smith et al., 2009). Issues that led to identification of unknown specimens included trawl damage to the specimen and/or individuals too small to accurately identify. Damaged samples had either a significant amount of tissue missing or only part of the body was found. Individuals with a total standard length of less than 2 mm were too small to accurately identify in some cases.

Data Analysis
Two diversity measures were applied to the larval fish assemblages. Species richness (S) is commonly used to represent total number of species per sample but here we estimated T F as the number of families present in each sample. Similarly, Shannon diversity (H') was based on diversity at the family level following the equation where n is the total number of individuals and f i is the number of individuals for each family. Diversity measures T F and H' were used for statistical testing, with each station consisting of both surface layer and mixed layer samples. A two-way analysis of variance (ANOVA) was used to examine effects of station location and month-year date with separate models developed using TF and H' as dependent variables. A Bonferroni adjustment was used to account for multiple testing to decrease chances for a Type I error. Twoway ANOVAs were also used to examine inter-and intra-annual differences in both TF and H' for surface layer, mixed layer, and combined samples. Tukey's honestly significant difference (HSD) test was used to test for post hoc differences among means. All statistical analyses were run using R (version 3.4.2) with alpha set at 0.05 (Wood, 2011;Oksanen et al., 2017).
Generalized additive models were used to examine the influence of oceanographic factors for varying time periods (month, year) on T F and H'. Explanatory variables used in GAMs were month, year, SST, SSH, distance to Loop Current boundary, salinity (SAL), dissolved oxygen (DO), depth, invertebrate biomass, and Sargassum biomass. GAMs are extensions of general linear models and allow fixed effects to be modeled by 3 http://www.ngdc.noaa.gov/mgg/gdas/gd_designagrid.html using a smoothing function (Guisan et al., 2002). General GAM construction follows the equation: Where E[y] equals the expected values of the response variable (T F or H'), g represents the link function, β 0 equals the intercept, x represents one of k explanatory variables, and S k represents the smoothing function of each respective explanatory variable. In addition to oceanographic data collected at each station described earlier, remotely sensed data (SSH, distance to Loop Current) were included as explanatory variables in GAMs. Spatial autocorrelation was not deemed to be an issue for fixed stations given the dynamic nature of oceanographic conditions across our sampling corridor. A manual procedure was used to identify influential variables on TF and H' , and the final model for each diversity measure was based on minimizing the Akaike information criterion (AIC). Collinearity among variables was examined using Spearman's test and variance inflation factor (VIF), (ρ > 0.60 and VIF > 5); collinearity was not found to be an issue thus all environmental variables were tested. A manual backward stepwise selection process was used to remove explanatory variables that did not influence T F or H' based models. Stepwise selection ended when all remaining variables were significant (p > 0.05) or the AIC value started to increase when non-significant variables were removed. Percent deviance explained (DE) was calculated for each model to examine overall fit. Once the final model was selected, each variable was removed individually to see the response in AIC, and DE in order to assess the relative importance of each predictor variable following Rooker et al. (2012).

Assemblage Composition
A total of 17,091 larvae (N = 9,551 in 2015 and N = 7,540 in 2016) comprising 99 families were collected over 2 years of sampling in the NGoM ( Table 1). The top five families by percent composition in 2015 from the surface layer accounted for over 70% of the larvae collected: carangids (jacks), clupeids (herrings), exocoetids (flyingfishes), scombrids (mackerels and tunas), and istiophorids (billfishes). For the mixed layer in 2015, myctophids (lanternfishes), scombrids, carangids, gonostomatids (bristlemouths), and gobiids (gobies) were the dominant families by percent composition (Table 1). General trends in catch percent composition persisted in 2016 and numerically dominant families in the surface layer were carangids, exocoetids, scombrids, istiophorids, and hermiramphids (halfbeaks), with carangids alone accounting for nearly 40% of the larvae collected. Deep-pelagic taxa again dominated the mixed layer from 2016 with 41% of the catch comprised of myctophids, gonostomatids, and bregmacerotids (codlets) larvae. A small percentage of the fish larvae collected (∼6%) could not be positively identified because of damage or the larvae were too small. Distinct differences in the percent composition of certain families were observed between months for both surface and mixed layer samples; albeit certain taxa were consistently high in both June and July (Table 1) Table 2).
Of the 99 families collected, the percent frequency of occurrence by station of 44 families was greater than 10% in either surface or mixed layer samples in 2015 or 2016 (Table 1), representing high diversity across stations. In the surface layer, exoceotids, carangids, scombrids, and hemiramphids were relatively common and present at the majority of stations sampled across both years ( Table 1). Several families were also common to the mixed layer with percent frequency of occurrence from combined years being over 50%, including myctophids, bregmacerotids, and scombrids. Certain taxa were common in 1 year but conspicuously less common in the other year sampled, in particular scombrids and carangids.

Biodiversity: T F and H'
Taxonomic richness in the surface layer and mixed layer varied between the 2 years surveyed (ANOVA, p < 0.001), with mean T F per station being higher in 2015 (6.3 ± 2.8) than 2016 (4.6 ± 3.2) for the surface station, and similarly, with mean T F per station being higher for the mixed layer in 2015 (12.4 ± 4.6) than 2016 (10.7 ± 4.7; Figure 2). Mean T F per station was similar between June (5.9 ± 2.6) and July (5.2 ± 3.5) surveys (ANOVA, p > 0.05) for surface stations. Mean T F per station in the mixed layer was significantly higher in June (12.9 ± 4.2) than July (10.4 ± 4.9) surveys (ANOVA, p < 0.001).
Both T F and H' varied spatially in the NGoM, with the most pronounced horizontal trend occurring between the north and south sampling transects and in areas impacted by MARS (Figure 3), where salinity was lower. In general, mean T F and H' was higher along the northern transect (28.0 • N) across all months and years sampled (Figure 2). In both 2015 and 2016, the northern transect had higher mean T F (10.5 and 9.1) and  Percent total of top families of the raw data by net type, year, and month are presented.
Frontiers in Marine Science | www.frontiersin.org Responses plots from GAMs indicated that T F and H' for fish larvae in the surface layer were higher at high sea surface temperatures (>28 • C), lower sea surface heights (0.3-0.5 m), lower salinity, higher invertebrate biomass, farther from the Loop Current, and at lower Sargassum biomass (Figures 4, 5).  Table 3). Responses plots from GAMs indicated that T F and H' for fish larvae in the mixed layer were higher at SSTs above 28 • C, lower sea surface heights (0.3-0.5 m), lower salinity, higher invertebrate biomass, and farther from the Loop Current (Figures 6, 7).

DISCUSSION
Larvae of epipelagic and mesopelagic species were collected throughout our sampling corridor in both surface and mixed layer samples. Common epipelagic fishes (e.g., carangids, exocoetids, and scombrids) accounted for almost half of the fish percent composition assemblage in surface waters, and the dominance of larval taxa that inhabit the epipelagic zone as adults has also been reported in ichthyoplankton surveys of the Straits of Florida (Richardson et al., 2010), tropical Atlantic Ocean (Katsuragawa et al., 2014), and the Pacific Ocean (Vilchis et al., 2009). Densities of these taxa were markedly higher than any mesopelagic taxon collected (e.g., myctophids∼ 0.1 larvae 1000 m −3 ) in the surface layer. In contrast, mesopelagic fishes, most notably myctophids, bregmacerotids, and gonostomatids, dominated percent composition collections in the mixed layer, with myctophids alone accounting for nearly one quarter of the larval fish assemblage in the upper 100 m and present at high densities (>40 larvae 1000 m −3 ). These results are consistent with the findings that mesopelagic larval fishes dominated the catch composition in the mixed layer of other regions in the Atlantic Ocean, including the Mediterranean Sea (Alemany et al., 2010;Torres et al., 2011). Direct comparisons with other studies are limited because the majority of surveys using comparable sampling gears focused on specific taxa rather than the entire family level ichthyoplankton assemblage (e.g., Rooker et al., 2013;Kitchens and Rooker, 2014;Randall et al., 2015); however, an earlier study by Richards (1993) characterized the ichthyoplankton assemblage in the NGoM using bongo net tows to 200 m with an observed T F of 100, of which our study is highly similar (T F = 99). They also reported that myctophids, carangids, and gonostomatids were commonly collected in the upper 200 m, supporting our observation that the mixed layer represents important habitat of mesopelagic fishes during the early life period.
Mesopelagic fish larvae, particularly myctophids, bregmacerotids, and gonostomatids, were numerically dominant by percent composition in samples from the mixed layer. At night, juveniles and adults of these taxa are known to   migrate from the mesopelagic zone to the epipelagic zone (D'Elia et al., 2016); Rodriguez et al. (2006) also reported high catches of mesopelagic fishes in the epipelagic zone of the Canaries-African Coastal Transition Zone. All samples were collected during the day and their presence in the upper 100 m of the water column suggests that the earliest life stages remain in the epipelagic zone in the daytime hours, and diel vertical migration between the zones commences later (Moku et al., 2003). Several midwater taxa, including species within these three families, hatch in the epipelagic zone and begin migration as they transition from larvae to juveniles (Watanabe et al., 2002). A large fraction of the individuals collected in our surveys from these families were relatively small (<5 mm SL) with many specimens appearing to be recently hatched, likely accounting for the high numbers of mesopelagic taxa in our bongo net collections from the mixed layer. Taxonomic richness and Shannon diversity (H') varied across the sampling corridor, with estimates of both diversity measures generally higher along the northern transect (28.0 • N). It is possible, even likely, that T F and H' were higher along the northern transect because this region borders the outer continental shelf, and thus both continental shelf and oceanic species are likely present in this region, with mixed communities leading to higher diversity. Many of the families of fish larvae collected along the northern transect in this study were indicative of continental shelf assemblages (McEachran and Fechhelm, 2010), and a greater presence of continental shelf species was often found at stations impacted by freshwater inflow (green water, lower salinity, and higher turbidity). That said, the northern transect stations were seaward of the continental shelf in slope waters where fish larvae of oceanic taxa (e.g., exocoetids, istiophorids, and scombrids) are known to occur. While the northern transect was essentially a mixed assemblage of both continental shelf and oceanic taxa, nearly all of the stations in the southern transect (27.0 • N) were in oceanic waters, which explains the high abundances of exocoetids and scombrids. As a result, the larval fish assemblage was primarily comprised of oceanic species with limited contribution of continental shelf species, leading to lower overall diversity or reduced T F and H' relative to stations in the northern transect.
Assemblage diversity also varied temporally and both T F and H' were generally higher in June than July in both sampling years. In the surface layer, exocoetids, mullids, and clupeids comprised a significantly higher percentage of the assemblage in June for both years, while carangids and scombrids were higher in July. In the mixed layer, myctophids and bregmacerotids dominated the June assemblage while carangids and scombrids comprised a greater proportion of the catch in July. Temporal shifts in the abundance and assemblage composition of larval fishes are often attributed to seasonal patterns of spawning in the Gulf of Mexico (Sanvicente-Añorve et al., 1998), equatorial Atlantic Ocean (Mourato et al., 2014), and inland waters off Australia (King et al., 2016), but other factors such as the position of mesoscale features or oceanographic conditions are also known to influence presence and spatial distribution of fish larvae (Cowen et al., 2000;Randall et al., 2015;Cornic et al., 2018). Results of the present study are consistent with other studies conducted in the NGoM that indicate higher numbers of exocoetids in June (Randall et al., 2015) and higher numbers of scombrids in July (Cornic et al., 2018), with both studies attributing seasonal patterns in larval abundance to temporal variation in spawning activity. Carangids, myctophids, and bregmacerotids are known to display variable spawning throughout the year (Moku et al., 2003;Ditty et al., 2004;  Namiki et al., 2007), and this may have contributed to observed temporal shifts in the presence of certain taxa in our collections. Inter-annual variability in the abundance and diversity of larval fishes are common and often associated with temporal shifts in the location of mesoscale features (Richardson et al., 2010;Rooker et al., 2013). In 2015, a higher northward penetration of the Loop Current corresponded with higher T F and H', while the northward penetration of the Loop Current was reduced in 2016 and lower T F and H' were observed. This suggests that diversity of the larval fish assemblage in this region is related to the northward extension of the Loop Current, perhaps due to physical convergence, as the Loop Current water itself is highly oligitrophic, and these results are consistent with previous studies (Rooker et al., 2012;Cornic et al., 2018).
The intrusion of the MARS also affects the distribution and abundance of fish larvae in the NGoM (Govoni et al., 1989;Grimes and Finucane, 1991), and a primary physicochemical indicator of MARS intrusion is salinity. In the present study, salinity was an important explanatory variable in both T F and H' GAMs, indicating that assemblage diversity for larval fishes may be highly dependent on the spatial configuration of lower salinity intrusions from MARS. Freshwater discharge from MARS in the spring creates a salinity gradient in the NGoM that ranges from the river delta to the continental shelf over the summer months (Schiller et al., 2011;O'Connor et al., 2016). Stations with highest diversity of larval fishes were often found in areas with lower salinity, suggesting that areas impacted by freshwater inflow may serve as habitat for a wider range of taxa, including both continental shelf and oceanic species. We observed that both T F and H' were higher in low salinity areas because both continental shelf taxa (serranids, lutjanids, and sciaenids) and oceanic taxa (exocoetids, scombrids, and istiophorids) were often present in collections, leading to higher diversity. Generally, the MARS plume is larger in area and outflow in June relative to July as the greatest amount of freshwater is discharged in the spring (Aulenbach et al., 2007). Results from this study showed higher diversity of larval fishes in our June surveys for both 2015 and  2016, suggesting that the influx of freshwater from the MARS impacted the assemblage composition and the location of areas with higher T F and H'. Moreover, 2015 had significantly higher diversity measures than 2016, which also appears associated with the MARS plume, as there was a greater freshwater discharge in 2015 (896,600 ft 3 s −1 ) than in 2016 (539,150 ft 3 s −1 ) 4 . MARS freshwater inflow into the oceanic ecosystem is also associated with an influx of nutrients that increase primary and secondary productivity (Lohrenz et al., 1997;O'Connor et al., 2016) and likely increases food opportunities for larval fishes. Thus, areas impacted by MARS may represent favorable habitat that supports growth and survival during early life (Grimes and Finucane, 1991), which may have also contributed to higher T F and H' observed at stations influenced by MARS. This is consistent with findings that show physiochemical processes, such as salinity, have been shown to influence larval fish distribution and the interaction between larval fish and the surrounding environment ultimately determines survival and recruitment success (Fogarty et al., 1991;Bruce et al., 2001).
Spatial variability in SSH and SST were also important drivers of T F and H' in this study. GAMs indicated that diversity increased in areas with lower SSH (cold-core eddies) and mid-level water temperatures (28-30 • C). Cold-core eddies are associated with upwelling, as cold, nutrient-rich waters in these features support higher primary productivity (Biggs et al., 1997), and assemblage diversity has been shown to increase in areas of elevated productivity in both marine and terrestrial ecosystems (Waide et al., 1999;Cardinale et al., 2002). Convergent zones where two mesoscale features meet are also responsible for aggregating plankton and are therefore favorable conditions for the survival of fish larvae in the GoM (Bakun, 2006) as well as several other marginals seas including the Mediterranean Sea (Alemany et al., 2010), Caribbean Sea (Erisman et al., 2017), and Gulf of California (Avendaño-Ibarra et al., 2013), potentially leading to the increased diversity of larval fishes along these features. In addition to the fronts physically transporting larvae to convergent zones, these zones also increase feeding opportunities for larvae, leading to higher survival rates (Bakun, 2006;Acha et al., 2018;Sato et al., 2018). Results from recent studies in the NGoM of pelagic larval fishes yield similar results, with billfishes, dolphinfishes, and tunas being associated with frontal features and convergent zones (Rooker et al., 2013;Kitchens and Rooker, 2014;Cornic et al., 2018).

CONCLUSION
Baseline estimates of biodiversity and assemblage structure are critical understanding the impacts of anthropogenic stressors on marine ecosystem, and this study serves as reference point for assessing the impacts of changing conditions onlarval fish assemblages in oceanic waters of the NGoM. Biodiversity hotspots of fish larvae in the NGoM were located in areas where continental shelf and oceanic communities coalesced, with T F and H' highest along the northern transect due to 4 https://waterdata.usgs.gov/ the influence of both MARS and mesoscale oceanographic features (Loop Current), confirming that biodiversity hotspots for larval fishes (high T F and H') in the NGoM will occur primarily in convergence zones (frontal features). As a result, factors that alter physicochemical conditions (i.e., freshwater inflow linked to MARS), or the geographic position of these oceanographic features (shifts in western boundary currents due to climate change; Chen et al., 2019) will ultimately influence assemblage diversity in the NGoM, possibly leading to broader changes in ecosystem structure and stability. Additionally, our assumption that fish larvae of numerically dominant families that use the epipelagic zone as adults (istiophorids, carangids, scombrids, and exocoetids, etc.) account for the majority of ichthyoplankton in the surface layer, while the mixed layer will have a significant contribution of mesopelagic taxa, was also supported. Mesopelagic families, particularly myctophids and gonostomatids were an important component of the mixed layer assemblage, and this finding highlights the ecological connectivity that occurs between epipelagic and deep pelagic zones.

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee.

AUTHOR CONTRIBUTIONS
CM, KC-S, MC, TS, and JR contributed to the design, collection, identification, and analysis of biological samples. CM and JR wrote the manuscript. All authors contributed to the article and approved the submitted version.