Monitoring Climate Change Impact on the Genetic Population Structure: The Case of the Fivebeard Rockling (Ciliata Mustela, Linnaeus, 1758) In Its Southern Limit of Distribution

Ciliata mustela is a marine inshore fish which occurs from central Portugal to northeastern Norway. We studied the population structure of this species using cytochrome b gene and the first intron of the nuclear S7 ribosomal protein gene and samples ranging from central Portugal to Gullmars Fjord, Sweden. We tested the following alternative hypotheses: 1) is the Portuguese population of the fivebeard rockling self-sustainable? or 2) is this population dependent on migrants from the north? We found no detectable subdivision among locations. Genetic diversity indices did not change along the study area. The data support persistence during one or several glacial cycles and a rapid expansion about 10 thousand years ago. In Portugal the populations of this species are strongly affected by climatic oscillations with severe reductions in warm years (bellow detection) and recoveries in cold years. We found that the percentage of private haplotypes is consistently lower in Portugal than in other locations. Our results support the hypothesis that the Portuguese population is mainly dependent on migrants from northern locations. We discuss the possibility of using species as C. mustella to monitor short term effects on the genetic structure of populations and their relation with climate change. *Corresponding author: Sara M Francisco, Eco-Ethology Research Unit, ISPA – Instituto Universitário. Rua Jardim do Tabaco 34, 1149-041 Lisboa, Portugal, Tel: +351 218 811 700 (ext. 348); Fax: +351 218 860 954; E-mail: sara_francisco@ispa.pt Received October 04, 2013; Accepted October 25, 2013; Published November 01, 2013 Citation: Robalo JI, Lima CS, Francisco SM, Almada VC, Almada F, et al. (2013) Monitoring Climate Change Impact on the Genetic Population Structure: The Case of the Fivebeard Rockling (Ciliata Mustela, Linnaeus, 1758) In Its Southern Limit of Distribution. J Phylogen Evolution Biol 2: 123. doi: 10.4172/2329-9002.1000123 Copyright: © 2013 Robalo JI, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Introduction
Many phylogeographic studies of west European marine fish focused on the relationships between the Atlantic and the Mediterranean, and on the role of the Mediterranean as a potential glacial refugium [1]. In recent years, however, there has been an increasing effort to study both boreal species and those occurring along west Europe [2]. These studies led to the identification of many potential glacial refugia well to the north of the Mediterranean. In a recent review [3], considered the following refugia: Mediterranean, southwest Europe, waters around southwest Ireland and Southwest Britain, the Hürd Deep (a saline lake that existed near the entrance to the North Sea during glacial times) and areas near north Norway, the Faeroes (a group of islands located between Norway and Iceland), and Iceland.
Patarnello et al. and Hickerson and Cunningham [1,4] have discussed factors that may affect the population's structure of a species and, therefore, their hypothetical glacial refugia. Some of the factors referred were: water temperature needed for adults' survival, reproduction and larval development, habitat in which reproduction occurs, adult dispersal, planktonic larval duration, larval dispersal and behaviour, existence of parental care, etc. The importance of each of these factors can only be evaluated using as many species as possible and exploring how their contrasting traits may have affected their persistence and post glacial recolonization.
In this paper we present data on the population structure of the fivebeard rockling, Ciliata mustela (Linnaeus, 1758) and discuss phylogeographic patterns amongst the locations sampled. C. mustela is a demersal cold water marine fish, distributed in the Northeast Atlantic from Lisbon north to Finnmark, around the British Isles, in the Skagerrak and Kattegat and Iceland [5]. This species was reintegrated in the family Gadidae after a complex history of varying family assignments [6]. It is an inshore species, frequent in tide pools and subtidal habitats when shelter is available. It can reach a maximum length of 25 cm but the most common length is of about 17 cm. The eggs and larvae are planktonic and no parental care is provided. Henriques et al. [7] demonstrated that the abundance of the species is quite sensitive to climatic conditions in Portugal. While in years with cold winters and increased southwards flow they can be sampled with some consistency, in years with warm winters and flow from the south the species ceases to be detectable in underwater fish censuses and tide pool surveys south of Lisbon. In a study on fish fauna changes in Portuguese marine waters, Cabral et al. [8] noted that, while some species of tropical affinities like Diplodus bellottii are being detected in increasing numbers, C. mustela is becoming a very rare fish at least from Tagus (Central Portugal) southwards. Therefore, this species is probably very unstable in the southern part of its distribution likely suffering serious range reductions if water temperatures continue to increase. This fact makes it particularly suitable to evaluate the impact of climate changes in the genetic diversity of a species.
In this study we aim to compare the genetic diversity of C. mustela in Portugal with that in other parts of the species range. The fact that in Portuguese waters the population is very likely reduced in warm years lead us to the following alternative hypotheses: 1) is the Portuguese population of the fivebeard rockling self-sustainable?, or 2) is this population dependent on migrants from the north? In the present paper we use one mitochondrial and one nuclear marker to test these hypotheses with samples ranging from Central Portugal to Norway and Sweden.

