Rapid selection against inbreeding in a wild population of a rare frog

Populations that are small and isolated can be threatened through loss of fitness due to inbreeding. Nevertheless, an increased frequency of recessive homozygotes could increase the efficiency of selection against deleterious mutants, thus reducing inbreeding depression. In wild populations, observations of evolutionary changes determined by selection against inbreeding are few. We used microsatellite DNA markers to compare the genetic features of tadpoles immediately after hatch with those of metamorphosing froglets belonging to the same cohort in a small, isolated population of the threatened frog Rana latastei. Within a generation, the inbreeding coefficient (FIS) decreased: at hatch, FIS was significantly >0, whereas FIS was <0 after metamorphosis. Furthermore, heterozygosity increased and allelic frequencies changed over time, resulting in the loss of genotypes at metamorphosis that were present in hatchlings. One microsatellite locus exhibited atypically large FST values, suggesting it might be linked to a locus under selection. These results support the hypothesis that strong selection against the most inbred genotypes occurred among early life-history stages in our population. Selective forces can promote changes that can affect population dynamics and should be considered in conservation planning.


Introduction
Small, isolated populations are particularly subject to inbreeding, and are thus intrinsically vulnerable to extirpation (Keller and Waller 2002). Both laboratory and field studies have shown that increased inbreeding associated with small and isolated populations can strongly reduce both the variance and mean of traits that are important components of fitness, including survival, fertility and individual growth rate (Saccheri et al. 1998;Keller and Waller 2002;Rowe and Beebee 2003;Frankham 2005;Wright et al. 2008). The reduced survival of inbred individuals can further reduce population size and increase susceptibility to stochastic processes, such that mutual reinforcement among environmental, demographic and genetic stochasticity creates an 'extinction vortex' (Gilpin and Soulé 1986). However, if enough genetic variability is retained by the population, less inbred individuals, or individuals not carrying deleterious alleles can be favoured through selection. If inbreeding depression is mainly caused by the effect of strongly deleterious, recessive alleles, natural selection may be particularly effective against them when they occur in the homozygous state (Wang 2000). Two processes can thus occur in small populations. First, the loss of inbred individuals can cause a further reduction of population size, leading to cycles of exacerbated inbreeding and drift in subsequent generations (referred to as 'mutational meltdown'; Lynch et al. 1995). On the other hand, selection against inbred genotypes can purge inbreeding depression (Hedrick 2001). In this case, populations may be able to recover demographically despite initial demographic losses by removing inbred individuals and reducing the frequency of deleterious alleles. Selection can be most effective at countering inbreeding depression under particular conditions, such as when the fitness disadvantages of deleterious mutations are particularly strong, inbreeding does not occur randomly, and migrants are not (re)introducing deleterious alleles (Glemin 2003;Edmands 2007).
Detecting the ongoing population genetic processes that lead to mutational meltdown or recovery can provide key insights for the management of threatened species at greater risk of inbreeding. However, detecting either of these patterns in wild animal populations can be challenging. Under natural settings, measuring the relative performance of individuals with respect to their inbreeding status usually requires long-term capture/mark/recapture analyses over multiple generations and the reconstruction of pedigrees. For this and other reasons, the measurement of inbreeding depression in the wild has been done primarily in larger vertebrates where parentage is easier to determine and the number of offspring per female is small (Keller and Waller 2002;Kruuk and Clutton-Brock 2008). For taxa where such parameters are not so amenable to measurement, two other approaches are used. Inbreeding depression can be measured under laboratory conditions, however inbreeding depression is frequently stronger under natural (stressful) conditions, and results obtained by laboratory studies may not reflect the complexity of natural conditions (Bijlsma et al. 1999;Joron and Brakefield 2003;Halverson et al. 2006). Alternatively, averaged population responses (e.g., F IS ) can be used, by comparing populations with different inbreeding histories rather than measuring fitness of individuals (Saccheri et al. 1998).
Molecular genetics can be a powerful approach for the study of microevolutionary changes over a range of temporal scales, from a few generations to thousands years. In anuran amphibians, typically 80-95% of mortality occurs between hatch and metamorphosis (Vonesh and De la Cruz 2002) and inbreeding depression may be most readily detected through the measurement of relevant population genetics parameters during these periods of heavy mortality (Rowe and Beebee 2003;Halverson et al. 2006;Chapman et al. 2009). In this study, we used microsatellite DNA polymorphisms to study temporal variation of genetic features in small, isolated population of a threatened frog suffering strong inbreeding depression (Ficetola et al. 2007). We analyzed the genetic features of the same cohort at two subsequent life-history stages, immediately after hatch and at metamorphosis; our sampling protocol was therefore designed to compare the genetic features of a cohort before and after the effects of potentially strong natural selection. We hypothesized that strong selection against the most inbred genotypes might cause strong shifts in genetic features of the population. Evolutionary signatures of selection against inbreeding include an increase in heterozygosity, a decrease in the inbreeding coefficient F IS , and loci exhibit-ing outlying values of F ST (Beaumont 2005). Moreover, the reduced frequency or even elimination of specific alleles all may occur if selection is removing deleterious alleles expressed in the homozygous condition.

