Sampling on a seasonal basis (March, June, September, and December) was conducted in three consecutive years (2017, 2018, and 2019) at 20 sites (Table 1 and Figure 1) across marine and transitional waters along the eastern Adriatic coast. The specific ecological characteristics of each site are presented in Table 1. All sites were selected based on knowledge where juvenile fish occurrence had been previously described (Dulčić et al., 1997, 2005, 2007; Matić-Skoko et al., 2007) over a similar biotope characterised mainly by sandy (S1, S3, S5, S6, S9, S10, S12, S15, S18, and S20), rocky (S11 and S16), or mixed sandy-rocky-pebbles (S2, S4, S7, S8, S13, S14, S17, and S19) substrata covered by photophilic algae alternating with patches of Cymodocea nodosa seagrass beds. At each site, three replicates separated by tens of meters were performed with a specially constructed small shore seine (L = 25 m; minimum mesh size 4 mm). To ensure comparability, the same net was used during the entire study period at all sites. Hauls were also performed in the same bathymetric range, from 0 to 2.2 m depth. A total of 36 samples per site were taken over a 3-year period. Sites were divided into four groups of five locations based on water type (marine or transitional) and geographical position (north or south): transitional north (T_N), transitional south (T_S), marine north (M_N), and marine south (M_S). The border between the northern and southern sites is Cape Ploča (43°29.648′N; 15°58.184′E) as the most prominent point of land along the eastern Adriatic coast. This cape represents the geographic and climatological boundary between the northern and the southern Adriatic, and is often characterised by strong sea currents and swells as weather systems from the north and the south come in contact (Gloginja and Mitrović, 2021). Surface and bottom temperature (°C) and salinity (‰) measurements were performed using a multi-parameter probe model YSI 85. All fish individuals were identified to the species level according to Jardas (1996) and Dulčić and Kovačić (2020), measured (total length to the nearest 0.1 cm) and weighed (total body weight to the nearest 0.1 g). Sampling covered all stages of fish, but only juveniles were selected for further analysis. Further on, the species Atherina sp. were excluded, since they were the most dominant and most common species in all samples, and the majority of specimens were adults, while Sardina pilchardus and Engraulis encrasicolus were caught only sporadically in high abundance, thus skewing the results. The following species were excluded from further analysis as only adults, and no juvenile specimens, were caught: Lichia amia, Pseudocaranx dentex, Sciena umbra, Uranoscopus scaber, and Zeus faber.

Total abundance (A), as the number of individuals (N) per species, and species richness (R), as the number of species (S), were counted, while the mean number of individuals and species per haul were recorded for each site and each period.

All analyses were performed using PRIMER-E software (Clarke and Gorley, 2015) with the add-on package PERMANOVA+ (Anderson, 2001; Anderson et al., 2008). Graphs were prepared using SigmaPlot (v. 13.0; Systat Software Inc., San Jose, CA, United States).

Multivariate statistical testing of abundance data for juvenile fish assemblages used a four-factor design [Year (Ye)/Season (Se)/Geographic position (Ge)/Type of water (Ty)], within one factor Period (1–12): 2017 (1–4), 2018 (5–8), and 2019 (9–12) constructed for cyclicity purposes. Factors Ye and Se were considered fixed while factors Ge and Ty were randomly nested in Se. The additional factor Period_Type of water_Geographic position was combined. Data sets for each site were imported in the PRIMER workspace and combined into a single matrix. Variability in the numbers of individual species in the samples was used to carry out dispersion weighting for each species. The dispersion-weighted data was transformed by a milder square-root transformation, as demonstrated by Clarke et al. (2014) for fish communities, particularly those in estuaries and nearshore coastal waters, where the prevalence of juvenile and small schooling species is high. Bray–Curtis similarity was then calculated on the dispersion weighted, transformed data.

Metric Multidimensional Scaling (mMDS) ordination used data at the same level of sample averaging for each Period, Geographical position, and Type of water combination to visualise the extents to which ichthyofaunal composition differed across Period_Type of water_Geographic position. All plots were constructed by calculating the distances between each pair of group centroids, i.e., the relevant average in the “Bray–Curtis space” of all samples (Anderson et al., 2008).

