Negative frequency‐dependent selection maintains shell banding polymorphisms in two marine snails (Littorina fabalis and Littorina saxatilis)

Abstract The presence of shell bands is common in gastropods. Both the marine snails Littorina fabalis and Lttorina saxatilis are polymorphic for this trait. Such polymorphism would be expected to be lost by the action of genetic drift or directional selection, but it appears to be widespread at relatively constant frequencies. This suggests it is maintained by balancing selection on the trait or on a genetically linked trait. Using long time series of empirical data, we compared potential effects of genetic drift and negative frequency‐dependent selection (NFDS) in the two species. The contribution of genetic drift to changes in the frequency of bands in L. fabalis was estimated using the effective population size estimated from microsatellite data, while the effect of genetic drift in L. saxatilis was derived from previously published study. Frequency‐dependent selection was assessed by comparing the cross‐product estimator of fitness with the frequency of the polymorphism across years using a regression analysis. Both studied species showed patterns of NFDS. In addition, in L. fabalis, contributions from genetic drift could explain some of the changes in banding frequency. Overdominance and heterogeneous selection did not fit well to our data. The possible biological explanations resulting in the maintenance of the banding polymorphism are discussed.

for generations (Hartl & Clark, 2014). For example, genetic diversity can be maintained by a balance between migration and mutation and genetic drift, in particular, when the effective population size is sufficiently high and the trait's effect on fitness is either weak or absent (e.g., Cook, 1992;Hoffman et al., 2006;Oxford, 2005). Therefore, for years, there was a debate on whether polymorphic traits were the result of stochastic or deterministic processes (Hartl & Clark, 2014;Kimura, 1991;Wright, 1977).
One of the most powerful and intuitive deterministic mechanisms by which polymorphisms can be maintained in natural populations is negative frequency-dependent selection (NFDS; Ayala & Campbell, 1974;Fitzpatrick et al., 2007;Svensson & Connallon, 2019;Trotter & Spencer, 2008). Additional mechanisms are spatially or temporally heterogeneous selection (or multiple niche selection model; Levene, 1953) and overdominance. Importantly, the different mechanisms explaining polymorphisms are not mutually exclusive but can occur simultaneously.
The banding polymorphism of the marine flat periwinkle Littorina fabalis (Turton, 1825) and the rough periwinkle Littorina saxatilis is ideal to study these matters as it tends to comply with most of the characteristics cited in the first paragraph: It has been shown to have a simple genetic basis (Ekendahl & Johannesson, 1997;Kozminsky, 2011); it is pervasive (Figure 1; Johannesson & Butlin, 2017); and it is easily scored (Figure 1). Furthermore, these F I G U R E 1 Banded individuals of Littorina fabalis and Littorina saxatilis and their frequency across years. (a) Examples of L. fabalis individuals and sample frequency over a period of 7 years in the Spanish population of Abelleira; photographs by Dr. Daniel Estévez. (b) Examples of Littorina saxatilis individuals and average frequency over a period of 21 years for Swedish populations recolonized with banded (blue) or white and/or red (red) individuals; photographs by Fredrik Pleijel two species are good models for evolutionary research (reviewed in Rolán-Alvarez et al., 2015). Both L. fabalis and L. saxatilis are intertidal rocky shore snails present on North Atlantic shores. They are gonochoric and polygamous and show internal fertilization. Littorina fabalis is an annual species; that is, overlap between generations is minor (Reid, 1996;Rolán-Alvarez et al., 2015;Williams, 1992), while for L. saxatilis, we follow the assumption by Johannesson and Butlin (2017) that the Swedish wave ecotype of this species (from which the data of this study are derived) is reasonably annual as well. Accordingly, the effects of the different evolutionary forces can be estimated year by year. L. fabalis is oviparous (laying egg masses from which juvenile snails hatch) and associated with brown algae, for example, Fucus spp., while L. saxatilis is ovoviviparous (females carry the offspring during the larval development and give birth to miniature snails) and can be mostly found on rock surfaces. The two species have low dispersal ability, and migrations among localities separated by a few km are occasional, presumably by rafting (Johannesson & Butlin, 2017;Johannesson & Johannesson, 1995;Reid, 1996).
The data for this study were obtained from two independent research efforts. For L. fabalis, phenotypic data on banding frequency were obtained from a study on shell color and mate choice (Estévez et al., 2020), although the data on banding are presented here for the first time. For L. saxatilis, we used published data from a study by Johannesson and Butlin (2017;https://doi.org/10.5061/ dryad.427p0). For estimation of genetic drift in L. fabalis, we used an estimation of effective population size performed from microsatellite data, whereas for L. saxatilis, we followed the same interpretation as previous authors. As frequency was captured along a time series for both datasets, we could estimate candidate evolutionary forces affecting the banding pattern polymorphism. In particular, we aimed at following a similar methodology to that shown by Martínez-Rodríguez et al. (2019). Thus, the objectives of this study were as follows: (a) to determine the relative contribution of genetic drift in the maintenance of the banding polymorphism in L. fabalis and (b) to estimate whether a mechanism of NFDS contributes to the polymorphism in both species.

| Original sampling
Mating pairs and unmated individuals of L. fabalis were collected from the Fucus spp. canopy in a rocky shore located in Abelleira, A Coruña, NW Spain (42°48′0.30″N-9°1′14.87″W), for a period of 7 years from 2011 to 2017 (see further details in Estévez et al., 2020). The density of the population is relatively large compared with other populations, with tens of individuals per algae. All individuals were collected at the lower part of the intertidal to avoid collecting individuals of the sibling species L. obtusata, which inhabits the upper part (Williams, 1990). Additionally, the inspection of the sexual anatomy (see Reid, 1990Reid, , 1996 and a molecular analysis using species-specific microsatellite loci (see Carvalho et al., 2016), performed with the 2014 sample, indicated that L. obtusata individuals accounted only for 1.4% of juveniles and 1.0% of females, both difficult to classify. This indicates that the amount of error incorporated due to misidentified individuals was negligible. Identified individuals of L. obtusata were discarded before the analyses. The presence/absence of bands was recorded for a total of 4,717 specimens across the 7 years. Shell band frequencies in this population were never artificially manipulated.
Frequencies of banded L. saxatilis individuals from eight Swedish skerries were obtained by Johannesson and Butlin (2017) from a sample of 11 years out of a period of 21 years to evaluate the contribution of different evolutionary forces to the maintenance of rare color morphs. The skerries were recolonized with individuals of a single color morph coming all from a larger island (58°49′18.6″N; 11°02′19.2″E), after a toxic algal bloom wiped out the original populations in 1992 (Johannesson & Johannesson, 1995). A group of skerries (W1, W2, W3, and R) were recolonized with white and red snails, while another group of skerries (B1, B2, and B3) were recolonized with banded snails at the start of the experiment (Johannesson & Butlin, 2017). In subsequent samplings, there was no manipulation of morph frequencies. Notice that similar results in both kinds of populations would suggest that manipulation itself had no effect in our fitness estimate. One of the skerries did not have enough data to evaluate the NFDS pattern and consequently was removed from the analysis (W4). Original results from Johannesson and Butlin (2017) pointed out that either NFDS or overdominance was possible mechanisms for the persistence of low frequencies of banded snails in most populations of the area.

| Estimation of genetic drift
The effective population size (N e ) obtained by Estévez et al. (2020) was used here to test whether changes in allele frequencies for the polymorphism of bands in L. fabalis could be explained by genetic drift. The genotypes of the microsatellites analyzed in this study were used to estimate the N e using the linkage disequilibrium method as implemented in the software NeEstimator v.2.01 (Do et al., 2014). We assumed a model of one locus and two alleles with banded dominant over unbanded pattern (following Ekendahl & Johannesson, 1997;Johannesson & Butlin, 2017;Kozminsky, 2011).
Then, the observed changes in frequencies between years, assuming that each year represents a generation (see Williams, 1992), were compared with the expected values under a null model of genetic drift, where variance in the shift of allele frequency in one generation is given by (Fisher, 1922)

| Estimate of fitness
Fitness was estimated based on the frequency changes between 2 consecutive years (see, e.g., Cook et al., 1999;Martínez-Rodríguez et al., 2019). This interannual fitness estimate represents the response to selection under a quantitative genetic approach, and assuming a Mendelian genetic basis with no environmental effects, the response is equivalent to the total selection differential (Walsh & Lynch, 2018). Fitness was estimated by the cross-product estimator, which is the most suitable for qualitative traits (W; Knoppien, 1985, Rolán-Alvarez & Caballero, 2000, Martínez-Rodríguez et al., 2019 and which can capture long-term fitness better than other methods (see Brommer et al., 2004). One of its major advantages is its frequency independence, which makes it very suitable for

| Analysis of deterministic balancing mechanism
A linear regression of the percentage of banded individuals on interannual fitness [expressed as the log 10 (W)] was carried out to find patterns of frequency-dependent selection. A single linear regression ANOVA was used to assess the significance of the relationship between fitness and frequency in the L. fabalis and L. saxatilis populations across the sampling years.
In L. saxatilis, the same hypothesis, namely that fitness is a function of frequency of the allele for banding, was evaluated in each population. Complementarily, the same hypothesis was estimated using all populations by two different methods: (a) Using the average of β, the parametric value of the regression coefficient, from all populations, a t test was employed to evaluate the null hypothesis β = 0; and (b) Fisher's test was used for combining the probabilities, which takes the p value from each independent test and follows the equation: where k is the number of different tests (Sokal and Rohlf, 1995). The significance of this test is evaluated against a chi-square distribution with 2k degrees of freedom. Butlin (2017), was discarded as the explanatory force in most cases, although the contribution of genetic drift to the frequency change was not excluded in remaining models.
The linear model between frequency and fitness is represented in Figure 3. The regression line showed a negative relation between fitness and percentage of banded individuals (notice that both fitness and percentage are relative values, and hence, the same results will be observed for unbanded individuals). The relationship between frequency (using always the starting frequency, i.e., that of year + 0 in all intervals) and fitness was significant (F 1,4 = 27.31, p =0.006).

| Littorina saxatilis
Fitness estimates for L. saxatilis in the seven populations (population W4 not included) are shown in Table 2. Significant estimates were observed in both decreases and increases in frequency across year intervals. These were particularly high in cases where frequency changes were high: Fitness increased much more from low starting frequencies and vice versa.
Johannesson and Butlin (2017) studied banded frequency in eight populations. Our reanalysis presented here shows the existence of a frequency-dependent pattern, as fitness of banded specimens varies negatively with its frequency, and this occurred with the population being manipulated either for color or for banding at the start of the experiment (Table 2; Figure 3). Furthermore, notice that these patterns represent a period of 21 years in all populations, supporting the stability of the mechanism and its persistence in spite of the initial frequencies. The regression was significant only in one population although four more were marginally significant (Table 3). However, when the whole subset of populations is tested simultaneously, the result is significant: A t test on the average of β coefficients indicated an inverse relationship between fitness and frequency, being such average equal to −0.61; additionally, Fisher´s test combining the probabilities of the seven ANOVAs was strongly significant (Table 3). Following these two results, the observed pattern adjusts to what it would be expected under NFDS (Figure 3).

| D ISCUSS I ON
The maintenance of genetic polymorphisms in the absence of balancing selection (i.e., under mutation-drift-migration equilibrium) has been documented in a number of cases, with examples including color polymorphisms in the leopard frog Rana pipiens (Hoffman et al., 2006) and in the periwinkle Littoraria pallescens (Cook, 1992;Parsonage & Hughes, 2002). Even polymorphisms in Cepaea snails, where there is a long history of research on selective mechanisms, have been attributed to some extent to equilibriums between migration and drift (Bellido et al., 2002). Some of these species have large population sizes, and consequently, neutral processes could be associated with the transitory changes of these genetic polymorphisms in the absence of selective forces. Contrarily, genetic drift is expected to reduce genetic diversity in L. fabalis at Abelleira as the species presumably suffers strong bottlenecks during winter (Rolán-Alvarez et al., 1995;Williams, 1996; although other populations outside NW Spain can be more dense and stable; see Reid, 1996). The average percentage of banded individuals in the L. fabalis population, Abelleira, was 8.3% ± 1.9 SD. These observations suggest that although genetic drift could explain subtle changes in frequency for the polymorphism (Figure 2), when a change occurs over a certain threshold, a balancing mechanism protects the polymorphism from being lost. The evidence in the related species, L. saxatilis, is similar. On the one hand, it is a low-dispersal species that shows intermediate density populations along the rocky shore and is expected

TA B L E 1 Estimates of interannual (W ± SD) fitness for banded
Littorina fabalis individuals using the cross-product estimator (W) across a period of 7 years. The total number of individuals sampled (N) and the percentage of banded individuals correspond to the first year in each interval ( to have low effective population size in certain areas (Fernández et al., 2005, but see Reid, 1996). On the other hand, genetic diversity in L. saxatilis can rapidly recover after strong bottlenecks, probably because one single female is able to recolonize one population hav- B2-42.5% ± 12.3 SD; B3-51.7 ± 7.6 SD) and the ones manipulated for other colors showing the lowest (W1-3.2% ± 2.5 SD; W2-4.9% ± 2.6 SD; W3-1.9% ± 1.3 SD; R-2.8% ± 1.5 SD). The banded shell pattern in this species persisted over time, and drift could not account for all observed changes in banding frequency (Johannesson & Butlin, 2017). These trends in both species are similar to what has been proposed in candy-striped spider (Oxford, 2005). In this example, color polymorphism was found to be almost neutral within a specific range of frequencies with selection regimes being very weak. However, when drift pushed frequency of a color morph over this range, evidence was found that some form of natural selection pushed back the frequency to previous values.
Interannual fitness of banded individuals was found to decrease with frequency, thus fulfilling the necessary condition for a mechanism of NFDS (for a similar approach, see Cook et al., 1999;Martínez-Rodríguez et al., 2019). This was observed in both species ( Figure 3; Table 2) and independently of the initial or average frequency of the phenotype. It could be argued that the initial frequencies obtained for L. saxatilis could have a profound effect on the evolutionary outcome, especially whether drift had a significant contribution. However, in the period of 21 years, a balancing mechanism was observed for populations either with low or higher frequencies ( Figure 3). Johannesson and Butlin (2017) were supportive of a mechanism of NFDS maintaining rare polymorphisms in L. saxatilis but suggested that overdominance could also be a candidate mechanism. However, although overdominance has been observed to maintain polymorphisms in a number of cases (Kellenberger et al., 2019;Schou et al., 2017), it is somewhat unlikely to produce the patterns observed in Figure 3. On the contrary, one argument supporting overdominance would be that both types of experimental populations in L. saxatilis (manipulated for banding or not) would maintain the observed equilibrium frequencies irrespectively of the initial frequencies. Notice that our observations cannot definitively support or reject overdominance, since it could produce a trend similar to a frequency-dependent pattern, at least for a limited number of generations. Importantly, associative overdominance where the polymorphism is controlled by a supergene or a chromosomal inversion is gaining momentum in explaining adaptation as more of these regions are being found in various taxa (Wellenreuther & Bernatchez, 2018). In fact, several chromosomal inversions have been found in L. saxatilis (Faria et al., 2019). Deleterious mutations can accumulate in haplotypes containing inversions causing homozygotes to have lower fitness. Such haplotypes could be maintained in heterozygotes, which would be in low frequency (see Jamie & Meier, 2020, and references therein). Thus, this mechanism cannot be rejected given our current observations and deserves future consideration.
Another possible balancing mechanism is heterogeneous selection, which could theoretically maintain a trend similar to that F I G U R E 3 Observed negative relation between annual fitness (log 10 W) and frequency of banded Littorina saxatilis (B1, B2, B3, R, W1, W2, and W3) and Littorina fabalis (Abelleira) individuals in different sampled populations (notice that this frequency refers to the phenotype, not the allelic frequency). Notice that a random pattern does not necessarily produce such a trend as it may cause several changes in frequencies in the same direction (causing a direct association between fitness and frequency) TA B L E 2 Estimates of interannual (W ± SD) fitness for banded Littorina saxatilis individuals for populations manipulated for conspicuous (W1, W2, W3, and R) and banded (B1, B2, and B3) individuals. Notice that since the W estimator is a relative measurement of fitness, fitness for unbanded individuals is the inverse. NA = impossible to calculate

TA B L E 3
Linear regression reanalysis of fitness (log 10 W, dependent variable) over the percentage of banded individuals (predictor) in the 8 populations studied by Johannesson and Butlin (2017). One of the populations (3)  observed in Figure 1 for L. fabalis. However, fitness in this case would be dependent on the frequency of exploitable habitats rather than other polymorphism variants' frequencies, and thus, more erratic variation than that observed in Figure 1 would be expected. Such erratic variation would indeed be expected in L. fabalis due to the relative complexity of its habitat: The brown algae canopy presents yellow, brown, and green color variations, and its density fluctuates along the year with strong reductions during winter (Rolán-Alvarez et al., 1995). On the other hand, L. saxatilis from the studied dataset inhabited a relatively simple rocky shore that remains identical over the year. Overall, we follow the same argument exposed by Brisson (2018) that any particular rare variant would not increase its fitness unless its most suitable background is encountered frequently. Therefore, the most plausible interpretation is that there is a true frequency-dependent natural selection mechanism involved although we cannot neither prove nor reject heterogeneous selection based on our observations. Frequency-dependent, predator-driven selection (apostatic selection) has been suggested for the band polymorphism in the closely related species L. obtusata (Smith, 1976). This kind of selection has also been proposed in guppies (Olendorf et al., 2006), salamanders (Fitzpatrick et al., 2009), and land snails (Tucker, 1991), among others. Most of these studies convey in hypothesizing that apostatic selection is a result of search image employment by predators, which adjusts to the frequency of the different colors or patterns (Bond, 2007). This in turn would imply that the polymorphism evolves in a "fine-grained" environment, that is, homogeneous environments where all polymorphic variants are equally cryptic (Surmacki et al., 2013). The main predator of both L. fabalis and L. saxatilis is the shore crab Carcinus maenas (Ekendahl, 1998;Williams, 1992); however, given the poor vision of crabs, intertidal fish or birds seem to be better candidates for visual predation (see Johannesson & Ekendahl, 2002;Reimchen, 1989). A possible combination of both apostatic selection and heterogeneous selection could occur in the sense that the environmental complexity could determine the strength of apostatic selection and crypsis (Bond & Kamil, 2006 In summary, we have found that banded pattern in L. fabalis and L. saxatilis shows a NFDS pattern of change in polymorphism frequencies along years. As several other intertidal gastropods show variation in shell banding (such as Nucella lapillus and L. obtusata), we believe that similar explanations could be invoked in the evolution of this polymorphism, especially on those showing direct development, and consequently being more prone to local adaptation. For example, a common predator or group of predators could be driving local adaptation in coexisting snail species of similar sizes causing apostatic selection. Further investigations on the ultimate biological mechanism causing this pattern will be needed, as well as conclusive work that can elucidate the relative role of either apostatic or heterogeneous selection.

ACK N OWLED G EM ENTS
We thank Mary Riádigos for administrative help. We also thank Prof. Kerstin Johannesson for a critical assessment of methods and advice to improve the manuscript. This work was supported by

CO N FLI C T O F I NTE R E S T
We hereby declare that no author of the research presented here has any conflict of interest. Writing-review & editing (equal).

DATA AVA I L A B I L I T Y S TAT E M E N T
Information on frequencies of banded individuals across temporal and spatial scale, and estimates of fitness can be found in Dryad, https://doi.org/10.5061/dryad.2280g b5rb.