Study species and area
The Italian agile frog, Rana latastei Boulenger, is endemic to Northern Italy and adjacent countries and is classified as vulnerable to extinction by the IUCN (Andreone et al. 2008). Habitat requirements are river washes, forest brooks and ponds associated with lowland forest cover, primarily along the Po River drainage in Italy. In this area, habitat has been subject to intense urbanization and agriculture. In particular, the western part of the distribution of R. latastei is broken into a limited number of small and isolated forest fragments. Many breeding sites are characterized by small census size and geographic isolation (Ficetola et al. 2007) and this combination of isolation and small population size has increased the risk of extinction of populations (Ficetola and De Bernardi 2004). Further, the part of the species range from which our study population comes from is characterized by a 3· reduction in genetic variability when compared to the eastern part of the range (Garner et al. 2004) and experiments have shown that reductions in genetic variability correlate well with reduced fitness (Pearman and Garner 2005;Ficetola et al. 2007). Breeding occurs in late winter; females lay a single egg mass per year (Barbieri and Mazzotti 2006), therefore the number of clutches represents both annual population reproductive output and the number of breeding females in a population. Multiple paternity has never been observed in R. latastei. Some studies have described multiple paternity in related species, but the offspring of polyandrous mating is a small proportion. For instance, in Rana temporaria, clutch piracy affected 84% of clutches, but pirate males fertilized only 26% of pirated clutches, and in these clutches 24% of embryos were fertilized by pirates, i.e., 5% of all eggs were sired by a secondary male (Vieites et al. 2004). Similarly, in Rana dalmatina only 4% of all eggs are sired by a secondary male (Lodé et al. 2005). Our results are robust even with the removal of 5% of larvae from the analyses, therefore the occurrence of multiple paternity at rates similar to the ones observed in related frogs would have minimal influence on the outcome of this study (see also Ficetola and De Bernardi 2009;Ficetola et al. 2010 for further details).
We studied one population located close to the Alserio Lake, in Northern Italy (population AL in the map in Ficetola et al. 2007). At this location, we detected a total of 43 clutches during the 2003 breeding season, most of which were found in a single pond (hereafter: study pond; 32 egg masses in 2003). The study pond is small (surface area 60 m 2 ), shallow (maximum water depth: 50 cm) and with clear water; clutches are deposited onto submerged tree roots and sticks below the water surface, thus, they are easily detected. Clutch counts therefore provide a reliable standard census method for agile frogs (Salvidio 2009). A few clutches were found in a small satellite wetland within 50 m of the study pond; clutches from the satellite pond were not considered here. The study population is isolated: the next available breeding location was located 3 km away, and the intervening distance was composed of unsuitable habitat, lacking both suitable breeding locations and contiguous forested habitat; a motorway acts as a further barrier.

