Predictive models of the preferential distribution of demersal fish larvae in the southern part of the California Current ; Modelos predictivos de la distribución preferencial de larvas de peces demersales en la parte sur de la corriente de California

1 Instituto Politécnico Nacional-Centro Interdisciplinario de Ciencias Marinas, Apdo. Postal 592, La Paz, B.C.S. 23096, Mexico. (MTP-A) E-mail: tere_peiro@hotmail.com. ORCID iD: https://orcid.org/0000-0003-3550-8456 (RF-R) (Corresponding author) E-mail: rfunes@ipn.mx. ORCID iD: https://orcid.org/0000-0002-6139-0861 (RG-A) E-mail: rarmas@ipn.mx. ORCID iD: https://orcid.org/0000-0003-4027-6143 (SHT) E-mail: strujil@ipn.mx. ORCID iD: https://orcid.org/0000-0002-2533-3915 2 Universidad Autónoma de Baja California-Facultad de Ciencias Marinas, Apdo. Postal 453. Ensenada, B.C. 22860, Mexico. (RD) E-mail: rdurazo@uabc.edu.mx. ORCID iD: https://orcid.org/0000-0002-3290-8633 3 CONACYT-Centro de Investigaciones Biológicas del Noroeste, La Paz, B.C.S. 23096, México. (ROM-R) E-mail: raul.martinez.rincon@gmail.com. ORCID iD: https://orcid.org/0000-0002-1204-716X


