A Simulation Study of the Ecological Speciation Conditions in the Galician Marine Snail Littorina saxatilis.

In the last years, the interest in evolutionary divergence at small spatial scales has increased and so did the study of speciation caused by ecologically based divergent natural selection. The evolutionary interplay between gene flow and local adaptation can lead to low-dispersal locally adapted specialists. When this occurs, the evolutionary interplay between gene flow and local adaptation could eventually lead to speciation. The L. saxatilis system consists of two ecotypes displaying a microhabitat-associated intraspecific dimorphism along the wave-exposed rocky shores of Galicia. Despite being a well-known system, the dynamics of the ecotype formation remain unclear and cannot be studied from empirical evidence alone. In this study, individual-based simulations were used to incorporate relevant ecological, spatial, and genetic information, to check different evolutionary scenarios that could evolve non-random mating preferences and finally may facilitate speciation. As main results, we observed the evolution of intermediate values of choice which matches the estimates from empirical data of L. saxatilis in Galician shores and coincides with previous theoretical outcomes. Also, the use of the mating correlation as a proxy for assortative mating led to spuriously inferring greater reproductive isolation in the middle habitat than in the others, which does not happen when directly considering the choice values from the simulations. We also corroborate the well-known fact that the occurrence of speciation is influenced by the strength of selection. Taken together, this means, also according to other L. saxatilis systems, that speciation is not an immediate consequence of local divergent selection and mating preferences, but a fine tuning among several factors including the ecological conditions in the shore levels, the selection strength, the mate choice stringency, and cost to choosiness. The L. saxatilis system could correspond to a case of incomplete reproductive isolation, where the choice intensity is intermediate and local adaptation within the habitat is strong. These results support previous interpretations of the L. saxatilis model system and indicate that further empirical studies would be interesting to test whether the mate choice mechanism functions as a similarity-like mechanism as has been shown in other littorinids.

Introduction

Speciation caused by ecologically based divergent natural selection can occur at small spatial scales. Microgeographic adaptation occurs when neighbor groups of individuals have adaptively diverged (Richardson et al., 2014) and this adaptation at a small spatial scale may happen despite the high potential for mixing within the dispersal area. This idea has produced controversy between theoretical studies that show that natural selection can overcome the effects of gene flow (Dieckmann and Doebeli, 1999; Gavrilets, 2004; Barluenga et al., 2006; Hey, 2006; Bolnick and Fitzpatrick, 2007; Rolán-Alvarez, 2007; Debarre and Gandon, 2011; Weissing et al., 2011; Savolainen et al., 2013; Butlin et al., 2014; Richardson et al., 2014; Getz et al., 2016; Rettelbach et al., 2016; Servedio and Boughman, 2017; Kopp et al., 2018) and those that admit that this is feasible but that the conditions necessary for it to occur make it difficult to observe in the real world. (Jiggins, 2006; Babik et al., 2009; Meyers et al., 2012; Martin, 2013; Savolainen et al., 2013; Martin et al., 2015; Fernández-Meirama et al., 2017a; Foote, 2018).

There are various plausible scenarios for the evolutionary interplay between gene flow and local adaptation. They may give rise to a monomorphic population, or to generalists adapted to different habitats, or to polymorphic subpopulations of locally adapted specialists, which, in case of low gene flow could lead to speciation (Coyne and Orr, 2004).

However, for local adaptation to occur, a gene flow-reducing mechanism, such as non-random mating due to mate choice, may be necessary in addition to selection. Furthermore, to detect mate choice, it is important to take into account the true spatial scale at which mate choice occurs to avoid a version of the so-called Simpson’s paradox in which the mating pattern existing in some data sets will disappear or be reversed when the groups are combined. Thus measuring mate choice at the correct scale is the key to understanding this evolutionary process (Rolán-Alvarez et al., 2015b; Estévez et al., 2018).

Local adaptation and speciation through non-random mating are expected to occur when the traits involved in the divergent selection and mate choice are the same, the so-called magic trait (Gavrilets, 2004; Thibert-Plante and Gavrilets, 2013; Servedio and Boughman, 2017; Richards et al., 2019).

Because of the long time scale necessary for local adaptation and mate choice evolution under gene flow, simulation has become an essential part of speciation research. Recent efforts on the study of the evolution of reproductive isolation focused on the complex interactions that emerge in the presence of local adaptation, mate choice, and sexual selection under different genetic structures, migration models, etc. (Sadedin et al., 2009; Thibert-Plante and Hendry, 2009; Thibert-Plante and Hendry, 2011; Cotto and Servedio, 2017; Sachdeva and Barton, 2017; Perini et al., 2020).

The present work is based on L. saxatilis as a model organism and used individual-based simulation to identify some relevant parameters, such as local selection intensity, hybrid zone environment, and demography that can influence the evolution of local adaptation and non-random mating preferences. The L. saxatilis model case relies on the existence of two main ecotypes that have evolved by natural selection to adapt to distinct microhabitats in different geographical regions in parallel (see next section). The two ecotypes have evolved different average sizes, morphologies, behavior, and physiologies, and may potentially show partial reproductive isolation from each other due to different microhabitat and mating preferences (Rolán-Alvarez, 2007; Johannesson et al., 2010; Rolán-Alvarez et al., 2015a). This system is of particular interest in the Galician populations where these ecotypes live at different shore levels, while at the mid-shore they meet and partially hybridize producing a hybrid zone that has been suggested as particularly suitable for ecological speciation in presence of gene flow (Rolán-Alvarez, 2007; Butlin et al., 2014; Boulding et al., 2017). Despite the number of empirical studies that have estimated biological and ecological parameters for this model system (reviewed in Pérez-Figueroa et al., 2005; Johannesson et al., 2010; Rolán-Alvarez et al., 2015a), several questions remain unanswered (see below). Here, we perform spatially-explicit computer simulations for studying the L. saxatilis Galician system incorporating the evolution of mate choice and reproductive isolation and using the parameter information available from empirical studies. Noteworthy, there were previous efforts of modeling another L. saxatilis system, the Swedish one (Sadedin et al., 2009; Perini et al., 2020) which has also two ecotypes although under a different spatial and microhabitat distribution (see next section). Thus, although building on the Swedish 2009 modeling attempt, our model and simulations have different demographic and evolutionary settings and utilize a different mating preference function for better matching of the Galician demographic and evolutionary conditions (Rolán-Alvarez, 2007; Carvajal-Rodríguez and Rolán-Alvarez, 2014; Rolán-Alvarez et al., 2015a). The effect of introducing a cost of choosiness that was not considered in previous simulation studies was also included in our simulations. However, to validate our program we also replicate the 2009 Swedish model (see Program validation section).

