Speciation reversal in European whitefish (Coregonus lavaretus (L.)) caused by competitor invasion.

Invasion of exotic species has caused the loss of biodiversity and imparts evolutionary and ecological changes in the introduced systems. In northern Fennoscandia, European whitefish (Coregonus lavaretus (L.)) is a highly polymorphic species displaying adaptive radiations into partially reproductively isolated and thus genetically differentiated sympatric morphs utilizing the planktivorous and benthivorous food niche in many lakes. In 1993, Lake Skrukkebukta was invaded by vendace (Coregonus albula (L.)) which is a zooplanktivorous specialist. The vendace displaced the densely rakered whitefish from its preferred pelagic niche to the benthic habitat harbouring the large sparsely rakered whitefish. In this study, we investigate the potential influence of the vendace invasion on the breakdown of reproductive isolation between the two whitefish morphs. We inferred the genotypic and phenotypic differentiation between the two morphs collected at the arrival (1993) and 15 years after (2008) the vendace invasion using 16 microsatellite loci and gill raker numbers, the most distinctive adaptive phenotypic trait between them. The comparison of gill raker number distributions revealed two modes growing closer over 15 years following the invasion. Bayesian analyses of genotypes revealed that the two genetically distinct whitefish morphs that existed in 1993 had collapsed into a single population in 2008. The decline in association between the gill raker numbers and admixture values over 15 years corroborates the findings from the Bayesian analysis. Our study thus suggests an apparent decrease of reproductive isolation in a morph-pair of European whitefish within 15 years (≃ 3 generations) following the invasion of a superior trophic competitor (vendace) in a subarctic lake, reflecting a situation of "speciation in reverse".

Introduction

Reproductive isolation is a process that impedes the exchange of genes between the members of distinct populations [1]. In ecological speciation, reproductive isolation evolves as a result of adaptive divergence between different sub-populations in response to divergent natural selection within a population [1], [2]. An important promoter of divergent natural selection in ecological speciation is the availability of various ecological opportunities. This provides alternative ecological niches that may favour different behavioral and morphological adaptations which eventually may lead to genetic divergence and reproductive isolation of different morph types [3]–[5]. The reproductive isolation may include both pre-zygotic and post-zygotic isolation mechanisms (Table 1). Pre-zygotic isolation mechanisms mainly involve spatial and temporal isolation (e.g. different spawning sites and time) and sexual selection [6]. The temporal and spatial pre-zygotic isolation mechanism may be important as it restricts gene flow among populations [7]. Sexual selection also represents a potential driver of reproductive isolation (mate choice) [6], [8]–[11]. In fish, premating reproductive isolation mechanisms are usually weak and may permit occasional gene flow among populations [12], [13]. However, the strength of pre-zygotic reproductive isolation mechanisms depends heavily on the environmental and ecological conditions which the isolation barriers are built upon [14]. For example in cichlids and Alpine whitefish, eutrophication has changed the reproductive strategies of species [11], [15], whereas the American signal crayfish has altered the reproductive behavior of males of three-spined sticklebacks [16], resulting in “speciation reversal” [17]. Post-zygotic isolation mechanisms may include extrinsic and intrinsic post-zygotic mechansims. Extrinsic mechansims may involve ecological inviability, whereas intrinsic mechanisms include hybrid inviability and hybrid sterility [1]. Such mechanisms have been documented in three-spined sticklebacks (Gasterosteus aculeatus), lake whitefish (Coregonus clupeaformis), and Alpine whitefish (Coregonus spp.) [18]–[20].

Table 1

Putative reproductive isolation mechansims in European whitefish (Coregonus lavaretus) and other post-glacial fish.