The RELATE statistic (P) was employed to test whether the pattern of temporal catch composition change conformed to cyclicity. If there is no tendency to cyclicity, then P will be close to zero (Clarke and Warwick, 2001). Separate two-way crossed analysis of similarities ANOSIM was used to interpret the relative size of the overall spatio-temporal effects on fish compositions, using the same resemblances as for the PERMANOVA tests. The species that mainly contribute to the separation between the factor Period_Type of water_Geographic position were determined by the SIMPER procedure, which shows the average contribution of each species to the dissimilarity between the groups. The results were presented for the specific season along all three consecutive years, characterised by most significant changes in communities.

To explore the data further, we used canonical variate plot (CAP) analysis. CAP is a routine for performing canonical analysis by calculating principal coordinates among groups of samples to predict group membership, positions of samples with another single continuous variable, or finding axes having maximum correlations with another set of variables (Anderson and Willis, 2003). This was initially run separately for each of the two factors: “Geographic position” and “Type of water” and then merged into a scatter plot. Particularly, CAP was used to estimate the accuracy of fish composition and environmental variables (temperature and salinity).

The mean monthly variation in water temperature and salinity averaged over four groups (M_S, M_N, T_S, and T_N) are shown in Figure 2. The observed sea surface temperatures followed a seasonal cycle, with the highest temperatures in June-September (particularly in 2018) and the lowest in December in all three years (Figure 2A). Observing the differences between marine and transitional waters, less variation was observed among waters in the south, regardless of the type of water (Supplementary Table 1). The mean seasonal temperature in T_N was always lower than in the remaining groups, though this is especially evident in December when higher and more similar results were recorded in the other three groups compared to T_N. Namely, northern transitional waters had the lowest mean values recorded in December 2017 (10.12 ± 1.81°C) and December 2018 (10.44 ± 1.39°C). It is evident that the mean temperature in December 2019 was higher for the T_N group than in the two previous years, and this was much more uniform among the groups (Figure 2A). Further, the mean temperature in June 2019, regardless of water type, was 2.78°C lower than in 2017 and 4.78°C lower compared to 2018. Analysis of variance revealed significant differences in temperature between years (PERMANOVA: Pseudo-F = 10.371; P = 0.0004), seasons (Pseudo-F = 50.181; P = 0.0013), interactions of years and seasons (Pseudo-F = 10.624; P = 0.0002), and among groups (M_S, M_N, T_S, and T_N) in relation to seasonal changes in temperature (Pseudo-F = 4.6813; P = 0.0008) (Table 2).

As expected, the measured salinity values showed greater variation between water types such that marine areas showed less seasonal fluctuations and continuous higher values (Figure 2B and Supplementary Table 1). Lower salinity values were observed for northern transitional waters (mean salinity value across the sampling period was 23.53 ± 6.23‰) than southern transitional waters. Two negative peaks of mean salinity values were observed in June and December 2019, regardless of water type. The mean salinity was always higher in 2017 and 2018 in the south, regardless of water type, while an inversion occurred in 2019 and mean salinity values were higher in the north in both types of waters. Also, every September in all three years, the mean salinity values were highest and most similar among the groups following the summer drought (Figure 2B). However, analysis of variance did not reveal significant differences in salinity between years, seasons, type of waters, and geographic position and their interactions, except for strong statistically significant annual variations among water type/geographic position groups (M_S, M_N, T_S, and T_N) in relation to seasonal changes in salinity, i.e., as nested in season (Pseudo-F = 1.9354; P = 0.0001) (Table 3).

