Population structure of European sprat ( Sprattus sprattus ) in the Greater North Sea ecoregion revealed by otolith shape analysis

A successful discrimination of fish populations is essential for sustainable management and assessment. Otolith shape analysis has been used on several species to reveal their population structure. Our aim was to use the otolith shape of European sprat ( Sprattus sprattus ) to investigate large-and small-scale geographical variability across the Greater North Sea ecoregion. The otolith shape was extracted from digitalised images and transformed into Wavelet coefficients to be analysed with multivariate statistics. Otolith shape was observed to follow the genetic population structure recently defined for the region, supporting the latest revision of the stock boundaries. Four main groups were identified based on phenotypic variability: (i) Norwegian fjords; (ii) North Sea and offshore Skagerrak – Kattegat; (iii) coastal Skagerrak – Kattegat; and (iv) Uddevalla fjord. However, 4-fold cross-validations based on Linear Discriminant Analysis resulted in low accuracy limiting at the moment the ability to use otolith shape analysis for population identification at an operational basis. Our results show the importance of coastal areas, which might be inhabited by distinct populations of sprat that are currently not acknowledged in the management and assessment.


Introduction
An accurate understanding of the population structure underlying a fish stock is of vital importance for sustainable management.Misconceptions of the relative contribution of different components may lead to wrong estimations of the fishing pressure that the stock may sustain and, consequently, to suboptimal utilization of the resource (Kerr et al., 2016).A prerequisite for investigating population structure in marine fish species is a successful discrimination of populations, either based on genetic or phenotypic differences.In recent years, genetic studies have revealed detailed population structure even for species with high gene flow like pelagic fish (Gíslason et al., 2020;Martinez Barrio et al., 2016).Traditionally, population discrimination was achieved by applying a variety of phenotypic methods such as growth analysis (Gröhsler et al., 2013), meristic characteristics (Meng and Stocker, 1984;Misra and Bowering, 1984;Turan, 2004), or otolith analysis investigating microchemistry (Geffen et al., 2011), microstructure (Brophy and Danilowicz, 2002), shape (Agüera and Brophy, 2011;Lee et al., 2018;Libungan et al., 2015a) or stable isotopes (Schade et al., 2019).Also, tagging experiments (Block et al., 2005), assuming different behavioural traits can be applied to discriminate populations.Among phenotypic methods, otolith shape analysis is a particularly powerful tool to investigate population structure, because the otolith shape is controlled by both genetic and environmental factors (Berg et al., 2018;Cardinale et al., 2004;Vignon and Morat, 2010), and thus highly variable between species and populations (Lombarte and Lleonart, 1993).Along with genetic affiliation, body size and age (Campana and Casselman, 1993), temperature (Gagliano and McCormick, 2004), diet (Mille et al., 2016), salinity (Berg et al., 2018), and habitat depth (Lombarte and Lleonart, 1993) have been shown to be important variables affecting otolith shape, making it a successful marker for stock discrimination in several fish species (Libungan et al., 2015a;Stransky et al., 2008).
The European sprat (Sprattus sprattus, Linnaeus 1758) is a small-bodied, pelagic schooling clupeid (Peck et al., 2012) inhabiting the Baltic, the Northeast Atlantic from Northern Norway down to Morocco, the Adriatic Sea in the Mediterranean, and the Black Sea (Debes et al., 2008).As other main forage fish, sprat plays an ecologically important role in the ecosystem as both planktivore and prey, and it is one of the main food sources in the diet of numerous species including demersal fish, seabirds and marine mammals (Lundström et al., 2010;Lynam et al., 2017).The production of pelagic eggs, and the opportunistic vagrant behaviour, suggest a potential for high dispersal and gene flow (Limborg et al., 2009), as found for instance in the Baltic Sea where the stock is composed of a single large population (Shvetsov et al., 1995).Sprat throughout the Greater North Sea ecoregion (Fig. 1), including the Skagerrak and Kattegat (ICES division 3a) and with the exception of the Norwegian fjords, are currently considered and managed as one stock (ICES, 2018), also comprising a single genetic population (McKeown et al., 2020;Quintela et al., 2020).Quintela et al. (2020) highlighted the presence of two additional genetic populations occurring in the Norwegian fjords and the Baltic Sea.The geographic boundaries with the neighbouring Norwegian and Baltic stocks remain uncertain (ICES, 2018) and potential mixing or hybridization cannot be excluded.Patterns of hybridization have been previously demonstrated by Limborg et al. (2009), where the genetic differentiation among sprat was reported to reflect the gradient in salinity between samples from the Kattegat and the Baltic Sea.
Despite the marine environment may seem to offer fewer physical barriers compared to terrestrial environments, many fish species have been proven to display population structures reflecting obstacles to gene flow over relatively small geographic distances (Limborg et al., 2012).The complexity of the Scandinavian fjord system is a potential source of physical barriers between spatially close populations (Libungan et al., 2015b), and local aggregations of sprat in the Skagerrak and Kattegat being morphologically or ecologically adapted are known (Lindquist, 1968;Molander, 1940).One location of particular interest is the Uddevalla fjord system along the Swedish west coast (hereafter referred to as "Uddevalla fjord").This is a peculiar environment composed of several different fjords separated by small straits and shallow sills (Hansson et al., 2013), with two sills 8-and 20 m deep where the two main branches meet the Skagerrak (Gustafsson and Nordberg, 2000).Sprat from this area show distinct size and age at maturity, slower growth and high rates of parasite infestation ("black spot disease" transmitted by the trematode Cryptocotyle lingua; Vitale et al., 2015).Furthermore, genetic divergence occurs between sprat found in the Uddevalla fjord and the Kattegat (Quintela et al., 2020), corroborating the existence of a relatively isolated fjord population.
To our knowledge, it has not previously been studied whether the otolith shape of sprat reflects their geographic distribution.In addition, the availability of a genetically validated population structure constitutes a strong baseline to investigate differences in otolith shape Fig. 1.Greater North Sea ecoregion with sampling locations of European sprat used for otolith shape analysis.Main map shows locations of samples targeted during 2014-2018 (dataset I).Inset shows locations of samples targeted during 2003-2004 (dataset II) (dataset II).Rectangles represent samples from fjords, triangles from coastal areas and circles from offshore areas.Skagerrak (green), North Sea (black), Norwegian fjords (pink), Kattegat (blue), Gullmar fjord (yellow), and Uddevalla fjord (red).
F. Saltalamacchia et al. between geographic areas.This study, therefore, aims to (a) evaluate whether otolith shape analysis supports the genetic-based interpretation of population structure which recently led to a revision of the sprat stock boundaries with the North Sea, Kattegat and Skagerrak into a single stock area, (b) investigate otolith shape differentiation on a small geographical scale by comparing coastal and offshore areas in the Skagerrak and Kattegat off Sweden.To address these aims we applied otolith shape analysis on sprat from the Greater North Sea ecoregion, also including samples from Norwegian fjords.

