Genetic divergence of tanaidaceans ( Crustacea : Peracarida ) with low dispersal ability

In this study, the phylogeographic patterns of nuclear, ribosomal and mtDNA gene fragments of five tanaidacean species (Zeuxo, Tanaidae) from the Atlantic, Pacific and Mediterranean Sea were investigated. We aimed to interpret results in the framework of current hypotheses on the distribution of small invertebrates with very limited dispersal ability. Evidence for a surprisingly high genetic divergence was found for intertidal tanaidaceans from the North Atlantic. This is a result of poor dispersal potential, as tanaidaceans have direct development, no pelagic stage, and very limited swimming capacity. However, lower genetic divergence was found between an intertidal tanaid species from the North Atlantic and two from the North Pacific, which suggests a scenario of recent colonization following the last glacial maximum. The species Zeuxo normani was found to be a species complex consisting, at least, of Z. normani (California), Z. cf. normani (Japan), Z. cf. normani (Australia), Z. sp. A (Korea), and Z. holdichi (Spain and France). Our results showed that traditional species identification underestimates tanaidacean diversity and that what have been previously perceived as reliable diagnostic morphological characters, are, however, variable and unreliable.


INTRODUCTION
Small invertebrates have a large potential for analysis of spatial patterns of diversity, providing ideal candidates to test for models in biogeography and macroecology (Witman andRoy 2009, Tuya et al. 2010).However, un-derstanding the diversity of such organisms is hampered by lack of comprehensive data on species distribution and by a difficult, and unreliable, taxonomy (Larsen 2005).A partial solution to this problem is provided by analyses of mtDNA, which have previously proven useful in obtaining reliable data (Larsen and Froufe 2013).
Pleistocene glaciations are known to have played a major role in shaping the phylogeographic patterns of many organisms, particularly during the last glacial maxima (Väinölä 2003, Hewitt 2004, Maggs et al. 2008).Recently, gene flow between the Pacific and the Atlantic has been detected for several marine taxa; it is generally accepted that dispersal often goes in the Pacific-Atlantic direction (Teeter 1973, Meehan et al. 1989, Palumbi and Wilson 1990, Palumbi and Kessing 1991, Taylor and Dodson 1994, Rawson and Hilbish 1995, van Oppen et al. 1995, Väinölä 2003).Gene flow related to inter-glacial periods is suggested for a number of invertebrates, including littoral bivalves (Väinölä 2003), sea urchins (Addison and Hart 2005), crabs (Harrison and Crespi 1999) and even peracarids, such as mysids (Audzijonytė et al. 2005).These findings show that the entire northeast Atlantic region has played an important role in the diversification of the marine biota, especially in the case of cold-water species.Some authors, however, suggest that the situation is a bit more complicated, including dispersal and gene flow going in both directions (Hardy et al. 2011).
There is current consensus that species from temperate regions had to contract their distribution ranges into warmer southern areas during glacial periods, whereas during inter-glacial eras, organisms were able to recolonize warming northern areas (Hewitt 1996).Such changes have resulted in a latitudinal gradient of genetic divergence, with higher levels of diversity at southern locations, which might have functioned as glacial refugia (Luttikhuizen et al. 2008, Kettle et al. 2010, Xavier et al. 2012).For example, in the eastern Atlantic, the Macaronesian archipelagos, the Iberian Peninsula, and the Atlantic coasts of North Africa have been suggested as such glacial refugia areas (Coyer et al. 2004, Maggs et al. 2008, Kettle et al. 2010).
Tanaidaceans are among the most ecologically important and diverse benthic micro-crustaceans in the marine realm.The order currently contains more than 1000 species; the discovery of new species has been prominent in the last decade (Larsen 2005, Blazewicz-Paszkowycz et al. 2012, Bamber and Błazewicz-Paszkowycz 2013).Yet, phylogenetic affinities of tanaidaceans are still poorly understood (Larsen andWilson 2002, Bird andLarsen 2009).Tanaidaceans can occur at high densities, often exceeding 10000 individuals per m2 , with some reports of over 140000 individuals per m 2 (Delille et al. 1985).Despite the obvious ecological importance displayed by such population densities, tanaidaceans are neglected in most ecological surveys, though several ecologists have acknowledged their relevance (Sokolova 1959, Pires 1980, Tuya et al. 2010), for example to monitor environmental impacts (Riera et al. 2011).Species of the genus Zeuxo (Tanaidacea: Peracarida) are broodcaring, non-swimming, and semi-sedentary animals for which long-distance dispersal is thought to be possible only by rafting, human transport via ballast water, or fouling (Larsen 2005).Therefore, their low dispersal potential makes them a good model for studying past vicariance events associated with glacial survival.
The identification of cryptic, or sibling, species is a well-known problem in systematics (Knowlton 1993).Traditional classification has focused mainly on morphological comparisons of organisms, sometimes resulting in conflicting results when the homology of morphological characters was difficult to establish.For example, the morphological differentiation among species of tanaidaceans is notoriously difficult, because of factors such as sexual and ontogenetic polymorphism, as well as the presence of cryptic (and sometimes even sympatric) species (Larsen 2001).This is quite apparent for tanaidaceans, since even family-level diagnostic characters might change with development and gender (Larsen andWilson 1998, Larsen 2001).This paper aims to address these problems of species identification of cryptic species complexes based on morphological characters by combining these with molecular information.
The general purpose of this study was to investigate the phylogeographic patterns of certain species of Zeuxo, and in particular what has been thought to be a single species, Z. normani (Richardson 1905), which has an almost cosmopolitan distribution.This species is recorded from all over the world, including Australia (Edgar 2008), Bermuda (Greve 1974), California (Richardson 1905, Sieg 1980), Japan (Shiino 1951), Sri Lanka (Stebbing 1905) and Spain (Cacabelos et al. 2010).Such a cosmopolitan distribution is highly unlikely for organisms with low dispersal potential, such as tanaids (Larsen 2005), so the presence of a species complex is expected.To compare genetic divergence between species of the genus from different parts of the world, we include analysis of Z. holdichi (Bamber, 1990), Z. exsargasso (Sieg, 1980) and two additional new species from Korea and Turkey.In particular, we compare diversity at the species level by coupling traditional taxonomy and DNA analyses.To accomplish this objective, we examined DNA sequence variation using three different genetic markers, one mitochondrial (COI), one nuclear large subunit ribosomal DNA sequence (28S rRNA), and one nuclear marker (H3).