Therefore, the main objective of this study is to try to answer the following questions which have been explicitly studied for the Galician system, although they can be extended to similar ecological speciation scenarios. Thus, given the conditions outlined in the Galician L. saxatilis model system (see details in Model and Methods but briefly, low dispersal ability, disruptive ecological adaptation, a magic trait, and an area of sympatry), our questions are as follows:

1. Can we expect to find locally adapted specialists in the presence of gene flow?

2. Can mate choice (choosiness) and so full reproductive isolation, be evolved as a side effect of ecological adaptation?

3. Can the intensity of mate choice be increased at the area where the ecotypes meet because of reinforcement or other related mechanisms?

4. What is the impact of considering a cost in mate choice in relation to patterns of natural and sexual selection observed in nature?

Following, we briefly review the L. saxatilis species as model organisms; then, we explain the simulation model adapted to the Galician L. saxatilis scenario to answer the questions mentioned above.

Materials and Methods

Model Organism

L. saxatilis is a marine snail gastropod living on North Atlantic rocky shores. The species crawl directly on rocks feeding on microalgae, diatoms, and detritus. The species is gonochoric, shows internal fertilization and its ovoviviparous, showing direct development, i.e., with the female having a brood pouch with dozens to hundreds of embryos that are born as miniature snails. Due to these characteristics, L. saxatilis shows low dispersal ability which also facilitates its capability to adapt to local conditions, and so local ecotypes have been frequently described on different localities and regions (Reid, 1996; Rolán-Alvarez, 2007; Johannesson et al., 2010; Rolán-Alvarez et al., 2015a). In particular, there is a striking sympatric polymorphism on Galician rocky shores (Northwest Spain), where two ecotypes adapted to two different shore levels and microhabitats, coexist partially in sympatry at the mid-shore (i.e., can frequently meet and mate partially assortatively, Rolán-Alvarez, 2007). A polymorphism determined by similar ecological forces (wave exposition and crab predation) is known in Britain and Sweden although the ecotypes are rarely (Britain) or never (Sweden) found in sympatry there (Rolán-Alvarez et al., 2015a). In Galicia, an upper-shore ecotype [named “Crab” in (Butlin et al., 2014)] lives on the barnacle belt, showing a larger and more robust shell, colored with alternate ribs and black lines, that protects against crab predation (Rolán-Alvarez, 2007; Butlin et al., 2014; Boulding et al., 2017). This ecotype also shows good resistance to desiccation, osmotic, and sun stress. In addition, the lower-shore ecotype [named “Wave” in (Butlin et al., 2014)] appears associated with the mussel belt and has a smaller, softer, and smoother shell, with a bigger aperture to accommodate a massive muscular foot. This strong foot is necessary to maintain attached snail to the substratum due to the strong waves that commonly impacts at first instance on the lower shore. There are no crabs that can predate on this ecotype on the lower shore, and so, corresponding upper and lower shore microhabitats, represent a differential selective regime for the ecotype subpopulations. At the mid-shore, where the above microhabitats overlap forming a patched distribution, both ecotypes meet and hybridize producing a number of intermediate morphological forms, although they show partial reproductive isolation with each ecotype mating preferentially with specimens of the same ecotype. These intermediate forms may usually resemble genetically to one or another ecotype though they occasionally can be truly intermediate genotypes (Galindo et al., 2013; Kess et al., 2018). The partial premating reproductive isolation seems to be linked to the size differences existing between these two ecotypes and could be a side effect of the size assortative mating typically observed in the species, being considered as a typical magic trait (Fernández-Meirama et al., 2017a; Boulding et al., 2017).

Simulation Model

There are few simulation studies on the interaction of natural selection with gene flow and the evolution of local adaptation in Littorina (Boulding, 1990; Pérez-Figueroa et al., 2005; Sadedin et al., 2009; Westram et al., 2018; Perini et al., 2020) and only the one from Sadedin and co-workers incorporates the evolution of mate choice and reproductive isolation.

The following model is a modification of the Sadedin et al. (2009) one, attending to the specific parameter values and spatial structure of the Galician L. saxatilis compared to the Swedish case. In the Swedish case, the local morphological and behavioral adaptation occurs within islands in the Swedish archipelago in a horizontal wave exposure gradient through cliff and boulder habitats interspersed along the shore (Sadedin et al., 2009). On the contrary, the Galician ecotypes vary along a vertical, within locality microgeographic wave exposure gradient with different spatially varying selection favoring large sizes in upper-shore (wave-sheltered) and small sizes in lower-shore (wave-exposed). Unlike the Swedish case, the Galician ecotypes jointly with the intermediate forms can be observed at approximately similar frequencies at the mid-shore (Rolán-Alvarez, 2007; Galindo et al., 2013; Kess and Boulding, 2019). The migration rates of these snails are relatively low (less than 2 m from the released point after 1 month), and the migration vectors of transplanted snails show a significant direction towards their native tidal height (Erlandsson et al., 1998).

The Galician simulation model began with a pregnant female well adapted to the sheltered habitat arriving at the most upper deme in the shore, which is a realistic scenario of how L. saxatilis colonize new habitats (Janson, 1987). This founder female produced offspring having optimum phenotype for the wave-sheltered shore and a random mating value for the mating trait (see below). The model was spatially explicit in one dimension with three shore levels (upper, intermediate, and lower).

Individuals were diploids with separate sexes; each individual was constituted by two different additive quantitative traits and a set of eight microsatellite loci. The first quantitative trait is an ecological magic trait x i.e., it defines the individual fitness in a specific shore level and also is the target trait for the mate choice when mating is not random. An example of this kind of trait in L. saxatilis could be size (Boulding et al., 2017). The second trait c is a mate choice trait and is related with the sign of the preference for mate choice (positive or negative assortative) and also with the male choosiness. The choosiness is defined as the absolute value of the linear function C = 2c−1, while the sign of the preference depends on the sign of C. If C = 0 (c = 0.5) there is random mating; negative values of C (c < 0.5) imply negative assortative mating, while the positive values of C (c > 0.5) generate positive assortative mating (the absolute value, choosiness, varies between 0 and 1).

Each quantitative trait (x and c) was constituted by additive unlinked bi-allelic loci (see Table 1). The alleles may have value of 0 or 1. The trait value is scaled between 0 and 1 by summing the alleles and then dividing by L. Since we did not consider environmental variation, the genotypic and phenotypic values were the same.

TABLE 1

Model parameters.

	Parameter	Symbol	Value
Demography	Generation number	T	20,000
	Carrying capacity exposed deme	K	15,000
	Intermediate and sheltered demes		3,750
	Mean offspring number	b	50
	Between deme migration probability	m
		0 demes	0.75
		1 deme	0.15
		2 demes	0.1
