Together with genetic predispositions, environmental constraints are the primary causes for shape variability within and among populations. Phenotypic plasticity is the ability of a single genotype to produce a range of different forms in response to its environment (Bradshaw, 1965; West-Eberhard, 1989; Sultan and Stearns, 2005). This trait, if adaptive, is generally perceived advantageous for an organism or population because it facilitates resistance to a range of different stressors allowing similar genotypes to populate heterogeneous habitats (Chevin and Hoffmann, 2017). Therefore, changes in morphology can be strong indicators for increasing or changing environmental stressors, knowledge about which may be important when monitoring the consequences of anthropogenically induced global climate change (Telesca et al., 2018). So far, studies have focused predominately on population and community composition responses (Ashton et al., 2017) or genetic markers to study the effect of global change. However, an understanding of phenotypic consequences was also found to be crucial to forming a complete picture of biological responses to a world in change (Hellberg et al., 2001; Cross et al., 2018; Telesca et al., 2018; Morley et al., 2019; Telesca et al., 2021).

Climate change associated with rising atmospheric carbon dioxide levels, increasing terrestrial and seawater temperatures, ocean acidification (OA) and more frequent extreme weather events is an emerging global threat to most of the Earth’s biosphere (Collins et al., 2013). In the oceans, annually increasing sea surface temperatures and OA pose strong stressors for marine organisms (Byrne, 2011; Fitzer et al., 2015). Species with calcium carbonate (CaCO ₃) shells or skeletons are expected to be especially impacted by global climate change (Waldbusser et al., 2013; Waldbusser et al., 2014; Przeslawski et al., 2015; Suckling et al., 2015; Thomsen et al., 2015). The reduction in seawater pH associated with rising atmospheric CO ₂ concentration results in potentially harmful effects (e.g., shell dissolution), making an organism more susceptible to predation, parasites, or infectious diseases (Gaylord et al., 2011; Viotti et al., 2019; Barclay et al., 2020). Marine calcifiers are reported to be especially vulnerable at the larval stage to the synergistic effects of OA and ocean warming (Przeslawski et al., 2015), as this can lead to premature death or severe retardation of shell formation (Gazeau et al., 2011; Thomsen et al., 2015; Waldbusser et al., 2015). To overcome these impacts, CaCO ₃ bearing organisms are required to spend a larger proportion of their metabolic energy on calcification to maintain the status quo (Melzner et al., 2011), which could have far-reaching implications for their energy balance and competitiveness.

There is a consensus that calcification is linked to CaCO ₃ substrate abundance and speciation, i.e., to the CaCO ₃ saturation state Ω _{CaCO 3} in seawater, although its exact role in the acidified ocean is still debated (Bach, 2015; Cyronak et al., 2015; Waldbusser et al., 2016). The Ω _{CaCO 3} is a function of both calcium ion [Ca ²⁺] and bicarbonate [ C O 3 2 − ] availability as well as the thermodynamic solubility product K_sp of the respective CaCO ₃ polymorphs (for molluscs, most important are aragonite and calcite) at a specific pressure, salinity and temperature (Zeebe and Wolf-Gladrow, 2001) (eq. 1).

In seawater, both [Ca ²⁺] and the dissolved inorganic carbon concentrations (DIC) are tightly coupled to seawater salinity, whereby the speciation of DIC in seawater is pH-dependent. This means that both a reduction in the salinity and the pH can reduce the saturation state of seawater. At Ω _{CaCO 3} < 1, CaCO ₃ dissolution becomes thermodynamically favourable; however, as pointed out by Waldbusser et al. (2016), low Ω _{CaCO 3} states likely pose more of a kinetic rather than a thermodynamic constraint on marine calcifiers considering the rapid rate of calcification observed especially at the larval stage (Waldbusser et al., 2016). Irrespective of whether thermodynamic or kinetic bottlenecks are responsible, continuously low Ω _{CaCO 3} states can have serious implications for shell-forming organisms (Doney et al., 2009; Watson et al., 2012; Waldbusser et al., 2014). In recent years, many in-situ and laboratory studies have added to our understanding of the effect of climate change on marine calcifiers (Byrne, 2011; Byrne and Przeslawski, 2013; Kroeker et al., 2013; Parker et al., 2013; Fitzer et al., 2015; Kroeker et al., 2016; Cross et al., 2018). However, relatively short-term experiments (Hofmann et al., 2010), heterogeneous study designs, and partly conflicting results among studies have complicated interpretations and highlight the complexity of possible organism responses to a changing environment. While long-term, multi-generational experiments are key to studying phenotypic and genetic adaptions to environmental change (Hofmann et al., 2010), they require many months to years and thus pose a serious time constraint (Clark et al., 2019a; Clark et al., 2019b). This is especially true if organisms with long lifespans and/or a slow reproductive cycle are to be studied. Thus different approaches to investigating the effect of climate change in a natural setting that also account for both phenotypic and genetic adaptions receive increasingly more attention (Cross et al., 2018; Telesca et al., 2019; Bullard et al., 2021; Telesca et al., 2021). These types of studies focus on organism responses in the temporal domain by comparing archival collections with modern specimens collected from a strongly constrained spatial range (Fisher et al., 2009; Cross et al., 2018; Telesca et al., 2021) or in the spatial domain by studying the effect of environmental constraints on an organism over an extensive geospatial range, made possible by recent advances in ecological modelling (Telesca et al., 2018; Telesca et al., 2019).