Taxon sampling
i) Thirteen specimens of Zeuxo (this species is here termed 'A' because we believe it is a new undescribed species and to separate it from sp 'B' from Turkey), collected by scraping from buoys and harbor walls by Mr. Ho Sung Hwang from three different locations in Korea were examined and six (two from each locality) were successfully sequenced (Table 1).
ii) Six specimens of Zeuxo holdichi (Bamber, 1990) were collected from rock pools by Dr. F. Arenas from one location in Northern France.Two were successfully sequenced but only for the 28S gene.Forty-one additional specimens were collected from a variety of algae by snorkelling at 1-1.5 m depth by K. Larsen and E. Froufe from one location in northern Spain.All specimens were examined morphologically and six were successfully sequenced (Table 1).
iii) One hundred and nineteen specimens of Zeuxo exsargasso (Sieg, 1980) were collected by K. Larsen from a variety of algae and eel-grass by snorkelling at 0.5-1.0m depth from two locations in the Cape Verde archipelago and by F. Tuya from a variety of algae by snorkelling at 1.5-4 m depth at one location at Gran Canaria Island.All specimens were morphologically examined and four individuals from Cape Verde and five from Gran Canaria were successfully sequenced (Table 1).iv) Fifty specimens of Zeuxo sp.B collected by Mr. Fikret Öndes from one location in Turkey were examined; three were successfully sequenced (Table 1).v) Over 200 specimens of Tanais dulongii Audouin, 1826 were collected by K. Larsen, F. Tuya, and E. Froufe from a variety of algae by snorkelling at 0.5-1.5 m depth at Mindelo, northern Portugal.Multiple specimens were successfully sequenced.Additionally, three specimens from Helgoland (collection date missing) were successfully sequenced (Table 1).

