Genetic Impact of a Severe El Niño Event on Galápagos Marine Iguanas (Amblyrhynchus cristatus)

The El Niño-Southern Oscillation (ENSO) is a major source of climatic disturbance, impacting the dynamics of ecosystems worldwide. Recent models predict that human-generated rises in green-house gas levels will cause an increase in the strength and frequency of El Niño warming events in the next several decades, highlighting the need to understand the potential biological consequences of increased ENSO activity. Studies have focused on the ecological and demographic implications of El Niño in a range of organisms, but there have been few systematic attempts to measure the impact of these processes on genetic diversity in populations. Here, we evaluate whether the 1997–1998 El Niño altered the genetic composition of Galápagos marine iguana populations from eleven islands, some of which experienced mortality rates of up to 90% as a result of El Niño warming. Specifically, we measured the temporal variation in microsatellite allele frequencies and mitochondrial DNA diversity (mtDNA) in samples collected before (1991/1993) and after (2004) the El Niño event. Based on microsatellite data, only one island (Marchena) showed signatures of a genetic bottleneck, where the harmonic mean of the effective population size (Ne) was estimated to be less than 50 individuals during the period between samplings. Substantial decreases in mtDNA variation between time points were observed in populations from just two islands (Marchena and Genovesa). Our results suggests that, for the majority of islands, a single, intense El Niño event did not reduce marine iguana populations to the point where substantial neutral genetic diversity was lost. In the case of Marchena, simultaneous changes to both nuclear and mitochondrial DNA variation may also be the result of a volcanic eruption on the island in 1991. Therefore, studies that seek to evaluate the genetic impact of El Niño must also consider the confounding or potentially synergistic effect of other environmental and biological forces shaping populations.


