Mate Choice Contributes to the Maintenance of Shell Color Polymorphism in a Marine Snail via Frequency-Dependent Sexual Selection

Natural color polymorphisms are widespread across animal species and usually have a simple genetic basis. This makes them an ideal system to study the evolutionary mechanisms responsible for maintaining biodiversity. In some populations of the intertidal snail Littorina fabalis, variation in shell color has remained stable for years, but the mechanisms responsible are unknown. Previous studies suggest that this stability could be caused by frequency-dependent sexual selection, but this hypothesis has not been tested. We analyzed shell color polymorphism in mating pairs and surrounding unmated individuals in two different populations of L. fabalis to estimate sexual fitness for color, as well as assortative mating. The estimated effective population size from neutral markers allowed us to disregard genetic drift as the main source of color frequency changes across generations. Shell color frequency was significantly correlated with sexual fitness showing a pattern of negative frequency dependent selection with high disassortative mating for color. The results suggested a contribution of male mate choice to maintain the polymorphism. Finally, the implementation of a multi-model inference approach based on information theory allowed us to test for the relative contribution of mate choice and mate competition to explain the maintenance of color polymorphism in this snail species.

Natural color polymorphisms are widespread across animal species and usually have a simple genetic basis. This makes them an ideal system to study the evolutionary mechanisms responsible for maintaining biodiversity. In some populations of the intertidal snail Littorina fabalis, variation in shell color has remained stable for years, but the mechanisms responsible are unknown. Previous studies suggest that this stability could be caused by frequency-dependent sexual selection, but this hypothesis has not been tested. We analyzed shell color polymorphism in mating pairs and surrounding unmated individuals in two different populations of L. fabalis to estimate sexual fitness for color, as well as assortative mating. The estimated effective population size from neutral markers allowed us to disregard genetic drift as the main source of color frequency changes across generations. Shell color frequency was significantly correlated with sexual fitness showing a pattern of negative frequency dependent selection with high disassortative mating for color. The results suggested a contribution of male mate choice to maintain the polymorphism. Finally, the implementation of a multi-model inference approach based on information theory allowed us to test for the relative contribution of mate choice and mate competition to explain the maintenance of color polymorphism in this snail species.