DNA extraction, amplification and sequencing
To locate and evaluate the taxonomic value of the new samples relative to other Zeuxo species, we included the only sequences available in GenBank, provided by Drumm (2010), corresponding to one individual of Z. normani (for accession numbers see Table 1).For this reason, we sequenced the same genes, i.e.COI, 28S and H3.
Total genomic DNA was extracted from 27 Zeuxo spp.specimens (Table 1) using the JETQUICK Tissue 1 ♀, unknown DNA kit, following the method described in Larsen and Froufe (2010), with the final elution step done with 30 µl.Additionally, two individuals of Tanais dulongii (Audouin, 1826) were also sequenced for this manuscript to be used as outgroups.All the molecular procedures, i.e.DNA extraction, PCR conditions, primers used and sequencing, were the same as those described below for Zeuxo spp.Primers used in both amplification and sequencing were LCO1490 and HCO2198 (Folmer et al. 1994) for the COI gene, 28S-RD1.3fand 28S-rD4b (Whiting 2002) for the 28S gene, and H3AF and H3AR universal primers (Colgan et al. 1998) for the H3 gene.
The PCR conditions (25-µl reactions) were as follows: each reaction contained 2.5 µl 10× Invitrogen PCR Buffer, 0.5 µl 10 mM of each primer, 1.5 µl 50 mM MgCl 2 , 0.5 µl 10 mM dNTP's, 0.1 µl Invitrogen Taq DNA Polymerase, 3µl DNA template and water to make up the 25-µl reaction.The cycle parameters were initial denaturation at 94°C for 3 min, denaturation at 94°C (30 s), annealing at 52°C (45 s) for both coding genes, 56°C (45s) for 28S and extension at 72°C (45 s) repeated for 38 cycles, with a final extension time of 5 min at 72°C.Amplified DNA templates were enzymatically purified (ExoSap) and sequenced using the ABI PRISM BigDye Terminator protocols.Sequences were read on an ABI-310 and are available for each species in GenBank (Table 1).

Phylogenetic analyses
Chromatograms were checked manually using ChromasPro 1.41 (technelysium.com.au).Additionally, one individual from Paratanais sp.(for accession numbers see Table 1) from Drumm (2010) and the two sequenced Tanais dulongii were incorporated in the data set, to be used as outgroups.The sequences of 28S were aligned in MAFFT (Katoh et al. 2005) and the sequences of the coding genes, COI and H3, were aligned with ClustalW using Bioedit v. 5.0.9 (Hall 1999).All the resulting alignments were then adjusted manually, but were not found to require additional editing.
The final combined data set alignment was then analysed using the maximum likelihood (ML) and Bayesian inference (BI) methods.The best-fit model of nucleotide substitution evolution under the corrected Akaike information criterion was estimated using JModelTest 0.1.1 (Posada 2008).The generalized time-reversible model plus gamma (GTR+G) was the most appropriate model of evolution for this dataset (G= 0.4320) and was therefore used in the phylogenetic analyses.ML trees were built in PhyML (Guindon and Gascuel 2003) with 1000 bootstrap replicates and searching for the best-scoring ML tree.Phylogenetic BI was performed on MrBayes version 3.1.2(Ronquist and Huelsenbeck 2003).Sequences were partitioned according to genes.Each partition was allowed independent parameters of sequence evolution under the model GTR+G.Analyses started with program-generated trees, with four heated Markov chains with default incremental heating; two independent runs 1.5×10 7 generations long were sampled at intervals of 100 gen-erations, producing a total of 100000 trees.Burnin was determined upon convergence of log likelihood and parameter estimation values using Tracer (Rambaut and Drummond 2007).
Parameters of genetic divergence, i.e. total number of polymorphic and parsimony-informative sites, haplotype (Hd) and nucleotide diversity (Nd), were calculated with DnaSP v.5.10.01 (Librado and Rozas 2009) for each gene separately, for each species and/or for the whole datasets.The estimates of evolutionary divergence over sequence pairs between groups and the average evolutionary divergence over sequence pairs within groups (all using uncorrected p distances) were calculated using MEGA5 (Tamura et al. 2007).Standard error estimates were obtained by a bootstrap procedure (1000 replicates).

Morphological analysis
Between four and six specimens were fully dissected from each locality, including both adult males and females, using chemically sharpened tungsten wire needles, and were examined under a high-powered compound scope.Full descriptions of both genders of the new species as well as re-descriptions of the known species will be published separately in a specialized taxonomical journal.All specimens used for molecular analysis were examined for diagnostic characters.This involved dissections of mandibles, pleopods and uropods, but not the full array of appendages.