The dog whelk, Nucella lapillus (Linnaeus, 1758), is an important predatory gastropod inhabiting the rocky intertidal zone of the North Atlantic (Burrows and Hughes, 1989; Hughes and Burrows, 1994). Nucella lapillus preys predominately on habitat-forming foundation species such as blue mussels or barnacles and thereby exerts a strong top-down control on rocky-shore ecosystems (Trussell et al., 2003; Benedetti-Cecchi and Trussell, 2014; Contolini et al., 2020). Its spatial distribution ranges from the northern tip of Norway to southern Portugal, spanning three climate zones. Nucella lapillus is a direct-developing gastropod without a larval stage (Moore, 1938). Its locomotor activity is limited to movements between prey (usually no more than a few centimetres) on which individuals are reported to feed for hours, if not days (Fretter and Graham, 1962). Nucella lapillus is only occasionally found below the tide marks and does not voluntarily crawl over sandy or muddy areas. The absence of a larval stage and the limited locomotor range of N. lapillus limits its dispersal so that distinct phenotypes exist within different allopatric populations, forming numerous discrete gamodemes (Berry and Crothers, 1974). Its calcareous shell consists of two layers, made of irregular calcite on the outside and crossed lamellar aragonite on the inside, which may be separated by a transitional spherulitic layer (Mayk, 2021). Site-specific exposure to wave action and predators are the most studied natural causes of shape variability among N. lapillus populations (Hughes and Elner, 1979; Vermeij and Currey, 1980; Palmer, 1990; Pascoal et al., 2012). On exposed shores, N. lapillus usually have short squat shells with a larger aperture size to accommodate a larger foot, which improves adherence to substrata. Squatter shells and better substratum adherence reduce the risk of dislodgment and are thus interpreted to be a morphological adaption to increase survival (Kitching et al., 1966; Guerra-Varela et al., 2009). Nucella lapillus on protected shores, where individuals are more regularly exposed to predators, tend to produce stronger and thicker shells (Currey and Hughes, 1982; Appleton and Palmer, 1988) and have increased growth rates (Burrows and Hughes, 1989). Individuals from protected shores also tend to exhibit relatively lighter soft tissue body weight as thicker shells provide comparably less space to accommodate the internal organs (Kitching et al., 1966). Site-specific exposure regimes have also been associated with genetic variability (karyotypes). Nucella lapillus show a pronounced chromosomal polymorphism; first reported in populations from the coast of Roscoff, France (Staiger, 1953). Known chromosome numbers range from 2n = 26 to 2n = 36 between populations (Bantock and Cockayne, 1975; Pascoe and Dixon, 1994). It was suggested that karyotypes correlate with local conditions, as the 2n = 26 karyotype appeared to be more common on exposed shores, whereas the 2n = 36 karyotype appeared to be more common on protected shores (Kirby et al., 1994). However, other studies did not confirm this finding, suggesting that shell shape in N. lapillus is a direct response to its physical and biological environment, expressed within a single generation and without the need for a change in karyotype (Pascoe, 2002; Pascoe, 2006). In any case, N. lapillus shows pronounced phenotypic adaptions to its habitat allowing this intertidal muricid to become a dominant top-down predator. New and continuous stressors, as brought about by OA and ocean warming, could, however, endanger its position, potentially jeopardising the trophic system on large parts of the North Atlantic rocky shore. While a wealth of studies exists in the literature that examines local effects on N. lapillus’ shell morphology, comparatively little is known about large-scale shape trends, despite its wide distribution. However, this would be of particular interest to understand how this important intertidal predatory gastropod adapts to different environments and has the potential to provide new insights and understanding of this species’ capacity to respond to a world in change. This paper aims to close this gap by investigating response patterns of N. lapillus across its entire latitudinal range with a particular focus on identifying significant relationships between environmental gradients and phenotypic responses which may indicate sensitivities of this intertidal muricid to environmental change. To do so, we analysed shape variance among N. lapillus populations, sampled over a 30 degrees latitudinal range, from the northern tip of Norway to the south of Portugal, and used multivariate ecological models to identify significant relationships of shell shape to global change relevant environmental parameters.

