The Effect of Recurrent Floods on Genetic Composition of Marble Trout Populations

A changing global climate can threaten the diversity of species and ecosystems. We explore the consequences of catastrophic disturbances in determining the evolutionary and demographic histories of secluded marble trout populations in Slovenian streams subjected to weather extremes, in particular recurrent flash floods and debris flows causing massive mortalities. Using microsatellite data, a pattern of extreme genetic differentiation was found among populations (global F ST of 0.716), which exceeds the highest values reported in freshwater fish. All locations showed low levels of genetic diversity as evidenced by low heterozygosities and a mean of only 2 alleles per locus, with few or no rare alleles. Many loci showed a discontinuous allele distribution, with missing alleles across the allele size range, suggestive of a population contraction. Accordingly, bottleneck episodes were inferred for all samples with a reduction in population size of 3–4 orders of magnitude. The reduced level of genetic diversity observed in all populations implies a strong impact of genetic drift, and suggests that along with limited gene flow, genetic differentiation might have been exacerbated by recurrent mortalities likely caused by flash flood and debris flows. Due to its low evolutionary potential the species might fail to cope with an intensification and altered frequency of flash flood events predicted to occur with climate change.


Introduction
Climate change poses a serious threat to species persistence. Many species will experience selection in new directions and at new intensities, and the degree to which species respond adaptively to dynamic and uncertain futures will determine their survival over the coming decades and millenia [1]. Intensification of weather extremes and associated catastrophic disturbances is emerging as one of the most important aspects of climate change and research is advancing from studying the impacts of changes in mean climate values (trend effects) to those produced by changes in the magnitude or frequency of extreme events (event effects) [2]. Catastrophes are characterized by their statistical extremeness combined with a short duration relative to the life cycle of the organisms affected; they can disrupt ecosystem, community or population structure and change resources, substrate availability, or the physical environment [2][3][4]. Evidence suggests that the frequency and intensity of extreme weather events (i.e. floods, droughts) is increasing in many regions in response to global climate change [5][6][7]. Despite the urgent need to advance research on extreme events and catastrophic disturbances, their evolutionary consequences have largely been unexplored [8].
Freshwater salmonids are commonly subject to substantial environmental variability in the form of changes in stream flow at different time scales [9][10][11]. Extreme events such as floods, droughts or landslides play an important role in the regulation of population dynamics in salmonids, to the extent that in highdensity streams, fish demography and persistence might be largely shaped by extreme flow events [12][13][14][15]. The direct and short-term effects of floods are largely a result of high-water velocities and sediment movement that cause the displacement and death of fish. The impact of floods is expected to be higher in the coming years. According to IPCC predictions, an increase of rainfall is expected in the next 50 years along with an intensification and altered frequency of catastrophic disturbances [5]. Recent advances in the statistical theory of extreme events suggest that large flood events will be more frequent, with return times markedly shorter than expected even in the absence of climate change [16].
Our model system is the stream-dwelling salmonid marble trout Salmo marmoratus, an endemism of the Adriatic Sea and its tributaries, currently restricted to Northern Italy, former Yugoslavia (Slovenia, Croatia, Bosnia-Herzegovina) and Albania [17]. Marble trout is considered to be one of the most endangered species in the Adriatic basin [18]. For decades, massive restocking practices have been conducted throughout its natural range by means of introduction of exotic brown trout. Hybridization between marble and brown trout has been so extensive that hybrids now dominate most rivers [19]. Molecular data confirm a high level of introgression in the Po river in Northern Italy [20,21], the Soca river system in the Italian/Slovenian border [18,22,23], and recently, the Adige river system in South Tyrol [24,25]. However, all studies also reported pure populations of marble trout in headwaters of all river systems, persisting above natural barriers (i.e. waterfalls). While those barriers prevent the upstream migration of conspecifics and consequent hybridization, the secluded nature and the small size of the remnant marble trout populations makes them extremely vulnerable to stochastic factors, including environmental events such as floods, droughts or landslides.
A project for the conservation and rehabilitation of marble trout in Slovenia started in 1993 [17]. Only seven remnant pure (nonintrogressed with brown trout; [23]) populations of marble trout remain in the Adriatic basin of Slovenia, persisting in totally isolated headwaters without predators or fishing [17]. Recurrent major floods and debris flows are the most important risk factor for the viability of marble trout populations in Slovenian streams [15,26]. The trigger of debris flows in the area is extreme precipitation, with records of intense flows going back to the 18 th century, and with a presumable occurrence interval of approximately 50 years (i.e. in the 20 th century major floods were recorded in 1926, 1954 and 1990; [27]). However, an intensification of the frequency of flash floods has been observed in the last decade, with four important flood events recorded in the 1999-2010 period ( Table 1).
The aim of this study is to explore the impact of catastrophic weather events on the genetic composition of isolated marble trout populations from the Adriatic basin of Slovenia ( Figure 1). As a result of intense precipitation and associated flash flood and debris flows, which causes mass mortalities but does not alter connectivity among isolated streams and basins, low levels of genetic diversity and high genetic substructuring at local geographic scale are expected. In this study, we were particularly interested in the comparison of pre-and post-flood samples, testing shifts in genetic composition and genetic diversity and estimating the possible existence of bottleneck/population declines.