INTRODUCTION
Negative frequency-dependent selection (NFDS) operates when the relative fitness of a morph within a natural population is an inverse function of its relative frequency, i.e., the more abundant a morph is, the lower its fitness (Wright, 1939). Although NFDS has been considered the most intuitive and powerful force maintaining polymorphisms (Ayala and Campbell, 1974;Fitzpatrick et al., 2007;Trotter and Spencer, 2008;Svensson and Connallon, 2019), relatively few examples exist. Thus, the biological causes for this pattern in natural populations remain largely unknown.
Patterns of NFDS can arise in several distinct, non-exclusive scenarios: driven by predation (apostatic selection), where less frequent morphs have a survival advantage (Tucker, 1991;Olendorf et al., 2006;Fitzpatrick et al., 2009;Weir, 2018); by sexual conflict, where less frequent morphs of one sex obtain fitness benefits in detriment or against fitness interests in the other sex (Svensson et al., 2005;Gosden and Svensson, 2007); by mate competition, where the less frequent morph within one sex obtains a greater number of matings than the common one when competing for them (Sinervo and Lively, 1996;Bleay et al., 2007;Rios-Cardenas et al., 2018); and by mate choice, in cases when rare or dissimilar morphs are preferred and therefore have relatively more matings (Hughes et al., 2013;Chouteau et al., 2017;Valvo et al., 2019). It becomes obvious that the main difference between these scenarios is the component of fitness that is exploited by the less frequent morph. Although the simultaneous study of all these components in one species would be a cumbersome task, both mate choice and mate competition scenarios can be easily studied together as they relate to the same component of fitness. In this study, we refer to these two scenarios as sexual fitness (following Hedrick, 1983, p.163;Bell, 2008, p.167), a component of fitness defined by the capacity of individuals to find mates.
Mate choice is the process by which certain traits in one sex lead to non-random mating (Kokko et al., 2003;Edward, 2015). This can lead to NFDS when there is a tendency by the choosy sex to mate with the less frequent or rare morphs, which consequently enjoy greater sexual fitness (Hughes et al., 2013;Valvo et al., 2019); or when there is a fixed preference to mate for dissimilar morphs (Chouteau et al., 2017). A usual proxy for inferring such a choice is disassortative mating. Whenever observed, a straightforward consequence of disassortative mating is that the polymorphism is protected from fixation or loss, independently of other findings (Takahashi and Hori, 2008;Hedrick et al., 2016Hedrick et al., , 2018Chouteau et al., 2017). Disassortative mating is, however, a rare event in nature (Jiang et al., 2013). It has therefore not been extensively studied, though a recent review suggests it may be more frequent in marine invertebrates than previously assumed (Olsen et al., 2020). A common explanation for the evolution of disassortative mating is the avoidance of inbreeding and the consequent maladaptation due to the accumulation of deleterious recessive alleles (Pusey, 1987;Tregenza and Wedell, 2000;Duthie and Reid, 2016;Olsen et al., 2020), similar to the evolution of selfincompatibility in plants (reviewed by Delph and Kelly, 2014). Another common explanation comes from studies in vertebrates, where disassortative mating has been shown to increase MHC heterozygosity and overall offspring fitness (Consuegra and Garcia de Leaniz, 2008;Huchard et al., 2013).
Mate competition, the contest of polymorphic individuals within one sex to carry out a greater number of matings (Andersson, 1994), can also lead to patterns of NFDS. Perhaps the most well-known case is that of Uta stansburiana lizards. In this species, males of different color-related mating strategies increase their fitness when they are less frequent, following a rock-paper-scissors dynamic (Sinervo and Lively, 1996;Bleay et al., 2007). Each of the three existing colors is associated with a mating strategy that is in advantage or disadvantage in relation with the others. Another example of NFDS between different mating strategies comes from swordtail fish (Rios-Cardenas et al., 2018). Two strategies were observed in this species: males with large blue bodies and yellow fins that display a great courting behavior; and sneaker males, completely yellow or completely blue, that mate females without consent using coercive behavior. The higher predation faced by courting males results in an interaction between opposing mating success and survival, resulting in NFDS.
Littorina fabalis (formerly known as Littorina mariae; Sacchi and Rastelli, 1966;Reid, 1996) is a small intertidal gastropod distributed across NE Atlantic rocky shores ranging from Portugal to the Barents Sea (Reid, 1996;Rolán-Alvarez et al., 2015a). It can be found grazing upon epiphytes on brown algae of the genus Fucus (Williams, 1990). It is gonochoric and has direct development (oviparous), with females laying a mass of fertilized eggs on the algae where they develop until they emerge as crawlaways (i.e., miniature adult snails). The species show a higher interpopulation genetic differentiation than the sibling species L. obtusata, presumably due to frequent bottlenecks during winter (more frequent storms in this season detach the thalluses of the Fucus sp. algae from the rock causing great mortality) and therefore pointing to a relatively low effective population size in the most exposed shores (Rolán-Alvarez et al., 1995). Other important selective agents in L. fabalis are crabs (e.g., Carcinus maenas), intertidal fish (e.g., Lipophrys pholis) and marine birds (Reimchen, 1979(Reimchen, , 1989Rolán-Alvarez et al., 1995). Among these, intertidal fishes are more likely to be visual predators than either crabs or birds (see Ekendahl, 1998 andEkendahl, 2002 for experiments of predation). Mating occurs mostly during low tide when the snails are more active. The male actively follows the mucus trails of the females (Johannesson et al., 2010;Saltin et al., 2013) and upon encountering one, it crawls over the shell of the latter, by its right side, and inserts the penis into the pallial cavity.
Populations of L. fabalis, especially those with the highest densities, often show a striking shell color polymorphism (yellow, brown, orange, and olive among other colors, see Reid, 1996;Rolán-Alvarez et al., 2015a). Predation is a popular explanation for the maintenance of this polymorphism and has been thoroughly investigated in the past with attention to crypsis and apostatic selection. For example, Reimchen (1979) found that the blennid Lipophrys pholis (previously named as Blennius pholis) actively predates those colors that are conspicuous in relation to the background. In L. saxatilis and L. obtusata this was further supported (Johannesson and Ekendahl, 2002) while it was shown that crabs could not account for visual predation (Ekendahl, 1998). Crypsis and apostatic selection have also been thoroughly studied in mangrove littorinids with contrasting evidence of crabs being visual predators (Parsonage and Hughes, 2002;McKillup and McKillup, 2008).
One population from NW Spain, Abelleira, has been monitored for more than 20 years and shows a rather stable shell color polymorphism over this period (Rolán-Alvarez and Ekendahl, 1996;Rolán-Alvarez et al., 2012, 2015b. Furthermore, disassortative mating between yellow and brown individuals and sexual selection favoring brown males have been observed in this population. However, these studies did not explicitly test if a frequency-dependent selection process is behind this stable shell color polymorphism (see Figure 1), or if genetic drift alone could explain the fluctuations in color frequencies across years. Disassortative mating caused by mate choice is expected to generate frequency-dependent sexual selection (Pusey, 1987;Takahashi and Hori, 2008;Hedrick et al., 2016Hedrick et al., , 2018. In addition, disassortative mating can also be generated as a mechanism to avoid inbreeding and its consequent negative effects. Thus, disentangling between these two possible causes of disassortative mating is crucial to understand the role of frequency-dependent sexual selection in the observation of NFDS patterns. Here, we investigate the mechanisms behind the maintenance of color polymorphism in L. fabalis by analyzing data collected across 7 years from one Spanish population (Abelleira), complemented by 1-year data from two more septentrional populations located in the White Sea. More specifically, we pursue the following objectives: i) to assess the contribution of genetic drift to explain the observed frequency changes, as theoretical models indicate that the maintenance of polymorphisms could be due to stochastic processes; ii) to test for a NFDS pattern attending to the sexual component of fitness; iii) to use a multimodel inference approach based on information theory to investigate whether mate choice, mate competition, or an interaction of both could explain the observed mating frequencies; and iv) to evaluate whether inbreeding avoidance could explain the evolution of disassortative mating in this species.