Classification of putative reproductiveisolation mechanisms	Components	References
Pre-zygotic isolation mechansims
1. Temporal isolation	Different spawning time	[22], [37]
2. Habitat isolation	Different spawning ground or microhabitat segregation	[17], [22], [44], [45]
3. Mate choice	Size based assortative mating	[46]–[48]
Post-zygotic isolation mechansims
Extrinsic
Ecological inviability and incompatibilities	Lack of ecological niches for intermediates, reduced feeding efficiency dueto non-optimal feeding apparatus, asynchronous hatching time, emergence	[19], [20]
Intrinsic
Hybrid inviability	Lethal or partial developmental abnormalities or embryonic mortalitiesin hybrid offsprings due to the genomic incompatibilities	[18], [49]–[51]

Sympatric morph-pairs with partial reproductive isolation frequently occur among post-glacial freshwater fish species such as whitefish Coregonus spp. and Arctic charr Salvelinus alpinus (L.) [18], [21], [22]. Due to the young geological age, these systems are considered to constitute an early phase of the speciation process [22], [23]. However, anthropogenic impacts like climate change, habitat destruction, and introductions of non-native species may cause instability of native ecological conditions that the pre-mating isolation mechanisms between morphs or species-pairs are dependent upon [15], [17], [24]–[27]. This may result in breakdown of weak pre-zygotic reproductive isolation, inducing hybridization between closely related species or morph-pairs [4], [17]. Increased hybridization between sympatric species following environmental disturbances has been reported in both plants and animals [26], [28], [29].

In northern Fennoscandian lakes, a morph-pair of whitefish referred to as the densely–rakered (DR) and large sparsely–rakered (LSR) whitefish, commonly occur in sympatry [30]–[32]. The DR whitefish is characterised by numerous long, thin and densely packed gill rakers, whereas the LSR whitefish has fewer, shorter, and more sparsely packed gill rakers and usually a larger body size. The DR whitefish generally occupies the pelagic zone feeding predominantly on zooplankton, whereas the LSR whitefish predominantly feeds on zoobenthos in littoral habitats [30]. This use of divergent niches has repeatedly been observed to promote phenotypic and partial reproductive isolation, thus resulting in genotypic differentiation between the two sympatric whitefish morphs across lakes in northern Fennoscandia [21], [32]–[34].

Vendace Coregonus albula (L.) is a pelagic, highly specialized zooplankton feeder known to be competitively superior to whitefish in zooplankton foraging [35]–[38]. Vendace were intentionally and repeatedly translocated into the tributaries of Lake Inari in northern Finland in 1964–66 [39], [40] and by the 1980’s, a large population of vendace was established in the lake. From Lake Inari, vendace migrated downstream into the Pasvik watercourse where the species was observed for the first time in 1989 [38]. In 1993, a few vendace were recorded in Lake Skrukkebukta in the lower part of the watercourse [38]. Since then, the vendace population in L. Skrukkebukta has been steadily growing [41]. Long-term studies of fish and zooplankton following the vendace invasion have demonstrated a strong impact of the invader on the size-distribution and abundance of zooplankton [42], which in turn have resulted in a diet shift of the DR whitefish from zooplankton to zoobenthos combined with a competitive relegation from the pelagic to the littoral habitat [35], [38], [43].

Despite the fact that sympatric DR and LSR whitefish are partially reproductively isolated and genetically differentiated [21], [34], the competitive displacement of the DR whitefish from the pelagic into the littoral habitat, which is inhabited by the LSR whitefish, may have enhanced the probability of weakening pre-zygotic isolation between DR and LSR whitefish. Moreover, the overall abundance of the DR whitefish has in some lakes been reduced to around 10% of its pre-invasion value [35]. This raises the hypothesis that the vendace invasion in the watercourse may have triggered breakdown of reproductive isolation between the DR and LSR whitefish. Here, we aimed at testing this hypothesis by contrasting multilocus genotype data from the DR and LSR whitefish in Lake Skrukkebukta at the arrival (1993) and 15 years after (2008) the invasion of vendace. We predicted that while the DR and LSR whitefish would be phenotypically and genetically differentiated prior to the vendace invasion, a breakdown of reproductive isolation and species collapse may have occurred following the invasion (but see [15], [17]).

