Chromosomal inversions have been shown to play a major role in local adaptation by suppressing recombination between alternative arrangements and maintaining beneficial allele combinations. However, so far, their importance relative to the remaining genome remains largely unknown. Understanding the genetic architecture of adaptation requires better estimates of how loci of different effect sizes contribute to phenotypic variation. Here, we used three Swedish islands where the marine snail Littorina saxatilis has repeatedly evolved into two distinct ecotypes along a habitat transition. We estimated the contribution of inversion polymorphisms to phenotypic divergence while controlling for polygenic effects in the remaining genome using a quantitative genetics framework. We confirmed the importance of inversions but showed that contributions of loci outside inversions are of similar magnitude, with variable proportions dependent on the trait and the population. Some inversions showed consistent effects across all sites, whereas others exhibited site-specific effects, indicating that the genomic basis for replicated phenotypic divergence is only partly shared. The contributions of sexual dimorphism as well as environmental factors to phenotypic variation were significant but minor compared to inversions and polygenic background. Overall, this integrated approach provides insight into the multiple mechanisms contributing to parallel phenotypic divergence. 

Snails were collected along the Swedish West coast, where Crab and Wave habitats frequently come into contact (Figure 1B), in 2013 and 2014. The sites were Ramsö (58°49’27.8”N11°03’45.3”E), here CZA, Inre Arsklovet (58°50’00.5”N 11°08’19.6”E), here CZB, and Yttre Arsklovet (58°49’51.3”N11°07’59.0”E), here CZD. At each location snails were collected along a transect that started in the wave habitat (exposed bedrock) continued to the crab habitat (boulder fields in the centre of a bay) and extended into the next wave habitat at the other side of the bay. Exact positions of all individuals were recorded in three-dimensional space using a Total station (Trimble M3). Positions along the transects were then converted to a one-dimensional distance measurement along a least-cost path that minimised path length and constrained movement to areas of high snail density. The first snail collected along the transect path (starting in the Wave habitat) was assigned position 0. Total path lengths were: 362.47 m (CZA), 257.00 m (CZB), 270.17 m (CZD). Environmental variables were recorded to describe the habitat at 1000-2000 points for each transect. These included substrate type (bedrock/boulder), presence of barnacles as indicator of wave exposure and presence of fucoid seaweed (indicator of weak wave action and more sheltered habitats). These measurements were then summarised as “habitat PC” by a principal component analysis. Sex of each snail was determined by dissection. Weight (wet weight of snails including shell), shell length and thickness were recorded. Thickness was measured with a thickness gauge (NeoteckDTI Digital Dial Indicator Probe, 0.001mm resolution) at the widest point of the aperture and the average of three measurements per individual was used. Foot area was measured from photos of snails moving inside seawater-filled Falcon tubes using the program ImageJ (http://imagej.nih.gov/ij/). Foot area (in cm2) was divided by shell length to obtain relative foot area. For describing shell shape and aperture characteristics, several size-independent parameters (height-growth, width-growth, aperture-position, aperture-size, aperture-shape) were used, based on a growth model developed by Larsson et al. (Larsson et al. 2020; J. R. Soc. Interface 17.). They describe whether the shell shape is more elongated with a high spire (small height- and width-growth) or globose (large growth parameters). Small values for aperture-position and large values for aperture-shape indicate narrow apertures, wide apertures have a large aperture-position and small values for aperture-shape. Boldness behaviour was measured as Bold Score, which is the log-time until an individual crawls out of its shell after disturbance (higher Bold Score means individuals are less bold). Each individual was measured three times and the average was used as a boldness proxy. For inversion genotypes a PCA on all SNPs within previously described inversion regions (see Faria et al. 2019, Mol. Ecol. 28:1375–1393), was applied which resulted in distinct clusters (see Faria et al. 2019 for details of this method). Inversion genotypes are coded as 0 and 2 for homozygotes and 1 for heterozygous.