At all sites, a total of 720 small beach seine samples were taken across three consecutive years (Tables 1, 4). A total of 24,807 fish juveniles belonging to 23 families and 92 species were included in this survey. The most abundant species were Sarpa salpa (2,805 individuals in total), Pomatoschistus marmoratus (2,028), Pagellus acarne (1,806), Chelon ramada (1712), Diplodus annularis (1,683), Chelon auratus (1,550), and Diplodus vulgaris (1,074), which together comprised 51.03% of the total sample. A high number of species were sporadically present throughout the 3-year study period (though with <30 individuals). Some species, such as Symphodus roissali and Chelidonichthys lucerna, were completely absent in 2019. Although the collection effort (as the total number of samples per year, month and site) was the same, the total number of species and individuals clearly fluctuated across sites (Table 4, Figure 3, and Supplementary Table 2) and across groups. The number of species varied among sites from 1 (at S18 and S19 in June and September 2019, respectively) to 24 species (at S16 in June 2018) and among the groups from 1 (T_N in December 2019) to 25 species (M_N in June 2018). The number of species varied less over the 3-year period at marine stations than at transitional ones, with the largest interannual changes observed in the T_N group. Moreover, PERMANOVA analysis of the number of juvenile fish species found significant differences along three consecutive years (Pseudo-F = 1.7459; P = 0.0362), across seasons among years (Pseudo-F = 1.4222; P = 0.0406) and among groups [M_S, M_N, T_S, and T_N, as Ge(Se) × Wa (Se) interactions] in relation to seasonal changes in the number of species (Pseudo-F = 2.0571; P = 0.0002) (Table 5).

Similar to the number of species, abundance fluctuated considerably less in marine waters and more in transitional ones, especially in the north. The abundance of individuals was the highest in 2018, and the highest in M_N group (4067 individuals). On the contrary, abundance was lowest in 2019, especially in the M_S (1,317) and T_N groups (1,407). Although some fluctuations in abundance were evident, especially the 52.5% reduction in the M_N group in 2019 (1,930), PERMANOVA analysis of juvenile fish abundance for all sites found no significant differences for the three consecutive years, seasons, type of water and geographical position together with their interactions, except among groups in relation to seasonal changes in the number of individuals (Pseudo-F = 1.9354; P = 0.0001) (Table 6). Regarding species, the abundance of Boops boops, Pagellus acarne, Pagellus erythrinus, Symphodus tinca, and S. roissali was drastically decreased in 2019 compared to previous years (Table 4). Additionally, P. acarne was only present in southern marine waters in 2017 and 2018.

On the metric MDS ordination plot, derived from the distance among centroid matrices, four groups were defined, clearly separating both geographical position and type of waters (Table 3 and Figure 3) ordinating groups from transitional to marine and then from north to south (T_N, T_S, M_N, and M_S), indicating that the greatest differences were found between the groups T_N and M_S (Figure 4). The cyclicity of abundance data, obtained from the RELATE test, following a seasonal pattern and was less visible for T_N. It seems that there were similar patterns for marine groups, namely between M_S in 2017 (Rho = 0.621, P = 0.839) and 2018 (Rho = 0.414, P = 0.665) and M_N in 2017 (Rho = 0.828, P = 0.335) and 2018 (Rho = 0.828, P = 0.333) while in the final year of the study period, similarity dropped drastically. Interestingly, there was no similar pattern between the three consecutive years at T_S (Rho = 0.014, P = 0.365). Species composition in T_N was governed by a similar pattern in 2018 (Rho = 0.414, P = 0.679) and 2019 (Rho = 0.414, P = 0.035). Two-way crossed ANOSIM tests showed and confirmed that the fish fauna composition was significantly related to both factors, types of water (R = 0.106, P = 0.001) and geographic position (R = 0.188, P = 0.001), with the latter factor more prominent. Moreover, one-way ANOSIM for the interaction type of water_geographic position confirmed that the magnitude of difference was lower between groups in 2017 (R = 0.419; P = 0.001) and 2018 (R = 0.458; P = 0.001), and highest in 2019 (R = 0.899; P = 0.001).

We additionally ran a CAP analysis and related the environmental data to the community composition information. Temperature and salinity values for the four water types/geographic position groups of factors were plotted as distances among centroids based on community data (Figure 5), revealing that temperature had a prominent influence on community composition in T_N, less on T_S, while salinity had a positive effect on M_S and M_N. A separate CAP analysis was run for each of the two factors (“type of waters” and “geographical position”), which also gave successful discrimination. In particular, 84.03 and 78.99% of samples from marine and transitional waters, respectively, were correctly allocated based on the community composition information. Further, 72.88 and 71.67% of samples were correctly allocated to north or south, respectively. These results suggested that transitional waters are more sensitive to interannual differences, especially in the north (Figure 6).