Mean plots of environmental parameters revealed clear gradients from north to south ( Figure 1 ). As expected, mean sea surface temperature increased gradually towards the south. Maximum difference in average temperature between sampling sites was 11.49°C. Seawater salinity showed an increasing trend towards the south with pronounced local variability ( Figure 1 , sampling sites 7-9). Salinity variations between sampling sites were up to 3.14 PSU. Seawater pH decreased between 70°N to 50°N, but sampling sites south of 50°N exhibited higher seawater pH ( Figure 1 ). The maximum seawater pH range encompassed between our sampling sites was 0.13. Minimum seawater Ω _Ca and Ω _Ar values at the sampling sites were 3 and 2 respectively. Seawater Ω _Ca states showed no apparent trend over the whole studied latitudinal range but a significant increase in saturation south of 50°N.

Relative contributions of the first five shape-PCs to the total morphospace range are shown in Figure 2B . Overall, shape variance explained by shape-PC1 (53.17%) was a widening of the shell. Higher shape-PC1 values described wider, more globular shells with almost rectangular shoulders and proportionally larger body whorls. Lower values described a narrowing of the shell with proportionally smaller body whorls and flatter shoulders. Shape-PC2 (14.76%) described a rounding of the body whorl. Higher shape-PC2 values described an overall rounding of the body whorl with a less contoured siphonal canal. Shape-PC3 (8.65%) described an “offsetting” of the body whorl to the spire, which gave the impression of tilting the shell’s longitudinal axis. Shell shape-PC4 (5.94%) and shape-PC5 (4%) only explained minor shape variations. Shape-PC4 described a slight widening and shift of the body whorl, especially around the shoulder, and shape-PC5 described contouring of the shell. A plot of among-population variation described by shape-PC1 and shape-PC2 showed a clear separation among populations, especially along the shape-PC1 axis ( Figure 3 ). Mean plots of selected shell descriptors (shape-PC1 to shape-PC5) also showed a clear effect of locality and significant latitudinal trends. Linear regression analyses confirmed significant shell shape changes in relation to sampling site latitude as expressed by shape-PC1, shape-PC2 and shape-PC3. Shape-PC4 and shape-PC5 did change significantly with latitude ( Figure 2 ).

The best-supported model ( Table 2 ) revealed a significant relationship between the selected shape descriptors and environmental gradients. The model comparison revealed that sea surface temperature, A_T , SIR, shell height, and sampling site ID (as a random effect term) best described the observed changes in shape-PCs. The model identified highly significant non-linear relationships of shape-PC1 with sea surface temperature, A_T , SIR and shell height ( Figure 4A ). Shape-PC2 showed a non-significant relationship with sea surface temperature but highly significant relationships with A_T and SIR. Shell height had only a minor effect on shape-PC2 ( Figure 4B ). Shell shape variance described by shape-PC3 was also not associated with sea surface temperature, but showed significant correlations with A_T , SIR and shell height ( Figure 4C ). Mean effect size estimates ( Figure 5 ) revealed differences in the relative contribution of environmental predictors on observed changes in shape. The SIR and shell height showed a significant effect on shell shape-PC1 to shape-PC3, and in addition, A_T exhibited a significant effect size for shape-PC2. Overall, effect size estimates primarily linked changes in the carbonate system as represented by A_T and especially the SIR to changes in shell shape. Most notably, we identified increasing SIR values resulted in the formation of wider, more globular shells with a more pronounced shoulder. Effect size estimates for sea surface temperature were non-significant for any of the three shape descriptors.

