Ecological causes of fluctuating natural selection on habitat choice in an amphibian.

We estimated natural selection targeting three traits related to habitat choice in a frog (Pseudacris maculata) breeding in pools on the rocky shores of Isle Royale, Michigan, over 16 years. Our aim was to identify the form and ecological causes of annual variation in directional and correlational selection as expressed in the survival and growth of tadpoles. We found directional selection favoring early breeding, but pool choice was under weak stabilizing selection. However, the form of stabilizing selection and the position of the optimum trait value shifted among years with the severity of disturbance and the intensity of biotic interactions. In years when wave wash and pool desiccation were severe, selection shifted to favor tadpoles in habitats where these risks were less pronounced. If predatory dragonfly larvae were abundant, selection favored tadpoles in small pools where dragonflies did not occur. When intraspecific competition was strong, selection favored early broods within a broader range of pool types. The agents of selection in this study-biotic interactions and disturbance-are common to many ecological systems and frequently exhibit temporal variation; this suggests that fluctuating selection may be widespread in natural populations.

Fluctuating selection plays a key role in evolutionary models of dispersal, phenotypic plasticity, bet‐hedging, the genetic load, and maintenance of genetic variation (Lynch and Gabriel 1987; McPeek and Holt 1992; Lande and Shannon 1996; Ghalambor et al. 2007; Simons 2011; Rees and Ellner 2019). Over longer time scales, fluctuating selection may strengthen stabilizing selection or create the appearance of stabilizing selection (Travis 1989; Lande and Shannon 1996; Tiffin and Rausher 1999; Lande 2007), and this has important consequences for interpreting macroevolutionary patterns of stasis in the fossil record (Lieberman and Dudgeon 1996; Gingerich 2001; Estes and Arnold 2007; Hunt and Rabosky 2014; Voje et al. 2018). Taken together, these ideas motivate attention to the empirical issue of how prevalent fluctuating selection is in nature.

Although it is often argued that natural selection varies in space and time (Lieberman and Dudgeon 1996; Grant and Grant 2002; Bell 2010; Calsbeek et al. 2012; Voje et al. 2018), empirical support for variable or fluctuating selection has been elusive. A database of long‐term selection studies compiled by Siepielski et al. (2009) contained little evidence for temporal fluctuations in directional selection after accounting for variation caused by sampling error in the original estimates (Morrissey and Hadfield 2012). Empirical evidence for spatial variation in the strength of directional selection is also surprisingly weak (Siepielski et al. 2013). However, these results contradict many arguments and observations about natural selection and variability in the environment (Thompson 2005; Estes and Arnold 2007; Bell 2010; Futuyma 2010; Calsbeek et al. 2012; Hunt and Rabosky 2014). This has led some to conclude that we simply do not yet have enough data of sufficient quality to detect the variability in natural selection that must be there (Kingsolver and Diamond 2011; Morrissey and Hadfield 2012; Siepielski et al. 2013; de Villemereuil et al. 2020).

General insight into the importance of fluctuating selection could emerge from a better understanding of ecological mechanisms that cause variability in selection (Calsbeek et al. 2012; Morrissey and Hadfield 2012). For example, several long‐term studies connect annual variation in the availability or timing of food resources with changes in individual performance and natural selection (Grant and Grant 2002; Charmantier et al. 2008; McAdam et al. 2019). Although these studies did not experimentally manipulate food resources, they are nevertheless convincing because they demonstrate the causes of fluctuating selection in addition to establishing the broader ecological context of the evolutionary process. They also provide insight into the general conditions under which we should expect to observe fluctuating selection. In this report, we take a similar approach to explore how annual variation in ecological conditions drives fluctuating selection on habitat choice in an amphibian.

This article reports data from a 16‐year field study of a frog population at Isle Royale, MI, USA. In each year, we recorded the habitats chosen by ovipositing adults and evaluated the fitness consequences of those choices by measuring natural selection on their tadpoles. There was strong annual variation in natural selection on habitat choice, affecting directional, quadratic, and correlational aspects of selection. Moreover, we discovered that annual variation in weather conditions and the importance of species interactions were the main ecological causes of fluctuating selection. Our data support an unusually detailed portrait of the temporal dynamics of natural selection and how it arises from an ecological context that is variable and unpredictable, but nevertheless has a clear causal interpretation.

Methods

THE HABITAT OCCUPIED BY TADPOLES

Our study area was located on North Government Island (2.9 ha), in Lake Superior, just off the northeastern coast of Isle Royale (maps and photographs in Fig. 1 and Fig. S1). Adult chorus frogs (Hylidae: Pseudacris maculata) overwinter in coniferous forest and lay eggs between May and July in small unvegetated pools along the exposed bedrock shores of small islands (Figs. 1C, 2A; Smith 1983, 1990). The exposed shore is 34 m from lake to forest, and it slopes gently upward to 8 m higher than the lake at the forest edge. By searching the rock pools every two days, we were able to detect broods of tadpoles soon after they appeared, recapture and measure them repeatedly while they developed, and record the incidence of ecologically relevant factors such as predation, competition, and disturbance.

Figure 1

Maps and photographs of our study area on North Government Island at the northeastern end of Isle Royale in northwestern Lake Superior. Gray shading represents freshwater lakes and the purple shading in (B) indicates open rock shores. The photograph in (C) shows a view toward the southwest from our study area: a 75‐m section of the shoreline, about 35 m broad and containing 121 pools, toward the northeastern end of North Government Island. Amplexing pairs of Pseudacris maculata frogs oviposit and tadpoles develop in the small pools that form in depressions in the rock. Photo credit: J. Van Buskirk.