Data collection
Stored dry otoliths of sprat were available and retrieved from collections belonging to the Department of Aquatic Resources of the Swedish University of Agricultural Sciences, the National Institute of Aquatic Resources of the Technical University of Denmark, and the Norwegian Institute of Marine Research in Bergen.Otoliths were collected in the early 2000 s and the mid-2010 s, covering different areas of the Greater North Sea ecoregion (Fig. 1).After measuring total length the fishes were sexed and the sagittal otoliths extracted.The otoliths were used independently by two scientists for individual age determination (winter rings counting) and stored dry in plastic bags until image capturing.
Since the sampling was inconsistent between the two periods across different areas, the two decades were analysed separately and with different intent, to avoid the potential confounding effect of time on the otolith shape.The first dataset included otoliths of sprat collected between 2014 and 2018 from Norwegian fjords (NW), North Sea (NS), offshore and coastal Skagerrak (SK and CS, respectively), and offshore Kattegat (KA) (Fig. 1), and it was used to address the first aim and to investigate the large-scale stock structure.This dataset comprised of 548 individuals ranging from age 1-3 (Table 1).The range of age classes was selected with the intent to minimise ontogenetic influences on shape while maintaining a sufficiently large sample size for analytical purposes.Many samples could not be retrieved over spawning season, and therefore this fish may represent both local and external components.
The second dataset consisted of samples collected in 2003-2004 from Skagerrak andKattegat (Fig. 1Fig. 1, inset) and was used to investigate the small-scale population structure of sprat in this region.These samples were collected along the coastal Skagerrak (CS), coastal Kattegat (CK), offshore Kattegat (KA), and in the fjords of Uddevalla (UF) and Gullmar (GF).In total, 860 individuals of age from 1-3 were analysed (Table 1).In this case, sampling was performed between late January and early July, thus including the spawning time of the species in the region (Vitale et al., 2015).