Results
A total of 18 out of 24 loci were polymorphic, showing from 2 to 8 alleles per locus ( Table 2). All locations showed low levels of genetic diversity as evidenced by heterozygosities ,0.25 and allelic richness ,2 (Table 3). Many loci showed a discontinuous allele distribution, with missing alleles across the allele size range, and few or no rare alleles (Table S1). A total of 31 location-specific private alleles were found, representing 49% of the total number of alleles sampled. Across locations, Studenc showed the highest genetic diversity, with intermediate levels found in Zadlascica and the lowest genetic variability found in Lipovscek and Sevnica (Table 3). When comparing pre-and post-flood samples, a moderate drop in genetic diversity was observed in all four locations in terms of H o and H e and allele richness (AR), mostly due to the loss of rare alleles (Table 3; Table S1). However, statistic comparisons between pre-and post-flood values of genetic A highly significant extreme overall genetic differentiation was found among samples (F ST = 0.716; p,0.001). All pairwise F ST and genetic distances between samples from different locations were high (F ST = 0.649-0.779; D CE = 0.592-0.839) and statistically significant (p,0.001; Table 4). By contrast, comparison of preand post-flood samples from the same location showed low nonsignificant F ST and genetic distances with the exception of Studenc (F ST = 0.041, p = 0.020; D CE = 0.021, p = 0.035). Accordingly, comparison of allele frequencies between pre-and post-flood samples using an exact test showed no temporal differences at Lipovscek (p = 0.993), Sevnica (p = 0.450) and Zadlascica (p = 0.985), while significant differences were found at Studenc (p,0.001) caused by loci CA059136 and Str151.
Accordingly, an AMOVA analysis partitioned genetic differentiation significantly among locations (F CT = 0.713; p,0.001) but not among (pre-and post-flood) samples within locations (F SC = 0.011; p = 0.067) ( Table 5). We conducted a Multidimensional Scaling analysis using D CE at all loci ( Figure 2). When plotting the values of the first and second principal components, samples clustered according to location, the Zadlascica location (from a tributary of the Soca river) appeared very different from the rest of locations (from tributaries of the Idrijca river), but so did Studenc despite being 20 km apart from Sevnica. In this sense, when testing Isolation-by-Distance, no correlation was found between genetic and geographic (shortest waterway) distances (r = 0.715; F = 4.19; p = 0.110). Assignment analysis using STRUCTURE confirmed the extreme genetic differentiation among locations and suggested a scenario with K = 4 groups as the most likely (p,0.001), corresponding to the four sampled locations, Zadlascica, Lipovscek, Sevnica and Studenc.
The MSVAR procedure for assessing past demographic history strongly supported a decline in all marble trout populations. All sampled points of log 10 (r) were substantially below zero in all eight samples, suggesting that the past population size was larger than the current population size ( Figure 3). All sampled points of log 10 (h) were negative, pointing to small current effective population sizes. Using a conservative mutation rate of m = 10 24 ,  (Table 7) using the sofware BOT-TLENECK and M_P_VAL showed some indications of population contraction at all locations. Garza and Williamson [28] suggested that values of M lower than 0.7 would indicate evidence of a bottleneck, while values above 0.8 would denote no bottleneck history. In our data set, all Zadlascica and Studenc samples plus the pre-flood Lipovscek sample showed M values between 0.518-0.661. The two Sevnica samples plus the post-flood Lipovscek sample showed M values between the 0.7-0.8 limits.