Ovipositing frogs at Isle Royale can choose among pools arrayed along three gradients of habitat variation. The first is the spatial location of the pool along an axis from the edge of the lake to the edge of the forest. This axis represents a gradient in disturbance, because lower pools are more frequently swept clean by storm waves from Lake Superior (Smith 1983). Wave disturbance in turn creates a gradient in the risk of predation and habitat permanence for tadpoles. Larval dragonfly predators (Aeshna juncea) are sensitive to disturbance, and therefore reach high density only in large pools that are more than about 60% of the distance between the lake shore and the forest (Fig. 3A). Mortality due to wave wash is important only in the lower half of the shore (Fig. 3B). Tadpoles that are washed into Lake Superior are effectively dead, because the lake contains fish and is too cold for tadpoles to complete metamorphosis.

Figure 3

Distribution of dragonfly larvae and the probabilities of drying and wave wash. (A) Average density of Aeshna dragonflies was highest in pools that are large and close to the forest. (B) The proportion of Pseudacris maculata broods that were washed into the lake by storm waves was highest for pools that are close to the lake. (C) The proportion of broods that dried before reaching metamorphosis was highest in small pools. Location on the shore is the fraction of distance between the edge of Lake Superior and the forest. Small black symbols represent pools within which no Aeshna occurred (A) or for which the source of mortality was never observed (B, C). Panel (A) depicts all pools; panels (B) and (C) include only pools in which at least three broods of P. maculata were recorded.

The second habitat axis is the size of the pool, which ranges from 0.03 m2 to 17 m2 (median 0.38 m2). Predation risk and habitat permanence also correlate with pool size, because dragonfly larvae are confined to large pools and drying primarily affects small pools (Fig. 3A and C). A tadpole will perish if the pool dries before it reaches metamorphosis.

The timing of breeding represents a third habitat axis (Fig. 2B). Broods that appear in the pools early in the summer face a different set of conditions than those that arrive late. Water temperatures rise steadily during the season, from an average of 13 °C in late May to about 20 °C in mid‐July (Fig. S2A). This provides a more favorable environment for growth later in the season. But the advantages of late arrival may be offset by increasing levels of competition for food, because the density of older tadpoles in the pools climbs from under 1·m−2 during late May to about 40·m−2 in late June, until density declines somewhat after the onset of metamorphosis in July (Fig. S2B). Moreover, earlier work indicates that late oviposition is associated with late metamorphosis in this population, which is strongly correlated with reduced survival to reproduction at age two years (Smith 1987).

Figure 2

Spatial and seasonal distribution of Pseudacris maculata broods on the study area. (A) Distribution of broods with respect to pool surface area and location. Contour intervals represent the estimated number of broods per day that appeared in the pool. Red symbols are pools that were used at least once by ovipositing frogs; white symbols are pools that were never used. Location on the shore is the fraction of the distance between the edge of the lake and the forest. (B) Frequency distribution of hatching dates for 864 broods over 16 years from 1983 to 1998. (C) The number of broods detected in each year is listed along the left side, shaded regions are periods during which we were present on the study area, and vertical ticks indicate estimated hatching dates.

All three of these axes reflect features of the habitat over which adult frogs can presumably exert some choice. An adult female can choose the date at which she arrives on the shore to oviposit, and once there she can visit pools of any size and at any distance from the forest.

MONITORING ADULT HABITAT CHOICE

We recorded the habitats chosen by ovipositing frogs by noting the appearance of eggs and hatchlings during searches of all pools at two‐day intervals, and counting tadpoles within a week after they hatched (stage 25; Gosner 1960). The appearance of eggs or a brood of tadpoles in a pool was registered as a choice by a pair of frogs of that particular habitat, defined by pool location, surface area, and hatching date. This method of scoring habitat choice is subject to error if ovipositing pairs visit more than one pool, or if more than one pair visits a pool on the same night. Both of these events do occur. Nocturnal observations of adult frogs suggest that up to 15% of the records may involve more than one pair ovipositing in the same pool on the same night, and this is especially common in large pools (Smith 1990). Therefore, we may have underestimated the number of times large pools were chosen for oviposition. An error in the opposite direction will occur if pairs visit more than one pool while in amplexus, because we will overestimate the number of distinct broods appearing on that date. An upper limit on the frequency of this error comes from examining the spatial distribution of broods that appeared simultaneously. Averaging over years, 10.5% of broods were recorded within 24 hours of another brood in a different pool within 3 m. These broods could have been produced by the same pair moving between nearby pools.

For the 16‐year period from 1983 to 1998, we were always present on the study area during June, but our coverage of earlier and later periods varied among years (Fig. 1C). On average, we initiated fieldwork on 26 May (extremes 12 May–8 June) and left the study area on 16 July (extremes 24 June–1 August).

ESTIMATING TADPOLE FITNESS

The fitness consequences of habitat choice by adult frogs were assessed by tracing the fate of the tadpole broods that resulted from their choices. The practice of estimating selection on parental characters by assigning components of offspring fitness to their parents can produce biased estimates (Clutton‐Brock 1988; Grafen 1988; Wolf and Wade 2001). But the approach works well when offspring performance is determined strongly by the phenotype of the parents and when offspring traits that influence their own performance are genetically uncorrelated with the parental characters under analysis (Grafen 1988; Wolf and Wade 2001). These conditions are likely to be upheld in our system. Habitat choice by adults can place tadpoles at risk of drying, wave wash, competition, or predation; these circumstances are not under the control of the offspring. Moreover, properties of a tadpole that influence its performance—related to characters such as feeding behavior and tail morphology (Smith and Van Buskirk 1995)—are probably genetically uncorrelated with habitat selection by adult females. Thus, in this case, the fitness consequences of parental habitat choice may be safely estimated by the performance of their offspring.

We estimated the distribution of offspring fitness using capture and recapture of broods separable by body size and time. Each brood was captured, counted, and measured on at least two occasions separated by an average of 16 days. On each capture, we removed all tadpoles during repeated searches until no new individuals were encountered. The body size distribution of the brood was estimated from a haphazard sample of 15 individuals, for which we measured the length of the body excluding the tail to the nearest 0.5 mm. These captures were used to estimate survival at the level of the tadpole even though individuals were not uniquely identifiable. For example, if there were 10 tadpoles in a brood on the first capture and five on the second capture, then five were scored as surviving and five were scored as dying. Daily proportional growth rate in body length was assigned to the survivors assuming that survival was independent of body size and rank‐order size did not change. We did not make a second capture during the 1984 field season, so there are no data on tadpole performance in that year.