Methods Sampling
A total of 90 samples of C. mustela were collected from 5 locations along West Europe, from Portugal to Sweden ( Figure 1): Cascais (near Lisbon, Portugal, LI), Ferrol (Northern Galicia, Spain, FE), Roscoff (France, RO), Helgoland (Germany, HE) and Gullmars fjord (Sweden, GU). Samples from the North Sea (Gullmars fjord and Helgoland) came from beach-seining surveys by the Institute of Marine Research (Norway) and Alfred Wegener Institute (Germany) respectively. The samples from Roscoff were provided to us as fin clips by Station Biologique de Roscoff (France -ASSEMBLE project). Samples from Spain and Portugal were collected by hand-netting in tide-pools. These individuals were kept out of water for a very brief period (less than 30 seconds), a small portion of dorsal fin was clipped and they were released immediately to their pool of origin. Specimen collection complied with the current laws of each country.

DNA analysis
All sequences were edited with CodonCode Aligner v.4.0.1. (CodonCode Corporation, MA, USA) and aligned using Clustal X2 [11] using default settings. ARLEQUIN software package V.3.5 [12] was used to estimate the genetic diversity within each sample, to assess potential population differentiation and to perform neutrality tests. It was also used to perform analysis of molecular variance (AMOVA) and to compute pairwise F ST s. Significance levels of all multiple statistical tests were corrected using the false discovery rate (FDR) approach [13] implemented in QVALUE package of software R [14]. In the case of the S7 intron the analyses were also run in ARLEQUIN, after allowing the program to reconstruct the genotypes present, using the ELB algorithm. For the S7 gene deviations from Hardy-Weinberg Equilibrium were assessed using F is as implemented in ARLEQUIN. The χ2-test developed by Salicru et al. [15] was used to access the significance of the differences of haplotype diversity values among populations. This method takes into account differences in the size of the samples to be compared. Because the sample of Helgoland was very small relative to all others (6 individuals) we also compared genetic diversity among locations by randomly sub sampling, with repetition, as implemented in the software PAST v.2.17 [16]. This program does not implement the estimation of gene diversity and nucleotide diversity indices. However, it estimates the Margalef index which can be also used with genetic data as an haplotype richness index. The Shannon H was also used as another indicator of genetic diversity. In addition, a rarefaction analysis was performed assuming a sample size equal to the smallest sample. In both cases, comparisons were based on the inspection of the confidence intervals generated by the program. For algorithms and computational procedures see the reference manual of PAST (and references therein).
To estimate nucleotide pairwise genetic distances between locations we used distance [17], as implemented in ARLEQUIN, to correct for inherited ancestral polymorphism. Haplotype networks were constructed for both markers using TCS v. 1.18 [18]. If the analyses performed failed to detect differentiation between sampling sites the corresponding sequences were pooled as likely members of a single population. Mismatch analysis [19], Fu's Fs [20] and Tajima's D [21] tests were performed to test for possible departures from neutrality.
The appropriate model of sequence evolution for each fragment was determined using the jMODELTEST program [22,23], under the Akaike Information Criterion (AIC). Past population demography of C. mustela was inferred with the cytb data using both mismatch distributions and the extended Bayesian skyline plot (BSP) model as implemented in BEAST v.1.6 [24] employing the Bayesian MCMC coalescent method, a HKY+I+G model of substitution, and a strict clock. This substitution model was chosen as it is the closest model available to the model selected by jMODELTEST. Also, BEAST documentation advises against the use of more complex models (unless they are explicitly selected) to avoid over-parameterization, slower mixing and longer runs. The adoption of a strict clock follows the recommendations of BEAST authors that indicate that it is preferable to start with a strict clock unless inconsistency among different runs of the programs raised doubts on the validity of this model, and a relaxed clock should be attempted. The Bayesian distribution was generated using results from five independent runs of 150 million MCMC steps obtaining effective samples sizes (ESS) of parameter estimates over 200, with a burn-in of 10%. The sampling frequency was every 1000 chains. The time to most recent common ancestor (t MRCA ) and the median and corresponding credibility intervals of the BSP were depicted using Tracer v.1.4 [25]. Effective population size, exponential growth parameter and age of populations were also estimated by the MCMC approach implemented in LAMARC V. 2.1.3 [26], using 10 runs of 12 short chains of 1,000 steps and five long chains of 50,000 steps, with a burn-in of 10,000. The sampling frequency was every 1000 chains. In order to compute estimates of effective population size, their changes with time, and the age of populations, we used the following divergence rates: 2% for cytb and 0.46% for S7. Mutation rates in fish are widely variable [27] and in the absence of a clock calibration for these markers in C. mustela, we address the uncertainty by assuming the mutation rates followed by Bowen et al. [28] and Bernardi and Lape [29].

