Evidence for Carbonate System Mediated Shape Shift in an Intertidal Predatory Gastropod

Phenotypic plasticity represents an important first-line organism response to newly introduced or changing environmental constraints. Knowledge about structural responses to environmental stressors could thus be an essential measure to predict species and ecosystem responses to a world in change. In this study, we combined morphometric analyses with environmental modelling to identify direct shape responses of the predatory gastropod Nucella lapillus to large-scale variability in sea surface temperature and the carbonate system. Our models suggest that the state of the carbonate system and, more specifically, the substrate inhibitor ratio ( [ H C O 3 − ] [ H + ] − 1 ) (SIR) has a dominant effect on the shell shape of this intertidal muricid. Populations in regions with a lower SIR tend to form narrower shells with a higher spire to body whorl ratio, whereas populations in areas with a higher SIR form wider shells with a much lower spire to body whorl ratio. These results indicate that a widespread phenotypic response of N. lapillus to continuing ocean acidification can be expected, potentially altering the phenotypic response pattern to predator or wave exposure regimes with profound implications for North Atlantic rocky shore communities.


INTRODUCTION
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 3 ) 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 2 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 3 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 3 substrate abundance and speciation, i.e., to the CaCO 3 saturation state W 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 W CaCO 3 is a function of both calcium ion [Ca 2+ ] and bicarbonate ½CO 2− 3 availability as well as the thermodynamic solubility product K sp of the respective CaCO 3 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 2+ ] 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 W CaCO 3 < 1, CaCO 3 dissolution becomes thermodynamically favourable; however, as pointed out by Waldbusser et al. (2016), low W 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 W 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 topdown 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.

Shell Collection
Specimens used in this study were sampled from the extensive Crothers collection at the Natural History Museum, London, United Kingdom (NHM) ( Table 1). These samples were collected by John Crothers and colleagues during the 1970s-1980s from numerous sites to study shell shape and colour changes with a special focus on wave exposure levels (Crothers, 1971;Crothers, 1979;Crothers, 1980;Crothers, 1981;Crothers, 1982;Crothers, 1983;Crothers, 1985a). From this collection, we sub-sampled a total of 1575 specimens from 19 sites spanning nearly the entire latitudinal range of N. lapillus on the European North Atlantic and Arctic coast ( Figure 1). We only included specimens with intact apices that showed no extensive signs of erosion and that had likely reached sexual maturity during their lifetime (Galante-Oliveira et al., 2010).

Environmental Descriptors
To study the effect of environmental constraints on shell morphology, we used data from the global oceanin-situ reprocessed carbon observations product (Copernicus Marine Service). This comprehensive global data set contains a compilation of in-situ observational gridded data collected from 1950 until 2020 for key parameters such as water temperature, salinity, alkalinity and pH from two up-to-date carbon and biogeochemistry data products, namely the Surface Ocean Carbon Atlas (SOCATv2020) and the Global Ocean Data Analysis Project (GLODAPv2.2020) that underwent rigorous data quality and bias checks (Bakker et al., 2016;Olsen et al., 2020). The data has a 1x1 degree spatial resolution at 33 different water depths. For this study, we used only sea surface data to reflect the conditions in the intertidal zone. Seawater conditions used in the further analysis were inferred from the nearest data grid to each sampling site. To have access to all carbonate system parameters, we calculated missing values using pH and alkalinity as input parameters to the function carb() in the R package seacarb v3.3.0 (Lavigne et al., 2008). The default dissociation constants (K1 and K2) for the respective temperature and salinity ranges were used in the calculation (Waters and Millero, 2013).