Using a combined EFA and GAMM approach, we identified significant shape changes in the intertidal N. lapillus that followed latitudinal ( Figure 2 ) and environmental trends ( Figures 4 , 5 ). We identified a significant narrowing of shells (i.e., a reduction in the shell aspect ratio) with relatively smaller body whorls from north to south. While we did not identify microclimate effects of individual sites on shell formation, which likely contributed to the observed considerable shape variation between individuals from the same population ( Figure 2 ), our large scale trends are so clear that the variability of mean environmental conditions among sampling sites must have been significant enough to influence the shell formation in the studied specimens. Although the intertidal zone is known to exhibit pronounced temporal and spatial variability in the carbonate system (Krause-Jensen et al., 2015), our results suggest that a shift in the mean has a significant influence on shell formation. Effect size estimates of environmental predictors on shape-PCs revealed that differences in the carbonate system had the strongest influence of all environmental predictors on shell shape. In contrast, sea surface temperature showed a non-significant effect size for any of the three shape-PCs ( Figure 5 ). Despite overall high Ω _Ca and Ω _Ar throughout the sampling sites, shell formation of N. lapillus appears to be remarkably sensitive to variations in the carbonate system. Our models indicate a significant link between shell formation of N. lapillus to variations in the “reactant” [ H C O 3 − ] and “inhibitor” [H ⁺] ratio (SIR) (Jokiel, 2011a; Jokiel, 2011b; Jokiel, 2013), which is expected to decrease in the coming years and centuries under increasing OA conditions. This suggests an overall shift in N. lapillus shell shape can be expected as a response or consequence to changing water chemistry constraints.

There is only a handful of studies explicitly investigating shell shape responses of marine gastropods to changes in the carbonate system, and these have produced partly conflicting results. A study of the shell aspect ratio in the gastropod Phorcus saucciatus collected from a CO ₂ vent and a control site reported a similar shape response to those here, as specimens from the vent site tend to form narrower and more elongated shells (Viotti et al., 2019). However, seemingly conflicting results were reported from laboratory studies where the combined effect of increasing seawater temperature and decreasing pH resulted in the formation of more globular shells in the intertidal gastropod Littorina littorea and in juvenile N. lapillus (Melatunan et al., 2013; Rühl et al., 2017). This discrepancy among studies is not easily resolved but highlights the complexity of comparing species responses from laboratory and field studies. One of the biggest issues in comparing gastropod shape studies is the lack of a common vocabulary and methodology to describe shape responses. One prominent example here is the interchangeability of the words “globular”, “wide”, “spherical”, “squat”, “round” to describe seemingly similar shape changes for which a visual representation of morphospace extremes ( Figures 2B and 5 ) would provide a better basis of comparison.

The discrepancy between study results also raises the interesting question of whether a universal shape response to changing water chemistry can be expected or if different species will respond in different ways to OA. A crucial point when shell shape responses to a variable carbonate system are to be compared is the careful separation between kinetically, thermodynamically, and stress-induced shape alterations, or in other words, if the phenotypic response is adaptive or non-adaptive. At the simplest level, changes in the seawater carbonate system as brought about by OA result in a reduction of the bio-calcification substrate accompanied by increases in [H ⁺] abundance that may constrain crystal formation at the site of calcification through lowering Ω _{CaCO 3} (Waldbusser et al., 2013; Cyronak et al., 2015). This means that less material or substrate is available to the individual for shell formation and that calcification likely becomes metabolically more expensive as the individual is required to spend more energy on [H ⁺] extrusion (Waldbusser et al., 2013; Tresguerres, 2016). While the latitudinal shape trends observed in our study are purely kinetically induced (in accordance with the kinetic–energetic hypothesis) (Waldbusser et al., 2016) and are therefore a direct indicator for variations in substrate and inhibitor abundance, we hypothesise that shape variation in OA laboratory studies is likely a mix of kinetic-, dissolution-, and stress-induced phenotypic responses precluding any direct comparisons.