Sampling protocol and analyses
We subsampled all 32 viable clutches deposited in the study pond in March 2003, taking a single tadpole per egg mass for DNA extraction. This sample therefore represented the entire hatching cohort. In late June-early July, we returned to the study pond and we sampled the pond banks to collect metamorphic froglets (stage 39-42; Gosner 1960). We captured active individuals; we also used a net to collect leaf litter and cryptic individuals; overall, we obtained 29 froglets. We collected a small tissue sample (remnant tail tip) from each metamorph for DNA extraction. We restricted our sampling to these semiaquatic stages to prevent sampling any possible metamorphs that may have come from the adjacent woodland.
Overall, 78 microsatellite loci developed either for R. latastei or for other Rana species have been tested on R. latastei samples. Thirty-one of these loci were monomorphic, and 41 did not amplify or showed presence of null alleles, therefore only six polymorphic loci are actually available for this species, a proportion significantly lower than in other frog species (Garner and Tomio 2001;Primmer and Merilä 2002; T.W.J.G. unpublished data); this is probably caused by the colonization history of this small range species combined with the extreme landscape fragmentation (Garner et al. 2004;Ficetola et al. 2007). DNA was extracted using a Chelex-based protocol and amplified at variable microsatellites following Garner and Tomio (2001). [ 33 P]-labeled fragments were electrophoresed through 6% polyacrylamide gels and fragments were detected using autoradiography. Alleles were visually scored against a pGemÒ-3zf+ (Promega, Madison, Wisconsin, USA) sequence reference marker. Further details on genetic analyses are provided in Ficetola et al. (2007).
We used Genepop 3.4 (Raymond and Rousset 1995) to test for linkage disequilibrium among the microsatellite loci that proved to be variable, and to test for Hardy-Weinberg equilibrium in larvae and metamorphs. To compare the genetic variability of hatchlings and of metamorphs, we estimated the percentage of polymorphic loci, allelic richness, expected and observed heterozygosity for both developmental stages. We also tested for genetic differences among developmental stages by calculating respective F ST values, and permuting these values 10 000 times to test for significance.
We calculated F IS for hatchlings and for metamorphs as a measure of the level of inbreeding in a cohort (Wright 1969); we estimated 95% confidence intervals for each coefficient by bootstrapping 10 000 times, and tested for significant differences in the inbreeding coefficient among developmental stages. We estimated F IS separately for hatchlings and metamorphs, therefore differences in F IS might be in principle due to differences in the reference population considered for F IS estimation. To address this, we used simulations to evaluate whether random mortality or sampling, and changes in the reference population, may have caused the observed changes in F IS . We simulated 1000 individuals with the same F IS of hatchlings and the same allele frequencies detected in the hatchlings. Subsequently, we randomly selected 29 individuals from the simulated population to form the metamorphic sample, and we estimated F IS for these individuals. This process was repeated 1 000 000 times.
We used Fdist2 (Beaumont and Nichols 1996;Beaumont 2000) to identify loci showing the potential signature of selection. This approach is based on the principle that loci under selection or tightly linked to a locus under selection should exhibit outlying F ST values (outliers) (Beaumont 2005). We generated a null distribution based on 30 000 simulated loci; the mean F ST was calculated on the basis of empirical distribution and then iteratively modified as outliers were identified and removed (Beaumont and Nichols 1996). To evaluate the robustness of our results to different approaches in the estimation of expected F ST , this analysis was repeated by calculating F ST without the iterative removal of outliers, which yields more conservative results. This second analysis produced very similar results (not shown). Results from Fdist2 may be misleading when analyzing multiple populations grouped in phylogenetic hierarchies (Latta 2008). For analysis with Fdist2, we treated the larval cohort as the first 'population', and the metamorphic cohort as the second 'population'. Therefore, phylogenetic hierarchies are absent from our data.
Finally, we used a Correspondence Analysis (CA) to detect differences in genotype distributions between tadpoles and metamorphs by graphically projecting all the individuals on the space defined by factors extracted on Rapid selection against inbreeding Ficetola et al. the basis of the similarity of their allelic states. This last analysis was completed using Genetix 4.05 (Belkir 2004).