Morphometric Analysis
We measured standard dimensions, namely shell height, shell width, aperture height and aperture width, using digital Vernier callipers. Every dimension was measured three times per shell. Aperture size was estimated from height and width measurements using the equation for the area of an ellipse. Continuous shell shapes (outlines) were obtained using elliptic Fourier analysis (EFA) (Giardina and Kuhl, 1977;Kuhl and Giardina, 1982) using an adaptation of the methodology laid out elsewhere (Telesca et al., 2018). For this, shells were individually positioned on an illuminated platform with the aperture facing downward. Images were taken from a fixed distance with a DSLR Camera (Nikon 3000, Sigma 105 mm macro lens) mounted on a photo stand in a dark room. This produced images of the shells as black silhouettes with clear outlines, ideal for automatic outline tracing. We decided to analyse shell shape variability with the shell apertures facing downward in a natural position, rather than in the more common but unnatural aperture upwards position, to study shape changes that are more meaningful for N. lapillus in nature. Shell silhouette images were imported to R (R Core Team, 2020), and outlines were digitised and subsequently turned into closed shape polygons (defined by x-y coordinates) using the R package Momocs (Bonhomme et al., 2014). Shape outlines were smoothed using a running average smoothing algorithm to remove digitisation noise before each shell was individually aligned along its longitudinal axis. This was done by defining a secant line that connected the two furthest apart points (i.e., the shell anterior with the posterior) within the polygon shape. The shape was then rotated until the slope of the secant line was zero. Thereafter, shape polygons were centred and scaled with regard to orientation and size. An EFA with 10 harmonics encompassing 99% of the total harmonic power was computed on the outlines. Principal component analysis (PCA) was used to identify the axes of most variance among individuals (hereafter shape-PCs). For further analyses, we used the first five shape-PCs that encompassed more than 86% of the total shell shape variance among specimens as representative shape descriptors. A graphical representation of the shell shape analysis pipeline can be obtained from the Supplementary Material. To visualise the transition between extremes of the morphospace, we generated deformation grids (Thompson, 1992) and iso deformation lines (Bookstein, 1992) for each shape-PC. Individual relationships of the first five shell shape descriptors with sampling site latitude were analysed using simple linear regression analyses (Figure 2).

Data Analysis and Model Formation
The effect of key environmental parameters on shell shape variability was investigated using generalised additive mixed models (GAMM, mgcv) (Wood, 2011) to accommodate the hierarchical structure and non-linear relationships in our dataset. Initial data exploration was done following Zuur et al. (2010). No significant outliers were detected among shape-PCs. Boxplots of shape-PCs by sampling sites showed heterogeneous variances between shape-PCs which was an expected consequence of the PCA. Since the aim of this work was to study the effect of global change on shell formation, we considered a range of associated environmental parameters of the physical and chemical domain as model predictors. As stated earlier, W CaCO 3 states are indicators for spontaneous carbonate precipitation and, as such, have been shown to be good predictors of bio-calcification (Sanders et al., 2021). However, apart from the theoretical framework, there is little evidence supporting its mechanistic role in molluscan shell formation (Bach, 2015). There is growing evidence that hydrogencarbonate concentration ½HCO − 3 rather than ½CO 2− 3 is likely the primary carbon species used in bio-calcification (Herfort et al., 2002;Allemand et al., 2004;Mackinder et al., 2010;Stumpp et al., 2012). Thus, we also included the substrate inhibitor ratio (SIR) along with the other parameters as possible predictors to describe variability in the carbonate system that is likely more mechanistically relevant to biological calcification (Jokiel, 2013;Bach, 2015;Cyronak et al., 2015). The SIR is a ratio of ½HCO − 3 concentration by proton availability (H + ) in mol μmol -1 (eq. 2) and, as such, is a representative of the inorganic carbon substrate availability in relation to precipitation conditions.
The predictors included in the initial model were: temperature, salinity, total alkalinity (A T ), pH, ½HCO − 3 , W Ar , W Ca , SIR, sampling site ID and shell height. We selected the first three shape-PCs that exhibited a significant correlation with latitude ( Figure 2) to be included in the GAMM. Given that all three selected shape-PCs are shape-representatives of the same shell, normalised shape-PC1 to shape-PC3 descriptors were included simultaneously within the same GAMM to allow for a holistic description of the shell shape variability. The number of knots per descriptor was manually defined to constrain A B FIGURE 2 | (A) Latitudinal trends of shell shape descriptors (shape-PC1 to shape-PC5) of the intertidal gastropod Nucella lapillus. Points represent the mean at every sampling site, and error bars denote the standard deviation. Linear regression significances are reported for each shape-PC (Significance levels: n.s. p > 0.05, *p < 0.05, **p < 0.001, ***p < 0.0001). (B) Deformation grids of shape-PC1 to shape-PC5 visualising the transition between extremes of the morphospace represented by each shape-PC (blue: mean -3s, red: mean + 3s). Percentage contributions to total shape variance between individuals are reported for each shape-PC. unnatural degrees of wiggliness, and cubic regression splines were used. The dredge() function in the R package MuMIn (Bartoń, 2020) was used to analyse all possible combinations of predictors, and the Akaike Information Criterion (AIC) was used as the principal tool in the model optimisation process. To overcome potential complications with collinearity in the model caused by correlated independent variables only models with an acceptable VIF factor < 3 of covariates were considered. The sampling site ID was included as a random effect term to account for the dependency of specimens from the same collection site. The best-supported model was then fitted using the restricted maximum likelihood approach (REML) and validated by inspection of standardised residual patterns to verify model assumptions. Variogram inspection suggested no significant spatial autocorrelation. The final model was of the form: where shapePC ij represents the j th observation for each shape-PC (j = 1-3) at every sampling site i (i = 1-19). f() * fshape denote the smoothing functions of temperature, alkalinity, SIR and shell height expressed by individual smoother for each level of j and ∈ i is the normally distributed random intercept. Effect size estimates of standardised predictors were calculated from the GAMM following Telesca et al., (2018). The standardisation of model predictors was done by subtracting the predictor's mean from its standard deviation. 95% confidence intervals were constructed along mean effect sizes to compare the magnitude and significance of environmental predictors on shell shapes. All exploratory data analyses and modelling were performed in R (R Core Team, 2020) and all data and scripts are publicly available (see data availability statement).