Sampling Design
Sampling involved collecting random mating pairs along the four closest unmated individuals of L. fabalis (within 10 cm of diameter; Figure 2). Collections were performed within the Fucus sp. canopy area in the lower fringe of the intertidal, following the approach implemented in previous studies (Rolán-Alvarez et al., 2012, 2015b. The mating pair and the closest unmated individuals are defined in this study as a microarea (number of matings from Although not conclusive, finding similar results would suggest that the same evolutionary forces are acting upon the color polymorphism across the species distribution range. The Spanish population was sampled once a year over a period of 7 years (2011)(2012)(2013)(2014)(2015)(2016)(2017), while the Russian population was sampled in the two localities in 2015. All collections were performed in summer, since the number of adults and mating pairs increases during this season (Rolán-Alvarez and Ekendahl, 1996). Once in the laboratory, the snails were sexed and the shell color was determined by visual inspection. A printed color model was used to classify them into four colors (brown, olive, yellow and orange) following previous studies (Rolán-Alvarez et al., 2015b). In subsequent analyses we grouped these four colors into two categories, dark (brown, olive) and light (yellow, orange), as the raw spectral reflectance data from L. fabalis shells from Abelleira showed consistent similarities within these categories (Rolán-Alvarez et al., 2015b). Olive and orange were the least frequent in Abelleira, while orange and brown were the least frequent in White Sea.

Estimation of Effective Population Size and Inbreeding
The samples collected in Abelleira (NW Spain) in 2014 (N = 564; both the mating pair and unmated individuals) were genotyped for 8 microsatellite loci developed for L. fabalis (Carvalho et al., 2015(Carvalho et al., , 2016. DNA extraction was performed following Galindo et al. (2013) and details are given in the Supplementary Material (see Supplementary Table 1).
The effective population size (N e ) was estimated from the microsatellite genotypes using the linkage disequilibrium method implemented by NeEstimator v.2.01 (Do et al., 2014) with a minor allele frequency of 0.05. This method estimates the correlation between each pair of allele frequencies using an unbiased estimation of linkage disequilibrium (Weir, 1979), and then the average correlation across loci to estimate N e (see Waples and Do, 2008). The N e value was subsequently used to test whether changes in allele frequencies for shell color could be explained by genetic drift alone. This assumed a model of one locus and two alleles with dark color being dominant, as suggested for the related species L. saxatilis by Johannesson and Butlin, 2017). The test was done by estimating the variance in the shift of allele frequencies in one generation, Wright (1922), where p 0 and q 0 are the allele frequencies at generation 0. Since L. fabalis is an annual species (see Williams, 1992) and collections in Abelleira were done across 7 years, this estimate was calculated six times, one for each annual interval. Notice that we refer here to the population of Abelleira where the species is annual, unlike snails from Russian populations, whose cycle is longer. The null hypothesis (i.e., frequency changes in a trait are caused by genetic drift) would not be rejected if the observed values of the allele frequencies across the studied years in Abelleira fell within the 95% confidence interval of the In order to test whether disassortative mating evolved as a response to inbreeding, molecular coancestry coefficients (molecular similarity sensu Walsh and Lynch, 2018, p. 693) were calculated between all pairs of individuals with Metapop2 (López-Cortegano et al., 2019). The average coancestry between individuals involved in mating pairs was compared with the average pairwise coancestry between all six individuals collected within the microarea (mating pair + 4 unmated individuals, Figure 2), and the significance of the difference between these was evaluated by a randomization one-way ANOVA (9999 permutations; Peres-Neto and Olden, 2001). If the mating pair coancestry was significantly lower than the microarea estimate, this would be an indication that disassortative mating evolved to avoid inbreeding.

Estimation of Assortative Mating
Assortative mating was calculated using the I PSI estimator, which has been shown to be robust in cases of populations with variable frequencies of the traits or datasets with small sample sizes (Rolán-Alvarez and Caballero, 2000;Pérez-Figueroa et al., 2005). The I PSI was calculated for shell color (L-light, D-dark) with the number of mating pairs from the sampled microareas as where each PSI is the number of observed mating pairs divided by the expected number under random mating for each possible combination. I PSI takes values from −1 (maximum disassortative mating) to 1 (maximum positive assortative mating) with 0 equaling absolute panmixia. The estimation of assortative mating in the wild may be biased due to the scale of choice effect (SCE). This effect biases the value of assortative mating (typically toward positive values) due to a mismatch between the sampling scale and the real scale of mate choice in nature (Rolán-Alvarez et al., 2015b). Therefore, to minimize SCE bias, the estimates of I PSI were obtained by pooling mating pairs into homogeneous color groups based on the light color frequency within their microareas: 1) Light group (≥ 66%), 2) Mixed group (33-66%), 3) Dark group (≤ 33%). The I PSI values from these groups were weight-averaged (I PSI ) each year to reduce the effect of samples' size differences. Significance for I PSI values was assessed through bootstrapping (10000 iterations) using JMATING V. 1.0.8 (Carvajal-Rodríguez and Rolán-Alvarez, 2006), while significance for I PSI was estimated using a one sample t-test (null hypothesis I PSI = 0).