Results
Three loci (RlatCa17, RlatCa21 and Rt2Ca9) were monomorphic across the entire sample: therefore all subsequent analyses were performed using the remaining three polymorphic loci (RlatCa9, RlatCa27 and RlatCa41). We did not detect evidence of linkage disequilibrium among any of these loci (all P > 0.08). In the earlier (larval) lifehistory stage, analysis of Hardy-Weinberg equilibrium showed a nonsignificant trend for heterozygote deficiency (P = 0.07). When loci were examined individually, only locus RlatCa27 exhibited a significant heterozygote deficit in larvae (P = 0.004). Conversely, in the later, postmetamorphic life-history stage we observed a significant global heterozygote excess (P = 0.013), despite the lack of any effect when any one locus was analyzed independently (all P > 0.1). Percentage of polymorphic loci, allelic richness and expected heterozygosity were similar in both stages. Observed heterozygosity of the cohort increased from 0.404 at hatch to 0.547 in metamorphs; nevertheless, standard errors were large (Table 1). At hatch, the inbreeding coefficient was significantly >0 (F IS = 0.183, 95% CI: 0.027 to 0.305). In contrast, the inbreeding coefficient of metamorphic individuals was negative (F IS = )0.169, 95% CI: )0.027 to )0.337, Fig. 1) and significantly different from the hatch value (P < 0.001). This decrease in the inbreeding coefficient over life-history stages was strongly influenced by one locus (RlatCa27) but was evident at all three of the variable loci (Fig. 1). Exclusion of locus RlatCa27 from the analysis resulted in a marginally non-significant trend (P = 0.057). Simulations showed that the observed changes in F IS could not be explained by random mortality/sampling, as the observed F IS of metamorphs was lower than the first percentile of the distribution of simulated data (Fig. 2). Genotype frequencies differed significantly among lifehistory stages (F ST = 0.101, P < 0.001), but pairwise F ST was not significant if locus RlatCa27 was removed from the analysis (F ST = 0.010, P = 0.169). Locus RlatCa27 exhibited an atypically large F ST value among life-history stages (P = 0.001, Fig. 3). This result is particularly striking because Fdist2 is expected to have relatively low power for detecting outliers when the number of loci included in the analysis is not large (Beaumont and Nichols 1996;Nielsen et al. 2006).
Both the significant decrease in the inbreeding coefficient and change in genotype distribution detected among life-history stages was primarily driven by the disappearance of certain homozygotes in the metamorphic sample that were present in the hatch data set. Specifically, 38% of hatchlings were homozygous for the same allele (154 bp) at locus RlatCa27, while none of the metamorphosing froglets were homozygous for this allele at this locus (Fig. 4). The frequency of this allele dropped from 0.518 in larvae to 0.054 in the metamorphs. One further homozygote was eliminated from the genotype pool as development progressed (166 bp) and no new homozygotes occurred in the metamorphic data pool that were not previously detected in the earlier life-history stage (Fig. 4). The first component extracted from the CA explained 24.3% of the total variation of genotypes and showed clearly how the genotypes of metamorphs constituted a small subset of the hatchling genotypes (Fig. 5). Hatchling scores along the first CA component ranged from +2 to )1, whereas all metamorph scores were <0.5. Genotype variability along the first component significantly decreased from hatch to metamorphosis (Levene's test for homogeneity of variances: F 1,59 = 29.0, P < 0.001, Fig. 5).

Discussion
Our study showed that within a cohort of R. latastei F IS decreases over time while observed heterozygosity increased slightly. One locus was predominantly responsible for atypically large differences among life-history stages (F ST ) that are considered to be evidence of selection at a locus (Beaumont and Nichols 1996). Further, only a subset of the genotypes present among clutches persisted through to metamorphosis and the frequency of a common allele, present only in hatchlings in the homozygous state, was strongly reduced in the metamorphic cohort. We present arguments below for why we believe these observed patterns are indicative of strong selection against inbreeding.
Microsatellites are usually considered to be neutral loci. Inferring selective processes from changes in neutral markers requires the elimination of the alternative explanations, such as immigration, drift, nonrandom mating and mutation. In our study, nonrandom mating and mutation can be ignored, as we sampled a single, nonreproductive cohort. Immigration had no effect on our sample because of pond isolation and the fact we sampled before metamorphs from nearby wetlands could have emigrated. Stochasticity may have a role in some of the changes observed; for instance, the frequency of some but not all the homozygotes decreased in metamorphs (Fig. 4). However, stochasticity cannot be the sole cause for the observed changes, as simulations showed that it is highly unlikely that the strong decrease of F IS we have observed could be attributed to stochastic processes (Fig. 2). Thus, we conclude that natural selection was the most likely driver of the temporal variation in genetic frequencies we observed; the strong decrease of F IS suggests Rapid selection against inbreeding Ficetola et al. that selection against inbreeding is the most likely explanation of the observed genetic changes. Nevertheless, it should be remarked only locus RlatCa27 was not in Hardy-Weinberg equilibrium in larvae; other possibilities therefore exist for the pattern at that locus. For example, allele 154 might be associated with mate choice in adults, causing an excess of homozygotes in eggs, then selected against in larvae. We are not the first to observe a relationship between measures of genetic diversity at microsatellite loci and measures of fitness (Chapman et al. 2009). Furthermore, previous research has identified specific patterns in the correlation between marker variability and fitness that are matched by our data set. For example, in red deer, reed warblers, harbor seals and sea lions, heterozygosity at multilocus microsatellite genotypes was strongly correlated with components of fitness, but in all these studies single microsatellite loci were shown to have a significantly greater contribution to interindividual variability and locus-specific correlations with fitness components (Pemberton et al. 1999;Hansson et al. 2004;Acevedo-Whitehouse et al. 2006). Balloux et al. (2004) showed that, in the absence of strong linkage disequilibrium and complex mating systems, it should be necessary to sample large numbers of microsatellite loci to detect a true fitness/heterozygosity correlation. This has not prevented the relationship between diversity estimated at microsatellites and fitness from being detected frequently in studies using a handful of markers. In our study species, six loci were sufficient for detecting significant variation in genetic diversity across the species range (Garner et al. 2004) and this variation in diversity correlated strongly with fitness measures (Pearman and Garner 2005;Ficetola et al. 2007). The fact that we detected an effect with only half that many loci reflects the pattern of genetic variability inherent to the western part of the species range. Fixation for three of six published loci used in this study was previously detected in the Lombardy Region because of the joint effect of postglacial colonization and human-driven habitat fragmentation (Garner et al. 2004;Ficetola et al. 2007); during the course of PCR primer development and crossamplification with primers developed for other Rana species, another 31 primer sets were characterized as monomorphic in R. latastei (Garner and Tomio 2001;Primmer and Merilä 2002). Thus, based on a random sample of dinucleotide repeat regions that amplify reliably, the true number of polymorphic loci may be as low as 3 of 37.
F IS of larvae was significantly higher than zero (Fig. 1), which at first glances is puzzling, as F IS is a measure of inbreeding as nonrandom mating (Keller and Waller 2002), and we considered a single population breeding in one small pond, where it is difficult to invoke substruc-turing (Wahlund effect). Nevertheless, large F IS values are often observed in small, isolated populations of anuran amphibians (e.g., Hitchings and Beebee 1997;Andersen et al. 2004). Increased levels of polygyny are commonly reported for anurans, and have also been described for the study population (Ficetola et al. 2010); polygyny can result in inbreeding coefficients greater than expected under random mating; this effect can be exacerbated in small populations (Balloux et al. 2004). No matter the process originating positive F IS in larvae, the change in F IS was more striking than its absolute values: F IS strongly decreased during larval development, indicating than the less inbred genotypes were favoured by selection.
Linkage between microsatellites and loci under selection can be particularly strong in small, inbred populations carrying a heavy genetic load and after bottlenecks, as in these populations larger portions of the genome are expected to be under balancing selection, and probability is greatest for loci linked to overdominant loci to remain polymorphic after bottleneck events (Ohta 1982;Balloux et al. 2004;Van Oosterhout et al. 2004;Zartman et al. 2006). The few microsatellites retaining heterozygosity can therefore constitute a nonrandom sample of the genome (Santucci et al. 2007). The study population corresponds well to this situation: it is small, isolated and passed through a demographic bottleneck some years before the season in which we sampled for this study (Ficetola and De Bernardi 2005). From this, the relative importance of three loci in our species for describing genetic variability is likely greater than in other, more genetically variable species. A limited number of variable microsatellites is typical of endangered species suffering loss of genetic diversity (Gebremedhin et al. 2009a). These species are often the ones for which studies are more urgent; lack of variable microsatellites should not hamper performing analyses that can provide important management information (Richards-Zawacki 2009).
There is debate if the relationship between inbreeding and fitness is due to the effects of many loci of small effect or the effects of a few loci, or even one, of marked effect (Chapman et al. 2009). The limited panel of available markers available for this study does not allow us to distinguish between a multilocus or single locus model. Even so, although all our polymorphic loci contribute to the observed changes in F IS and heterozygosity, the strong effect of locus RlatCa27 suggests that certain allele(s) at this locus may be linked to lethal or semi-lethal allele(s) at a locus involved in early development. The observation that not all loci decreased their frequency in the homozygote state (e.g., 162 bp: Fig. 4) supports the hypothesis that not all of the microsatellite allele(s) is/are associated with deleterious mutations. Such linkage patterns between microsatellites and genes exposed to strong selection have been inferred and even detected several times Acevedo-Whitehouse et al. 2006;Nielsen et al. 2006;Mäkinen et al. 2008;Amos and Acevedo-Whitehouse 2009), and might be more detectable in small, isolated populations (Balloux et al. 2004;Van Oosterhout et al. 2004).
An increased frequency of recessive homozygotes through inbreeding is expected to increase the efficiency of natural selection against deleterious mutants and reduce inbreeding depression through purging (Crow 1970;Boakes and Wang 2005). If inbreeding depression is mainly caused by deleterious mutations having a strong effect, selection may rapidly eliminate deleterious alleles, and population fitness rebounds (Wang et al. 1999;Wang 2000) while the inbreeding levels of populations are reduced (Moore 2007). However, whether purging actually occurs to an extent that enables the recovery of populations has been strongly debated, primarily due to the equivocal nature of both field and experimental evidence (Leberg and Firmin 2008;Chapman et al. 2009). The disappearance of one common allele occurring primarily in the homozygous state among life-history stages (154 bp at locus RlatCa27) fits the signature of strong selection against inbreeding. Theoretical models suggest that this type of selection can be particularly effective in species with large reproductive ability (Wang 2000). In many anuran species, including R. latastei, fecundity is on the order of thousands of eggs. Moreover, larval mortality is high and density dependent, creating ideal conditions for selection to take effect early during development. If a sufficient number of females reproduce, the death of a large proportion of tadpoles may not threaten population persistence in the short term (Vonesh and De la Cruz 2002). This suggests that purging could be more effective in anuran populations than in taxa with limited individual reproductive effort, such as larger vertebrates (Jamieson et al. 2006;Chapman et al. 2009).
Under these conditions, selection may affect population dynamics by reducing inbreeding depression and the risk of an extinction vortex. Relationships between evolutionary and population dynamics can have important consequences for management and in population viability analyses, and when alternative options (e.g., translocations versus in situ management) are compared (Ficetola and De Bernardi 2005;Ficetola et al. 2007). For instance, actions allowing for a rapid demographic rebound (Wang et al. 1999), such as habitat reclamation, can be particularly effective for the recovery of populations. Our knowledge of the dynamics of inbreeding depression in the wild is far from complete, however, and more studies are required to improve decision making that account for the genetic responses that may or may not be associated with population recovery. Long-term studies of amphibian populations are increasingly available, and will be necessary to evaluate the relative importance of selection and inbreeding in determining population dynamics, and whether selection is strong enough to rescue isolated populations from inbreeding depression (Bensch et al. 2006). Our analysis considered only one population. Populations often experience different histories, conditions and selective forces; analyses including a larger number of populations are better for describing general patterns (Halverson et al. 2006).
To date, many studies on temporal variation of genetic features in wild, nonmodel species have been performed over long temporal scales. Our analysis showed a striking change in genetic features during a very short time, which might affect population dynamics. The fate of the population will be therefore determined not only by the interaction between environment, demographic processes and inbreeding depression. Selective forces can act over the functional variability among individuals, promoting changes influencing the future of endangered species, and should be considered in the analysis of population dynamics (Bensch et al. 2006;Saccheri and Hanski 2006). The short time scale over which these forces can operate may increase the complexity of interactions between environment, demography and genetics; molecular genetics can help to elucidate these processes in the wild, also at short temporal scales. The panel of tools available for conservation genetics is quickly broadening. The possibility to detect the signature of selection in natural populations using both neutral and non-neutral markers may increase our understanding of evolutionary processes in threatened populations (Gebremedhin et al. 2009b), providing a widening range of management options for conservation.