Presence/absence data were obtained from over 1,200 sampling hauls conducted during several acoustic surveys in the Adriatic Sea. Specifically, we scrutinized data collected on board R/V G Dallaporta and R/V BIOS DVA performed in early summer and late summer between 2004 and 2021, in the framework of national (ECHOADRI and PELMON) and Mediterranean International Acoustic Survey (MEDIAS) projects (Leonori et al., 2021), as shown in Table 1 . The survey design consists of systematic parallel line transects perpendicular to the coastline or adapted to the geomorphology of the survey area, conducted during the daytime in bathymetry between 10 m and 200 m, as shown in Figure 1 (MEDIAS, 2021).

In parallel to acoustic sampling, routine biological sampling was carried out to enable identification of echo traces of fish schools and to obtain a biological sample. Therefore, pelagic trawl hauls were made on-demand at a vessel speed of ~ 4 knots for 30 min at different locations within the survey area, adapting the sampling intensity to the occurrence of echo traces and fish schools in each area. Biological sampling was performed by the two research vessels equipped with the same type of four-seam pelagic trawl net ~ 15 m horizontal and ~ 11 m vertical opening with an 18-mm stretched mesh size cod-end suitable for catching the small pelagic fish schools of juveniles and adults that are the survey’s target species. On the R/V G. Dallaporta, however, the net behavior is currently monitored in real time owing to the SIMRAD FX80 trawl sonar system, whereas on the R/V BIOS DVA, the net positions are monitored using the SIMRAD ITI system. Once on board, the total catch is weighted, and in the event of very abundant catches, a subsample is collected for further analysis.

In this study, any haul holding at least 10 sprat specimens was designated as a sprat presence location (1); otherwise, it was deemed absent (0). Several thresholds were attempted during preprocessing scrutinization. A threshold of 10 specimens was selected based on the expertise of the authors of the present work and the available dataset dealing with the minimum sample size needed for the analysis (Hernandez et al., 2006), yielding a data matrix encompassing 215 presence and 1,035 absence points (see map in Supplementary Figure S1 ).

The Adriatic Sea is an enclosed basin characterized by an extensive continental shelf that, with the exception of Pomo/Jabuka Pit, does not exceed the 200-m isobath in the north-mid-part as shown in Figure 1 (Russo and Artegiani, 1996). The shallow waters, together with the characteristic cyclonic (anticlockwise) circulation, cause high-seasonal and annual temperature variations, while the nutrient load is strongly influenced by the Po runoff in the northern part of the basin, which could lead to strong variations affecting the primary production. The freshwater inflow determines more pronounced eutrophic conditions for the western Adriatic Sea compared to the eastern Adriatic (Giani et al., 2012).

Despite the fact that the small pelagic fish biomass and distribution in the Mediterranean Sea can be affected by the Atlantic Multidecadal Oscillation (AMO) and North Atlantic Oscillation (NAO) (Tsikliras et al., 2019), the most important environmental drivers for the sprat life cycle in semi-enclosed basins, like the Adriatic and Baltic seas, are temperature, primary production, depth, salinity, and oxygen concentration (Ojaveer and Kalejs, 2010; Tsikliras et al., 2019; Schickele et al., 2021). The feeding success of sprat larvae can be affected even by wind and cloud coverage, but temperature emerged as the main driver for ontogenesis and sprat growth (Petereit et al., 2008; Voss et al., 2008; Frisk et al., 2015). For these reasons, in this work, we focused on the factors that most commonly affect the distribution of small pelagic fish that may be relevant to the study area and the target species (Fernández-Corredor et al., 2021). Specifically, satellite data derived chlorophyll-a concentration (CHL in mg m⁻³), sea surface water potential temperature (SST in °C), bottom seawater potential temperature (BST in °C), dissolved oxygen (O₂ in mmol m⁻³), net primary production (NPP in mol m⁻³), phytoplankton concentration (PHYC in mol m⁻³), sea surface salinity (SSS in PSU), and sea-level anomaly (SLA in m) were downloaded from the Copernicus Marine Services (Escudier et al., 2020; Clementi et al., 2021; Cossarini et al., 2021), in the period 2004–2021, while depth in meters was downloaded from the MARSPEC database (Sbrocco and Barber, 2013). The POLCOM-ERSEM model (Kay, 2020) of total chlorophyll-a concentration (CHLT in mg m⁻³) and sea surface water potential temperature (SSTF in °C) under IPCC RCP 4.5 and 8.5 scenarios were downloaded from the Climate Data Storage for the contemporary period (2006–2021), short-term projections (2029–2049), mid-term projections (2054–2074), and long-term projections (2079–2099). All environmental variables were downloaded at the maximum resolution, and only the layer at 5 m depth was retained for computations except for BST and SLA (see Supplementary Table S1 for more details).

