Genetic and phenotypic differentiation of lumpfish (Cyclopterus lumpus) across the North Atlantic: implications for conservation and aquaculture

Aims

There were three aims of this study. First, given their limited larval dispersal and evidence for homing, we hypothesised that lumpfish might display genetic isolation by distance, with populations closer together being more genetically similar than those further apart (Rousset, 1997). By sampling across the whole range, we aimed to estimate the level of genetic and phenotypic differentiation, and determine the patterns of gene flow across the species’ range, at spatial scales relevant to management. Secondly, we wanted to know to what extent lumpfish translocations could pose a potential genetic risk to local populations. For this, we examined if Icelandic and Norwegian lumpfish, the two most common sources of lumpfish in aquaculture, were genetically distinct from lumpfish populations present in other salmon farming areas (i.e., Ireland, Canada, Scotland and the Faroe Islands) with which they might interbreed. As lumpfish aquaculture is very recent, our estimates of genetic structure of wild lumpfish populations might serve as genetic baselines against which the impact of farm escapees might be gauged, as has been done for Atlantic salmon (Gilbey et al., 2017). Finally, as some lumpfish populations may be endangered, we provide estimates of effective population size, and test for the existence of genetic bottlenecks to better understand their conservation status.

Collection of samples

Fin tissue was obtained from 410 lumpfish originating from 15 sites across the species’ range (Table 1) and were stored in 96% ethanol at −20 °C until analysis. Sites located within an 80 km radius (the estimated maximum range of dispersal Kennedy et al., 2015) were pooled together to minimise the risk of spatial pseudo-replication. Samples were pooled from the Faroe Islands (Klasvík and Kollafjørður, c. 20 km), Denmark (Køge Bay and Mosede Havn, c. 13 km) and Sweden (survey hauls from Bornhölm to Öland, and from Gotland to Gotska Sandön). Pooled groups were named after the site contributing the largest number of samples. Biometric data on length (mm) and weight (g) were available for eight of the 15 sites (Table 1). Swansea University, College of Science Ethics Review Committee, provided full approval to Benjamin Whittaker for the genetic analysis of tissue samples already collected by other researchers (STU_BIOL_90920_181018132845_1).

DNA extraction and amplification