Genome Structure	Number of microsatellites		L = 8
	Ecological magic trait (L loci)	x	L = {4, 8}
	Male mating trait (L loci)	c	L = {4, 8}
Selection	Selection strength	σs	{0.15, 0.45, 1}
	Selection in middle habitat	θ	{no, yes}
Mating	Evaluated females per male	N f	10
	Mating preference tolerance	σa	{0.05, 0.1}
	Mate choice cost		{no, yes }
Mutation	Neutral mutation rate per locus	μ0	10–3
	Trait mutation rate per locus	μ	10–5

The generations were discrete non-overlapping and the simulation events occurred in the following order: birth, viability, mating, migration, reproduction, and mutation. The parameter values can be consulted in Table 1.

Spatial Distribution

The spatial distribution of individuals in the simulation tried to resemble the Galician rocky shore where L. saxatilis inhabits (Rolán-Alvarez, 2007). The population was simulated as a one-dimensional array with a total of 26 demes. The first 20 demes corresponded to the upper shore which is the sheltered habitat in the Galician coast, the next two demes were the intermediate habitat, and the last four corresponded to the lower shore which is the exposed habitat. Each deme had a maximum capacity depending on the carrying capacity of its habitat. A deme in the lower-shore exposed habitat had a carrying capacity of K = 15,000, while the intermediate and sheltered had K = 3, 750 each. This distribution of demes per habitat tries to match the true habitat distribution and densities observed in Galicia (Rolán-Alvarez, 2007). However, the symmetric models (with equal number of demes and population density) were also run in order to interpret the consequence of this asymmetry (see below). Only pregnant females were considered as migrants. A migrant female can move through one or two demes, the offspring will be born at the arrival deme. Each individual develops its whole life within a deme.

Viability

The survival condition of each individual depended on the distance between its ecological trait x and the optimum in the habitat where it was born. The optimum was 0, 0.5, and 1 for the exposed, intermediate, and the sheltered shore, respectively. We also considered the case of no selection in the intermediate habitat so that under this scenario a fraction of individuals may reproduce in a spatial region (two demes) that is neutral with respect to ecological speciation (Cotto and Servedio, 2017). The fitness w i of an individual i with phenotype x i was given by (Sadedin et al., 2009) where x i is the value of the ecological trait, θ is the optimum in the habitat where the individual was born, and σs is the inverse of selection strength.

The offspring survival was density-dependent following a Beverton-Holt model (Beverton and Holt, 1957). Thus, the individual viability depended on the total number of local juveniles N, the mean number (b = 50) of offspring of each individual, and the carrying capacity (Sadedin et al., 2009) where w i is the individual fitness (Equation 1) and K is the carrying capacity of the deme.

Mating

In L. saxatilis, mate choice is performed by males (Rolán-Alvarez et al., 2015a). In our simplest model (without cost), each male always mate with one of 10 females randomly chosen from its own deme. The probability of mating between the male and any of the 10 females is proportional to the FND ψ function (Ψ , see Eq 3). where Ψ depends on the value of the male mating trait c and the distance D = |x m −x f| between the ecological traits of the male and the female candidate. The parameter C is a linear function of c, C = 2c−1, and D is the maximum D (1 in this case), and σ is the tolerance or inverse strength of the male mating preference (Carvajal-Rodríguez and Rolán-Alvarez, 2014).

From (3), when c > 0.5 (C > 0) the mating is positive assortative, when c < 0.5 (C < 0) the mating is negative assortative, and it is random when the phenotype c = 0.5 (C = 0).

The formulae in (3) is a modification of the original model of (Sadedin et al., 2009) to correct for the known problem that, when phenotypic values are equal to 0.5, positive assortative mating occurs regardless of the choice value (Carvajal-Rodríguez and Rolán-Alvarez, 2014).

As already explained, the choosiness is the absolute value of C while the sign of the preference depends on the sign of C. Thus, provided that some choosiness value evolved (|C| > 0), the resulting mate choice mechanism corresponds to a matching rule model combined with a magic trait (Kopp et al., 2018).

Under the mating without cost, each male always mates and produces offspring, independently of how choosy he is, so there will be mating even if all Ψ values are very low. This is achieved by defining the mating probability as the FND function for a given male and female pair relative to the total FND functions for all females that the given male encounters.

The following is an example of how mating is implemented without cost. Consider a situation where a male i encounters female 1 and they have Ψ = 10–3, and for any of the remaining females in the deme, the pair of male i with any female j ≠ 1 has Ψ = 10–4. The sum of Ψ values for this male with the ten females he encounters is S = 10–3 + 9 × 10–4 = 0.0019. In the model without cost, the probability of mating for male i is Pr (i × 1) + 9Pr (i × j) = Ψ /S + 9Ψ /S = S/S = 1 where we divided by S due to the condition that the male must mate. Thus, we see that the probability of male i mating with female 1 is approximately 0.53 (i.e., 10−3/S) which is one order of magnitude greater than with any other female, so it is proportional to the Ψ values.

Mate Choice Cost

Given the value Ψ of the preference function, we considered also the case with cost to choosiness (Bolnick, 2004). In the model with cost, the mating probability is equal to the FND function for a given male and female pair (the coefficient of proportionality is set to 1) and therefore an overly choosy male might not find a mate. Following our previous example, male i randomly encounter female 1 plus nine females from its deme and will try to mate with all of them but now the mating probability coincides with the Ψ value and so, the probability of mating for male i is Pr (i × 1) + 9Pr (i × j ≠1) = Ψ + 9Ψ = S = 0.0019. Therefore, there is more than a 99.8% chance that this male will not mate. If no mating occurs the male is discarded and the next male is considered. If no mating takes place at all, the deme will be empty.

Migration

The life cycle occurs within the deme and ends after reproduction. The pregnant females can disperse 0 (no migration), 1 or 2 demes in either direction with probability 75, 15, and 10%, respectively. After migration, the female produces the offspring within the new deme and dies. These numbers were chosen in order to simulate a scenario of high gene flow under sympatry, as this was probably the starting point of the system before any local adaptation was evolved (Rolán-Alvarez, 2007).

Reproduction

After mating and migration, the number of offspring was obtained from a Poisson with mean b (b = 50 see Sadedin et al., 2009). Recall that males are the choosy sex so each male randomly encounters some females in the deme and successful mating will occur depending on the FND function i.e., on the phenotypes of male and female and the choosiness of the male. If a specific mating is successful then the offspring number per mated male-female pair is sampled from the Poisson distribution with a mean of 50. The genetic composition of the offspring depended on the independent segregation of the trait loci and microsatellites from the parental gametes.

Mutation