When available, daily values encompassing years and months in which the surveys were conducted, (June, July, and September from 2004 to 2021) were processed in the Q-GIS environment using the SAGA plugin to obtain seasonally averaged regular grids (Conrad et al., 2015). Early summer raster files were derived from daily values for June and July, while late summer values were obtained by averaging September values. Subsequently, satellite environmental maps were interpolated through a cubic-spline interpolation and aligned to obtain a finer grid of 0.031° E (approximately one nautical mile) for each explanatory variable.

The environmental variables and the response variable were treated as late summer (data from the western Adriatic Sea from 2004 to 2014 and the eastern Adriatic Sea from 2009 to 2021) and early summer (data from the western Adriatic Sea from 2015 to 2021), merging data from different years and areas. Each dataset was randomly split into training (70% of the total hauls) and validation datasets (30% of the total samples), ensuring a balance was maintained to avoid any bias in model performance (Fielding and Bell, 1997). Since CHL and PHYC had shown a non-normal distribution, they were logarithmically transformed to achieve a normal distribution of the data. Depth was similarly transformed using a cubic root. Before applying the models, we assessed the correlation between independent environmental variables by means of the Pearson coefficient. Any variables exhibiting a significant correlation, characterized by a Pearson correlation r greater than 0.7 and less than −0.7, were independently selected (Bosch et al., 2018). Accordingly, from any possible combination of uncorrelated variables, four and three predictor sets were defined for early summer and late summer, respectively (see Table 2 ).

Using presence/absence data as response variables and the predictor sets as explanatory variables, four modeling techniques were fitted in the “biomod2” package in R (Thuiller et al., 2023), including two regression models, GAM and GLM, and two machine learning models, RF and GBM. This package has the advantage of being able also to compute an ensemble model consisting of the average of the fitted individual models. As is customary for presence/absence data, the two regression models were fitted to the training dataset assuming a binomial error distribution and a logit link function selecting the best model on the basis of Akaike’s information criteria (AIC), and deviance explained (Marra and Wood, 2011). GAM models were computed using the “GAM_mgcv” algorithm assuming a binomial error distribution and a logit link function, while the degree of smoothing, expressed as the number of knots k, was automatically selected (Wood, 2017). The two machine learning methods were set using the default settings of tuning and fitting parameters given in the biomod2. Finally, an ensemble model was computed through the weighted mean technique (EMwmean), where the probabilities derived from the chosen models are assigned weights based on their respective evaluation true skill statistic (TSS) scores.

The combination of two datasets, predictor sets, four modeling techniques, and ensemble models led to 70 fitted models. Firstly, we evaluated the variable importance through the method used in the biomod2 package. It adopts a machine-learning approach after the models have been trained. It randomly perturbs one of the variables in each permutation and calculates a correlation score between the standard prediction and the new prediction. This score is regarded as an estimate of the variable’s significance within the model. The higher the score, the greater the importance of the predictor variable in the model. It is important to note that this method does not consider interactions between variables that have been previously assessed through a correlation matrix (Thuiller et al., 2023). Only environmental variables that reached a variable importance score of 0.3 in at least one of the four models were retained for the following steps. The obtained refined predictor sets were fitted for a second run of the individual models. Next, an ensemble model was performed for each refined predictor set. The predictive performance of each model was validated by validation with the independent test dataset and quantified based on the TSS, the area under the curve (AUC), and accuracy (sensitivity and specificity) values. In particular, the TSS value was used to choose the model with the best performance, characterized by the lowest level of potential prediction error for sprat-habitat suitability projection. The model cannot accurately discriminate if the AUC is less than 0.5; an AUC of 0.7 indicates satisfactory discrimination; an AUC of 0.8 indicates good discrimination; and an AUC of 0.9 or more indicates very good discrimination (Hosmer et al., 2013). TSS is a measure of agreement between the predicted and observed values. Values between 0.0 and 0.4 indicate slight to fair model performance; values between 0.4 and 0.6 indicate moderate model performance; values between 0.6 and 0.8 indicate substantial model performance; and values between 0.8 and 1.0 indicate almost perfect model performance (Somodi et al., 2017). Accuracy is determined by sensitivity and specificity, which measure the model’s ability to properly identify the values correctly assigned as presence (positive) and absence (negative), respectively. We calculated the aforementioned metrics by comparing the predicted values to the observed values using a threshold criterion based on equal sensitivity and specificity with the exception of AUC, which is a threshold-independent measure.

