Monitoring Skeletal Anomalies in Big-Scale Sand Smelt, Atherina boyeri , as a Potential Complementary Tool for Early Detection of Effects of Anthropic Pressure in Coastal Lagoons

: Mediterranean coastal lagoons are increasingly affected by several threats, all concurrently leading to habitat degradation and loss. Methods based on ﬁsh for the assessment of the ecological status are under implementation for the Water Framework Directive requirements, to assess the overall quality of coastal lagoons. Complementary tools based on the use of single ﬁsh species as biological indicators could be useful as early detection methods of anthropogenic impacts. The analysis of skeletal anomalies in the big-scale sand smelt, Atherina boyeri , from nine Mediterranean coastal lagoons in Italy was carried out. Along with the morphological examination of ﬁsh, the environmental status of the nine lagoons was evaluated using a method based on expert judgement, by selecting and quantifying several environmental descriptors of direct and indirect human pressures acting on lagoon ecosystems. The average individual anomaly load and the frequency of individuals with severe anomalies allow to discriminate big-scale sand smelt samples on the basis of the site and of its quality status. Furthermore, a relationship between skeletal anomalies and the environmental quality of speciﬁc lagoons, driven by the anthropogenic pressures acting on them, was found. These ﬁndings support the potentiality of skeletal anomalies monitoring in big-scale sand smelt as a tool for early detection of anthropogenic impacts in coastal lagoons of the Mediterranean region.


