Assignment of Gilthead Seabream Sparus aurata to Its Origin through Scale Shape and Microchemistry Composition: Management Implications for Aquaculture Escapees

: This study tests the suitability of the gilthead seabream scales as a proxy for origin selection in wild and anthropogenically pressured environments. Scale morphology and microchemistry were used to discriminate the habitat selection of two wild, farmed and wild farm-associated populations where landmark and outline-based scale morphometrics, trace-element chemistry and scale microstructure characteristics were analysed. The morphometric techniques successfully di ﬀ erentiated between the farmed and wild origin scale phenotypes. Reduced discrimination sensitivity between the wild and wild farm-associated origin was, however, reported. The discrimination based on microchemistry (B, Ba, Mn, K, Sr and Zn) classiﬁed the scales with high accuracy according to their origin (wild vs. farmed vs. wild farm-associated) and sampling locations, thus proving itself as a powerful tool in provenance study of gilthead seabream. Disparity in scale microstructure characteristics accounted for radii, circuli and inter-circulus spacing, hence unveiling the di ﬀ erences in growth and environmental conditions between the wild and farmed ﬁsh. In brief, scale shape was found to be a potent exploration tool for farmed ﬁsh identiﬁcation, whereas scale microchemistry yielded a good resolution in identifying gilthead seabream membership among di ﬀ erent habitats. Considering the importance of this species in aquaculture and ﬁsheries throughout the Mediterranean, more research is needed to assess the usefulness of scales as nonlethal biogeochemical tags.


Introduction
Population connectivity and spatial structure information provide a basis for understanding fish population dynamics and play a key role in conservation and fisheries management [1]. In a spatially heterogeneous environment, high gene flow may actually favour the evolution of increased phenotypic plasticity over adaptive genetic divergence between populations [2]. Individual phenotypical adaptation of morphological, behavioural and physiological traits to new conditions may occur within one or two generations, allowing them to persist in novel environments. Several studies have activities in Croatia. Two farmed populations having Atlantic brood-stock origin were sourced from the above-mentioned locations, from two of the largest commercial grow-out farms (FS1, n = 56; FS2, n = 100), each with an annual production capacity of >100 tonnes. Selected wild sampling locations (WS1, n = 92; WS2, n = 100) were distant from all aquaculture facilities (70 km distance for WS1, 35 km for WS2), where escapees were not expected ( Figure 1). Sampling of wild mature adults was primarily conducted in collaboration with small-scale fisheries and all fish collections used in the study were sampled during the spawning season, from December 2015 to February 2016. Total length (cm) and weight (g) were recorded for each individual. All specimens were adults with average total lengths (TL ± SD) of 26.7 ± 4.2 cm, 26.5 ± 1.8 cm, and 26.4 ± 2.9 cm for the wild, farmed and wild farm-associated collections, respectively. No significant difference was found for TL between sampling locations (p > 0.05). Table 1. Summary of sampling locations from the middle Adriatic Sea, along with the origin and number of individual gilthead seabream Sparus aurata used in the present study. The number of fish scales per population (n) used in the landmark and outline-based morphometrics, micro-structure and trace element analyses are also presented.