Based on the validation process, the best ensemble models were used to predict the habitat suitability map of sprat in early summer and late summer for the entire period (from 2004 to 2021) in the Adriatic Sea. Specifically, we plotted the probability of occurrence of the species over a grid mesh size of 0.031°. Rare, occasional, preferential, and persistency habitat allocation maps were created by applying a threshold to the mean probability values: areas with a probability of occurrence< 0.25 were set as rare; areas with a probability of occurrence between 0.25 and 0.5 were set as occasional; areas with a probability of occurrence > 0.5 and< 0.75 were designated as preferential; areas with a probability of occurrence > 0.75 were set as persistency (Giannoulaki et al., 2017). The average seasonally mean environmental values of each selected variable were used to predict the annual habitat suitability map over the time series and seasons, providing information on the distribution of the species during seasons and areas where acoustic surveys had not been conducted.

Since there is often a noticeable shift in the scale of the IPCC-style model variables when compared to observed satellite oceanographic conditions, we compared the temperature and chlorophyll-a concentration between the hindcast (Copernicus) and model (POLCOM-ERSEM) in the overlapping period (2006–2021). Hence, we calculated the differences between the two datasets and adjusted the model-derived variable following the delta model bias correction approach described in Hare et al. (2012). The difference between hindcast T° (SST) and model T° of the overlapping period (SPTP) (ΔT_ijk ) was computed for each year (i), season (j), and RCP (k). Next, we computed the mean delta of the overlapping period ( Δ ¯ _jk ) where n is the number of years.

Accordingly, the difference between satellite CHL and modeled CHLTP (ΔCHL _ijk ) was computed through the ln difference.

Successively, the mean values for the future SST( SPTF ¯ j k ) and CHL( CHLTF ¯ j k ) were computed across three distinct time spans to encompass short-term periods (2029–2049), mid-term periods (2054–2074), and long-term periods (2079–2099). Finally, the two variables were corrected for each period (l) as follows:

After the treatment of the future environmental dataset, the best-fit ensemble models identified through the validation process for predicting the sprat-habitat suitability (hindcast) were employed for forecasting the sprat-habitat suitability under the new environmental conditions in (i) early and late summer, (ii) RCP 4.5 and 8.5, and (iii) short-, mid-, and long-term periods. Finally, we quantified the difference in the extension of the habitat suitability map of the species comparing current end future projections in terms of percentage gain and loss of area through the “BIOMODRangesize” function.

The environmental factors describing the primary production of CHL, PHYC, and NPP were significantly correlated with each other (r > 0.9; p-value< 0.0001) in both seasons. Moreover, NPP has a strong negative correlation with SSS and a strong positive correlation with BT (r > 0.75; p-value< 0.0001). Therefore, only CHL was retained as a proxy for primary production during the computations of the single and ensemble models GLM, GAM, RF, GBM, and EMwmean (the correlation matrix can be found in Supplementary Figures S2, S3 ). Four sets of predictor variables for early summer and three for late summer were used, as listed in Table 2 . Since SLA and O₂ consistently had a score< 0.3, these variables were discarded, and a second run for each single and ensemble model was performed using only the selected variables from Table 2 .

The best ensemble model performance was obtained with predictor set 1 (P1) in both seasons, as shown in Figure 2 . Therefore, the variables used in the final EMwmean for early and late summer present, past, and future projections are SST, CHL, and depth. Individual model evaluation metrics shown in Table 3 indicated that the model performs better for early summer predictions than for the late summer model. Nevertheless, in both seasons, the sensitivity scores ≥ 0.90 indicate a very low probability of obtaining false absence values, and the specificity ≥ 0.80 proves a high accuracy in detecting presence values. Moreover, the machine learning methods RF and GBM gave more accurate results than regression models GAM and GLM. Overall, SST and CHL resulted as the main explanatory variables for sprat, followed by depth in both seasons, as shown in Figure 3 .