Introduction
The big-scale sand smelt Atherina boyeri is a small, short-lived euryhaline fish species, widely distributed throughout the Mediterranean and Black Sea [1] and in the eastern Atlantic Ocean [2]. It usually spends its life cycle in coastal and brackish environments [3] but also can dwell in freshwater habitats [4] and hypersaline conditions [5]. The big-scale sand smelt is a gregarious species and carnivorous, feeding on small crustaceans, worms, molluscs, and fish larvae [6]. For this reason, it is considered as one of the main links between primary benthic and planktonic consumers and the higher trophic level in the food-webs of the aquatic habitats it inhabits [7]. Because of the morphological variability [3,8,9] and polymorphism [10][11][12][13] of the species, the taxonomy of the genus Atherina is still controversial. Several studies based on morphological [12], biochemical [13], and molecular investigations [14,15] have identified three distinct groups that have been proposed to be elevated to the rank of species [16]: one corresponding to the lagoon type (A. lagunae, [16]) and two of marine type: non-punctuated (A. boyeri, [16]) and punctuated (A. punctata, [16]). However, the terminology is not validated, and many studies, especially ecological investigations where no morphological and genetic analyses are performed, still refer to the A. boyeri complex.
A. boyeri is a typical resident fish in most Mediterranean coastal lagoons, along with other small sized fish such as the pipefish (Sygnatus spp.), the killifish (Aphanius fasciatus), the blenny (Salaria pavo), and many Gobiid species. Resident fish are numerically dominant in lagoon fish assemblages [17,18], and although generally without economic value, they play an important role as key links in the lagoon food web. A. boyeri is the only resident fish exploited by Mediterranean lagoon fisheries, that mostly target marine seasonal migrants, such as many species from Sparidae, Mugilidae, and Moronidae families, and the catadromous species Anguilla anguilla [1,[19][20][21].
In the Mediterranean region, lagoon fisheries have always been important from both the socio-economic and the cultural points of view, and still represent a main provisioning ecosystem service delivered by lagoon ecosystems [22][23][24][25]. These include biodiversity conservation, as well as cultural and recreational services, and regulating and supporting services such as protection by coastal erosion, flood control, and waste and pollution assimilation. Lagoon ecosystem services provision relies on the strategic location in the coastal area as well as on the ecological features of these transitional water bodies [21,26]. Due to their intermediate position between the continental and marine domain, coastal lagoons show significant ecological gradients that make these environments very dynamic, at the temporal and spatial scale, and highly resilient and productive. On the other hand, precisely because of their position between land and sea, and along with the increase of human activities and consumption patterns, lagoons are increasingly affected by several threats [27,28]. These can directly or indirectly lead to changes in the structure and functioning of these ecosystems, as well as to habitat degradation and loss. Therefore, the need to assess the overall quality of coastal lagoons considering anthropogenic impacts affecting their ecological status has become a priority.
Fish have been widely used to assess habitat alterations and describe the characteristic of the aquatic ecosystem [29] based on a number of reasons (e.g., sensitivity to stressors, longevity, complex habitat requirements) [30][31][32][33]. The long-standing tradition of ecological, physiological, and ecotoxicological research on a large number of fish species [34][35][36] has allowed to provide insight into environmental degradation at different levels. Since the issuing of the Water Framework Directive (WFD, 2000/06/EC) [37], the interest has shifted towards a global evaluation of the ecological status of aquatic ecosystems. The WFD implementation now also foresees the use of fish as Biological Quality Element (BQE), and this has driven towards the development of suitable methodologies to include fish in the evaluation of the quality status of aquatic ecosystems [38][39][40]. For transitional waters, and specifically for lagoons, methods based on fishes are still under validation. Potential methods rely on the use of multi-metric indices that take into account the fish assemblage as a whole, and supplementary information (e.g., environmental parameters and human activities) to evaluate anthropogenic pressures acting on the lagoon [40][41][42].
On the other hand, [43] have been recently pointing out that an approach based on the use of a single fish species as a biological indicator, in comparison with the use of fish assemblages as a whole, remains valid. Information from single fish populations addressing health, condition status or contamination might be used as complementary information on habitat environmental conditions. Particularly for transitional water bodies, especially wetlands and coastal lagoons, resident fish are considered the most suitable [43] indicator species of lagoon environmental conditions [36,43,44]. They spend their entire life cycle within the water basin [1], and they can integrate the effects of natural environmental stress of the ecosystem [20] and those of pressures deriving from human activities [18,45], with responses at the different levels of the biological organization [46].
Within the many levels of investigation supporting the use of fish as bioindicators of environmental degradation (see [29] for a review), morphological anomalies and skeletal anomalies in fish have been proposed [47][48][49] as potential easily detectable indicators of a disturbed development possibly related to environmental disturbance. Skeletal anomalies are taken into account in one of the first proposed multi-metric index based on fish fauna to evaluate anthropogenic pressures acting on water bodies, the Index of Biotic Integrity (IBI) [50]. The IBI includes among indicators the presence of DELT anomalies, stated as the ratio of fish specimens with deformities, fin erosion, lesions, and tumors [51]. The development of skeletal anomalies is a process that can occur along the entire life cycle of a fish (larva, juvenile, and adult). Moreover, the skeleton is a system characterized by high plasticity, whose anatomy, histology, mechanical properties, and number of elements depend on environmental and genetic factors [52]. The influence of environmental, genetic, and epigenetic factors on the skeleton is further modulated by the ecology, biology, and physiology of fish. Morphological anomalies are rare in natural populations from undisturbed ecosystems, but many authors have hypothesized a relationship between skeletal anomalies and environmental conditions [47,50,53,54]. As a result, skeletal anomalies in fish have been used as biomarkers of pollution stress [55][56][57].
Skeletal deformities (lordosis, scoliosis, etc.) have been detected in the big-scale sand smelt and related to contamination by pollution [58][59][60] and environmental factors [61]. This fish is described as very sensitive to water quality, responding rapidly to its deterioration [1]. Several studies identify the big-scale sand smelt as key fish species to reveal potential contamination (physiological biomarkers, biochemical and genotoxic biomarkers: [62]) and to evaluate the overall environmental quality (body condition [43]) in transitional water ecosystems.
In the present study, skeletal anomalies in A. boyeri complex from nine Mediterranean coastal lagoons in Italy are examined, and their prevalence, frequencies, typologies, and severity evaluated on a comparative basis. The nine study sites have been chosen on the basis of their typology, and of the environmental and management setting. Their quality status is assessed by a methodology based on expert judgement approach that relies on several environmental descriptors of direct and indirect human pressures acting on lagoon ecosystems. The results of the skeletal anomalies analysis in the nine big-scale sand smelt populations are compared with the results obtained from the quality status assessment of the nine lagoons. The aims of the study were: (i) to analyze and describe the skeletal condition of the A. boyeri complex in different populations (ii) to verify if the analysis of skeletal anomalies in this resident species allow s for detecting the presence of anthropogenic pressures on coastal lagoons, thus contributing to the assessment of the ecological status in Mediterranean transitional waters.

Study Sites
The study sites are nine coastal lagoons along the coasts of Italy (Figure 1 Table 1.

Environmental Assessment
The lagoons' overall environmental quality was assessed using a method based on expert judgement [64], whose background, rationale and methodology are detailed in [65]. The method relies on the preliminary choice, by subjective expert judgement, of several indicators of direct and indirect human pressure acting on the lagoon ecosystems, which are quantified on a range scale (from 0-not present to 5-very high). These concur to the objective scoring of values of three Category Pressure Indexes, CPIs, related to three categories of anthropogenic pressure: "Changes in morphology and hydrology" (CPI-Morpho), "Use of landscape and lagoon resources" (CPI-Use) and "Water quality" (CPI-Qual). CPIs scores concur in turn to the calculation of a Final Pressure Index, FPI. For a detailed description of the choice of pressure indicators for the lagoons under study, their attribution to categories and their quantification, data sources, and formulas for the CPIs and FPI calculation, see Section S1 in Supplementary Material file and [65].

Fish Samplings and Preparation
Fish samplings were carried out between 2014-2017. In all sites, big-scale sand smelts are exploited by commercial fisheries, except two (Fogliano and Caprolace). Therefore, 97-135 individuals were taken from commercial catches or scientific surveys at each site. A total of 1023 fish were painlessly euthanized by lethal anaesthesia with 2-phenoxy ethanol 4% and fixed in buffered formalin 10% (pH 7.2; 0.1 M).
Individual standard length (SL, mm) was measured on digital images using the software ImageJ (https://imagej.nih.gov/ij/download.html). Specimens were whole mount stained for calcified tissues with Alizarin red S (modified from [66]) and stored in 100% glycerol. Individuals were examined for skeletal anomalies detection using a Leica MZ12 Stereo Zoom Microscope.

Analysis of Skeletal Anomalies
Skeletal anomalies were classified using an alphanumeric code following [57,67]. For the complete list of skeletal anomalies and the relative binomial key, see Table 2. The anatomical terminology is according to [68]. Table 2. List of the considered skeletal anomalies types. The letter indicates the body region; the number indicates the anomaly type, and the symbols *, # indicate a variation (or a sub-type) of the anomaly type. Severe anomalies are shown in bold.  The data obtained from the monitoring of the skeletal anomalies were used to build a matrix (RM, i.e., raw data matrix). The RM was transformed into a binary matrix (BM, i.e., presence of each type of skeletal malformation = 1; absence = 0). The RM was used to calculate the frequencies (%) of each type of anomaly on the total number of anomalies. The BM was used to calculate the frequencies (%) of the individuals affected by each type of anomaly in each group. For each fish sample, the following metrics were calculated: (i) frequency (%) of individuals with at least one anomaly; (ii) average anomaly load, as number of observed anomalies/number of individuals with anomalies; iii) number of observed types of anomalies; (iv) frequency (%) of individuals with at least one severe anomaly (axis deformation, cephalic anomalies, vertebral centra fusions and deformations); (v) average severe anomaly load, as number of severe anomalies/number of individuals with severe anomalies; and (vi) percentage of severe anomalies, as the ratio of the number of severe anomalies over the number of total anomalies.

Data Analysis
To interpret and compare results of the lagoons environmental assessment, a Principal Components Analysis (PCA) [69] was run on the dataset relative to scores of all pressure indicators in each site. Fish standard lengths (SL) were compared with the Kruskal-Wallis test followed by Mann-Whitney pairwise test with the Bonferroni correction. Differences in SL among samples were considered significant for a p-value < 0.05.
Correspondence analyses (CA) [70] were carried out on the RM and other matrices laid out with sub-set of data to highlight the possible role of anomalies types in discriminating fish samples by the site of origin. As a final step, in order to explore potential relationships among the categorical variables taken into account, i.e., (a) the skeletal anomalies presence in fish, (b) the anthropogenic pressures in the nine study sites, (c) the site of origin of the samples, a Multiple Correspondence Analysis (MCA) [71] was carried out. Based on anomalies' absolute contribution to the CA ordination model, six anomalies were retained for the MCA, i.e., B4z, B5, C6T, E11L, G11, M8*. The dataset of active variables brought together (a) the binary codes (0/1) for the presence of the six selected anomalies) and (b) the anthropogenic pressure indicators. For the latter, each variable (single pressure indicator) was coded with several columns, relative to each possible score (from 0 to 5) for each of the 13 indicators, and removing variables with no variance (all 0, all 1), thus obtaining a total of 45 columns for the environmental variables. For example, the Water renewal time indicator (EXC) is inflated in 5 variables, i.e., EXC/0 . . . EXC/5 corresponding to 5 possible scores. The final binary matrix consisted of 842 rows, i.e., fish presenting at least one of the six selected anomalies, x 51 columns, 6 relative to the presence/absence of each retained anomaly and 45 to the occurrence or not of pressure indicators score levels. One supplementary variable was added to represent the study sites from which fish were sampled.
All statistical procedures were performed with Past v.4.03 software [72], except the MCA that was carried out in R 3.6.0 using the FactorMine package.

Environmental Assessment
The results of the lagoons environmental assessment are reported in Figure 2, where the scores attributed to the 13 pressure indicators are given in radar plots for each lagoon under study. In Table 3  The PCA results, in which the first two PCs account for 68% of the total variation of the dataset, allow to differentiate the nine lagoons based on the relative influence of the 13 pressure indicators. The ordination of the lagoons in the two-dimensional plot ( Figure 3) shows SAB, GOR and GRA with positive loadings for the first PC (40% of total variation), associated, among the pressure indicators, to Aquaculture as the leading variable, Dissolved Inorganic Nitrogen and Chlorophyll-a. ORB, COM, CAB, LES, CAP, and FOG share negative loadings for the first PC, associated with pressure indicators belonging to the CPI-Morpho, such as Inlet and Freshwater supply. The second PC (accounting for 28% of the overall variability) loadings further discriminate these lagoons, in a first group (ORB, COM and CAB) with positive loadings and associated to pressures such as Fixed Barrier and Water Exchange, and a second group (LES, CAP and FOG) with negative loadings and associated to FW Supply, Banks status and Landscape activities.

Skeletal Anomalies Analysis
The average SL of fish from all samples was 60 (±10.2 SD) mm, ORB samples show the smallest size observed (mean SL 45 mm; range 38-90 mm), and GOR shows the largest size (mean SL 74 mm; range 54-95 mm). Differences in SL among samples were significant for a p-value < 0.05 except for CAB, FOG, SAB, LES, and COM.
All fish from all the lagoons under study presented skeletal anomalies, with 100% of specimens affected by at least one anomaly except for a single individual in CAP ( Table 4). The average anomaly load (average number of anomalies per malformed specimen) ranges between 6 and 19, in the samples from ORB and CAB respectively, with intermediate values in samples from the other sites. No lordosis, kyphosis, deformed jaws and deformed opercular plates were observed. Four sand smelts showed anomalies detectable at an external examination: scoliosis in two fishes (equal to 1.4% of individuals) from FOG, one (0.8%) from CAP and one (1%) from ORB (See Tables S1 and S2 in Section S2 of Supplementary Materials). Other severe anomalies were partially or totally fused and deformed vertebral centra (Figure 4). Severe anomalies affected a percentage ranging from 1% (GRA) up to 13% (COM) of fish in the samples (Table 4). Only three sites (CAP, LES, and GRA) showed less than 6% of individuals affected by severe anomalies.   The number of observed types of anomalies in the different samples ranged between 25 (GRA) and 44 (CAB). The most frequent anomalies were those affecting the hemal arches (viz. C6 and C6T) (Figure 4b) of the hemal vertebrae (GRA: 50%; LES: 48%; GOR: 36%; COM: 33% of total anomalies). Furthermore, the 13-21% of the observed anomalies are those affecting the interdorsals (M8*) (Figure 4c) in fish from CAP, FOG, COM, CAB, and SAB. Anomalies involving the epurals (e.g., G10*) and hypurals (e.g., G19) of the caudal fin (range 27-37% of anomalies), and those affecting abdominal vertebral centra (average 11%) are observed in all samples (see Tables S1 and S2 in Section S2 of Supplementary Materials).

Correspondence Analysis
A series of CAs have been conducted first on RM (1023 fish × 63 variables), and then on matrices that consider subsets of data. The results of these statistical analyses are illustrated in the Supplementary Materials, Section 2 and Figure 5. The CA performed on the RM (1023 fish × 63 variables, i.e., 62 anomalies and the variable ABS, to take into account specimens devoid of anomalies) gave an ordination model in which 20.7% of the variance was explained by the first three axes ( Figure S1 in Section S2 of Supplementary Materials). Given the low variance obtained by this ordination model and the unclear distributions of the variables (i.e., skeletal anomalies) with respect to the sites of origin, a second CA was performed on a subset of data. The CA was performed on a matrix (1023 specimens and 22 variables, i.e., 21 anomalies and the ABS variable) that were retained, excluding the rarest anomalies (frequency < 4% in all the specimens) ( Figure S2 in Section S2 of Supplementary Materials). This CA gave an ordination model in which 43.4% of the total variance was explained by the first three axes. Despite the higher explained variance, the ordination model was difficult interpret, because of the weight of the ABS variable, which forces the sites of origin to aggregate close to the origin of the axes.
A further CA performed on the matrix including all fish (1022 in total, the single specimens without any anomalies being excluded) and 21 anomalies (all the anomalies except the rarest) gave an ordination model in which 30% of the total variance was explained by the first three axes. The obtained ordination model fairly discriminates the sites based on the 21 types of skeletal anomalies represented in the corresponding fish samples. In Figure 5

Multiple Correspondence Analysis (MCA)
In order to highlight potential relationships among type of anomalies, sites and anthropogenic pressures, an MCA was carried out. The MCA returned a cumulative explained inertia for the first two dimensions amounting to 38%. Therefore, a two-dimensions MCA solution was considered the most adequate. The results are shown in Figure 6a-c, separately for each categorical variable. A comparison of plots for the two variables, relative to anomalies types ( Figure 6a) and site origin (Figure 6b) evidence along Dimension1 two groups for each variable: one is for negative values of Dim1, with anomaly C6T and sites LES, GRA, and GOR, and one is for positive values (all the other sites, and the other anomalies). Along Dimension2, E11L, B5, B4z, and M8 locate in the same quadrant of ORB, SAB, and CAB; C6T plot in the negative semi-axis of both Dim1 and Dim2, together with GRA and LES. The site GOR plots in the second quadrant, and the sites COM, CAP, and FOG, along with G11, plot in the fourth quadrant. Figure 6c shows the plot of the variables relative to anthropogenic pressure indicators, whose position can be compared with the plot relative to the "site of origin" (Figure 6b). Indicators (Figure 6c) related to the category of land and resources use (CPI-Use), such as   (Figure 6c), those belonging to the CPI-Qual, in particular, Contaminants, Chlorophyll and Dissolved Inorganic Nitrogen with high scores (CON5, CHL4; DIN5) and, among CPI-Use, Aquaculture (AQU5), along with the site GOR. A third group, for negative values of the two Dimensions, brings together variables relative to pressure indicators of all categories presenting intermediate scores, as well as GRA and LES among the sites (Figure 6b).

Discussion
In this study, the skeletal condition of A. boyeri complex was examined in populations from nine Mediterranean coastal lagoons along the coasts of Italy. The study sites differ in typology, ecological features, and management frameworks. The lagoons' quality status was also evaluated based on several environmental descriptors of direct and indirect human pressures acting on the ecosystems. The results allow for some considerations on different aspects, which are discussed in relation to all the available studies on skeletal anomalies in coastal lagoons wild fish.
A first, unexpected, result is that all A. boyeri from all sites presented at least one skeletal anomaly, with only one specimen (from CAP) free of any anomaly. This investigation on skeletal anomalies in big-scale sand smelt has been carried out following the same approach used in a previous study on wild juveniles of three Mugilid species from three sites, two of which were Goro and Lesina [57]. In that study, 16% of juvenile mullet (on a total of 873 juveniles) had at least one skeletal anomaly in Goro, and 15% (on a total of 910) in Lesina. In both lagoons, the same frequency of juvenile fish with severe anomalies (12%) was found, but a higher anomaly load was observed in mullets from Goro (5.3) with respect to those from Lesina (1.9). In the present investigation, carried out with the same methodology but targeting a resident species at the adult stage, sand smelts from Goro showed a worse general condition as well with regard to skeletal anomalies than fish from Lesina, as shown by the average anomaly load (14 vs. 9), the number of anomalies types (34 vs. 30), the occurrence of severely deformed fish (6% vs. 2%), the frequency of severe anomalies (1% vs. 0%) and severe anomaly load (2 vs. 1). In juvenile mullets, severe anomalies such as those affecting jaw and opercular plate, kyphosis, and lordosis were detected [57], while none of these anomalies was found in adult sand smelt from any of sites considered in the present study.
Differences in skeletal anomalies frequencies (both in total and considering only severe anomalies) and in the skeletal anomaly patterns detected in different species dwelling in coastal lagoons, i.e., three Mullet species [57] and the A. boyeri complex (present study), could be ascribed to the different life stages under examination: mullets were all early juveniles while sand smelt were adults. Morphological anomalies are usually rare in natural fish populations, especially in adult fish, because depending on the severity of the anomaly, larvae and juveniles are progressively eliminated from the population due to natural mortality, predation, and competition [29].
Although 100% of sand smelts in almost all sites object of this study were affected by skeletal anomalies, some of them were very rare (affecting ≤ 4% of the total fish) or extremely common (affecting ≥ 75% of total fish). Only one anomaly (G19) was present in the 75% of the total fish analyzed (considered as background noise), while the other 41 affected very few individuals (≤4% of total fish). However, only four sand smelts showed anomalies detectable at an external examination: scoliosis in two fish (equal to 1.4% of individuals) from FOG, one (0.8%) from CAP and one (1%) from ORB.
Such anomalies can be expected in fish species living in fluctuating environments which vary consistently over space and time [1,3], especially in adult big-scale sand smelt, a species with high morphological variability and plasticity [3,8,73,74], and elevate phenotypic polymorphism [10,11]. On the other hand, some metrics, such as the average anomaly load of anomalies and the frequency of individuals with severe anomalies, allowed us to discriminate samples based on the lagoons of origin.
The comparison with published surveys carried out on sand smelt in other sites evidenced the presence of lordosis and/or scoliosis affecting the 3.6% of 1479 sand smelts from the Nerevta river estuary, in the middle Adriatic [59], a single specimen found in the Homa Lagoon, Turkey [60], and 9.7% of 2175 A. lagunae from the Tunis North Lake [58]. All these sites are described as heavily polluted. The absence of lordosis and kyphosis in big-scale sand smelts from the nine lagoons surveyed in this study is then in line with the results achieved with the environmental assessment, that highlight that in all the study sites the anthropogenic impact is of medium level, and with the fact that in none of the study sites severe specific pollution is reported. Many authors have hypothesized a relationship between skeletal anomalies in fish living in degraded ecosystems and environmental conditions [47,50,[53][54][55][56][57]. Anomalies are taken into account in the IBI [51], which incorporates DELT anomalies (including skeletal deformities), and returns a classification of good environmental status for a water body if less than the 2% of examined fish is affected, degraded if the frequency ranges between 2 and 5%, and heavily degraded if it exceeds the 5%. Following this approach, and considering only the severe anomalies in fish, the lagoons of COM, ORB, CAB, FOG, GOR, and SAB should be considered heavily degraded. The severe anomalies were the above mentioned scoliosis and fusions (total or partial) and deformations of the vertebral centra, that contributed to the scores of some general metrics (i.e., frequency of individual with severe anomalies; frequency of severe anomalies; average severe anomaly load). Therefore, the absence of spinal deformations, and the presence of severe anomalies, give a partially contradictory evaluation. More in-depth analyses, based on larger samples from a larger number of sites, would be needed to settle the issue.
However, based on results of the present study and taking into account the twosummary metrics related to skeletal anomalies in fish samples, i.e., the frequency of individuals with severe anomalies and the average individual anomaly load, the former better evidence a relationship with the lagoon overall quality status as evaluated by the expert judgement approach applied in this study. Based on the type and severity of the anthropogenic pressures acting on the nine lagoons under study, in most of them, a medium level of anthropogenic impact was highlighted, with different scores of the Final Pressure Index (higher in GOR, COM, SAB and lower in others: GRA, ORB, FOG, CAB), and only two (CAP and LES) were less exposed to anthropogenic impacts. Accordingly, lagoons with higher FPI values show higher frequencies of severely affected individuals. On the reverse, the less impacted lagoons show lower frequencies of individuals affected by severe anomalies (LES). An exception is GRA, which also shows a low frequency of individuals with severe anomalies, albeit showing a medium level of anthropogenic impact. The average individual load of anomalies seems less representative, showing lower values only in LES and ORB fish.
Some conclusive remarks can be made on the basis of results attained in the present study. With some distinctions, the results obtained with the two different approaches, one based on the assessment of lagoons quality status using anthropogenic pressure indicators and the other which considers the morphological analysis of the skeletal malformations of a resident fish species, return two similar evaluations for the nine coastal lagoons taken into consideration.
The analysis of skeletal anomalies in A. boyeri complex from nine Mediterranean coastal lagoons has allowed us to describe in detail anomalies and their occurrence in this lagoon-resident fish. The skeletal anomalies monitoring in big-scale sand smelt evidenced some potentiality as a tool for early detection of anthropogenic impacts if some overall metrics and specific anomalies are taken into consideration. Some typical anomalies characterize big-scale sand smelt from individual lagoons, while the frequency of individuals with severe anomalies and the average individual anomaly load allowed us discriminate among fish from different coastal lagoons and other anthropogenic impacts.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073 -4441/13/2/159/s1, The following Supplementary Material contains additional information to MM and Results sections and is organized in two sections: Section S1. Environmental quality assessment; S2. Skeletal anomalies. Section S1. Environmental quality assessment. Table S1. Single pressure indicators that contribute to three categories of anthropogenic pressures. Table S2. Pressure indicators, description and degree levels of change (impact score) proposed for the quantification of the anthropogenic pressures. Section S2. Skeletal anomalies. Table S1. Frequencies (%) of individual affected by each anomaly. Table S2. Frequencies (%) of each type of anomaly on the total number of anomalies. Figure S1. Correspondence analysis (1023 fish x 63 variables). Figure S2. Correspondence analysis (1023 specimens and 22 variables). Figure S3. R output for the analysis of the variable contribution for the matrix 1022x2.

Conflicts of Interest:
The authors declare no conflict of interest.