Comparative phylogeography and demographic history of five sibling species of Pseudocalanus (Copepoda: Calanoida) in the North Atlantic Ocean

Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

Contents lists available at ScienceDirect

Journal of Experimental Marine Biology and Ecology journal homepage: www.elsevier.com/locate/jembe

Comparative phylogeography and demographic history of five sibling species of Pseudocalanus (Copepoda: Calanoida) in the North Atlantic Ocean Ole Nicolai Staurland Aarbakke a,⁎, Ann Bucklin b, Claudia Halsband c, Fredrika Norrbin a a b c

Department of Arctic and Marine Biology, UiT- The Arctic University of Norway, 9037 Tromsø, Norway Department of Marine Sciences, University of Connecticut, Storrs, CT 06269, USA Akvaplan Niva AS, Fram Centre, 9296 Tromsø, Norway

a r t i c l e

i n f o

Article history: Received 4 August 2014 Received in revised form 2 October 2014 Accepted 4 October 2014 Available online 21 October 2014 Keywords: Zooplankton Phylogenetic inference Evolutionary neutrality Demographic inference Bayesian skyline

a b s t r a c t Species of the genus Pseudocalanus (Copepoda: Calanoida) are widespread in the North Atlantic Ocean and provide an opportunity to explore ocean-basin-scale genetic relationships of closely related holozooplanktonic species. A phylogeographic analysis of five species of Pseudocalanus, using two mitochondrial DNA markers and one nuclear DNA marker, was conducted. The results show that there are differences in the overall genetic structure among the five species. Three species, P. minutus, P. moultoni and P. elongatus displayed no or little genetic structuring, while P. acuspes and P. newmani showed comparatively higher genetic structure. These results, combined with tree-based and demographic-history analyses, may indicate that the five copepod species investigated herein could represent two diverging evolutionary branches. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Planktonic copepod species exhibit a broad range of phylogeographic patterns, from genetically homogenous populations throughout their ranges to marked geographic structuring. In some cases, sibling-species pairs of copepods have been shown to exhibit very different patterns of distribution of genetic variation, with one species showing panmixia while another shows marked differentiation of geographic populations (Blanco-Bercial et al., 2011; Durbin et al., 2008). The evolutionary and ecological processes that underlie and determine the patterns of distribution of copepods and other planktonic marine species can be inferred from analysis of the patterns of genetic variation (Knowles, 2009). Phylogeographic analysis (i.e., the study of geographic patterns of genetic variation) may be used to recognize and distinguish among the different processes that may contribute to the present-day observable patterns of distribution and abundance. Many putative cosmopolitan species may comprise genetically distinct populations (e.g., Durbin et al., 2008) or cryptic species with discrete biogeographical distributions (e.g., Lee, 2007). Phylogeographic studies of holozooplankton have shown that changes in ocean current circulation patterns over geological time ⁎ Corresponding author at: Dept. of Arctic and Marine Biology, Faculty of Biosciences, Fisheries and Economy, UiT – The Arctic University of Norway, PO-Box 6050 Langnes, 9037 Tromsø, Norway. Tel.: +47 77 64 60 32. E-mail address: [email protected] (O.N.S. Aarbakke).

http://dx.doi.org/10.1016/j.jembe.2014.10.006 0022-0981/© 2014 Elsevier B.V. All rights reserved.