The newborn underwent mutation after reproduction. The per individual haploid locus mutation rate was 10–3 for the stepwise mutation model of microsatellites (Slatkin, 1995) and 10–5 for the additive locus (Sadedin et al., 2009).

The life cycle was allowed to run for 20,000 generations. We performed 20 replicates for each parameter combination and 48 different cases of parameter combinations for the Galician model. The number of runs for some cases was extended to 80,000 generations to evaluate long-term results, as this number of generations corresponds to an upper threshold for Galician L. saxatilis based on a generation time of 6 months and estimates of divergence time between ecotypes (Quesada et al., 2007).

Initial Conditions

As already indicated in the simulation settings, the experiment began with a founder event in which the most upper deme of the sheltered habitat was colonized by the offspring of a single female adapted to this habitat (x = 1). The offspring number was 50 (initial population size). All the individuals were homozygous for the quantitative traits and had the optimal ecological phenotype for this habitat (x = 1). The male mating trait was initially fixed at c = 0.5 (random mating). The microsatellite genotype was heterozygote for every locus and individual.

The simulations were performed on C++11 and the source code is available at https://acraaj.webs.uvigo.es/software/FernandezMeiramaetal_sourcecode.zip.

Program Validation

To validate our implementation, we replicated two other ecological speciation simulation models, namely, the model in (Sadedin et al., 2009) concerning the Swedish L. saxatilis ecotypes and the model in (Gavrilets et al., 2007) concerning the speciation process of cichlids in a crater lake (Barluenga et al., 2006). This replication effort comes in the spirit of exploring the robustness and reproducibility of individual-based models in order to augment the credibility and consistency of computational modeling and simulation and thereby improving the evolutionary ecology theoretical field in general and speciation in particular (Massol and Debarre, 2015; Thiele and Grimm, 2015; David et al., 2017).

Regarding the replicated models, we have implemented only the mating system by similarity (matching rules sensu Kopp et al., 2018) with asymmetric resource availability. Also, regarding the cichlids model, we have implemented only the two patches space structure. The replication of such models allowed us to validate our implementation, while at the same time it permitted us to interpret the sensitivity of the evolutionary outcome to some of the assumptions in the model, as the type of preference function and the presence or absence of mate choice cost.

Analysis

We analyzed the data obtained in the simulations using different variables. For example, the occurrence of ecological adaptation to the exposed habitat (lower shore) was recorded when the ecological trait changed its mean value by more than 75% (x < 0.25) as it is an arbitrary but reasonable indicator that most individuals have adapted to the new habitat and similar thresholds have been used in the past (Sadedin et al., 2009). Similarly, reproductive isolation evolution was recorded when the mean choosiness value was higher than 0.1 or, that is the same, the mean value of c was appreciably different from 0.5 (c < 0.45 or c > 0.55). In addition, we measured the genetic differentiation between habitats by means of the F ST (Wright, 1951) and Q ST indices (Leinonen et al., 2013). The fixation index F ST was computed for the microsatellites as F ST = 1−(((H low + H mid + H up)/3)/Htotal), where H shore was the heterozygosity at the lower, middle or upper shore levels, and H total was the total heterozygosity without considering the shore levels. Similarly, the quantitative genetic analogue of F ST was computed for the ecological trait as Q ST = V between/(V between + 2*V within) where V is the variance between or within habitats. Hence, if high mating discrimination happened in the middle habitat, c ≥ 0.9, jointly with the occurrence of ecological adaptation, Q ST ≥ 0.9, then this was considered as reproductive isolation between the two ecotypes and can be interpreted as a positive case of speciation caused by ecologically-based divergent natural selection i.e., ecological speciation.

The simulation studies should be statistically analyzed in order to be properly understood and gain scientific credibility (Lotterhos et al., 2018). To evaluate the different simulated scenarios in the Galician model, we investigated various simulation outcomes (dependent variables) that could be influenced by a priori different combination of simulation parameters (independent variables) under a classical multifactorial ANOVA (see studied variables in Table 2). First, our independent variables were combined using an orthogonal design (24 parameter combinations ×20 replicates per combination). Notice that we do not use ANOVA as classical way for hypothesis testing, but rather as an exploratory tool to compare the relative contribution of different factors (and interactions) into the dependent ones. The relative contribution of each factor to the overall variation was measured by the eta-squared coefficient (Sokal and Rohlf, 1981). The eta-squared is finally expressed as % with respect to the sum of all etas in the ANOVA (in order to allow comparison among different ANOVAs). In this way, we can identify the most interesting relationships (those with higher eta-squared i.e., higher proportion of variance of the dependent variable explained), which can be detailed a posteriori in the following sections. The per habitat average of the quantitative trait phenotypic values, x and c, were used as dependent variables under this ANOVA approach, and also the within habitat population size (N, considering all demes 1per habitat), as well as the F ST and Q ST estimates described above. Additionally, we used the Pearson correlation (r) of x values among partners in mates within habitat, as a way of estimate the assortative mating for the ecological (magic) trait.

TABLE 2

Independent (factors) and dependent (response) variables for the ANOVA.

Independent variables	Selection strength (σs)
	Mating strength (σa)
	Number of loci (L)
	Selection in the mid-shore (θ)
Dependent variables	Ecological trait (x)
	mating trait (c)
	Mating correlation for x (r)
	Population size (N)
	Differentiation measures (F ST and Q ST)

Results

When developing the Galician model we adapted it to replicate the 2 × 1 model (two habitats, 1 deme per habitat) in (Gavrilets et al., 2007) concerning the speciation process of cichlids in a crater lake (Barluenga et al., 2006). We obtained the same results regarding speciation events, choice trait, and F ST values, as reported from the original article. We also reproduced the similarity scenario for the Swedish L. saxatilis model ( ) using their same mating preference function and obtained qualitatively the same results regarding adaptation, neutral genetic differentiation, ecotype formation, and non-random mating evolution.

Finally, we repeated the simulations with the Swedish model but using the preference function defined in (3) instead of the original one (Sadedin et al., 2009). We obtained some new cases of evolution of negative assortative mating and fewer invasion and adaptation cases. These results are expected since it is known that the original function in (Sadedin et al., 2009) has an anomalous negative assortative mating behavior (Debarre, 2012) implying less negative assortative mating cases than expected (Carvajal-Rodríguez and Rolán-Alvarez, 2014). However, the main conclusions in Sadedin remained, namely, ecotypes may persist indefinitely with moderate genetic differentiation without speciation. As Sadedin et al highlights both ecological divergence and the evolution of reproductive isolation, are crucial to ecological speciation but since they are driven by different forces that can even conflict, the prediction about ecological speciation requires highly detailed data jointly with a quantitative simulation model specially designed for the problem at hand as we did with what we have called the Galician model.

Galician Model

