Inorganic Elements in Mytilus galloprovincialis Shells: Geographic Traceability by Multivariate Analysis of ICP-MS Data

The international seafood trade is based on food safety, quality, sustainability, and traceability. Mussels are bio-accumulative sessile organisms that need regular control to guarantee their safe consumption. However, no well-established and validated methods exist to trace mussel origin, even if several attempts have been made over the years. Recently, an inorganic multi-elemental fingerprint coupled to multivariate statistics has increasingly been applied in food quality control. The mussel shell can be an excellent reservoir of foreign inorganic chemical species, allowing recording long-term environmental changes. The present work investigates the multi-elemental composition of mussel shells, including Al, Cu, Cr, Zn, Mn, Cd, Co, U, Ba, Ni, Pb, Mg, Sr, and Ca, determined by inductively-coupled plasma mass-spectrometry in Mytilus galloprovincialis collected along the Central Adriatic Coast (Marche Region, Italy) at 25 different sampling sites (18 farms and 7 natural banks) located in seven areas. The experimental data, coupled with chemometric approaches (principal components analysis and linear discriminant analysis), were used to create a statistical model able to discriminate samples as a function of their production site. The LDA model is suitable for achieving a correct assignment of >90% of individuals sampled to their respective harvesting locations and for being applied to counteract fraud.


Introduction
Fish products are a valuable source of nutrients for humans. Global consumption of seafood increased steadily from 1961 to 2017, at an annual rate of 3.1%, and aquaculture became a fundamental source of fish production, contributing up to 46% of the total fish market. Among marine animal organisms, bivalve mollusks represent an inexpensive source of protein, characterized by high biological value, and also providing essential minerals, trace metals, and vitamins [1]. Mussels, as a result of their filter-feeding habit and sessile lifestyle, tend to accumulate substances from the surrounding environment; hence, they may be considered as natural "devices" enabling pollution biomonitoring [2].
Only after 10 days were the valves perfectly clean and easily crunchable, yielding a fine homogenous powder.
Microwave sample digestion experiments were performed to obtain a clear mineralized solution. Initially, shell samples were digested using nitric and hydrofluoric acid (according to the EPA Method 3052, [24]), but a solid precipitate appeared in the solution (probably the low soluble CaF 2 salt). Other acidic mineralization mixtures were tested without satisfactory results, especially regarding calibration linearity and repeatability. The best conditions were those obtained following EPA Method 3051A [25].

Chemometrics
Input data are organized in datasets (matrices), whose rows correspond to sample replicates (objects), while the columns correspond to the measured variables. Each column corresponds to a metal concentration (overall data matrix is available as supplementary material S1).
A Pearson correlation matrix for the farmed shells, reported in Table 1, describes the associations between variables. Calcium was highly and directly correlated to Ni (0.79). There were other significant correlations, positive between Co and Cu (0.56) and negative between Mg and Ca (−0.50), and between Mg and Ba (−0.62). Such correlations may indicate the coexistence of Ca and Ni in mussels; perhaps these metals may be taken up using an analogous mechanism by the organism. Mg, instead, seems to compete with Ca as the structural material of the mussel shell. The correlation matrix in Table 2 shows the correlations and anti-correlation between elements in the shells collected from natural banks.
The strongest anti-correlation was between Sr and Ni (−0.80). The correlations and anti-correlations between Zn and U (0.60), Ca and Co (−0.60), Mn and Al (0.57), Ca and Ni (0.56), Co and Mg (0.52), Cd and U (0.51), and Cd and Ba (0.50) were significant but less strong.

