Development and Validation of a Multi-Locus PCR-HRM Method for Species Identification in Mytilus Genus with Food Authenticity Purposes

Abstract

DNA-based methods using informative markers such as single nucleotide polymorphism (SNPs) are suitable for reliable species identification (SI) needed to enforce compliance with seafood labelling regulations (EU No.1379/2013). We developed a panel of 10 highly informative SNPs to be genotyped by PCR-High resolution melting (HRM) for SI in the Mytilus genus through in silico and in vitro stages. Its fitness for purpose and concordance were assessed by an internal validation process and by the transference to a second laboratory. The method was applicable to identify M. chilensis, M. edulis, M. galloprovincialis and M. trossulus mussels, fresh, frozen and canned with brine, oil and scallop sauce, but not in preserves containing acetic acid (wine vinegar) and tomato sauce. False-positive and negative rates were zero. Sensitivity, expressed as limit of detection (LOD), ranged between 5 and 8 ng/μL. The method was robust against small variations in DNA quality, annealing time and temperature, primer concentration, reaction volume and HRM kit. Reference materials and 220 samples were tested in an inter-laboratory assay obtaining an “almost perfect agreement” (κ = 0.925, p < 0.001). In conclusion, the method was suitable for the intended use and to be applied in the seafood industry.

1. Introduction