The results are first explored by factorial ANOVA to emphasize the main factors (and interactions) contributing the most (highest variance explained) to the dependent variables, and second the details of these effects are fully described in the next sections.

Parameter Effects: Orthogonal Analysis of Variance

The relative importance of the different simulated parameters, as depicted by the selection and mating strength, number of genes, etc, was summarized by the ANOVA from which the factors and interactions that explained the larger percentage of variance in the evolutionary outcome can be outlined (Table 3, also see Supplementary Material for the full table of the factorial ANOVA).

TABLE 3

ANOVA results showing % of variance (% partial eta-squared; see text) explained by each factor for the ecological trait (x) mate choice (C), mating correlation (r), population size (N), and differentiation measures (F ST and Q ST).

Factor	Response
	x	C	r	N	F ST	Q ST
σs	40*, 35*, 26*	26*, 20, 58*	46*, 64*, 45*	30*, 29*, 35*	48*	38*
σa	3, 1, 9*	1, 0, 1	5*, 9*, 0	1, 0, 1	2*	1
L	0, 9*, 0	9*, 9, 0	8*, 1, 18*	7*, 1*, 21*	7*	18*
θ	4, 3*, 1	12*, 9, 3	2, 0, 0	12*, 34*, 0	0	22*
L × σs	42*, 13*, 31*	16*, 16, 30*	33*, 15*, 30*	12*, 2*, 33*	31*	30*
θ × σs	2, 14*, 0	19*, 17, 0	2, 2, 0	20*, 28*, 2	0	31*

Only the two most important factor interactions are presented. The values within cells correspond to the % at lower, middle, upper shore, by this order, except for differentiation (F ST and Q ST) that are between shore levels. The asterisk indicates significance at the 0.001 level.

The selection strength σs was clearly the key element with highest effect in adaptation, choice evolution, population size, and differentiation. The number of genes had typically a significant although lower impact on most dependent variables. The mating strength (σa) and selection on the mid-shore (θ) showed in general low impact on most variables. On the other hand, the interaction of selection with the number of genes, and with less importance the θ × σs interaction, had also great importance for most variables. Some factors require special attention in certain variables. For example, the number of genes (L) shows a moderate impact in the mating correlation (r) and both F ST and Q ST variables. In the following, we present these outcomes more in detail.

Ecotype Formation: Habitat Colonization and Adaptation

The founder population had the ecological trait perfectly adapted (x = 1) to the sheltered (upper shore) habitat. As time passed some individuals migrated to the intermediate and lower habitats where the ecological conditions were different so that the optimal value for the exposed habitat (lower shore) was x = 0. In the previous section, we have emphasized that σs and L × σs were the most relevant factors determining the evolution of the ecological trait x (Table 3). We will now describe these patterns in detail. In Figure 1, we may appreciate the summary of the adaptation process after 20,000 generations, through the interplay between migration and selection within a habitat, averaged over the different simulated scenarios and demes.

FIGURE 1

Mean value per habitat of the ecological trait for individuals living in the sheltered (upper), intermediate and exposed (lower) habitats.

First, adaptation to the upper shore (Figure 1, white bars, x ≥ 0.75) happened in virtually all cases which are not surprising since the founder individual was perfectly adapted to this habitat.

Habitat colonization and adaptation to the lower shore happened under weak to intermediate selection ( ), irrespective of the value of the other parameters. Under strong selection, (σs = 0.15) adaptation occurred in cases with few loci (L = 4) and non-neutral hybrid zone (optimum 0.5) and also under neutral hybrid zone and mating strength σa = 0.1. For cases corresponding to scenarios under strong selection and more number of loci (L = 8), adaptation only occurred when the hybrid zone was non-neutral. Thus, we may appreciate that more than 80% of the simulated cases successfully colonized the new habitat and adapted to it having the optimal ecological phenotype (Figure 1, black bars, x ≤ 0.15, Figure 1).

The middle shore (Figure 1, grey bar) had some cases where generalists (0.4 < x < 0.6) evolved. The evolution of the generalists in the middle shore corresponded to scenarios with non-neutral hybrid zone and strong selection (σs = 0.15) with adaptation in the lower shore. Other cases under strong selection evolved ecological trait phenotypes close to the upper shore optimum (Figure 1, x > 0.75). In the latter cases, there was no colonization of the lower shore and they only differ from the cases in which generalists evolved in that the ecological trait was neutral in the hybrid zone. Thus, under strong selection in the upper and lower shores, the environmental conditions of the hybrid zone as defined by the selection in the mid-shore (θ) are key to the evolution of generalists and the colonization of the lower shore.

However, more than 60% of the cases in the middle habitat had an ecological trait value closer to the optimum for the lower habitat (Figure 1, x < 0.25). This asymmetry, i.e., ecological trait values in the intermediate habitat being closer to the lower shore than to the upper shore optimum can be explained by the higher carrying capacity of the lower habitat which implied a higher gene flow from lower to the middle than from upper to middle through the vertical axis of the Galician coast.

To confirm the cause of the asymmetric effect for adaptation in the middle habitat, we studied a specific symmetric case corresponding to L = 4 loci under intermediate selection (σs = 0.45) and neutral hybrid zone. We simulate 100 replicates of Galician model and also 100 replicates of a symmetric model with the same number of demes in the upper and lower shore and equal per deme carrying capacity in the three habitats (Supplementary Table S1). We found no differences in the adaptation to the lower shore and confirmed the effect of differential migration to adaptation to the middle shore where the mean ecological trait value was 0.15 ± 0.04 for the Galician model, while it was 0.48 ± 0.07 for the symmetric model (Supplementary Table S2).

The asymmetric effect is a significant difference with the Swedish model where intermediate zones are non-neutral (optimum 0.5) and their distribution is not localized just in the middle point of the vertical coast axis but in the circular horizontal borders between sheltered and exposed so before each sheltered and at the end of each exposed in a circular (as in an island coast) distribution.

In addition to the ecological trait values, we also studied the time to adaptation to the exposed habitat. Basically there were only two opposite outcomes; most scenarios produced fast adaptation in less than 1,000 generations but some did not reach adaptation in 20,000 generations (they appear as the right-most bar in Figure 2). As before, the strength of selection was the fundamental factor affecting the speed of adaptation. Most fast adaptation cases involve intermediate-low selection strength, while most cases of non-adaptation corresponded to strong selection.

FIGURE 2

Generations until adaptation to the lower shore (X < 0.25) for the different scenarios assayed.

There was also an interaction effect between the strength of selection and the number of loci (L × σs in Table 3) implying that under strongest selection (σs = 0.15) there was twice the probability of adaptation under the fewer loci scenario (L = 4) than under the greater one (L = 8).

Mating Trait