Estimation of Sexual Fitness
For each year, sexual fitness was estimated in Abelleira and both locations in the White Sea by comparing shell colors before (unmated individuals) and after (mated individuals) selection (Arnold and Wade, 1984;Rolán-Alvarez and Caballero, 2000;López-Cortegano et al., 2020). The individuals belonged exclusively to one of these two groups. Thus, mated and unmated individuals were statistically independent. We carried out two different statistical methods. First, we used logistic regression to study how variation in shell color (independent variable) could determine fitness (dependent variable). In this case, fitness was considered a binary variable with a value of 0 (unmated individuals) or 1 (mated individuals). This method is commonly employed when fitness is binary and the character under study is a quantitative trait (Brodie et al., 1995;Janzen and Stern, 1998). Significance was evaluated through the Wald test (Sokal and Rohlf, 1995, p. 770). Second, we used the cross-product fitness estimator (W), within sex and sample. This is the most suitable estimator for qualitative traits (Knoppien, 1985;Rolán-Alvarez and Caballero, 2000;Martínez-Rodríguez et al., 2019), and has been shown to estimate longterm 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 studying frequency dependent mechanisms (Martínez-Rodríguez et al., 2019). As an example, the W estimator for dark individuals (D) was calculated as where AS and BS are the frequency of individuals after and before selection (mated and unmated, respectively) for D (dark) and L (light) colors. The W-values and their significance were assessed through bootstrapping (after 10000 iterations) using the software JMATING, by resampling the observed color frequencies before (BS) and after (AS) selection independently. Linear regression of the frequency of a color on sexual fitness (given as the logarithm of W in base 10) was carried out to find potential patterns of NFDS (i.e., that fitness decreases with frequency). A single linear regression ANOVA was used to assess significance. Parametric statistical tests were implemented in SPSS Statistics v.22.0 (IBM Corp.).

Contribution of Mate Choice and Mate Competition by Multimodel Comparison
We used multimodel inference (see Burnham et al., 2011) to estimate the contribution of mate choice and mate competition to the observed patterns of shell color polymorphism. These analyses were carried out with the software InfoMating v 0.2 (Carvajal-Rodríguez, 2020), which estimates the mutual mating propensities or MMPs (i.e., the tendency to mate of each particular mating combination: LL, LD, DL and DD) expected from the frequency of mating pairs and unmated individuals under each specific model tested (Carvajal-Rodríguez, 2018). The Akaike information criterion corrected for low sample size (AICc) was used to determine the weight of each model to select the best one (Burnham et al., 2011). Since the MMPs of each model are defined by mate choice and/or mate competition parameters (see Figures 2-4 in Carvajal-Rodríguez, 2020), this allowed us to identify which type of parameter is contributing the most to the observed frequencies. We analyzed Abelleira and the White Sea pooling data from the different years and populations, respectively. Additionally, a second analysis was made for pools of homogeneous groups in each population. Given the number and complexity of the obtained models (see Supplementary  Table 2), we classified them into mate choice, mate competition and mixed models, attending to the parameters included in each. For example, mate competition parameters were assumed when the parameter is present in either the whole column or row from Supplementary Table 2, otherwise the parameter was assumed to be mate choice.

Role of Genetic Drift and Inbreeding
The estimated mean effective population size (N e ) in Abelleira was 600, with 95% confidence interval of 394 and 1076. These values were used to calculate the expected variance in the change of allele frequency over one generation (from 1 year to the next) assuming only genetic drift (Figure 3). Using the mean N e , the allele frequency change between some years could not be explained by drift alone (2011, 2013, and 2014). Similar results were observed with the 95% confidence interval values of N e (Figure 3).
The average molecular coancestry coefficient for all mating pairs was 0.345 ± 0.078 (mean ± SD) and the overall average for the pairwise comparisons within microarea (including mated and unmated individuals) was 0.346 ± 0.044. Thus, mated individuals did not show significant differences in relatedness compared to surrounding individuals (ANOVA, F 1 , 188 = 0.03; p = 0.859).