Image processing
Only otoliths that did not show any signs of damage were used in this study.The analyses were carried out in R (R Core Team, 2017) following the procedure illustrated by Libungan and Pálsson (2015).The otoliths were positioned sulcus down on a dark background, with the ventral side facing upwards.High-contrast images were taken with a LEICA DFC320 digital camera mounted on a LEICA stereo microscope MZ16 FA, using the software LEICA Application Suite 4.1.0(build 1264).The otolith outlines were extracted from the images using the "shapeR" package (Libungan and Pálsson, 2015).The detected outlines were overlaid on the original pictures and visually evaluated to ensure that each outline narrowly traced the edge of the otolith.When the outline was not accurate, the original picture was re-elaborated with a different threshold value, or manually edited with GIMP (GNU Image Manipulation Program v. 2.8.22).
All the outlines were smoothed to reduce pixel noise.To account for variations in size, rotation, and position of the otoliths, the outlines were automatically set with area = 1, rotated so that the longest Feret axis was positioned horizontally, and centred on their centroid.Radii were drawn at equal angular intervals from the centroid to several coordinates along the outline, and used to calculate 64 independent wavelet shape coefficients by Discrete Wavelet Transform.

Otolith shape analysis
The following analyses were performed for each of the two datasets individually.Directional asymmetry (left-right side), sex, and body size are known to be potential confounding factors when studying otolith shape (Castonguay et al., 1991;Ider et al., 2017), and were thus investigated.To account for interactions between fish size and geographic area on the otolith shape, an Analysis of Covariance (ANCOVA) was performed on the wavelet coefficients.Bonferroni correction was applied to compensate for the increased likelihood of Type I error due to multiple testing of the variables.The coefficients that showed significant interactions (P < 0.05) between body size and geographic area (grouping factor) were discarded.The remaining were automatically adjusted for allometric relationships with fish size by using a common within-group slope normalization technique (Lleonart et al., 2000).The standardized wavelet coefficients were visually inspected for normality before further statistical analyses.
When available (approximately 90% of the individuals), the right otolith was selected for the analyses.While no differences have been found between left and right otoliths of Baltic sprat (Aps et al., 1991), no studies investigated this matter on sprat from the Skagerrak-Kattegat area.Since the original datasets only included one otolith per individual, a new dataset of 160 otolith pairs (left and right) was built and used only for the purposes of this check.The test for variations in otolith shape between male and female individuals was performed on the original datasets.No significant differences were observed between left and right otoliths (t-test, P > 0.05), nor between sexes with respect to sampling location (ANOVA, P > 0.05).Consequently, left otoliths were flipped horizontally and included in the analyses in cases where the right one was not available, and the analyses were performed on merged data for

Table 1
Details of sprat samples used in the study for each dataset.ID represents the abbreviation of each sampling location.Note that not all months were sampled in all years.Mean total length (TL, cm) ± standard deviation (sd) as well as sampling size (n) per group are given.males and females.Canonical Analysis of Principal coordinates (CAP) using the capscale function (vegan package; Oksanen et al., 2013) was performed on the standardised wavelet coefficients to explore the relationships between otolith shape and geographical area.The canonical scores were further tested for significance (α = 0.05) using ANOVA-like permutation tests with 2000 permutations.

Classification
Individual assignment of the otoliths in respect to their geographical areas was evaluated using a Linear Discriminant Analysis (LDA) on the standardised wavelet coefficients by means of functions written in the caret package (Kuhn, 2008).LDA is a supervised method used to discriminate among predefined groups of individuals based on a sample of observations from the original groups.The classification success rate was estimated using 4-fold cross-validation with 1000 repetitions, which omits one-quarter of the observations at a time (test set) and then classifies them using a classification model built on the rest of the observations (training set).Since the discriminative power of LDA is known to be affected by highly unbalanced class sizes, a down-sampling procedure was applied to ensure that the number of individuals in each class matched the abundance of the least conspicuous one.

Norwegian fjords -North Sea -Skagerrak/Kattegat
Otolith shape compared among sprat from the Norwegian fjords (NW), North Sea (NS), offshore Skagerrak (SK), coastal Skagerrak (CS), and offshore Kattegat (KA) by CAP and ANOVA-like permutation tests showed significant differences between areas (P < 0.05, Fig. 2).The first two canonical axes explained approximately 90% of the shape variations between regions (71.7% and 17.9% for CAP1 and CAP2, respectively).The relative positions of the group centroids along CAP1 discriminated the NW samples from the rest of the samples, which remain closer to each other along this axis.No significant differences were found between NS vs. SK and SK vs. KA (P > 0.05), whereas the differences for NS vs. KA were significant (P = 0.01).CS differed slightly from all other areas (P < 0.05).An LDA with 4-fold cross-validation was performed to check whether otolith shape could be used to classify individuals back to their area of origin, resulting in an overall success rate of 43% (Table 2).