We have already emphasized that it was σs and σs × L the most consistent factors determining the evolution of c (Table 3). We will now describe how these effects occur in detail. The phenotypic value for the male mating trait, c, was initially 0.5 which implies random mating (i.e., no choosiness |C| = |2c−1| = 0). Values of c above 0.5 imply positive assortative, while below 0.5 imply negative assortative mating. To discard noisy variation around 0.5, we considered that positive mate choice evolved when the mean c values were above 0.55; while negative choice evolved when mean values were below 0.45. Thus, from Figure 3 we appreciate that positive assortative mate choice (c > 0.55) evolved in the three habitats in most cases although there were some cases (about 15% in the three habitats) with negative mate choice (c < 0.45, Figure 3).

FIGURE 3

Mean value of the mate choice trait for individuals living in the upper (sheltered), middle, and lower (exposed) habitats.

There was a clear pattern behind the negative assortative mating scenarios, the combination of few loci (L = 4) with strong selection (σs = 0.15) (Supplementary Table S3). The negative assortative mating only evolved under strong selection scenarios and the most favorable scenarios include few loci, low tolerance (σa = 0.05), and the selective-middle zone. It seems that under this scenario it was useful to evolve a preference for the different types in order to maintain the polymorphisms avoiding the fast fixation of sub-optimal genotype combinations. If fact, these negative assortative mating cases did not evolve fast, i.e., it took at least 1,500–10,000 generations until the x trait was below 0.25 in the exposed shore and after 20,000 generations the mean fitness still was sub-optimal (w < 0.9) in the three habitats.

When the selection was not so strong and/or the number of loci was high, the result was mostly positive assortative mating (70% of cases in the three habitats with average choosiness of 0.45, c = 0.725, see next section).

As with the ecological trait, we also studied the time required to evolve choice in the lower habitat. Interestingly, the time needed for this character to evolve is much longer than for the ecological trait.

Mating Trait Versus Mating Correlation

The increase of reproductive isolation in speciation scenarios is typically estimated using size (or ecologically related traits) assortative mating (Jiang et al., 2013; Janicke et al., 2019). However, the evolution of assortative mating can be only produced by means of mate choice (choosiness in our model). Therefore, it is of special interest to compare a genetic trait that has evolved, such as male mating trait (c) values, with the observed mating pattern, as measured by the correlation (r) between male and female phenotypes x (ecological trait) in mated pairs, since this is a common measure of assortative mating (Figure 4). This could help to identify potential situations in which demographic effects may give the wrong impression of mate choice evolution.

FIGURE 4

Mating trait value (c, abscissa axis) vs. the mating pattern as measured by the correlation (r, ordinate axis) for the ecological trait in mated pairs.

There was good agreement between the sign of the choice and the mating pattern described by the correlation. Regarding choosiness, higher |C| values implied higher absolute correlation values and vice versa. The agreement was quite good in the range 0.2 < c < 0.8 i.e., a choosiness below 0.6 (|C|< 0.6). With choosiness values above 0.6, (c ≤ 0.2 and c ≥ 0.8 values) there was a saturation of the correlation value reaching its maximum in several cases. Also noteworthy is that there were some cases with high choosiness that presented very low correlation (points close to the abscissa at the right-hand end of Figure 4). These latter cases corresponded to scenarios where the ecological trait was virtually fixed.

In general, about 1/3 of the cases, (33, 34, and 29% for the lower, middle, and upper shore, respectively) reached intermediate (0.4–0.6) positive choosiness values as a result of the trade-off between divergence and sexual selection thus favoring local adaptation under migration (Rettelbach et al., 2013; Cotto and Servedio, 2017; Sachdeva and Barton, 2017).

As with the ecological trait, we studied the possible effect of the asymmetry of the Galician model in the evolution of the mating trait and observed that in this case, the asymmetry in carrying capacity seems to have little or no effect on the evolution of the mating trait (Supplementary Table S2). The values obtained for the mating trait matchwell the empirical estimation of mate choice values for L. saxatilis and other littorinids (Fernández-Meirama et al., 2017a; Fernández-Meirama et al., 2017b).

Neutral (F ST) and Quantitative (Q ST) Genetic Differentiation

Regarding the relative contribution of the distinct factors to the population differentiation, the analysis of variance (Table 3) showed that the selection strength (σs) and its interaction with the number of loci (L × σs) had the highest effect on both neutral and quantitative differentiation. In addition, all the interactions in the table involving strength of selection (σs) as well as the number of loci (L) showed a large impact on quantitative differentiation (Q ST). This suggests the great importance of the strength of selection (as a factor and in interactions) in determining the variance explained in QST. However, the ANOVA did not investigate the relationship between F ST and Q ST.

The Q ST vs. F ST comparison (Figure 5) presented the typical pattern of local adaptation in the presence of gene flow, with high Q ST vs. low F ST values (Leinonen et al., 2013). The only exception occurs when adaptation is not reached within the 20,000 generation interval (e.g., strong selection with L = 8 loci and neutral hybrid zone, see the ecotype formation section), in that cases the Q ST values remained low.

FIGURE 5

Neutral vs. quantitative genetic differentiation.

Intermediate Habitat Effect

We have modeled two different ecological scenarios for the intermediate habitat (2 demes). In the selective-middle habitat scenario, there is positive selection for the ecological trait with the optimum in 0.5. In the neutral-middle habitat scenario, the ecological trait was neutral. The latter is a realistic scenario not usually considered in ecological speciation modeling studies (but see Cotto and Servedio, 2017). The kind of scenario (selective or not) had a minor impact (affecting mean trait values) regarding the ecological trait (x) in the middle shore but a moderate to high impact for c (at the lower shore), N (from low to mid-shore), and Q (Table 3).

Both under the selective and neutral-middle scenarios, there were in the middle shore, a higher percentage of cases with mean ecological phenotype closer to the lower zone optimum (47 or 64% for selective or neutral middle respectively) than to the upper zone optimum (0 or 11% for selective or neutral middle, respectively). This asymmetry was already evident in Figure 1, and was related to the Galician microhabitat configuration jointly with the higher carrying capacity of the lower habitat that provoked a higher number of migrants arriving from this habitat (Supplementary Table S2). This result matches the empirical observation regarding the size of hybrids in the mid-shore which even though they are genetically more heterogeneous than Crab (upper shore) or Wave (lower shore) ecotypes they tend to be closer to the lower shore phenotype (Galindo et al., 2013; Diz et al., 2021).

Regarding the impact of the middle scenario over the colonization of the lower shore, there were more colonization cases of the exposed habitat when there was a selection in the middle shore (97% of those cases) than when the middle shore was neutral (79% of those cases). The analysis of variance (Table 3) showed that the selection strength (σs) and its interaction with the middle scenario (θ × σs) had a significant effect for the x trait in the middle zone.