Extreme differentiation among marble trout populations
Freshwater fish species show a greater average degree of genetic differentiation among locations than marine species, resulting from the isolation of fish populations among drainages [29]. The marked zoogeography produced by the historical patterns of isolation among drainages is subjected to the continuous remodeling of the river drainage and to climatic fluctuations of which the glacial/interglacial periods of the Pleistocene played a crucial role [30]. Species respond actively to fluctuations in their natural range, with hydrographic networks being used for colonization and for retreat throughout geological times. Higher levels of population subdivision have been found in many salmonid studies reflecting complex genetic structures often resulting from isolation among drainages. In brook charr Salvelinus fontinalis, a global F ST of 0.37 was found among 26 populations from La Maurice National Park in Canada located 3-42 km apart [31]. In cutthroat trout Oncorhynchus clarkii, Taylor et al. [32] found an overall F ST of 0.32 and pairwise F ST values up to 0.45 in populations from British Columbia, while the recent study of Pritchard et al. [33] reported a global F ST of 0.41 for Rio Grande populations. In bull trout Salvelinus confluentus, Taylor and Costello [34] found an F ST = 0.33 among 20 North-West America coastal locations, while recent papers reported an overall F ST of 0.15 and pairwise F ST values up to 0.31 in Alberta [35] and pairwise F ST values up to 0.66 between lake samples in a 50 km area off Montana [36].
Collectively, the extreme level of genetic differentiation found among Slovenian marble trout populations (overall F ST = 0.716; pairwise F ST between drainages = 0.649-0.779; 49% private alleles sampled) exceeds the highest values reported in salmonid populations found above waterfalls and/or in small streams like the ones in this study. Such extremely high F ST is indicative of complete genetic isolation that has persisted over time. This is concordant with the previous study of Fumagalli et al. [23], which reported a global F ST of 0.66 in the same area using a smaller number of microsatellites. At present, the four locations in our study are completely isolated from each other and from the main river drainage by means of impassable natural barriers (waterfalls) that restrict dispersal abilities. Baloux and Lugon-Moulin [37] argued that even in the total absence of gene flow, extreme genetic differentiation of the magnitude found in our study is not expected due to the high mutation rate of microsatellites [38], with appearance of new variants counteracting within-population fixation of alleles. However, the low level of genetic diversity observed in all locations (low heterozygosities and low number of alleles) implies a strong impact of local genetic drift, and suggests  that, together with absence of gene flow, genetic differentiation among locations may have been exacerbated by recurrent mortalities likely caused by flash flood and debris flows. The combined effects of genetic drift and inbreeding following demographic bottlenecks may have thus contributed to alter distinctly the genetic composition of each population, so that over time different populations ended up having different alleles at each locus as well as different allele combinations from multiple loci, resulting in the observed pattern of extreme genetic differentiation. The MSVAR procedure clearly suggested a population contraction with a reduction in population size of 3-4 orders of magnitude to current effective densities of around 3-41 (CI: 1-396) individuals per location. Bottleneck signatures were observed in all samples. First, most polymorphic loci presented a biallelic pattern with few (or none) rare alleles. Second, a discontinuous allele range was found at many loci (Table S1). For instance, only alleles 104, 110, 118 and 124 were sampled across locations at locus Strutta58, which suggests that the non-observed alleles have been lost over time. By contrast, samples collected in tributaries of the river Po in Northern Italy showed a continuous allele distribution and a higher number of alleles at five microsatellite loci [39], suggestive of demographically stable populations with no signatures of bottlenecks. Loci Ssa85 and Str15, which were monomorphic in Slovenian locations (our study), showed 6-9 alleles in Italian locations. Loci Str85, Str543 and Str591, which presented three alleles each in our study, showed higher allelic richness in the Italian samples (7-17 alleles). Pujolar et al. In contrast with the extreme differentiation found at microsatellite loci, one single mitochondrial DNA control region haplotype (MA1) has been reported for marble trout in Slovenia [22,39]. This suggests that despite the current complete geographic isolation of marble trout populations in Slovenian streams, those populations were interconnected in a relatively recent past prior to the fragmentation of habitats. On the basis of microsatellite data, Fumagalli et al. [23] suggested that the pure populations within the Idrijca drainage belonged to an independent river system. Nevertheless, our analysis on a larger number of microsatellite loci does not support the division by river basins, as the two locations from the Idrijca basin in our sampling (Sevnica and Studenc) did not cluster together in the Multidimensional Scaling analysis, and appeared clearly differentiated despite being merely 20 km apart within the same river basin.