According to response curves shown in Figure 4 , areas characterized by shallow waters and elevated primary production, as measured by the chlorophyll-a concentration (mg m⁻³) in early summer, had a greater likelihood of harboring sprat. CHL is the main driver in late summer too, when there is a shift in the preferential depth of the species to deeper waters. It is interesting to note the discrete thermal niche of sprat in early summer, when the species in the Adriatic Sea is mainly found between temperatures of 23°C and 24°C, as opposed to late summer, when the probability of occurrence increases at temperatures of approximately 21°C ( Figure 4 ).

The habitat suitability maps depicted in Figure 5 show a persistency area located in front of the Po Delta in the eastern part of the basin during early summer, followed by a reduction in late summer. The overall suitable areas are concentrated in the northeastern coastal waters during early summer, with an occasional presence up to the mid-Adriatic Sea. In late summer, these areas split into two main preferential spots, one on the western side of the basin and the other around the early summer persistency area.

Along the time series, past predictions of the sprat-habitat suitability map reveal a fluctuating extent of suitability area. The area associated with a high probability of sprat presence extended further south and offshore in both seasons in 2007, 2010, and 2014 ( Supplementary Figures S4, S5 ). Early summer is characterized by a general increase in the suitability region across the time series, with the exception of 2017 and 2018, when the suitability range was limited to the north Adriatic Sea ( Supplementary Figure S4 ). An opposite trend was observed in late summer, when the last 6 years of the studied time series (2015–2021) were characterized by a shrinkage of suitable sea grounds for sprat, which is particularly evident in 2018, along with a reduction in the western Adriatic Sea in 2011 ( Supplementary Figure S5 ).

The goodness-of-fit results of the ensemble models under the long-term RCP 8.5 scenario presented in Table 3 underline a good to high accuracy and a moderate to good performance of the final ensemble model in forecasting the habitat suitability of sprat between 2079 and 2099. The robustness of the models under other scenarios is presented in Supplementary Table S2 . Our projections indicate in early summer a loss of habitat of up to 94% for the warming scenario RCP 8.5 and up to 51% for the stabilization scenario 4.5 in the long-term period (2079–2099), as shown in Table 4 and Figure 6 . Conversely, in the short-term period (2029–2049), we observed a spread of sprat-suitability habitat toward the Croatian waters, which are characterized by lower SST ( Figures 7 , 6 ). The same shift was detected also in the mid-term period (2054–2074), when, however, the habitat loss in the northeastern Adriatic Sea overtook largely the gain ( Figure 6 ; Table 4 ). Late summer is generally characterized by a predicted migration of the species from persistency areas to deep and offshore waters. In this season, we registered a severe drop in habitat suitability wideness in the long-term period, albeit smaller compared to early summer, with a peak for RCP 8.5 of −43% ( Table 4 ).

Habitat suitability modeling is an effective technique for predicting species distribution for Ecosystem-Based Management (EBM) purposes, as it enables the prediction of species distribution by examining the nonlinear interactions between species and environmental variables (Araùjo et al., 2019). The present study used an ensemble model to present a comprehensive analysis of the past, present, and future habitat suitability maps of a sentinel species, sprat, in the Adriatic Sea in line with the good practice in SDMs (Araùjo et al., 2019), unlike previous studies conducted on sprat at a regional scale (Voss et al., 2008; De Felice et al., 2021). One of the main advantages of this approach is the characterization of model uncertainty. Along with the modeling technique, the main strength of this study comprises the amount and quality of presence/absence data employed to train and test the model. For this purpose, we used two decades of data collected during acoustic surveys in the Adriatic Sea (Leonori et al., 2021). These comprise a higher number of presence values at the local scale compared to the publicly available dataset used by Schickele et al. (2021) to forecast sprat distribution in the Mediterranean and the Atlantic seas. The use of a fishery-independent dataset may reduce the uncertainty due to the fishery behavior and make available a real absence dataset instead of pseudo-absence data which are usually applied in basin scale studies and may affect model results (Pennino et al., 2016; Schickele et al., 2021). We obtained good to high performance and accuracy of the single models and the ensemble under current and future scenarios (see Table 3 ), which validates the robustness of the data. As demonstrated by other authors, machine learning methods resulted in being the most accurate ones (Quinci et al., 2022). Conversely, we obtained a similar evaluation score between regression models and machine learning methods, which means that both approaches are able to provide good predictions of sprat-habitat suitability (Čengić et al., 2020).