(e.g., Blanco-Bercial et al., 2011) and species-specific ecological differences (e.g., Goetze, 2005; Goetze and Bradford-Grieve, 2005) can cause genetic differentiation of populations and potentially lead to reproductive isolation and speciation. Behavioral patterns, such as diel vertical migration, may also cause genetic differentiation of geographic populations of open ocean planktonic copepods (Goetze, 2011). Climate variability (e.g., global warming), environmental impacts (e.g., acidification of the oceans), and species interactions (e.g., competition with other species) may cause population size changes over time in marine animal species. Information about ancestral population-size and growth parameters can be inferred from analysis of genetic variation using coalescent theory (i.e., tracing all alleles of a gene shared by all members of a population to the most recent common ancestor, Griffiths and Tavare, 1994; Hudson, 1990; Kingman, 1982). Persistent populations should exhibit high genetic diversity, with high rates of private haplotypes and show signs of demographic stability. Conversely, populations displaying low diversity, relatively fewer haplotypes and a demographic change signal could be recently established or have undergone a population expansion or reduction. Several statistical tests have been developed that examine occurrence of past demographic changes (e.g., Tajima's D (Tajima, 1989), Fu and Li's D and F (Fu and Li, 1993), and Fu's Fs (Fu, 1997)). Through a comparison of observed variability (measured as heterozygosity, haplotype or nucleotide diversity, and numbers of pairwise differences or segregating sites) these tests determine whether haplotype frequencies deviate from evolutionary neutrality (Kimura, 1985). Both extensive

480

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

gene flow and population range expansions can create similar patterns in mitochondrial sequence data. Thus, an observed lack of isolation by distance and/or spatial partitioning of molecular variance do not necessarily imply high gene flow and contemporary panmixia (Venkatesan et al., 2007). Because the neutrality tests explore slightly different properties of the data, careful examination of the combined results of these statistics allows elucidation of the most-likely demographic history. Bayesian skyline analysis is a method that uses contemporary genetic diversity within a given population to infer past changes in population size. Under coalescent theory, the population size at a given point in time is inversely proportional to the probability of two lineages coalescing during that generation (Drummond et al., 2005).

in the English Channel. Finally, Holmborn et al. (2011) found the genetic diversity of P. acuspes within the Baltic Sea to be unusually low, finding only two haplotypes in 83 individuals. This study assesses the phylogeography and demographic history of five species of Pseudocalanus in the North Atlantic using two mitochondrial markers and one nuclear marker. The findings are important for future phylogeographic studies of holozooplankton in the North Atlantic Ocean, as well as biogeographic studies of Pseudocalanus species.

1.1. Genus Pseudocalanus

Samples for this study were collected from 10 locations spanning the geographic distribution of five species of Pseudocalanus. Because Pseudocalanus species are difficult to identify, stations were chosen in known and suspected areas of sympatry, in order to obtain as many species as possible in several locations. The majority of the biological material analyzed in this study was obtained from the Census of Marine Zooplankton (CMarZ) archives located at the Department of Marine Sciences, University of Connecticut, USA (Table 1). All samples had been treated in accordance with the protocols described by Bucklin (Bucklin, 2000).

Calanoid copepods of the genus Pseudocalanus are abundant and widespread in the North Atlantic Ocean, where the genus is represented by at least five species (Corkett and McLaren, 1978; Frost, 1989). Despite extensive sympatry, overlapping species ranges (Fig. 1), and remarkably similar morphology, species of Pseudocalanus display a variety of life histories and ecological preferences. Among both species and individuals, there may be significant differences in timing of reproduction, number of generations, overwintering strategies and behavior (e.g., Carter, 1965; Halsband and Hirche, 2001; Lischka and Hagen, 2005; McLaren et al., 1989; Norrbin, 1991, 1994; Renz et al., 2007, 2008; Yamaguchi et al., 1998). Genetic differences between species of Pseudocalanus have been established for the glucose phosphate isomerase (GPI; Sévigny et al., 1989) locus and mitochondrial cytochrome oxidase I (COI) (Bucklin et al., 2003). Several studies have probed within-species genetic diversity of the genus. Sévigny et al. (1989) found significant differences in allelic frequencies at the GPI locus between Pacific and Atlantic populations of P. newmani and P. moultoni, but no significant differences among populations within ocean basins for P. newmani, P. moultoni, P. acuspes or P. minutus. Using COI, Aarbakke et al. (2011) showed divergence between Pacific and Atlantic populations of P. moultoni, and found smallbut-significant genetic distances between P. moultoni populations in the NW and NE Atlantic. Unal et al. (2006) reported differences in levels of intraspecific variation for COI among P. elongatus populations, with low variation in the Black Sea and comparatively high genetic diversity

2. Materials and methods 2.1. Sample collection

2.2. DNA extraction, PCR amplification, sequencing and sequence alignment From each sample jar, mature female (stage CVIF) Pseudocalanus spp. were picked using a dissecting microscope and placed in a separate vial with 96% ethanol. In instances when mature females were rare or absent, mature males (CVIM) or juvenile fifth copepodites (CV) were used. Individuals used in the analysis were selected at random from each vial (i.e., using a plastic pipette and no magnification). The specimen was then submerged in a large petri dish filled with MilliQ purified water to remove the alcohol prior to DNA extraction. Total genomic DNA was extracted using a DNeasy Blood and Tissue kit (Qiagen, Valencia, CA, USA) with a final elution of 90 μl. Three genes were selected for this study: mitochondrial cytochrome oxidase subunit I (COI), mitochondrial cytochrome B (CytB) and the nuclear internal transcribed spacer region 1 of ribosomal DNA (ITS-1). The consensus invertebrate primers LCO-1490: 5′-GGTCAACAAATCATAAAG

Fig. 1. Map of the North Atlantic Ocean showing the distribution of Pseudocalanus spp. (Aarbakke et al., 2011; submitted; Bucklin et al., 1998, 2001; Frost, 1989; Grabbert et al., 2010; Holmborn et al., 2011; Laakmann et al., 2013; Unal et al., 2006).

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488 Table 1 Overview of sampling locations of Pseudocalanus spp. with numbers of COI identified individuals of five species in parentheses. Species

Collection location (n)

P. acuspes P. newmani P. minutus

NE Atlantic: Balsfjord (52), Håkøy (1), Barents Sea (2) NW Atlantic: Gulf of Maine (18), Georges Bank (2), Nova Scotia (4) NE Atlantic: Norwegian Sea (36), Håkøy (5), Balsfjord (7), NC Atlantic: Iceland20 (51), Iceland21 (39) P. moultoni NW Atlantic: Georges Bank (30), Nova Scotia (1), NE Atlantic: Håkøy (9), Svartnes (10), NC Atlantic: Iceland (1) P. elongatus NE Atlantic: East English Channel (13), NC Atlantic: Iceland (8)

ATATTGG-3′ and HCO-2198: 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′ Folmer et al. (1994) were used to amplify COI. The PCR protocol was 94 °C for 30 seconds, 50 °C for 40 seconds and 72 °C for 40 seconds, with an initial annealing step of 3 minutes and a final elongation step of 10 minutes. COI was amplified for all 385 individuals for species identification. A second mitochondrial gene, CytB, was amplified using the primers 151 F 5′-TGTGGRGCNACYGTWATYACTAA-3′ and 270R 5′AANAGGAARTAYCAYTCNGGYTG-3′ (Merritt et al., 1998). The PCR protocol was 94 °C for 40 seconds, 50 °C for 45 seconds and 72 °C for 45 seconds, with an initial annealing step of 3 minutes and a final elongation step of 15 minutes. Total reaction volume was 25 μl and the PCR was run for 35 cycles. The nuclear marker, ITS 1, was amplified using the primers F1665-18S 5′-CCGTCGCTACTACCGATTGAACG-3′ and R73-5.8S 5′-GTGTCGATGTTCATGTGTCCTGC-3′ (primers made by R. J. Machida, University of Tokyo in Figueroa, 2011). The PCR protocol used for ITS-1 was identical to that for CytB. The total reaction volume and run duration for all PCR reactions was 25 μl and 35 cycles, respectively. For all three genes, the PCR products were purified with either PCR gel extraction kit or MinElute PCR purification kit (both Qiagen) depending on the strength of the bands when they were visualized on a 1% agarose gel. Cycle sequencing was performed on a ABI 3130 genetic analyzer using the BigDye ver. 3.1 cycle sequencing kit (Applied Biosystems, Inc., Carlsbad, CA, USA). For all three genetic markers the sequencing protocol was 96 °C for 10 seconds, 50 °C for 5 seconds and 60 °C for 4 minutes for a total of 35 cycles. All sequences were aligned in ClustalW (Thompson et al., 1994) and the raw files were compared to the sequences to verify the quality of the base calling. COI sequence variation was used to identify the species through BLAST searches of the GenBank database (Altschul et al., 1997). 2.3. Population genetic analysis 2.3.1. Population structure Hierarchical analysis of molecular variance (AMOVA; Excoffier et al., 1992), population pairwise genetic distances and Wright's FST (Weir and Cockerham, 1984; Wright, 1951) values were calculated in Arlequin ver. 3.5.1.2 (Excoffier and Lischer, 2010). To detect significant deviation from the hypothesis of no difference between populations, 50,000 haplotype permutations were performed and the corresponding P-values were subject to a sequential Bonferroni correction (Holm, 1979) at 0.0001, 0.01 and 0.05 significance levels. 2.3.2. Tree-based analysis Tree-based analyses were performed in Molecular Evolutionary Genetic Analysis (MEGA) Ver. 5 (Tamura et al., 2011). Statistical selection of the best-fit models for each dataset, among 88 candidate models, was done for Bayesian Information Criterion (BIC), Akaike Information Criterion (AIC), AIC corrected for sample size (AICc) and Decision Theory (DT) optimized for Maximum Likelihood using jModelTest (Posada, 2008). 2.3.3. Network analysis Statistical parsimony COI and CytB networks (Clement et al., 2000; Templeton et al., 1992) with a 95% cutoff were constructed for all

481

species to visualize within-species relationships of haplotype frequency and distribution. 2.3.4. Nucleotide diversity and tests of neutrality Haplotype polymorphism and nucleotide diversity were calculated in DNasP ver. 5 (Librado and Rozas, 2009), estimated according to Nei (1987) as the average number of nucleotide differences per site between two sequences. To test whether the haplotype frequencies deviated from evolutionary neutrality, a series of statistical tests was performed using coalescent analysis. These were Tajima's D (Tajima, 1989), Fu and Li's D and F (Fu and Li, 1993) and Fu's Fs (Fu, 1997). Through a comparison of observed variability (measured by heterozygosity, haplotype or nucleotide diversity, and number of pairwise differences or segregating sites), these tests determine whether haplotype frequencies deviate from evolutionary neutrality. Strobeck's S Strobeck (1987) was calculated, which compares the observed haplotype number to the expected number of haplotypes based on the frequency distribution derived from the inferred mutation rate. 2.3.5. Demographic inference Prior to the demographic analyses, the sequence data were tested for substitution saturation using the DAMBE (data analysis in molecular biology and evolution) software package ver. 5.3 (Xia and Xie, 2001). Bayesian skyline plots for each species were generated using all COI sequence data in BEAST ver. 1.7.5 (Drummond et al., 2012) and visualized in Tracer ver. 1.5 (Rambaut and Drummond, 2007). The substitution model used was HKY with 4 gamma categories applying the strict molecular clock, according to AICc score obtained from JModelTest. Three partition designs were used: no partition and two and three partitions of codon positions 1 + 2 and 3 and 1 + 2 + 3, respectively. The Markov Chain Monte Carlo (MCMC) length was set to 107 with a burn-in of 10%. To our knowledge, no direct estimate of mitochondrial mutation rate exists for calanoid copepods. Therefore, an estimate suggested by Blanco-Bercial et al. (2011) for Clausocalanus of 1.2–1.5% per million years was used. This estimate is based on a range of mutation rates from Knowlton and Weigt (1998) and Lessios (2008). 3. Results 3.1. Cytochrome oxidase I Mitochondrial COI sequence data from the five species of Pseudocalanus comprised 290 sequences and resulted in an alignment length of 473 base pairs (b-p). No insertions or deletions were detected. The alignment contained 177 polymorphic sites, of which 147 were parsimony informative, and 74 haplotypes were discerned (Table 2). The model selected for the phylogenetic analysis, based on the AIC and AICc scores from JModeltest, was Hasegawa, Kishino and Yano (HKY; Hasegawa et al., 1985) with gamma correction. The tree-based analysis showed that all five species displayed high levels of intraspecific genetic structure, a pattern most pronounced for P. acuspes and P. newmani (Fig. 2) Analyses of molecular variance (AMOVA) were done by grouping the populations of each species into Northeast, North Central and Northwest Atlantic. For P. minutus, P. moultoni and P. newmani, 100, 97.1 and 92.4% of the variance was attributed to differences within populations, respectively. Conversely, Pseudocalanus acuspes and P. elongatus showed higher levels of genetic structuring, having 68.3% (P b 0.05) and 16.4% (P b 0.0001) of the variance among populations. Minimum spanning networks were built using statistical parsimony with a 95% connection limit to visualize spatial distribution and relative frequency of the haplotypes within each of the five species (Fig. 5). Three species, P. newmani, P. elongatus and P. minutus, displayed a “star shape” with several less common haplotypes radiating from one major haplotype. The haplotype network for P. moultoni also showed one major haplotype (found in all sampling locations except Nova

482

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

Table 2 Summary statistics of COI and CytB sequence variation for populations of P. minutus, P. moultoni, P. elongatus, P. acuspes and P. newmani. N = number of sequences analyzed; BP = sequence length in base pairs; H = number of haplotypes; Hd = haplotype diversity; π = nucleotide diversity; TajD = Tajima's D statistic (Tajima, 1989); FLD and FLF = Fu and Li's D* and F* statistics (Fu and Li, 1993); Fs = Fu's S statistic (Fu, 1997). P b 0.05 = ** P b 0.0001 = **. Gene

Species

n

BP

H

Hd

π

TajD

FLD

FLF

Fs

COI

P. P. P. P. P. P. P. P. P. P.

138 51 21 55 24 33 34 18 10 10

473 473 473 473 473 163 163 163 163 163

19 22 6 19 8 7 9 4 2 3

0.69 0.87 0.56 0.82 0.77 0.58 0.76 0.40 0.36 0.51

0.0029 0.0044 0.0018 0.0094 0.0051 0.0065 0.0074 0.0027 0.0022 0.0034

−1.68* −1.55* −1.70* −2.45** −1.61* −0.03 0.01* −0.03* −0.02 −0.02

−3.22** −2.88* −1.24 −3.66** −2.11* −0.08 −0.08 −0.03 0.01 −0.06

−3.15** −2.98* −1.59 −3.83** −2.28* −0.02 −0.02* −0.01 −0.01 −0.06

−12.14** −27.39** −2.43* −5.85* −0.96 0.09 −0.09** 0.24 0.32 0.29

CytB

minutus moultoni elongatus acuspes newmani minutus moultoni elongatus acuspes newmani

Scotia), with 10 private haplotypes radiating from it; an additional 4 common haplotypes were separated from the major haplotype by 2 or 3 base substitutions, with private haplotypes radiating from them (Fig. 5). Compared to the other four species, P. acuspes displayed a more diverse pattern, with two major haplotypes and several private haplotypes separated from the two major haplotypes by a high number

of mutations (i.e. the two Barents Sea haplotypes). Three P. acuspes haplotypes (a, c and i; all from Balsfjord) were too divergent for the set 95% connection limit. Haplotype diversity ranged from 0.56 to 0.87 and nucleotide diversity from 0.0018 to 0.0094. While P. moultoni and P. acuspes had the highest haplotype diversities, P. newmani and P. acuspes displayed the

Fig. 2. Maximum likelihood tree based on COI sequence variation showing relationships among all COI haplotypes from five species of Pseudocalanus. Phylogenetic analysis used the HKY (Hasegawa et al., 1985) model, with gamma correction. Bootstrap values among 1000 sub-replicates are shown at branch points (Felsenstein, 1985). Haplotypes are given as letters corresponding to minimum spanning networks (Fig. 5) and GenBank Accession Nos. KF-991131-KF991202.

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

483

Fig. 3. Maximum likelihood tree based on CytB sequence variation showing relationships among all CytB haplotypes from five species of Pseudocalanus. Phylogenetic analysis used the HKY (Hasegawa et al., 1985) model, with gamma correction. Bootstrap values among 1000 sub-replicates are shown at branch points (Felsenstein, 1985). Haplotypes are given as letters.

highest nucleotide diversity (Table 2). Several tests of evolutionary neutrality were performed for all five species. Fu's Fs was negative for all species and the values were significant for all species except P. newmani. Fu and Li's D and F were also negative for all species and significant for four species, but not for P. elongatus. Tajima's D was significantly negative for

all species (Table 2) and Strobeck's S was high in all species (range 0.995–1.000). The Bayesian skyline analysis, which summarizes instantaneous estimates of effective population size, showed population expansions by approximately one order of magnitude for P. minutus and P. moultoni.

Fig. 4. Maximum likelihood tree based on ITS 1 sequence variation showing relationships among ITS 1 consensus sequences from five species of Pseudocalanus. Phylogenetic analysis used the HKY (Hasegawa et al., 1985) model, with gamma correction. Bootstrap values among 1000 sub-replicates are shown at branch points (Felsenstein, 1985). GenBank Accession Nos. KF991203-KF991207.

484

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

Fig. 5. (color) Map of the North Atlantic Ocean (A); the circles indicate collection locations for samples analyzed in this study; circle colors are reflected in haplotype network diagrams (B–F); empty circles indicate samples that contained no Pseudocalanus. B–F: COI gene genealogies for P. minutus (B), P. moultoni (C), P. elongatus (D), P. newmani (E) and P. acuspes (F). Each node represents one base change; node size corresponds to the number of individuals sharing the same haplotype (n). The haplotype names for each species are given in Fig. 2.

The timing of the expansion events differed between the two species, starting approximately 25,000 and 65,000 years ago for P. minutus and P. moultoni, respectively. The population of P. elongatus displayed no expansion and the results suggest that this population has been stable for the past 15,000 years (Fig. 6). The skyline analysis also showed that the

P. newmani population has been stable for 250,000 years. According to these results, the population of P. acuspes has been stable for several hundred thousand years, but underwent a bottlenecking event and reached its population low point approximately 50,000 years before present (YBP), after which the population has been increasing (Fig. 6).

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

485

3.2. Cytochrome b Cytochrome B sequence data from the five species of Pseudocalanus comprised 105 sequences and resulted in an alignment length of 163 b-p; no insertions or deletions were detected. The alignment contained 55 polymorphic sites, of which 53 were parsimony informative. Twenty-five haplotypes were recovered. The tree-based analysis displayed varying levels of intraspecific phylogenetic structure (Fig. 3), with the highest level of structuring observed for P. minutus. Minimum spanning networks constructed with CytB data mirrored the general pattern seen in the COI networks for P. elongatus, P. minutus and P. moultoni (data not shown). Similarities of CytB and COI networks for P. newmani and P. acuspes could not be ascertained due to fewer recovered haplotypes for CytB compared to COI (3 and 2 for P. newmani and P. acuspes, respectively). Haplotype and nucleotide diversities for CytB agreed with the configuration observed for COI with respect to P. minutus, P. elongatus, P. moultoni and P. newmani. Compared to the COI values, CytB haplotype and nucleotide diversities for P. acuspes were lower. The evolutionary neutrality tests based on CytB sequence data all centered around zero and most were slightly negative (Table 2). Tajima's D was significant for P. elongatus (− 0.03, P b 0.05) and P. moultoni (0.01, P b 0.01). Fu and Li's F (−0.02, P b 0.05) and Fu's Fs (−0.09, P b 0.01) were significant for P. moultoni (Table 2). 3.3. Internal transcribed spacer region I The nuclear marker, ITS-1, was sequenced for 98 individuals resulting in a 454 b-p alignment with five heterozygous sites. There was little intra- or interspecific variation and consensus sequences of the five species of Pseudocalanus were used to construct a phylogeny to assess the evolutionary relationship between the five species established by the mitochondrial markers (Fig. 4). The phylogeny displayed two clades: one comprising P. minutus, P. moultoni and P. elongatus; with P. acuspes and P. newmani constituting the other. 4. Discussion 4.1. Comparisons of Pseudocalanus spp.

Fig. 6. Bayesian skyline plots showing the median values and 95% highest posterior density (HPD) for: P. minutus (A), P. moultoni (B), P. elongatus (C), P. acuspes (D), and P. newmani (E). Only the highest posterior density (area shaded in grey) should be considered for robust inference of change in population size. Time on the x-axis is linear; population size on the y-axis is a logarithmic scale.

4.1.1. Pseudocalanus minutus The lowest levels of genetic structure in this study were found in Pseudocalanus minutus, with all variation explained within populations according to AMOVA analysis. Although all sampling locations except Håkøy revealed private haplotypes, the four most common COI haplotypes comprised 79% of the individuals and were approximately evenly distributed among the Iceland 1, Iceland 2 and the Norwegian Sea samples. The same pattern was evident in the CytB sequence data with three common haplotypes. The only previous study of genetic structure of P. minutus, by Sévigny et al. (1989), was spatially less extensive, and revealed no divergence between populations of P. minutus from Hudson Bay and Nova Scotia based on the glucose phosphate isomerase (GPI) locus. These results indicate that either contemporary gene flow is high across the sampled range for P. minutus, or that the species has experienced a genetically homogenizing event with insufficient time for population structure to reestablish. The tests of evolutionary neutrality showed that Fu's Fs was negative for both mitochondrial genes (CytB, ns), which occurs when excess rare haplotypes are present and is suggestive of either population expansion or genetic hitchhiking (Fu, 1997). Fu and Li's D and F were also negative for both mitochondrial genes. Negative D and F indicate an excess of recently-derived haplotypes, which can occur through population expansion or background selection. Tajima's D was also negative, which indicates population size expansion and/or purifying selection (Tajima, 1989). Strobeck's S was high (0.995–1.000). These results show deviation from neutrality expectations and indicate that the population has undergone either

486

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

selection or expansion. The Bayesian skyline analysis indicated that the population of P. minutus underwent an expansion some 25,000–10,000 YBP. This estimate coincides with the warmer period between the Last Glacial Maximum (27,000–19,000 YBP) and the Late Glacial Maximum 13,000–10,000 YBP). Bucklin and Wiebe (1998) hypothesized that Calanus finmarchicus may have experienced a range reduction of 75% and latitudinal displacement approximately 18,000 YBP. The distributional ranges for C. finmarchicus and P. minutus are similar and it seems possible that, as the ice retreated northward and sea levels rose following the LGM, P. minutus expanded its population accordingly. 4.1.2. Pseudocalanus moultoni Displaying the widest geographical range in this study, P. moultoni was distributed throughout the NW, NC, and NE Atlantic. Still, 97% of the genetic variation (by AMOVA) was explained by withinpopulation comparisons, and the most common haplotype was found in individuals from four locations (Georges Bank, Iceland, Balsfjord and Håkøy). The second most common haplotype was nearly as widely distributed. Compared to P. minutus, there were relatively more private haplotypes from all locations, except Iceland, indicating that P. moultoni has had more time to establish a genetic structure. The results for P. moultoni also displayed deviation from neutrality and even more strongly (compared to P. minutus) indicated that the population has undergone either selection or expansion. Fu's Fs was very negative for COI and also significantly negative for CytB. Fu and Li's D and F and Tajima's D yielded negative values for COI. Similar to the shape of the plot for P. minutus, the Bayesian skyline for P. moultoni displayed a population expansion that started around 65,000 YBP and leveled off approximately 35,000 YBP. These results suggest that P. moultoni began to expand and stabilize approximately at the onset of the glaciation culminating in the LGM. Furthermore, no evidence indicated that the population of P. moultoni was reduced during the LGM. This species has a large distributional range and is possibly oceanic and cosmopolitan; it may thus have the ability to endure a host of environmental conditions. It may be speculated that P. moultoni tracked its habitat during the LGM and was able to maintain a relatively large population. 4.1.3. Pseudocalanus elongatus The AMOVA results indicated strong spatial structuring between the two geographic groups of P. elongatus with 16% of the variance explained among populations. This species displayed clear genetic structuring when considering the two mitochondrial markers. Both COI and CytB showed that all individuals from the East English Channel belonged to one haplotype (N = 13), while only one Iceland individual shared that main haplotype. The remainder of the Iceland individuals displayed private haplotypes. Analysis including the P. elongatus COI sequences deposited by Unal et al. (2006) and Holmborn et al. (2011) revealed that Unal et al.'s main haplotype “H3” and the single sequence deposited by Holmborn et al., were identical to the P. elongatus haplotype “b” in the present study. The remainder of Unal et al.'s sequences, both Black Sea and remaining United Kingdom haplotypes, blended into the star shape pattern shown in Fig. 5, the greatest distance from haplotype “b” being two mutations. The Bayesian analysis indicated that the population has been stable for the past 15,000 years, implying that P. elongatus, like P. minutus and P. moultoni maintained population sizes during the Last Glacial Maximum. Of all the species of Pseudocalanus, P. elongatus exhibits the southern-most distribution and could have sustained its populations closer to the equator. 4.1.4. Pseudocalanus acuspes Variation in COI showed that P. acuspes had the highest amongsample variance of all species in this study. Compared to the other species of Pseudocalanus, P. acuspes exhibited high divergence of haplotypes within a single sample. Additionally, two private haplotypes in the Barents Sea sample were separated from the two most common haplotypes by 9–11 mutations and three Balsfjord haplotypes differed