Regarding the mating trait and mating correlation, the middle habitat scenario did not contribute to explaining any significant variance for the mating trait or the correlation over the simulations (Table 3 and compare rows within Supplementary Tables S4, S5). Also, the choosiness between the different habitats was quite similar (compare columns in Supplementary Table S4). However, there were important differences in the correlation between habitats. The mean mating correlation was 0.24 in the exposed area, 0.55 in the middle, and 0.34 in the sheltered area. However, these general mean values are affected by cases with negative choices and also by cases with poor adaptation. Obviously, when excluding the negative choice scenarios the correlation increases but the pattern with higher correlation in the middle habitat still holds. Specifically for the intermediate selection cases, the correlation values in the middle habitat can be above 0.85 (Supplementary Table S5). The effect was even higher when the symmetric model is considered (Supplementary Table S6) so that the higher correlation in the middle habitat was not an effect of the asymmetry in the Galician model.

Finally, regarding genetic differentiation, F ST was not affected by the presence or absence of selection in the intermediate habitat, contrary to quantitative differentiation (Q ST) which was 20% higher on average when the middle habitat was selective (Q ST = 0.93), compared to when the middle was neutral (Q ST = 0.74). The type of middle habitat scenario explained 22% of the variance in Q ST over the simulations (Table 3).

Ecological Speciation

We defined as a proxy for complete reproductive isolation the concurrence of high local adaptation (Q ST ≥ 0.9, see Figure 5) and high middle-shore choosiness (see Figure 4, c ≥ 0.9 i.e., positive assortative mating with choosiness |C| ≥ 0.8). This reproductive isolation caused by ecologically-based divergent natural selection is called ecological speciation. Thus, when taking jointly the evolution of adaptation and assortative mating, we observed that there was ecological speciation in about 5% of the simulations. All the ecological speciation cases happened under a few selective loci (L = 4) scenarios. However, they were not uniformly distributed for the different selection strengths with all speciation cases linked to intermediate or weak selection strength. We increased the simulation runs to 100 to check this result and observed an 11–15% percentage of speciation for intermediate and weak selection, while there was no speciation for the strong selection cases (σs = 0.15, Supplementary Table S7). These same results were also obtained in the long-term (80,000 generations).

As before, the relative importance of the different evolutionary scenarios was summarized by the ANOVA (Supplementary Table S8). The strength of selection was again the most important factor influencing the dependent variables. The factor L × σs, i.e., the interaction between the number of loci and strength of selection, had decreased importance. In general, many more factors were now non-significant in several variables, which suggests that with cost less adaptation and evolution of variables is obtained.

In the previous scenarios, there was no cost for being choosy. Adding a cost to the mate choice implies that the overly choosy males may remain unmated. As a result, there were no evolutionary outcomes in the form of negative choice or negative assortative mating (negative mating correlation). Also, the less positive choice was evolved in the three habitats (see Table 4).

TABLE 4

Mean choosiness |C| linked to positive assortative mating (negative values excluded) for the different habitats with selective or neutral-middle scenario and with or without mate choice cost.

	No cost			Cost
	Lower	Middle	Upper	Lower	Middle	Upper
Selective middle	0.42	0.42	0.44	0.34	0.34	0.34
Neutral middle	0.49	0.47	0.47	0.39	0.38	0.39

Values are averages for the different L and σa.

Finally, there were fewer colonization events in the presence of cost. Under the selective-middle habitat scenario, 83% of the simulations with cost underwent colonization, compared to 97% without cost. The same pattern occurred under the neutral-middle habitat scenario, 66% of colonization with cost vs. 79% without cost. However, when colonization occurred the fitness in the exposed habitat was higher under cost scenarios (mean fitness 99.7%) indicating that adaptation was better in these cases. Finally, we do not find ecological speciation when mate choice cost was included but in just one run of the case with weak selection (σs = 1), low number of loci (L = 4), higher tolerance (σa = 0.1), and selective-middle scenario. This case implies a 0.2% compared to 5% speciation cases when the cost was absent.

Discussion

We have analyzed the evolution of local adaptation and mate choice in a system displaying extreme microhabitat environmental heterogeneity. This system resembles the L. saxatilis microhabitat-associated dimorphism along the wave-exposed rocky shores of Galicia. L. saxatilis dimorphism is a good example of microgeographic adaptation and incipient-speciation which occurs despite high gene flow based on the species expected levels of dispersal (Rolán-Alvarez 2007; Richardson et al., 2014). In our individual-based simulations, we investigated the dynamics of the ecotype formation through the inter-deme and habitat spatial and temporal scales. This deme-based fine spatial scale model showed how the strength of selection and the evolution of choosiness may interact in order to produce local adaptation and reproductive isolation (Cotto and Servedio, 2017).

Our results are consistent with the key role of natural selection in the ecological speciation processes (Barton, 2010). The strength of selection on its own or through interaction with other factors explains most of the variance in the dependent variables over the scenarios simulated for both with and without mate choice cost. That is, most of the differences in ecological and mating trait values are explained by the strength of selection.

We have seen how the population, as the individuals disperse through inter-habitat demes, evolve different phenotypic values for habitat specialization. Given the high environmental heterogeneity between the sheltered and exposed habitats, we obtained specialists for these habitats. The outcome was slightly different in the middle habitat, where the generalist or specialist evolved depending on how the environmental heterogeneity was defined. When the middle habitat was neutral this facilitated the co-existence of specialists from both shore levels although with a higher % from the lower shore due to the higher number of migrants coming from this shore. On the contrary, when the middle habitat was non-neutral with an optimum in-between both extreme shore levels, which implies a less-heterogeneous environment, more intermediate phenotypes evolved. These results are consistent with previous studies showing the positive relationship between local adaptation and environmental heterogeneity (Berdahl et al., 2015; Svardal et al., 2015).

Concerning our question number 1 about the pattern of local adaptation under gene flow, which expects high quantitative genetic divergence for certain traits, without a general divergence at the genome level, we confirmed the expected pattern and observed high quantitative trait genetic divergence (Q ST) jointly with low neutral genetic differentiation (F ST).

We found that the time-scale for ecological adaptation and the occurrence of ecological speciation were influenced by the strength of selection. If the selection was strong the time for adaptation was longer. Moreover, weak or intermediate selection and few selective loci were favored under ecological speciation. However it should be noted that we modeled independent loci, i.e., a few unlinked loci may represent a relatively wide genomic region depending on the linkage relationships of the species. Obviously, the fewer the loci, the smaller the genomic region. Our results match the previous one for the Swedish model (Sadedin et al., 2009) where ecotype formation was slower when the number of loci involved was larger, and also recent results for the L. saxatilis Swedish system that shows the importance of adaptive recombination suppression within inversions to facilitate ecotype evolution and adaptation (Koch et al., 2021). Moreover, recent work (Kautt et al., 2020) has highlighted that simple genetic architecture, like the one we modeled, does not necessarily lead to speciation even with magic traits which are also consistent with our findings. It seems, however, that polygenic selection might be more efficient in driving sympatric speciation under some circumstances however this setting was out of the scope of the present study.