We calculated measures of growth and survivorship from the brood‐level data to obtain the expected fitness for individual tadpoles. The date of first capture (DT1) and average proportional growth rate (GR) for each brood were used to estimate the date on which its average body length was 3.5 mm (DH), which we used as a measure of hatching date. The estimated date of hatching was: DH = DT1 – log(SZ1/3.5)·GR−1, where SZ1 is the mean body size of the brood on DT1. For 146 broods, we did not have a measure of early growth rate, and therefore we estimated their growth using data from 452 broods that were first measured at <4 mm. These data showed that GR depends on initial tadpole density (DN; nr.·m−2) and average water temperature (T;°C) as follows: GR = 0.0359 – 1.467e‐05·DN + 8.530e‐04·T. Metamorphic condition is determined by growth and development during the larval stage. Size at metamorphosis (SM; mm) increases with daily proportional growth rate as follows: SM = 10.34 + 28.84·GR, and date at metamorphosis (DM; days) is related to date of hatching (DH) and growth rate: DM = 64.52 + 0.75·DH – 186.5·GR.

Tadpole growth and development provide a reasonable measure of individual quality because field studies of amphibians have repeatedly demonstrated that adult survival or breeding condition is related to metamorphic condition (Smith 1987; Semlitsch et al. 1988; Berven 1990, 2009; Scott 1994; Altwegg and Reyer 2003). We used data from the Isle Royale P. maculata population to project the expected survival to age two years based on estimated time and size at metamorphosis: logit[survival to age 2] = −1.969 + 0.601·SM – 0.198·DM, where SM and DM are both expressed in SD units (Smith 1987). In other words, this means that froglets emerging at small size or late in the season are less likely to return to breed. This equation gives the expected survival to age 2, not the actual measured survival, and therefore does not include variance that would be added moving forward during the juvenile phase.

We used the relationships above to connect data on hatching date and growth during the interval between the first and second recaptures with the probability of survival to adulthood. This probability was calculated for individual tadpoles that survived to the second capture; individuals that died were assigned a value of zero. Survival to age 2 was then transformed to a quantity that we call expected fitness by relativizing it (dividing by the mean), multiplying by 10, and then rounding to the nearest integer. The rounding procedure enabled us to model variation in expected fitness assuming a zero‐inflated distribution, which was necessary because 44.2% of individuals had a fitness value of zero (mean expected fitness = 9.996, SD = 10.704; Fig. S3). Error introduced by rounding was small: the correlation between expected fitness and survival to age 2 predicted by regression was r = 0.9998.

In all, 1051 separable broods containing 55,537 tadpoles appeared in the pools, of which 62 broods were disqualified from analysis because they included tadpoles washed in from nearby pools by rain storms. Of the remaining 989 broods, 264 were not measured a second time due to time constraints, 46 were detected so late that we could not estimate their date of hatching, and 35 vanished for no known reason. This left 644 broods containing 42,172 tadpoles for which all data were available and their eventual fate was known.

STATISTICAL ANALYSES

We estimated linear, quadratic, and correlational selection gradients acting on pool choice of adult frogs on the three habitat axes using multiple regression of expected fitness, w, for each tadpole against habitat values for the pool that it occupied. The full model was (Lande and Arnold 1983): where L is the pool location on the shore, S is pool size (log of m2), D is the date of hatching, α is the intercept, the βs are linear selection gradients, and γs are the quadratic and correlational selection gradients. The three habitat axes were standardized before analysis (mean = 0, SD = 1), and the SD of the quadratic terms was scaled to 1. The random part of the model included intercepts for pool and year, and random slopes for heterogeneity among years against all fixed effects to estimate year‐to‐year variation in the linear and nonlinear selection gradients. Directional selection gradients and fluctuating selection on pool location, pool size, and hatching date came from a reduced version of equation (1) without higher order effects, as recommended by Lande and Arnold (1983, p. 1218). We used the full regression model to estimate quadratic and correlational selection and their heterogeneity among years.

The analysis employed a zero‐inflated Poisson distribution to accommodate the observed distribution of w (Fig. S3). The model therefore assumed that expected fitness was determined by two potentially independent, sequential events: a logistic process reflecting survival of larvae in the pools prior to metamorphosis (called “Tadpole Survival”) and a Poisson process reflecting the influence of larval growth and development time on postmetamorphic survival to age two years (called “Tadpole Quality”). Equation (1) was fitted using maximum likelihood via the Template Model Builder implemented in the R package glmmTMB (Brooks et al. 2017). The importance of fixed effects was judged by inspecting the 95% Wald confidence intervals of the parameters, and the importance of random effects was evaluated using likelihood ratio tests.

Equation (1) reveals whether selection gradients fluctuate significantly, but it does not indicate whether the direction of selection differs among years. We addressed this question by estimating the position of the phenotypic optimum relative to expressed trait values in each year, using a method proposed by Chevin et al. (2015). The model expresses fitness (W) of individuals in year t as a Gaussian function of individual trait values (z): where W max,t is the maximum fitness in year t, θ t is the trait value conferring maximum fitness in year t, and ω describes the width of the fitness function representing stabilizing selection. When trait values are scaled such that SD = 1, the directional selection coefficient (β) equals (θ – z)/(ω2 – 1) (Lande 1976). Hence, this model clarifies that fluctuating selection can arise from variation in the optimum, the mean trait value, or both. We estimated θ t, separately for each of the three traits (habitat axes), using general linear models derived by log‐transforming equation (2) (following Chevin et al. 2015). The random part of the model included random intercepts for pool and year, and random slopes for heterogeneity among years in linear and quadratic selection. The three models were fitted within a Bayesian framework using Stan, with the brms interface in R version 4.0.2 (Burkner 2017; Stan Development Team 2019). We used default uninformative priors, fit four chains each with 1000 warmup iterations and 1000 sampling iterations, and found that convergence was excellent (Ȓ was ≤ 1.01 for all parameters). Fluctuating selection was visualized by comparing the posterior distribution of estimated θ with the observed trait values in each year.