Principal Component Analysis
A PCA model was computed for the matrix containing data from farmed mussel shells (153 rows × 14 columns). The first two PCs explain 36.9% of the total variance. Figure 1a shows the scores plot in the plane PC1 vs PC2. Figure 1b shows the loadings in the same space of the first two PCs. In the scores plot, letters identify the sampling sites and colors the sampling areas.  Figure 1b shows the loadings in the same space of the first two PCs. In the scores plot, letters identify the sampling sites and colors the sampling areas.
The influence of the different parameters can be evaluated from both the scores plot and the loadings plot. Except for three objects (three replicates of the same samples) in sampling sites "T" and "M", both characterized by high levels of Sr and Ba, Figure 1a shows a good grouping of the sampling sites. Shells characterized by high levels of Ni and Ca are located to the right side of the graph; samples with high levels of Mg, Cu, Co, Mn, and Zn are located to the left side of the graph. However, it can also be observed that sampling sites close to each other are often placed on the opposite sides of the scores plot, as for "E" and "D", both close to the city of Fano.  The influence of the different parameters can be evaluated from both the scores plot and the loadings plot. Except for three objects (three replicates of the same samples) in sampling sites "T" and "M", both characterized by high levels of Sr and Ba, Figure 1a shows a good grouping of the sampling sites. Shells characterized by high levels of Ni and Ca are located to the right side of the graph; samples with high levels of Mg, Cu, Co, Mn, and Zn are located to the left side of the graph. However, it can also be observed that sampling sites close to each other are often placed on the opposite sides of the scores plot, as for "E" and "D", both close to the city of Fano.
The factors affecting the grouping are unknown and could be related to the different water depths, or the different locations of the production plants with regards to the coast, including the relative position with respect to a river mouth, which may strongly influence the inorganic element water loading.
Another PCA model was created for the dataset containing mussels from natural banks (63 objects × 14 variables).
For the natural bank shells, Figure 2a shows the distribution of the objects and Figure 2b the distribution of the variables in the space PC1 vs. PC2. The first two PCs explain 43.9% of the total variance.
including the relative position with respect to a river mouth, which may strongly ence the inorganic element water loading.
Another PCA model was created for the dataset containing mussels from banks (63 objects × 14 variables).
For the natural bank shells, Figure 2a shows the distribution of the objects and 2b the distribution of the variables in the space PC1 vs. PC2. The first two PCs 43.9% of the total variance.
In this case, samples from the same production site are not always well groupe is particularly true for the samples coming from the Pesaro area (red samples in 3a), where no clear grouping is highlighted. This could probably be the result of m growing on natural cliffs, which makes their shell's inorganic composition more in geneous, within the same sampling site, rather than what happens in controlled fa which leads to greater homogeneity.

Linear Discriminant Analysis
A classification method (LDA) was applied to the data in order to further hi the eventual clusters that emerged from the PCA.