Methods

Study Area and Sample Collection

This study was carried out in the dimictic and oligotrophic Lake Skrukkebukta situated in the Pasvik watercourse, Finnmark county, northern Norway (Figure 1). Lake Skrukkebukta (69° 33′ N, 30°7′ E; 21 m a.s.l.) has an area of 6.6 km 2, and a maximum depth of 37 m. In addition to whitefish, the lake harbours perch (Perca fluviatilis L.), pike (Esox lucius L.), brown trout (Salmo trutta L.), burbot (Lota lota (L.)) and nine-spined stickleback (Pungitius pungitius (L.)). The dataset used for this study consists of gill samples from two sampling years; at the arrival of vendace (1993) and after its invasion and establishment as a superior pelagic feeder (2008). Fish were sampled in the pelagic, littoral and profundal zones using 40 m long gill nets with eight sections of 5 m each, with different mesh sizes of 10, 12.5, 15, 18.5, 22, 26, 35 and 45 mm (knot to knot). In the pelagic zone, 6 m deep floating nets were used, whereas in the littoral and profundal zones, 1.5 m deep bottom nets were employed. The whitefish were classified in the field into either DR or LSR whitefish based on the overall gill raker and body morphology [30]–[32], [38]. In a study using a long term dataset, Siwertsson et al [52] revealed that the gill raker numbers is a stable phenotypical trait in European whitefish (Coregonus lavaretus) populations in subarctic lakes and constitutes a reliable means to classify morphs of European whitefish. The gill arches were removed and preserved in 96% ethanol for further analysis. The gill rakers on the first left branchial arch were counted under dissection microscope. For the genetic analyses we used 93 fish from year 1993 and 96 fish from year 2008.

Figure 1

Map of Pasvik watercourse, Finnmark county, northern Norway.

Black arrow indicates the location of Lake Skrukkebukta in Pasvik water course.

Fish were euthanized by means of cerebral concussion prior to sample collection. A fishing permission is required from the fishing right owner, which on Government land in Finnmark county is the County Governor of Finnmark with legal authority through LOV 1992-05-15 nr 47, §13. Accordingly, we obtained permissions for the gill net fishing in Lake Skrukkebukta from the County Governor. No ethical permission is required from the Norwegian Animal Research Authority for collection with gill nets and the associated sacrifice of fish (FOR 1996-01-15 nr 23, the Norwegian Ministry of Agriculture and Food).

Analysis of Gill Raker Number

For the present dataset, the number of gill rakers represents the only phenotypic trait that could be associated with individual genotypes. The number of gill rakers is an important trait involved in trophic utilization via the feeding efficiency [30], [32], [53], [54] and in coregonids, it has been shown to be heritable [55], [56] and influenced by diversifying selection [23], [57]. Moreover, since sympatric divergence in this system is promoted by ecological opportunity [54] and number of gill rakers contribute to the local adaptation [23], this trait could directly or indirectly impact on the extent of reproductive isolation. Thus, Kahilainen et al [54], demonstrated that increasing number of gill rakers facilitates the foraging of small prey, particularly zooplankton, and that fish with a high number of gill rakers generally are associated with a pelagic, zooplanktivore niche, whereas fish with lower gill raker numbers are associated with benthic habitat and diet utilization. Densely rakered gills are more likely to be clogged by bottom sediments, whereas sparsely rakered gills may not be able to filter small zooplankton efficiently [54]. In a field experiment demonstrating extrinsic post-zygotic isolation mechanism, three-spined stickleback hybrids showed slower growth rate as compare to the parental ones as a consequence of inheritance of intermediate phenotypic traits (including gill rakers) which likely resulted in lower fitness [19]. We used this trait to examine the potential breakdown of morphological differentiation by exploring the modality of the gill raker number distribution over the two sampling years (1993 and 2008). Kernel density estimation was used to produce smoothed gill raker distributions. A kernel smoother estimates the underlying probability density function for the given dataset and identifies the number of distributions. This was achieved in R using the data smoothing function density [58]. Secondly, we employed a model-based clustering method for R (MCLUST version 3) [59] to determine the number of morphological groups represented in the two sampling years. The data were fitted to models with one or a mixture of up to three Gaussian distributions. The MCLUST functions for univariate data consist of only two models with equal (denoted E) or varying variance (denoted V), respectively. The best model was selected based on the Bayesian Information Criterion (BIC). We calculated ΔBIC for each population by estimating the difference between BIC for the best and the next best model. Following Kass and Raftery [59], ΔBIC >10 was adopted as very strong support for the best model, 6< ΔBIC <10 as strong support, 2< ΔBIC <6 as moderate support, and ΔBIC <2 as equivalent support for the best and the next best model.