We investigated four ecological risks that may cause annual variation in linear and nonlinear selection gradients: pool drying, wave wash, predator abundance, and tadpole density. These risks were suggested by earlier findings on the importance of disturbance and biotic interactions for amphibian larvae (Smith 1983, 1990; Woodward 1983; Smith and Van Buskirk 1995; Werner et al. 2009). For each potentially causative agent, we fit a model similar to equation (1) except that the fixed effects included the ecological risk and interactions between it and all linear, quadratic, and correlational terms. These interaction terms tested for variation in linear and nonlinear selection gradients among years differing in the causative factor. Pool was included as a random effect. The four ecological risks were represented by the annual values of (1) the proportion of pools that dried between 26 May and 15 July; (2) the average wave height, as a proportion of the distance between the lake and the forest, measured every two days between 26 May and 15 July; (3) the total number of Aeshna dragonfly larvae on the study area, estimated by nocturnal mark‐recapture surveys (Van Buskirk 1993); and (4) the total number of tadpoles detected on the study area. The number of dragonflies was not available for the first three years (1983–1985). One of the six pairwise correlations among risks was significant: pool drying was less common in years with high wave wash (r = −0.55, N = 16, P = 0.027). Year‐to‐year variation in risks was largely unpredictable (Fig. S4): temporal autocorrelation was detected only for Aeshna numbers (r = −0.545 at a lag of two years; P < 0.05). These models were fit in package glmmTMB assuming a zero‐inflated Poisson distribution.

Results

The distribution of broods between 1983 and 1998 indicated that ovipositing pairs used nearly all pools at least occasionally, except for those at the edge of the lake (Fig. 2A). Broods were most frequent in large pools of intermediate height on the shore. Hatching dates ranged from May 24 to July 25, with 81% of broods hatching during the month of June (Fig. 2B).

SELECTION ON HABITAT CHOICE

There was weak natural selection on habitat choice, pooled over the entire study, partly because the strength and form of selection varied greatly among years (Table 1; Fig. 4). Stabilizing selection affected pool location—tadpoles survived best in pools at intermediate location on the shore (Fig. 4A). There was no significant natural selection on pool size (Table 1). Hatching date experienced directional selection favoring frogs that oviposited early in the season, caused mostly by tadpole survival. The form of selection on date was partly stabilizing because tadpole quality declined somewhat in the earliest and latest broods (Fig. 4C).

Table 1

Results of hierarchical linear models testing for directional, nonlinear, and correlational selection acting on the three axes of habitat choice: pool location, pool size, and date of breeding

	Fixed Effects	Random Effects
Source	Estimate	Source	Var. Comp.	LR‐stat
(A) Tadpole survival (logistic part) †
Intercept	−0.840 (−2.308, 0.629)	Pool	4.7285	3964.9
Location	−0.354 (−1.117, 0.409)	Year	5.0019	8458.9
Location2	0.700 (0.158, 1.243)	Location (year)	1.6325	802.1
Pool size	–0.498 (−1.043, 0.048)	Location2 (year)	0.2653	378.3
Pool size2	0.206 (−0.216, 0.627)	Pool size (year)	0.5979	856.2
Date	1.125 (0.654, 1.596)	Pool size2 (year)	0.1338	122.1
Date2	0.330 (−0.335, 0.995)	Date (year)	0.8209	1321.8
Location × pool size	–0.293 (−0.747, 0.162)	Date2 (year)	1.4378	275.2
Location × date	–0.133 (−0.354, 0.089)	Location × pool size (year)	0.2823	297.8
Pool size × date	0.225 (−0.163, 0.613)	Location × date (year)	0.1327	116.8
Pool size × date (year)	0.4783	192.4
(B) Tadpole quality (Poisson part)
Intercept	2.949 (2.744, 3.153)	Pool	0.0727	6946.2
Location	–0.065 (−0.165, 0.035)	Year	0.0583	18780.5
Location2	0.076 (−0.081, 0.233)	Location (year)	0.0279	4160.0
Pool size	0.020 (−0.061, 0.100)	Location2 (year)	0.0730	702.2
Pool size2	–0.011 (−0.143, 0.121)	Pool size (year)	0.0154	2601.0
Date	–0.053 (−0.215, 0.109)	Pool size2 (year)	0.0514	1465.2
Date2	–0.376 (−0.776, 0.024)	Date (year)	0.1017	2787.6
Location × pool size	0.068 (−0.021, 0.156)	Date2 (year)	0.5630	972.2
Location × date	0.251 (−0.093, 0.596)	Location × pool size (year)	0.0153	677.3
Pool size × date	–0.101 (−0.205, 0.004)	Location × date (year)	0.4336	442.2
Pool size × date (year)	0.0300	554.7

In part A, the model estimates the probability that mortality arises from the zero‐inflation process (Lambert 1992); positive coefficients mean that the independent variable was associated with lower survival. The response is expected fitness as defined in the text, which showed a zero‐inflated Poisson distribution (Fig. S3). Directional selection gradients and their variation among years came from a reduced version of equation (1), as described in the text. To calculate stabilizing/disruptive selection gradients, the quadratic regression coefficients shown here must be doubled (see equation (1); Stinchcombe et al. 2008). For fixed effects, covariates were standardized and the table reports the mean of the estimate with the 95% Wald confidence interval in parentheses. For random effects, the table shows the estimated variance component and the likelihood ratio (LR) statistic from nested models with and without the random effect. All random effects were significant at P < 0.0001. Boldface highlights significant effects. Sample size is 81 pools, 15 years, and 42,172 individual tadpoles.