Farms and Natural Banks
The first analysis was applied to all samples (dataset with 216 objects × 14 var and considering the mussel origin (farmed or natural banks) as an a priori catego LDA model showed an encouraging NER of 93.5%, with only three farmed objec classified as natural and 11 natural objects misclassified as farmed (5 from samplin "V" and 6 from "Y"). The LDA model gave a 98.0% correct assignment for farme 82.5% for natural growing mussels. Table 3 shows the confusion matrix, sensitivity ificity, and NER for such an LDA model. In this case, samples from the same production site are not always well grouped. This is particularly true for the samples coming from the Pesaro area (red samples in Figure 3a), where no clear grouping is highlighted. This could probably be the result of mussels growing on natural cliffs, which makes their shell's inorganic composition more inhomogeneous, within the same sampling site, rather than what happens in controlled farming, which leads to greater homogeneity.

Linear Discriminant Analysis
A classification method (LDA) was applied to the data in order to further highlight the eventual clusters that emerged from the PCA.

Farms and Natural Banks
The first analysis was applied to all samples (dataset with 216 objects × 14 variables), and considering the mussel origin (farmed or natural banks) as an a priori category. The LDA model showed an encouraging NER of 93.5%, with only three farmed objects misclassified as natural and 11 natural objects misclassified as farmed (5 from sampling point "V" and 6 from "Y"). The LDA model gave a 98.0% correct assignment for farmed, and 82.5% for natural growing mussels. Table 3 shows the confusion matrix, sensitivity, specificity, and NER for such an LDA model. minant plot in Figure 3a, two samples collected in "B" sampling sites were misclassified as "C", one sample from "F" as "E", one sample from "H" as "R", two samples from "M" as "G", two samples from "P" as "Q", and one sample from "R" as "I". Figure 3b shows the corresponding LDA-loadings plot, from which it can be seen that the most important variables for the discrimination were Co, Sr, Ni, Mn, and Zn, those furthest from the origin.
The values of sensitivity and specificity (Table 4), which were always higher than 66.7% (one misclassified sample may weigh for 10-20% if a class carries 9 or 6 samples), confirmed the good performance of the classification model.   Next, data were divided into "farmed" and "natural" samples, in order to check, by LDA, if a stronger discrimination than that observed by PCA was present between sampling sites.
Farmed Samples A LDA model was calculated for the farmed samples (dataset with 153 objects × 14 variables). The sampling sites were considered the a priori category. The classification method applied to this dataset gave an overall NER of 94.1%: only nine samples out of 153 were misclassified, assigning them to different sampling sites. As shown in the discriminant plot in Figure 3a, two samples collected in "B" sampling sites were misclassified as "C", one sample from "F" as "E", one sample from "H" as "R", two samples from "M" as "G", two samples from "P" as "Q", and one sample from "R" as "I".  Figure 3b shows the corresponding LDA-loadings plot, from which it can be seen that the most important variables for the discrimination were Co, Sr, Ni, Mn, and Zn, those furthest from the origin.
The values of sensitivity and specificity (Table 4), which were always higher than 66.7% (one misclassified sample may weigh for 10-20% if a class carries 9 or 6 samples), confirmed the good performance of the classification model. Samples from Natural Banks A LDA model was created from the dataset of natural bank mussel shells (dataset with 63 objects × 14 variables). Similarly, in this case, the a priori category was given to the sampling sites. Figure 4 and Table 5 describe the quality of the classification model. The NER reached the value of 98.4%: only one sample out of 63 was wrongly assigned to another class (an object belonging to the "V" site was assigned to "Y"). The discriminant plot in Figure 4a confirms the good discrimination between all classes, with only a slight overlap between "K" and "U" classes (but without misclassification of samples between them). The LDAloadings plot in Figure 4b shows that Co, Ni, Mg, and Sr are the most important variables for such a discrimination.

All Samples
A further LDA model was computed on the entire dataset, with both farmed and natural samples (dataset with 216 objects × 14 variables). The aim was to evaluate if the sampling sites could be discriminated, regardless of their production method (farmed or natural banks). Thus, in this case, the a priori categories were also the sampling sites. The LDA NER was again very good: 93.5%. Only 14 objects out of 216 were misclassified. Moreover, only three objects from the natural bank samples were misclassified into farmed shells, and only one object from the farmed shells was misclassified as a natural bank samples. The sensitivity and specificity (Table 6) for all 25 classes were higher than 62.5%. The discriminant plot in Figure 5a shows, besides the good discrimination of class "Z", several overlaps between classes, thus, particularly in this case, it is important to evaluate the numerical results of LDA (Table 6) to assess the model performance. The LDA-loadings plot in Figure 5b shows that the most important variables for discrimination are Co, Ni, and Zn. object belonging to the "V" site was assigned to "Y"). The discriminant plot in Figure 4a confirms the good discrimination between all classes, with only a slight overlap between "K" and "U" classes (but without misclassification of samples between them). The LDAloadings plot in Figure 4b shows that Co, Ni, Mg, and Sr are the most important variables for such a discrimination.   62.5%. The discriminant plot in Figure 5a shows, besides the good discrimination of class "Z", several overlaps between classes, thus, particularly in this case, it is important to evaluate the numerical results of LDA (Table 6) to assess the model performance. The LDA-loadings plot in Figure 5b shows that the most important variables for discrimination are Co, Ni, and Zn.  As described above, the implemented statistical model enabled a good classification performance, with a high percentage of correct assignment and efficient discrimination of farmed samples from wild ones. Mussels from natural banks, indeed, grow in an environment significantly different from farmed ones, thus also strongly influencing the aspect and composition of shells. First, the lifetime scale of natural mussels is usually longer and less standardized than farmed mussels, enabling the shells to grow thicker and more robust. Moreover, mussels from natural banks, being closer to the coast and having a larger nutrient availability, share their growing habitat with several other organisms, which also colonize their shells, unavoidably influencing their composition.
In Table 7 we compare our results with the ones obtained in other studies, based on ICP and dealing with the geographic traceability of bivalves by analyzing the mineral composition of their shells. In the table are reported: the species analyzed, the sampling sites and areas, the elements detected, the statistical methods applied, and the models correct assignment score (%). It emerged that the model developed in the present work had a good assignment performance (LDA: 98.4% for NB mussels, 94.1% for farmed mussels), and better than most of the other studies reported in Table 7. Comparable, and even slightly better, were the results obtained by M. Bennion et al. (2019) [18], even if the latter, to achieve a 100% correct assignment, needed to combine statistical modeling with the chemical composition of three structures (clean shell, foot, and periostracum) instead of analyzing only the clean shell, as we did.