INTRODUCTION
In ocean ecosystems, habitat characterization provides predictive information on the spawning habitat and distribution of species, which is important for understanding the factors that regulate spawning and could be used for assessing ongoing changes in the pelagic ecosystem (Ibaibarriaga et al. 2007, Bernal et al. 2007, Bautista-Romero et al. 2018. The preferred distribution reflects positive habitat selection either by breeders or through higher survival of eggs and larvae (Asch and Checkley, 2013). In this regard, knowledge of the three-dimensional characteristics of the water column (physical process and productive waters) is particularly important for predicting spawning and distribution of some species (Weber andMcClatchie 2010, Lynn 2003). Additionally, fish assemblages vary according to midwater oxygen content, sea surface temperature and advection by the California Current, which in many cases covary with basin-scale indices (Koslow et al. 2013). Hence, determining habitat quality is the key aspect for explaining recruitment variability (Ibaibarriaga et al. 2007, Kitchens andRooker, 2014).
The habitat characterization of fish larvae facilitates the identification of the key set of variables governing the permanence and distribution of organisms. It is known that in demersal fish species, larval abundance responds to seasonal variability of environmental forcing, which is closely related to reproductive periods, geographic location and potential habitat for the expansion or shrinkage of the distribution range (Bailey 1981, Funes-Rodríguez et al. 2009, Peiro-Alcantar et al. 2016). On the western coast of the Baja California Peninsula, larval fluctuations are related to the predominance of sub-Arctic water in winter-spring and part of summer, in addition to the warming effect of surface and subsurface counter-currents from late summer to autumn to the south of the peninsula (Durazo and Baumgartner 2002, Durazo 2015, Valle-Rodríguez and Trasviña-Castro 2017. Additionally, other oceanic processes that promote changes in thermohaline structure and hydrostatic stability occur across spatial scales. These include coastal upwelling, eddies, jets, filaments, fronts, and El Niño Southern Oscillation (ENSO), all of which reflect the transitional nature of water masses (Roden 1991, Lynn and Simpson 1987, Durazo 2015. The demersal fish larval community in the southern portion of the California Current system is highly diverse and splits into two larval assemblages around Punta Eugenia and Bahía Sebastián Vizcaíno (28°N) (Peiro-Alcantar et al. 2016). One assemblage congregates to the north, characterized by the predominance of species of sub-Arctic−transitional affinity (winter-spring). The other appears to the south, com-posed mostly of species of subtropical affinity during summer-autumn (Funes-Rodríguez et al. 2009, Peiro-Alcantar et al. 2016, including some species of fishing importance (e.g., Merluccidae, Scorpaenidae, Pleuronectidae) The objectives of this study are to characterize the habitat of the dominant demersal fish larvae across the distribution range from data collected in the southern part of the California Current (1997-2000 and 2006-2010) and to assess their seasonal response to environmental variables. Knowing how environment and ecological patterns regulate these populations will allow a better understanding to predict the optimal and limiting factors across their geographic distribution range.
The sampling plan included 85 stations in waters off the Baja California Peninsula as part of the Mexican California Current Research Program (IMECOCAL), similar to the arrangement of stations of the California Cooperative Oceanic Fisheries Investigations (Cal-COFI) program in Mexican waters. The distance between stations was ~35 km along transects that were set perpendicular to the coastline, separated ~70 km from each other (Fig. 1). At each station, vertical casts were conducted using a factory-calibrated CTD (Sea-Bird) fitted with dual conductivity and temperature sensors. The explanatory variables included sea surface temperature and salinity at 50 m depth, absolute dynamic height (ADT; dyn cm; http://marine.copernicus.eu), mixed layer depth, water column stratification (φ, J m 3 ; Simpson et al. 1990) at a 0-50 m range across the water column, and spatial and temporal information [geographic coordinates (latitude-longitude) and season]. We used ADT data because these allow a wider view of the surface ocean circulation during each cruise than the spatial scope of in situ dynamic height calculations; however, the two were similar.
Plankton samples were collected with Bongo nets (0.61 m mouth diameter; 505 μm mesh size) fitted with digital flowmeters to measure filtered water volume. Nets were towed obliquely to the surface from a maximum depth of ~200 m. Samples were preserved in 4% borate-buffered formaldehyde in seawater. Ichthyoplankton was sorted from plankton samples and larvae were identified to species level (Moser 1996). Fish larval abundances were standardized to larvae per 10 m 2 of marine sea surface for comparison with other studies in the California Current; detailed descriptions of gear, sampling procedures, and standardization are given in Smith and Richardson (1977). The most frequent and abundant demersal taxa were selected to characterize their habitat (Peiro-Alcantar et al. 2016), including species subjected to large-scale commercial exploitation and recreational fishing (Merluccius productus, Sebastes sp.), and others of lesser commercial interest (Citharichthys stigmaeus, C. xanthostigma, Symphurus atricaudus and Synodus lucioceps). Sebastes is a speciose genus characteristic of sub-Arctic to warm temperate zones. The Sebastes sp. larvae used here belong to the most frequent and abundant morphotype; their larvae were somewhat more abundant off central California and declined south of the Southern California Bight (Smith and Moser 2003).
Generalized additive models (GAMs) were used to investigate the influence of spatial, temporal and environmental factors on demersal fish larvae on the western coast of the Baja California Peninsula. The use of GAMs in fish habitat modelling is valuable for exploring species-environment relationships, seeking to predict the habitat preference of adults in the reproduction stage and the potential habitat of larvae (Weber andMcClatchie 2010, Kitchens andRooker 2014). All GAMs were performed in R (R Core Team 2017) using the mgcv package (Wood 2011). Since response variables are abundance data (continuous and positive values), the fitted GAMs used the Gaussian distribution family with log as a link function. The log link function was used to reduce the high overdispersion of data. The level of smoothing (i.e. the number of bases (k) in smoothing functions) for each predictor variable was restricted to 4 to avoid overfitting. Collinearity among variables was examined using the Pearson correlation coefficient, so highly correlated variables (R>0.7) were removed from the models. Since temperature at 50 m depth and dynamic height were highly correlated (R>0.7), only one of these two variables was used in the models, depending on its contribution to the model. Thus, we used up to seven predictor variables with a low correlation (Table 1) in the final models. This technique has been used successfully to describe species-environment relationships and to predict habitat preferences or potential habitat distribution (Austin 2007, Elith andLeathwick 2009). The mathematical representation of the GAM applied to abundance of selected species can be expressed as follows: log (abundance) = α + f 1 X 1 + f 2 X 2 + f 3 X 3 + f 4 X 4 + + f 5 X 5 + f 6 X 6 + X 7 where α and X n denote intercept, longitude and latitude (interaction term), temperature, dynamic height, salinity, mixed layer depth, water column stratification and season (as a parametric term), respectively; and f n are smoothing functions (thin plate regression spline) (Wood 2011).

Surface ocean circulation
Seasonal climatology maps of the ADT (1993ADT ( -2013 show the cyclonic eddy of Southern California clearly visible in winter, summer and autumn ( Fig. 2A, C, D). There is also an increase in the number of mes-  oscale structures related to an anticyclonic flow, with an elongated meander being more pronounced during summer and autumn (Fig. 2D). Sea surface circulation is mostly towards the equator during spring and summer. This flow, which corresponds to the California Current, is driven by the prevailing southward winds. In spring, when winds along the coast are stronger, upwelled waters near the coast have lower temperature and consequently display the lowest elevation along the entire peninsula (Fig. 2B). Upwelled water is evident in the sea level minimums (ADT) and is reduced significantly in summer and autumn, especially in the region south of 28°N (Fig. 2C, D). Also during spring, mesoscale structures are reduced and the cyclonic gyre of Southern California is non-existent.
In the southern region, a coastal flow towards the pole is evident in summer and increases in autumn (Fig. 2C, D), as indicated by counter-currents. During summer-autumn, connectivity may be established between waters of tropical origin from the southern margin of the peninsula and the Pacific anticyclonic gyre (subtropical waters) from the eastern margin, as indicated by current lines (Fig. 2C, D).

Habitat characterization
Performances of the best-fitted models are shown in Table 2. Effect plots of the best-fitted models (Fig.  3) were used to describe the effect of environmental, temporal and spatial factors on selected demersal  fish larvae on the west coast of the Baja California Peninsula.
Preference in winter-spring M. productus is expected to occur principally in spring in coastal and offshore waters in the northern area (29-31°N), associated with lower dynamic height (<0.9 m) and lower salinity (<33.5 at 50 m depth), in relatively unstratified waters (Phi<50 m) over a wide range of mixed layer depths (40-100 m). The higher abundance of Sebastes sp. is predicted to occur during the spring in coastal and offshore waters in the northern area, but the model suggests that Bahía de Sebastian Vizcaíno (29°N) is also important. This species is estimated to occur at lower temperatures (<16°C at 50 m depth) and salinity (33.4-34.0 at 50 m depth) at a shallow mixed layer depth (<60 m). A higher abundance of C. stigmaeus occurs mainly in the northern area (29-31°N), associated with lower temperature and salinity (<16°C and <33.5, respectively), in areas with a shallow mixed layer (<20 m) and unstratified water (Phi<30 m).
In general, M. productus and Sebastes sp. showed a preference for coastal and oceanic waters of the northernmost region in spring. These species of northern distribution, including C. stigmaeus, were recorded in waters where temperature and salinity are lower, in relatively unstratified waters and at a shallow mixed layer. The exception is M. productus, which was found over a wide range of mixed layer depths, and is expected to occur when the dynamic height is lower. This suggests a preference of these species for northern coastal areas, where the winds along the coast are stronger and unstratified waters are consequently related to newly emerged waters along the peninsula, during the main flow of the California Current sea surface circulation towards the equator.

Preference in summer-autumn
A higher abundance of C. xanthostigma is expected to occur in the northern zone, with a core of peak abundance just off Punta Baja (30°N). Although dynamic height has little effect on this species, it is predicted to occur when dynamic height is higher (>0.95 m), associated with lower salinity (<33.5) and unstratified waters (Phi<30 m) in areas with shallow mixed layers (<20 m). In the case of S. atricaudus, the highest abundance occurs in coastal waters across the study area during the autumn. The model suggests a preference for lower dynamic height (<0.80 m) in areas with shallow mixed layer (<40 m), with low-to-high water column stratification (<20 m and >100 m depth) and salinities (33.0 and 34.0). Higher abundances of S. lucioceps should occur during the autumn in coastal waters centred on the wide continental shelf, mainly in Bahía Sebastián Vizcaíno (28°N) and to the south. This species is associated with higher dynamic height (<1.0 m) and salinity (>34.0) in unstratified waters (Phi<30 m).
These species of widespread distribution across the study area are centred on the wide continental shelf in Bahía Sebastián Vizcaíno (28°N) and to the south, except for C. xanthostigma, which occurs mainly in the northern zone. Although C. xanthostigma is estimated to occur at relatively lower salinities (contrasting with S. atricaudus and S. lucioceps, which are reported over a relatively broader range), none of the species was related to temperature. Except for S. atricaudus, these species are expected to occur associated with peak sea levels in relation to the countercurrent flow south of 28°N during summer and autumn.

DISCUSSION
The preferred habitat distribution of fish larvae was analysed with the use of GAMs in order to identify the most significant environmental and geographic predictors in the southern waters of the California Current. The analysis allowed a set of variables to be identified according to habitat requirements that are usually specific, rather than the influence of a single variable. In this regard, it has been pointed out that combinations of ocean variables may directly or indirectly serve as predictors of spawning areas (Weber and McClatchie 2010, Asch and Checkley 2013, Checkley et al. 2017.

Preference in winter-spring
M. productus, Sebastes sp., and C. xanthostigma belong to the so-called "Northern Complex", distributed in the north-south plane across the area of influence of the California Current, but probably close to the coast when their distribution reaches the southern waters off the peninsula (Moser et al. 1987). C. stigmaeus and C. xanthostigma display a particular affinity for warm water, but are included in the "Northern Complex" due to the high affinity with the northern anchovy, Engraulis mordax; C. stigmaeus is the most temperate taxon in this group (Moser et al. 1987). The distribution of the slope/transitional larval fish assemblage is primarily related to the bathymetric range of adults.
Larvae of M. productus and Sebastes sp. show the Southern California Bight as their preferential distribution centre, while larvae of Sebastes sp. and Citharichtys stigmaeus are distributed across the California Current , Bautista-Romero et al. 2018. In contrast, C. xanthostigma, Synodus lucioceps, and Symphurus atricaudus are associated with fauna in warm to temperate areas (Moser et al. 1987). Their distribution is limited to the Baja California peninsula, mainly in the vicinity of Punta Eugenia (28°N), although Symphurus spp. and C. xanthostigma also occur with a low frequency in southern California .
The highest larval abundances of M. productus (Bailey et al. 1981, Funes-Rodríguez et al. 2009) and Sebastes spp. (Moser andBoehlert 1991, Moser et al. 1993) are consistent with those reported for coastal and slope waters north of the peninsula, closely related to the habitat used by adults and larvae. A second increase in M. productus to the south of the peninsula was also observed, evidenced as latitudinal differences in abundance, which seemed to break at Punta Eugenia (28°N). This distribution is probably related to the stock of dwarf hake, similar to the findings reported by Funes-Rodríguez et al. (2009). Additionally, the abundance of M. productus decreases off Punta Baja (~30°N), apparently bordered by waters transported towards the coast from the outer margin of the California Current around Isla Guadalupe. This finding could be related to an anticyclonic flow around an elongated meander, such as the one described by Durazo and Baumgartner (2002) and Soto-Mardones et al. (2004). These observations also explain the transport of mesopelagic larvae to the coast (Funes-Rodríguez et al. 2011)  M. productus and Sebastes sp. occur mainly in winter-spring and are scarce in other seasons of the year (Moser and Boehlert 1991, Funes-Rodríguez et al. 2009). In this study, their larvae occurred principally in spring under sub-Arctic waters characteristic of the California Current. The current is defined by more homogeneous flow off Baja California, intense winds that favour coastal upwelling, low temperature and salinity from the surface to 100 m depth off California, and maximum coast-to-ocean slope at the sea surface, mainly in spring (Lynn and Simpson 1987, Soto-Mardones et al. 2004, Durazo 2009, 2015. The most significant environmental predictors of both species were similar and defined mostly by the seasonal variation in the strength of the California Current southward flow along the west coast of the Baja California Peninsula. Although larvae of M. productus are positively related to relatively lower SST, other studies (Bjorkstedt et al. 2002, Bautista-Romero et al. 2018) recorded a preference for low salinity, low dynamic height and unstratified water in a wide mixed layer depth. This preference explains the broad distribution of M. productus larvae in waters south of the peninsula. Similarly, larvae of Sebastes sp. are related to unstratified waters and low temperatures, these parameters being the most significant predictors of its larval distribution, as evidenced in localities with low SST, high geostrophic flow or lower ADT values in other studies (Weber and McClatchie 2010, Ash and Checkley 2013, Kosllow et al. 2013. This would explain the relationship with other species that take advantage of inshore areas in spring, characterized by highly productive and relatively stable conditions in the Southern California Bight (Weber and McClatchie 2010).
Upwelling areas may exert major effects in the planktonic stages, derived from the advection of eggs and larvae outside of breeding areas and consequently affecting the recruitment of stocks (Bailey 1981, Sakuma andRalston 1997). Similarly, the scarcity of larvae of coastal species is associated with the region of upwelling and offshore Ekman transport, suggesting that adults avoid spawning in upwelling areas, thus avoiding offshore transport of their larvae into the California Current (Doyle et al. 1993). Although coastal species with benthic eggs tend to be relatively scarce during the upwelling season in Angola Bay, three species of flatfishes showed positive responses to upwelling events. These were associated with the highest microzooplankton biomass, suggesting that the spawning of these species is timed to maximize food availability (Pattrick and Strydom 2014).
To compensate for the effect of the offshore Ekman transport, the spawning season in M. productus and Sebastes sp. occurs in early winter in California ), when offshore transport shows a seasonal decrease (36-27°N, Bailey 1981). Additionally, most late eggs and larvae of M. productus are located near or within the thermocline (Ahlstrom 1959, Cass-Calay 1997, Moser et al. 1997) limited by the mixing layer, in a vertical distribution range below the photic layer (50-100m depth) (Sakuma et al. 2007). Larvae of Sebastes sp. are captured mainly within the mixed layer (Sakuma et al. 2007) and increase abruptly in abundance near upwelling fronts (Moser andBoehlert 1991, Bjorkstedt et al. 2002). This finding is likely associated with surface convergence related to the increased thermocline slope, so the distribution of fish larvae may change along the pycnocline slope associated with fronts (Bjorkstedt et al. 2002). In other regions, larvae of flatfishes showed positive responses to upwelling events, suggesting that spawning is timed to maximize food availability for larvae in Algoa Bay. However, larvae hatching from benthic eggs tend to be larger and more developed upon hatching than those from pelagic eggs, and would therefore have some enhanced ability to resist transport (Pattrick and Strydom 2014).

Preference in summer-autumn
Larvae of C. xanthostigma and C. stigmaeus occur mainly in the second half of the year, the latter throughout winter , Bautista-Romero et al. 2018. During the second half of the year, a reduced core of sub-Arctic water occurs as far south as the tip of the peninsula (28°N). This leads to a progressive narrowing of the California Current near the surface as a result of mixing with warmer and more saline waters along its southerly path (Durazo 2015), governed by a warm countercurrent on the southern margin of the peninsula (Valle-Rodríguez and Trasviña-Castro 2017). Both species exhibited the best fit to the model at low salinities and a mixed layer depth. Dynamic height had little or no effect on these species, which occur mainly associated with the eastern margin of the Pacific gyre in a poorly stratified water column. Similarly, other studies have reported higher catches of Coryphaena hippurus larvae near fronts and eddies in the northern Gulf of Mexico (Kitchens and Rooker 2014).
Citharichthys spp. larvae may migrate vertically through the thermocline and pycnocline, which could favour recruitment to settlement areas near the coast Ralston 1997, Sakuma et al. 1999). In this study, flatfish larvae preferred coastal shelf waters, with some species (e.g. S. atricaudus and S. lucioceps) preferentially distributed across a wide continental shelf in Bahía Sebastián Vizcaíno. This finding is consistent with other regions where the narrow continental shelf reduces the potential spawning habitat (e.g. sardine), with clear limits in terms of temperature and bottom depth in European Atlantic waters (Bernal et al. 2007). Additional examples are the high larval density and number of coastal and shelf species on the coastal side, apparently limited by the shelf-slope front off the northern Catalan coast (Sabatés and Olivar 1996). However, in other areas the largest number of early stages occur in shelf waters of high productivity influenced by large rivers, such as the Bay of Biscay (Ibaibarriaga et al. 2007) and the Columbia River plume (Doyle et al. 1993). This is not the case of the subtropical Baja California Peninsula, a region with negligible freshwater input, favourable winds and upwelling events over the whole year, which peak from spring to early summer (Durazo 2015).
S. atricaudus and S. lucioceps are part of the socalled "Southern Coastal Complex", distributed mainly in the coastal zone around Punta Eugenia (28°N) (Moser et al. 1987). Despite their similar geographic distribution, the preferred habitat of each of these two species differs. S. lucioceps thrives mainly in relatively homogeneous waters with high dynamic height and salinity, even when their larvae occur throughout the year with peak abundance in autumn ).
For its part, S. atricaudus occurs in areas with shallow mixed layer depth and unstratified waters. This is consistent with its peak seasonal abundance (summerautumn), as it remains within cyclonic eddies around Punta Eugenia and thus avoiding offshore transport, while settlement occurs in late autumn and winter (Charter and Moser 1996).
In summary, oceanic circulation and environmental forcing have significant consequences (extension and delimitation) in the distribution of demersal fish larvae in the southern part of the California Current. Using GAMs to explore the preferred habitat allowed us to identify common environmental variables that influence larval distribution (temperature, salinity, dynamic height and degree of column stratification). This is possible since larval associations reveal that different species share similar life cycles despite the dispersive environment. The fact that these occur sympatrically (eggs and larvae) may contribute to their being affected by the same mechanisms . However, the set of variables and ranges where habitat preference is maximized is usually species-specific, and this allows a better description of the most important environmental characteristics that are closely related to the habitat used by adults and larvae.