Assortative Mating
All I PSI estimates obtained from the studied populations were negative (except the dark homogeneous group in location 2 of the White Sea; Table 1), indicating a higher frequency of mates between individuals with different shell color (heterotypic) than  = 0.667) and males (dotted; R 2 = 0.004). Notice that, since the W estimator for a given morph is relative to the other morph, using it for one or other morph will yield the same results.
with the same color (homotypic) in both populations. Estimates of I PSI for each of the three homogeneous groups were variable, with values ranging from 0.00 to −0.59. Only 9 (6 in Abelleira and 3 in Russia) out of 27 I PSI estimates were high in absolute value and significant under bootstrapping (p ≤ 0.05). This is likely due to the substantial reduction in the number of individuals when they are grouped in homogeneous units. After this grouping, the estimates in Abelleira varied from −0.18 to −0.40 with an overall mean of −0.29, which was significant (t 2 = 7.07; p = 0.019), pointing to consistent disassortative mating in this population. Average values in the White Sea were not significant but were very similar to those obtained in Abelleira, thus revealing the same pattern for unrelated populations. This suggests that disassortative mating may be widespread across the two areas.

Sexual Fitness Estimation
Estimates of sexual fitness using the cross-product fitness estimator (W) and the Wald test for each sex are shown in Table 2. Sexual fitness for light males was always lower than that for dark males (with the exception of White Sea 1), but the difference was only significant for Abelleira in 2011, 2012, and 2014. In the remaining years, light males' fitness was only slightly lower than that of dark males. On the contrary, light females' sexual fitness tended to be larger compared to that of dark ones, with the exception of White Sea 2 and Abelleira in 2012 (P = 0.002). The cross-product estimator had more statistical power and showed more significant cases than the Wald test under logistic regression. The results were, however, similar from both approaches. Figure 4 shows the regression of sexual fitness (as log 10 of W) of light color on the percentage of light color individuals. Notice that the result for the dark color is the inverse but would yield the same observations and conclusions. Importantly, a significant pattern of NFDS was observed for females (F 1 , 7 = 14.00; P = 0.007) but not in males (F 1 , 7 = 0.03; P = 0.872; Figure 4).

Mechanism Inference by Multimodel Comparison
The relative weight of the tested models, and the mutual mating propensities for the best-fit model for each population are presented in Figure 5 (Abelleira and White Sea). In Abelleira, models including a parameter of mate choice and mate competition accounted for 66% of the total weight, while models including only mate choice accounted for 29%. Within homogeneous groups these results were very similar except in the light group where models with only mate choice parameters accounted for 69% of the total weight. In the White Sea, both types of models (mate choice and mixed) presented similar weights. While dark and mixed homogeneous groups showed mate choice models to account for a major part of the total weight (62 and 77%, respectively), mixed models accounted for a 66% of the total weight in the light group.
The best-fit model in Abelleira included both mate choice and mate competition parameters. The estimate of the mate choice parameter was 0.85 ± 0.069 (Figure 5). Thus, mating between homotypic pairs (DD or LL) had 15% less propensity than mating between heterotypic pairs (DL or LD). This preference decreased to 0.55 ± 0.01 within homogeneous groups (a 45% reduction; Figure 5). The estimate of the mate competition parameter was 1.31 ± 0.106 (Figure 5), which implies a sexual advantage of 31% for dark males over light ones, and 1.45 ± 0.23 within homogeneous groups (Figure 5). Results for the light homogeneous group differed somewhat as no mate competition parameters were defined in the best-fit model, but rather a choice component was found in favor of heterotypic mating pairs and against LL ones.
The best-fit model for the White Sea was defined by mate choice only. The estimated mate choice parameter for this model was 0.77 ± 0.06 (a reduction of 23% in the mating propensity of homotypic pairs). The best-fit models within homogeneous groups were somewhat different. The dark group showed the same pattern as in Abelleira, with a mate choice parameter of 0.42 (a reduction of 58% in the mating propensity of homotypic pairs). The mixed group showed two different mate choice parameters. The first parameter reduces the mating propensity of DD mating pairs, while the second increases it in mating pairs of the DL kind (first letter corresponds to the male color). The light group showed that dark males had 51% less propensity to mate and there was a choice component against LL mating pairs.

DISCUSSION
Littorina fabalis color polymorphism has been shown to be stable across years in the current and previous studies (Figure 1; Rolán-Alvarez and Ekendahl, 1996;Rolán-Alvarez et al., 2012). This is in spite of frequent bottlenecks during winter, at least in NW Spain (Rolán-Alvarez et al., 1995), which could aggravate the effects of genetic drift causing allelic fixation. Here we have presented several examples that support the contribution of a mechanism of NFDS for the maintenance of this polymorphism across the species distribution. Our results indicate that such a mechanism depends on precopulatory mate choice. We have also found some results pointing to directional sexual selection in dark male color. However, such directional selection would also cause a fixation of alleles in the long run and thus its contribution to polymorphism maintenance is not likely. Possible alternative scenarios and implications of the present study are discussed below.