Figure 4

Annual variation in natural selection on three components of habitat choice in Pseudacris maculata frogs. Fitness is the expected survival of offspring to age two years, relativized across the entire sample. Curves show the relationships between fitness and aspects of habitat choice predicted from the model in equation (1). The heavy black curve is based on the population‐level fixed effects, orange curves are model predictions for the separate years, and tick marks show the distribution of observations (N = 42,172 individuals).

All aspects of selection varied among years (Fig. 4). This included linear, quadratic, and correlational selection, targeting all three habitat dimensions and involving both tadpole survival and tadpole quality (Table 1). For example, tadpoles appearing early in the season generally enjoyed higher survival, but there was considerable heterogeneity in directional and stabilizing selection on hatching date such that some years differed substantially from the long‐term average (Fig. 4C). Correlational selection involving all three pairwise relationships between the habitat choice traits also showed highly significant year‐to‐year variation (Table 1).

Some of the exceptional years visible in Figure 4 can be explained by specific events that we observed in the field. The clearest example occurred in 1994, when a powerful storm on June 14 washed away 93.4% of the 3013 tadpoles that were alive at the time, and only 163 tadpoles appeared in new broods after the storm. This event is visible in Figure 4A as a year with low relative fitness except in pools near the edge of the forest.

Much year‐to‐year variation in selection was caused by changes in the estimated phenotypic optimum (θ t in equation (2)). Observed values of pool location, pool size, and hatching date (purple in Fig. 5) were somewhat variable among years, but estimates of θ (yellow in Fig. 5) were much more variable. As a result, the mean phenotype expressed in the population was in some years too large and in others too small. For example, directional selection favoring the upper pools occurred in 1988 when θ for pool location moved toward the forest; selection for use of low pools was observed in 1997 and 1998 when θ moved closer to the lake (Fig. 5A). In both cases, adult frogs bred in pools at similar location but fluctuating selection was created by movement of the optimum. This was less true for hatching date because the optimum was nearly always earlier than the observed value, except in 1994 when the early broods were destroyed by waves (Fig. 5C).

Figure 5

Fluctuating selection on three components of habitat choice was caused primarily by annual variation in the position of the optimum. In yellow are the median and 80% credible intervals of the posterior distribution of the estimated optimum (θ t in equation (2)); in purple are the means and 80% quantiles of the observed trait values for each year.

ECOLOGICAL CAUSES OF VARIATION IN NATURAL SELECTION

Annual variation in selection on pool choice was significantly associated with four ecological risks to tadpoles: pool‐drying, wave wash, predation risk, and intraspecific competition (Table 2). In dry years, selection favored frogs that chose pools in the lower part of the shore, whereas tadpoles in upper pools were favored in wet years when few pools dried (Fig. 6A). In years with heavy wave wash, selection shifted strongly to favor tadpoles farther from the lake in the upper part of the shore (Fig. 6B). When dragonfly predators were abundant, broods performed better if they appeared in lower pools below the region where Aeshna are found (Fig. 6C). Years with high tadpole density created stabilizing selection favoring tadpoles at intermediate location (Fig. 6D). In most cases, ecological effects on fitness surfaces were caused by changes in both tadpole survival and tadpole quality (Table 2).

Table 2

Summary of four mixed‐effects linear models testing for effects of ecological risks on year‐to‐year variation in selection acting on the three axes of habitat choice

	Ecological Risk
	Proportion Pools Drying	Mean Wave Height	No. of Dragonflies	No. of Tadpoles