Genetic consequences of serial bottlenecks
The analysis of pre-and post-flood samples in our study allowed us to explore the genetic consequences of population bottlenecks. The effects of bottlenecks are directly related to the increase of Table 5. Results from AMOVA analysis partitioning genetic differentiation among locations (4 locations: Zadlascica, Lipovscek, Sevnica and Studenc) and among samples within locations (pre-and post-flood samples).  Table 4. Pairwise F ST estimates (above diagonal) and genetic distances (D CE , below diagonal) between samples, labeled as in Table 3. stochastic events associated with small population sizes, leading in most cases to a loss in genetic variability [40]. Recently, Bouzat [41] suggested that the potential genetic outcomes of demographic bottlenecks can only be assessed when considering replicated bottlenecked populations. In our study, we explored four separate locations that experienced episodes of massive mortalities with consequent reductions in population size ranging from 56 to 77%. Collectively, and taking into the consideration the magnitude of the demographic decline, only a limited genetic effect was observed. A moderate drop in genetic diversity was found in all locations, both in terms of heterozygosities and allelic richness, but all comparisons were statistically not significant. However, a total of 9 alleles were lost locally, in all the cases representing rare alleles with frequencies ,0.05 prior to the flood that were no longer sampled after the flood events. Moreover, bottlenecks resulted in alteration of allele frequencies, particularly secondary alleles, which in some cases experienced a drop of 20-30% in frequency. This was notably observed at Studenc, where genetic composition was significantly different in pre-and post-flood samples, but not at the other populations. The moderate reduction of allelic diversity and heterozygosity might be attributable to the particular demographic history of those populations. In his recent review, Bouzat [41] emphasized the potential role of population history in determining the outcome of demographic bottlenecks. Specifically, one can expect that populations experiencing recurrent bottlenecks might have had their genetic pool already eroded over time, which would decrease the effectiveness of both purifying selection and random allele elimination. This holds particularly true for the Slovenian marble trout populations in our study, which have been repeatedly impacted by severe flood events, with a presumable occurrence interval of 50-100 years [27]. Bayesian demographic analysis using the MSVAR procedure is concordant with a scenario of serial bottleneck episodes that might have been occurring for 150-1340 years.
On the question of how do the Slovenian population still maintain some polymorphism despite experiencing flood events for hundreds of years, one possible explanation is that our sampling is not representative of the entire population, and other compartments of the population might keep some genetic variation. One hypothesis could be that young fish (age 0+ and 1+) might be less sensible to flows because usually they are not found in the main stream but in more protected small tributaries, repopulating the stream even if all adults have been removed. However, field observations suggest that eggs and young fish are more vulnerable to floods than adults and that there is no recruitment at all when floods occur after reproduction (Crivelli, unpublished information). Alternatively, an upstream part of the population could be preserved. This hypothesis is supported by tagging data (Crivelli, unpublished information) that show nontagged individuals in post-flood samples, whereas all pre-flood samples had been tagged, which suggests that the newly-arrived non-tagged individuals came from a compartment of the population located upstream, flushed by the flood.