by N 3%. These results differ from Sévigny et al. (1989), who found no significant difference among samples from Bedford Basin, Resolute Bay and the Baltic Sea. Holmborn et al. (2011) found that mitochondrial diversity was exceptionally low in Baltic Sea P. acuspes, reporting a haplotype diversity of 0.024 and discerning only two haplotypes in 83 individuals taken from the Gulf of Finland in the north to the Arkona Basin in the south. Inclusion of the COI sequence from Holmborn et al. (2011, GenBank accession no. HM770074.1) in the TCS analysis revealed that it was identical to haplotype “e” (Figs. 2 and 5). The populations sampled by Sévigny et al. (1989) and Holmborn et al. (2011) may represent relicts remaining after the last glaciation, and the pattern of low genetic diversity in the Baltic Sea has been observed for several taxa Holmborn et al. (2011). The tests of evolutionary neutrality in this study yielded significant negative values and the COI sequence data for P. acuspes indicated that the population underwent a bottleneck event some 200,000–50,000 YBP, with slight population decrease and subsequent population increase. 4.1.5. Pseudocalanus newmani The COI sequence data displayed relatively high levels of intraspecific sequence divergence for P. newmani, compared to the other species. The network analysis showed private haplotypes for all three sampling locations where P. newmani was found, and there were several theoretical haplotypes not recovered among the sequenced individuals. This finding agrees with Sévigny et al. (1989), who found P. newmani to be the most variable among the Pseudocalanus species. The tests of neutrality yielded negative values, and the Bayesian skyline analysis indicated that the sampled P. newmani population has remained stable in size for the past 250,000 years. 4.2. Possible evolutionary divergence within Pseudocalanus Owing to results from the tree-based analyses, which revealed similarities in intraspecific genetic structuring, interspecific phylogeography and demographic history of the five species of Pseudocalanus analyzed in the present study, it is conceivable that these species constitute two evolutionary clades: a high-dispersal oceanic clade consisting of P. minutus, P. moultoni and P. elongatus and a coastal/shelf clade constituted by P. acuspes and P. newmani. Despite having the greatest geographic distance between sampled populations, P. minutus, P. moultoni and P. elongatus displayed the lowest levels of genetic structuring. This can be explained by either recent migration or continuously high dispersal. Conversely, P. acuspes and P. newmani were only found in locations spatially much closer, but showed relatively higher levels of genetic structuring. All three genetic markers used in this study displayed the same basic evolutionary relationship among the five Pseudocalanus species. Moreover, the phylogenetic trees for COI and ITS 1 formed two clades, with P. minutus, P. moultoni and P. elongatus in one, and P. acuspes and P. newmani in the other. Previous records of the biogeography of Pseudocalanus species, combined with our sampling, suggests that the species within each clade share traits related to dispersal: Pseudocalanus minutus, P. moultoni and P. elongatus are widely distributed and appear to be able to maintain their populations across the North Atlantic Ocean. Conversely, P. acuspes and P. newmani seem to be associated with continental shelves and coastal habitat. This distinction may have been missed previously due to failure of earlier studies to discriminate the species, and a lack of consensus regarding Pseudocalanus habitat preferences. While Frost (1989) designates P. moultoni as a temperate coastal species, McLaren et al. (1989) describes the same species as offshore, comparing it to the more coastal P. newmani. Aarbakke et al. (2011) found high connectivity among P. moultoni in samples from Georges Bank, northern Norway and the Svalbard archipelago, and a possible genetic break between Atlantic and Pacific populations of the same species. Laakmann et al. (2013) recorded P. moultoni in the southern North Sea. Corkett and McLaren

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

