Unveiling the Relationship Between Sea Surface Hydrographic Patterns and Tuna Larval Distribution in the Central Mediterranean Sea

Thunnus thynnus (Atlantic bluefin tuna, ABT) and other tuna species reproduce in the Mediterranean Sea during the summer period. Despite the Central Mediterranean Sea, the Strait of Sicily in particular, being a key spawning site for many tuna species, little is known on the effects of oceanographic variability on their larval distribution in this area. The abundance and presence-absence of larval specimens for three tuna species (ABT, bullet tuna and albacore) were modeled in order to examine their relationships with environmental factors, by analysing historical in situ information collected during seven annual surveys (2010–2016). The results revealed that most tuna larvae for the three species were found in the easternmost part of the study area, south of Capo Passero. This area is characterized by a stable saline front and warmer nutrient-poor water, and it has different environmental conditions, compared with the surrounding areas. The models used to investigate the presence-absence and abundance of the three species showed that ABT was the most abundant, followed by bullet tuna and albacore. The presence and abundance data collected are comparable with those of other spawning areas in the Mediterranean. Regarding biological and physical parameters, the results suggest that temperature, salinity, and day of the year are the key factors for understanding the ecological mechanisms and geographical distribution of these species in this area. Temperature affects the presence of ABT larvae and salinity, which, with a physical barrier effect, is a key factor for the presence-absence of bullet and albacore and for albacore abundance.