The Stochastic Scenario
Genetic variation is expected to be quickly eroded under stochastic forces when the populations have low effective sizes (Cook, 1992;Hoffman et al., 2006). L. fabalis is a species which suffers strong bottlenecks during winter, possibly due to the detachments of Fucus thalluses during storms (Williams, 1996). This could cause severe reductions in effective size (Rolán-Alvarez et al., 1995) and therefore also presumably in color frequencies. However, in the 7 years of study in Abelleira, frequency fluctuations do not indicate such severe reductions (Figure 1). Therefore, we checked whether the observed changes in allele frequencies could be just explained by genetic drift without the need to claim for any selective mechanism.
Our analysis of genetic drift could be questioned by the facts that the genetic basis of color in Littorinids, including L. fabalis, is not yet completely resolved (Reimchen, 1989;Ekendahl, 1995;McQuaid, 1996;Kozminsky, 2014;Johannesson and Butlin, 2017), and because we have combined all the observed colors in just two categories. Our decision for this latter was based on phenotypic similarity (see Rolán-Alvarez et al., 2015b) and because orange and olive colors had a low frequency in all years of study. However, our assumption that shell color is FIGURE 5 | Mutual mating propensities obtained for the best models in Abelleira and the White Sea in both total numbers and homogeneous groups, pooled across all years (Abelleira) and populations (White Sea). Male in columns and female in rows. Values within boxes are the best model mutual mating propensities for the four different mating combinations. Pie charts represent the sum of weights for mate choice, mate competition, mixed and random mating models (see Supplementary  Table 2). coded by just one locus is conservative, because the variance in gene frequency across years would be further diminished considering two or more loci under different inheritance models (see simulations by Estévez, 2018, Chapter 3). Our results show that, although genetic drift cannot be discarded when fluctuations in color frequencies are small, genetic drift alone cannot explain the maintenance of the polymorphism over time in our study population. This is in agreement with findings in L. saxatilis where genetic drift cannot explain the persistence of rare colors (Johannesson and Butlin, 2017). In contrast, previous results in Littoraria pallescens showed genetic drift to be the main driver of polymorphism maintenance in the absence of selection (Cook, 1992). However, in other mangrove littorinids, heterogeneous natural selection has gained support as a more likely explanation for such a maintenance (Parsonage and Hughes, 2002). Our results also suggest that there is a frequency threshold over which selection acts to prevent allele fixation. That is, when a morph's frequency surpasses a boundary, selective mechanisms occur that outrun the relative strength of drift. This has also been proposed in the candy-stripe spider (Oxford, 2005) and points to the alternating roles of selection and drift in maintaining polymorphisms (see O'Hara, 2005).
Overall, this suggests that genetic drift is less important than previously thought in this population and species. There may be two reasons for this. First, these populations are dense and stable with a limited population reduction during winter (personal observation). Second, in the related species L. saxatilis, quick recovery of the populations' genetic diversity has been observed after bottleneck episodes (Johannesson and Johannesson, 1995;Piñeira et al., 2008). Surprisingly, the rarest colors in L. saxatilis persisted over time in spite of strong population reductions (Johannesson and Butlin, 2017). The fact that pregnant females in this species carry offspring from multiple males may prevent allele loss by genetic drift (Panova et al., 2010). We have no direct evidence of such a process in L. fabalis but its evolutionary proximity to L. saxatilis, and the fact that females mate multiple times, suggests that they may spawn multiple sired offspring. This has in fact been studied in the sibling species L. obtusata where between 4 and 6 sires were found to contribute to females' broods (Paterson et al., 2001). Therefore, it is likely that each L. fabalis female spawns a genetically diverse brood, and consequently, this can contribute to maintaining a larger effective size despite adult population fluctuations.

