Blue mussels of the Mytilus edulis species complex from South America: The application of species delimitation models to DNA sequence variation

Smooth-shelled blue mussels, Mytilus spp., have a worldwide antitropical distribution and are ecologically and economically important. Mussels of the Mytilus edulis species complex have been the focus of numerous taxonomic and biogeographical studies, in particular in the Northern hemisphere, but the taxonomic classification of mussels from South America remains unclear. The present study analysed 348 mussels from 20 sites in Argentina, Chile, Uruguay and the Falkland Islands on the Atlantic and Pacific coasts of South America. We sequenced two mitochondrial locus, Cytochrome c Oxidase subunit I (625 bp) and 16S rDNA (443 bp), and one nuclear gene, ribosomal 18S rDNA (1770 bp). Mitochondrial and nuclear loci were analysed separately and in combination using maximum likelihood and Bayesian inference methods to identify the combination of the most informative dataset and model. Species delimitation using five different models (GMYC single, bGMYC, PTP, bPTP and BPP) revealed that the Mytilus edulis complex in South America is represented by three species: native M. chilensis, M. edulis, and introduced Northern Hemisphere M. galloprovincialis. However, all models failed to delimit the putative species Mytilus platensis. In contrast, however, broad spatial scale genetic structure in South America using Geneland software to analyse COI sequence variation revealed a group of native mussels (putatively M. platensis) in central Argentina and the Falkland Islands. We discuss the scope of species delimitation methods and the use of nuclear and mitochondrial genetic data to the recognition of species within the Mytilus edulis complex at regional and global scales.