Coping with environmental change
Our findings strongly suggest that the evolutionary and demographic history of marble trout populations living in secluded Slovenian streams has been shaped by massive mortalities caused by disturbance factors such as flash flood and debris flows. Overall, genetic analysis revealed an extremely high genetic differentiation among remnant pure populations, together with much low levels of genetic diversity in all populations. Effective population sizes estimated using the MSVAR procedure ranged from 3 to 41 (CI:1-396) individuals per population. The low genetic variability found does not seem to affect the viability of the populations so far, as they have persisted up to the present despite recurrent and unpredictable disturbances. To mitigate their ecological impact, fish populations might exhibit several adaptations that result from trade-off among growth, reproduction and survival. According to life-history theory, recurrent floods can have important evolutionary consequences by selecting for life-histories that are synchronized to either avoid or exploit the direct and indirect effects of extreme flows [42]. However, due to the predicted intensification of extreme weather leading to increased floods, heatwaves, droughts and rainfall in the next 50 years [5,43], marble trout may not have the sufficient genetic potential for the evolution of life-histories able to mitigate their impact. Recent studies by the European Commission's Joint Research Centre (JCR) have confirmed rising temperatures (around 1.5uC in the last 35 years) Figure 3. Plots of the simulated points from the marginal posterior distributions of log 10 (r) (x-axis) and log 10 (t f ) (y-axis). Samples labeled as in Table 3. Collectively, while the Slovenian marble trout populations have coped successfully with recurrent yet unpredictable disturbance events over time, the low genetic variation found in the species makes it difficult to assess its evolutionary potential since a genetically depauperate population might fail to adapt to future environmental change. Moreover, all populations are closed, and there is no potential for spontaneous colonization of new habitats or re-colonization after local extinction. It has been proposed that an effective population size of 500 individuals is large enough to maintain genetic diversity for key history traits [45]. The estimation in our study of effective population sizes of 3-41 individuals is one-two orders of magnitude lower than the threshold of 500 individuals, plus the upper-bound CI of population size does not overlap with the threshold value. This suggests that the Slovenian marble trout populations are beyond the critical population size for genetically secure populations. The future development of an integrated genetic/ecologic model that explores different demographic scenarios with varying degree of intensity and frequency of flood events might help disentangling the role of catastrophic disturbances in determining the genetic structure and genetic diversity of marble trout.  Table 6. Summary statistics of MSVAR across all locations, including the genetic parameter h = 2N 0 m, which was calculated using a mutation rate of m = 10 24 , r = N 0 /N 1 (N 0 = current effective population size; N 1 = effective population size at the time of population expansion or decline), and t f = t a /N 0 (t a = number of generations since the beginning of the expansion/decline).

Microsatellite amplification
Minute sections of tissue from ethanol-preserved finclips were digested in a lysis buffer containing 100 ml TE Buffer, 7 ml 1 M DTT (dithiothreitol) solution pH 5.2 (diluted in 0.08 M NaAC) and 2 ml proteinase K solution (20 mg/ml) for at least 8 hours at 56uC. After incubation at 96uC for 10 min, samples were centrifuged at 13,000 rpm for 11 min, and the supernatant was stored at 220uC.
PCR products were obtained in a GeneAmp PCR System 2700 Thermocycler (Applied Biosystems) using the QIAGEN Multiplex PCR Kit. PCR reactions consisted of 2 ml template DNA, 5 ml QIAGEN Multiplex PCR Master Mix, 0.2 ml 10 mM forward and reverse primers, and water up to 10 ml. PCR conditions were as follows: 3 min at 95uC, 35 cycles of 30 sec at 94uC, 90 sec at 57uC and 1 min at 72uC, and final elongation for 5 min at 60uC. PCR products were visualized in 1.8% agarose gels and screened for microsatellite polymorphism using an ABI 3130 AVANT automatic capillary sequencer (Applied Biosystems). Alleles were sized according to a Liz500 (50-500 bp) marker.