Sample Collection and Preparation
Scales were removed from an area above the lateral line, between the posterior insertion of the dorsal fin and the anterior insertion of the anal fin and were stored in manila envelopes. All equipment in contact with scale samples was plastic, acid-washed in 10% nitric acid (TraceMetal Grade, Seastar Chemicals Inc., British Columbia, Canada), triple rinsed in Milli-Q water (Millipore, Billerica, MA, USA), dried in a laminar flow hood, and stored in sealed plastic bags prior to use. In the laboratory, the fish scale cleaning process included sonication for 5 min in 3% ultrapure hydrogen peroxide to loosen organic debris followed by brushing with a toothbrush and washing with Milli-Q water. The procedure was repeated twice.
A portion of the decontaminated scales was used for elemental analyses, while a single undamaged scale from each fish was selected for geometric morphometric analysis. The cleaned and dried scale was examined and photographed under 40x objective lens of the Dino-Lite AM311 digital microscope (AnMo Electronics Corp., Hsinchu, Taiwan) for further image-processing analysis ( Figure 2). In addition, three to five scales from each fish were cleaned, microscopically examined and the number of annuli was counted along the longest antero-posterior axis on the scale for fish age determination. Farm fish were excluded from the age analysis. Scale surfaces were also examined by scanning electron microscope (FEI Quanta 650F, SEM, Thermo Fischer Scientific, Waltham, MA, U.S.A.) to provide more information on the morphology and biomineralized structures present in fish scales. U.S.A.) to provide more information on the morphology and biomineralized structures present in fish scales.  Table 1.  Table 1.
Water 2020, 12, x FOR PEER REVIEW 4 of 17 U.S.A.) to provide more information on the morphology and biomineralized structures present in fish scales.  Table 1. . Scale features as primary radii (R), inter-radial circuli (IRC) and area used for circuli counts and distance measurement (CC) in the rostral field are also presented.

Morphometric Analysis
Landmark-based and outline-based approaches were used to identify scale morphometric relationships among populations of different origins and geographic sites. No allometric adjustments were conducted since all specimens used in the present study originated from the same age (2 years) and length classes.
Landmark-based geometric morphometrics (GM). The total dataset was used for GM analysis based on seven-landmark configuration, following the Ibáñez [13] landmark digitalization scheme ( Figure 2). Shape variation was analysed in the software MorphoJ 2.03b [19] using a full Procrustes fit to remove variation in scale, position and orientation. Principal components analysis (PCA) was applied to the Procrustes coordinates while canonical variate analysis (CVA) was performed to test group pairwise differences and produce ordinations. The jack-knife cross-validated classification table was performed in PAST ver. 3.0 software [20].
Elliptic Fourier analysis (EFA). Scale images were converted into black and white silhouettes for the procedure of outline extraction in the "Momocs" R package [21], where they were transformed into a list of (x, y) coordinates for each pixel around the contour of a given scale shape, and turned into quantitative descriptors by Elliptic Fourier analysis (EFA). Harmonics were added until at least 99% of the variance in the scale shape was reconstructed [22] and Fourier descriptors were obtained for further multivariate statistical analysis. To demonstrate classification of individuals into their sampled origin and location, a linear discriminant analysis (LDA) was performed.

Trace Element Analysis
Trace elements in scales from 96 selected wild, farmed and wild farm-associated seabream (Table 1) of similar total length (25 ± 0.8 cm) and weight (250 ± 73 g) were analysed using inductively-coupled plasma mass spectrometry (ICP-MS) at the environmental laboratory Maxxam Analytics (Montréal, Québec, Canada). The concentrations of Ba, B, Cu, Fe, Mg, Mn, K, Sr, Ti, Zn, Hg and Pb were measured in paired samples. Quality control was performed regularly (blank, certified reference material, duplicate, and fortified blank). The correlation coefficient for duplicates was 0.999. Results obtained for fortified blanks were within the limits set by the laboratory. Only those elements whose ICP-MS concentrations were above the respective detection limits were considered suitable for statistical treatment. Hg and Pb were not detected in all samples and were therefore omitted from the reported dataset. Finally, a total of 10 trace elements was selected for further analysis.
Fourth-root transformation of elemental concentration data was performed to stabilize the homogeneity of variance. To investigate the extent to which element concentrations in scales varied by origin, Kruskal-Wallis one-way analysis of variance range tests were used followed by Dunn's post hoc multiple comparison tests. Two-way permutational multivariate analysis of variance (PERMANOVA) was further applied to the Bray-Curtis resemblances of elemental concentration data to detect differences in the elemental compositions of scales in individuals of different origin and location. Origin (three levels: wild, farmed and wild farm-associated) and Location (two levels: Site1 and Site2) were fixed and orthogonal factors. To inspect variance contribution sources, sizes of the multivariate variance components obtained using mean squares from PERMANOVA were compared between the main factors [23]. Significant terms (p < 0.05) of the main PERMANOVA test were examined further by appropriate post-hoc pair-wise comparisons. Finally, canonical analysis of principal coordinates (CAP; [24]) was used to visualize multivariate patterns of scale trace element composition and to classify individuals according to their natural origin. Classification accuracies were evaluated through the percentage of correctly classified individuals using the leave-one-out cross-validation method, where an allocation success of 12.5% would be expected by chance alone (six groups). Elements contributing most to differences among groups were identified by examining Spearman rank correlations between the elemental concentrations and primary CAP axes where size and direction of the vector corresponds to the strength and sign of the relationship, respectively [23]. Relationships between the individual element and canonical axes with a correlation coefficient (r) > 0.30 are displayed on the CAP ordination as a vector overlay. In addition, combined trace elements and scale shape based on Fourier descriptors data was obtained from 96 individuals of known microchemical signatures to evaluate allocation success when two different methodologies were combined. The combined data was normalized prior to CAP analysis.