Most of the simulations evolved intermediate values of choosiness as a result of the trade-off between the effect of natural and sexual selection. Similar results have been obtained previously under some theoretical models and specific parameter range (see e.g., Rettelbach et al., 2013; Cotto and Servedio, 2017; Sachdeva and Barton, 2017 and references therein) and it is interesting that we have obtained the same result with a model for the specific case of L. saxatilis in the Galician coast using empirical estimates of the parameters (when available) and that the intermediate choosiness obtained in the simulations matches quite well empirical estimates of L. saxatilis and other littorinids (Fernández-Meirama et al., 2017a; Fernández-Meirama et al., 2017b; López-Cortegano et al., 2020).

Therefore, regarding our question number 2 on mating choice evolution as an effect of ecological speciation we may say that yes, under L. saxatilis Galician-like conditions, choosiness evolved during the process of ecological adaptation. Most of the evolved choosiness corresponded to positive assortative mating although negative assortative mating evolved under the combination of few loci and strong selection. As a general pattern, intermediate selection and few selective loci favored the evolution of choosiness. However, this does not mean that the evolution of choosiness will irremediably cause ecological speciation. Actually, choosiness was typically maintained at intermediate values.

This was caused because there was an interplay between the evolution of assortative mating and local selective pressures. We modeled the evolution of choosiness as an ecologically neutral trait, although having as the target trait a magic trait, with realistic migration estimated for L. saxatilis in the Galician coast (Rolán-Alvarez et al., 2015a). These conditions imply that the evolution of choosiness is constrained by the maintenance of the trait polymorphism required for evolving local adaptation which is in turn affected by the current choosiness. If choosiness is strong, then rare, non-optimal, phenotypes would mate between themselves. Therefore, during the process of local adaptation, the alleles coding for strong choosiness are eliminated by viability selection. Under this setting, the highest amount of divergence occurs at intermediate choosiness values for a wide range of parameter values, which has been also observed in other studies (Rettelbach et al., 2013; Cotto and Servedio, 2017; Sachdeva and Barton, 2017). As far as we know, this is a new outcome regarding the evolution of mate choice for the L. saxatilis systems and is consistent with our intermediate estimates of assortative mating ( Johannesson et al., 1995; Rolán-Alvarez et al., 1999; Rolán-Alvarez et al., 2004; Cruz et al., 2004 ) and mate choice (Fernández-Meirama et al., 2017a; Fernández-Meirama et al., 2017b) for the Galician system. Positive assortative mating evolved in the three simulated habitats.

We also studied the correspondence between the preference based on the c trait and the assortative mating correlation measure. We already know that both measures matched well if the sampling is done at the correct scale where choice occurs (Estévez et al., 2018). However, we found important differences when comparing the choosiness values (recall that choosiness is a linear function of the c trait) with the value of the assortative mating correlation measure, especially when interpreted in the context of reinforcement. Reinforcement is defined in a broad sense as the evolution of traits that minimize hybrid formation in response to selection (Servedio et al., 2013; Pfennig, 2016). The effect of selection can be indirect i.e., reinforcement occurs due to selection against hybrids that indirectly results in prezygotic isolation, or on the contrary, it involves the effect of natural or sexual selection acting directly on the prezygotic phenotypes involved in mate attraction. In the latter case, reinforcement is better called reproductive interference (Butlin and Ritchie, 2013; Hochkirch, 2013; Shaw and Mendelson, 2013).

We found a higher mating correlation in the middle habitat than in the lower and upper ones. Therefore, if we use mating correlation as a proxy for assortative mating (Jiang et al., 2013; Janicke et al., 2019) we would say that we are detecting reproductive interference.

However, choosiness in the middle habitat, where hybrids can be formed, was very similar to choosiness in the other habitats, and so actually we cannot conclude that reproductive interference was evolved in our scenarios.

This putative reinforcement-like pattern cannot be due to a character displacement for the ecological trait (x) since there were no more extreme phenotypes in the middle area than in the other two. However, there was higher variation for the ecological trait in the middle habitat than in the others (coefficient of variation was at least one order of magnitude higher). Therefore, we should infer a higher sensitivity of the correlation measure to outliers (Pernet et al., 2013) which on the other hand would also explain the observed saturation effect of r over c. Actually, estimates from L. saxatilis empirical data indicate that mate choice has not increased in populations where the two ecotypes meet compared with those with only one ecotype (Fernández-Meirama et al., 2017a). All of the above corroborates that evolving high mate choice levels is difficult in presence of high gene flow. Therefore, regarding question number 3 in the introduction, we may say that mate choice does not increase in the middle habitat under the simulated conditions of the Galician L. saxatilis system.

Finally, we considered question number 4 about the impact of adding a cost to the mate choice and we observed in concordance with previous studies (Kopp and Hermisson, 2008), fewer reproductive isolation cases and just residual ecological speciation when the cost was added. As expected, when choosiness evolved its value was lower, but still relevant to have an ecological impact. In addition, there was increasing importance on how the middle habitat was modeled and its impact on adaptation. When the middle habitat was neutral regarding the ecological trait, there were twice more adaptation failures than when the middle habitat was selective because when selection favors an intermediate phenotype this probably facilitated the transition to the optimal phenotypes for the lower shore (Schneider and Burger, 2006).

In summary, it seems that the L. saxatilis model system may correspond well to an incomplete ecological speciation case with intermediate choosiness and strong within habitat adaptation. The results obtained in our model system suggest that the evolution of ecological speciation is not an immediate consequence of local divergent selection and mating preferences, but a fine-tuning among the environmental conditions in the microhabitat and the hybrid middle zone, the genetic basis of the traits, the selection intensity, and the mate choice stringency and cost. Similar results which have been obtained for other ecological speciation model cases (Safran et al., 2013; Cotto and Servedio, 2017) suggest that incomplete ecological speciation may not be rare.

These results may guide the search for new empirical data providing direction for future models to assist our understanding of the L. saxatilis system and ecological speciation in general. For example, it would be interesting to check the particular condition that could facilitate higher values of choosiness and so putatively complete the ecological speciation process in this model system or even test whether the mate choice mechanism functions as a similarity-like mechanism as has been shown in other littorinids (Lau et al., 2021).