(A) Tadpole survival (logistic part) †
Location	0.060 (−0.411, 0.532)	0.011 (−0.475, 0.496)	−0.008 (−0.412, 0.396)	−0.036 (−0.525, 0.454)
Location2	0.326 (−0.004, 0.655)	0.344 (0.004, 0.683)	0.522 (0.235, 0.810)	0.354 (0.012, 0.696)
Pool size	0.143 (−0.354, 0.640)	0.062 (−0.449, 0.572)	−0.001 (−0.427, 0.425)	0.106 (−0.404, 0.615)
Pool size2	0.237 (−0.031, 0.505)	0.170 (−0.104, 0.444)	0.294 (0.053, 0.534)	0.164 (−0.109, 0.438)
Date	0.510 (0.480, 0.540)	0.070 (0.037, 0.102)	0.355 (0.322, 0.388)	0.462 (0.431, 0.494)
Date2	0.031 (−0.004, 0.065)	0.192 (0.160, 0.224)	0.056 (0.022, 0.091)	−0.037 (−0.066, −0.008)
Location × pool size	−0.018 (−0.268, 0.231)	−0.037 (−0.292, 0.219)	0.046 (−0.174, 0.266)	−0.003 (−0.261, 0.255)
Location × date	−0.068 (−0.098, −0.037)	−0.072 (−0.100, −0.043)	−0.154 (−0.187, −0.122)	−0.061 (−0.089, −0.033)
Pool size × date	0.130 (0.099, 0.161)	0.105 (0.077, 0.132)	0.187 (0.156, 0.218)	0.107 (0.080, 0.134)
Risk	–0.392 (−0.431, −0.353)	0.408 (0.367, 0.449)	−0.041 (−0.083, 0.000)	0.163 (0.124, 0.201)
Risk × location	0.150 (0.116, 0.184)	−0.033 (−0.063, −0.002)	−0.029 (−0.061, 0.003)	−0.114 (−0.145, −0.083)
Risk × location2	0.178 (0.145, 0.211)	−0.132 (−0.166, −0.098)	0.077 (0.043, 0.111)	0.200 (0.169, 0.231)
Risk × pool size	−0.133 (−0.166, −0.101)	−0.014 (−0.045, 0.017)	0.114 (0.081, 0.147)	0.020 (−0.011, 0.051)
Risk × pool size2	−0.006 (−0.043, 0.030)	0.009 (−0.026, 0.045)	0.169 (0.134, 0.205)	−0.129 (−0.164, −0.094)
Risk × date	0.683 (0.644, 0.722)	−0.389 (−0.424, −0.354)	0.240 (0.208, 0.273)	0.125 (0.096, 0.154)
Risk × date2	0.002 (−0.039, 0.044)	0.074 (0.037, 0.110)	−0.186 (−0.225, −0.147)	0.165 (0.132, 0.198)
Risk × location × pool size	0.053 (0.023, 0.084)	–0.015 (−0.047, 0.017)	0.049 (0.017, 0.081)	−0.102 (−0.130, −0.074)
Risk × location × date	−0.006 (−0.041, 0.030)	−0.346 (−0.377, −0.314)	0.391 (0.357, 0.424)	0.053 (0.025, 0.081)
Risk × pool size × date	0.175 (0.141, 0.209)	−0.093 (−0.124, −0.062)	−0.040 (−0.071, −0.008)	0.053 (0.024, 0.083)
(B) Tadpole quality (Poisson part)
Location	−0.053 (−0.127, 0.020)	−0.020 (−0.108, 0.069)	−0.076 (−0.154, 0.003)	−0.036 (−0.109, 0.037)
Location2	0.054 (0.002, 0.105)	0.030 (−0.031, 0.092)	0.043 (−0.012, 0.098)	0.014 (−0.037, 0.066)
Pool size	0.030 (−0.045, 0.105)	0.041 (−0.049, 0.132)	0.026 (−0.054, 0.106)	0.003 (−0.072, 0.077)
Pool size2	0.001 (−0.041, 0.044)	−0.023 (−0.074, 0.028)	0.017 (−0.030, 0.064)	−0.009 (−0.051, 0.034)
Date	−0.103 (−0.107, −0.098)	−0.054 (−0.059, −0.049)	−0.032 (−0.038, −0.027)	−0.107 (−0.112, −0.102)
Date2	−0.083 (−0.088, −0.078)	−0.092 (−0.096, −0.087)	−0.139 (−0.144, −0.133)	−0.048 (−0.052, −0.043)
Location × pool size	0.016 (−0.025, 0.057)	−0.004 (−0.052, 0.045)	0.018 (−0.026, 0.062)	−0.002 (−0.043, 0.040)
Location × date	0.039 (0.034, 0.043)	0.008 (0.004, 0.013)	0.020 (0.015, 0.025)	0.033 (0.029, 0.037)
Pool size × date	0.025 (0.020, 0.030)	−0.016 (−0.021, −0.011)	0.002 (−0.003, 0.008)	0.022 (0.017, 0.026)
Risk	0.015 (0.010, 0.021)	−0.005 (−0.012, 0.001)	0.023 (0.017, 0.029)	−0.097 (−0.102, −0.092)
Risk × location	−0.071 (−0.075, −0.066)	0.060 (0.055, 0.065)	–0.022 (−0.026, −0.017)	0.035 (0.030, 0.039)
Risk × location2	0.022 (0.017, 0.027)	−0.023 (−0.028, −0.018)	0.002 (−0.003, 0.006)	−0.001 (−0.006, 0.003)
Risk × pool size	−0.022 (−0.027, −0.018)	0.024 (0.019, 0.029)	−0.057 (−0.062, −0.052)	0.010 (0.005, 0.015)
Risk × pool size2	0.020 (0.014, 0.026)	0.016 (0.010, 0.022)	−0.007 (−0.012, −0.001)	0.020 (0.015, 0.026)
Risk × date	–0.083 (−0.088, −0.079)	0.067 (0.062, 0.072)	−0.123 (−0.127, −0.118)	0.057 (0.052, 0.062)
Risk × date2	−0.004 (−0.009, 0.001)	−0.016 (−0.021, −0.011)	0.055 (0.049, 0.060)	−0.040 (−0.045, −0.035)
Risk × location × pool size	0.025 (0.020, 0.030)	−0.005 (−0.010, 0.000)	0.040 (0.035, 0.044)	0.004 (0.000, 0.009)
Risk × location × date	0.049 (0.045, 0.054)	−0.036 (−0.040, −0.031)	−0.028 (−0.033, −0.024)	−0.052 (−0.056, −0.047)
Risk × pool size × date	0.048 (0.043, 0.053)	−0.053 (−0.058, −0.048)	−0.052 (−0.057, −0.047)	−0.040 (−0.045, −0.035)

The logistic part of the model estimates the probability that mortality arises from the zero‐inflation process; positive coefficients indicate that the independent variable was associated with lower survival.

The response was expected fitness, which showed a zero‐inflated Poisson distribution (Fig. S3). Trait values and the four ecological risks were standardized (mean = 0, SD = 1). Random intercepts for pool were included in both parts of the model. The table reports the mean (95% Wald CI in parentheses). Boldface highlights estimates for which the CI did not overlap zero. The ecological risks are the proportion of pools on the study area that dried in each year; the average wave height measured at two‐day intervals, as a proportion of the distance from the lake to the forest; the log of the number of two‐ and three‐year‐old dragonfly larvae on the study area; and the log of the number of tadpoles on the study area in each year. N = 80 pools and 42,172 individuals (for dragonflies, N = 77 pools and 38,465 individuals).

Figure 6

Ecological causes of annual variation in natural selection on habitat choice with respect to pool location by ovipositing Pseudacris maculata frogs. Surfaces depict fitted values from the models shown in Table 2. Red symbols represent the trait values of all individuals (N = 42,172; in panel C, N = 38,465 individuals). Selection favored lower pools near the lake shore during years that were relatively dry and had few waves or competitors.