INTRODUCTION
Atlantic bluefin tuna (Thunnus thynnus, Linnaeus, 1758) is one of the most valuable tuna species in the world. Exploited by humans since ancient times, it has long been a valuable resource whose meat can now reach one of the highest prices in the international market (Mylonas et al., 2010;Sun et al., 2019). The International Commission for the Conservation of Atlantic Tunas (ICCAT) manages this species and considers two different stocks: the western stock whose main spawning ground is in the Gulf of Mexico (although according to new evidence spawning also occurs in the Slope Sea), and the eastern stock that spawns in the Mediterranean Sea (Muhling et al., 2017). Within the Mediterranean Sea, three main areas are well recognized, the western area located around the Balearic Islands; the central area in the Tyrrhenian Sea and around Sicily, northeast of the Gulf of Sirte, Malta, and Tunisia; and the Eastern area that comprises the Ionian Sea and the area around Cyprus (Karakulak and Yıldız, 2016;Muhling et al., 2017).
There is substantial discussion about the sub-population structure of Atlantic bluefin tuna (ABT) within the Mediterranean Sea and the occurrence of spatio-temporal mixing during spawning. The major spawning fraction is formed by large migrant adults that enter from the Atlantic Ocean following the Atlantic waters through the Strait of Gibraltar from April to May. Spawning occurs from June to July (Alemany et al., 2010;Reglero et al., 2012;Abascal et al., 2016). Some studies suggest that there is also a resident fraction of adults that remain in the Mediterranean Sea throughout the entire year or at least for a greater part of it (Cermeno et al., 2015;Livi et al., 2019). These individuals feed in different locations during the winter (Cermeno et al., 2015) depending on their ontogenetic stage (Sarà and Sarà, 2007). Although there is no conclusive scientific evidence, this resident fraction is believed to spawn in the Eastern Mediterranean and may also be in the Central Mediterranean (Cermeno et al., 2015). Migrants and residents may mix during reproductive events (Livi et al., 2019), and adult tagging data show that some individuals that remained in the Mediterranean Sea move to the Central Mediterranean during the spawning period where large migrant adults can also be found (Block et al., 2005;Fromentin, 2010). Studies to investigate the lesser-known spawning areas, such as the central Mediterranean Sea, are crucial to increase the understanding of the dynamics of the Mediterranean stock of this species. As such, questions that remain unanswered warrant attention.
The presence of larvae in the Western Mediterranean, one of the best studied spawning grounds, suggests ABT to be an environmental spawner tightly linked to the less saline Atlantic waters that enter through the Strait of Gibraltar and flow northward. There, they encounter the more saline resident water mass around the Balearic Islands, a high dynamic mesoscale area with a marked front (Alemany et al., 2010;Reglero et al., 2012;Alvarez-Berastegui et al., 2014). The Atlantic waters then enter the Strait of Sicily, flowing toward Tunisia and the South of Sicily, from west to east (García Lafuente et al., 2002;Cuttitta et al., 2004;Bonanno et al., 2014). A consecutive sampling of tuna larvae over a 3-year period has shown that the presence of tuna larvae in Tunisia is also related to this water mass (Koched et al., 2012(Koched et al., , 2013(Koched et al., , 2016Zarrad et al., 2013). In contrast, little is known about the south of Sicily, an area characterized by high mesoscale variability, vortices, and meanders (García Lafuente et al., 2002) because of incoming Atlantic freshwater and warmer resident water (Robinson et al., 2001). These environmental data have not yet been analyzed with tuna larval data from oceanographic surveys in the study area.
Another common characteristic in the Balearic Islands and Tunisia is that ABT larvae are commonly found together with two other tuna species, bullet tuna (Auxis rochei, Risso, 1810) and albacore (Thunnus alalunga, Bonnaterre, 1788), suggesting that the three species have overlapping spawning grounds (Alemany et al., 2010;Koched et al., 2013). Evidence from gonad sampling suggests that the timing of spawning occurs around the same time in the Western and in the Central Mediterranean, between June and July (Heinisch et al., 2008), a pattern that is corroborated by temporal changes in the abundance of larvae (Alemany et al., 2010).
Considering these basic knowledge gaps and the apparent need for in-depth analysis, this study aims to better analyze the lesser-known central Mediterranean Sea spawning area. The objectives are to better understand the distribution and abundance of the three species and relate them to the biotic and abiotic factors of the area, thus verifying their main ecological drivers. This can enhance the knowledge of the Central Mediterranean, expanding the ecological knowledge of these species within the basin. We explored the environmental drivers affecting ABT, bullet and albacore tuna at the larval stage by analyzing presence-absence and abundance data, focusing on the early life stage.

Field Sampling
The samples are from ichthyoplankton surveys carried out in the Strait of Sicily by the Italian National Research Council (CNR) during the summertime (June-August) from 2010 to 2016. A regular sampling grid was used (1/10 • × 1/10 • along the continental shelf, and 1/5 • × 1/5 • offshore). The number of sampled sites investigated in each year, as well as the sampling period, is summarized in Table 1. The sites are located along transects perpendicular to the southern Sicilian coast (Figure 1).
Ichthyoplankton samples were collected 24 h a day, with oblique tows of a Bongo 40 net with 200 µ mesh size. It was equipped with a flow meter and towed at a ship speed of 2 knots, with a descent speed of 0.75 m/s and an ascent speed of 0.33 m/s. The net was lowered from the seabed to the surface, but sampling depths did not exceed 100 m.
The collected samples were immediately stored onboard in 70% ethanol and later processed in the land-based laboratory for identification of tuna larvae at the highest possible taxonomic level according to Rodríguez et al. (2017). Following larvae sorting, the remaining zooplankton was collected to obtain dried zooplankton weight. Water column temperature Albacore 0 0 0 ( • C), and salinity (PSU) measurements were gathered at all sites using the multi-parameter CTD SBE 11 plus probe (underwater unit).

Data Exploration and Statistical Modeling
Mean larval densities at stations with different bottom depths were explored to identify potential spatial preferences regarding bathymetry. The distribution of observed larval presence and absence was explored using density plots. Generalized additive models (GAMs) allowed us to investigate the relationships between the spatial distribution of tuna larval specimens and environmental factors. Presence-absence data, from all the stations, and abundance data, only from positive stations, were used as dependent variables in the models. This helped determine the main covariates that control larval occurrence (calculated as the probability of presence) and the density of tuna larvae per unit area. The binomial distribution family was selected for presence-absence data and the Gamma family for abundance. To transform the catch effort for abundance, we standardized the number of larvae with the following formula: where CPUA is the "capture per unit area", larvae/m 3 is the number of larvae divided by the filtered water volume in m 3 , and BD is the depth in meters reached by the bongo net during sampling (Alvarez-Berastegui et al., 2018). Considering the volume of filtered water and the maximum depth of the sampling gear, we were able to obtain an unbiased and normalized expression of larval density per m 2 . The first model building phase was dedicated to selecting variables that best describe the system and that were potentially useful for understanding the tuna larvae spatio-temporal distribution. This step also considered information from similar studies performed on other fish species. We started by including all the available and suitable data for the analysis, such as geographic coordinates and dissolved oxygen, but excluding stratification of the surface water, which is a typical condition characterizing the summer season in the southeastern sector of the study area, where most of the tuna larvae were found; and terrestrial runoff, because of limited drainage basin in an island such as Sicily. Finally, tidal excursion was excluded as this is very limited in the waters of the study area.
By Pearson' s correlation, the collinearity among the variables was analyzed. In case of a significant correlation between two variables, we selected the one that fitted best with the dependent variables and excluded the other. For example, longitude was positively and significantly correlated to temperature anomaly, mostly because of the upwelling in the northwestern zone (Supplementary Material), which, in some years, was stronger and produced lowest temperature anomaly values.
The last covariates considered in the models were: "Year" as a factor and with random effect, in order to focus the analysis on environmental relationships and not on interannual fluctuations; Frontiers in Marine Science | www.frontiersin.org surface "Salinity" (PSU) as recorded by the CTD probe at 5 m; "Day of the year, " to verify the presence-absence and abundance peak, due to the temporal choice for spawning as dictated by the adults; "Zooplankton density" (standardized for the volume of filtered water, dry weight mg/Vol), to understand the role of possible food availability; "Residual of superficial temperature" ( • C), calculated as the residual of a GAM, where temperature (at 5m) was fitted to the factor "Day of the year", in order to avoid the correlation between these parameters and allow for the different timing of annual surveys.
The final formulas used for the (i) presence-absence model and (ii) larval abundance model were as follows: A step-wise procedure allowed the progressive removal of nonsignificant covariates (p > 0.05), ensuring that only independent variables were used and non-significant variables were excluded. This had the effect of the most significant variables being easier to understand, removing interferences in the additive model, dictated by less significant variables. All simulations, statistical analyses, and plotting were performed using the R statistical software (RStudio Inc, 2016). In particular, GAMs were fitted using the mgcv library (Wood, 2011).

Spatio-Temporal Distribution
While the three species under investigation [Atlantic bluefin tuna (ABT), bullet tuna and albacore] are commonly found in the samples (Table 1), other tuna and tuna-like species are sporadically found, such as Sarda sarda, Euthynnus alletteratus, and Katsuwonus pelamis. ABT larvae were constantly more abundant in the eastern side of the study area (between 14.5 and 15.5 • E longitude) and occurred every year between Capo Passero (the southernmost tip of Sicily) and Malta Island (Figure 1). Similar annual abundance was found, but in 2016 there was higher average abundance ( Figure 1A). Bullet tuna larvae had a geographical distribution similar to that of ABT annually, although they were more numerous closer to the coast and more prevalent east of longitude 15 • E ( Figure 1B). Regarding larval densities, low larval densities were found in the western part of the Strait, while higher densities were found along the eastern end of the Sicilian coast and Malta ( Figure 1B).
Lastly, albacore, the tuna species found to have the lowest presence probabilities and abundance at positive stations, was found between 2011 and 2015 ( Figure 1C) and in well-defined areas at the times of samplings (the shelf area between Capo Passero and Malta), with rare presence between longitudes 12 and 14 • E ( Figure 1C). Salinity data showed saltier water masses in the eastern zone for all the years (Figure 2). A thermohaline front was found in the same area, every year, beyond 15 • E of longitude.

Hydrographic Conditions
The Zooplankton data showed oligotrophic waters every year, except for a few points with higher biomass in 2010 and general greater abundance in the western sector (Figure 2).

Larval Data and Environmental Factors
The mean larval density related to the bottom depth showed that the three species were distributed in areas with different bathymetric profiles (Figure 3). ABT was the species distributed more offshore, in areas where the bottom depth is deeper than 100 m. Bullet tuna distribution was more coastal, showing higher probabilities of presence in shallower stations (0-100 m), although individuals were found in the entire bathymetric range. Albacore was present evenly in bathymetric range deeper than 50 m. Regarding temperature, ABT larvae were present at temperatures higher than 21 • C. Similarly, bullet tuna was found in waters warmer than 20.9 • C, and albacore was found in waters warmer than 20 • C (Figures 4A-C). With reference to salinity, the presence of the three species showed two peaks at around 37.5 and 38.5 (Figures 4D-F). The GAM analysis explored the relationship between larval distribution and putative covariate factors that may drive the presence and abundance of the three species ( Table 2).
The presence of Atlantic bluefin tuna (ABT) was positively related to "Day of the year, " "Temperature residuals, " and "Year" (p < 0.001 for all the three variables) (Table 2, Figures 5A1,A2). The presence of ABT increased during the spawning season, reaching a maximum around day 200 (July 19th). The larvae were found in warm waters, and the probability of their presence was higher when temperatures exceeded the average temperature trends. The presence of bullet tuna was positively related to "Day of the year" (p < 0.001) with higher probability of finding larvae as the season advances. Regarding "Salinity" (p < 0.05), bullet tuna showed a flat trend of up to 38 PSU and increased from that point onward (Figures 5B1,B2). As seen for bullet tuna, the presence of albacore increased with the "Day of the year" (p = 0.01) and with "Salinity" (p < 0.001) (Figures 5C1,C2).
Regarding the model used to analyze the abundance in positive stations, ABT was found to be negatively related to "Day of the year" (p < 0.01) (Figure 6A). For the bullet model, even with the stepwise procedure, none of the explanatory variables were significant. The albacore model for abundance in the positive stations showed an increasing linear trend with "Depth, " a dome-shaped relation with "Salinity" with a maximum at intermediate levels; and with "Zooplankton abundance, " the trend decreased until intermediate values remained constant at higher zooplankton abundances (p < 0.05) (Figures 6B1-B3).
For this species, it should be considered that in the eastern sector where the larvae are more abundant, in the 2 years with greatest zooplankton abundance (2010 and 2016), we did not find   albacore larvae. It is also necessary to consider that there are very few events of large zooplankton abundance and that from the low to medium abundance that has been detected, the effect on albacore abundance decreases (Figure 6B2).

DISCUSSION
The results indicate that the central Mediterranean Sea, in particular the Strait of Sicily, is a significant larval habitat for the three tuna species considered in this study. ABT and bullet were generally more abundant than albacore, the latter having not been caught in two of the years during the sampling time series. Therefore, there is a temporal-spatial overlap during the early life stages of these species. Temperature, salinity, and day of the year emerged as key variables to better understand the spatio-temporal distribution of larval habitats of these species in the area. Tuna larvae were found in the eastern side of the Strait, where surface water temperatures were found to be consistently warmer than those of the surrounding areas. The GAM results on the presence-absence of ABT larvae show that the day of the year is highly significant, with a dome-shaped relationship reaching a maximum at around day 200, corresponding to the beginning of the second half of July. Tuna species are found in waters above 22 • C, and the timings are consistent with studies carried out in the Western Mediterranean (e.g., Alemany et al., 2010;Reglero et al., 2012) and south of the Central Mediterranean (e.g., Koched et al., 2013;Zarrad et al., 2013). The results reinforce the idea that tuna eggs and larvae are mainly found in oceanic waters and are constrained by temperatures above 20 • C, at which they can develop and survive . In the study area, we also found some larvae in the August period, suggesting that some adults may have spawned outside the common spawning period, when temperatures are still suitable for egg-hatching.
There was a higher probability of finding larvae at stations where the temperature was higher than the average temperature for that day. Adults can choose where to spawn and may choose to spawn when and where temperatures are favorable for the development of their larvae. ABT have very high fecundity and are multiple batch spawners (Ciannelli et al., 2015). While they are in the spawning grounds, they have a very fast ovarian development and oocyte maturation that is propelled by the presence of increasing warming water, which poses an evolutionary constraint on the spawning locations (Ciannelli et al., 2015). The results suggest that the moment in which adults choose to spawn increases up to mid-July and then decreases. Interestingly, in the study area, we also found some larvae during August. This indicates that some adults have spawned outside the common peak period. The reproductive success of these events that occur later in the year is not known. These smaller ABT that spawned later in the year might be found by larger conspecific larvae that spawned during previous spawning events. If these larger larvae have reached the piscivorous phase, they might feed on their conspecifics and carry out cannibalism (Reglero et al., 2011). As Bongo-40 nets that predominantly catch small larvae were used, statistical analyses on the difference in size distribution across the seasons and/or stations were not possible.
Regarding salinity, the minimum, maximum, and mean values recorded in the stations are typical of the study area, characterized by the confluence of Atlantic waters with Eastern Mediterranean waters in the Ionian Sea. In fact, we have detected a strong variation in salinity values, corresponding to a key oceanographic structure, the saline front in the south of Capo Passero, to the east of the study area. This structure appears to be stable during the summer period, confirming previous observations (Patti et al., 2010;Bonanno et al., 2014;Cuttitta et al., 2016). Across the front the larval occurrence is higher, and this is evident for the three species that show two presence peaks related to salinity. Larvae originating from the spawning activity in the southeastern sector of the study area accumulate with larvae advected by the Atlantic waters from the offshore western area, causing the higher concentrations observed in the retention area for the three species investigated. In the same area, this advection mechanism has also been verified for other species in previous studies (García Lafuente et al., 2002;Cuttitta et al., 2018;Patti et al., 2018). So, definitively the dynamics of this frontal area seems to play a key role for larval retention in the eastern zone. This hypothesis has been confirmed by the GAM models, as bullet and albacore tuna showed a higher probability of presence at higher salinity. This is not dictated by physiological constraints but by a purely physical factor. The different densities of the two water masses where we find higher larvae presence and the presence of the saline front every year suggest that this oceanographic structure works as a physical barrier, concentrating the tuna larvae. Geostrophy and associated surface currents play a key role in determining the fate of the early life stage, and Capo Passero assumes an important function as a retention area for eggs and larvae, not only for tuna but also for other species that gather from the northwest, resulting from along-shore advection Patti et al., 2020).
Atlantic bluefin tuna (ABT) larvae were absent in the western sector. This area is characterized by coastal upwelling zones due to mistral winds that induce offshore Ekman transport of surface waters, causing the upwelling of cold and nutrient-rich water masses (Patti et al., 2010;Torri et al., 2018). In general, tuna larvae are found in oligotrophic areas, probably as a strategy to encounter fewer predators (Bakun, 2013). Despite the eastern area being less productive than the western area, local enrichment processes such as cold filaments from the upwelling zone together with the thermohaline front and eddies could concentrate prey, improving feeding opportunities (Cuttitta et al., 2016Torri et al., 2018). These oceanographic processes could create patches of higher food abundance that, together with scarcity of predators, could increase the chances of larval survival.
It was not possible to depict the reproductive peak of bullet and albacore tuna in this study, since the data suggest that the spawning window of bullet is longer than the timing of the surveys (Allaya et al., 2013). This is demonstrated by the "Day of the year" results. The bullet tuna reproductive strategy is based on a protracted spawning behavior, compared with ABT. Concerning bullet tuna abundance, the lack of significant variables in the applied model highlights how the larvae of this species are heterogeneously distributed within the study area ( Figure 1B). We chose to implement a unified model approach for the three species, both for presence-absence and abundance, but clearly, each species has its own peculiarities, and the results suggest the need for further species-specific studies. Moreover, because of the complex nature of the ecological data, low values of explained deviance are common, also considering that other environmental variables, not available for this study, could potentially be useful to improve the knowledge of larval distribution.

CONCLUSION
The study area located in the center of the Mediterranean Sea, with direct influence of Atlantic and Levantine waters, creates unique environmental conditions selected as optimal spawning grounds by different tuna species as has been observed for other species (Falcini et al., 2020;Mangano et al., 2020;Patti et al., 2020). In this study, similarities with other well-characterized spawning grounds in the Central and Western Mediterranean were highlighted, where the three selected tuna species also coexist (Alemany et al., 2010;Koched et al., 2013). This study targets the lesser-known central Mediterranean Sea spawning area, revealing that temperature, salinity, and day of the year were the main variables driving tuna larval geographical distribution. However, the lack of a long time series in the region and the different sampling methods in the area compared with the other areas made the detection of significant temporal trends in the observed patterns difficult. As such, this study has given some interesting insights and paves the way for future studies on the impact of environmental variability on the interannual fluctuations in the abundance of these species and in their reproductive success in the Central Mediterranean Sea.
More research in the Central Mediterranean tuna spawning grounds is needed to understand its role compared with other areas in the Mediterranean Sea. It is also possible that there is a spatio-temporal mixing in spawning grounds between the two ABT Mediterranean sub-populations. The crucial issue of the ABT Mediterranean stock structure remains inconclusive. The results could indicate a reproductive event or multiple reproductive events, resulting from Atlantic or resident ABT, or both. We still need to know more about this area, as we still do not know the larval origin and if the larvae come from the resident adult fraction or not. By studying these lesser-known spawning areas, we wanted to take the first step, laying the foundations, to fill critical gaps in the understanding of the unique dynamics of the Mediterranean stock for these species that, to date, have remained unanswered.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
SR, AC, and GS contributed to the conception and design of the study. MT, BP, and AC carried out the sampling. SR, AC, and PR taxonomic identification. SR, MT, and BP organized the database. SR, MT, DAB, and PR performed the statistical analysis. SR wrote the first draft of the manuscript. GS and AC secured funding. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
Fondo Sociale Europeo Sicilia 2020 supported SR's PhD research presented in this study. Data collection was mainly supported by the Italian National Research Council (CNR) through the USPO office, and by the FAO Regional Project MedSudMed Assessment and Monitoring of the Fishery Resources and the Ecosystems in the Straits of Sicily, funded by the Italian Ministry MIPAAF and co-funded by the Directorate-General for Maritime Affairs and Fisheries of the European Commission (DG MARE). Sampling activity was also supported by other research projects, namely, Grandi Pelagici, funded by the Sicilian Regional Government, and SSD-PESCA and the Flagship Project

ACKNOWLEDGMENTS
We thank all the CNR and IEO-COB technicians for the support during species identification, especially Melissa Martin Quetglas and the PhD students (Catalina Mena Oliver and Daniel Ottmann Riera) for the help during preliminary data exploration, Dr. Luigi Giaramita and Carlo Patti for their support during oceanographic surveys and sampling collection. We thank Francisco J. Alemany, Walter Ingram, Lorenzo Ciannelli, and Hanane Elyaagoubi for their helpful comments on the preliminary study. We thank Mr. Emanuele Gentile, master of the R/V Urania and R/V Minerva Uno, and his crew for their work in support of plankton sampling during oceanographic cruises.