The SIMPER routine revealed that 14 species (Table 7) in varying order of percentage contribution were responsible for the most (>50%) dissimilarities in the interactions between the composition of juvenile fish sampled in June among the three consecutive years. It was evident that five species that gave a significant contribution to T_N in 2017 (e.g., Sparus aurata, Diplodus vulgaris, and Symphodus roissali) were not present in 2018, while other species appeared (Sarpa salpa and Chelon ramada). SIMPER revealed that the average dissimilarity was highest between M_S and T_N (90.51%) while T_S and M_N were more similar (79.41%).

Generally, the determined fluctuations in sea temperature and salinity during the study period followed the well-known seasonal cycle specific to the Adriatic Sea (Grbec and Morović, 1997). More pronounced interannual fluctuations in water temperature were found for transitional waters, especially those situated in the north. According to the Croatian Meteorological and Hydrological Service, mean monthly air temperatures for most of 2018 were above the average of the reference period 1961-1990, while in May 2019, temperatures were below the reference average along the entire coast. Specially in the area gravitating toward T_N, temperatures in May 2019 were nearly 5°C below average (DHMZ, 2019). Sites within this group are related to the transitional waters of the Zrmanja and Krka Rivers, whose water systems are more exposed to the mainland due to relief characteristics (heavy rains, snow melt from nearby mountains and cold north-easterly winds in spring) and the maritime influence is less pronounced. Usually, salinity fluctuations are associated with freshwater river flows, though the significantly lower temperatures of these waters, especially in the spring months are due to their origin and flow through extremely mountainous terrain. On the contrary, the transitional waters of the Neretva and Cetina Rivers have wider mouths toward the open sea, especially the Neretva Delta, and thus are more influenced by the gentle maritime environment (Krvavica and Ružić, 2020). Unlike temperature, there were no prominent interannual fluctuations in salinity. However, reduced salinity was evident in the northern transitional waters for the same reasons as mentioned above. In late spring 2019, reduced salinity was recorded along the entire coast as a consequence of an extremely rainy period. For example, the amount of precipitation in May (120 mm) and November (210 mm) 2019 in the middle of study area (S8 and S9) was twice the average cumulative amount of precipitation in the reference period 1961–1990 for those months (DHMZ, 2019). In 2017 and 2018, the average salinity in the south, regardless of water type, was always higher than in the north, confirming the unusually high rainfall in the north (DHMZ, 2019). As expected, statistically significant interannual differences in mean seasonal water temperatures also led to a change in the species occurrence among the three consecutive years. Thus, June 2019 was characterised by markedly low salinity and low temperature after a cold and rainy winter and spring, which likely caused the absence of the following species: B. boops, B. podas, P. acarne, P. erythrinus, S. tinca, S. roisalli, and C. auratus. All these species, except C. auratus, spawn in spring (Bartulović et al., 2011; Dulčić and Kovačić, 2020) and their settlers likely delayed settlement due to the prevailing conditions. Interestingly, a significantly reduced abundance of B. boops and B. podas were observed in southern transitional waters (T_S) in the same period. These species also spawn in early spring (Dulčić et al., 2004; Matić-Skoko et al., 2007; Dobroslavić et al., 2017; Zorica et al., 2020) and apparently their juveniles did not yet inhabit the nursery areas by June. Late spawning and a longer retention in the nursery area could be the reason for increased abundance of L. mormyrus in September 2019. Increased abundance is most often a positive response to average warmer habitat conditions (Jones et al., 2009). It seems that the longer the spawning season for a population, the higher its probability of finding suitable conditions for larval survival, resulting in a higher recruitment level (Sale, 2006). This may imply that if environmental conditions change significantly, species with a shorter spawning period may experience a more pronounced shift in spawning time and consequently in settlement onset. Interannual variability in temperature has been linked with changes in fish year-class strength for a number of species (Planque and Fox, 1998; Laurel et al., 2017; Barbeaux and Hollowed, 2018). Temperatures affect the onset of spawning and spawning period, but also the ecology of early developmental stages of fish, thus defining the moment of their migration or entry into benthic nursery areas (García-Rubies and Macpherson, 1995; Dulčić et al., 1997, 2007; Biagi et al., 1998). Also, positive relationships between temperature and impacts on settlement and recruitment are generally observed in cold-temperate regions, whilst negative relationships are observed in warm-temperate regions (Planque and Fox, 1998; Barbeaux and Hollowed, 2018; Downie et al., 2021). Environmental changes at higher latitudes that are under greater coastal influence appear to be more drastic or more rapidly reflected in changes at the community, species or individual level.