Samples and Datasets
Mussels samples belonging to the Mytilus galloprovincialis species were harvested from 25 sites (18 farms and 7 natural banks) along the Marche region coast (central Italy), in the Adriatic Sea, in the framework of a regional surveillance plan. Natural banks are cliffs where mussels grow naturally, while farms are off-shore breeding sites, located roughly two nautical miles away from the coast. Figure 6 shows the geographical distribution of the sampling sites, and Table 8 summarizes the sampling periods and the number of samples for each site. European regulations request periodical monitoring of mollusk production areas, analyzing marine biotoxins, and microbiological and chemical contaminants [4][5][6]. The samples described here were selected from those of commercial size (>4 cm). Shell samples were taken between June and October 2018, and each site was generally sampled 3 times, for a total of 72 samples. Each one of the 72 samples was analyzed in triplicate, and 14 elements (Al, Cu, Cr, Zn, Mn, Cd, Co, U, Ba, Ni, Pb, Mg, Sr, Ca) were determined by ICP-MS, selecting those recognized from the bibliography as the most interesting to characterize geographical origin [22]. The results (72 samples × 3 replicates) were used to create a dataset whose columns were the concentration of the 14 elements (variables) and whose rows were the analysis replicates; therefore, the dimensions of the datasets were: 216 × 14 for the whole data, 153 × 14 for the farmed dataset (153 rows = 51 samples × 3 replicates, 14 variables), and 63 × 14 for the natural bank dataset (63 rows = 21 samples × 3 replicates, 14 variables).

Sample Pretreatment and Preparation
The mussel samples were mechanically opened, the soft tissues removed, and the shells rinsed with water. To analyze a sufficiently homogeneous and representative quantity of sample, 4 valves were submitted to pre-treatment by immersion in H 2 O 2 30% solution for 10 days [22]. The dried shell was then homogenized to obtain a powder. Microwave sample digestion was accomplished following the Environmental Protection Agency procedure 3051A: 0.25 g of homogenized shells was treated with 9 mL HNO 3 67-69% and 1 mL of HCl 37% in a high-pressure microwave system (Milestone-Ethos1-HPR1000), using a temperature/pressure programming procedure [25].
The thus obtained digested solution was diluted with ultrapure water, and adjusted to a final volume corresponding to 50 ± 1.0 g in polypropylene vessels. Further dilutions were performed according to the elements to be analyzed.

Instrumental Analysis
The analyses of Al, Ba, Cd, Ca, Co, Cr, Cu, Mg, Mn, Ni, Pb, Zn, Sr, and U were performed by ICP-MS, on an ELAN DRC II (Perkin Elmer) equipped with an ASX-520 autosampler (Perkin Elmer). The sample was introduced by a peristaltic pump, a Meinhard quartz concentric nebulizer, and a cyclonic quartz spray chamber. The internal standard solution was automatically delivered through a mixing block. Argon (purity 99.999%) was used as the plasma torch, nebulizer, and auxiliary gas, while ammonia (purity 99.99%) was the reaction gas in DRC mode. The ICP-MS parameters were (Table 9): Table 9. Optimized ICP-MS acquisition method. The calibration curves used for quantification were matrix-matched from 0.1 µg L −1 to 500 µg L −1 , choosing for each element the most appropriate range.

ICP-MS PARAMETERS
Each calibration curve was created on at least 3 points other than the matrix sample. The best fit curve was obtained by interpolating the response factors (counts sec −1 element/counts sec −1 internal standard) and the respective concentrations, and using the least-squares algorithm.
The method performances were assessed by evaluating linearity, precision (repeatability), trueness, and quantification limits (LOQ). The last figure of merit was calculated based on signal-to-noise ratio, using the criterion 6σ [26]. The experimental values for LOQ were: 0.002 mg/kg for Cd, Pb, U; 0.020 mg/kg for Ba, Co, Cr; 0.040 mg/kg for Cu; 0.20 mg/kg for Mn, Ni, Sr; 0.50 mg/kg for Al, and Zn; 20 mg/kg for Ca, Mg.

Quality Assurance/Quality Control
The laboratory background contamination level was monitored by analyzing reagent blank samples; accuracy was assessed by batch certified reference material (CRM Dolt5) and spiked sample analysis.
The shell matrix was shown to be highly inhomogeneous, therefore batch-to-batch repeatability checking was necessary, and all the samples were analyzed in triplicate in order to have robust results.