Nevertheless, some bias may persist due to the consistency between present and future environmental data and variables selection. The delta model bias correction approach is considered one of the best methods for bias correction of coastal data, dealing with any possible shift in scale between climate models and current environmental variables (Hare et al., 2012). However, the choice of bias correction method could influence the spatial and temporal variability of the future projection (Pozo Buil et al., 2023). The nine environmental parameters considered in the present study may not adequately explain the whole variability in sprat-habitat suitability since other authors demonstrated the influence of other parameters such as currents, distance from the coast, wind, and cloud coverage (Petereit et al., 2008; Voss et al., 2008; Frisk et al., 2015; Tsikliras et al., 2019). However, we ended up with seven sets of parameters composed of up to four environmental variables and final ensemble models with three main explanatory variables, which explained over 40% of the model’s deviance according to GAM and GLM.

CHL was assessed as the main driver of sprat distribution, as suggested for other clupeid species (Fernández-Corredor et al., 2021). Primary production has seldom been linked to sprat presence. Nevertheless, De Felice et al. (2021) demonstrated a significant positive correlation in the Gulf of Lion as well as for SLA in the Adriatic Sea, in contrast to this work, where SLA did not show any significant influence on sprat distribution. This strong positive correlation may be mostly attributable to the diet of sprat, which consists primarily of copepods that are particularly abundant when CHL concentration levels are high (Tičina et al., 2000; Chust et al., 2014). Therefore, a clear link through food web between primary producers (CHL as a proxy of phytoplankton abundance) and secondary producers (zooplankton) as sprat’s prey and sprats is an important factor in sprat’s habitat suitability. SST, as expected, is the second key factor in determining sprat population dynamics since it is a cold-favoring species (Petereit et al., 2008; Peck et al., 2012). SST is an important factor affecting the main pelagic fish species worldwide (Lanz et al., 2009; Giannoulaki et al., 2013; Ebango et al., 2020). It affects the sprat adult growth rate and fecundity, as well as the rate of larval development and survival (Peck et al., 2012; Frisk et al., 2015). Our findings indicate an optimal thermal niche between 23°C and 24°C in early summer, which decreases to less than 22°C in late summer. Depth is the third factor that affects sprat-habitat suitability. Our results indicate that sprat are more likely to be found at depths between 35m and 65 m but not deeper than 125 m, as already demonstrated by other studies (Dänhardt et al., 2007; Voss et al., 2008). The importance of depth in determining the sprat’s habitat suitability might be related to the cooler and more dense water mass present in deeper areas, below the thermocline, and eventually to the species’ feeding behavior, which involves feeding on zooplankton and small crustaceans found in the mid-water layers (Tičina et al., 2000). In fact, the most northern part of the Adriatic Sea is an area where North Adriatic Dense Water (NAddW), the densest Mediterranean water, is generated by extreme surface layer cooling by winter cold winds, eventually driving the thermohaline circulation, ventilating the deep layers, and changing the biogeochemical properties of the Adriatic Sea (Pranić et al., 2024).

The habitat suitability maps presented here indicated that the sprat distribution area is concentrated in the northern Adriatic Sea, while the suitability decreased toward the southern regions. Moreover, our investigation reveals a west-east migration from the shallow feeding area near the Po delta in June and July to deeper waters in September for the reproductive season, where a further drop in temperature during the winter months will determine the start of the spawning season for the species in the Adriatic Sea (Dulčić, 1998). The temperature shift could be responsible for the migration from coastal to offshore and deeper and colder waters (see Supplementary Figures S7, S9 ) that we recorded in Figure 5 . The specific cyclonic circulation of the Adriatic Sea, influenced by river inflow and strong seasonal winds, could lead to a breakdown of the thermocline and mixing of the water layers in the late summer (Raicevich et al., 2015). The unique characteristics of the Adriatic Sea and the studied period could also justify the nonsignificant influence of salinity and water circulation dynamics considered in this study through SSS and SLA, unlike other studies conducted in the Mediterranean Sea and the Baltic Sea (Voss et al., 2008; Tsikliras et al., 2019).