(1978) discussed the distribution of Pseudocalanus species with respect to distance offshore, and noted that all the areas of highest abundance of Pseudocalanus species are less than approximately 400 km offshore. However, they cite authors who refer to Pseudocalanus species as “typically oceanic” in the Bering Sea (Motoda and Minoda, 1974) and abundant further offshore in the northern North Atlantic (Østvedt, 1955). Corkett and McLaren (1978) concluded that Pseudocalanus species are primarily coastal, but allowed that the neritic-oceanic gradient is unrefined and permits many exceptions. According to Corkett and McLaren (1978), Pseudocalanus species generally occur in the epipelagic zone, but have been found as deep as 1,350 m (Corkett and McLaren, 1978). Similarities in behavior, and especially patterns of diel vertical migration (DVM), may retain P. acuspes and P. newmani in coastal habitat, while P. minutus, P. moultoni and P. elongatus may be more readily dispersed. Studies of DVM behavior that discriminated Pseudocalanus species have only been conducted for P. newmani, which shows tremendous behavioral plasticity (e.g., Frost and Bollens, 1992; Ohman, 1990). Ohman (1990) showed that P. newmani expressed three types of diel vertical migration within a single fjord, while Frost and Bollens (1992) demonstrated the phenotypic behavioral plasticity of P. newmani in response to potential predation by visual and non-visual predators. 5. Conclusions Sequence analysis of the COI and CytB mitochondrial DNA fragments from five species of Pseudocalanus across the North Atlantic Ocean revealed that there are marked differences among Pseudocalanus species with respect to genetic structure, diversity and demographic history. Three species, P. minutus, P. moultoni and P. elongatus show relatively low genetic structuring and short-term demographic stability, following population increases in P. minutus and P. moultoni after and prior to the Last Glacial Maximum, respectively. Two species, P. acuspes and P. newmani display relatively higher levels of genetic structuring, despite being spatially less extensive compared to the three other species, and were found only in coastal samples of the Northeast and Northwest Atlantic, respectively. Tree-based analysis of COI and ITS-1 variation for the five species distinguished two clades, which may represent an evolutionary divergence based on characteristic distribution, habitat preference, and/or behavior. Acknowledgements We gratefully acknowledge Nancy Copley (Woods Hole Oceanographic Institution, USA) and Elvire Antajan (Institut Francais de Recherché pour L'exploitation de la Mer, France) for sample collection assistance. We thank Leocadio Blanco-Bercial (Dept. Marine Sciences, University of Connecticut, USA) for advice and assistance with the data analysis. This study is a contribution from the Census of Marine Zooplankton (CMarZ, see www.CMarZ.org), an ocean realm field project of the Census of Marine Life.[SS] References Aarbakke, O.N.S., Bucklin, A., Halsband, C., Norrbin, M.F., 2011. Discovery of Pseudocalanus moultoni Frost, 1989 in Northeast Atlantic waters based on mitochondrial COI sequence variation. J. Plankton Res. 33, 1487–1495. Altschul, S.F., Madden, T.L., Schaffer, A.A., et al., 1997. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 25, 3389–3402. Blanco-Bercial, L., Alvarez-Marques, F., Bucklin, A., 2011. Comparative phylogeography and connectivity of sibling species of the marine copepod Clausocalanus (Calanoida). J. Exp. Mar. Biol. Ecol. 404, 108–115. Bucklin, A., 2000. Methods for population genetic analysis of zooplankton. In: Harris, R.P., Wiebe, P., Lenz, J., Skjoldal, H.R., Huntley, M. (Eds.), ICES Zooplankton methodology manual. Academic, London, pp. 533–570. Bucklin, A., Bentley, A.M., Franzen, S.P., 1998. Distribution and relative abundance of Pseudocalanus moultoni and P. newmani (Copepoda: Calanoida) on Georges Bank using molecular identification of sibling species. Mar. Biol. 132, 97–106.

487

Bucklin, A., Guarnieri, M., McGillicuddy, D., Hill, R.S., 2001. Spring evolution of Pseudocalanus spp. abundance on Georges Bank based on molecular discrimination of P. moultoni and P. newmani. Deep Sea Res. II Top. Stud. Oceanogr. 48, 589–608. Bucklin, A., Frost, B.W., Bradford-Grieve, J., Allen, L.D., Copley, N.J., 2003. Molecular systematic and phylogenetic assessment of 34 calanoid copepod species of the Calanidae and Clausocalanidae. Mar. Biol. 142, 333–343. Bucklin, A., Wiebe, P.H., 1998. Low mitochondrial diversity and small effective population sizes of the copepods Calanus finmarchicus and Nannocalanus minor: possible impact of climatic variation during recent glaciation. J. Hered. 89 (5), 383–392. Carter, J.C., 1965. The ecology of the calanoid copepod Pseudocalanus minutus Kroyer in Tessiarsuk, a coastal meromictic lake of northern Labrador. Limnol. Oceanogr. 10 (34), 345–353. Clement, M., Posada, D., Crandall, K.A., 2000. TSC: a computer program to estimate gene genealogies. Mol. Ecol. 9 (10), 1657–1660. Corkett, C.J., McLaren, I.A., 1978. The biology of Pseudocalanus. Adv. Mar. Biol. 15, 1–231. Drummond, A.J., Suchard, M.A., Xie, D., Rambaut, A., 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29 (8), 1969–1973. Drummond, A.J., Rambaut, A., Shapiro, B., Pybus, O.G., 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. Mol. Biol. Evol. 22, 1185–1192. Durbin, A., Hebert, P.D., Cristescu, M.E., 2008. Comparative phylogeography of marine cladocerans. Mar. Biol. 155 (1), 1–10. Excoffier, L., Smouse, P., Quattro, J., 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics 131, 479–491. Excoffier, L., Lischer, H.E.L., 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Res. 10, 564–567. Felsenstein, J., 1985. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39, 783–791. Figueroa, D.F., 2011. Phylogenetic analysis of Ridgewayia (Copepoda: Calanoida) from the Galapagos and of a new species from the Florida keys with a reevaluation of the phylogeny of Calanoida. J. Crust. Biol. 31 (1), 153–165. Folmer, O., Black, M., Hoeh, W., et al., 1994. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotechnol. 3, 294–299. Frost, B.W., 1989. A taxonomy of the marine calanoid copepod genus Pseudocalanus. Can. J. Zool. 67, 525–551. Frost, B.W., Bollens, S.M., 1992. Variability of diel vertical migration in the marine planktonic copepod Pseudocalanus newmani in relation to its predators. Can. J. Zool. 49 (6), 1137–1141. Fu, Y.-X., 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147 (2), 915–925. Fu, Y.X., Li, W.H., 1993. Statistical tests of neutrality of mutations. Genetics 133, 693–709. Goetze, E., 2005. Global population genetic structure and biogeography of the oceanic copepods Eucalanus hyalinus and E. spinifer. Evolution 59 (11), 2378–2398. Goetze, E., 2011. Population differentiation in the open Sea: insights from the pelagic copepod Pleuromamma xiphias. Integr. Comp. Biol. 51 (4), 580–597. Goetze, E., Bradford-Grieve, J., 2005. Genetic and morphological description of b i N Eucalanus spinifer T. Scott, 1894 (Calanoida: Eucalanidae), a circumglobal sister species of the copepod E. hyalinus ss (Claus, 1866). Prog. Oceanogr. 65 (1), 55–87. Grabbert, S., Renz, J., Hirche, H.-J., Bucklin, A., 2010. Species-specific PCR discrimination of species of the calanoid copepod Pseudocalanus, P. acuspes and P. elongatus, in the Baltic and North Seas. Hydrobiologia 652, 289–297. Griffiths, R., Tavare, S., 1994. Simulating probability distributions in the coalescent. Theor. Popul. Biol. 46 (2), 131–159. Halsband, C., Hirche, H.J., 2001. Reproductive cycles of dominant calanoid copepods in the North Sea. Mar. Ecol. Prog. Ser. 209, 219–229. Hasegawa, M., Kishino, H., Yano, T., 1985. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22 (2), 160–174. Holm, 1979. A simple sequentially rejective multiple test procedure. Scand. J. Stat. 6, 65–70. Holmborn, T., Goetze, E., Põllupüü, Põllumãe, A., 2011. Genetic species identification and low genetic diversity in Pseudocalanus acuspes of the Baltic Sea. J. Plankton Res. 33 (3), 507–515. Hudson, R.R., 1990. Gene genealogies and the coalescent process. Oxf. Surv. Evol. Biol. 7 (1), 44. Kimura, M., 1985. The neutral theory of molecular evolution. Cambridge University Press. Kingman, J.F., 1982. On the genealogy of large populations. J. Appl. Probab. 19A, 27–43. Knowles, L.L., 2009. Statistical phylogeography. Annu. Rev. Ecol. Evol. Syst. 40, 593–612. Knowlton, N., Weigt, L.A., 1998. New dates and new rates for divergence across the Isthmus of Panama. Proc. R. Soc. Ser. B Biol. Sci. 265 (1412), 2257–2263. Laakmann, S., Gerdts, G., Erler, R., Knebelsberger, T., Martínez Arbizu, P., Raupach, M.J., 2013. Comparison of molecular species identification for North Sea calanoid copepods (Crustacea) using proteome fingerprints and DNA sequences. Mol. Ecol. Res. 13 (5), 862–876. Lee, C.E., 2007. Global phylogeography of a cryptic copepod species complex and reproductive isolation between genetically proximate “populations”. Evolution 54 (6), 2014–2027. Lessios, H.A., 2008. The great American schism: divergence of marine organisms after the rise of the Central American Isthmus. Annu. Rev. Ecol. Evol. Syst. 39, 63–91. Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25 (11), 1451–1452. Lischka, S., Hagen, W., 2005. Life histories of the copepods Pseudocalanus minutus, P. acuspes (Calanoida) and Oithona similis (Cyclopoida) in the Arctic Kongsfjorden (Svalbard). Pol. Biol. 28 (12), 910–921.

488

O.N.S. Aarbakke et al. / Journal of Experimental Marine Biology and Ecology 461 (2014) 479–488

McLaren, I.A., Laberge, E., Corkett, C.J., Sevigny, J.M., 1989. Life cycles of four species of Pseudocalanus in Nova Scotia. Can. J. Zool. 67 (3), 552–558. Merritt, T.J., Shi, L., Chase, M.C., Rex, M.A., Etter, R.J., Quattro, J.M., 1998. Universal cytochrome b primers facilitate intraspecific studies in molluscan taxa. Mol. Mar. Biol. Biotechnol. 7, 7–11. Motoda, S., Minoda, T., 1974. Plankton of the Bering Sea. In: Hood, D.W., Kelley, E.J. (Eds.), Oceanography of the Bering Sea with emphasis on renewable resources, Institute of Marine Sciences, University of Alaska. Occasional Publication No. 2., pp. 207–242. Nei, M., 1987. Molecular evolutionary genetics. Columbia University Press, New York. Norrbin, M.F., 1991. Gonad maturation as an indication of seasonal cycles for several species of small copepods in the Barents Sea. Polar Res. 10 (2), 421–432. Norrbin, M.F., 1994. Seasonal patterns in gonad maturation, sex ratio and size in some small, high-latitude copepods: implications for overwintering tactics. J. Plankton Res. 16 (2), 115–131. Ohman, M.D., 1990. The demographic benefits of diel vertical migration by zooplankton. Ecol. Monogr. 60 (3), 257–281. Østvedt, O.J., 1955. Zooplankton investigations from weather-ship M in the Norwegian Sea 1948–49. Hvalrådets skrifter 40, 1–93. Posada, D., 2008. jModeltest: phylogenetic model averaging. Mol. Biol. Evol. 25, 1253–1256. Rambaut, A., Drummond, A., 2007. Tracer. Version 1.5, http://tree.bio.ed.ac.uk/software/. Renz, J., Peters, J., Hirche, H.J., 2007. Life cycle of Pseudocalanus acuspes Giesbrecht (Copepoda, Calanoida) in the Central Baltic Sea: II. Reproduction, growth and secondary production. Mar. Biol. 151 (2), 515–527. Renz, J., Mengedoht, D., Hirche, H.J., 2008. Reproduction, growth and secondary production of Pseudocalanus elongatus Boeck (Copepoda, Calanoida) in the southern North Sea. J. Plankton Res. 30 (5), 511–528. Sévigny, J.-M., McLaren, I.A., Frost, B.W., 1989. Discrimination among and variation within species of Pseudocalanus based on the GPI locus. Mar. Biol. 102, 321–327.

Strobeck, C., 1987. Average number of nucleotide differences in a sample from a single subpopulation: a test for population subdivision. Genetics 117, 149–153. Tajima, F., 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 133, 693–709. Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., Kumar, S., 2011. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol. Biol. Evol. 28, 2731–2739. Templeton, A.R., Crandall, K.A., Sing, C.F., 1992. A cladistics analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132, 619–633. Thompson, J.D., Higgins, D.G., Gibson, T.J., 1994. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, positionspecific gap penalties and weight matrix choice. Nucleic Acids Res. 22 (22), 4673–4680. Unal, E., Frost, B.W., Armbrust, V., Kideys, A.E., 2006. Phylogeography of Calanus helgolandicus and the Black Sea copepod Calanus euxinus, with notes on Pseudocalanus elongatus (Copepoda, Calanoida). Deep Sea Res. II Top. Stud. Oceanogr. 53, 1961–1975. Venkatesan, M., Westbrook, C.J., Hauer, M.C., Rasgon, J.L., 2007. Evidence for a population expansion in the West Nile virus vector Culex tarsalis. Mol. Biol. Evol. 24 (5), 1208–1218. Weir, B.S., Cockerham, C.C., 1984. Estimating F-statistics for the analysis of population structure. Evolution 38, 1358–1370. Wright, S., 1951. The general structure of populations. Ann. Eugen. 15, 323–354. Xia, X., Xie, Z., 2001. DAMBE: software package for data analysis in molecular biology and evolution. J. Hered. 92 (4), 371–373. Yamaguchi, A., Ikeda, T., Shiga, N., 1998. Population structure and life cycle of Pseudocalanus minutus and Pseudocalanus newmani (Copepoda: Calanoida) in Toyama Bay, southern Japan Sea. Plankton Biol. Ecol. 45 (2), 183–193.