Data analysis
Within-sample genetic diversity statistics were assessed by observed (H o ) and expected (H e ) heterozygosities per locus using GENETIX version 4.05 [54] and allelic richness (AR) using FSTAT [55]. Differences in genetic diversity among samples were tested by one-way ANOVA using STATISTICA version 6.0 (StatSoft Inc.). Deviations from Hardy-Weinberg Equilibrium (HWE), linkage disequilibrium and differences in allele and genotype frequencies among samples were tested using GENE-POP version 3.4 [56].
Neutrality of the markers was tested using the software LOSITAN [57], which implements a F ST outlier detection approach. This method evaluates the relationship between F ST and expected heterozygosity (H e ), describing the expected distribution of Wright's inbreeding coefficient F ST vs. H e under an island model of migration with neutral markers. This distribution is used to identify outlier loci that show excessive high or low F ST compared to neutral expectations. Such outlier loci are candidates for being subject to selection. We used the 0.99 criterion in order to minimize false positives as suggested by the authors.
Population structure was explored using non-hierarchical and hierarchical F-Statistics [58] calculated using ARLEQUIN [59]. First, overall and pairwise F ST values were calculated. Genetic differentiation was also partitioned among locations (4 locations: Zadlascica, Lipovscek, Sevnica and Studenc) and among samples within locations (pre-and post-flood samples). Significance tests were assessed with 10,000 permutation tests. In all cases, significance levels were corrected for multiple comparisons using Bonferroni [60]. Pairwise multilocus comparisons between samples were calculated by Cavalli-Sforza and Edwards [61] chord distance (D CE ) and graphically represented by Multidimensional Scaling (MDS) analysis using STATISTICA version 7.0 (StatSoft). Isolation-by-Distance (IBD) was tested using a Mantel test implemented in GENETIX, by correlating linearized genetic distance (F ST /(12F ST )) vs. geographic distance (shortest waterway distance measured along the streams between pairs of samples).
A model-based clustering algorithm as implemented in the software STRUCTURE [62] was used in order to infer the most likely number of populations in the data. The software organizes individuals into a predefined number of clusters (K) with a given likelihood, which may represent putative groups. The analysis was performed with 1,K,8 to account for population substructuring within species, using the admixture model and without a population prior. The most likely K was determined using the criterion of Evanno et al. [63] and then used to assign each individual. A burn-in length of 10 5 iterations followed by 10 6 additional Markov Chain Monte Carlo (MCMC) iterations were performed. A minimum of 5 runs were performed for each K to check repeatability of results.
Historical demographic changes were inferred using the Bayesian coalescent-based approach implemented in MSVAR [64]. Using a strict stepwise mutation model, the procedure provides distributions of the exponential growth rate r = N 0 /N 1 (where N 0 is the present effective population size, N 1 is the effective population size at the time of population expansion or decline), the time since the population started to expand or decline t f = t a /N 0 (where t a is the number of generations since the beginning of the expansion/decline) and the genetic parameter h = 2N 0 m. A mutation rate of m = 10 24 was used [65,66]. Rectangular priors were chosen for all parameters, with limits of (29, +5) for log 10 (r), log 10 (t f ) and log 10 (h) and starting values of 1 for h at each locus, r and t f . The analysis was performed for an exponential model of population change. We used 50,000 thinned updates and a thinning interval of 50,000 steps, with an initial 10% discarded as burn-in. Convergence was assessed using Tracer v1.4 [67].
Finally, we tested for a recent genetic bottleneck episode with two different approaches. First, we used the software BOTTLE-NECK [68], based on the principle that after a recent reduction of effective population size, number of alleles (k) decreases faster than heterozygosity (H e ) at polymorphic loci. Thus, in a recently bottlenecked population, the observed gene diversity is higher than the expected equilibrium gene diversity (H eq ) which is computed from the observed number of alleles (k), under the assumption of a constant-size (equilibrium) population [69]. We used the Multiple-step Stepwise (TPM) model [70], which consists of mostly one-step mutations but a small percentage of multi-step changes, and is the recommended model for microsatellite data sets rather than the Infinite Alleles (IAM) or Single-step Stepwise (SMM) models [71]. The proportion of singlestep mutation events was set to 90% (variance = 12%). Observed and expected heterozygosities were compared using a Wilcoxon sign-rank test as suggested by Piry et al. [68]. Second, we calculated M, the mean ratio between number of alleles (k) and range in allele size (r), assuming that during a bottleneck episode k decreases faster than r (M_P_Val; [28]). Hence, the value of M decreases when a population is reduced in size. Average M was calculated across loci and compared with the critical value M crit estimated after 10,000 simulations and assuming the population to be at equilibrium. In all simulations, three different values of h were used (5, 10 and 20). A range of mutation models were examined and conservative values were used for p s (frequency of one-step mutations) and D g (average size of non one-step mutations), p s = 0.90 and D g = 3.5, respectively.

Supporting Information
Table S1 Allele frequencies at all loci. Samples labelled as in Table 3. (DOC)