Introduction
Species delimitation, i.e. the act of identifying boundaries at the species level [1,2], is necessary for systematics, ecology and evolution, and fundamental to the accurate assessment of a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 biodiversity as well as for implementing conservation policies. Often this act can be relatively easy owing to allopatry or prezygotic barriers to reproduction, but in many cases species delimitation is made difficult by the presence of cryptic variation and the limitations of many species concepts to effectively recognise such entities [3,4]. In recent years, there has been an increase in the number of methods for delimiting species [5][6][7][8][9][10], some of which attempt to provide a robust theoretical and repeatable analytical framework for species identification. Species delimitation methods usually examine molecular variation data (usually DNA sequences) and utilise phylogenetic reconstructions to statistically identify a threshold that distinguishes species from populations. The accuracy of the method depends largely on the rate of speciation, the population size and the genetic variation used (number of informative sites or loci). These delimitation approaches remove the subjectivity associated with an individual researcher's view of what constitutes a species, and have been helpful in solving taxonomic problems, particularly those of closely related species (e.g., Lemer et al. [11]), but are not yet widely applied. Marine bivalve mussels of the genus Mytilus are naturally distributed throughout all continents of the world excluding Antarctica [12], and as common members of the intertidal community they play important roles in energy transfer from the pelagic to the benthic realm [13,14]. The smooth-shelled blue mussel Mytilus edulis species complex consists of three widely recognised and closely related species: Mytilus edulis Linné, 1758, M. galloprovincialis Lamarck, 1819 and M. trossulus Gould, 1850 [13]. However, there has been longstanding debate about the existence of other smooth-shelled blue mussel species, dating back over 100 years (refer to [15][16][17] for reviews and to Hilbish et al. [14]; Gérard et al. [18] for more recent molecular interpretations). This is true for all of the Southern hemisphere, but particularly so for the Atlantic and Pacific coasts of South America, despite a large body of work addressing the subject in regions such as Chile (e.g., [18][19][20][21][22][23][24][25][26][27][28] and references therein). To the best of our knowledge, models of species delimitation have not been applied to the Mytilus edulis species complex problem, although they may have the ability to resolve the global taxonomy of this widespread and important group.
In South America, blue mussels occur naturally from Dichato, Chile (36˚32 0 S; 72˚56 0 W) on the Pacific coast, around Cape Horn (55˚58 0 S; 67˚17 0 W), and extend north along the Atlantic coastline to a northern limit at Punta del Este, Uruguay (34˚58 0 S; 54˚57 0 W) [29]. The distribution of these animals also includes the Patagonian and Atlantic Ocean islands (e.g. Tierra del Fuego and the Falkland Islands). These blue mussels are important members of the benthic fauna of the South American coast [30], and are an important resource for aquaculture in the region, in particular, in Chile [31]. The Atlantic coast blue mussel was originally described by d'Orbigny [32], based on morphometric grounds, as M. platensis. Only eight years subsequently, again based on morphometric grounds, the Pacific coast blue mussel was described by Hupé [33] as M. chilensis. Since the advent of genetic markers, a range of marker types (with or without associated morphometric analyses) have been applied to native South American mussels to address the question of their taxonomy. Variously, South American mussels have been described as M. edulis-like [13], as M. edulis platensis [24], as M. platensis [26,34], as a Southern hemisphere lineage of Mytilus galloprovincialis also found in other Pacific Ocean regions such as New Zealand and Australia [22,23,35], as Mytilus edulis chilensis [36,37], and as Mytilus chilensis [28,[38][39][40][41][42][43]. Several authors have noted that different marker types and analyses of different genomes (i.e., mitochondrial versus nuclear DNA) may provide different answers, and that newer generations of marker types, increased genomic coverage, and also wider geographic sampling may be required to definitively answer the question of which species occur where (e.g., McDonald et al. [13], Oyarzún et al. [41], Larraín et al. [43]). In addition, specifically relating to the situation on the Pacific coast of South America, it was noted (e.g., Borsa et al. [24], Oyarzún et al. [27]) that at the time the name M. chilensis had no formal standing, despite its widespread usage. This situation has changed recently, and WoRMS now lists M. chilensis Hupé, 1854 as having "accepted" species status. For native mussels on the Atlantic Ocean coast of South America, WoRMS continues to recognise the valid status of M. platensis d'Orbigny, 1846. This fluidity of taxonomic status highlights the challenge of working with a species complex, the advances in taxonomy that new molecular markers can provide, and the difficulty of achieving a taxonomic outcome that is both biologically relevant and accepted by workers in the field.
The development of species-specific nuclear and mitochondrial DNA RFLP assays to blue mussels from many parts of the world has substantially improved our understanding of the phylogeography and specific status of members of the smooth-shelled blue mussel complex [19,23,39,[44][45][46][47][48][49][50][51]. However, the situation in the Southern hemisphere is still much less clear than that in the Northern hemisphere, and in need of further attention for three reasons. First, delimitation of species allows clarification of the taxonomy in a stable and consistent way. This is of importance because the species is the fundamental unit in biology and is used for activities such biodiversity classification, biosecurity monitoring and breeding programmes [52-55]. Second, clarification of the taxonomy and systematics permits a better understanding of the conservation threats posed by invasive mussels and the consequences of hybridisation and introgression on the genetic integrity of native mussels [27,34,43]. Third, species delimitation may have important implications in other areas, such as food production (e.g. aquaculture), where strict regulations exist to protect consumer rights and for reasons of traceability around which species may be grown, moved within a country or between countries, sold on the local market, or exported (i.e. European Normative, Regulation (CE) N˚104/2000 and N˚2065/ 2001 -[39,43,56-58]. In this context, blue mussels form the basis of mussel aquaculture industries in many countries [59,60] so that accurate species determination is important, both in terms of traceability and marketing [43,61,62].
The aim of the present research was to address the question of the taxonomic status of blue mussels from the Atlantic and Pacific coasts of South America, and thereby to better understand the phylogeography of the taxa identified. For context, such work needs to be carried out with reference to other Mytilus species elsewhere. We use both nuclear and mitochondrial DNA sequence variation to delimit putative native and introduced taxa in the full distributional range of Mytilus spp. in South America-Chile, Argentina, Uruguay and the Falkland Islands. The uncertain taxonomy of the Mytilus spp. system is used to test speciation hypotheses with different species delimitation methods, in an attempt to move beyond qualitative assessments of population-specific or lineage/taxon-specific genetic differences. We view this first application of species delimitation models to the South American blue mussels as a case study which, if successful, may then be applied to the global situation using existing and/or new genetic data sets.

Ethics statement
This study was carried out in accordance with the principles of the Basel Declaration and recommendations of Universidad Austral de Chile committee. The protocol was approved by the biosafety institution of the Universidad Austral de Chile (No 008/17). The samplings were carried out according to the authorisation granted by the Subsecretaria de Pesca y Acuicultura (SUBPESCA, Rex No. 2898/2015). The animals involved in this study were minor invertebrates (Mollusca: Bivalves: Mytilidae).

Mussel collection
Mussels were collected from the intertidal region by hand or from the shallow subtidal zone by SCUBA divers. Samples (8 to 20 mussels per site, mean = 17.4) were collected between 2002 and 2016 at 20 sites in Chile, four in Argentina, one in Uruguay and one from the Falkland Islands (Fig 1; Table 1). In total, 348 mussels were sampled within a size range of 27.3 to 65.9 mm shell length. In addition, reference samples (mantle tissue) from the Northern hemisphere including Canada (Bras d' Or lake), UK (Cornwall) and Spain (Moaña-Pontevedra) were used.

DNA extraction and visualisation
Each mussel was dissected and a small piece of mantle edge tissue was fixed in 95% ethanol and stored at 4˚C. A subsample of~30 mg of mantle edge tissue from each individual was used for total genomic DNA extraction using a DNA kit according to the manufacturer's instructions (Geneaid 1 ). Sizes of amplified fragments were estimated from a 50 bp DNA ladder (Invitrogen TM ) on gels stained with SYBR1 Safe DNA.
Standard PCR amplifications were carried out in 25 μL reaction mixtures containing 2 μL DNA template, 0.2 mM total dNTP, 2mM MgCl 2 , 0.4 mM each primer, 1U of Taq (Invitrogen TM ), the manufacturer's PCR buffer, and sterile distilled water. The PCR conditions involved an initial denaturation at 94˚C for 5 min, followed by 30 cycles of 94˚C for 30 s, annealing at 55˚C for 30s and elongation at 72˚C for 1 min, followed by further elongation at 72˚C for 5 min.
Mytilus F-type mtDNA COI and 16S haplotypes are known to possess a more informative phylogenetic signal than M-type mtDNA haplotypes [18,47,66,67]. Because of this, only female mitochondrial sequences were considered in the present study. Amplicons were purified and sequenced by Macrogen (South Korea). Both sequence directions were determined, using the individual primers from the original reaction. DNA sequence was edited using Geneious1  (Table 2). On the other hand, mytilids sequences used in previous publications were downloaded from GenBank ( Table 2). All nucleotide sequences were aligned using MAFFT v.7 [68] under the iterative method of global pairwise alignment [69], and default settings were chosen for all parameters.

Phylogenetic analyses
Phylogenetic trees were constructed using Maximum Likelihood (ML) and Bayesian inference (BI): (i) with species of the family Mytilidae (S1 Fig) and (ii) with species of the genus Mytilus. Evolutionary models and partitioning strategies were evaluated with PartitionFinder v2.1.1 [79], which identified the best partition using the Bayesian Information Criterion, BIC [80]. A ML tree was inferred using GARLI v2.0 [81] with branch support being estimated by nonparametric bootstrap (BS) (200 replicates). Bayesian analyses were performed using MrBayes v3.2 [82]. Each Markov chain was started from a random tree and run for 5.0x10 7 generations with every 1000 th generation sampled from the chain. Stationarity was checked as suggested in Nylander et al. [83]. All sample points prior to reaching the plateau phase were discarded as "burn-in", and the remaining trees were combined to find the a posteriori probability of phylogeny. Analyses were repeated four times to confirm that they all converged on the same results. Posterior probability (PP) values >0. 90 were taken as statistical support for a clade being present on the true tree [84]. The support values (PP and BS) were located in each node (Figs 2 and 3).

Species delimitation within the Mytilus species complex
We employed five different methods for species delimitation, including four separate single locus (COI) analyses and a combined multilocus (COI + 16S + 18S) approach. We recognise that these sorts of models are only as good as the data that are used to test them, and for this reason we are interested to determine if a single marker or a multilocus approach is most powerful. Firstly, we used the Generalized Mixed Yule Coalescent model (GMYC single), a method specifically developed for only one mitochondrial locus [5], and for when the majority of phylogenetic signal is found in mtDNA. This algorithm estimates the number of ''species" by classifying the branching rates of a phylogram as being the result of interspecific or intraspecific lineage branching patterns (sensu Pons et al. [5]). After removing duplicate sequences (COI) because they may cause problems with downstream GMYC analyses [85], the best-fitting substitution model was chosen with the help of PartitionFinder v2.1.1 [79]. Using unique haplotypes, we built an ultrametric phylogenetic tree (Fig 2) in BEAST v1.8.1 [86]. We ran phylogenetic analysis under a lognormal relaxed clock set to an evolutionary rate of 9.51 x 10 −8 [70] considering a coalescent tree with constant population size, using a random starting tree, and with 1 x 10 8 Markov Chain Monte Carlo (MCMC) generations sampled every 1,000 th generation. We implemented two independent runs and combined results using LogCombiner v1.8.1 [86], burning the first 25% of the samples and then using Tracer v1.5 [87] to check for minimum adequate Effective Sample Size (ESS values > 200) and to visually inspect stationarity and convergence by plotting likelihood values. A consensus was built with TreeAnnotator 1.8.1 [86] using the maximum clade credibility method. This tree was used as input to estimate GMYC single in the package SPLITS (SPecies LImits by Threshold Statistics- [88]) using R v3.0.1 [89].
Secondly, we used a Bayesian version of this model (bGMYC), which addresses the uncertainty in the trees by sampling over a posterior distribution of sampled trees [90]. bGMYC analyses were performed by running the eponymous R package on the 100 trees sampled during the MCMC in BEAST v1.8.1 (we discarded the first 90% trees as 'burn-in'). We ran each tree for 50,000 generations, discarding the first 40,000 generations as burn-in and using thinning intervals of 100 generations (as recommended by the authors). The threshold parameter priors (t1 and t2) were set at 2 and 96, and the starting parameter value was set at 25. In the case of bGMYC analyses, the convergence of the MCMC was assessed by checking the evolution graph of the posterior probability against the number of generations, as advised in the bGMYC tutorial. Thirdly, the Poisson Tree Processes (PTP) model [8] was employed to infer molecular clades based on our inferred molecular phylogeny. The PTP method estimates the mean expected number of substitutions per site between two branching events using the branch length information of a phylogeny and then implements two independent classes of Poisson processes (intra and inter-specific branching) before clustering the phylogenetic tree according to the results (sensu Zhang et al. [8]).
Fourthly, we used bPTP, which is an updated version of the original maximum likelihood PTP (maximum likelihood PTP search result is part of the bPTP results). It adds Bayesian support (BYS) values to delimited species on the input tree. A higher BYS value on a node indicates all descendants from this node are more likely to be from one species. The two analyses   (PTP and bPTP) were conducted on the web server for PTP (available at http://species.h-its. org/ptp/) using the MrBayes topology as advocated for this method [8,91]. Fifthly, we used the program BPP [10] to compare different models of species delimitation and species phylogeny in a Bayesian framework, accounting for incomplete lineage sorting due to ancestral polymorphism and gene tree versus species tree conflicts [6,92]. This program conducts multilocus, coalescent-based analyses requiring a guide tree and specification of two priors involving population size and divergence time. Thus, we used the A10 mode, which delimits species using a species tree estimated in BEAST v2.4.8 [93,94]. The population size parameters (θs) were assigned the gamma prior G (9, 100). The divergence time at the root of the species tree (τ0) was assigned the gamma prior G (8, 1000), whilst the other divergence time parameters were assigned the Dirichlet prior [6] (Eq 2). Each analysis was run at least twice to confirm consistency between runs.  Table 2. https://doi.org/10.1371/journal.pone.0256961.g003 The taxonomic index of congruence (Ctax) between pairs of species delimitation methods was estimated, following Miralles and Vences [95]. To identify the most congruent species delimitation approaches, the mean Ctax value for each method was also estimated (S1 Table).

Analysis of broad spatial scale genetic structure
To test the robustness of the species delimitation approach we employed an independent analytical approach to identify population clusters of blue mussels in South America. In principle, if the DNA-based species delimitation approach is accurate, then when applied to the population genetic structure of mussels this should be reflected by different spatially explicit clusters of mussels in South America, for example on the Atlantic and Pacific coasts. This methodology also has the benefit of being able to identify the occurrence and location of regions of interbreeding between two (or more) different clusters of mussels. In addition, the presence of non-native species, such as Northern hemisphere M. galloprovincialis, may be identified. The spatially explicit Bayesian clustering program Geneland 3.2.4 [96] (an extension of program R 3.1.2. R Development Core Team) was used to investigate spatial genetic structure using COI data (Table 1). We converted variable base sites into bi-allelic (allele-like) data, so that the COI input file was a binary file. We ran ten independent runs, where the parameters for possible populations were K = 1-18, and the number of MCMC iterations was 4,000,000, saving every 100 steps (S2 Fig).
After comparing the results of the analyses, we selected the run with the highest posterior probability and post-processed it for graphical display. A burn-in of 10,000 generations (20%) was trimmed from the posterior in the post-processing. A contour map of the posterior mode of population membership was drawn to visualise genetic substructure within South America.

Phylogenetic relationships within the genus Mytilus
Two alignments were performed. First, we aligned the COI data for a total of 600 sites. Testing for the best fit substitution model resulted in the selection of models GRT+I+G (1 st , best fit), TVM+G (2 nd ) and TIM+G (3 rd ) for COI codon position (Fig 2). Then we aligned the three DNA markers for a total of 2638 sites: 441 were variable and 206 were phylogenetically informative. Two of three markers corresponded to mitochondrial dataset with a total of 1068 nucleotide sites, of which 339 were variable and 186 were phylogenetically informative. The evolutionary models and the partitioning strategy obtained in Partitionfinder were SYM+I (1 st ), F81 (2 nd ), HKY+I (3 rd ) (COI), HKY+G (16S) and K80+I (18S).
Known phylogenetic relationships of species within the genus Mytilus based on COI sequence variation were generally well captured. Mytilus californianus Conrad, 1837 (Fig 2,  group a) and Mytilus coruscus Gould, 1861 (Fig 2, group b) formed a divergent clade from the other species of Mytilus, but with relatively low support (BS = 51%, PP = 0.59). M. trossulus (Fig 2, group c), M. edulis (Fig 2, group i) and M. galloprovincialis (Fig 2, group h) all formed reasonably well-supported clades, with M. edulis as the sister species of M galloprovincialis. Mussels from Chile (Fig 2, group d) formed a well-supported clade, different from the species of the Northern hemisphere. The M. chilensis clade appeared to be monophyletic with respect to native Mytilus spp. from Australia (Fig 2, group e) and New Zealand (Fig 2, group f), although with low bootstrap support and posterior probability support (BS = 60%, PP = 0.88).
The placement of samples within the clades from our different geographic collecting sites was generally consistent with expectations, based on published data. The M. chilensis clade (presumptive native mussel diversity on the Pacific coast of South America) contained native mussels from Chile, but also included mussels from the Atlantic coast of Argentina and Uruguay, as well as from the Indian Ocean location of the Kerguelen Islands. The M. galloprovincialis clade contained mussels from the two Mediterranean Sea sites (Italy and Greece), but also contained mussels from Australia, South Africa and Chile. The M. edulis clade contained mussels from Sweden and the Atlantic coasts of the USA and Canada, but also included mussels from Argentina and the Falkland Islands (Fig 2). A final group, of mixed geographic origin, was also recognised, containing mussels from the Atlantic coast of the USA and Finland (Baltic Sea) (Fig 2; AY130046 and AF242015, see detail in Table 2). This geographically anomalous group was described as clade "g". However, this clade showed low support.

Species delimitation analyses
The most congruent result amongst single-and multi-locus analyses recognised eight monophyletic lineages (including the outgroup) as different species (Fig 2 -PTP and bPTP; mean Ctax = 0.89, see all Ctax values in S1 Table). The Mytilus chilensis and M. trossulus lineages were recovered in phylogenetic analyses and were also supported in the Bayesian tree of concatenated sequences (Figs 2 and 3; BS >90; PP>0.99).
In the bGMYC the distribution of ratios of the Coalescence to Yule rates was above 0, and without negative values (S3 Fig), indicating that the model is a good approximation of the reality of the data. Bayesian GMYC analyses detected Mytilus californianus, M. coruscus, M. trossulus, M. galloprovincialis and M. edulis in the Northern hemisphere. Three species were detected in the monophyletic clade (PP = 1; BS = 98%) of the native mussels of the Southern hemisphere (d = Mytilus chilensis from Chile, Argentina, Uruguay and Kerguelen island; e = Mytilus sp. from mainland Australia and Tasmania; and f = Mytilus sp. from NZ- Fig 2). The results from the PTP, bPTP, GMYC and bGMYC analyses differ in their ability to differentiate between species in Australia versus New Zealand. However, the difference in Ctax values between PTP and bGMYC was marginal (PTP mean Ctax = 0.89 and bGMYC mean Ctax = 0.82, see S1 Table), indicating that this clade would include three species (d, e and f- Fig 2).
Results of the BPP species delimitation analyses are shown in Fig 3. In this multilocus analysis five putative Mytilus species were identified, with a speciation rate = 43.480; coalescent rate = 869.620; null logl = 717.926; max logl = 807.662; and P-value < 0.001. As expected, both M. coruscus and M. californianus were recognised as species. M. trossulus (from Canada and USA) was also recognised as a distinct species. Although there was clear separation of M. edulis (containing native mussels from Canada, Sweden, the UK and Argentina) from M. galloprovincialis (containing native mussels from Spain plus introduced mussels from Chile and South Africa) within the tree, the BPP did not recognise these two groups as separate species. The fifth and final group recognised by the BPP was a mixed Southern hemisphere (putatively M. chilensis) clade, containing native mussels from the Auckland Islands (New Zealand Southern Ocean), Chile and Uruguay.
Analysis of broad spatial scale genetic structure. Based on COI sequence variation, the software program Geneland, which takes into account spatial information, indicated that K = 3 groups was the most likely structure in South America (Fig 4A). The three groups were geographically clustered, with one group being centred around northern Chile (Northern hemisphere M. galloprovincialis-highlighted in yellow Fig 4B), the second group being centred around the coast of Argentina and the Falkland Islands (native southern Atlantic Ocean mussels- Fig 4C), and the distribution of the third group spanning southern Chile, the Straits of Magellan, and the southern part of Argentina (native M. chilensis- Fig 4D). This latter group also included the mussels from Uruguay. The assignment probabilities of individuals to their respective clusters were high, at � 0.90 (Fig 4).

Discussion
For the first time, species delimitation models have been applied to both single locus (COI) and multilocus (COI+16S+18S) datasets to examine in a robust, repeatable and objective manner the taxonomic relationships within the genus Mytilus. Our primary focus has been to examine the diversity of the genus Mytilus from South America, including the Falkland Islands. In this test case study, blue mussels collected from sites around South America have been examined in the context of the global diversity of the genus. In addition, genetic sequence variation (COI) that has been used in previous studies where they addressed the evolutionary history of the Mytilus edulis complex (i.e. [18,67,70,71]) was analysed for the South American mussels in a geospatial framework (Geneland).

Species delimitation within the Mytilus edulis species complex
Our samples cover the distribution of Mytilus in South America, with the exception of Brazil where the introduction of M. galloprovincialis has recently been reported [97].  [13], Borsa et al. [24]). On the other hand, our analyses corroborate the monophyly of native mussels on the Pacific coast of South America (Mytilus chilensis Hupé, 1854). In the Southern hemisphere, M. chilensis, along with two distinct evolutionary lineages from New Zealand (putative M. aoteanus Powell, 1958) and Australia (putative M. planulatus Lamarck, 1891) coexist with Mytilus galloprovincialis (the Mediterranean mussel), which has been introduced to parts of Chile (Bahía de Concepción), Brazil, Australia and New Zealand [18,20,22,35,41,67,[97][98][99][100]. Although the most congruent unilocus analysis (PTP: mean Ctax = 0.89) delimited in one species the blue mussels native to Argentina, Chile, Uruguay, Kerguelen Islands, New Zealand, Australia and Tasmania (see Fig 2 -groups d, e, f), the recovery of species from the General Mixed Yule Coalescent model (mean Ctax = 0.82) separated this clade into three species (d = Chile, part of Argentina, Uruguay and Kerguelen Islands; e = Australia; and f = New Zealand). Many previous studies have reported pronounced genetic differentiation between mussels from Australia/Tasmania and New Zealand (e.g., Gérard et al. [18], Sanjuan et al. [101], Pickett et al. [102]). Therefore, the SDM results support the taxonomy described by Lamy [103] and Powell [104] who separated the Southwest Pacific Ocean mussels into two species: New Zealand animals were classified as Mytilus aoteanus Powell, 1958 (currently synonymous with M. galloprovincialis-WORMs), and those of Australia as Mytilus planulatus Lamarck, 1819.
Overall, the GMYC model can be considered suitable for the data. The bGMYC provides reliable results when the branching rate of the coalescent process is higher than the branching rate under a Yule process. For our data the distribution of the ratio of the coalescence rate to the Yule rate is between one and two, with no negative values (S3 Fig). Nevertheless, we recommend studying mussels native to Australia and New Zealand to test the specific hypothesis within the framework of species delimitation.
Generally, the analysis of broad spatial scale genetic structure based on the COI sequence variation used in the spatially explicit Geneland analysis supported the species delimitation models, and identified three main groups within South America. Invasive Northern hemisphere M. galloprovincialis was observed in northern Chile, in the vicinity of Concepcion, consistent with earlier reports of its occurrence here [20,35,105]. A second group of mussels was identified on the Atlantic Ocean coastline of Argentina and also in the Falkland Islands, corresponding to putative M. platensis (e.g., Zbawicka et al. [34]). However, according to the species delimitation analyses, these mussels correspond to clade i (Mytilus edulis-Fig 2). The third and largest group in terms of area of distribution was recorded for the southern Pacific Ocean coastline of South America, south of Punta Lavapié (37˚20'), and included all of the Straits of Magellan and extended north onto the Atlantic Ocean coast of Argentina. This putative species is consistent with M. chilensis (e.g., Larraín et al. [43]). Interestingly, and perhaps somewhat surprisingly, the Uruguay population of Punta del Este, which is near the northern limit of Mytilus sp. on the Atlantic coast of South America, was also identified as belonging to this third group. We do not know if this population is a genuine member of the group or if perhaps it represents a localised introduction of putative M. chilensis into Uruguay. This clearly warrants further examination, but is beyond the scope of this study. The recent report of introduced Northern hemisphere M. galloprovincialis in southern Brazil [97] highlights how quickly the situation can change as records of invasive species establishment change from year to year.

Advantages and limitations of the species delimitation approach
Species delimitation models have strengths and weaknesses and the outcome (interpretation of species-specific differences) is only as good as the data set being used (single versus multiple locus or marker; mitogenome versus nuclear genome; rapidly evolving versus conserved loci) and also will very much depend on the analytical approach (the model and its assumptions) employed.
The four different single locus (mtDNA) species delimitation tests identified nine (GMYC and bGMYC) or seven (PTP and bPTP) species, whereas the one multilocus (mtDNA + nDNA) test (BPP) identified only five species within the genus Mytilus (excluding the outgroup). GMYC [5] uses as input an ultrametric tree estimated from a single locus. The method models the transition point between cladogenesis and allele coalescence by assuming that the former will occur at a rate far lower than the latter. This results in a shift in the rate of branching of the genealogy that reflects the transition between species-level processes and population-level processes (taken from the review by Cartens et al. [2]). On the other hand, bGMYC takes into account phylogenetic uncertainty gene tree estimates using a Bayesian approach [90]. Both implementations of the GMYC are likely to delimit well-supported clades of haplotypes as independent lineages and as such may be prone to over delimitation (sensu Cartens et al. [2]). This probably explains the erroneous delimitation of clade g (see Fig 2).
The BPP method implements a reversible jump Markov chain Monte Carlo search of parameter space, which includes population divergence and estimated distributions of gene trees from multiple loci [6]. The method uses sequence data, and the user is asked to define the topology of the species tree [10]. Then, the posterior probability of the proposed nodes of the species tree is calculated. Whilst inaccurately specified guide trees can lead to false-positive delimitations, the accuracy of BPP is not dependent on its ability to estimate gene trees (sensu Cartens et al. [2]). Cartens et al. [2] indicated that the validation approaches (such as BPP) are often given more weight by empirical investigations because they explicitly model the process of lineage diversification. However, the results should always be interpreted with caution as the methods are not perfect. In our case, the validation of M. chilensis as the predominant species on the Pacific Ocean coast of South America is robust, whereas the validation of M. platensis on the Atlantic Ocean coast is not robust.

Use of sequence data for the identification of species within the genus Mytilus
Many mitochondrial protein coding genes have been used to study phylogenetic relationships amongst species [106] due to the high rate of substitution occurring in the third codon position. However, the nuclear rDNA (e.g. 18S) is one of the most highly conserved DNA regions and has been used to reconstruct phylogenies from phyla to orders [107,108]. Whilst the joint use of nuclear and mitochondrial genetic data has been very useful for recovering a phylogeny in the Mytilidae [109], the 18S gene has a better phylogenetic signal in analysis between genders (e.g., Distel [65], Owada [110], Liu et al. [111]), whereas mitochondrial genes may have more utility in studies of species complexes (e.g. [112][113][114][115]). Because of this concatenation of DNA sequences from different genes (indeed, different genomes with very different properties including size, evolutionary rate, and mode of inheritance) with different strengths or depths of evolutionary signal, it is likely that the multilocus analysis (BPP) failed to recover the species of the Northern hemisphere (M. edulis and M. galloprovincialis), although they were evident in the unilocus analyses (PTP, bPTP and bGMYC). However, in both analyses of species delimitation the specific status of Mytilus chilensis is confirmed for animals that inhabit southern Chile, part of Argentina, Uruguay and the Kerguelen Islands.

Comparison of species delimitation results with recent SNPs analyses
The species delimitation tests were performed on DNA sequence data, some of which has been used for decades to assess phylogenetic and phylogeographic variation (e.g. Wares and Cunningham [70], Riginos et al. [71]). The most recent analyses of Southern hemisphere Mytilus population genetic variation have employed highly variable nuclear DNA markers such as single nucleotide polymorphisms (SNPs). These markers have produced new insights into the taxonomy and phylogeography of native and introduced mussels from New Zealand [116], Chile [43], Argentina [34], Brazil [97] and Australia [100] and offshore islands [117]. Generally speaking the interpretation of the phylogeography, and therefore indirectly the taxonomy, of blue mussels from these Southern hemisphere regions has tended to be based on analyses such as STRUCTURE [118,119], AWclust [120], CA [121,122] and DAPC plots [123] that identify distinct genetic clusters, and by inference some degree of genetic isolation and therefore putative taxonomic identity. SNPs data sets usually provide the greatest level of detail (definition) of all genetic data sets used to date, but are not designed to work in a phylogenetics setting because their high mutation rates are generally not suited to such an approach. This is why taxonomic problems are usually approached from the phylogenetic and coalescence perspective. For the Southern hemisphere mussels, greatest congruence was observed between the GMYC single analysis of the COI data set and the SNPs-based phylogeographic and taxonomic interpretations [34,43,100,116,117]. Overall, the SNPs analyses indicate that in the Southern hemisphere there are clear genetic clusters, that is, distinct genetic entities, in the different major geographic areas surveyed: Chile (putative M. chilensis), Argentina and Uruguay (putative M. platensis), New Zealand (putative M. aoteanus), Australia (putative M. planulatus), plus individuals of mixed ancestry in offshore regions, such as the Falkland Islands (southern Atlantic Ocean), the Kerguelen Islands (southern Indian Ocean), and the Auckland/Campbell Islands of New Zealand (southern Pacific Ocean). The broad-scale phylogeography of Southern hemisphere blue mussels has recently been reviewed by Gardner et at. [124]: our SDM and Geneland results are consistent with interpretations presented by many authors that are described in this review. The congruence between the SNP and GMYC results is probably due to the fact that the precision of the GMYC model increases when a marker with a small effective population size and a high mutational rate is used [125]. Clearly, an important next step will be to test the species delimitation models on the SNPs data sets [34,43,116,117] to compare the objective results from the models against the subjective, and at times controversial, interpretation of the researchers in question in describing the taxonomy of Southern hemisphere blue mussels.

Conclusions
One of the longstanding debates about the taxonomy of the M. edulis species complex, in particular in the Southern hemisphere, has revolved around some of the subtle differences that exist (for example, between native mussels from the Pacific Ocean coast of Chile versus the Atlantic Ocean coast of Argentina, or between South America and Australasia, or on remote offshore islands such as the Falkland Islands and the Kerguelen Islands), and whether these are minor differences that reflect an affinity to Northern hemisphere species (i.e., sub-species status), or larger evolutionary differences that may reflect distinct species. Ultimately, this interpretation may depend on the definition that is applied of what constitutes a biological species [4,126] or to what has been called individual researcher ' . . . taxonomic preconception.' (Gérard et al. [18], p. 84). Regardless, the importance of reaching a standardised and universally agreed taxonomy is apparent, in particular in terms of food labelling, biosecurity and conservation [27,43,59]. The unusual basis of mtDNA inheritance in blue mussels, the close evolutionary relationships amongst the taxa, and the potential for hybridisation between any pair of co-occurring taxa makes the taxonomic classification and phylogenetic reconstruction for this group extremely challenging, using either morphological or molecular data [13][14][15][16][17][18]127]. It is these same challenges that make this group of mussels ideally suited for testing using species delimitation models, in particular with new markers such as SNPs. Whilst the application of five different species delimitation models in this case study focussed on South American mussels has not identified the full range of putative species identities, the approach nonetheless shows great promise if applied to more informative markers such as SNPs as well as to a larger data set with greater geographic coverage.