Results
The cytochrome b (cytb) data set consisted of an amplicon with 605 nucleotides, with 29 polymorphic sites (90 individuals, Genbank Accession Numbers JX127057-JX127146), representing a total of 25 haplotypes. The amplification of the first intron of the S7 gene yielded a fragment of 619 nucleotides long, with 25 variable sites (88 individuals, Genbank Accession Numbers JX126969-JX127056).
The cytb haplotype network ( Figure 2) showed a relatively shallow star-like network, with a central haplotype (identified by TCS as the ancestral one) representing approximately half of the sampled fish (46 of 90 specimens). All other haplotypes are connected to this one, the majority being only one mutation away. The most derived haplotype is only three mutations apart. There was no clear geographic structure in haplotype distribution, as all derived clades except one are shared between the Atlantic (Cascais, Ferrol and Roscoff) and the North Sea (Gullmars fjord and Helgoland). Fish from Lisbon displayed the smallest proportion of private haplotypes (even after rarefaction) ( Table  1). The S7 network is more complex than that of cytb which was to be expected given the fourfold effective population size of nuclear genes, when compared to mitochondrial markers. Five alleles, all closely related (with a central one identified by TCS as the ancestral), represent the vast majority of the sequences, without geographic structure in haplotype distribution. The more derived allele is only four mutations apart.
For the cytb fragment and S7 intron, the genetic diversity indices for each population are summarized in Table 1. The diversity indices were moderate and did not vary significantly among locations. For the cytb data, haplotype diversity values were not significantly different among populations according to the χ 2 test developed by [15] (χ 2 =10.253, p=0.0683). For both cytb and S7, comparisons of the diversity indices based on resampling with repetition (performed with PAST) did not yield significant results as all confidence intervals showed some overlap (results not shown). No pair of confidence intervals was separated according to the two methods (bootstrap and permutations). Rarefaction analysis did not change this situation even when considering the individual rarefaction as low as the smaller sample size (6 individuals for cytb and 5 sequences for S7).
The pairwise F ST s and corrected average pairwise differences yielded only non significant differences between sample pairs, for cytb (results not shown). For the S7 gene the only population pairs  Figure 1). In the case where haplotypes are shared among regions, shading is proportional to the frequency of the haplotype in each region. The area of the circles is proportional to the frequency of each haplotype. In both networks the haplotype considered the ancestral by TCS is the central one.  =0.1560, p=0.05181).