Microsatellite Analysis

The total DNA was extracted from gill lamellae using E.Z.N.A E-Z 96 Tissue DNA Kit (Omega-Biotek) following the manufacturer’s protocol. Individuals were genotyped at 20 polymorphic microsatellite loci: BFRO-018 [60], BWF1, BWF2 [61], Cla-Tet01, Cla-Tet03, Cla-Tet06, Cla-Tet09, Cla-Tet10, Cla-Tet13, Cla-Tet15, Cla-Tet17, Cla-Tet18 [62], Cocl-Lav04, Cocl-Lav06, Cocl-Lav10, Cocl-Lav18, Cocl-Lav27, Cocl-Lav49, Cocl-Lav52 [63], C2-157 [64]. The loci were co-amplified in four 2.5 µl PCR multiplex reactions as described earlier [65] using the QIAGEN Multiplex PCR kit following the manufacturer’s protocol. The PCR products were denatured in Hi-Di Formamide, containing LIZ-500 internal size standard (Applied Biosystems) and separated using an ABI-3130×l Genetic Analyzer (Applied Biosystems). The alleles were scored using the automatic binning function, with predefined bins, as implemented in Gene-Mapper 3.7 software (Applied Biosystems). All alleles were subsequently verified visually by two independent persons. In addition, the included replicate and blank samples were manually verified to ensure the validity of the data. Finally, we manually verified all identified private alleles to ensure that they are not an artefact from inconsistent scoring.

Statistical Analyses

Micro-Checker 2.2.3 [66], employing 1,000 bootstraps, was used to check for null alleles and genotyping errors. Loci showing null alleles were removed from the dataset before subsequent analysis.

Number of individual clusters (populations) at each of the two sampling years were inferred by STRUCTURE 2.2.3 [67]. We calculated the number of populations represented (K) and the individual admixture coefficients (q) of individuals belonging to the pure parental populations of respectively LSR, DR, or genetic hybrids of the two whitefish morphs. We used an admixture model with correlated allele frequency with burn-in period of 10,000 and number of iterations of 50,000 in five replicates. The statistical significance for the likely number of clusters in each of the two sampling years was tested by Kruskal-Wallis non-parametric test employed in SPSS [68] by comparing the ln Pr (X|Κ) values for each of 5 runs for variable K. We applied threshold values for q at 0.2 and 0.8 to categorise the individuals as hybrids or pure, with individuals having q-values between 0.2 and 0.8 being considered as hybrids and the others as pure individuals (q ≤0.2 and q ≥0.8). These threshold q-values have been shown by simulation studies to provide the most efficient detection of hybrids [69], as increasing the threshold q-value will misclassify the hybrids as pure individuals.