The corresponding analyses for pool size demonstrated that selection was strongly altered by the four ecological risks (Table 2; Fig. S5). During dry years, survival declined in small pools but the quality of surviving tadpoles increased; the resulting pattern showed stronger stabilizing selection for pools of somewhat larger size (Fig. S5A). With increasing numbers of dragonflies, natural selection favored frogs that used smaller pools. In years with high competition, both tadpole survival and tadpole quality improved in larger pools where densities were lower. Selection on hatching date was also strongly modified by all four ecological risks, with early breeding especially favored in years with many drying pools, few waves, and many predators (Table 2; Fig. S6).

The strength and configuration of correlational selection was also modified by the four ecological risks (Table 2). For example, enhanced risk of pool drying twisted the shape of the bivariate expected fitness surface on hatching date and pool location such that earlier breeding in pools closer to the lake was favored (Fig. 7A and B). Increased wave wash caused a shift from maximum fitness early in the season low on the shore to a sloping surface favoring later breeding in upper pools (Fig. 7C and D). The influence of wave wash on correlational selection played out in opposite directions for tadpole survival and tadpole quality (Table 2). The causal effect of wave wash through survival (reflecting improved survival later in the season, when storms are less frequent, and in upper pools, which are rarely reached by storm waves) was strong enough to overwhelm the effect of wave wash through tadpole quality (reflecting reduced competition in pools that were cleared out by storm waves). Years with high dragonfly numbers caused a transition from a correlational saddle that favored later breeding in higher pools to directional selection favoring lower pools where dragonflies do not occur (Fig. 7E and F). Finally, when tadpoles were abundant, changes in both survival and tadpole quality reorganized the selection surface to create strong stabilizing selection for breeding at intermediate dates and lower on the shore (Fig. 7G and H). Correlational selection involving the other two pairs of habitat axes was also modified by the four ecological risks, sometimes in complex ways (Table 2; Figs. S7, S8).

Figure 7

Ecological causes of annual variation in correlational selection on habitat choice with respect to hatching date and pool location by ovipositing Pseudacris maculata frogs. Surfaces depict fitted values from the models shown in Table 2, with high relative fitness indicated by the warm orange color. For panels on the left side, the surface was estimated using a value of the ecological agent equal to the average of the three lowest years; for panels on the right side, the value of the agent was the average of the three highest years. These average values are given in parentheses above each panel. Red symbols are the values of broods occurring in years during which the ecological agent was low (left side) or high (right).

Discussion

This study documents pervasive fluctuating selection on choice of breeding habitat in a frog. There was annual variation in directional and quadratic selection for all three habitat traits, and in correlational selection involving all combinations of traits. Fluctuating selection was associated with annual changes in the phenotypic optimum. These results confirm the view that fluctuating selection probably occurs frequently (Stanley and Yang 1987; Lieberman and Dudgeon 1996; Estes and Arnold 2007; Bell 2010; Hansen 2012; Hunt and Rabosky 2014). However, empirical evidence for fluctuating selection is surprisingly thin (Morrissey and Hadfield 2012). In the best‐studied cases, naturally occurring variation in a known agent of selection can be tied to variation in selection coefficients (e.g., Grant and Grant 2002; Charmantier et al. 2008; McAdam et al. 2019; de Villemereuil et al. 2020). In our study, annual changes in the fitness surface arose from variation in competition, predation, and disturbance— agents that modify which pools are favorable for tadpole survival and growth. To be sure, causal interpretations would be more convincing if the agents of selection were manipulated (Wade and Kalisz 1990), but we argue that correlations over many years establish a good case for causation and support general conclusions about the ecological context of fluctuating selection.

Spatial and temporal habitat variation relevant for ovipositing frogs is clearly apparent in the Isle Royale rock pools. A striking property of the three main habitat axes—pool location, size, and time of the season—is that they are associated with ecological or physical factors that tend to exert opposite and compensatory influences on tadpole performance (Fig. 3, Fig. S2). Pool location on the shore represents a tradeoff between predation risk up high and probability of wave wash down low, so that frogs avoiding one risk factor suffer increased exposure to the other (Smith 1983). Pool size represents a tradeoff between predation risk and probability of pool drying. Time of breeding defines a tradeoff between cold temperature early in the season and increasing competition for food as the summer progresses (Fig. S2).

The same factors that affected performance along the three habitat axes also created fluctuating selection on habitat choice. A good example was predation risk, represented by the number of Aeshna dragonfly larvae on the study area. In years with many predators, selection favored frogs that chose smaller pools low on the shore—this is the habitat where Aeshna rarely occurs. Selection for small lower pools is interpretable in terms of the known effects of dragonflies on larval anurans, including direct mortality, suppressed foraging activity, and reduced growth and development rate (Smith and Van Buskirk 1995; Benard 2004). But in years with few dragonflies, selection shifted and frogs did better if they bred somewhat later in the season and in large, upper pools (Fig. 6, Figs. S5, S6). The rock surface in upper pools supports more algae, so enhanced tadpole growth in the upper habitat probably explains the change in natural selection when predators were scarce.

Variation in the other three ecological agents also created fluctuating selection in interpretable ways. In years with many drying pools, selection favored frogs that bred early in the season, well before temperature increased and rainfall declined in July. Annual variation in wave wash was strongly associated with selection on pool location and hatching date, because heavy waves tend to occur early in the season and primarily destroy broods that are close to the lake. In years with high tadpole numbers, which indicate the strength of competition for food, survival, and growth rate were reduced and selection shifted to favor frogs that bred earlier and used larger pools of intermediate height on the shore. This was because resource competition is reduced in large pools (which tend to be less crowded) and early in the season (before other tadpoles have arrived). These mechanistic connections between annual changes in natural selection and variation in identifiable causal agents help explain the shifts in the phenotypic optimum that we observed (Fig. 5).