Kattegat/Skagerrak
CAP and ANOVA-like permutation tests on small-scale perspectives along the coastal Skagerrak (CS), coastal Kattegat (CK), offshore Kattegat (KA), and the two fjords of Gullmar (GF) and Uddevalla (UF) show that otolith shape varies significantly among areas (P < 0.05, Fig. 3).Around 90% of the variability is explained by the first two canonical axes (69.1% and 19% respectively), with the greater divergence occurring between the UF and all the other samples (Fig. 3).KA otoliths differed significantly from UF, GF, and CS (P < 0.05).Pairwise comparisons showed no significant differences between CS and CK (CAP, P > 0.05), nor between CS and GF.The overall classification accuracy based on LDA with cross-validation was 26% (Table 3).

Discussion
Otolith shape analysis of European sprat showed clear patterns of differentiation among geographic locations in the Greater North Sea ecoregion, suggesting that shape variability has the potential to be used to explore the population structure of this species.However, the accuracy obtained by LDA on our data was too low to allow successful discrimination of individuals based on the shape coefficients only.Four distinct groups were detected: (i) Norwegian fjords; (ii) North Sea and offshore Skagerrak-Kattegat; (iii) coastal Skagerrak-Kattegat including Gullmar fjord; and (IV) the Uddevalla fjord.Although the underlying causes of the shape variability were not examined, these results are congruent with recent genetic studies (McKeown et al., 2020;Quintela et al., 2020), and show that otolith shape follows the genetic population structure defined for this area.Similar agreements between otolith shape and genetic methods have been observed for several fish species  (e.g.Kristoffersen and Magoulas, 2008;Libungan et al., 2015a;Valentin et al., 2014).An exception seems to be represented by the sprat from the Swedish coastal areas, that have not yet been investigated by other methods but whose otolith shape seems to depart from the observed "North Sea -Skagerrak -Kattegat" clustering pattern.Among the sampled locations, the largest differentiation was detected between the Norwegian fjords and the North Sea -Skagerrak -Kattegat.Previous genetic studies on microsatellite DNA loci (e.g.Glover et al., 2011) observed that sprat from the Norwegian fjord are highly different from the North Sea, Celtic Sea, and Baltic Sea populations, indicating little or absent movement of individuals into and from the Norwegian fjords.The occurrence of reproductively isolated components of sprat within the western Norwegian fjords had already been suggested by Naevdal (1968) based on the observation of hemoglobin and serum proteins, and small-scale patterns of differentiations have been reported for other fish species such as cod and herring (Knutsen et al., 2003;Bekkevold et al., 2005).
The samples from the North Sea and offshore Skagerrak-Kattegat did not significantly differ, showing complete accordance between morphological and genetic methods.In this respect and considering the recent scientific advice (ICES, 2018), corroborating evidence from independent methodologies is highly valuable to support the revision of boundaries among assessment units, which should mirror as much as possible the underlying biological components.In this study, the similarities between the North Sea and Skagerrak-Kattegat were assessed within the age classes 1-3, since these were the only classes available from the North Sea samples.It would therefore be advisable to conduct further studies on otolith shape variations between the two regions by sampling older ages in the North Sea and more locations in the offshore Skagerrak.
As expected, given the lack of known physical or environmental barriers that could prevent gene flow, no differences were found among coastal samples from Skagerrak and Kattegat, with the exception of the Uddevalla fjord.Otolith shape was in fact observed to be homogeneous along the Swedish coast, while some differentiation was found between the offshore and coastal areas both in Skagerrak and Kattegat.This suggests that a coastaloffshore differentiation might exist along the Swedish west coast, similarly to what has been inferred by genetic studies for southern Norway (Quintela et al., 2020).Additionally, it would be of interest to assess whether otolith shape reflects the mixing between adjacent populations that has been observed by means of genetics throughout the southern Kattegat and the Belt-Baltic Sea (Limborg et al., 2009;Quintela et al., 2020).To date the level of such mixing, as well as its seasonal variability, remains unknown.
A clear differentiation was observed between the Uddevalla sample and the rest of the coastal Skagerrak-Kattegat, confirming what was shown by previous studies based on morphometrics (Lindquist, 1968), life history traits (Vitale et al., 2015), and genetics (Quintela et al., 2020).Due to their geographical characteristics, the Uddevalla fjord constitutes a sheltered environment relatively isolated from the Skagerrak (Nordberg et al., 2001).This area is highly influenced by low saline Baltic waters (Björk et al., 2000), and it is characterised by lower water temperature during summer and local conditions of hypoxia and anoxia (Vitale et al., 2015).By comparing morphometric and meristic characters among samples from Skagerrak and Kattegat, Lindquist (1968) observed a lower number of vertebrae in the Uddevalla sprat, similar to the number reported for the Black Sea and Adriatic Sea populations.The work carried out by Quintela et al. (2020) showed a clear genetic divergence between sprat in the Uddevalla fjord and the outer sea.In contrast, the otolith shape of sprat from the inner part of Gullmar fjord (just a few kilometres north along the Swedish Skagerrak) did not differ from the rest of the coastal samples.Compared to the Uddevalla fjord, Gullmarn shows a less complex structure, with a deeper sill at − 42 m (Nordberg et al., 2000) separating the fjord from the Skagerrak; nonetheless, local mechanisms of differentiation in the inner part of the fjord have been demonstrated for other species such as cod, for which a separate spawning component was identified (Øresland and André, 2008).In contrast to cod, which are more stationary, sprat could still migrate and interbreed with individuals from outside the fjord.
The observed phenotypic variability in the otolith shape of sprat seems to capture population structure at both broad and small geographical scales.However, the low accuracy shown by the LDA needs to be considered, since it poses limits to an operational use of the method.There are a number of factors that could have negatively  impacted on the discriminative power.An operational algorithm for stock discrimination relies on clean baselines of stock-specific otolith shapes built from individuals of known origin (Hüssy et al., 2016).One of the basic assumptions in morphometric studies is the use of samples collected during spawning time and, possibly, of spawning (or close to-) individuals (Campana and Casselman, 1993), to make sure that the underlying population structure is not masked or confused by the mixing of different biological components.However, the fish used for this study was collected over different seasons, and maturity stages could not be considered.This may bring potential detrimental effects on the strength of the signal in the data.Moreover, the combination of various age classes, even though based on a small range and with data corrected for fish size, may reduce the model ´s ability to capture peculiarity between different locations (Mapp et al., 2017).It is likely that the discriminatory power could be improved by extending the analyses to a larger sample set better conformed to these assumptions.Finally, the possibility has to be considered that otolith shape of sprat, even though showing clear patterns of differentiations, does not vary enough to provide clear gaps between groups that could be successfully used by classification algorithms.Despite LDA has shown high accuracy rates for both cod and herring, other discrimination methods like machine learning algorithms could have led to slightly higher assignments rates (Smoliński et al., 2020).
To evaluate the use of geometric morphometric methods in population studies, it is always important to validate the results by means of other independent techniques and, when available, genetic references are preferred (e.g.Afanasyev et al., 2017;Hüssy et al., 2016).The general patterns shown by this study appear consistent with biogeographical considerations, previous knowledge on sprat population structure (Limborg et al., 2009(Limborg et al., , 2012)), and the most recent genetic analyses (Quintela et al., 2020), showing that this method has a potential to be applied on European sprat population studies and deserves further attention to assess its operational purposes.Considered the various limits posed by the nature of the samples, these results appear promising despite the low classification accuracy performed by the LDA.Otoliths are routinely collected for age reading and can be stored for many years, allowing the analysis of long time series.Moreover, otolith shape analysis is a non-destructive, relatively quick and inexpensive method; if proven operatively successful in terms of discrimination power, this technique could be easily applied to the determination of geographical boundaries between stocks and to address issues of stock mixing, as done for several other species such as herring (Libungan et al., 2015b), anchovy (Zengin et al., 2015), mackerel (Turan, 2006) and red snapper (Sadighzadeh et al., 2014).Its implementation in a comprehensive, multidisciplinary approach, should therefore be evaluated.

Fig. 2 .
Fig. 2. Canonical scores on the first two discriminating axes for each sampling location in the Greater North Sea ecoregion sampled in 2014-2018 (dataset I).Coastal Skagerrak (CS), offshore Kattegat (KA), North Sea (NS), Norwegian fjords (NW), offshore Skagerrak (SK).The intervals around the group centroids represent the 95% confidence interval for the standard error of the mean.

Table 2
Cross-validated classification matrix of Linear Discriminant Analysis between the sampled locations in the Greater North Sea ecoregion sampled in 2014-2018 (dataset I).Norway (NW), North Sea (NS), offshore Skagerrak (SK), offshore Kattegat (KA), coastal Skagerrak (CS).The values are given as percentages of the total.Bold values indicate the percentage of correct classification on the total count of individuals.

Table 3
Cross-validated classification matrix of Linear Discriminant Analysis between the sampled locations from Skagerrak-Kattegat sampled in 2003-2004 (dataset II).Coastal Kattegat (CK), coastal Skagerrak (CS), Gullmar fjord (GF), offshore Kattegat (KA), Uddevalla fjord (UF).The values are given as percentages of the total.Bold values indicate the percentage of correct classification on the total count of individuals.