The populations defined from the q-values were used for further statistical analyses (see Table S1 and Table S2). Expected heterozygosity (He) was estimated using GENEPOP 4.1 [70], [71]. The exact test [72] implemented in GENEPOP 4.1 was used to test for deviations from Hardy-Weinberg and linkage equilibrium within each population, and the corresponding P (probability value) was obtained by following Markov chain conditions as 10,000 dememorisation, 200 batches and 10,000 iterations. Estimates were corrected for multiple tests using sequential Bonferroni corrections as described by [73]. Allelic richness (AR) was estimated accounting for differences in sample sizes using the rarefaction procedure [74] as implemented in HP-RARE [75]. The estimates calculated here were based on the smallest number of samples (N = 27 in Sb93_Hybrids). The R function cor.test [58] was used to test the association between the gill raker numbers and admixture values (q-values).

Results

Analysis of Gill Raker Numbers

The kernel density plots of gill raker numbers over the 15 years revealed a bimodal distribution with modes growing closer with time since invasion (Figure 2A and 2B). The MCLUST analysis of frequencies of gill raker numbers between the LSR and DR whitefish followed a bimodal distribution before and after the vendace invasion (2008: ΔBIC = 9; 1993: ΔBIC = 8.4).

Figure 2

Frequency distribution of gill raker numbers in DR-LSR whitefish morph-pair samples collected in 1993 (A) and 2008 (B) from Lake Skrukkebukta.

The solid line on the histogram indicates the kernel density function.

Genetic Analysis

Four loci showed departures from Hardy-Weinberg equilibrium after Bonferroni corrections (Cla-Tet10, Cla-Tet17, Cla-Tet15 and Cocl-lav52), most likely due to null alleles or stuttering as indicated in the Micro-Checker analysis. These loci were, therefore, discarded from subsequent analyses. Summary statistics for the microsatellite loci and populations are listed in Table S1 and Table S2.

The STRUCTURE analysis of temporal data confirmed the presence of two populations in the 1993 samples. Mean ln Pr (X|Κ) was significantly smaller for K = 2 (mean ln Pr (X|Κ) ± S.D. = −4467.58±9.91, Kruskal-Wallis test, P<0.05) as compared to either K = 1 (mean ln Pr (X|Κ) ± S.D. = −4505.46±0.76) or K = 3 (mean ln Pr (X|Κ) ± S.D. = −4567.96±57.94). However, in 2008, the most likely number of population was K = 1 (mean ln Pr (X|Κ) ± S.D. = −4706.78±0.55, Kruskal -Wallis test, P<0.05) rather than K = 2 (mean ln Pr (X|Κ) ± S.D. = −4824.6±19.28) or K = 3 (mean ln Pr (X|Κ) ± S.D. = −4812.42±17.18). The proportion of individuals assigned to one of the parental categories based on the threshold q-values (0.2/0.8) was 66% in 1993, mostly with non-overlapping 90% credible intervals (Figure 3A). In contrast, in 2008 all the individuals revealed q-values within the range from 0.2 to 0.8 with overlapping 90% credible intervals (Figure 3B).

Figure 3

Frequency distribution of individual admixture proportions of DR-LSR whitefish morph-pair upon arrival (1993, A) and after (2008, B) a biological invasion.

The q-values indicate admixture proportions with the associated 90% credible intervals as estimated by STRUCTURE 2.2.3 [67]. The solid horizontal lines indicate the threshold q-values for identifying hybrids and pure whitefish morphs.

The scatter plot of individual gill raker number and individual admixture coefficients (STRUCTURE q-values) showed two genetically and morphologically well differentiated groups of individuals in 1993 (Figure 4A). The individuals sampled in 2008, however, displayed no clear pattern of correlation between the individual admixture coefficients and gill raker number (Figure 4B). Support for this scenario was given by the correlation test for which a significant correlation between gill raker number and admixture values was observed in 1993 (Pearson correlation coefficient, r = 0.73, P<0.05), but not in 2008 (Pearson correlation coefficient, r = −0.177, P>0.05).

Figure 4

Scatter plots comparing genetic and morphological data of DR-LSR whitefish morph-pair at arrival (1993, A) and after (2008, B) a biological invasion.