INTRODUCTION
The El Niñ o-Southern Oscillation (ENSO) refers to a complex set of ocean-atmosphere interactions that take place throughout the Pacific basin [1]. El Niñ o events, which represent one phase of ENSO, are characterized by an accumulation of warm surface water in the central and eastern regions of the tropical Pacific. The resulting high sea-surface temperatures (SST) feed back on atmospheric circulation patterns worldwide, causing a range of environmental changes, from severe droughts in parts of Asia and the Western Pacific, to harsh winter conditions and flooding in North America [2][3].
During periods of acute El Niñ o warming, there is widespread mortality in aquatic organisms in the eastern Pacific [4]. The suppression of ocean upwelling and the rise in SST cause a marked decrease in primary productivity with ecological consequences for the entire marine food web [5]. Paleoclimate data depicts a 15,000 year history of El Niñ o, making it part of the natural evolutionary process of many populations and ecosystems. However, climate models suggest that the strength and frequency of El Niñ o events have recently increased, and will continue to do so, due to global warming [6]. In fact, the El Niñ o events of 1982-1983 and 1997-1998 were the strongest recorded in the last century and perhaps the last 400 years [7][8]. Thus, in the visible future, human activity may push environmental conditions to new extremes, and the capacity of populations to respond to these changes remains largely unknown.
The Galápagos archipelago ( Figure 1) lies in the primary region of ENSO activity, experiencing the large increases in SST and rainfall that are characteristic of El Niñ o events [9]. In the absence of ocean upwelling, food abundance is low, causing starvation of many marine organisms [10]. During the most recent, severe El Niñ os, population crashes of 77% (1982)(1983) and 65% (1997)(1998) were recorded in the endemic Galápagos penguin [11], while nearly 100% of Galápagos fur seal yearlings and large males perished [12]. Perhaps the best-studied example of the impact of ENSO on natural populations is the Galápagos marine iguana (Amblyrhynchus cristatus), a species that is only found in the archipelago. Marine iguanas are known to inhabit the shoreline of all the major islands in Galápagos, where they forage almost exclusively on algae from intertidal and (nearshore) subtidal zones.
The digestion of algae is made possible by a community of bacterial micro-symbionts that exist in the hindgut of the iguanas [13]. During severe El Niñ o events, when SST is elevated, the coastal environment is dominated by algae species that the iguanas cannot digest. On some islands, this has led to extreme undernourishment of iguanas, causing lowered body condition, increased stress hormone levels, significant changes in breeding behaviour, and mortality rates as high as 90% [14][15][16][17]. During La Niñ a, the other phase of ENSO, preferred algae are abundant and marine iguana populations have been shown to rebound quickly in terms of both census numbers and body condition [14].
Sharp population declines, or bottlenecks, like those seen during the recent El Niñ os, may translate into losses of genetic variation that can lead to increased rates of inbreeding and the fixation of deleterious alleles, and hinder the ability of populations to adapt to future changes in the environment [18][19][20]. In addition, there are a number of studies which have shown that decreased genetic diversity is associated with lowered immunocompetence, making populations more vulnerable to disease (see [21] for a recent example).
Although long-term demographic studies from several islands in Galápagos reveal the drastic effects of El Niñ o on marine iguana census numbers and reproductive behaviour, it is difficult to estimate the effective population size (N e ) from such data, and methods based on demographic information tend to overestimate N e . Conversely, indirect methods of N e estimation based on genetic data have been shown to be particularly informative, making them an important component in evaluating the consequences of potential bottleneck events such as El Niñ o [22].
Only a few studies have examined the impact of ENSO on genetic diversity [23][24][25], but there have been no attempts to use genetic data to quantify N e in natural populations that have experienced a severe El Niñ o event. When samples from two or more generations are available, N e estimates based on the temporal variance of allele frequencies [26][27] have been shown to be the most reliable method for detecting a genetic bottleneck in both simulated [22] and empirical datasets [28]. This is because levels of genetic drift are inversely proportional to population size, causing rapid changes in allele frequencies in small populations, and the magnitude of these changes between two time points indicates the severity of the bottleneck [29].
In this study, we estimated the harmonic mean of N e for marine iguana populations from 11 islands between 1991/1993 and 2004, a period that covers the severe 1997-1998 El Niñ o. Effective population size estimates were derived from the temporal variance in allele frequencies of 13 microsatellite loci representing more than 800 individuals collected before (1991/1993) and at least one generation after (2004) the warming event. In order to corroborate our results, we also analysed the complete mitochondrial control region for the same populations and most individuals for signs of temporal genetic changes. The large number of individuals and genetic loci used in our analysis, combined with the temporal sampling design, makes this one of the most comprehensive studies to date trying to measure the short term genetic effects of a severe El Niñ o event on natural populations.

Sampling
Approximately 800 marine iguanas from 11 islands in the Galápagos archipelago were sampled at two different time points before (in 1991/1993) and after (in 2004) the 1997-1998 El Niño event (see Figure 1 and Table S1 for a detailed description of taken samples). Sample sizes for each population and time point ranged from 20-58 individuals. We sampled mainly semi-adult and adult marine iguanas with sex classes evenly distributed across samples. Approximately 1-2 ml of blood were taken from the caudal vein of each individual and placed in storage buffer (100 mM Tris, 100 mM EDTA, 2% SDS). Total genomic DNA was extracted using the QIAamp 96 DNA Blood Kit (Qiagen).