Chemometrics
Pearson correlation matrix was calculated to evaluate correlation or anti-correlation between variables (metals concentration). A direct correlation takes place when the correlation index is positive, otherwise, the variables are anti-correlated [27]. When the correlation index between two variables is less than |0.3|, the two variables are considered weakly correlated. The correlation is considered strong when the correlation index is >|0.7|. When the index is between |0.3| and |0.7| the correlation is considered significant. This choice was empirically made by the authors, according to an already published method [28].
Principal components analysis is a well-known chemometric technique used for the exploratory analysis of a data structure [29]. It consists of projecting the data into a subspace whose axes are called principal components (PCs). The PCs are obtained from a linear combination of the original variables and contain a decreasing amount of information, when moving from the first to the last PC [30]. The most important PCA outputs are the scores and the loadings plots, which show respectively the objects (rows of the dataset) and the variables (columns of the dataset) in the PC space, and allow studying the relationships between them. These graphical models are used to verify whether the objects form any clusters, which would be an indication of possible discrimination power in the variables, and whether the clusters agree with the eventual assignment of samples to a priori categories. When PCA shows clustering, it makes sense to go on attempting the creation of multivariate classification models, in a mathematical and not only graphical form.
Linear discriminant analysis (LDA) is a classification approach based on the maximization of the ratio of the between-class variance to the within-class variance in the dataset, thereby pursuing the maximal discrimination between known sample classes [31]. It assumes that all variables have a Gaussian distribution and are statistically independent. A first validation of the model is carried out by cross-validation (CV) [32]. While PCA is an unsupervised method that allows defining a representation of the data in an orthogonal space without any hypothesis about the belonging of the objects to a given class, LDA is a supervised method that, starting from a priori known classes, aims to create a classification model based on a certain number of independent variables.
Several parameters can be used for the quality estimation of the classification models [33]. All indexes can be derived from the confusion matrix, which is a square matrix with dimensions G × G, where G is the number of a priori classes. The diagonal elements n gg represent the correctly classified objects, while the off-diagonal elements represent the objects erroneously classified.
The non-error rate (NER) is a measure of the quality of a classification model: it represents the percentage of correctly assigned objects. The NER is defined as follows (Equation (1)): where n is the total number of objects. It is also called accuracy, or classification rate. The sensitivity (Sn g ) describes the model's ability to recognize objects belonging to the g-th class and is defined as (Equation (2)): where n g is the number of objects belonging a priori to the g-th class. The specificity (Sp g ) is the ability of the g-th class to reject the objects of all other classes and is defined as (Equation (3)): in which n g ' is the total number of objects assigned to the g-th class [27]. NER, sensitivity, and specificity can assume values between 0 (no class discrimination) and 1 (perfect class discrimination) [34].
Person's correlation matrices and the PCA were calculated by the software CAT [35], based on the R environment (R Core Team, Vienna, Austria). LDA was performed using the R software.

Conclusions
A method enabling the simultaneous ICP-MS analysis of 14 elements (Al, Cu, Cr, Zn, Mn, Cd, Co, U, Ba, Ni, Pb, Mg, Sr, Ca) in mussel shells after microwave digestion was developed and optimized. The method was applied to mussel shells collected in 25 sampling sites located along the central Adriatic Sea coast, and the obtained results were processed by chemometric tools, in order to discriminate mussels as a function of their production method and site.
Among the here applied statistical methods, LDA models were able to achieve a high percentage (>90%) of correct assignments of the samples to a specific sampling point. Therefore, this model could be applied as a useful tool to ensure mussel traceability, counteract fraud, and guarantee product quality.
Moreover, this work indicates that the contents of Al, Cu, Cr, Zn, Mn, Cd, Co, U, Ba, Ni, Pb, Mg, Sr, and Ca in the shells of Mytilus galloprovincialis have discriminating power concerning the production method and sampling site. In particular, LDA showed a potentially interesting discriminating power for Co, Ni, Zn, and Sr, due to the high LDA loadings shown in all the computed models. However, to obtain a robust multivariate screening method, further work is needed: the dataset will be implemented with a higher number of samples for each site and for each sampling period, whose eventual influence on results has not yet been explored.
Supplementary Materials: The following are available online. Table S1: Sample data matrix.

Data Availability Statement:
The data presented in this study are available as supplementary material.
Acknowledgments: Special thanks go to Palombo Paolo from Istituto Zooprofilattico Sperimentale dell'Umbria e delle Marche "Togo Rosati"-Ancona-Italy for analytical support in experimental work and to Valeria Melai from Istituto Zooprofilattico Sperimentale dell'Abruzzo e del Molise "Giuseppe Caporale"-Teramo-Italy for Ca ICP-OES confirmation analysis.