Molecular analysis
The final alignment of the combined data set yielded 1226 bp, of which 426 positions were variable and 286 parsimony-informative (98 and 63 for H3; 206 and 168 for COI; and 122 and 55 for 28S) for 27 individuals of Zeuxo species from seven localities (Table 1), and the three outgroup sequences.Both ML and the BI analysis recovered the same tree for the combined H3+COI+28S data that is shown in Figure 1A.Three highly supported clades were retrieved within Zeuxo with a large genetic divergence of 20% for the combined data set (Fig. 1A).The only Z. cf.normani available in GenBank was collected in Japan, no voucher material exists and we doubt the validity of this identification.The sequences from Japan seem almost identical to the new sequenced individuals from Korea, and cluster with another very well-supported sub-clade of Z. holdichi from the Atlantic coast of Spain and France (Clade 1. Fig. 1A).All three Zeuxo specimens from Turkey grouped in Clade 2. The nine specimens from the Macaronesian archipelagoes of Cape Verde and Canaries (Z.exsargasso) grouped together in a single clade (Clade 3) that had strong internal support.The relationships between the three major clades were poorly resolved in both analyses.
Additionally, as the two Z. holdichi individuals from France (Table 1) were only successfully analysed for 28S, the same phylogenetic analyses were performed for this data set alone.The final 28S alignment yielded 564 bp, of which 122 positions were polymorphic and 55 parsimony-informative.Within Zeuxo spp.specimens, there were 11 haplotypes depicted (differing by 1-2 bp), and relatively high levels of nucleotide variability were found (Hd=0.590,Nd=0.063).The Tamura-Nei plus Gamma model (TrN+G) was the most appropriate model of evolution for this dataset (G=0.5690).The results (both ML and BI analyses recovered the same tree) are shown in Figure 1B.Again, three highly supported clades were retrieved within Zeuxo spp: Clade 1 is a politomy and includes all the individuals collected in Korea (Z.sp.A, Korea), France and Spain (Z.holdichi) and also the Z. normani (hereafter termed Z. cf.normani) registered in GenBank (Drumm 2010) from Japan; Clade 2 is composed of all the Z. exsargasso individuals sequenced from both the Canaries and Cape Verde archipelagos, which cluster well into two sub-groups; and Clade 3 includes the three specimens (Z.sp.B) from Turkey.
The Bayesian analysis revealed the same three highly supported monophyletic clades retrieved with 28S, but in the H3 phylogeny (see supplementary material, Fig. S1) Clade 3, which included individuals from Turkey, clustered with 84% support with Clade 2 (Zeuxo exsargasso individuals from Gran Canaria and the Cape Verde).When 28S was used, Clade 3 came basal to the other two clades (Fig. 1B).
H3.The final H3 alignment was 273 bp long and included the Z. cf.normani previously published (Gen-Bank HM016171) and a total of 24 new Zeuxo spp.individuals (Table 1) in a total of 28 sequences.It contained 98 polymorphic and 63 parsimony-informative sites for the total data set, and within Zeuxo spp.specimens only eight haplotypes were depicted and high levels of nucleotide variability were found (Hd=0.759,Nd=0.062).The tree for the H3 gene is given in supplementary material Figure S1.
COI.The final COI alignment was 390 bp and it consisted of 28 sequences (Table 1), which contained 206 polymorphic and 168 parsimony-informative sites.Within Zeuxo spp.specimens, there were seven haplo-types depicted and high levels of nucleotide variability were found (Hd=0.850,Nd=0.033).The tree for the COI gene is given in supplementary material Figure S2.
The overall success rate for DNA extraction was 75%.Remaining specimens are currently kept at the collection in LMCEE.Zeuxo sp.A and B will be formally described in another paper and the remaining material (types) will be returned to natural history museums in their respective countries of origin.