The individual admixture values (q-values) estimated by Structure 2.2.3 [67] are plotted against the gill raker numbers. The q-values indicate admixture proportions of individuals. The solid horizontal lines indicate the threshold q-values to identify hybrids and pure whitefish morphs.

Discussion

The results of this study suggest a breakdown of reproductive isolation between DR and LSR whitefish within a 15-year period following the invasion of a superior trophic competitor, vendace. Most noteworthy, a Bayesian admixture analysis of the genetic data failed to separate the DR and LSR whitefish in the 2008 sample, whereas a clear separation of the two populations was evident in 1993. Furthermore, population-pure individuals, i.e. genetically assigned as either LSR or DR whitefish, were nearly absent in 2008 compared to a frequency of approximately 66% pure individuals revealed in the 1993 sample. Also, the association between gill raker numbers and genotypic data was strong and significant in 1993, whereas a much weaker and non-significant correlation was observed in 2008. This reflects a breakdown of the association between the genotype and the phenotypic trait. The high frequency of hybrids between the LSR and DR whitefish in 2008 suggests that hybridization may have contributed to a decrease in the abundance of the DR whitefish observed after the vendace invasion [35].

Increased hybridization should presumably promote a collapse in the divergent pattern of gill raker numbers in the two morphs since hybrids would be expected to consist of intermediate phenotypes in respect to this heritable trait. In lakes harbouring reproductively isolated LSR and DR whitefish, the gill raker numbers generally follow a distinct bimodal distribution pattern [30], [32]. In the present study, the LSR and DR whitefish showed a bimodal distribution pattern of the gill rakers in both sampling years, even though all individuals in 2008 were assigned as hybrids according to the genetic analyses. A likely explanation for this apparent contradiction may be the few generations passed since the first vendace observation (i.e. approximately three generations over the 15 year time), not allowing enough time for intermediate phenotypes to accumulate (but see the discussion part on association between gill raker numbers and genotypes). Notwithstanding, in the 2008 sample the two modes in the gill-raker distribution representing DR and LSR whitefish, are less distinct and located closer to each other compared to the 1993 sample. This may suggest an on-going process towards a collapse in the bimodal gill-raker number distribution. An alternate, albeit non-exclusive, explanation could be selection against intermediate phenotypes as the gill raker number is a trophic trait influenced by natural disruptive selection ( [23], but see [76]). The gill raker number has a polygenic basis [77] and is expected to follow a complex pattern of inheritance and expression relative to other more discrete traits. The breakdown of the bimodal distribution pattern of gill raker numbers in hybridizing whitefish may, therefore, take longer than the rapid collapse in the divergent distribution of neutral microsatellites, as observed in the present study. This is because the homogenizing effect of hybridization may be delayed by the counteracting effect of disruptive natural selection. Such a mechanism was also suggested from a study involving two ecotypes of rainbow smelt (Osmerus mordax), in which natural selection maintained the adaptive phenotypic differences (including gill raker counts) in spite of the presence of high gene flow [78].

A genetic analysis (five microsatellite loci) combined with analysis of a phenotypic trait (body shape) spanning over a 25 year period, revealed a collapse of a limnetic-benthic species-pair of three-spined sticklebacks into a hybrid swarm [17]. Using STRUCTURE analysis as a proxy for departures from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD), the study also showed a decline in linkage disequilibrium across years following the collapse. This decline was not evident in the pairwise loci comparisons nor within years. Likewise, the authors did not identify any decrease in observations of HWE after the collapse. In the present study involving 16 microsatellite loci, no deviations from HWE and presence of LD among loci were identified. However, the STRUCTURE output and scatterplot of q-values corroborate assumptions of a weakening of genetic differentiation between the whitefish morphs after the vendace invasion. Thus, this observation supports our expectation of loss of genetic differentiation between the LSR and DR whitefish.