Our mechanistic interpretation of causal agents also suggests that fluctuating selection is probably widespread. The specific biotic interactions and physical processes that create tradeoffs in our study may be particular to the Isle Royale rock pools, but compensatory effects of the physical and biotic environment occur generally on ecological gradients (Grime 1973; Wilbur 1980; Wellborn et al. 1996; Hopcraft et al. 2010). Species distributions along such gradients are frequently associated with functional traits that influence performance with respect to physical factors such as disturbance and biotic interactions such as resource competition or predation (Woodward 1983; Smith and Van Buskirk 1995; Werner et al. 2007; Weiher et al. 2011; Wong et al. 2019). Environmental variation shifts the balance among controlling factors in many systems (Chamberlain et al. 2014; Piovia‐Scott et al. 2017; Carvajal‐Endara et al. 2020), and this is likely to cause fluctuating selection in general, just as it does in the Isle Royale rock pools.

Opposite and compensatory risks along habitat axes should average out to create stabilizing selection (Travis 1989; Tiffin and Rausher 1999). Indeed, the overall pattern of selection on all three habitat axes was weakly stabilizing. Fitness surfaces contained ridges or oblong peaks under certain conditions, but the general picture was that pool location, pool size, and linear combinations of the two were frequently under stabilizing selection (Figs. 4, 7, Fig. S7). Quadratic and correlational selection gradients were of roughly the same order of magnitude as linear gradients (Tables 1, 2). Previous studies comparing estimates of natural selection reveal that quadratic and correlational selection are typically weaker than directional selection (Kingsolver et al. 2001; Kingsolver and Diamond 2011). However, methodological issues may cause underestimation of nonlinear selection (Blows and Brooks 2003; Hereford et al. 2004; Kingsolver and Diamond 2011; Haller and Hendry 2014). Analyses specifically tailored to estimate the joint magnitudes of directional, quadratic, and correlational selection indicate that nonlinear selection can be widespread and strong (Simms 1990; Schluter and Nychka 1994; Blows and Brooks 2003; McGuigan et al. 2011).

One motivation for estimating the strength of nonlinear selection is that stabilizing and correlational selection influence the magnitude of the evolutionary load. A population experiences load whenever individuals express phenotypes that are not optimal. The total evolutionary load is the sum of L G, the variance load (caused by the deviation of individual phenotypes from the optimum due to genetic and environmental variance), and L θ, the lag load (due to deviation of the population mean phenotype from the optimum, θ) (Maynard Smith 1976; Lande and Shannon 1996; Chevin 2013). Load is highest when stabilizing selection is very strong (so that small deviations from θ reduce fitness substantially) or when the environment changes rapidly. For example, an analysis of temporal variation in selection on breeding date in 12 bird species and seven mammals estimated that lag load caused by annual variation in the optimum date was about 0.20 for birds and 0.40 for mammals (de Villemereuil et al. 2020). These values are high because the estimated width of the stabilizing selection curve was narrow. McAdam et al. (2019) also reported high L θ in a study of litter size in squirrels (≈0.25), but in this case load was caused by annual changes in the environment associated with a masting food source. In our study, optimal trait values fluctuated strongly among years but the stabilizing fitness function was broad; hence, estimates of load were low (Fig. S9). Annual values of the total load associated with pool location and pool size were between 0.013 and 0.067, with L θ accounting for about 60–90% of the total. Stabilizing selection on hatching date was even weaker still, and therefore load caused by that trait was negligible. Taken together, these studies illustrate how data on the strength of nonlinear selection, combined with information on θ and the distribution of phenotypes, can reveal differences among populations in the extent of maladaptation. Comparisons over time can then reveal the causes. So far, the conclusion is that maladaptation can be pervasive in populations that undergo fluctuating selection, and this is often related to changes in food resources (Grant and Grant 2002; Charmantier et al. 2008; McAdam et al. 2019; de Villemereuil et al. 2020).

Phenotypic selection for early reproduction, similar to what we observed, has been noted previously in birds and plants (Rowe et al. 1994; Munguia‐Rosas et al. 2011; de Villemereuil et al. 2020). This pattern is commonly attributed to a negative correlation between breeding time and body condition or nutrition, so that high‐quality individuals breed early and produce high‐quality offspring (Price et al. 1988; Austen et al. 2017). If body condition is primarily determined by the environment, no genetic response will occur in spite of persistent phenotypic selection (Price et al. 1988; Rausher 1992). An alternative explanation is suggested by the observation that tadpole density, and hence resource competition, increases during the season at Isle Royale (Fig. S2B). In years when crowding is severe, selection may disfavor frogs that oviposit after the others. Priority effects such as this—which are well known between species of larval amphibians (Alford and Wilbur 1985; Murillo‐Rincon et al. 2017)—could create frequency‐dependent selection favoring reproduction at an earlier time than competitors, regardless of the actual date (Haller and Hendry 2014; de Villemereuil et al. 2020). These explanations cannot be differentiated with the data at hand, but could be tested using information on the body condition of ovipositing frogs.

If strong directional and stabilizing selection is common in nature, then it is germane to ask how phenotypic and genetic variation is preserved. Explanations for the maintenance of habitat breadth based on fitness trade‐offs are poorly supported empirically (Futuyma and Moreno 1988; Sexton et al. 2017). But theory suggests that fluctuating selection could be important by increasing phenotypic variation (Slatkin and Lande 1976; Bull 1987) and genetic variance (Lynch and Gabriel 1987; Ellner and Hairston 1994; Svardal et al. 2015; Rees and Ellner 2019). Thus, even in the presence of generally stabilizing selection on habitat choice, the never‐ending movement of the optimum may prevent specialization on the habitat favored by long‐term average selection.

AUTHOR CONTRIBUTIONS

JVB designed the study, obtained funding, helped collect data, analyzed data, and wrote the manuscript. DCS designed the study, obtained funding and permits, collected the data, and contributed to manuscript preparation.

DATA ARCHIVING

The data and analysis scripts are available on dryad: https://doi.org/10.5061/dryad.3ffbg79hx.

Associate Editor: A. Siepielski

Handling Editor: A. McAdam

Supplementary figures

Click here for additional data file.