Assortative mating, sexual selection and their consequences for gene flow in Littorina

When divergent populations are connected by gene flow, the establishment of complete reproductive isolation usually requires the joint action of multiple barrier effects. One example where multiple barrier effects are coupled consists of a single trait that is under divergent natural selection and also mediates assortative mating. Such multiple-effect traits can strongly reduce gene flow. However, there are few cases where patterns of assortative mating have been described quantitatively and their impact on gene flow has been determined. Two ecotypes of the coastal marine snail, Littorina saxatilis, occur in North Atlantic rocky-shore habitats dominated by either crab predation or wave action. There is evidence for divergent natural selection acting on size, and size-assortative mating has previously been documented. Here, we analyze the mating pattern in L. saxatilis with respect to size in intensively-sampled transects across boundaries between the habitats. We show that the mating pattern is mostly conserved between ecotypes and that it generates both assortment and directional sexual selection for small male size. Using simulations, we show that the mating pattern can contribute to reproductive isolation between ecotypes but the barrier to gene flow is likely strengthened more by sexual selection than by assortment.

term 'magic trait' is used for a subset of multiple-effect traits where the trait under 48 divergent selection also contributes to assortative mating (Servedio et al. 2011). Trimble total station (as in Westram et al. 2018). "Reference" snails (used as mating 199 partners, see below) were sampled at a fifth island ("ANG" in Westram et al. 2018; N 200 58°52'15.14", E 11°7'11.88") in Crab and Wave habitats away from the contact zone (in 201 total 2 0 0 individuals from each habitat per test shore). Both reference and transect 202 snails were sexed prior to mating experiments based on observation of the male penis. 203 If no penis was observed, individuals were assumed to be females. If the penis was 204 underdeveloped, individuals were considered sexually immature and excluded from the 205 mating experiments. Dissections of transect snails followed all experiments in order to 206 confirm initial sex determination and check whether females were immature or 207 parasitised. Trials involving immature or parasitised transect individuals, or individuals 208 whose sex had been determined incorrectly, were discarded. 209 Size was measured for both reference and transect snails as the maximum distance 210 between the top of the apex and the base of the aperture of the shell. Shape was 211 determined only for the transect snails and summarized as the first relative warp from a 212 landmark-based geometric morphometrics analysis, which captures the Crab-Wave axis 213 of variation (Ravinet et al. 2016;Westram et al. 2018). Shell shape of the reference 214 9 only whether or not a mating occurred in each trial. A positive outcome was recorded if 245 the male was in the mounting position for more than 1 minute (Saur 1990). 246

DATA ANALYSIS 248
For each mating pair, we had information about whether a mating event was observed 249 or not, the island where the transect snail was collected (CZA, CZB, CZC or CZD), 250 transect snail shape, the ecotype of the reference snail (Crab, Wave), the sex of the 251 transect snail (and therefore of the reference snail) and the sizes of the two snails, 252 which were used to calculate the ratio between the female and male size for each 253 included size ratio effects, with the square of size ratio being the strongest effect, but 273 varied in the other factors that entered the model (Table S1, S2). 274 We then fitted a model to the observed data in order to describe the relationship 275 between mating probability and the ratio of female to male size. We selected a function 276 to account for the decline in probability away from an optimum size ratio (the effect of 277 the square of size ratio in the GLM). This model allowed us to estimate parameters for 278 the mating pattern that we then applied to size distributions in nature to infer the 279 assortative mating and sexual selection generated by the mating pattern. The 280 parameter estimates were also used to simulate the barrier effects of size-assortative 281 mating and sexual selection (see below, CLINE SIMULATIONS). Initial trials showed 282 that a symmetrical Gaussian model that is commonly used to describe sexual selection 283 and assortative mating (Lande 1981;Gavrilets 2004) could not account for our 284 observations because the mating probability declines asymmetrically, more rapidly for 285 males larger than females than for males smaller than females. Therefore, the binary 286 outcome of the mating experiment (mated or non-mated pair) was fitted using logistic 287 regression to a skew normal function of the size ratio. Specifically, we expressed the OR) also depends on the parameter ߙ , which controls the amount of skew (Fig. 1). The 296 OR was estimated by taking the first derivative of Eq. (1) using Wolfram|Alpha (access 297 October 19, 2018) and finding its root using the function uniroot()in R version Our initial data exploration using generalized linear models (see above) suggested that 326 the relationship between mating probability and size ratio might vary according to island 327 and ecotype (or snail shape). Furthermore, although the unit of replication was the 328 transect-reference pair, it remained to be tested whether individuals differences in 329 shape between transect snails could have explained part of the variation in mating 330 probability. We tested the impact of these variables by fitting hierarchical models in 331 Stan. In these models, we replaced one or more of the parameters in Eq. (1) by a 332 'hyperparameter' that was a function of the island from which the transect snail was 333 sampled, the transect snail shape, the reference snail ecotype and the sex of the