Microsatellite analysis
Thirteen species-specific microsatellite loci were amplified in 806 individuals and scored for alleles as previously described [30][31]. We used the program NeESTIMATOR version 1.3 [32] and its implemented TM3 algorithm [33] to estimate effective population size (N e ) by measuring the temporal variance of microsatellite loci allele frequencies between time points using a likelihood-based approach. For the analysis of temporal changes in allele frequencies for N e estimation, the program NeESTIMATOR requires the user to define a reference time point of a population (generation 0) and at least one subsequent generation (in our case generation 1). Island populations sampled in 1991/1993 were entered as generation 0 and those sampled in 2004 as generation 1. Observed changes in allele frequencies and therefore estimates of N e only reflect changes that occurred from the reference point (generation 0 in 1991/1993) to the subsequent point (generation 1 in 2004). Since the TM3 program requires a maximum value for N e , we used the estimates of the maximum census size (N c ) of a respective island population [14] as the upper bound of N e , as in most published studies, N c exceeds or equals N e [29]. N e was then estimated as the maximum likelihood value of 10,000 updates within that range.
F ST differentiation between pre-and post-El Niñ o samplings was estimated with the program ARLEQUIN version 3.0 [34], including the test for significance (a = 0.05; see Table 1). Locus specific heterozygosity of samples from both time points was estimated with the program ARLEQUIN version 3.0 [34], and differences between time points for each island was tested with the non-parametric Wilcoxon signed-rank test using the Analyze-itH package for Microsoft ExcelH for significance (a = 0.05; see Table 1 and Table S2).

Mitochondrial analysis
Complete mitochondrial CR sequences (1183 base pairs) were generated for 838 marine iguanas (most of which were also analyzed for microsatellite loci). PCR was carried out on total genomic DNA using the primers IguanaCytb3 (59-ACCAGTA-GAACACCCMTTCATC-39) and 12s1984 [35]. PCR was performed for 35 cycles with an annealing temperature of 57uC and an extension time of 90 seconds. PCR products were purified using the Qiaquick PCR Purification kit (Qiagen) and sequenced

Island
Year of sampling Although the H-test was originally based on the number of segregating sites in a population (S), DNASP uses g in place of S. We chose to run the additional set of simulations with h because it should be less sensitive to low frequency migrants, which can greatly elevate the value of g, and give a false signal of lower than expected Hd.

Microsatellite analysis
The average heterozygosity of the 13 microsatellite loci calculated for each population and time point ranged from 0.637 for Pinta in 1993 to 0.815 for Fernandina in 2004. Island specific F ST values calculated between pre-and post-El Niñ o samples were not significantly different from zero for most comparisons, but was 0.008 and significant (p,0.05) for Marchena ( Table 1). The same level of F ST differentiation (0.008) was also found between time point samplings on Floreana, but in this case it was not significant at the 5% level. Based on the Wilcoxon signed-ranks test, locus specific heterozygosity was not significantly different (a = 0.05) between pre-and post-El Niño samples from a specific island, with the exception of Marchena and Santa Cruz (see p-values in Table 1 and the explicit locus specific comparisons in Table S2). The population on Marchena showed a significant decrease in locus specific heterozygosity (p = 0.032) from 1993 to 2004, while there was a significant increase (p = 0.04) on Santa Cruz from 1991 to 2004. Figure 2a shows the Bayesian estimates of N e and its associated confidence intervals (CI) obtained from microsatellite data using the temporal variance method. N e values for each population represent the harmonic mean of the effective population size over the time interval between pre-and post-El Niñ o samplings. A low N e value (,50 individuals) was only found for the population on  (Fig. 2a and Figure S1; see Material and Methods for details). Figure 2B depicts the low N e estimate of around 40 individuals (CI of 21-86 individuals) between 1993 and 2004 experienced by the population on Marchena (for comparison, see Figure S1 for plots of other populations). The large confidence intervals associated with N e estimates on the majority of islands shows that our study approach, given the number of individuals and microsatellite loci used, is not able to provide consistent estimates for large values of N e . This indicates that the effective population size must have remained fairly high on these islands during the time interval between samplings and that N e did not fall below the critical size of 50 individuals.