The Mediterranean Sea is one of the most threatened basins in the world due to multiple anthropogenic pressures and its specific characteristics (Azzurro et al., 2022). It is characterized by high biodiversity and highly variable oceanographic conditions, but it is also the fastest-warming marine region in the world (Pastor et al., 2020). It is well known that abiotic changes in marine ecosystems, driven by climate changes, affect the spatial distribution of marine fauna, as demonstrated by Andrews et al. (2019) in the Barents Sea and by Ouled-Cheikh et al. (2022) in the Mediterranean Sea. Since the Adriatic Sea is a semi-enclosed basin located in the most northern part of the Mediterranean Sea, the boreal pelagic species, which are migratory and capable of covering vast distances, are unable to escape northward toward suboptimal conditions due to the semi-enclosed nature of the basin, driven by global warming. This leads to a complex orography characterized by a widely extended continental shelf with a predominant depth of 50 m in the northern region (Russo and Artegiani, 1996). These conditions expose, even more, the Adriatic Sea to the negative effects of climate change (Denamiel et al., 2020), which might be a climatic “dead end” for the Mediterranean boreal marine species if global warming continues (Ben Rais Lasram et al., 2010). Along with overfishing, this threatens the Mediterranean Sea’s most productive region (FAO, 2023).

The application of SDMs provides the possibility to assess population dynamics along the time series and under future scenarios, which can be particularly valuable as part of EBM’s purpose (O’Higgins et al., 2020). We found a noticeable contraction in the extent of the sprat-habitat suitability map under high SST and low CHL, such as in 2018, which was the warmest month of the time series, especially in late summer. Conversely, we detected an expansion of sprat-suitable areas during colder years such as 2007, 2010, and 2013 (see Supplementary Figures S6–S9 ). These findings are in accordance with the sprat’s biology and ecology, considering that the sprat is a boreal species favoring cooler habitats.

We argue that the sprat population in the Adriatic Sea is likely to be negatively impacted by the temperature increase due to global warming. An increase in seawater temperature could negatively affect the spawning success of boreal species (Peck et al., 2012; Frisk et al., 2015). The habitat suitability projections presented in the results highlight a severe decline in suitable areas for sprats under increasing carbon dioxide emission scenario 8.5. We detected a continuous decline during early summer from short- to long-term projections, which led to an almost extinction of suitable areas between 2079 and 2099, while during the late summer, we foresaw a reduction of almost 50%. The projections under the “stabilization” scenario 4.5 in late summer reveal a general shift from coastal warm-temperature waters to deeper and colder waters, underlined by a gain of new habitat in the period between 2029 and 2074. Nevertheless, albeit minor compared to RCP 8.5 projections, we observed a strong shrinkage of sprat-habitat suitability in both seasons for a long-term period, as shown in Figure 7 . Climate change is affecting fish stocks around the world, and the recently observed poleward shift in cold-favoring species is leading to a change in the pelagic communities’ composition in favor of warmer water species (Cheung et al., 2013; Kanamori et al., 2019). The northward spread of sprat has already been foreseen in the Baltic and North-Atlantic Oceans due to an increase in SST, unlike in the Adriatic and the Mediterranean Sea, where the physical boundaries will prevent the northward spreading of the species (Schickele et al., 2021). The reduction in the sprat’s suitable habitat area found in our simulations is likely to have negative effects on future recruitments of the sprat population and significant ecological and economic consequences in the Adriatic Sea, since the sprat is an important forage species in the pelagic ecosystem for many predators, including seabirds and commercially important pelagic fish such as swordfish, tuna, and mackerel (Fasola et al., 1989; Knudsen et al., 2009; Navarro et al., 2017). The gradual decline in sprat biomass is therefore likely to have cascading effects throughout the ecosystem, both in terms of bottom-up control of predators and top-down control of zooplankton abundance. Consequently, other small pelagic species, such as anchovy in the Mediterranean Sea, may be more suited to forthcoming local conditions, competing with S. sprattus for the same ecological niche, being less affected by climate change (Bourg et al., 2015; Schickele et al., 2021). The trophodynamic role of sprat in the Adriatic Sea may be gradually replaced by the warm-temperate Clupeid species round sardinella, Sardinella aurita, which is already spreading its area of distribution northward, from the south-central Mediterranean Sea to the north Adriatic Sea, where its presence is currently considered uncommon (Schickele et al., 2021).