Microscopic Analysis
From the selected dataset (Table 1), intraspecific variation in scale microstructures of the rostral region was examined from the images obtained by a digital microscope, focusing on number of primary radii (R), the number of circuli and the distance between them in the inter-radial area. For each scale, circuli counting and measuring distances were restricted to the area between two central grooves that radiate from the focus to the scale anterior margin, including only 0.5 millimetres of the scale radius ( Figure 2). This scale-edge circuli approach also enabled farm fish to be included in analysis, since the majority of farmed fish had regenerated nuclei and deformed circuli in the focus vicinity [15]. Analysis was based on 350 selected wild, farmed and wild farm-associated seabream scales (Table 1) using ImageJ 1.52 software (http://imagej.nih.gov/). The Kruskal-Wallis test was employed to compare scale characteristics among tested fish origins. Location comparisons were omitted from analysis due to the highly unbalanced number of samples between sites. The analysis was followed by Dunn's multiple comparison tests if a significant difference (p < 0.05) was obtained.
To qualitatively compare regional scale morphology, SEM imaging was conducted on several scales per population. Freshly fractured surfaces of scales were fixed to pin stubs using conductive double-sided adhesive carbon tabs and then carbon coated. This ensured the best possible imaging conditions with the aim of providing structural characteristics of fish scales. Imaging was performed using FEI Quanta 650 apparatus installed at the Department of Earth Sciences, University of Geneva, Switzerland. The analytical setting included the SEM operation in a high vacuum at an acceleration voltage of 5 kV and a beam power that did not exceed 1 nA. Working distances were as low as 4 mm, thus providing optimal imaging resolution. All imaging was performed using a solid-state concentric backscatter (CBS) detector.

Scale Shapes
Landmark-based geometric morphometrics (GM). Observed scale shape variation was captured by the first two principal components (PC), which cumulatively described 76.7% of the total variance among landmarks. CVA output based on fish origin partially segregated individuals according to the specified groups (W, WF, F) with no significant difference in shape of wild and wild farm-associated groups as measured by both Procrustes and Mahalanobis distances (Supplementary Table S1). The most pronounced group differences were found in the pairwise comparisons of farm and wild or farm and wild farm-associated shape data. Shape changes related to variation along the first canonical variate (CV1, 96.9% of cumulative variation) showed that the scales of the farmed group were relatively wider and shorter, while scales of the wild and wild farm-associated groups were elongated, with a larger distance between the focus and anterior edge, and a straighter posterior edge as opposed to the convex edge observed in the farmed group ( Figure 3a). The cross-validated discriminant analysis using Procrustes shape data correctly classified 49.7% individuals to their origin (Table 2a). The highest classification score was recorded in the farmed group (67.7%) and the lowest in the wild farm-associated group (38%). When individuals were divided according to origin and sampling location, overall classification success was much lower and accounted for 35.8%. A common misallocation pattern was observed between sites of the same fish origin than between fish origins (Supplementary Figure S1a, Supplementary Table S2a). farmed and wild farm-associated individuals this was based on Ba and Mg, and to a lesser extent B ( Figure 4). The CAP model filled with combined trace elements and scale shape data was slightly less sensitive, having an overall allocation success of 96.8% (93/96). In addition to the misallocation of one individual from the WFS2 group, two farmed individuals from Site 1 were also misallocated to the farmed group from Site 2 and to the wild farm-associated group from Site 1 (Supplementary Table  S3b).   Table 1.
Elliptic Fourier analysis (EFA). In terms of overall scale shape, results of the EFA analysis were in line with the shape pattern observed with GM analysis. Namely, the largest proportion of variation in the Fourier coefficients among seabream origins was recorded in the exposed posterior scale portion (Supplementary Figure S2). Significant differences across samples were detected (F = 28.78, p = 0.001) and the LDA clearly separated the farmed group from wild and wild farm-associated ones (Figure 3b). The overall classification score was higher (66.7%) in comparison to the GM score (49.7%). The degree of classification success varied between origins: farmed group (85.16%), wild (57.29%) and wild farm-associated origin (61.5%) (Table 2b). An overall correct classification rate of 43.1% was recorded at the population level, where farmed populations were clearly segregated from both wild and wild farm-associated ones, whereas a high amount of overlap was observed between the WFS1-WS1-WFS2 populations (Supplementary Figure S1b, Supplementary Table S2b). Table 2. Classification success of 548 individuals of gilthead seabream Sparus aurata to their origin using discriminant function analysis, based on (a) Procrustes shape coordinates, and (b) Fourier shape coefficients. Group description codes are given in Table 1

Trace Element Chemistry
Concentrations of all trace elements in fish scales, except Cu, showed statistically significant variations among fish origins from both locations (Table 3). While concentrations of Fe, K, Sr, Ti and Zn were higher (p < 0.05) in scale samples from wild fish, conversely, B, Mn and Mg were higher in farmed fish (p < 0.05). For wild farm-associated fish, concentrations of all elements fell within the concentration range detected between wild and farmed groups, except for Mn, where the concentration was significantly lower in contrast to the other groups. Multivariate analysis revealed that the trace element composition of gilthead seabream scales differed significantly among fish origins and locations (PERMANOVA test; F(Origin) 2,95 = 506.4, p < 0.01; F(Location) 1,95 = 289.3, p < 0.01) where all pair-wise group comparisons were statistically different (p < 0.001) for both factors. The square root values of the multivariate pseudo-variance component from PERMANOVA, that can also be expressed as the percentage difference between assemblages, revealed much greater variation of elemental concentrations composition at the level of origin (6.9%) than location (3.4%). Discriminant analysis (CAP, Figure 4) for individuals grouped by origin and location was highly significant (p < 0.001), with the CAP model showing a leave-one-out allocation success of 97.9% (94/96). Two misallocated individuals from WFS1 were allocated to the equivalent farm-associated habitat WFS2 (Supplementary  Table S3a). CAP axis 1 separated individuals according to origin, where wild farm-associated groups showed a more similar elemental composition pattern than farmed ones (Figure 4). CAP axis 2 separated individuals according to the sampling location. Vector overlays indicated that seven elements (B, Ba, Mg, Ti, K, Sr, Zn) were responsible for group separations, where the classification of wild individuals was dominantly based on differences in K, Sr, Ti and Zn, while for farmed and wild farm-associated individuals this was based on Ba and Mg, and to a lesser extent B (Figure 4). The CAP model filled with combined trace elements and scale shape data was slightly less sensitive, having an overall allocation success of 96.8% (93/96). In addition to the misallocation of one individual from the WFS2 group, two farmed individuals from Site 1 were also misallocated to the farmed group from Site 2 and to the wild farm-associated group from Site 1 (Supplementary Table S3b).   Vector overlays indicate the correlation of the individual element with primary axes. Group codes are given in Table 1.

Scale Microstructures
The number of primary radii extending from the focus to the scale margin differed significantly by fish origin (Table 4, Kruskal-Wallis test, p < 0.05). Pairwise comparisons revealed that scales from farmed fish had a significantly higher number of radii than in other groups. Namely, radii counts ranged from 9-23 in farmed fish, 8-12 in wild and 7-14 in wild farm-associated fish. In addition, the number of scale-edge circuli significantly differed between groups (Kruskal-Wallis test, p < 0.05), where scales of farmed fish had 18% more circuli in the 0.5 mm restricted scale edge area than in the other two groups (Table 4). A higher density of circuli was accompanied with a narrower inter-circular space in the scales of farmed fish, with an average circulus distance of 15.4 ± 1.6 µm, as opposed to the average distances of 21.2 ± 2.0 µm and 20.3 ± 2.1 µm recorded for wild and wild farm-associated scale circuli, respectively (Figure 5a,b). The tooth-like structures (lepidonts) on the inter-radial circuli were numerous and dorsally compressed and functioned in the attachment of scales to the skin. Lepidontal damage was more pronounced in the anterior region of farmed fish scales than in the wild or wild farm-associated scale types (Figure 5c,d).   Table 1.

Scale Microstructures
The number of primary radii extending from the focus to the scale margin differed significantly by fish origin (Table 4, Kruskal-Wallis test, p < 0.05). Pairwise comparisons revealed that scales from farmed fish had a significantly higher number of radii than in other groups. Namely, radii counts ranged from 9-23 in farmed fish, 8-12 in wild and 7-14 in wild farm-associated fish. In addition, the number of scale-edge circuli significantly differed between groups (Kruskal-Wallis test, p < 0.05), where scales of farmed fish had 18% more circuli in the 0.5 mm restricted scale edge area than in the other two groups (Table 4). A higher density of circuli was accompanied with a narrower intercircular space in the scales of farmed fish, with an average circulus distance of 15.4 ± 1.6 μm, as opposed to the average distances of 21.2 ± 2.0 μm and 20.3 ± 2.1 μm recorded for wild and wild farmassociated scale circuli, respectively (Figure 5a,b). The tooth-like structures (lepidonts) on the interradial circuli were numerous and dorsally compressed and functioned in the attachment of scales to the skin. Lepidontal damage was more pronounced in the anterior region of farmed fish scales than in the wild or wild farm-associated scale types (Figure 5c,d).

Discussion
The ability of scale shape and/or scale trace element composition analysis to resolve fish origin or movements between distinct environments has been in the focus of research for several decades, especially for freshwater and anadromous species and references thereof [25,26]. Phenotypic plasticity and local adaptations are important considerations in delineating the population structure of marine fish. Recent studies have shown that scale shape can be used as a good discriminator between congeneric species [27], sympatric phenotypes [28], populations [29] or geographic variants [30]. Here, we explored morphological and microchemistry patterns of gilthead seabream scales from different population origins, i.e., wild, wild farm-associated and farmed, sampled over a relatively small spatial scale that were previously characterized genetically and morphologically (whole-body shape) [6,8]. Based on the results presented above, multiple lines of inference are reported herein. Clear divergence of farmed fish scale shape from the shapes of wild and wild farm-associated fish was shown, with divergence patterns analogous to those previously documented based on the analyses of genetic and whole-body shape from the same population dataset [6,8]. The typical shape generated by farming condition was characterized with a more elongated lateral axis of the exposed portion of the scale and convexed posterior scale edge, enabling the highest classification origin score (85.2%), with outline analysis displaying a 17.41% more correct classification rate than GM (see Table 2). On the contrary, wild and wild farm-associated fish origins shared a similar shape pattern characterized by a more elongated antero-posterior axis and shorter posterior scale margin. Taking into consideration the continuous scale loss as a common feature amongst sea-cage farmed fish because of handling and high fish densities, we can speculate that the farmed shape represents an adaptive response of the animal triggered by farming conditions. Namely, a wider posterior scale field enables greater coverage of the skin surface and serves more efficiently as a protective cover from physical injuries commonly recorded at farms, while the increased number of radii possibly provide better flexibility and anchoring of the scale at high stocking density [31]. It should be noted that LDA displayed a significantly lower correct classification rate in GM that in outline analysis ( Table 2). The small number of distinguishable homologous anatomical points offered by scale shape features (i.e., seven landmarks) limits the amount of shape information that could be extracted from the given configuration. In contrast, EFA analysis considered the whole outline of the scale, conserving all information related to shape with a higher correct classification rate in relation to GM. From a management perspective, automated outline data collection presents another advantage over GM since it is not impacted by the measurer's level of skill and expertise, which can crucially effect analysis results [32]. Therefore, more research is needed to evaluate and compare morphometric methods of distinguishing fish origin based on scale morphology to develop easy to use, fast and accurate tools for assessing the influence of escapees in fisheries landings. Namely, several ecological implications may arise from the presence of escapees in the wild, i.e., competition for resources, spread of diseases and genetic interactions with wild counterparts [6,12]. Moreover, escaped individuals may reduce the value of the species on the market and/or mislabelling of fish origin may mislead the final consumer [16].
Regarding the wild fish origins from the present study, Žužul et al. [6] observed a slight genetic divergence between wild and wild farm-associated gilthead seabream resulting from the use of different spawning grounds. While scale shape failed to clearly distinguish these two fish origins, whole-body shape gave high classification success [8]. Such disparity among the morphological characters could be explained by the fact that whole-body shape is more impacted by the conditional status of a fish, while scale shape is presumably less sensitive to short-term environmental effects and less dependent on the fish's condition [30,33]. This was seen in wild farm-associated fish where both body-morphological and physiological traits were affected because of the permanent food source and shelter offered around tuna farms [7,8].
Genetic and environmental factors may both contribute to intraspecific variability in scale shape in fish [33]. However, because of the overlapping pattern of scale shapes observed among different locations (Site1, Site2) and origins (W and WF), successful classification of individuals back to their population was not accomplished, supporting the idea that environmental factors, i.e., fish density and altered diet as in the case of farmed fish, play a significant role in the abovementioned scale adaptation. Such a farming footprint on scale shape is regularly accompanied by other morphological deviations, enabling very accurate discrimination of fish of farm or wild origin, i.e., the presence of regenerated scales as a result of scale loss, different circuli spacing and scale texture, and lateral line deformities [15,16,34,35]. High fish densities within cages and handling during farming cycle enables physical collisions and friction among individuals [15]; such conditions lead to the higher levels of scale loss, lepidontal alterations and breakage of the circuli on the scales, as found in the present study. A recent molecular study pointed out that subtle variation in scale shape, including disparity across the body, has a traceable genetic basis and that the same major signalling pathway is involved in both early patterning and later shaping of the same phenotype. The Fgf receptor fgfr1b was identified as a prime target in shaping cichlid scales [36]. Still, further research is required to elucidate the impacts of different environmental conditions on the molecular nature of major signalling pathways involved in scale development.
In contrast to scale shapes, which can be mediated both by genetic and environmental influences, scale chemistry can equally serve as a biogeochemical tag [17,18,26,37] and a pollution bioindicator of the habitat [38], as the hydroxyapatite or apatite-(CaOH) layer of scale can accommodate a wide range of elements directly from seawater [39]. In addition, several authors have noted the temporal stability and/or low rates of post-depositional alterations of scale-derived apatite of anadromous species [18,26,37,40], which allows correct fish classification back to their original environment.
In this study, canonical analysis of principal coordinates of scale chemistries correctly classified 97.9% of sampled gilthead seabream to their origin and sampling location. Even though the dataset was limited, seven elements (B, Ba, Mg, TI, K, Sr, Zn) were useful in separating groups. Farmed origin was characterized with high B and Mg contents, wild origin with high Sr, Fe, K and Zn contents, while wild farm-associated origin with low Mn levels. By comparing trace-element concentrations in otoliths to discriminate the wild and farmed seabream, Arechavala-Lopez et al. [41] used an analogue set of elements to successfully assign fish to their origin (96.6% of individuals correctly classified). Still, only a few elements showed a consistent relationship between scales and otoliths. While an increase in Mg and Ba as well as decrease in Sr concentrations in farmed fish in contrast to wild ones were recorded in both biological structures, the opposite pattern was seen for elements such as K, Fe and Zn. Lower concentrations were recorded in scales and higher in otoliths in farmed fish. Such inconsistent results have been reported elsewhere [41] and are coupled with greater scale element variabilities. This was attributed to the fact that the crystal lattice of apatite-group minerals is remarkably tolerant to structural distortion and chemical substitutions, thus leading to extremely diverse compositions [42,43], which is not the case for aragonite [17]. Furthermore, the scale composition is controlled by the bloodstream element supply [44] as opposed to otoliths whose chemistry reflects a more regulated pathway via an endolymph fluid [45].
Recent studies have shown that, for both marine and freshwater systems, the elemental ratio on the surface of fish scales is closely related to aquatic environment in which the fish grew up [38,46]. This is particularly significant for trace-elements such as Sr, Ba and Mn. Several authors have recognized Sr as the most important element in classifying the capture site of fish [18,37,38,46,47]. Wang et al. [38] observed the concentrations of Sr to be approximately two times higher in seawater than in freshwater, parallel with the higher slope of the Sr/Ca ratio in seawater fish scales than in freshwater fish scales. In this study, the abundance of Sr was lower in farmed seabream originating from hatcheries situated along the Atlantic coast of France compared to the Sr content in the scale of the wild origin seabream ( Table 3). The concentration of Sr in seawater does not show significant variations with depth and salinity and is usually within the range of 7.85 ± 0.03 ppm [48]. It is controlled by a range of factors, such as river discharge (i.e., runoff), groundwater runout, ocean crust-seawater interaction, and diagenetic reflux of Sr into the sea following sediment recrystallization [49,50]. The latter is more important in marginal shallow seas, such as the Adriatic Sea, which has been the prime environment of carbonate deposition since the Triassic and whose diagenesis is rather complex, releasing considerable amounts of Sr in the microenvironment during dissolution of carbonate metastable phases [51,52], the best example being the aragonite-calcite transformation [53]. Strontium content of Adriatic seawater therefore reflects its carbonate-rich bedrock geology [54], while the bedrock geology of open oceans may be devoid of a substantial presence of carbonates as in the case of the French Atlantic coast [55]. This readily explains the higher concentration of Sr in the scales of wild origin (i.e., Adriatic) seabream compared to Sr content in the scales of the farmed seabream originating from French hatcheries.
While the observed differences in the scale element profile of farmed and wild fish could be expected considering that the fish originated from different basins and different rearing conditions, the differences observed among wild and wild farm-associated fish are intriguing. Namely, both wild origins were sampled from the same basin, the Eastern Adriatic Sea, where the element composition in the seawater system is relatively stable and unaffected by weather [38]. It could be argued that the main difference among origins lies in fish diet preferences related to habitat use. As an opportunistic feeder, gilthead seabream adapts its diet to available food in the habitat [56]. Thus, wild fish are more oriented to a diet consisting of bivalves, arthropods and gastropods [57,58], while wild farm-associated seabream utilized tuna waste fish feed composed of small pelagic fish [7]. The strong attraction of marine fish species to farms fed on waste aqua feed has been well explored [59,60]. A previous study highlighted the impacts of diet on otolith composition in the marine species Pomatomus saltatrix [61] and thus a certain dietary effect can be expected on the elemental signature of other calcified structures, such as scales [62]. A shellfish-rich diet is a good source of elements such as K, Zn, Fe and Cu [63] and these elements had concentrations almost twice as high in wild fish than in wild farm-associated origin. Still, coasts are increasingly subjected to human pressures, and anthropic water discharges could also have impacted the elemental composition of wild fish scales. Such an impact has already seen elsewhere via otoliths usage for the identification of bioavailable contaminants [64]. To better understand the element profiles of fish scales, both water and fish feed chemistries should also be included in studies.

Conclusions
In summary, the results of the present study suggest that scale shape alone does not allow for reliable determination of interspecific population membership of gilthead seabream sampled from different habitats. Instead, it is useful as an explanatory tool for the identification of farmed fish because of its sensitivity to strong environmental impacts, such as high fish density or fish manipulation. Given that farmed fish will continue to escape into the wild, scale shape pattern may serve as a discrimination tool in controlling recent escape events in the framework of the national monitoring programme. Escaped fish can have ecological and socioeconomic impacts on the environment and local fisheries, due to competition for resources, spread of disease, and alteration of genetic diversity from hybridisation [12,58,65,66]. Several authors have pointed out that the incidence of escaped fish in fishery catches may lead to potential mislabelling of fish origin, selling a farmed fish as a wild one at much higher prices, constituting consumer fraud [16,67]. Summarizing results from 60 studies using different techniques to discriminate farmed from wild fish, Arechavala-Lopez et al. [68] concluded that the external appearance and morphometry of fish are useful markers for rapid assessments, especially for sea bream. Thus, as the most accessible structure with non-lethal sampling requirements, scale shape analyses can be used in monitoring ecologically important areas, such as Natura 2000 zones and marine protected areas.
In the present study, the application of scale microchemistry as a technique provided good results in identifying fish membership among different habitats. Such an ability to discriminate biologically differentiated populations may be helpful in fish traceability and stock delineation for European fisheries management. However, more research is needed to examine the temporal stability in scale hydroxyapatite of marine species, and in the water chemistry and fish diet effects on scale microchemistry, evaluating the performances of scales as a nonlethal biogeochemical tag.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4441/12/11/3186/s1, Figure S1: Plots of (a) canonical variate analysis (CVs) of overall scale shape variation along the first 2 canonical axes based on landmark-based geometric morphometrics analysis and Procrustes shape coordinates, and (b) linear discriminant analysis (LDA) of scale outlines variation along the first 2 axes based on elliptic Fourier analysis and Fourier shape coefficients of 548 individuals of gilthead seabream Sparus aurata from the eastern Adriatic Sea. Group codes are given in Table 1, Figure S2: Scale outlines comparisons of different gilthead seabream Sparus aurata origins based on Elliptic Fourier analysis, Table S1: Difference in scale shapes among wild (W), farmed (F) and wild farm-associated (WF) gilthead seabream Sparus aurata presented by Procrustes distances (a) and Mahalanobis distances (b), Table S2: Classification success of 548 individuals of gilthead seabream Sparus aurata to their origin and sampling location using discriminant function analysis, based on (a) Procrustes shape coordinates, and (b) Fourier shape coefficients. Group description codes are given in Table 1, Table S3: Classification success of from 96 individuals of gilthead seabream Sparus aurata to their origin and sampling location using (a) trace elemental composition of scales, and (b) combined elemental composition and shape coefficient data of scales. Group description codes are given in Table 1.