A Mitochondrial cytochrome b gene
For the pooled data, mismatch distributions for cytb and S7 were compatible with both the sudden demographic and the spatial expansion models for the European population of C. mustela. For the ctyb gene the spatial expansion model had, however, a slightly better fit (Lower SSD, Table 2). Visual inspection of the mismatch distributions (data not shown) indicated that, for the cytb, most pairs of sequences showed 0 differences, while 1 difference prevailed for S7. Figure 3 shows the Bayesian Skyline Plot generated by BEAST for the cytb dataset of C. mustela. The plot revealed that fivebeard rockling, as a whole, persisted for more than 70ky at low number, and experienced a rapid population growth in the last 10 thousand years (after the end of the Last Glacial Maximum, LGM), reaching an Effective Population Size (N ef ) of 700 thousand fish in the present day. The Times to the Most Recent Common Ancestor (t MRCA ) estimated based on the cytb and the S7 of C. mustela are shown in Table 2 and are of 19.181 ky for cytb and 129.960 ky for S7.
Demographic results estimated with Lamarc for the cytb and S7 of C. mustela are also shown in Table 2 and confirmed the BSP in revealing a long persistence for at least 227 thousand years. Given the uncertainty associated with the divergence rate adopted and the large errors involved in estimations the values of N ef obtained with the two methods are not discordant.

Discussion
The present data did not reveal any subdivision between Atlantic and North Sea samples. We are aware that our sampling in Helgoland was very limited. Nevertheless the population of Gullmars Fjord, which is located well to the north of Helgoland, had a much larger number of individuals and its close affinity with Atlantic sites was, in our view, well documented.
The genetic pattern found for C. mustela reveals similar values for the genetic diversity indices across the sampled area and shallow star-like haplotype networks, both for cytb and S7 ( Figure 2 and Table  1). This pattern is quite different from that postulated to explain the low diversity at the northern limit of a species. Traditionally, this low diversity is presented as a consequence of successive post-glacial foundations of new populations in an advancing front ( [30]; see [31] for a likely example with Pomatoschistus microps). Interestingly, several studies of temperate fish detected no variation of genetic diversity with latitude (e.g. Sprattus sprattus [32], Taurulus bubalis [2]).
In what concerns population history of fivebeard rockling our results suggest coalescence extending back to the middle of the last glaciation or more, with rapid expansion soon after it ended ( Figure  3 and Table 2). Its ability to live subtidally could facilitate survival in the vicinity of the ice as intertidal temperatures are considerably lower when compared to deeper waters. If, after surviving glaciations in a broader area, cold tolerant species were the first to recolonize northern locations, they likely suffered fewer bottlenecks, keeping substantial diversity and gene flow across their entire range.  Table 2: Demographic parameters of Ciliata mustela based on cytb and S7. Mismatch distributions: τ (time in generations; upper and lower bounds of 95% CI in parenthesis), t (time in years), θ 0 (theta0), N 0 (female effective population size before the expansion), θ 1 (theta1), N 1 (female effective population size after the expansion), SSD (sum of square deviation) and Hri (Harpending's Raggedness index). Estimates of population parameters with LAMARC: θ (theta), N f (female effective population size), G (growth rate), and N 1% (age of population assessed as the age at which N f drops below 1%). Estimates of t MRCA (time to most recent common ancestor) were assessed with BEAST.

Mismatch distribution
The occurrence of C. mustela in cold years in Portugal and its virtual disappearance in warm years [7] strongly suggest that the population is not self-sustained in this area and many recruits result from larval transport from the north. Indeed, the smaller number of private haplotypes in Portugal (cytb, 17%, Table 1), when compared to the remaining populations, provides evidence that supports this interpretation. Nevertheless, genetic diversity shows similar values from Lisbon up to Gullmars Fjord, our location more to the North of the species range. Additionally, the existence of private haplotypes in Lisbon is consistent with a scenario where some survivors of an old local population are mixing with a much greater contingent of northern migrants brought to the South by larval transport. Taken together, the present results support our second hypothesis, with the Portuguese population of the fivebeard rockling being mainly dependent on migrants from northern locations.
Caution should be taken while interpreting the present findings. First, the area of the British Isles is unsampled and there is evidence pointing to an important colonization route of the North Sea via Britain [33,34]. Second, if our divergence rates prove to be unrealistic the interpretations must be re-evaluated. We had no hard calibration points for this fish group and used values to the lower end of the published values for fish. We gave preference to the values of cytb because nuclear fragments, even the relatively fast evolving introns, have evolutionary rates that are typically much lower than those of mitochondrial DNA [35]. Third, in order to distinguish a consistent trend from a sporadic change from a year to the next, this analysis must be repeated in futures studies. Finally, sample sizes must be increased. C. mustela must be studied in more depth as it may prove to be a very interesting species to monitor the population variations between Portugal and the rest of Europe in relation to water temperature and other ecological and oceanographic conditions related to climate change.