Morphological results
A number of currently considered diagnostic characters for species identification within Zeuxo were found to be unstable, some even within a single individual.1) The number of uropod articles (Fig. 2).The uropodal article number varies between four and five in Z. sp.A from Korea (Fig. 2A, B), and between five and six in Z. exsargasso (Fig. 2C).2) The pleopod setation (Fig. 3).In this study, the pleopod endopod had only one inner seta in Zeuxo exsargasso, which is otherwise described as having two (Sieg 1980: 219.Fig. 61).The same setal character was shown to vary also in Z. holdichi (Fig. 3A, B) and in both this species and Z. exsargasso the spiniform armament of the pereopod carpus also seems to vary.
3) The carapace pigmentation (Figs 4 and 5).The specimens of Z. sp.A from Korea, while originating from three populations, displayed no variation in the carapace pigmentation patterns (Fig. 4A).The same situation was found in specimens of Z. sp.B from Turkey, although originating from one population only (Fig. 4B).However, the pigmentation pattern of Z. exsargasso from the Cape Verde (two populations) and Gran Canaria (one population) displayed large variation in the carapace pigmentation patterns within and between populations (Fig. 5A, B and Fig. 5C, D, respectively), as well as between the archipelagos (Fig. 5A, B versus Fig. 5C, D).Also, the pigmentation pattern for Z. holdichi varies between the French (Bamber 1990:1589, Fig. 1B) and the Spanish populations (Fig. 4C).