Balancing Mechanism Maintaining the L. fabalis Color Polymorphism
Even though there are some studies showing evidence of disassortative mating (Takahashi and Hori, 2008;Hedrick et al., 2016Hedrick et al., , 2018, this is a rare occurrence in natural populations and has been associated with type I errors (Jiang et al., 2013). Previous studies hypothesized that such a lack of cases could be due to the absence of spatial and temporal dynamics incorporated in the analyses (Rolán-Alvarez et al., 2015b;Fargevieille et al., 2017). Here we have accounted for the scale bias effect (SCE; see Rolán-Alvarez et al., 2015b) and have shown that disassortative mating was observed in all the studied years and in populations at both extremes of the species distribution. Preliminary laboratory evidence also indicates disassortative mating in this species and for the same colors from Abelleira (I PSI = −0.29 ± 0.102; p = 0.0054; unpublished results). In addition, two other NW populations of L. fabalis from Spain show a high and negative assortative mating for shell color (unpublished results). Thus, L. fabalis provides a good example of consistent disassortative mating for shell color.
Negative mate choice for color in L. fabalis (see below) could be the cause of disassortative mating, as observed by experimental trials in various study systems. For instance, in the butterfly Heliconius numata, experimental trials showed a tendency for disassortative mating which was hypothesized to maintain wing color polymorphisms in spite of positive frequency dependent selection by predators (Chouteau et al., 2017). Another well-known example is the white-throated sparrow, where disassortative mating is almost absolute in the wild (Hedrick et al., 2018) and experimental lab evidence supports negative mate choice (Houtman and Falls, 1994). However, as has been observed in several studies (see for example López-Cortegano et al., 2020, and references therein) results from laboratory trials cannot be readily extrapolated to natural populations. Therefore, one advantage of our study is that it shows what is going on directly in the field, and so it has an evolutionary relevance. On the contrary, a problem of field studies is that they can be affected by several source of bias (Alatalo et al., 1986;Crespi, 1989;Taborsky et al., 2009). Nevertheless, we provided estimates corrected for the SCE bias and we also estimated the strength of mate choice and sexual selection effects by independent methods, and both confirm the present interpretation.
First, our best fit models ( Figure 5) indicated negative mate choice in almost all cases (the diagonal for heterotypic pairs showed greater values than that of homotypic pairs). Second, the trend observed in Figure 4 for female sexual fitness strongly showed a pattern of NFDS, which is an expected outcome under male mate choice with preference for the dissimilar color (see Kokko et al., 2007;Wellenreuther et al., 2014). This is in accordance with existing evidence which shows that mate choice in Littorina is indeed preferentially exerted by males (Erlandsson and Kostylev, 1995;Ng et al., 2013Ng et al., , 2019Saltin et al., 2013). Thus, for L. fabalis there seems to be a pattern of NFDS in which the most frequent female color mate less often. Taken together, these pieces of evidence indicate the occurrence of mate choice for dissimilar color in L. fabalis. Mate competition could be affecting our results but our multimodel inference accounted for its contribution on the observed mating propensities. The multimodel inference results did point to a greater mating advantage in dark males, at least in Abelleira (up to a 31%), suggesting that dark males outcompete light males. Such occurrence, however, only explains why there is a greater frequency of mating pairs where a dark male is involved but not why there is a greater frequency of heterotopic pairs. In other words, that dark males have greater sexual fitness does not explain the evidence that points to negative mate choice.

Potential Causes for Observed Patterns
A common reason to mate preferentially with dissimilar morphs is to avoid inbreeding, which can cause fitness disadvantages in the offspring (Keller and Waller, 2002;Charlesworth and Willis, 2009;Duthie and Reid, 2016). Thus, patterns of disassortative mating are not rare in cases where homozygotes are at selective disadvantage (Tregenza and Wedell, 2000;Mays and Hill, 2004). Here, we tested if relatedness of mating pairs were significantly lower than that of unmated individuals from the same microarea. We did not find significant differences and our results do not support that overall preference for dissimilar morphs have evolved to avoid inbreeding, at least in the studied population. Previous studies in L. saxatilis yielded the same observation (Ng and Johannesson, 2016). However, the level of relatedness at some microsatellite loci does not necessarily reflect the relatedness for the whole genome. Mating can therefore respond to differences in certain genomic regions. For example, in the white throated sparrow, homozygotes for the chromosome allele 2 m suffer low fitness because of the exposition of deleterious alleles in that chromosome (Tuttle et al., 2016). Additionally, in this species, disassortative mating between morphs is almost absolute (Hedrick et al., 2018) and explains the polymorphism maintenance. Both pieces together indicate that negative choice may have evolved to reduce the number of less fit homozygous 2 m types in this bird species. In an analogous way to the whitethroat sparrows, the observation that light mating pairs is less frequent than the remaining mating combinations could suggest that homozygosity for light color may have deleterious effects in L. fabalis. A thorough analysis of the genetic basis of shell color in L. fabalis is fundamental for studying this hypothesis.
Shell color might not be the trait under direct choice but rather a correlated trait. Association between color and other traits is typical in nature (McKinnon and Pierotti, 2010). Sacchi (1969Sacchi ( , 1972 proposed several such associations including salinity and desiccation tolerance. Association between these traits can be due to genetic inversions which can originate linkage disequilibrium between locally adapted loci. Indeed, L. fabalis studies supported the existence of an inversion which affected size, the allozyme Ark Ginase, and a RAPD locus genotype (Johannesson and Mikhailova, 2004;Kemppainen et al., 2011). This inversion would create a locally adapted supergene codifying for several traits. In this scenario, the occurrence of disassortative mating would benefit the occupation of a wider variety of habitats (van Nes, 1998). Of the possible traits associated with color, variation at loci determining the immune system is a promising candidate. In invertebrates in general, color is usually linked to immune responses (see Scheil et al., 2013 and references therein). In turn, mating and immune characteristics have been found to be greatly linked (Lawniczak et al., 2007). In L. fabalis, male mate choice for different color could have evolved as a side effect of mate choice for individuals with dissimilar immune characteristics, in the same fashion as the MHC complex in vertebrates (e.g., Consuegra and Garcia de Leaniz, 2008;Huchard et al., 2013). In this case, the benefit would be to produce more fit heterozygous offspring for immunity genes. It is still uncertain whether male choice can evolve to improve offspring quality. While this has been suggested in house mice, Mus domesticus (Gowaty et al., 2003), other studies testing this hypothesis have found either a negative (Swierk and Langkilde, 2019) or absence of correlation (Bertram et al., 2017).
Male mate choice can also evolve to exploit other fitness benefits. In many cases, males prefer those females that denote greater fecundity (e.g., Edward and Chapman, 2012;Jaworski et al., 2018). In Littorina, an honest indicator of fecundity is female size and males tend to mate with larger females (Saltin et al., 2013, and references therein). Color can also be an indicator of female condition and thus inform males to increase fitness (Amundsen and Forsgren, 2001;Weiss, 2006). We have no evidence that male mate choice for dissimilar colored females is related to increase fecundity in L. fabalis. To explore this, an experimental approach could be performed where numbers of eggs resulting from homotypic and heterotypic mating pairs are compared. A greater number of eggs from heterotypic mating pairs would indicate male mate choice to increase fecundity in this species.
The observed greater sexual fitness in dark males may be driven by temperature and risks to desiccation, given the known relationship between shell color and sunlight incidence (see Heath, 1975;Cook and Freeman, 1986). Shell color in Littorina could be in fact associated to thermal physiology (reviewed by Rolán-Alvarez et al., 2015a). For example, evidence that sunlight is acting differently on dark and light shells has been claimed in the sibling species L. obtusata (Sergievsky, 1992;Phifer-Rixey et al., 2008) and other gastropods (Cepaea nemoralis, Heath, 1975;Littoraria pallescens, Cook and Freeman, 1986;Nucella lapillus, Etter, 1988). An increment in temperature in dark males may increase their sexual fitness. Elevations in sexual selection intensity due to temperature have in fact been shown in other organisms (Kvarnemo, 1996;Allen et al., 2012). Alternatively, an increase in temperature may trigger behaviors to avoid desiccation (see Miller and Denny, 2011). In our study system, it may cause dark males to move more actively for shelter and thus increase the rate of encounter with females. Interestingly, such sexual fitness advantage is observed in dark males in Abelleira (best models in total results in Figure 5), which represents the southern limit of L. fabalis distribution and where the temperatures are higher, but not in the White Sea, in agreement with the temperature hypothesis.

CONCLUSION
Shell color polymorphism in L. fabalis is persistent despite expected strong genetic drift. We have otherwise found that color polymorphism in this species follows a NFDS pattern when studying the sexual component of fitness. In particular, this pattern was only present in females. This, together with the occurrence of disassortative mating, is indicative of male mate choice, which is common in species of the Littorina genus. Since these patterns result from observations in the wild, we used multimodel inference to assess the possible role of mate choice and/or mate competition. The inference strongly supports mate choice but also indicates mate competition in males. Several hypotheses have been presented here that can account for either mate choice or mate competition. However, the connection, if any, between them is still to be elucidated and needs further research. In particular, understanding the underlying genetics of color and what is the correlation between different traits, could cast light into these mechanisms and explain the underlying biological causes. Experimental trials in the lab for both mate choice and mate competition can provide further information concerning these biological causes. This study has confirmed a rare pattern, disassortative mating, to occur in very distant populations of L. fabalis, and opens the door for further research into understanding the biological mechanisms behind polymorphism maintenance and mate choice in marine systems.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: doi: 10.6084/m9. figshare.13050494.v1 and doi: 10.6084/m9.figshare.13050509.v1.

AUTHOR CONTRIBUTIONS
DE performed samplings in Galicia, analyzed the data, designed the tables and figures, wrote and reviewed the content of the manuscript, and designed it for publication. EK performed samplings in the White Sea, provided the data from that location, and reviewed the content of the manuscript. AC-R provided software resources for the analyses of data and reviewed the content of the manuscript. AC contributed to the analysis of genetic drift and reviewed the manuscript. RF reviewed the manuscript, added valuable criticism toward methods, and contributed to improve the writing. JG and ER-A made equal contributions to the manuscript, performed samplings in Galicia, contributed to the writing of the manuscript, and reviewed it for publication. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by Xunta de Galicia (ED431C 2016-037), FONDOS FEDER (unha maneira de facer Europa) and Ministerio de Economía, Industria y Competitividad, Agencia Estatal de Investigación (CGL2016-75904-C2-1-P). JG was funded by a JIN project (Jóvenes Investigadores sin vinculación o con vinculación temporal, Ministerio de Ciencia, Innovación y Universidades, RTI2018-101274-J-I00). RF was financed by the FEDER Funds through the Operational Competitiveness Factors Program -COMPETE and by National Funds through FCT -Foundation for Science and Technology though a research contract within the scope of the project "Hybrabbid" (PTDC/BIA-EVL/30628/2017 and POCI-01-0145-FEDER-030628). UVigo Marine research centre is one of the five research centers from all fields of knowledge recognized as "Galician Singular Research Centres" by the Galician Regional Government through the "Excellence in Research (INUGA)" Program from the Regional Council of Culture, Education and Universities. The Centre's activity is co-funded by the European Union through the ERDF Operational Program Galicia 2014-2020. The research work of EK was carried out as part of the State Task RF (no. AAAAA19-119022690122-5).