Mitochondrial control region analysis
Mitochondrial control region (CR) data was generated for nearly the same set of individuals and populations used in the microsatellite analysis (see Table S1 for details). Among all the islands, the population from Marchena exhibited the largest decrease in haplotype diversity (Hd) and nucleotide diversity (p)   Table S3). Although less severe, a loss of mitochondrial variability was also detected on Genovesa (DHd = 20.232, t = 2.11, upper tail p,0.025; see Table S3). All remaining populations did not exhibit significant changes in levels of mtDNA diversity (see Table S3 and

Study feasibility
The temporal sampling of genetic data can be a powerful means of estimating N e between two closely situated time points, as is the case in our study. However, there are certain requirements of such methods, particularly the temporal variance approach, which must be met in order to achieve valid N e estimates and a significant statistical power. Simulations based on the temporal variance method have shown that a bottleneck with only 34 surviving individuals can be detected with a probability of higher than 80% after one generation when using 12 moderately polymorphic microsatellite loci (with mean of 4.3 alleles per locus) and at least 30 individuals per population [22]. Since loci with higher levels of polymorphism increase statistical power [22], the use of 13 highly polymorphic microsatellite loci (with a mean of 8 alleles per locus; see Table 1) in our study improves our ability to detect a critical reduction of N e beyond that of the simulation.
In addition, the temporal variance approach, as well as other tests comparing populations before and after a population reduction (e.g. tests for decreases in heterozygosity, allelic diversity, haplotype and nucleotide diversity, etc.), require that a population be sampled before and at least one generation after the potential bottleneck event. In our study, we collected marine iguana specimens from eleven different islands in 1991 or 1993, before the 1997-1998 El Niñ o, and in 2004, approximately six years after the end of the warming period. Since we sampled mainly adults and sub-adults in 2004, along with some juveniles, it is possible that some of these individuals were born prior to the El Niñ o period. However, if a given population was reduced to a critically small size (N e #50) during the 1997-1998 warming period, it is very likely that iguanas sampled in 2004 would primarily be individuals born in the first generation after the bottleneck event and would therefore reflect the gene pool of the reduced population. On the other hand, if the population size remained fairly high during the warming period, the gene pool reflected in 2004 would not deviate significantly from its pre-El Niñ o composition. In the first scenario (N e #50), given the number and variability of microsatellite loci employed in this study, we would be able to detect a critical reduction of N e with a rather high probability as pointed out above. In the second case, however, we would not expect to see any major changes in allele frequencies between time points, since the population makeup had not been significantly altered.
An additional concern is whether the six years following the 1997-1998 El Niñ o were sufficient for the development of a new generation into adult and sub-adult classes. Female marine iguanas usually reach the age of first reproduction in their fourth or fifth year; however, data from the severe 1982-1983 El Niñ o show that the age of first reproduction of female marine iguanas may be three years or less following ENSO-induced population crashes [16]. This is partially a result of rapid growth rates that are common during La Niñ a periods following El Niñ o [16]. Males, on the other hand, may not gain access to mates until they are at least 12-15 years old, due to strong sexual selection for largebodied individuals [39][40]. However, if a population is far below carrying capacity, which is often the case during strong El Niño events, males can reach maximum size in less than six years, as was found on the island of Santa Fe after the 1982-1983 El Niñ o [15]. Moreover, large, territorial males appear to suffer the highest rate of mortality during El Niñ o [15], increasing the likelihood for younger males to occupy mating territories and achieve mating success after strong El Niñ os. Thus, it is plausible that a new generation could have been established on many islands in the six years following the 1997-1998 El Niño.

The genetic impact of the 1997-1998 El Niñ o on Galá pagos marine iguanas
There have only been a few attempts to measure the impact of severe El Niñ o events on genetic diversity in natural populations, with mixed results. Randomly amplified polymorphic DNA (RAPD) markers have shown that heterozygosity in intertidal kelp (Lessonia nigrescens) populations along the Chilean coast fell by 50% compared with unaffected populations following the 1982-1983 El Niñ o [23]. Also, Galápagos penguins (S. mendiculus), whose populations have been severely diminished by ENSO, have significantly lower observed heterozygosity compared to the common and unaffected Magellanic penguin (Spheniscus magellanicus) from South America, possibly the outcome of a series of strong El Niñ o events [24]. Conversely, a study on butterflies (Arhopala epimuta) in Borneo found that El Niñ o-induced forest fires did not alter the temporal genetic structure of microsatellite loci (5 loci) or mtDNA (control region) in populations that were sampled before and after the 1997-1998 El Niñ o [25].
In this study, we focus on the impact of ENSO on multiple populations of Galápagos marine iguanas, which are also known to experience large population reductions (up to 90%) during intense El Niñ o periods [14]. Although high levels of mortality were recorded on a number of different islands during the 1997-1998 El Niñ o, the microsatellite data indicate that only one of the 11 populations examined (Marchena) showed signatures of either a critically small population size during the time interval in which the El Niñ o event occurred or a reduction in genetic diversity following the warming period. Based on the temporal variance approach, there is strong support for a N e below 50 on Marchena during the period between 1993 and 2004 (Figures 2a and 2b), indicating that the population must have become quite small during this time. F ST values, which are also based on differences in allele frequencies between two samples, support this finding, since a statistically significant level of differentiation between time points was only found for Marchena (Table 1). Lastly, only the population from Marchena exhibited a significant decrease in locus-specific heterozygosity values between the first and second samplings (Table 1 and Table S2).
The results of the mtDNA data were similar, where Marchena exhibited the largest decrease in haplotype diversity and nucleotide diversity from 1993 to 2004 (Figure 2a; Table 2). However, unlike the microsatellite data, the population from Genovesa also showed a significant decrease in Hd between samplings indicating a possible population reduction. Additionally, coalescent simulations showed that observed Hd values were significantly lower than expected under neutral evolution for both Marchena and Genovesa in 2004 but not in the first sampling predating the El Niño (see Table 2).
The genetic signature of a population decline detected between 1991 and 2004 on Genovesa likely reflects the population crash, from 15,000 to 900 individuals, that was observed during the 1997-1998 El Niño event, and there is no evidence for other incidents which may have drastically effected the demography of this island during this time period [14 and Martin Wikelski pers. comm.]. Conversely, for Marchena, marine iguanas may also have been highly impacted by a strong volcanic eruption in 1991, where lava flows continued for at least 40 days. This eruption is likely to have caused mortality of marine iguanas in both the terrestrial environment (as shown for Galápagos tortoises on the island of Isabela [41]) and the marine environment, where water was described as ''too hot to touch'' and was observed to have killed many fish and other aquatic organisms [42]. Although the actual eruption occurred two years before the initial sampling on Marchena in 1993, a new generation is not likely to have emerged until some point between 1993 and 2004, causing its genetic effects to be confounded with that of the 1997-1998 El Niño. In this case we can imagine three scenarios where populations were reduced due to the eruption, to El Niño, or the combined effect of both. The last is an intriguing option given the magnitude of the population bottleneck on Marchena compared to the other populations.
Additionally, there are two other forces which may have influenced marine iguana populations between 1991/1993 and 2004: the long term, but weaker El Niño period extending from 1990-1995 [1,11], and an oil spill that occurred off the island of San Cristóbal in 2001 [43]. The former cannot be ruled out, but the genetic impact of the oil spill could not be reflected in our data since the three years between the spill and the 2004 sampling are not sufficient for the establishment of a new generation of marine iguanas, which is an important requirement of the temporal variance method. Thus, this does not rule out the possibility that the oil spill impacted the effective population size of the iguanas, but rather that, with our study design, we are simply not able to detect it.
Under random mating, heterozygosity is lost from a population at a rate of 1/2N e (1/N ef in the case of mtDNA) per generation. Thus, genetic diversity is lost more rapidly in small populations. Yet, during brief bottleneck events, only a small amount of the total heterozygosity is lost, even after a considerable reduction in population size (e.g. 2% after the first generation of a bottleneck of N e = 50) [44]. However, the loss of diversity follows an exponential decay process, where long-term decreases in the number of individuals can have a detrimental impact on N e . Therefore, it is not surprising that, with the exception of Marchena, a single, strong El Niñ o event did not stimulate a loss of heterozygosity in the populations, even if mortality rates were high. The temporal variance approach, on the other hand, relies on shifts in allele frequencies rather than loss of diversity, and has been shown to be particularly sensitive in estimating N e in small populations [28,45]. Since, as discussed above, our study design offers sufficient statistical power to detect a bottleneck, the large confidence intervals surrounding the majority of N e estimates suggest that the effective size of most populations of marine iguanas were so large that they underwent only slight shifts in allele frequencies. Such temporal stability has been observed in a number of other studies [46,47] and is quite plausible given the large census size estimates of marine iguanas recorded of many islands [14].
Recommendations for future study design in assessing the genetic impact of El Niñ o Although we did not find a strong influence of El Niño on genetic diversity in marine iguanas, the data show that the genetic impact of a single, intense environmental challenge may vary even within a single species and depends on the specific history of a population. The low N e estimate for Marchena, calculated over an interval in which a severe El Niño event and a volcanic eruption occurred, suggests that it is important to consider the potential synergistic relationship between El Niño and other phenomena. While we saw little evidence of a bottleneck in the majority of our populations, future experimental designs must account for multiple natural (e. g. long-term, weaker El Niño periods) and human-induced (e.g. oil spill) disturbances, or else the causes of population reductions may be confounded.
Evolutionary theory predicts that the long-term effective size of a population is the harmonic mean of N e over many generations [44]. Under this scenario, long-term N e is heavily influenced by generations in which the effective population size is small. Therefore, the high mortality rates of marine iguanas associated with El Niño warming may compound over multiple ENSO cycles. Since climate models predict that the strength and frequency of El Niño events will continue to increase, a true understanding of the long-term impact of ENSO on population persistence may only come from experiments designed to measure the genetic impact of this phenomenon following a series of consecutive events. Table S1 Sampling localities by island and sample sizes for 13 microsatellite loci and mitochondrial control region sequences (1183 bp) for the two temporal samplings (1991/1993 and 2004). The first column lists sampling localities by island, specific sampling location (in parentheses), and geographical coordinates. Sample sizes for 13 microsatellite loci and mitochondrial control region sequences (1183 bp) are reported in separate columns for the two temporal samplings before the 1997-1998 El Niño (in 1991 or 1993, or both years for Santa Fé) and after the 1997-1998 El Niño in the year 2004. Found at: doi:10.1371/journal.pone.0001285.s001 (0.04 MB DOC)

Table S2
Locus specific heterozygosity values for 13 microsatellite loci for pre-and post-El Niño island samplings. Locus specific heterozygosities calculated using the program ARLEQUIN for marine iguana populations sampled before the 1997-1998 El Niño (in 1991 or in 1993 or as for Santa Fé in both years) and after the 1997-1998 El Niño in the year 2004. The first column (locus) shows names of microsatellite loci. The other columns report the heterozygosity for a given population for each time point. Differences in heterozygosity between time points were tested with a Wilcoxon signed ranks test and associated p-values are provided in the last line. Significant p-values (p,0.05) are marked with an asterik (*). The only significant decrease of locus specific heterozygosity from 1993 to 2004 was found for the population on Marchena, whereas Santa Cruz showed a significant increase during this period.  Figure S1 Graphical display of the TM3 results (estimates of the effective population size [Ne] generated from the temporal variance of allele frequencies using a likelihood-based approach) as provided by the program NeESTIMATOR. Graphical display of the TM3 results (estimates of the effective population size (Ne) generated from the temporal variance of allele frequencies using a likelihood-based approach) as provided by the program NeESTI-MATOR (see Methods). Point estimates for different values of Ne are shown (x-axis) with their corresponding log likelihood values (y-axis). For each island, the estimated Ne with the highest log likelihood, the calculated 95% confidence interval (CI), and the upper bound of the estimate are provided below the graphic. Found at: doi:10.1371/journal.pone.0001285.s005 (2.03 MB DOC)