The literature provides a wealth of information about local effects on N. lapillus shell shape. Shell aspect ratio variations in relation to local wave exposure regimes were first reported by (Cooke and Reed, 1895). The authors showed that N. lapillus on exposed shores exhibited smaller, squatter shells with a larger aspect ratio, whereas individuals from protected shores formed narrower, more elongated shells. This trend was later explained to be a response that improved fitness as smaller squatter shells and a proportionally larger foot of exposed individuals (obvious from relatively larger aperture sizes) conferred an advantage over more elongated shell forms with smaller foot sizes in high wave energy environments (Kitching and Ebling, 1967). However, the exposed morph showed a marked disadvantage on sheltered shores as animals of that shell form could barely retract into the shell, making them more vulnerable to predators and desiccation during low tide (Crothers, 1985b). Crothers (1973) found the shell height to aperture height ratio to be a convenient proxy for wave exposure trends for N. lapillus populations on Pembrokeshire, UK shores and thereby quantified the relationship between wave exposure and relative aperture size (Crothers, 1973). However, unlike the reported wave exposure induced shape changes, mean shape-PCs in this study showed no relationships with aperture size ( Figure 6 ). Likewise, we also did not observe a change in apertural teeth expression across studied populations which would indicate variations in wave exposure (Crothers, 1971) or predatory pressure (Appleton and Palmer, 1988) constraints on shell deposition in the studied N. lapillus populations. This suggests that shape variations, as explained by shape-PC1 to shape-PC3 (over these large spatial scales), were largely decoupled from local wave exposure and predation gradients.

While shape variations are an excellent proxy to identify new or changing environmental stressors in a population or species (Marshall et al., 2015), predictions of the long term consequences of plastic adaptions are not trivial - or to use the words of R. J. Berry and J. H. Crothers: “Natural populations of animals rarely behave like suspensions of inanimate matter in a frictionless fluid, responding to and predictable by simple mathematical theory” (Berry and Crothers, 1974). However, newly introduced stressors and associated phenotypic adaptions always represent a trade-off between different optima and can thus challenge an organism’s position in an ecosystem (Ghalambor et al., 2007). Our findings suggest that a trend towards narrower, more elongated shells can be expected under increasing OA conditions with, perhaps, significant implications for the predator-prey relationship on North Atlantic rocky shores. Although shape-PCs were largely decoupled from local wave exposure gradients, a continuing trend towards the formation of narrower and more elongated shells could have conceivable consequences for populations on (extremely) exposed coasts. In particular, since global climate change is associated with increasingly severe storm events IPCC 2014, the larger surface of elongated shells could lead to more frequent dislodgement and, as a consequence, could make N. lapillus populations more susceptible to predatory attacks. The herein observed shell shape response to OA could lead to a disintegration of the delicate balance between shape responses to wave exposure or predatory pressure extremes with the potential consequence of a realignment of rocky-shore community structures. However, our results and predictions of morphological responses to end century OA conditions should be viewed with caution for three reasons: First, the site-specific intra-population variability is many times larger than the average shift in shape observed in this study ( Figure 2 ), suggesting that there is a fair amount of tolerance involved in the shape response to environmental constraints. Second, a sensitivity analysis of shell morphology to different stressors is missing, which would be required to make robust predictions. Third, temporal and spatial constraints of environmental data used in this study limited fine-scale observations. While our models suggest a remarkable sensitivity of shell shape to variations in the carbonate system on a large geographic scale, it is not clear if shell morphology will be more sensitive to water chemistry or site-specific physical constraints under increasing OA conditions.

Irrespective of the direct consequences, our results show that changes in the carbonate system are likely to directly affect N. lapillus shell morphology. Although calcification conditions were, in classical terms, favourable at every sampling site, shell formation and resulting morphology were sensitive to changes in the SIR. Our findings uncovered a, so far unrecognised, large scale pattern in shell morphology of the intertidal N. lapillus, which highlights the coupling of phenotypic adaptions to variations in the carbonate system. However, although our findings suggest an adaptive trend in shell morphology under increasing OA conditions, further research is required to understand the phenotypic sensitivity to local and global constraints, which could present an exciting new research opportunity to understand shell formation constraints in a changing world.

The code and data that support the findings of this study are openly available on Zenodo at https://doi.org/10.5281/zenodo.6520504.

DM, LP, and EH conceived the original project and designed the study. DM performed the specimen collection, laboratory and modelling work, and analysed the data. DM wrote the first draft of the manuscript and all co-authors contributed substantially to revisions. All authors contributed to the article and approved the submitted version.

This study was supported by a NERC studentship awarded to DM (NE/L002507/1).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.