DISCUSSION
This study describes, for the first time, the genetic distribution pattern of species within the order Tanaidacea, which have low dispersal potential.Tanaidacean species can thus be prone to maintaining signatures of genetic drift caused by vicariant events, as seen for other organisms (Petit et al. 2003, Pelc et al. 2009).
The genetic analysis indicates that a 'Zeuxo normani' species complex exists (Fig. 1), spanning the Pacific and the North Atlantic.It furthermore indicates that there is a closer genetic relationship between the Pacific species and Z. holdichi than between Z. holdichi and other North Atlantic species.The Zeuxo normani species complex thus includes Z. normani California (Sieg 1980), Z. cf.normani, Japan (Drumm 2010), Z. cf.normani, Australia (Edgar 2008), Z. sp.A, Korea (this study), and Z. holdichi, Spain and France (this study).
The long-distance dispersal found in certain species with otherwise low dispersal abilities, such as Z. exsargasso, is thought to depend on rafting on floating substrata, particularly on algae (e.g. the genus Sargassum) (Sieg 1980).Thus, the general absence of genetic divergence throughout the distribution of Z. exsargasso (albeit not examined over the entire distribution of the species) could suggest that randomness and regional hydrographic patterns can promote high population connectivity via the Gulf Stream.It has been suggested (Thiel and Haye 2006) that there may be convergence zones where multiple rafts arrive, thus favouring high levels of gene flow between distant populations and erasing the genetic signatures of founder effects.Such 'arrival' zones for eastward dispersal are likely to be the Macaronesian archipelagos (Bamber 2012a).
Zeuxo holdichi (part of the 'Z.normani' species complex) may have retreated to refugia in the south (in northern Spain, Brittany and the English Channel) during the last glaciations, as was proposed for many cold-water taxa (e.g.Kettle et al. 2010).However, this species does not seem to be able to recolonize northern latitudes.Furthermore, it appears that a discontinuous interbreeding between the North Atlantic species (Z.holdichi) and the Pacific 'parent' species (Z.sp.A from  Korea and Z. cf.normani from Japan) allow for morphological observable separation of this species (Z.holdichi): this is plausible if we accept that gene flow is in the Pacific-Atlantic direction.In contrast, while the genetic COI divergence (see Table 2) between Z. cf.normani from Japan and Z. sp.A from Korea is on a similar level (10.7%) to that of Z. holdichi (11.2% and 12.5%, respectively) only a marginal morphological difference can be found between these species according to the description by Shiino (1951).Importantly, this shows that regardless of small genetic (COI) divergence, morphological difference can either be present or absent.
The molecular analyses showed, as expected, the same three highly supported clades within Zeuxo spp., although the relationships between them lacked strong support.It also showed that different species of Zeuxo appear to have a genetic divergence in the COI gene, ranging from 10.7% between Z. sp.A (Korea) and Z. cf.normani (Japan) to 28.2% to 28.5% between Z. sp.B (Turkey) and Zeuxo exsargasso (see Table 2).Although few genetic studies have been made on tanaidaceans, a 20% divergence has been reported between several species of Leptochelia (Drumm 2010, Larsen et al. 2012, Larsen and Froufe 2013).
The closer genetic similarity of Z. holdichi from the Atlantic coast of Spain to the Pacific Z. sp.A from Korea (11.2%) and Z. cf.normani from Japan (12.5%), compared with that of the geographically much closer Atlantic species, Z. exsargasso (28.2%-28.5%),indicates the presence of a Z. normani species complex spanning the Pacific and North Atlantic.This close genetic similarity of Z. holdichi to the North Pacific members of the Z. normani complex suggests a previous inter-glacial colonization.We then speculate that the populations of Z. holdichi from the Atlantic coast of Spain and France are remnants of such a circumpolar distribution.
Our results point out some issues on the current diagnostic characters for species of Zeuxo that were previously found to be unstable.First, it is well known that the number of uropodal articles (Fig. 2) can vary within species of Tanaoidea though it has still been accepted as one of the most important characters for species identification since the monograph by Sieg (1980).However, already Kudinova-Pasternak (1989) mentioned variations of this character in Z. beringi, and both Bamber (1990) and Edgar (2008) confirmed this observation in other species of Zeuxo.Similarly, this study, supported by genetic confirmation, proves that this character can be variable within a species and, indeed, even within an individual.A very thorough study on this character in the closely related genus Zeuxoides (Sieg, 1980) has been conducted by Bird (2008), demonstrating a clear correlation between body size (instars) and uropod article number.Our study has shown that variation can even exist between the left and right uropod of the same specimen (Fig. 2C), regardless of the body size.Therefore, the value of the number of uropodal articles as a character for species identification should be applied with caution.We have seen similar cases of uropod article variation within individuals of other tanaidaceans (Larsen et al. 2012, Larsen andFroufe 2013), and this is thus too often to just be considered a rare aberration.Second, the Parazeuxo subgenus diagnosis (Sieg 1980) depends mainly on the pleopod outer setation character and the number of spiniform setae on the pereopodal merus and carpus.This study has shown that the pleopod endopod setal character (Fig. 3) is variable in Z. exsargasso and Z. holdichi, thus making the subgenus diagnosis invalid.In both of these species, the spiniform armament of the pereopod carpus also seems to vary, and should be applied, again, with caution.Third, in recent years, some authors (Edgar 2008, Bamber 2012a, Kakui pers. comm.) have suggested that pigmentation could be used for species identification without the need for dissection.However, our results do not completely support this premise (Figs 4 and 5), at least for Z. exsargasso and Z. holdichi.Since both the morphological and genetic analysis confirmed the conspecificity of these respective specimens, this suggests that some species have variable carapace pigmentation patterns.No specific patterns were recorded relative to gender, while younger developmental stages appear to have the pigmentation concentrated on the anterior part of the carapace.
In summary, the conclusions of this study are the following: 1) The assembly of the northeast Atlantic comprises species that originated from both the North Pacific and the northwest Atlantic; 2) traditional species identification clearly underestimates diversity; and 3) what have been perceived as reliable diagnostic morphological characters can often be variable.
Hence, more detailed morphological, distributional and molecular surveys, with an emphasis on the intertidal assemblage, are needed to answer the questions of wide distributions records of non-swimming species without a dispersal phase.Clearly there are at least two proven possibilities: frequent rafting events (Larsen 2005, Bamber 2012a); and the presences of multiple, morphologically almost inseparable species known as species complexes (Larsen 2001, Bamber 2010, 2012b).

Fig. 1 .
Fig. 1. -A, maximum likelihood (ML) tree inferred using the TVM+G model of sequence evolution showing relationships of Zeuxo spp.from different locations.The tree is rooted using Tanais dulongii and Paratanais sp.Bootstrap support values above 60% for the ML analysis are shown below nodes, and posterior probability values for the Bayesian analysis above nodes (see Materials and Methods).When the three analyses have the same value, they are represented by an asterisk (*).B, analysis for 28S only, including specimens from France.

Table 1 .
-Details of Zeuxo material: collection locations, numbers of haplotypes, and gene sequences used in the present study

Table 2 .
-Mean genetic uncorrected p pairwise COI distances between the main Zeuxo spp.Standard error estimates are shown above the diagonal.