RESULTS
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 W Ca and W Ar values at the sampling sites were 3 and 2 respectively. Seawater W 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 nonsignificant 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 nonsignificant for any of the three shape descriptors.

DISCUSSION
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 nonsignificant effect size for any of the three shape-PCs ( Figure 5). Despite overall high W Ca and W 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" ½HCO − 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 2 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 nonadaptive. At the simplest level, changes in the seawater carbonate system as brought about by OA result in a reduction of the biocalcification substrate accompanied by increases in [H + ] abundance that may constrain crystal formation at the site of calcification through lowering W 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 FIGURE 3 | Scatter plot of the first two shape-PCs from PCA analysis performed on elliptic Fourier coefficients revealing the among-population variation of outline and shape features between Nucella lapillus specimens. Confidence intervals are shown for each origin group (ellipses with dashed linetype), and corresponding morphospace extremes are represented as grey shapes.
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 A B C FIGURE 4 | GAMM shape trend predictions of the best-supported model as expressed by (A) shape-PC1, (B) shape-PC2, and (C) shape-PC3 with selected descriptors of sea surface temperature, A T , substrate inhibitor ratio (SIR) and normalised shell height. Significance is reported for each predictor (Significance levels: n.s. p > 0.05, *p < 0.05, **p < 0.001, ***p < 0.0001). Colour gradients (red to blue) of significant predictor terms correspond to those in Figure 5. 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 A B C FIGURE 5 | (A) shape-PC1, (B) shape-PC2, and (C) shape-PC3 and the best-supported model's mean effect size plots, including sea surface temperature, A T , substrate inhibitor ratio (SIR) and shell height. Error bars represent the 95% CI, and the significance of either predictor is determined if the CI does not cross the zero line (p < 0.05). Asterisks mark significant predictors. Colour gradients (red to blue) correspond to the morphospace extremes visualised to the left of each effect size plot. 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.

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

AUTHOR CONTRIBUTIONS
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.

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