Insights on the origin of invasive copepods colonizing Basque estuaries; a DNA barcoding approach

Background

The introduction of NIS to estuaries and coastal embayment is of great concern. Commercial ships’ ballast water discharge and the northwards progression of species due to the ongoing climate change arise as the main factors explaining the rising occurrence of NIS species in Northern Atlantic waters. In this regard, regular monitoring of the plankton communities is paramount as to be able to respond properly to this potential issue.

Results

While monitoring the invasive copepod Acartia tonsa populations in the estuaries of Bilbao and Urdaibai (Basque country, Spain), we report here the Asian copepod Pseudodiaptomus marinus for the first time in the Iberian Peninsula waters. Individuals from both species were collected from July to October, 2013 for DNA sequencing of mitochondrial cytochrome c oxidase subunit I gene (MT-CO1). Phylogenetic analysis of MT-CO1 confirmed P. marinus identity.

Conclusions

Phylogeographic distribution of A. tonsa haplotypes in Europe along with the Bilbao port traffic patterns suggested a secondary invasion from an European source to Basque estuaries. The successful establishment of the A. tonsa population and the appearance of Pseudodiaptomus marinus confirm the need for regular plankton monitoring of estuarine and port waters. This applies also to nearby systems as these populations could represent a source of future dispersal.

The estuary of Bilbao (43.35° N, 3.03° W), once considered among the most polluted estuaries in Europe, initiated in the late 1980s a water quality restoration program that allowed the re-colonization of the zooplankton community. The surrounding neritic species successfully colonized and became autochthonous but with a polyhaline - euhaline habitat preference (Uriarte & Villate 2004, 2005; Albaina et al., 2009; Aravena et al., 2009). More in detail, the inner part of the system changed from nearly no plankton inhabitants to dominance of the opportunistic Acartia clausi Giesbrecht 1889, and, within a few years, to a diverse “Acartia species complex” fitting in salinity driven species-specific niches (Albaina et al., 2009; Aravena et al., 2009). However, since 2003 the invasive species Acartia tonsa Dana 1849 (first cited in 2001) displaced the rest of congeneric species becoming predominant in inner waters of the Bilbao estuary (Aravena et al., 2009). The combination of empty niches, suitable environmental conditions, and availability of proper vectors (high vessel traffic with massive ballast water exchanges, Butrón et al., 2011) made the estuary of Bilbao a perfect recipient for Non Indigenous Species (NIS) (Paavola et al., 2005). Today, this estuary is among the main ports in the North-East Atlantic (NEA) with connections to > 800 ports worldwide. Shortly after its dominance in the Bilbao estuary, the nearby estuary of Urdaibai (43.41° N, 2.69° W), where A. bifilosa Giesbrecht 1881 had been dominating since 1980s (e.g. Uriarte & Villate 2005), was also colonized by A. tonsa (Aravena 2009). And since its first appearance in 2003, this species has showed a huge increase in number, becoming a codominant species nowadays (data not shown).

The introduction of NIS to estuaries and coastal embayment is of great concern and regular monitoring of these particularly sensitive ecosystems for NIS and its ecological impact is paramount for a timely answer and any further mitigation/adaptation policy design. In the present study, we assessed zooplankton communities in the inner part of the estuaries of Bilbao and Urdaibai during the summer of 2013. Our objectives were to 1) examine the genetic information of the known NIS A. tonsa in both estuaries and gain insights into the invasion source by comparing mitochondrial DNA (mtDNA) sequences obtained from the native and introduced ranges (e.g. Downie 2002; Cristescu 2015 for a review) and 2) discover any novel invasive species to these Basque estuaries.