The whitefish morphs are differentiated based on the number and morphology of the gill rakers as well as their morphological appearance [30], [37], [52]. Microsatellite studies have repeatedly reported genetic differentiation between DR and LSR whitefish in various lakes in northern Norway [21], [34]. Hence, there is reason to assume that hybridization and gene flow between the sympatric whitefish morphs should promote a collapse of the association between genetic and phenotypic traits (in our case gill raker numbers) due to increased admixture. In the analysis of the three-spined stickleback species-pairs from Enos Lake, the hybrid swarms showed no or little association between body morphology and genetic scores [17]. Herein, the lack of associations between gill raker numbers and q-values in 2008 agrees with the above findings and supports our expectation of a breakdown of the association between morphology and genotype frequencies in the whitefish.

Nonetheless, information about spawning behaviour and the nature of reproductive isolation mechanisms between the DR and LSR whitefish in northern Fennoscandian lakes remain speculative at this time. However, Table 1 summarizes putative reproductive isolation mechanisms between DR and LSR whitefish morphs based on available information from whitefishes and other postglacial fish species. Time and/or place of spawning are considered as an important prezygotic isolation mechanisms between whitefish morphs [7], [46], similarly as for other postglacial fishes like Arctic charr [22]. However, these mechanisms do not seem to be similar across different study lakes [46], [79]. Local lake conditions such as fluctuations in water temperature, bottom substrate and depth may be important factors in determining spawning time and place [46]. Svärdson [37] suggested that the differences in spawning time or place are important reproductive isolation mechansims in whitefish in Swedish lakes. A similar mechanism has been reported in Alpine whitefish species, which spawn at different depths ( [7] and Table 1). In Lake Stuorajavri (northern Norway), the LSR whitefish comprised 98% of the spawning whitefish caught at two different spawning grounds (R. Knudsen, unpublished), suggesting that the two morphs at least have separate spawning locations.

The fact that sympatric adults of DR and LSR whitefish differ in size suggests that size-assortative mating may be an important mechanism maintaining reproductive isolation in the face of gene flow (Table 1). Size assortative mating has been described as a pre-zygotic isolation mechanism between kokanee and sockeye salmon (Oncorhynchus nerka), as well as in other whitefish sympatric forms (Coregonus spp) [37], [48]. Microhabitat segregation of spawning site may be another possible mechansims of premating reproductive isolation as suggested by various studies [17], [45], [80].

Information regarding extrinsic and intrinsic post-zygotic reproductive isolation in coregonids exists from studies of North American lake whitefish [50], [51] and whitefish (Coregonus spp.) from the Swiss Alps [20]. Furthermore, severe genomic incompatibilities and related differential embryonic mortality have been reported in the hybrids of parents belonging to two glacial races of lake whitefish ( [50], [51], Table 1). However, no such intrinsic hybrid inviability was reported from a study of Alpine whitefish although the hybrids showed discrepancies in hatching time [20]. Hence, the authors suggested that hatching time could be a powerful selective force against hybrids in Alpine whitefish [20]. Another extrinsic post-zygotic reproductive isolation mechanism relates to ecological selection against hybrids or ecological inviability. Ecological inviability is imputable to reduced foraging efficiency due to the intermediate gill rakers or absence of appropriate ecological niches [19]. Considering the young evolutionary age (11700-5800 years before present) of the northern Fennoscandian whitefish monophyletic lineage [21], [33] and previous results in Alpine whitefish [20], it is plausible that the speciation process has not proceeded far enough to develop a strong or complete reproductive isolation between the sympatric whitefish morphs. The on-going or incomplete speciation process is a key signature of young northern post-glacial systems as compared to systems of an older evolutionary age. We propose, therefore, that our results suggest a breakdown of reproductive isolation between sympatric whitefish morph-pairs, further leading to a hybrid swarm triggered by the competitive interaction between native whitefish and non-native vendace.