Mussels from Mytilus genus represented 20% of the worldwide mollusks production in 2018 [1]. Mytilus genus is one of the most cultivated and marketed bivalve, widely appreciated as a tasty and nutritious source of protein. Mussel aquaculture is a relevant economic activity for many coastal communities [2,3]. The related commercial species, Mytilus edulis (Linnaeus, 1758), Mytilus galloprovincialis (Lamarck, 1819), Mytilus chilensis (Hupé, 1854) and Mytilus trossulus (Gould, 1850) are taxonomically recognized in the World Register of Marine Species [4] and the Integrated Taxonomic Information System [5] together with other Mytilus spp. Nowadays, the European regulation on the common organization of the markets in fishery and aquaculture products [6] requests the declaration of the commercial designation of the species and its scientific name in the label. Each participant country has its official list of commercial designations and scientific names for fishery and aquaculture products (https://ec.europa.eu/fisheries/cfp/market/consumer-information/names_en, accessed on 21 May 2021). This specific traceability requirement is aimed at confirming the authenticity of the products. Species substitution can result in an inexpensive product being labelled as high-priced. It can also affect food safety via unnoticed consumption of allergens due to undeclared species [7,8]. Seafood mislabeling is well documented throughout history [9,10]; it impacts not only food authenticity [11] but also allows the trade in the markets of endangered species or products from illegal, unreported and unregulated (IUU) fisheries, threatening wildlife [12], hampering conservation and negatively affecting consumers decisions [13]. Nowadays, the breadth and depth of mislabeling are coming into sharper focus, thanks to DNA-based species verification methods [14].

DNA testing is currently reported as the regulatory tool of choice for seafood species identification (SI) [15,16,17]. Several DNA analysis techniques rely on Polymerase Chain Reaction (PCR) [18]; however, the informative potential of the targeted genomic regions should be rigorously tested before its adoption in a standardized system. Recently, PCR-based methods have been widely used in seafood authentication [19] and harvesting location detection [20].

For Mytilus SI in either fresh or processed products, different molecular markers and techniques have been developed [21,22,23,24]. However, the most common DNA marker used for this purpose targets the polyphenolic adhesive protein gene, highly conserved within, but polymorphic among species [22,24,25,26,27,28]. Targeting this region, Jilberto et al. [29] developed a high-resolution melting (HRM) assay that allowed them to identify M. chilensis, M. edulis and M. galloprovincialis and their hybrids. This technique is simple and affordable to most molecular analytical laboratories.

SI based in a single locus (“mono-locus”) is relatively easy to perform, but this approach can produce conflicting results among markers and/or when individuals from hybrid zones are analyzed. Single nucleotide polymorphism (SNP) multi-locus panels have allowed for identifying M. edulis, M. galloprovincialis, M. trossulus [30,31], M. chilensis, M. planulatus and M. platensis [32,33,34]. However, genotyping a high number of SNPs may not be practical in routine analysis; therefore, it is attractive to develop reduced panels selecting highly informative SNPs for SI. The identification of the most informative loci can be performed by different criteria such as F_ST outliers or minor allele frequency (MAF) [35,36].

Before the application of an analytical method in food analysis, its performance should be evaluated through compliance with quality criteria according to international guidelines [37,38,39]. Moreover, in the case of laboratories involved in import and export food testing or law enforcement, compliance with ISO/IEC Standard 17025 requirements is necessary [40]. The fit for purpose and the performance of a qualitative real-time PCR method can be first validated in-house, assessing the specificity, sensitivity, repeatability, reproducibility and practicability. Moreover, the transference of the method to a second laboratory can be performed [38].

To support the confidence of consumers, the food industry, business operators and regulators about seafood authenticity, and to avoid the mislabeling of mussels, it is essential to have available an affordable method whose results are internationally recognized.

In this study, we aim to develop a multi-locus SNPs method based on the PCR-HRM analysis, for the identification of species of the Mytilus genus (M. chilensis, M. edulis, M. galloprovincialis and M. trossulus), and to assess its fitness for purpose by an in-house validation process.

2. Materials and Methods

2.1. Method Development and SNPs Selection

To select highly informative SNPs for Mytilus SI (M. chilensis, M. edulis, M. galloprovincialis and M. trossulus) we analyzed two SNP datasets: (A) 49 SNP discovered by data mining from EST in M. edulis, M. galloprovincialis and M. trossulus, genotyped by Sequenom^® MassARRAY iPLEX [41,42]. These SNPs were genotyped in 311 individuals previously identified with the nuclear PCR-RFLP marker ME 15–16, as 189 M. chilensis, 27 M. edulis, 78 M. galloprovincialis and 17 M. trossulus. The second dataset consists of (B) 74 SNPs discovered by Araneda, 2016 [43] using RADseq in M. chilensis, along with 16 SNPs from the previous dataset (from Zbawicka et al. [41,42]) genotyped by the GT-seq method in M. chilensis, M. edulis and M. galloprovincialis. These SNPs were genotyped in 386 individuals previously identified with the nuclear PCR-RFLP marker ME 15–16, as 263 M. chilensis, 41 M. edulis and 82 M. galloprovincialis.

In each dataset, the SNPs were selected using two criteria: the F_ST outlier and the MAF. The first criterion (F_ST outlier) selects those SNPs showing F_ST values above the upper limit of the confidence interval (95%) of the distribution of F_ST obtained with the FDIST2 method [44] implemented in LOSITAN [45] and configured according to Araneda et al. [43]. The rationale of this selection criterion is based on the fact that loci showing F_ST outliers values are candidates for having undergone a selection and adaptation processes during speciation revealing high levels of genetic differentiation, increasing resolution for assigning individuals to species [45]. The second criterion, called MAF_MAX, selects loci that showed maximum differences in the MAF between the analyzed species, calculated by adding the absolute values of the differences between the MAFs for each pair of species. Those SNPs whose MAF_MAX values are over the 80th percentile in the MAF_MAX distribution were considered the most informative. Using this criterion, we select SNPs that showed a fixed allele in one species and the alternative allele present, either fixed or not, in the other species.

To select the highly informative SNPs for SI, four groups of SNPs were tested in each dataset: the first group contained the loci selected by the F_ST outlier criteria, the second group have the loci selected by MAF_MAX criteria, the third group contains loci selected by the F_ST outlier and the MAF_MAX criteria simultaneously. The fourth group contains SNPs selected either by the F_ST outlier or MAF_MAX criteria. The performance of each SNP group was tested, counting the individuals correctly assigned to species based on their multi-locus genotypes. Each individual was assigned to species, using the Training, Holdout and leave-one-out (THL) procedure (Anderson, 2010) implemented in GeneClass2 software [46], with the Bayesian method of [47]. We considered an individual assigned to a species when the probability of assignment was upper to 95% [48]; otherwise, the individual was considered a hybrid.

2.2. Sample Collection

Raw samples from natural populations were analyzed to test in vitro the performance of the selected SNPs. Individuals of the target species, M. chilensis, M. edulis, M. galloprovincialis and M. trossulus (34 per species, n = 136), were collected in four sites from putatively allopatric populations, avoiding described hybrid zones. The species of these samples was confirmed by PCR-RFLP Me15-16 AciI [49] and/or PCR-HRM PAPM (high resolution melting analysis on the polyphenolic adhesive protein gene of mussels) [29]. An individual was considered for further analysis when the two SI methods gave the same result. In addition, 20 individuals from each of the other mussel species commonly found in mussel products (Aulacomya atra and Choromytilus chorus) and species closely related to the Mytilus genus (Perumytilus purpuratus) were collected from natural populations (

Table 1. Experimental design including type and number of samples analyzed in each step of the develop and validation stages.

Species	n	Country of Origin	Product	Packing Media	Development	Validation
					D	A	P	LOD	R	T
Natural Populations
Mytilus chilensis	34	Chile	Raw	-	X	X	X	X	X	X
Mytilus edulis	34	Canada	Raw	-	X	X	X	X	X	X
Mytilus galloprovincialis	34	Spain	Raw	-	X	X	X	X	X	X
Mytilus trossulus	34	Canada	Raw	-	X	X	X	X	X	X
Aulacomya atra	20	Chile	Raw	-	X	-	X	-	-	X
Choromytilus chorus	20	Chile	Raw	-	X	-	-	-	-	X
Perumytilus purpuratus	20	Chile	Raw	-	X	-	-	-	-	-
Commercial Samples
ND	3	Chile	Canned	Brine	-	X	-	-	-	-
Mytilus edulis	3	Germany	Canned	Brine	-	X	-	-	-	-
ND	3	Spain	Canned	Brine	-	X	-	-	-	-
ND	3	Chile	Canned	Oil	-	X	-	-	-	-
ND	3	Chile	Canned	Oil	-	X	-	-	-	-
ND	3	Chile	Canned	Marinated	-	X	-	-	-	-
Mytilus edulis	3	Germany	Canned	Marinated	-	X	-	-	-	-
ND	3	Spain	Canned	Marinated	-	X	-	-	-	-
ND	3	Chile	Canned	Hot sauce	-	X	-	-	-	-
ND	3	Spain	Canned	Scallop sauce	-	X	-	-	-	-
Mytilus edulis	3	Denmark	Canned	Tomato sauce	-	X	-	-	-	-
ND	3	Chile	Frozen	Pre-Cooked	-	X	-	-	-	-

Abbreviations: A, applicability; D, diagnostic test performance (sensitivity and specificity of the method) and primer specificity; LOD, limit of detection; P, practicability; R, robustness; S, SNP selection; T, transferability test. X, Samples used to perform this step; ND: Species not declared.

; Appendix A—

Table A1. Sampling location of individuals from natural populations, collection site, country and geographic coordinates.

Species	Collection Site	Country	Latitude 1	Longitude 1
Mytilus chilensis	Hualaihue	Chile	−41.971389	−72.442778
Mytilus edulis	Prince Edward Island	Canada	46.100219	−63.144747
Mytilus galloprovincialis	Galicia	Spain	42.266208	−8.725078
Mytilus trossulus	Vancouver	Canada	49.292781	−123.127367
Aulacomya atra	Punta Patache	Chile	−20.808631	−70.209255
Choromytilus chorus	Piedra Azul	Chile	−41.548708	−72.590881
Perumytilus purpuratus	Zapallar	Chile	−32.552314	−71.466133

¹ Decimal degrees (DD).

). SI in these specimens was performed by morphological characters in the whole individual, before tissue sampling.

To determine the applicability of the method, canned and frozen pre-cooked commercial mussel products from Germany, Denmark, Spain and Chile were analyzed. Canned mussels with different packing media, “Brine”, “Oil”, “Marinated” and others with sauces (tomato, scallop and hot sauce), were included (Table 1). Of each can or package, three mussels (n = 3) were analyzed.

2.3. DNA Extraction

DNA was extracted from the mantle edge tissue of each individual by different methods: (i) phenol–chloroform method modified for Mytilus [50], (ii) cetyl trimethyl ammonium bromide (CTAB) modified by Arseneau et al. [51] and (iii) Chelex^® 100 method by Chai et al. [52]. The method (i) was used to extract individuals from natural populations and commercial samples; methods (ii) and (iii) were used to extract only the commercial samples. Methods (i) and (ii) include a RNAse treatment; however, as a rapid method, (iii) Chelex^® 100 does not. The quality and quantity of DNA were checked in 0.7% agarose gel in TBE buffer visualized with GelRed^® Nucleic Acid Gel Stain (Biotium Inc, San Francisco, CA, USA) and with a spectrophotometer NanoDrop 2000 (Thermo Fisher, Waltham, MA, USA) using only samples with absorbance ratio 260/280 nm > 1.8. All samples were diluted to a concentration of 20 ng/μL and stored at −18 °C for further molecular analysis.L and stored at −18 °C for further molecular analysis.

2.4. PCR-HRM Analysis and Primer Design

To genotype the selected informative SNP using PCR-HRM, primers were designed with Primer3 v4.1.0 [53] according to Jilberto et al. [29]. The PCR-HRM reaction was optimized at 96 °C for 2 min, followed by 35 cycles at 96 °C for 15 s and 59 °C for 30 s. The final melt curve analysis was performed from 55 °C to 95 °C, with increments of 0.1 °C per second. The reaction contained 0.1 mM of forward and reverse primers (Appendix A—

Table A2. Primer name, amplicon size, polymorphism and reference.

Primer Names	Amplicon Size (bp)	SNP	Reference
PAPM-F/PAPM-R	116–170	G/T	Jilberto et al. (2017)
L2-F/L2-R	106	A/T	This work
L3-F/L3-R	58	A/C	This work
L4-F/L4-R	51	A/G	This work
L5-F/L5-R	54	T/C	This work
L6-F/L6-R	132	A/G	This work
L7-F/L7-R	79	C/T	This work
L8-F/L8-R	54	G/T	This work
L9-F/L9-R	80	C/T	This work
L10-F/L10-R	57	A/G	This work

Jilberto et al. High-resolution melting analysis for identification of commercially important Mytilus species [29].

) and 20 ng of genomic DNA and 1× of HRM SensiFast™ kit (Meridian Biosciences^®, Cincinnati, OH, USA) in a final volume reaction of 8 μL. All analyses were performed in an ECO™ Real-Time PCR System (Illumina, Inc. San Diego, CA, USA) and a magnetic induction cycler quantitative PCR MIC qPCR Cycler (BioMolecular System, Brisbane Queensland, Australia) using the software provided with each equipment. The genotype of each locus was determined by comparing the shape of the melting curve of the query individual with the curve of a reference individual of each species. To be sure that the amplified fragment contains the target SNP, the amplicons of two individuals of each of the four target species were Sanger sequenced in Macrogen (Seoul, South Korea). For each SNP, the amplicon size and polymorphism position was confirmed.

2.5. Diagnostic Test Performance

We evaluated 216 individuals from natural populations (Table 1), calculating Sensitivity (S) and Specificity (E) as in Equation (1) and Equation (2) [54]:

$$\mathrm{S}=\frac{{\mathrm{N}}^{\underset{¯}{\mathrm{o}}} \mathrm{individuals} \mathrm{with} \mathrm{the} \mathrm{species} \mathrm{correctly} \mathrm{identified}}{{\mathrm{N}}^{\underset{¯}{\mathrm{o}}} \mathrm{individuals} \mathrm{from} \mathrm{target} \mathrm{species} \mathrm{analyzed}}$$

(1)

$$\mathrm{E}=\frac{{\mathrm{N}}^{\underset{¯}{\mathrm{o}}} \mathrm{individuals} \mathrm{correctly} \mathrm{excluded} \mathrm{from} \mathrm{the} \mathrm{species} }{{\mathrm{N}}^{\underset{¯}{\mathrm{o}}} \mathrm{individuals} \mathrm{from} \mathrm{non}-\mathrm{target} \mathrm{species} \mathrm{analyzed}}$$

(2)

2.6. Fitness for Purpose Assessment of the Multi-Locus PCR-HRM Method for SI in Mytilus Genus

The qualitative method parameters evaluated in the in-house validation process were applicability, practicability, specificity of primers, sensitivity (limit of detection—LOD), robustness and method transferability [37,38].

2.6.1. Applicability

The applicability or scope of application of the method [55] assessed the possibility of extracting DNA from frozen and canned commercial mussels that can be successfully analyzed by the SI HRM-PCR method (Table 1). Commercial samples were extracted by the phenol–chloroform, CTAB and Chelex methods, and the DNA quality was evaluated as described in Section 2.2. A DNA concentration ≥ 20 ng/µL was considered sufficient. The DNA that successfully met quality and concentration requirements was tested by the PCR-HRM multi-locus panel.

2.6.2. Practicability

The assessment of the ease of operation and throughput of the method [55] considered the minimal unit of analysis, equipment availability, cost, length of the laboratory analysis, occupational safety and staff training requirements [38].

2.6.3. Specificity of Primers

We evaluate the specificity of each primer pairs of the selected SNPs to ensure they amplify only in the target species. First, an in silico evaluation using Primer-BLAST was performed [56]. Secondly, we tested the specificity of each primer pair in vitro to detect any possible amplification in species close to the Mytilus genus (Perumytilus purpuratus) and other species commonly found in mussel products (Aulacomya atra and Choromytilus chorus). The false-negative rate (FNR) and the false-positive rate (FPR) were calculated as in Equations (3) and (4):

$$\mathrm{FNR}=\frac{ {\mathrm{N}}^{\underset{¯}{\mathrm{o}}} \mathrm{individuals} \mathrm{from} \mathrm{known} \mathrm{target} \mathrm{species} \mathrm{misclassified} \mathrm{as} \mathrm{non}-\mathrm{target} \mathrm{species} }{\mathrm{Total} \mathrm{individuals} \mathrm{of} \mathrm{known} \mathrm{target} \mathrm{species} \mathrm{tested}}$$

(3)

$$\mathrm{FPR}=\frac{ {\mathrm{N}}^{\underset{¯}{\mathrm{o}}} \mathrm{individuals} \mathrm{from} \mathrm{known} \mathrm{non}-\mathrm{target} \mathrm{species} \mathrm{misclassified} \mathrm{as} \mathrm{target} \mathrm{species} }{\mathrm{Total} \mathrm{individuals} \mathrm{of} \mathrm{known} \mathrm{non}- \mathrm{target} \mathrm{species} \mathrm{tested}}$$

(4)

2.6.4. Sensitivity (LOD)

To determine the limit of detection (LOD), the DNA concentration in ng/μL was expressed in Haploid Genome Equivalents (HGE) using the C-value. For M. chilensis, it was estimated as 1.56 pg, (an average among the species from Mytilus genus) and for the non-target species (Aulacomya atra), the C-value was 2.24 pg [57]. LOD correspond to the lowest concentration of the target species that can be detected by the method at least 95% of the time. A preliminary screening covering the range of 1000, 600, 200, 50, 10, 2 and 0.1 HGE was tested in the 10 SNPs. Following Broeders et al. [38] protocol, SNPs PAPM and L2 were analyzed in the DNA concentration range of 25, 20, 15, 12, 8, 5 and 1.25 ng/μL in 6-fold. The PCR-HRM was performed in the same conditions described in Section 2.3. The last dilution, where all six replicates gave a positive and specific amplification [58], together with higher and lower dilutions, were tested in 60 replicates [59] for SNP PAPM and L2. The lowest dilution level, at which all 60 replicates showed a specific positive amplification, was considered as the LOD [38].

2.6.5. Robustness

Robustness was evaluated in triplicate for SNPs PAPM and L2 through the fractional factorial design described in Youden and Steiner (1975) [60]. Seven factors or experimental conditions in the protocol were varied to test how small, deliberate modifications affect the results. We considered the following factors: DNA quality through 260/230 and 260/280 ratios, annealing temperature and time, primers concentration, the total volume of reaction and two different HRM kits. Robustness was evaluated either in DNA extracted by the phenol–chloroform and CTAB methods. The standard value of the technique, along with the lowest and highest values for each experimental condition and its combination per run are shown in

Table 2. Parameters used in eight trials for the robustness study *.

		Measurement Level (High–Low)
Factor	Standard Value	1	2	3	4	5	6	7	8
Annealing temperature (°C)	59	60	60	60	60	58	58	58	58
260/230 ratio	≈2.0–2.2	2.1	2.1	0.7	0.7	2.1	2.1	0.7	0.7
260/280 ratio	≈1.8–2.0	2.0	1.7	2.0	1.7	2.0	1.7	2.0	1.7
Annealing time (s)	30	40	40	20	20	20	20	40	40
Primer concentration (ng/μL)	10	15	5	15	5	5	15	5	15
Reaction volume (μL)	8	9	7	7	9	9	7	7	9
Reagent (HRM Kit)	-	S	Q	Q	S	Q	S	S	Q

HRM (High Resolution Melting) Kit: S (SensiFast™ HRM Kit, Meridian Biociences^®) and Q (qPCRBIO™ HRM Mix, PCRBiosystems). * Based on Youden and Steiner (1975) [60].

.

2.6.6. Transferability Test

The standardized operational procedure (SOP) to perform the SI analysis and 220 blind samples, containing 44 of each of the target species (M. chilensis, M. edulis, M. galloprovincialis, M. trossulus) and 44 individuals of the non-target species (Aulacomya atra and Choromytilus chorus), were transferred to a second independent laboratory with different analysts and thermocyclers. The concordance of the results between both laboratories was analyzed using Cohen’s kappa coefficient (κ), where p_o corresponds to the sum of observed concordances and p_e to the sum of concordances attributable to chance (Equation (5)) [61].

$$\mathsf{{\rm K}}=\frac{p_o-p_e}{1-p_e}$$

(5)

3. Results

3.1. Method Development and SNPs Selection

The group of SNPs that in silico analysis showed the best performance assigning individuals to species in dataset A were groups 2 (MAF_MAX) and 4 (F_ST outliers or MAF_MAX) (

Table 3. Assignment success in datasets A and B by criteria.

Group	Dataset A (n = 311)		Dataset B (n = 386)
	N° of SNP	Correctly Assigned Individuals (%)	N° of SNP	Correctly Assigned Individuals (%)
All SNP	49	309 (99.4)	90	363 (94.0)
(1) FST outliers	11	287 (92.3)	27	363 (94.0)
(2) MAFMAX	17	308 (99.0)	21	358 (92.7)
(3) FST outliers and MAFMAX	6	286 (92.0)	11	354 (91.7)
(4) FST outliers or MAFMAX	22	308 (99.0)	37	359 (93.0)

). Both panels correctly assigned 99.0% of the individuals to species, using only 17 and 22 SNPs, respectively. In dataset B, group 1 (F_ST outliers) showed the best performance with 27 SNPs, which correctly assigned 94.0% of the individuals to species. Consequently, we considered the 22 and 27 SNPs included in groups 4 and 1 in datasets A and B, respectively. Two SNPs were shared by both datasets, resulting in 47 most informative SNPs.

In these selected SNPs, all the sequences were checked in silico for any other polymorphism close to the SNP of interest. One SNP was removed because it can give unspecific curves during the HRM analysis. Another 14 were discarded to avoid SNP difficultly to genotype (A/T or C/G) or because it was not possible to design primers. From the remaining 32 SNPs, 12 we discarded because they gave redundant information (seven and five SNPs separate M. trossulus and M. chilensis, respectively, from the other species). Another three SNPs, having MAF values < 0.18 in one species and zero in the others, were considered not informative. Finally, two SNPs were discarded because no genotype was obtained for one of the species (Figure 1). On the other hand, one SNP A/T not considered by previous criteria was included because it was the only one whose allele frequencies were fixed and informative to separate M. edulis from M. galloprovincialis. After in silico selection, 16 SNPs were considered for the in vitro stage.

3.2. Sample Collection, DNA Extraction and Species Identification

The DNA obtained from raw individuals from natural populations satisfied the quality parameters for further analysis. In these samples, the species identified by mono-locus methods targeting the polyphenolic adhesive protein gene was confirmed by morphological traits. The results of the DNA extraction of twelve commercial samples of processed products are shown in Section 3.5.1 Applicability.

3.3. PCR-HRM Analysis and Primer Design

After the in vitro test, one SNP was removed because it amplified only in Mytilus chilensis. Another five SNPs were also discarded although they showed an observable amplification by qPCR, the obtained melting curves were not species-specific. Finally, the PCR-HRM method was successfully standardized in 10 of the most informative SNPs, named L2 to L10, also being the selected SNP amplified by PAPM primers from Jilberto et al. (2017), obtaining melting curves that allow to clearly distinguish each genotype (Figure 2).

Genotyping was successful for the 10 SNPs in all Mytilus individuals from natural populations. The amplicon sequencing of the 10 SNPs reveals that SNP loci were included in the amplified fragments. The 136 samples were correctly assigned to the expected species. The average assignment probability to each species was 100%.

3.4. Diagnostic Test Performance

All Mytilus individuals were correctly assigned to their species, and therefore all the individuals from other genus were correctly excluded. No false negative and false positive identifications were obtained; consequently, the specificity and sensitivity were maximum (S = 1.0, E = 1.0).

3.5. Fitness for Purpose Assessment of the Multi-Locus PCR-HRM Method for Species Identification in Mytilus Genus

3.5.1. Applicability

Starting with a reduced DNA extraction trial, testing two individuals from only four of the canned commercial samples (brine, oil, marinated and with hot sauce as packing media), we concluded that the Chelex method yielded low-purity DNA, evidenced in a 260/280 ratio of 1.18 ± 0.05 (Appendix A—

Table A3. DNA quality (concentration and purity with SD) of mussels from commercial samples obtained by different extraction methods.

DNA Extraction Method by Matrix	N° Cans or Packages	N° Individuals	Concentration (ng/μL)	Purity 260/280 Ratio
Phenol–Chloroform
Canned–Brine	3	9	596.26 ± 128.88	1.73 ± 0.09
Canned–Oil	2	6	766.27 ± 182.44	1.82 ± 0.01
Canned–Marinated	3	9	1345.68 ± 600.18	1.78 ± 0.15
Canned–Hot sauce	1	3	717.00 ± 457.03	1.88 ± 0.04
Canned–Scallop sauce	1	3	768.40 ± 224.99	1.73 ± 0.08
Canne–Tomato sauce	1	3	359.30 ± 76.83	1.68 ± 0.07
Frozen–Pre-cooked	1	3	846.10 ± 693.30	1.77 ± 0.13
		Average	771.29 ± 337.66	1.77 ± 0.08
CTAB
Canned–Brine	3	9	149.58 ± 69.05	1.85 ± 0.03
Canned–Oil	2	6	86.05 ± 38.61	1.81 ± 0.07
Canned–Marinated	3	9	174.94 ± 196.42	1.72 ± 0.63
Canned–Hot sauce	1	3	164.27 ± 37.13	1.88 ± 0.01
Canned–Scallop sauce	1	3	80.00 ± 71.70	1.70 ± 0.15
Canned–Tomato sauce	1	3	75.57 ± 50.70	1.86 ± 0.14
Frozen–Pre-cooked	1	3	323.50 ± 163.77	1.83 ± 0.16
		Average	150.56 ± 89.63	1.78 ± 0.26
Chelex(R) 100
Canned–Brine	1	2	489.00 ± 262.62	1.14 ± 0.08
Canned–Oil	1	2	971.90 ± 245.51	1.26 ± 0.04
Canned–Marinated	1	2	850.50 ± 48.79	1.25 ± 0.08
Canned–Hot sauce	1	2	275.25 ± 21.43	1.08 ± 0.01
		Average	646.66 ± 144.59	1.18 ± 0.05

). Consequently, we continued testing only the phenol–chloroform and CTAB methods in the commercial products listed in Table 1. The average purity was similar between the phenol–chloroform (1.77 ± 0.08) and the CTAB (1.78 ± 0.26) methods and was within the acceptable range (1.7–2.0). However, the DNA purity from mussels canned with tomato sauce, extracted by the phenol–chloroform method, was below the accepted range (1.68 ± 0.07). Therefore, this method was excluded from subsequent HRM analysis.

The average DNA concentration was 771.29 ± 337.66 ng/μL for the phenol–chloroform method and 150.56 ± 89.63 ng/μL for the CTAB method, higher than the established quality limit (20 ng/μL) in all the tested products (

Table 4. Species identification in commercial samples.

Commercial Samples	Packed in	Declared Species	n	Species Identified
				Mch	Me	Mg	Mt	NA
Canned in brine	Chile	ND	3	3
Canned in brine	Germany	M. edulis	3		3
Canned in brine	Spain	ND	3			3
Canned oil	Chile	ND	3	2				1
Canned oil	Chile	ND	3	3
Canned marinated	Chile	ND	3					3
Canned marinated	Germany	M. edulis	3	1				2
Canned marinated	Spain	ND	3					3
Canned hot sauce	Chile	ND	3	3
Canned scallop sauce	Spain	ND	3			2		1
Canned tomato sauce	Denmark	M. edulis	3					3
Frozen pre-cooked	Chile	ND	3	3

n: Number of individuals analyzed. ND: Species not declared. NA: Not assigned species. Mch: M. chilensis. Me: M. edulis. Mg: M. galloprovincialis. Mt: M. trossulus.

).

The CTAB method gave 100% of specific HRM amplifications in canned mussels with brine and hot sauce, and in frozen pre-cooked commercial products. Consequently, it is only applicable to mussels packed in such matrices. Canned mussels in oil and scallop sauce gave 66.7% of successful HRM genotyping, but in mussels marinated and in tomato sauce, the genotyping success was lower than 50%.

An example of the results of SI in the commercial samples used for the applicability step is presented in Table 4. This table only shows samples with eight or more SNP successfully genotyped, allowing SI. Species were identified in frozen products and those packed in brine, oil, hot and scallop sauces. In products packed marinated or in tomatoes sauce, the lower number of SNPs genotyped hampered SI. Although canned marinated products are not within the scope of the method, it is interesting to note that the only species that could be identified in this matrix does not match the one declared in the label (Table 4).

3.5.2. Practicability

Excluding the extraction step, the analysis time for a batch of seven samples (minimum unit of analysis) does not exceed 6 hours, considering one analyst. Real-time PCR equipment and HRM kits are usually available in molecular analysis laboratories. The cost of the necessary reagents to perform analysis per sample was estimated at USD 23. The reagents were not considered dangerous for the handler. The training time of an analyst is four hours. Consequently, the method is practical to apply under the described conditions.

3.5.3. Specificity

In the in silico test, five primer pairs (L6–L10) showed no amplification with any sequence included in the NCBI database. Two primer pairs (PAPM and L2) showed amplification in at least one of the targets species and were considered species-specific. The remaining three pairs of primers (L3, L4 and L5) may not be specific because they showed amplification in non-target species.

During the in vitro test, two primer pairs amplify in non-target species: L7 in C. chorus and L9 in A. atra and P. purpuratus. The Cq values (PCR cycle number at which the reaction curve gets over the threshold line) were statistically similar (Appendix A—

Table A4. Cq values (mean ± standard deviation) for SNPs L7 and L9 in target and non-target species.

SNP	Mytilus	A. atra	C. chorus	P. purpuratus
L7	22.51 ± 1.86	-	22.41 ± 1.48	-
L9	21.98 ± 1.55	23.30 ± 0.85	-	25.24 ± 1.34

) to the target Mytilus species according to the non-parametric Gao test [62].

3.5.4. Sensitivity (LOD)

In the preliminary screening, DNA samples of the target species M. chilensis ranged between 25 and 15 ng/μL and the 10 SNPs gave Cq values < 28. At these same concentrations, the non-target species Aulacomya atra gave Cq values > 30. This means that at the recommended DNA concentration (20 ng/μL), the method is sensitive, detecting the target species and excluding the non-target ones (Appendix A—

Table A5. Cq values in LOD preliminary screening in 10 SNPs in Mytilus chilensis and Aulacomya atra.

HGE	PAPM	L2	L3	L4	L5	L6	L7	L8	L9	L10
Mytilus chilensis
1000	23.117	22.526	23.620	22.943	22.596	23.179	23.101	23.698	23.873	23.217
600	25.539	25.468	26.757	25.462	25.756	25.794	25.788	26.081	25.921	25.821
200	26.940	26.373	27.422	26.758	26.729	27.335	27.129	27.452	26.623	26.838
50	29.306	29.338	30.566	29.374	29.008	29.704	29.605	29.282	29.237	29.184
10	32.327	31.782	34.827	31.972	31.574	32.835	32.149	32.963	31.807	31.385
2	-	34.144	-	33.374	33.213	34.606	33.514	-	34.679	33.464
0.1	-	-	-	-	-	-	-	-	-	-
Aulacomya atra
1000	-	31.124	29.976	26.763	-	-	30.502	34.784	24.858	-
600	-	32.012	28.861	27.868	-	-	31.999	-	26.143	31.696
200	-	-	29.697	-	-	-	34.869	-	30.475	-
50	-	-	-	-	-	-	-	-	31.612	-
10	-	-	-	33.398	-	-	-	-	33.735	-
2	-	-	-	-	-	-	-	-	34.788	-
0.1	-	-	-	-	-	-	-	-	-	-

-: Does not amplify in HRM-PCR.

). Considering the preliminary screening, when PAPM and L2 were genotyped in a narrower concentration range (25–1.25 ng/μL), the lowest dilutions at which all six repetitions gave a positive and specific amplification were 5 and 8 ng/μL, respectively (Appendix A—

Table A6. LOD₆ determination in PAPM and L2.

SNP	Dilution (ng/μL)	Cq ± SD 1	Positive and Specific Amplification (%) 1
PAPM	25	24.29 ± 1.19	94.44
	20	23.66 ± 1.83	77.78
	15	24.14 ± 0.94	94.44
	12	24.98 ± 3.57	77.78
	8	25.39 ± 2.20	83.33
	5	26.87 ± 3.72	77.78
	1.25	20.53 ± 7.99	11.11
L2	25	24.60 ± 0.16	100.00
	20	25.15 ± 0.10	94.44
	15	25.59 ± 0.24	100.00
	12	25.96 ± 0.31	100.00
	8	26.55 ± 0.20	100.00
	5	27.95 ± 0.23	100.00
	1.25	20.45 ± 7.19	50.00

¹ Average considering 3 days, with 6 replicates by tested dilution point per day.

). The LOD was 5 and 8 ng/μL for PAPM and L2, respectively, determined after testing the preliminary screening concentration and the immediately higher (8 ng/μL) and lower (1.25 ng/μL) concentrations (Appendix A—

Table A7. Cq values for the tested dilutions in the preliminary screening to determine LOD.

SNP	Dilution (ng/μL)	Cq ± SD 1	Positive and Specific Amplification (%) 1
PAPM	8	25.34 ± 0.26	100.00
	5 (LOD)	25.72 ± 0.25	100.00
	1.25	27.92 ± 1.28	76.67
L2	8	25.75 ± 0.57	100.00
	5 (LOD)	27.06 ± 0.94	83.33
	1.25	28.95 ± 0.80	6.66

¹ Average considering 60 replicates per dilution point.

). These values are much lower than the recommended DNA working concentration (20 ng/μL DNA), indicating that the method is sensitive in these conditions.

3.5.5. Robustness

The standard deviation of the method (S) does not exceed the acceptance criteria in each of the tests performed (Appendix A—

Table A8. Robustness. Differences (Δ = X-x), between each of the selected parameters by SNP and extraction method.

Parameter	SNP PAPM		SNP L2
	Phenol–Chloroform	CTAB	Phenol–Chloroform	CTAB
Annealing temperature (A-a)	0.105	0.167	0.222	0.505
260/230 ratio (B-b)	5.251	0.382	5.191	0.214
260/280 ratio (C-c)	0.828	0.848	0.311	0.067
Annealing time (D-d)	0.075	0.159	0.300	7.569
Primer concentration (E-e)	0.432	4.807	0.331	4.654
Reaction volume (F-f)	0.942	0.207	0.618	0.662
Reagent (G-g)	0.709	0.683	1.249	0.761
SDi	2.921	0.688	2.890	0.736
Acceptance criteria (S)	3.008	0.971	2.890	0.876

); therefore, the multi-locus PCR-HRM method is considered robust carried out under the conditions previously described.

3.5.6. Transferability Test

Cohen’s kappa coefficient was 0.925, indicating an “almost perfect agreement” in the results of both laboratories (Landis and Koch, 1977). In the 220 samples analyzed, only 15 individuals (6.8%) did not show agreement in the species identification between both laboratories (Appendix A—

Table A9. Observed and theoretical agreement analysis to evaluate the concordance between the results from two laboratories (A and B) applying the PCR-HRM method for mussel species identification.

		Observed Agreements
		Analyst B
		M. chilensis	M. galloprovincialis	M. trossulus	M. edulis	NA	Total Analyst A
Analyst A	M. chilensis	44	-	-	-	-	44
	M. galloprovincialis	-	52	-	-	1	53
	M. trossulus	-	-	43	-	1	44
	M. edulis	-	-	-	24	-	24
	NA	-	10	1	-	44	55
Total Analyst B		44	62	44	24	46	220
Theoretical Agreements Attributable to Chance (Landis and Koch, 1977)
		Analyst B
		M. chilensis	M. galloprovincialis	M. trossulus	M. edulis	NA	Total Analyst A
Analyst A	M. chilensis	8.8	-	-	-	-	44
	M. galloprovincialis	-	14.9	-	-	1	53
	M. trossulus	-	-	8.8	-	1	44
	M. edulis	-	-	-	2.6	-	24
	NA	-	10	1	-	11.5	55
Total Analyst B		44	62	44	24	46	220

NA: Not assigned species.

). The non-matching results arise because, in one laboratory, individuals showed two or more non genotyped loci, resulting in a GeneClass2 score under 95%, not allowing to assign the species. No individuals were assigned to different species by the participant laboratories.

4. Discussion

We developed a multi-locus assay for species identification in the Mytilus genus based on the 10 most informative SNPs from 123 analyzed. Using the multi-locus genotypes, the four target species, Chilean mussel (M. chilensis), Blue mussel (M. edulis), Mediterranean mussel (M. galloprovincialis) and the Bay mussel (M. trossulus), can be confidently identified. This method allowed us to discard the presence of other mussels not belonging to the Mytilus genus but commonly found in seafood products, such as A. atra and C. chorus.

4.1. Multi-Locus PCR-HRM Method Development

The F_ST outlier and MAF_MAX criteria were useful to select the most informative SNPs for SI that are intended to reveal variations among but not within species. The reduced 10 SNP panel was constructed including the SNPs selected by each of these two criteria.

The MAF_MAX criteria maximize the difference in allelic frequencies among species, selecting locus showing a maximum allelic frequency for one species, but a minimum for another (MAX_MAF) [36]. However, to use this approach it is necessary to sample several locations per species and with enough number of individuals, to rule out that the observed allele frequencies are population specific. This criterion is similar to the one used by Wilson et al. [30], these authors choose the more informative SNP using the “loading values” calculated after a discriminant analysis of principal components. Higher loading values reflects SNPs that capture more variance in the allele frequencies and therefore greater differences among species.

A locus shows high F_ST values when two or more groups are homozygous for contrasting alleles. Therefore, the F_ST outlier criterion has been used to select SNPs for geographic assignment because differences in their allele frequencies are probably a product of local adaptation [12,63]. In the Mytilus genus, species are restricted to specific geographic areas, especially the native species from the southern hemisphere, therefore, F_ST outlier SNPs are expected to be informative to differentiate them.

Wilkinson et al. [64] have shown a low correlation between both criteria selecting informative loci, a fact corroborated by our study. Only two loci (PAPM and L7) were selected by MAF_MAX and F_ST outlier at the same time. Therefore, it is important to considerer both criteria together when selecting the most informative loci for SI.

Searching for loci that can be technically genotyped, each of the potentially most informative SNPs must be revised to determine the feasibility to design primers. Moreover, to obtain clearly distinguishable melting curves, amplicons containing more than one polymorphism should be discarded. Other technical aspects of HRM analysis must be also considered [65]. However, optimal conditions scarcely happened, and to design a practical assay it is necessary to be flexible. For example, we included one SNP A/T because was informative to separate between M. edulis and M. galloprovincialis although these kinds of SNPs are not recommended as they are hard to genotype. We also removed some C/T and G/A SNPs because they gave redundant information.

The performance of the reduced 10 SNP panel was similar to the obtained with 49 SNPs (Dataset A) and 90 SNPs (Dataset B). Assignment probability, sensitivity and specificity were maximum for all four species tested here, indicating that the method performs well, correctly identifying the species and excluding individuals that do not correspond to species false positive. To our knowledge, the only reduced multi-locus panel to perform SI in the three European species was developed by Wilson et al. (2018) [30]. The performance of this 12 SNP panel, calculated from their raw data and excluding populations used for validation purposes, showed a maximum specificity for the three species. Sensitivity assigning M. trossulus and M. galloprovincialis was also a maximum; however, in M. edulis, this value reaches only 0.96. The panel developed here can detect one species more with fewer SNPs (10 vs. 12) and with higher performance. As Wilson et al. (2018) [30] have clearly shown, the most valuable aspect of multi-locus panels is the identification of putatively hybrid not detected by the mono-locus assays.

4.2. Fitness for Purpose Assessment of the Multi-Locus PCR-HRM Method

Before the commercial application of any analytical method, it must be validated following accepted quality criteria [37,38]. In laboratories involved in food testing, compliance with ISO/IEC Standard 17025 [40] requirements is mandatory. However, formal validation studies of qualitative techniques that use the PCR-HRM technique to identify marine species are scarce in the current literature. In the case of qualitative PCR-based methods, their scope or applicability relies on the capacity to obtain enough DNA with proper integrity and purity, and also in the absence of PCR interferents. In our study, the CTAB method allowed us to extract enough quality DNA in all tested food matrices. However, every testing laboratory must validate its own DNA extraction method before using this multi-locus assay.

The canning process involves the application of high temperatures (≈121 °C) to the product and its packing media, hermetically sealed in an anaerobic environment [66]. Is widely recognize that in highly processed food (canned or cooked), the DNA is fragmented below 300 bp, leading to extracting less and degraded DNA [9,67]. The HRM analysis use amplicons shorter than 300 bp [11] and the sensitivity of detection is enhanced in smaller amplicons [65]; therefore, it is ideal for genetic analysis of highly processed food. In our method, the size of the amplicons sizes were between 50 to 170 bp, appropriated to be successfully genotyped in processed food. However, the packing media may contain organic acids, ions, chelating agents and other ingredients that favor DNA fragmentation during the thermal process and/or interfere with the PCR reaction [67,68,69,70]. Although we obtained DNA of high quality from products packed in vinegar or sauces containing vinegar, amplification failed or the Cq values were larger than the threshold defined to obtain genotypes. It has been demonstrated that using commercial silica-based columns along with chaotropic solutions allow the removal of some PCR inhibitors in comparison with non-commercial methods such as the phenol–chloroform [71]. Our results pose the challenge to further testing to remove these interferents and/or using amplification enhancers [72,73], to expand the scope of our assay to include complex matrices.

The method was practical and easy to apply, and requires short training time. The HRM technology works in a closed tube, reducing human effort, errors and the risk of drag contamination [74,75]. In addition, it is amendable with small laboratories with low samples flow. The reagents cost are affordable (≈USD 23 per sample in Chile) and could be reduced in countries where it is not necessary to import reagents and consumables.

Although three primer pairs retrieve potential in silico amplifications during Primer Blast analysis, and two of them showed amplification during the in vitro assays in non-target species, the method was specific. The multi-locus approach means that a specimen is assigned to the species using the information of the 10 loci simultaneously. Therefore, the isolated amplification of one or two SNPs in the non-target species will not produce any species assignment. This panel avoid assigned errors derived from the analysis performed in a mono-locus way.

Although this method is qualitative, the DNA concentration that can be detected with 95% confidence (LOD_95% = 5 and 8 ng/μL) was 25–40% of the DNA concentration recommended in the SOP of analysis (20 ng/μL), showing adequate sensitivity for the intended use. The LOD of this qualitative method was appropriated on the base of the Cq values ranges, all positive and specific amplification observed showed Cq values below 29 [38,65]. The method was robust and stable in front of small variations in the eight common operational parameters tested in an orthogonal test as recommended by ISO 20813:2019 [76]. Transferability to a second laboratory allowed the evaluation of the method under other different condition (operator, equipment, etc.). As the method is multi-locus, when one lab fails to amplify some loci, the SI of the sample cannot be obtained. Therefore, no discrepancies in species assignment arise.

5. Conclusions

The multi-locus PCR-HRM method developed in this study was applicable to identify M. chilensis, M. edulis, M. galloprovincialis and M. trossulus mussels, fresh, frozen and canned with brine, oil and scallop sauce. As a limitation, the method is not 100% applicable or efficient in preserves that contain acetic acid (wine vinegar) and tomato sauce. The method is a practical, fast, sensitive and robust diagnostic tool to differentiate the target species of the genus Mytilus from others commonly found in seafood containing mussel.