The zooplankton sampling was carried out monthly, from July to October, in consecutive days during neap tides in waters with a salinity of ~30. Samples were obtained from horizontal tows (below halocline) of a 200 μm mesh plankton net (mouth diameter: 0.25 m) equipped with a flowmeter (Hydrobios). Net samples were transferred to the laboratory where subsamples were preserved either in 4 % borax buffered formalin seawater solution for taxonomic and abundance analysis or in absolute Ethanol for genetic analyses. The qualitative and quantitative analysis of zooplankton was carried out under a stereoscopic microscope (Olympus SZX12). Identification was made to species or genus level following mainly the ICES Identification Leaflets for Plankton (http://goo.gl/qfKg9V). In each sample, a minimum of 100 individuals of the most abundant taxa were counted before finishing sub-sampling; if not possible, the whole sample was identified. An Olympus IX70 microscope was used for Pseudodiaptomus marinus Sato 1913 identification confirmation. This species was identified following the keys in Grindley & Grice (1969) and Brylinski et al. (2012). Five to ten individuals per estuary (A. tonsa and P. marinus) were sorted for DNA sequencing of mitochondrial cytochrome c oxidase subunit I gene (MT-CO1).

DNA was extracted from single copepods using the Chelex method (following Chen & Hare 2008). Amplification of the MT-CO1 region was performed using Folmer et al. (1994)’s primers (LCO1490 and HC02198). Barcode region was PCR-amplified in a 25 μl reaction containing 2.5 μl of each primer (10 μM), 2.5 μl dNTP (0.2 mM of each dNTP), 2.5 units of GoTaq polymerase (5 u/μl; Promega), 2.5 μl MgCl2 (25 mM), 0.5 μl BSA (10 mg/ml), 5 μl 5X Buffer, 8 μl Milli-Q water (Merck) and 1 μl of DNA extract. Reactions were amplified trough 35 cycles (30s at 95 °C, 45 s at 57 °C and 45 s at 72 °C) followed by a final extension step of 7 min at 72 °C. PCR products were purified with Multiscreen HTS plates (Millipore) and further sequenced in both directions on an Applied Biosystems (ABI) 3139 capillary electrophoresis Analyzer, with ABI BigDye Terminator version 3.1 chemistry (Applied Biosystems). Base-calling from chromatograms was performed using Seqscape v2.5 (ABI). The BLASTN algorithm was used to verify that the target locus had been amplified and to discard the putative presence of pseudogenes and/or contamination. The herein generated sequences (GenBank accession numbers KM879190-KM879193 and KM879199-KM879216) were combined with 275 and 36 MT-CO1 sequences from A. tonsa and Pseudodiaptomus spp., respectively, obtained after an exhaustive GenBank search (Table 1). Both A. tonsa and Pseudodiaptomus spp. sequences were then aligned using Clustal W algorithm as implemented in Bioedit (Hall, 1999). The alignment length was 481 and 572 bp, respectively. From the MT-CO1 alignments, maximum likelihood (ML) trees were constructed with MEGA5 software (Tamura et al., 2011) following Hall (2013) guidelines. The best-fit model of evolution, T92 + G and HKY + G + I for, respectively, A. tonsa and Pseudodiaptomus spp. datasets (MEGA5 software), and 500 bootstrap replicates were selected. Phylogenetic trees were then graphically represented using FigTree v1.4.2 software. DnaSP v5 software was used to calculate haplotype and nucleotide diversities (Librado & Rozas 2009). We then inferred a haplotype network of MT-CO1 sequences using a median-joining algorithm implemented in Network v4 software (Bandelt et al., 1999).

Finally, as commercial ships’ ballast water exchange is generally considered as the main vector for NIS introduction to the marine environment, Bilbao port ship traffic data from the 2001–2013 period was obtained from the port authorities. The number of ships entering the Bilbao port and their previous port of origin were considered.

Copepod species composition and abundance patterns of summer 2013 were examined for two Basque estuaries (Table 2). While up to four Acartia species were recorded in the estuary of Urdaibai in summer 2013 (namely Acartia tonsa, A. bifilosa, A. clausi and Paracartia grani Sars 1904), A. tonsa was the sole representative in the estuary of Bilbao. It was the most abundant copepod in the latter except for July when Calanipeda aquaedulcis Kritchagin 1873 was. Pseudodiaptomus marinus, with no previous records in the estuary of Bilbao, was the third most abundant copepod in July 2013 and was recorded again, in a lower number, in September. This species was absent in the Urdaibai estuary where A. tonsa and A. bifilosa dominated the copepod community with comparable abundances.

A total of 18 and 4 MT-CO1 sequences were obtained for, respectively, A. tonsa and P. marinus. The latter showed a high level of similarity (98 %; BLASTN algorithm) to the sole GenBank sequence for P. marinus that was collected from Korea (Accession no. AY145436). All five sequences were clustered together in the phylogenetic tree for the Pseudodiaptomus genus (Fig. 1a). Apart from this, a total of three MT-CO1 haplotypes were recorded in the estuary of Bilbao, showing two to four variable sites among them and twelve against the sole available Asian haplotype (Fig. 1b).

Pseudodiaptomus marinus. a Maximum likelihood (ML) tree constructed with the HKY + G + I model of evolution. Forty MT-CO1 sequences (572 bp) from 8 different Pseudodiaptomus species were included in the analysis. The Bilbao estuary sequences (P. marinus AG39 – P. marinus AG42) clustered together in the lower branch (full black circle) along with the GenBank sole P. marinus record (grey circle). The specimen signalled by an asterisk that clustered within P. koreanus ones was identified as a new species (Pseudodiaptomus cf. inopinus Kyushu) in the original GenBank record (Accession no. AB576159) and might represent a misidentification. The scale bar indicates substitutions per nucleotide position. b Median-Joining (MJ) network analysis for the five P. marinus MT-CO1 sequences. Numbers correspond to the location of variable sites along the assembled region; white rhombuses represent missed haplotypes. Each circle represents one MT-CO1 haplotype; the size of the circles corresponds to the number of individuals sharing a haplotype (one except for a sole case with two individuals)

Full size image

Regarding A. tonsa, three MT-CO1 clades were evident in the ML tree where the Basque estuaries’ individuals were grouped with all the GenBank sequences (Fig. 2a). This included both the species’ native area (including samples along the Eastern USA coast) and European records (n = 293; Table 1). A total of 102 distinct MT-CO1 haplotypes were included in this analysis, eight of them with reported presence in European waters. Interestingly, all the European registers (n = 55) clustered together in clade-A. While the most common European haplotype (H-1) was shared by forty-four individuals, H-4, H-6, H-7 and H-8 haplotypes had a sole representative. Finally, two haplotypes were recorded in Basque estuaries. While H-1 was present in both estuaries, H-2 was only recorded at the estuary of Urdaibai (Fig. 2b).

Acartia tonsa. a Maximum likelihood (ML) tree, constructed with the T92 + G model of evolution, for a 481 bp MT-CO1 region including 293 A. tonsa sequences. Only A. tonsa specimens collected at European locations (n = 55) are labelled with both the GenBank accession number and the geographical location (NSEA and MEDIT stand for, respectively, North Sea and Mediterranean Sea). Basque estuaries’ specimens were labelled as BILBAO 1-10 and URDAIBAI 1-8. A sole A. hudsonica Pinhey 1926 sequence (EU274431, highlighted) was included as an out-group. The scale bar indicates substitutions per nucleotide position. b Median-Joining (MJ) network analysis for the European collected A. tonsa individuals. MJ network constructed as in Fig. 1b. Colours corresponded to the different geographical areas (legend superimposed). Numbers correspond to the location of variable sites along the assembled region

Full size image

In the summer of 2013 the copepod community of two Basque estuaries (namely the estuary of Bilbao and the estuary of Urdaibai) was monitored to discover any non-indigenous species and any changes of known NIS Acartia tonsa, with a goal to gain insights into NIS invasion patterns following a DNA barcoding approach (e.g. Downie 2002). Apart from the expected presence of abundant A. tonsa in both estuaries, we surprisingly recorded Pseudodiaptomus marinus for the first time in the estuary of Bilbao as well as the entire Iberian Peninsula waters. The MT-CO1 gene was selected because it is widely used as the standard DNA barcoding gene with a high species-discriminating power and a vast amount of publicly available sequences.

Both morphological characters and barcode sequences confirmed the presence of P. marinus (Fig. 1a), a species considered to be native to the Northwestern Pacific Ocean and that had been recently cited in NEA waters (see Brylinski et al., 2012). Interestingly, although only four individuals were sequenced, three distinct haplotypes were recorded (Fig. 1b). In comparison, we found one A. tonsa haplotype within the Bilbao (n = 10) and two in Urdaibai (n = 8) (Fig. 2b). Extremely low mtDNA haplotypic diversities are common in invasive species due to a strong founder effect (e.g. Holland 2000; Chandler et al., 2008). However, the degree of reduction depends on several factors such as the number (and scale) of introduction events and the diversity of the invasion origin (Cristescu 2015). However, the lack of P. marinus MT-CO1 records in the public databases prevents any discussion. Concerning this species, the first record in European waters was in 2007 in the north Adriatic Sea (De Olazabal & Tirelli 2011) and since then P. marinus has been cited in other locations along the Mediterranean Sea (Sabia et al., 2014 and references therein). Regarding North-East Atlantic (NEA) waters, it has been only reported in the southern North Sea in 2010 and in 2011, in two independent surveys (Brylinski et al., 2012; Jha et al., 2013). Taken into account an unpublished record at the Gironde estuary (Brylinski et al., 2012), to our knowledge, present work represents the fourth citation for this NIS in NEA waters and the first one in the Iberian Peninsula waters. However, the 243 ind. m^-3 abundance for P. marinus in July 2013 in the estuary of Bilbao is unprecedented as previous NEA records were always below 4 ind. m^-3. Moreover, the presence of females carrying eggs (5 out of 37 individuals sorted in July) suggested an established P. marinus population in the estuary of Bilbao. Finally, further analysis of archived samples (2005 to 2012 period; monthly surveys in waters with a salinity of ~30) allowed us to date the first recorded presence of this species back to October 2010 (10/30/2010, 28 individuals with no ovigerous female).

On the other hand, A. tonsa, native to the American and Indo-Pacific area, is a well-known NIS in European waters since the early 1900s (Brylinski 1981), with ships’ ballast water exchange proposed as the most likely introduction vector (e.g. Paavola et al., 2005). In summer 2013, A. tonsa was consistently more abundant (~2000 ind. m^-3, considering only adults; Table 2) in the estuary of Bilbao, where it was the only species of the genus Acartia, than in the estuary of Urdaibai, where this species was at the same abundance levels of A. bifilosa. The rich amount of MT-CO1 sequence data in GenBank for A. tonsa allowed us to compare the herein obtained MT-CO1 sequences against a large amount of representative sequences from both the species’ native and invaded area (Table 1, Fig. 2). Chen & Hare (2011), covering 20 estuarine systems along the whole USA Atlantic coast, reported three sympatric genetic lineages for A. tonsa (namely X, F and S) for both a mitochondrial and a nuclear DNA loci, respectively, MT-CO1 and ITS (Internal Transcribed Spacer). These three clades are also shown in Fig. 2a. Every European record, including a total of eight different haplotypes, clustered together in one clade (“A” in Fig. 2a) that corresponded to the “X” lineage in Chen & Hare (2011). The fact that, to date, only one of the three A. tonsa native lineages has been reported in European waters deserves further discussion. While an insufficient sequencing effort in European waters compared with North American ones could represent a putative explanation (Table 1), it is also likely that an undetermined physiological pre-condition such as a broader tolerance range in lineage “X” has allowed the colonisation of European waters by this lineage. In this sense, Chen & Hare (2011) reported a higher gene flow corresponding to the “X” lineage in North American estuaries including preference for a higher and larger range of, respectively, salinities and latitudes than the other two native lineages. This could represent a higher resilience capacity and an obvious advantage when facing longer displacements such as trans-oceanic ones and, more importantly, for successful settlement in the invaded systems. Environmental stress associated with ballast water environment and/or a better adaptation to European water conditions (either for adults or resting eggs) in those individuals could explain the reported pattern. In this sense, the lower haplotype diversity (1 haplotype out of 10 individuals) in Bilbao may be related to the selective effect of pollution, whereas higher haplotype diversity (2 out of 8) to the relatively pristine condition of Urdaibai, an UNESCO Biosphere reserve since 1984, with mainly recreational use and a relatively low urbanization pressure. However without further experiments, such as common-garden ones, this questions will remain unsolved.

Interestingly, the two haplotypes recorded in the Bilbao and Urdaibai estuaries (H-1 and H-2 in Fig. 2b) had been previously reported at both the native and invaded European locations. While H-1, the dominant haplotype in Europe (including the ten specimens of the estuary of Bilbao along with six out of eight from the estuary of Urdaibai; Fig. 2b), was also the predominant haplotype in the native area’s lineage “X”, the H-2 had been previously cited at three estuaries in the Eastern USA coast (Gang Chen, personnel communication) and at an undetermined Italian location (GenBank accession number HE647798). More interestingly, the haplotype network for European A. tonsa MT-CO1 records, showing closely related haplotypes in Basque waters but a well-diversified pattern in North Sea area, suggests the latter as putative center of origin (Fig. 2b). This fits also the chronology of A. tonsa’s records in NEA waters, with earliest records centered in the North Sea area (Figure 1 in Brylinski 1981). Although non-European origins of A. tonsa invaded to Basque estuaries cannot be ruled out, both the previously cited presence of H-1 and H-2 within European waters and the haplotype network shape points to a secondary introduction from an European source as the most likely process. Moreover, this is further supported by the Bilbao port maritime traffic patterns when considering the putative source of NIS. In this sense, during the 2001–2013 period, from an average of 3400 commercial ships per year that entered the port of Bilbao, only 3 and 2 % of the ships arrived from a North American and Asian (Indian and Pacific) port of origin, respectively, whereas 53 % arrived from Atlantic European ports.

Paavola et al. (2005) have cited that low species richness environments are more prone to NIS successful establishment and this is supported by present data showing the current preponderance of NIS in the Bilbao estuary following pollution abatement. Shortly after A. tonsa, P. marinus also invaded the system and both became settled in a few years’ time (Aravena et al., 2009, Uriarte et al., 2016). Both species had been recorded earlier in nearby NEA systems and, corresponding to a broad range of tolerance to environmental stressors such as salinity and pollution, are capable of proliferating in harbours and other anthropized environments (see Brylinski 1981 and Brylinski et al., 2012 for, respectively, A. tonsa and P. marinus). While data are still scarce, A. tonsa and P. marinus previous records’ pattern along with present data suggest an ongoing colonization process in NEA waters. More surveillance of plankton environment, especially at potential NIS sources such as harbours or marinas is critical to anticipate further spreading of these species. Finally, the development of multi loci markers such as SSRs or SNPs for these NIS species is desirable as they have the potential to determine a more precise geographical origin of marine invasions (e.g. Reusch et al., 2010).

• Pseudodiaptomus marinus was reported for the first time in the Iberian Peninsula waters.

• Previous records’ pattern along with present data suggest an ongoing colonization process in NEA waters.

• Phylogeographic distribution of Acartia tonsa, haplotypes in Europe along with the Bilbao port traffic patterns suggested a secondary invasion from an European source to Basque estuaries.

• The successful establishment of the A. tonsa population and the appearance of Pseudodiaptomus marinus confirm the need for regular plankton monitoring of estuarine and port waters.

DA’s work was supported by the ZabaldUz Program (PhD fellowship). SGIker technical and human support (UPV/EHU) is gratefully acknowledged. Thanks to T. Matellanes (Bilbao Port Authority) for providing maritime traffic data. Special thanks to the anonymous reviewers that greatly improved the different manuscript versions. The sampling program received funds from the University of the Basque Country (UFI11/37) and the Basque Government (GIC10/168).

Availability of data and materials

The dataset supporting the conclusions of this article is available in the GenBank repository (http://www.ncbi.nlm.nih.gov/nuccore/; Table 1).

Authors’ contributions

AA, FV and AE conceived and designed the work. AA, MA, DA, IU and AI designed the sampling and participated in laboratory work. MA, DA and IU collected the samples. AA performed the data analysis. AA wrote the manuscript that was read, edited and approved by the rest of authors.

Competing interests

The authors declare that they have no competing interests.

Ethics approval and consent to participate

Sampling copepods in Spain did not require any specific consent or permission. Animals did not suffer as they were transferred to Ethanol before sorting for DNA extraction.

Affiliations

1. Department of Genetics, Physical Anthropology & Animal Physiology, Faculty of Science and Technology, Laboratory of Genetics, University of the Basque Country UPV/EHU, Leioa, 48940, Spain
   • A. Albaina
   • , M. Aguirre
   • , D. Abad
   •  & A. Estonba
2. Department of Plant Biology and Ecology, Faculty of Pharmacy, Laboratory of Ecology, University of the Basque Country UPV/EHU, Gasteiz, 01006, Spain
   • I. Uriarte
   •  & A. Iriarte
3. Department of Plant Biology and Ecology, Faculty of Science and Technology, Laboratory of Ecology, University of the Basque Country UPV/EHU, Leioa, 48940, Spain
   • F. Villate
4. Environmental Studies Centre (CEA), Vitoria-Gasteiz city council, Vitoria-Gasteiz, 01008, Spain
   • A. Albaina

Corresponding author

Correspondence to A. Albaina.

Reprints and Permissions