The genetic architecture of the reproductive isolation mechanisms is largely unknown in European whitefish. In ecological speciation, the traits offering reproductive isolation must either be genetically correlated with traits under disruptive selection, or the traits offering the adaptive advantage should drive the reproductive isolation pleotropically. In such a situation, reproductive isolation may be achieved as a result of alleles that allow adaptation to a particular niche also will govern the mate choice [81]. Otherwise, the progress towards the evolution of reproductive isolation between resource specialists would require a build-up of non-random association between the set of genes optimal for each resource and the genes for assortative mating. In other words, reproductive isolation requires linkage disequilibrium between the genes affecting viability and those affecting assortative mating [81]. This suggests that the genes responsible for gill raker numbers in European whitefish must be in association with the genes responsible for intrinsic reproductive isolation. While this has not been investigated in European whitefish, a recent study on North American whitefish showed that covariation in expression of genes belonging to the same gene network was associated with phenotypic variation both in swimming behaviour and water depth preference and in number of gill rakers [82]. Since differential habitat preference (e.g. water depth) could promote reproductive isolation, and since both depth preference and number of gill rakers co-vary in association within the same gene network, these results on lake whitefish provide the first evidence that gill raker numbers could be involved in reproductive isolation between sympatric whitefish forms.

The observed hybridization and apparent re-admixture of DR and LSR whitefish following the vendace invasion in the Pasvik watercourse may represent an example of “speciation in reverse” [1], [4]. A few other examples of “speciation in reverse” in fishes have been described from natural systems inhabited by Alpine whitefish, three-spined sticklebacks, ciscoes and cichlids [4], [15], [17], [25,]. However, the reason for the collapse is not known in most of the reported cases and is only suggested in the collapse of morph-pairs in three-spined sticklebacks from Enos Lake [16], [17] and from species reversal in Alpine whitefish [15]. In the former case, the collapse of the species pair has been attributed to a coinciding establishment of American signal crayfish (Pascifasticus lenisculus) that may have altered the habitat use, mating behaviour and/or trophic distribution of sticklebacks in the lake. In the latter case, Vonlanthen et al. [15] suggested that anthropogenic induced eutrophication has partly destroyed the spawning grounds and caused relaxation of divergent selection thereby facilitated the collapse of native whitefish species-pairs in several Alpine lakes. Thus, there is reason to suggest that the relegation of DR whitefish into the habitat of LSR whitefish with a subsequent genetic homogenization may represents speciation reversal.

At the genomic level, the observed pattern of putative speciation reversal may be attributed to differential divergence at the local or global genomic scale, referred to as growth of genomic islands [83]. Genomic islands are genomic regions exhibiting greater differentiation than expected under neutrality. An ecological speciation model (with gene flow) of “genomic islands of divergence” predicts that the genes affecting the local adaptation and reproductive isolation may reside within the genomic clusters of islands [84]. Correlated response to selection particularly by virtue of linkage disequilibrium is most likely to reverse after the event of secondary contact. Recombination will break the larger genomic islands into smaller parts, but the genes under selection remains isolated in smaller linkage regions. In this scenario, the process may not be termed “speciation reversal” but would rather reflect the decay of genomic islands into smaller numbers and sizes. However, both scenarios are not necessarily exclusive as decay of genomic islands of divergence could actually be the key molecular mechanisms causing speciation reversal. This, however, remains to be investigated.

Conclusions

Our study suggests an apparent breakdown of reproductive isolation between two sympatric European whitefish morphs over a relatively short time span of 15 years (i.e. within three generations) following the invasion of a superior resource competitor (vendace) in a subarctic lake. Further studies are necessary to understand the mechanisms and dynamics that are weakening the reproductive isolation and the potential evolutionary consequences for the whitefish populations in the lakes of the Pasvik watercourse. In addition, a genome wide study would be necessary to shed more light on the effect and dynamics of admixture at the genomic level.

Summary statistics for each locus and four populations from 1993 and 2008.

(DOCX)

Click here for additional data file.

Details of population and basic genetic diversity measures in four populations from 1993 and 2008.

(DOCX)

Click here for additional data file.