The influence of environmental factors on species composition or occurrence is specifically evident at the fish community level for both species richness and abundance. Moreover, water temperature appears to have a greater impact on transitional waters that are often retracted deep into the mainland. This is even more relevant since more than 64% of fish juveniles representing 68 species and 28 families were recorded in estuarine nurseries, and most are species of commercial interest (Dulčić et al., 2007). On the other hand, marine waters appear to be a less volatile environment, and the establishment of the juvenile fish community was defined primarily by salinity. This is supported by the fact that the same species were recorded in marine waters in all three years. Dulčić et al. (1997) investigated several shallow marine coves on distant islands and found a similar temporal pattern and the prevalence of five dominant species, suggesting that only a low amount of variation in the juvenile fish abundance in the real marine environment can be explained by temperature and salinity. The absence of certain species (e.g., S. roissali and C. lucerna) in 2019 in both marine and transitional waters could be regulated by a factor that is more regional in nature. Similarly, Stagličić et al. (2011) reported that species abundance within littoral fish assemblages could be a consequence of factors occurring on a wider spatial scale (i.e., fisheries mismanagement, pollution, climate change). Contrary, C. saliens and S. salpa showed a higher abundance in northern marine waters (M_N) in 2019 and 2018, respectively, while S. cinereus was most abundant in T_S in 2018 in all four groups. However, there is the possibility that the sites where certain species were recorded are not essential or preferable nurseries for those species (i.e., poor site selection, low sample size, sampling gear selectivity). For example, P. erythrinus or T. ovatus were both lacking in southern marine areas in both 2018 and 2019. These species use other nursery areas, like marine coves on distant islands or at greater depths that are characterised by high and constant salinity (Dulčić et al., 1997; Biagi et al., 1998; La Mesa et al., 2011). In this study, P. erythrinus was found in 2017 exclusively at one station (S12) that is wide open toward the sea. Another typical marine species, P. pagrus, was also caught only at this site. The dispersal strategies of certain small migratory, shoaling species, such as Spicara sp. or B. belone, could also be the reason for their spatio-temporal absence (Hoare et al., 2000).

If transitional waters geographically located at higher latitudes are more susceptible to interannual change, then the assumption is that more vulnerable species will select those waters for their nurseries, especially the more abundant species in the northern-central Adriatic (Dulčić et al., 2004; Matić-Skoko et al., 2007; Dulčić and Kovačić, 2020). Firstly, species of economic interest such as D. labrax, S. salpa and S. aurata should be highlighted here. It is also important to mention the small, resident species that also have ecological significance in the benthic communities of transitional waters, e.g., Salaria pavo, S. ocellatus, and Z. ophiocephalus, as these herbivorous and omnivorous fishes seem to control both algae and invertebrates in rocky sub-littoral zones (Ruitton et al., 2000). S. ocellatus showed a significant association with an algal habitat, and its abundance can certainly be related to the strong seasonal differences in algal cover (Hinz et al., 2019). If species preferring transitional waters at higher latitudes also spawn in winter and/or early spring, such as D. vulgaris, D. sargus, S. aurata and D. labrax, this makes them particularly sensitive to interannual differences in environmental conditions. As mentioned earlier, they can potentially be affected by late spawning, causing delays in settlement, leaving nurseries, or even causing individuals to remain in nurseries year round. Particularly, this delay can lead to the fact that at the moment of entry there is no prey of suitable size or in sufficient quantity for all incoming settlers (Machado et al., 2017). Decreased growth can make individuals too weak to join adult populations in the open seas that season (Lanier and Scharf, 2007). Since many marine species undergo long-distance dispersal during the pelagic larval phase before settling in benthic habitats where they remain (Di Franco et al., 2012), and many species also undergo further ontogenetic changes in restricted nursery grounds, pre- and post-settlement emigration by juveniles can then have significant population-level consequences (Vasconcelos et al., 2013). Sensitivity to physical-chemical changes of the aquatic environment, lack of suitable prey, and predation are the main causes of juvenile fish mortality (Guidetti, 2001; Cuadros et al., 2018).