DNA was extracted using the Nexttec Isolation kit (NextTec, Wellingborough, UK) following the manufacturer’s protocol. The concentration of extracted DNA was quantified using a Nanodrop 2000 (Thermo Fisher Scientific Inc., Waltam, MA, USA) and diluted with DNA free water to 50 ng/µl where necessary. A 2 µl of sample DNA was used for amplification using a QIAGEN Multiplex PCR kit (QIAGEN, Manchester, UK) in a total reaction volume of 9 µl. Ten lumpfish specific microsatellite loci (Clu29, Clu34, Clu36, Clu45 and Clu12, Clu26, Clu33, Clu37, Clu40, Clu44 (Skirnisdottir et al., 2013) were genotyped in two separate multiplex reactions (Table S1). Amplification consisted of a single initial activation step at 95 °C for 15 min followed by eight cycles of touchdown PCR denaturation at 94 °C for 30 s, annealing from 64 °C or 60 °C to 56 °C in descending two-cycle steps of 2 °C and an extension at 72 °C for 90 s, 24 additional cycles with an annealing temperature of 56 °C and a single final extension at 60 °C for 30 min. An Applied Biosystems ABI3130xl Genetic Analyser (Applied Biosystems, UK) was used to resolve the fragments using GeneScan 500-LIZ(-250) as a size standard. Fragment length was established using GeneMapper v5.0 (Applied Biosystems, Foster City, CA, USA). Genotyping consistency was validated by repeating PCR, fragment analysis and scoring for 10% of samples.

Estimates of genetic diversity

We used Microchecker v2.2.3 (Van Oosterhout et al., 2004) to identify null alleles, allele dropout and stutter peaks, and Bayescan v2.1 (Foll & Gaggiotti, 2008) to test for loci neutrality. GENEPOP v4.2 (Rousset, 2008) was used to test for linkage disequilibrium, deviations from Hardy-Weinberg equilibrium, and to calculate allelic frequencies across populations. GeneAlEx v6.502 (Peakall & Smouse, 2012) was used to assess the number of alleles (N_A), effective alleles (N_E), private alleles (N_PA), expected (H_E) and observed heterozygosity (H_O), and to carry out a Mantel test of genetic isolation by distance.

Genetic differentiation, genetic structure and patterns of migration

We estimated the extent of genetic differentiation by calculating global and pairwise F_ST values between populations by 1,000 permutations with a Bonferroni correction (P <  0.00022) using Arlequin v3.5.2.2 (Excoffier & Lischer, 2010). We then conducted an Analysis of Molecular Variance (AMOVA) on pooled groups (separated more than 80 km) to partition genetic variation at three hierarchical levels (among populations, within populations, and among individuals), and estimated significance values using 1,000 permutations. To assess population structuring, a Bayesian cluster analysis was conducted in STRUCTURE v2.3.4 (Falush, Stephens & Pritchard, 2007; Hubisz et al., 2009; Pritchard, Stephens & Donnelly, 2000) to estimate the most likely number of genetic clusters (K) informed by individual genotypes. Admixture models with K values ranging from 2 to 15 were run using twenty iterations, a burn-in length of 10,000 and 50,000 Markov Chain Monte Carlo repeat simulations to estimate the likelihood of each K value. Results were fed into STRUCTURESELECTOR (Li & Liu, 2018) to identify the most likely number of clusters present based on the median of means (MedMeaK), maximum of means (MaxMeaK), median of medians (MedMedK) and maximum of medians (MaxMedK) criteria (Puechmaille, 2016). In addition, we used a Bayesian cluster analysis implemented in TESS v2.3.1 (Chen et al., 2007), which explicitly includes spatial information, to identify genetic discontinuities in our study area and estimates admixture proportions without assuming any predefined populations. For this, spatial coordinates were randomly generated for each sample to fall within one standard deviation of the maximum and minimum latitude and longitude of each site (Durand, Chen & François, 2008). Admixture models were run with 50,000 total sweeps, 10,000 burn-in sweeps, and 200 runs per K_max ranging from 2 to 15. The average Deviance Information Criterion (DIC) of the lowest 10 DIC values was calculated for each K_max to assess the most likely number of clusters. Runs within 10% of the lowest (DIC) for a given K_max were used for analysis. CLUMMP v1.1.2 (Jakobsson & Rosenberg, 2007) was used to average variation between repeated iterations for the most likely K values, and the resulting output was visualised using DISTRUCT v1.1.1 (Rosenberg, 2004). A neighbour joining tree was constructed with Populations v1.2.32 (Langella, 2002) to further analyse population differentiation using Nei’s standard genetic distance with 1,000 bootstraps per locus, and the resulting tree was visualised using TreeView (Page, 2003). Patterns of gene flow were estimated using div-Migrate, which defines a pool of simulated migrants based on the allele frequencies of each population pair and estimates asymmetric gene flow by measuring the genetic differentiation between each population and the migrants pool (Sundqvist et al., 2016). We used the Nm statistic (Alcala, Goudet & Vuilleumier, 2014) that combines elements of G_ST and D to calculate directional relative migration and assessed whether gene flow was significantly asymmetric between populations by running 5,000 bootstrap simulations.

Effective population size and evidence of genetic bottlenecks

Estimates of effective population size (N_e) for sites containing at least 19 individuals were calculated using the Linkage Disequilibrium Model (LDM) with a critical value of 0.02 in NeEstimator v2.1 (Do et al., 2014). Evidence of genetic bottlenecks was evaluated with Bottleneck v2.1 (Cornuet & Luikart, 1996) using 1,000 replicates under the Two-Phase (TPM) and the Stepwise (SMM) Mutation Models to assess heterozygosity deficiency.

Phenotypic variation

Variation in the Length–weight relationship between regions (West Atlantic, n = 30; East Atlantic, n = 65; English Channel, n = 60; Baltic Sea, n = 40), was examined by regression analysis on log-transformed data (R Development Core Team, 2013). We calculated relative weight (Wr) as the ratio of the observed weight divided by the predicted weight (from the regression obtained above) to obtain an index of body condition that is more appropriate for fish like lumpfish that have an unusual body shape (Nahdi, Garcia de Leaniz & King, 2016). The most plausible number of age classes represented in the samples, and the mean size at age (Macdonald & Pitcher, 1979) were calculated through mixture analysis of length-frequency data using PAST v3.17 (Hammer, Harper & Ryan, 2001). The Von Bertalanffy growth equation (Kirkwood, 1983) was fitted to estimate growth parameters in each region.

Population genetic diversity

All microsatellite loci were polymorphic. The mean number of alleles (N_A) ranged from 4.5 (Ro) to 6.8 (Kl, VB), mean expected heterozygosity (H_E) ranged from 0.592 (Öl) to 0.700 (Kl), and mean F_IS varied from −0.056 (Ro) to 0.110 (Öl) across all loci (Table 1). Initial analysis suggested that null alleles might be present at multiple loci (Clu34, Clu36, Clu12, Clu33, Clu37 and Clu40, Table S2 ). However, repeatedly removing each locus in turn showed little variation in F_ST values (Table S3–S8), and therefore all markers were retained for further analyses. No evidence of departures from neutrality or linkage disequilibrium was found after Bonferroni corrections for multiple tests (Rice, 1989). Deviations from Hardy-Weinberg equilibrium (HWE) were detected at five of the 15 sites (Table 1), but these involved only 12% of loci after Bonferroni correction (Table S9). The mean number of private alleles (N_PA) was relatively low, ranging from 0.00 to 0.40, with sites in the West Atlantic (FB = 0.30, WB = 0.40) and Baltic Sea (GS = 0.30) showing the highest values.

Population structure and gene flow

Global F_ST was 0.095 (P <  0.001) indicating a moderate but significant degree of genetic differentiation for a marine fish (Hartl, Clark & Clark, 1997; Ward, 2000; Hutchinson, Carvalho & Rogers, 2001). Results of AMOVA indicated that 83.5% of molecular variation was due to variation within individuals, 7% amongst individuals within populations, and 9.5% amongst populations (Table S10). Pairwise F_ST showed a significant level of genetic differentiation across most populations (Table 2), but populations closer together were genetically more similar after Bonferroni correction. On the basis of F_ST values, the strongest differentiation was found between West Atlantic and Baltic Sea populations. Results of a Mantel test support the existence of a significant, albeit weak, isolation by distance (R² = 0.1229, P = 0.01).

The most likely number of genetically distinct groups (K) ranged from K = 5 (MedMedK, MedMeaK) to K = 6 (MaxMedK, MaxMeaK) using STRUCTURESELECTOR (Fig. S1). Spatial cluster analysis using TESS suggested a K_max = 10 (Fig. S1), though only six of these genetic groups showed substantial representation, and four groups contributed only 3.3% to the genetic background. Distinct clusters were detected in the West Atlantic and Baltic Sea by both STRUCTURE and TESS, with a greater level of admixture across the East Atlantic (Figs. 1A–1C). Results were consistent in attributing a genetically unique pattern to the Mid Atlantic, English Channel clusters and a Norwegian site at Averøy. A neighbour joining tree (Fig. 2) showed similar patterns to that of the structuring analyses, highlighting the separation between the West Atlantic and Baltic Sea populations, and the higher degree of admixture within the East Atlantic group.

The effective number of migrants (N_m) ranged from 1.00 between sites in the English Channel to 0.03 between sites in the West Atlantic and Baltic Sea. The exchange of migrants was much higher within genetic clusters than among clusters (Table S11), with the highest levels of gene flow found within the East Atlantic and within the English Channel (Fig. 3). The only evidence of moderate asymmetric gene flow was from Norway towards the Faroe Islands (N_m = 0.507, P < 0.05).

Effective population size and evidence of genetic bottlenecks

Estimates of effective population size (N_e) based on a Linkage Disequilibrium Model (LDM) varied from 19 (Norway) to 70,148 (Denmark; Table S12). Sites with low N_e values (<75) were found across Iceland, Faroe Islands and Norway (Fig. 3). A significant deficiency of heterozygotes was identified in Ireland and Scotland using the Single Mutation Model (SMM) in Bottleneck (Wilcoxon signed-rank test, P = 0.0033 after Bonferroni correction), suggesting that these populations could have undergone a recent genetic bottleneck (Table  S13), but this was not detected by the Two-Phase Model of Mutation (TPM).

Phenotypic variation

The relationship between length and weight differed significantly between regions (F_4,192 = 917.2, P < 0.001; Fig. 4). Lumpfish in the Baltic Sea were heavier relative to their size than lumpfish in the East Atlantic and the English Channel (pairwise comparisons: Baltic—East Atlantic, estimate = −0.090 ± 0.036, t =  − 2.530, P = 0.012; Baltic—English Channel, estimate = −0.145 ±0.046, t =  − 3.171, P = 0.002), but were similar to those in the West Atlantic (pairwise comparison Baltic—West Atlantic, estimate = −0.094 ± 0.050, t =  − 1.891, P = 0.060). The relative weight of lumpfish differed between regions (F_3,191 = 2.841, P = 0.039) and was highest in the Baltic Sea and the West Atlantic, and lowest in the East Atlantic and the English Channel.

Mixture analysis identified multiple plausible age classes present amongst lumpfish sampled in the Baltic Sea (seven age classes), East Atlantic (four age classes) and English Channel (three age classes), but only a single plausible age class in the West Atlantic. Based on the parameters of the Von Bertalanffy Growth equation, the maximum age was estimated to be 6.0 years for Baltic populations, 5.7 yrs for populations in the East Atlantic and 7.5 yrs for southern populations spawning in the English Channel. Fitted growth equations differed significantly between regions (Table 3), with lumpfish in the Baltic Sea showing the slowest growth and those in the English Channel showing the fastest.

Conclusions and management implications

By 2020 c. 50 million lumpfish will be required by the salmon farming industry (Powell et al., 2018a; Treasurer, 2018) and most of these will come from the stripping of wild broodstock (Wittwer & Treasurer, 2018) caught in Iceland and Norway, and then shipped as eggs or larvae to salmon farms elsewhere. Information on lumpfish escapees is lacking but corkwing wrasse (Symphodus melops) deployed as cleaner fish in Norway have recently been found to escape and hybridise with local populations (Faust et al., 2018), and the same could happen with lumpfish. Efforts should thus be made to reduce the risk of lumpfish escaping from fish farms and interbreeding with local populations, as high propagule pressure associated with open-net pens is the single most important factor determining the impact of escapees (Consuegra et al., 2011).

Our study suggests that lumpfish translocations should be restricted within genetically homogenous groups to reduce the risk of genetic introgression between native and non-native populations. In this sense, lumpfish from some areas of Norway, and particularly from Iceland, may be ill-suited for deployment in Ireland, Scotland and the Faroe Islands, and vice-versa. Ultimately, closing the breeding cycle of the species in captivity, and producing sterile lumpfish for deployment in salmon farms, must be a research priority for both the conservation of the species and the cleaner fish industry (Powell et al., 2018a), as this will lessen dependency on wild broodstock and reduce the risk of genetic introgression.