MATING PATTERN CONSEQUENCES IN THE CONTACT ZONE 338
The parameters of the mating pattern were estimated from the observations in the 339 mating experiment, which was designed to investigate the probabilities of mating given  (Table S5). Generations were discrete and non-overlapping. The lifecycle was modelled in the 405 order: dispersal, recombination and mating, locally in each patch, then natural selection. 406 In the model, dispersal distance was Gaussian distributed with mean zero and standard 407 , so that, due to diploidy, overshooting of the local trait optimum was possible. Mating 419 was implemented according to five different models, one being random mating, and the 420 remaining four being different versions of the mating pattern based on the trait that was 421 also under natural selection (see below). In each model, we assumed that every female 422 produced a large (and the same) number of offspring (i.e., where one or both snails were inactive throughout the trial, 7 0 transect snails without 506 spatial information and 1 1 8 mating pairs with missing shell sizes. 507 508

MATING PATTERN IN THE LABORATORY TRIALS 509
Size-assortative mating acts as a barrier to gene flow when the probability of mating 510 between two populations of different sizes is reduced. To investigate this barrier effect, 511 the first step is to quantify how the probability of mating varies with respect to female 512 and male size distributions. The mating model, Eq. (1), was built for this objective and it 513 was fitted to the data from all four islands combined (Fig. 2). The probability of mating 514 followed a right-skewed distribution with maximum displaced from the center of the 515 distribution towards pairs where the female was 1 . 3 1 times larger than the male and 516 falling rapidly for pairs with other size ratios (Table 1; Fig. 2). As the size ratio between 517 the sexes increased/decreased, the mating function approached a probability close to 518 zero within the range of observed size ratios for males larger than females but not for 519 males smaller than females (

ASSORTATIVE MATING AND SEXUAL SELECTION 553
Clines in male and female size were observed on all four islands with centers close to 554 habitat boundaries (Fig. 3; Fig. S4). In all cases, sexual size dimorphism was greater in 555 Wave snails than in Crab snails, the variance in l n ሺ s i z e ሻ was also greater in Wave snails 556 and the variance increased in the centers of the clines (Fig. 3; Fig. S4; Table S5). 557 After generating virtual mating encounters using a custom script, we computed, for each where the size variance was smallest ( Fig. 3; Fig. S4). 566 Sexual selection was predicted to favor smaller males, and lower variance in male size 567 in all cases ( Fig. 3; Fig. S4). However, sexual selection was also predicted to vary along 568 the transects of the four islands in line with the size variance and difference between 569 female and male sizes. In some cases, the predicted effects were very small.

BARRIER TO GENE FLOW 594
In all five models we simulated (see illustrations in Fig. 4a, e, I, m, q) we found that, at 595 the end of the simulations, the average phenotype of females at the two habitat ends 596 matched their corresponding optimal phenotypes (Fig. 4b, f, j