A century of coping with environmental and ecological changes via compensatory biomineralization in mussels

Abstract Accurate biological models are critical to predict biotic responses to climate change and human‐caused disturbances. Current understanding of organismal responses to change stems from studies over relatively short timescales. However, most projections lack long‐term observations incorporating the potential for transgenerational phenotypic plasticity and genetic adaption, the keys to resistance. Here, we describe unexpected temporal compensatory responses in biomineralization as a mechanism for resistance to altered environmental conditions and predation impacts in a calcifying foundation species. We evaluated exceptional archival specimens of the blue mussel Mytilus edulis collected regularly between 1904 and 2016 along 15 km of Belgian coastline, along with records of key environmental descriptors and predators. Contrary to global‐scale predictions, shell production increased over the last century, highlighting a protective capacity of mussels for qualitative and quantitative trade‐offs in biomineralization as compensatory responses to altered environments. We also demonstrated the role of changes in predator communities in stimulating unanticipated biological trends that run contrary to experimental predictive models under future climate scenarios. Analysis of archival records has a key role for anticipating emergent impacts of climate change.


| INTRODUC TI ON
As global climate change accelerates (Kirtman et al., 2013), we need accurate predictions about biological responses if we are to prevent emergent damage to the biosphere (Nagelkerken & Connell, 2015;Urban et al., 2016). However, our ability to project changes to individual species and populations, let alone communities and ecosystems, in response to disturbances is limited (Kroeker et al., 2017;Urban et al., 2016). The overwhelming majority of models forecasting biotic change stem from experimental observations on relatively short timescales (Kroeker et al., 2013;Nagelkerken & Connell, 2015). Long-term data are lacking in natural systems incorporating the role of phenotypic plasticity, both within and across generations, in combination with long-term | 625 TELESCA ET AL. genetic adaption in populations, the processes identified as critical for responding to altered environments (Cross et al., 2018;Thomsen et al., 2017;Vargas et al., 2017). Thus, a mechanistic understanding of processes shaping biotic variation in functioning ecosystems, across relevant timescales to evaluate potential impacts of long-term acclimation and adaption, is critical for building the theoretical framework necessary to anticipate the scope of biological responses (Kroeker et al., 2017;Urban et al., 2016).
Compensatory mechanisms are important mediators of biotic dynamics with the potential to attenuate (Ghedini et al., 2015), or even reverse (Connell et al., 2017), the direct impacts of change on overall performance and ecological functions of organisms (Cross et al., 2019;Ghedini et al., 2015). Compensation can arise from a range of physiological, habitat, and ecological alterations acting at organismal, population (Ashton et al., 2017;Cross et al., 2019), or community level (Ghedini et al., 2015;Peck et al., 2015). These biological mechanisms have, however, seldom been examined in species with both high climate sensitivity and disproportionate ecological impact on healthy ecosystem functioning (McCoy & Ragazzola, 2014;Telesca et al., 2019;Urban et al., 2016), such as calcifying foundation species.
Climate change has caused major alterations in ocean physics and chemistry (Kirtman et al., 2013). These changes, which profoundly affect marine community structure and ecosystem dynamics, have been further exacerbated in nearshore habitats by increasing human pressures over the last century (Nagelkerken & Connell, 2015;Urban et al., 2016). Most environmental change research on marine calcifiers has focused on early life-history stages given the expectation, and largely the findings, that intensifying climate stressors act more strongly on larvae and juvenile stages than adults (Gazeau et al., 2013;Kroeker et al., 2013;Waldbusser et al., 2015). Such projections about calcifiers' susceptibility to change predominantly originate from short term to long term, from days to 1 or 2 years, laboratory (Kroeker et al., 2013;Nagelkerken & Connell, 2015), mesocosm and in situ manipulations (Ashton et al., 2017;Connell et al., 2017;Peck et al., 2015). These studies on simplified experimental 'communities', despite being long term for experiments, may not translate to very long-term population dynamics within real-world naturally complex ecological systems (Connell et al., 2017;Peck et al., 2015;Vargas et al., 2017). Recent in situ works have begun to address this problem (Ashton et al., 2017;Clark et al., 2019), but a rarely used approach to evaluate longer-term responses, from years to decades, and the target of potential population-level selection involves using archival collections (Cross et al., 2018;McCoy & Ragazzola, 2014;Pfister et al., 2016). Although historical studies on archival material are correlative in their nature and might, therefore, be constrained by the limited parallel records of environmental information at the corresponding timescale of the biological processes analysed, they are especially useful to (a) assess ecological effects of stressors as they are introduced over climate change relevant timescales (McCoy & Ragazzola, 2014); and (b) include evaluation of population responses via transgenerational phenotypic plasticity and possible long-term adaption, both keys to resistance (Thomsen et al., 2017;Vargas et al., 2017), complementing experimental data but not incorporated in any other approach (Cross et al., 2018;Vargas et al., 2017).
Marine mussels of the genus Mytilus are bed-forming foundation species throughout the world's eulittoral ecosystems. Owing to their economic value, ecosystem services provided, and projected vulnerability (Gazeau et al., 2013;Thomsen et al., 2017), mytilids have received much attention as key indicator species (Kroeker et al., 2016;Telesca et al., 2018Telesca et al., , 2019. However, because of the ubiquity of mussels, extensive museum collections of Mytilus from the last century are exceptionally rare. Calcareous mussel shells (Figure 1a,b) perform several vital functions including structural support and protection from predators. Because shell integrity determines survival, shell traits, and their resulting mechanical properties are subject to strong selection pressure, with functional success or failure a fundamental evolutionary driving force (Freeman, 2007;Watson et al., 2017). Thus, an understanding of potential compensatory mechanisms in shell biomineralization to withstand environmental insult and altered predation impacts is essential to inform realistic projections of calcifiers' persistence and ecological functions when selection pressures alter in changing environments.
Here, we describe temporal changes in biomineralization in the blue mussel Mytilus edulis using a unique museum collection of specimens sampled from intertidal breakwaters along 15 km of Belgian coastline at near decadal frequency between 1904 and 2016. By combining extensive long-term datasets of key environmental descriptors and predators in a naturally controlled system, across a centennial timescale incorporating long-term acclimation and adaption, we describe (a) compensatory trade-offs in shell biomineralization as potential mechanism for resistance to altered environments and species interactions; and (b) the role of changes in predator communities in stimulating unanticipated trends in biomineralization that run contrary to global-scale projections under future climate scenarios.
Specimens collected between 1904 and 1987 were obtained from archival shell collections of the Royal Belgian Institute of Natural Sciences (RBINS). This unique collection of a single species is composed of both wet (shells and body tissue preserved in 70% ethanol) and dry (shells only) specimens that were obtained from monitoring programmes during the 20th century. Only archival specimens with a well-preserved periostracum and no evident deterioration of internal calcareous shell layers, as well as with information on species identity, collection date and location were evaluated. Only general information on the tidal ranges of specimens (i.e. intertidal or subtidal) were reported in collection notes for most of the sampling sites; therefore, exact tidal height could not be obtained. To minimize the potential influence of differences in tidal habitats on mussels shell formation (i.e. differences in growth rate; Seed, 1969), we selected adult M. edulis with shell lengths of 39-66 mm, collected from intertidal stone breakwaters. Specimens of this size class were selected to (a) avoid bias in the analyses due to the use of shells with different size ranges for each collection site; and (b) control statistically for the effect of shell size variations on shell layers deposition. In all, 30 intertidal adults (shell length of 44-62 mm) were hand-collected in 2016 from a single location and substratum (intertidal stone breakwaters in Mariakerke). For each specimen, shell length was measured with digital calipers (0.01 mm precision) and used to control for within-year shell size differences (Telesca et al., 2018). Animal age was estimated by counting annual shell increments. Because of the limited control we had on the selection of specimens from specific intertidal heights, we tested for age differences in mussels of similar shell size between collection sites/years (indicative of potential differences in growth rate due to habitat heterogeneity).

| Shell preparation and layer measurement
We set left shell valves (n = 256) in polyester resin (Kleer-Set FF, MetPrep). Embedded specimens were sliced longitudinally along their axis of maximum growth (Figure 1a,b) using a diamond saw and then progressively polished with silicon carbide paper (grit size: P800-P2500) and diamond paste (grading: 9-1 µm). Shell thickness was measured on polished cross-sections photographed with a stereomicroscope (Leica M165 C equipped with a DFC295 HD camera; Leica). Because many shells had undergone evident environmental abrasion, which removed the periostracum and prismatic layer near the anterior shell side (umbo), we estimated whole-shell, prismatic, and nacreous layers thicknesses at the midpoint along the shell crosssection (Telesca et al., 2019). Periostracum thickness was measured at the posterior edge where it attaches to the external surface of the prismatic layer. In this way, we estimated the fully formed organic layer that was unaffected by decay and abrasion caused by physicalenvironmental agents during shell growth or potential deterioration (detachment) due to preservation (Harper, 1997; Figure 1b).

| Organic content analysis
We performed thermogravimetric analyses (TGA) to estimate the weight proportion (wt%) of organic matrix within the prismatic layer. M. edulis specimens were selected from 3 years, 1904, 1950, and 2016 (n = 14, 14, and 15, respectively), to explore differences in shell organic content. The periostracum was removed by sanding, and prismatic layer tiles (8 × 5 mm, n = 43) were cut along F I G U R E 1 Historical patterns of blue mussel shell thickness along the Belgian coast between 1904 and 2016. (a) Mytilusedulisshell morphology and structure. (b) Antero-posterior cross-section of a left valve showing the arrangement of shell layers. The shell consists of three layers: (i) the outer periostracum; (ii) the fibrous prismatic layer; and (iii) the nacreous layer. The organic periostracum is made of scleroproteins (quinone-tanned). The calcified fibrous prismatic and nacreous layers are mainly composed of calcium carbonate (CaCO 3 ) crystals of different mineral forms, calcite and aragonite, respectively, and minor amounts of inter-and intra-crystalline organic matrix. (c) Predicted temporal trends of shell layer thicknesses for the mean shell length of the sample (53.13 mm) between 1904 and 2016. Whole-shell and nacreous layer thickness increased linearly over time. Prismatic layer deposition significantly decreased between 1970and 1983, and increased in the 1911-1934period and after 1995. Periostracum thickness increased between 1945and 1962  the posteroventral shell margin, cleaned, air-dried, and then finely ground. We tested 10 mg of this powdered shell with a thermogravimetric analyzer (TGA Q500; TA Instruments). Samples were subjected to constant heating from ~25°C to 700°C at a linear rate of 10°C/min under a dynamic nitrogen atmosphere and weight changes were recorded (see Methods S1). We estimated the wt% of organic matter within the shell microstructure as weight loss percent during the thermal treatment between 150°C and 550°C ( Figure S2).

| Elliptic Fourier analysis of shell outlines
Shape analyses of M. edulis shell valves were carried out through a geometric morphometrics approach using an elliptic Fourier analysis (EFA) of outlines. Outline processing and EFA were carried out using the package Momocs (Table S2) in the software R (v3.5.2; R Core Team, 2018).
Outlines of left shell valve lateral views (n = 268) were digitized, processed, and analysed with an EFA following Telesca et al. (2018; see Methods S1). We chose seven harmonics, encompassing 99% of the total harmonic power ( Figure S3). Four coefficients per harmonic (28 descriptors) were extracted for each shell outline and used as variables quantifying geometric information.
A principal component analysis (PCA) was performed on the harmonic coefficients to define axes capturing most of the shape variation among individuals. Principal components (PCs) were considered as new shape variables. To understand the contribution of harmonic coefficients to shell shape, we reconstructed extreme outlines along each PC ( Figures S3 and S4). The first seven PCs (99% of outline variation) were analysed with a MANOVA to test for significant effects of collection year and shell length (size) on shape variances ( Figure S4). To visualize outline differences, we generated deformation grids and iso-deformation lines through mathematical formalization of thin plate splines analysis ( Figure S3).

| Environmental dataset
To understand temporal changes in marine conditions along the Belgian coast, we evaluated available long-term datasets of key environmental descriptors for the collection location: sea surface temperature (SST), sea surface salinity (SSS), chlorophyll-a (Chl-a) concentration (as a proxy for food supply; Thomsen et al., 2013) and dissolved oxygen (see Data S1). Key descriptors were selected based on known influences on mussel growth, their availability for the sampling location, and forecasted ocean alterations under climate change (Kirtman et al., 2013).  (Telesca et al., 2018(Telesca et al., , 2019. Other parameters and local process with a potential influence on mussel biomineralization could not be included due to their limited availability or absence of historical datasets for the collection sites.

| Mussel-predators dataset
To describe historical changes in predation regime and pressure on intertidal mussel beds and their shell traits, we analysed long-term abundance datasets of key intertidal predators with a potential influence on mussel shell traits: decapods (benthos and planktonic larvae), seagulls and dog whelks. We selected key predator taxa based on historical dataset availability for the study location, their predation strategy (i.e. durophagy and drilling), and known effects on shell biomineralization responses in Mytilus spp. (Freeman, 2007;Lowen et al., 2013;Sherker et al., 2017). The predatory starfish Asterias spp. was excluded from our study due to (a) the absence of historical datasets for the study location; (b) the reported changes in echinoderm abundances and ranges at a North Sea scale after the mid-1980s mainly regard echinoids and ophiuroids ); (c) their limited influence on intertidal mussel compared to subtidal mussel beds (Seed, 1969;Seed & Suchanek, 1992); and (d) the inducible defence response in Mytilus spp.
shells of this predator that pulls the shells apart, including development of stronger adductor muscles, increased byssal strength, and occasional trade-offs between growth and shell thickness with no net effect on shell deposition (shell weight ;Freeman, 2007;Lowen et al., 2013).

| Statistical analysis
Generalized additive mixed-models (GAMMs; Wood, 2017) were used to describe long-term variations in environmental descriptors and predator abundances, as well as temporal trends in mussel biomineralization and shell shape. GAMMs allowed us to (a) account for the hierarchical structure of the datasets (multiple M. edulis specimens from each collection site; n = 12-30); (b) specify flexible dependence structures of the response on the covariates; and (c) control for variations ('noise') among sampling locations (Bolker, 2015).
We carried out data exploration following the protocol of Zuur

| Time-series analysis
We used GAMMs to model long-term and seasonal changes in time series of environmental descriptors and the abundance of predators. Time series, consisting of discontinuous daily observations per year, were expressed as continuous monthly averaged measurements. Fixed covariates were month (seasonal variation, expressed as numeric 1, …, 12 indicator), year (trend), both fitted as smoothers, and month × year interaction that was defined through a smooth tensor product interaction (Wood, 2017). The tensor product interaction is used to represent functions of covariates, which are Boundaries and numbers of knots (limits and dimensions of the basis used for splines) were manually selected while effective degrees of freedom were estimated by the smooth function (Wood, 2017). A cubic regression spline basis was used for the trend smooth term, whereas a cyclic cubic spline basis was used for the seasonal smooth term to avoid discontinuity between January and December values (Wood, 2017). Models indicated significant within-year residuals autocorrelation, which required the use of different residual correlation structures (Wood, 2017). Autoregressive (for equally spaced time series) or conditional autoregressive (for unequally spaced observations) models, nested within each year, were fitted to the residuals to account for temporal autocorrelation.
A negative binomial generalized additive model with a log link function was used to model numbers of seagull breeding pairs as a function of year and to compare between the four species. The log link function ensures positive fitted values, and the negative binomial distribution is typically used for count data that are overdispersed with respect to a Poisson distribution. Covariates were year (smoother with a cubic spline), gull species (categorical, four levels: seagull species), and their interaction.

| Statistical modelling
Separate GAMMs were used to describe shell thickness (n = 256) and shape (n = 268) with respect to time, and to compare between shell measurements (thickness of whole-shell, periostracum, prismatic, and nacreous layers) and morphological traits (PC1 and PC2 from EFA coefficients). Fixed covariates were collection year and shell length (both fitted as smoothers). To estimate biologically meaningful intercepts, we centred the input variables. Collection year was centred to the year 1900, while length to sample mean shell length (53.13 mm). Shell length was included in all models to control for the effect of shell size variations on shell measurements and morphological traits. To incorporate dependency among observations from the same sampling site, we used site as a random intercept.
where Response ij is the jth observation for the response variable (thickness or shape variable) in site i (i = 1, …, 13), f is the smoothing function (cubic regression spline), and β 0 is an intercept. ε ij is the noise and Site i is the random intercept, which are assumed to be independently distributed with expectation 0 and variance σ 2 and 2 Site , respectively.
We modelled the wt% of organic matrix (n = 43) as a function of year of collection (categorical, three levels: 1904, 1950, and 2016) and prismatic layer thickness (continuous), to test for differences among time periods and association with shell deposition.
The response variable was coded as a value from 0 to 1; therefore, a generalized linear model with a beta distribution and a logistic link function was used. Pairwise contrasts with a standard Bonferroni correction were then used to test (α = 0.017) for differences in the wt% of organics among collection years. (1) To estimate potential interactions of local environmental conditions on mussel shell biomineralization, we first calculated descriptive statistics for individual factors per year (i.e. mean, median, standard deviation, 10th, 25th, 75th, and 90th percentiles, minimum and maximum annual values). To provide a first-order approximation of the water conditions during the lifespan of sampled specimens, we expressed the descriptive statistics as average values for the 6-year period prior to collection (4 years for individuals from 1904).
We then standardized all the descriptors to account for differences in measurements units by subtracting the sample mean from the variable values and dividing them by the sample SD Next, we performed separate PCAs on descriptive statistics for SST and SSS to create PCs, characterizing SST and SSS 'regimes' (Kroeker et al., 2016). We then used the first two PCs scores from each PCA as input variables (Table S3), along with PCs from EFAs (shape-PCs; Figure S4) included as proxies for food supply (Telesca et al., 2018), as fixed covariates (smoothers) in GAMMs to model variations in shell thickness measurements. Collection site was included as a random intercept.
To describe variations of shell traits under changing predation regimes, we modelled thickness of each shell layer as a function of predators' abundance (low/high: decapods and gulls) and presence (presence/absence: dog whelks). Thickness measurements were selected depending on the temporal availability of predator datasets.

| Model optimization and predictions
Models were optimized by first selecting the random structure and then the optimal fixed component. The principal tools for model comparison were the corrected Akaike information criterion (AICc) and likelihood ratio tests. Random terms were selected on prior knowledge of the dependency structure of the dataset.
Visual inspection of residual patterns indicated violation of homogeneity in most cases. This required the use of variance structures (generalized least squares) allowing the residual spread to vary with respect to predictors. To assess the presence of temporal dependence in model residuals, we used (partial) autocorrelation functions, for regularly spaced time series, or variograms, for irregularly space ones. When temporal autocorrelation was found, the best residual correlation structure was selected through comparisons of models with differing stochastic trend models in the residuals. The fixed component was optimized by rejecting nonsignificant interaction terms only that minimized the AICc value.
For environmental GAMMs with multiple smooth terms, selection was carried out using cubic regression splines with shrinkage (Wood, 2017). Final models were validated by inspection of standardized residual patterns to verify assumptions of normality, homogeneity, and independence. The proportion of variance explained by the models was quantified with conditional determination coefficients (cR 2 ).
Optimal models were used to describe changes in mussel shell characteristics over the last century. All predictions for temporal changes of shell thickness and morphological traits were computed controlling statistically for among-individual variations in shell length. We used a simulation-based approach to identify periods of statistically significant change in the time series (see Methods S1; Figure S5). This information was used to identify which parts of the predicted trend were significantly increasing or decreasing over time. All data exploration and statistical analyses were performed with the software R (Table S2).
Thermogravimetric analyses indicated the weight proportion (wt%) of organic content within the prismatic layer of mod-  Table S4). No significant difference in animal age (number of annual shell increments) between collection years for individuals of same shell size was detected (F 13,247 = 1.85, p = .28).

| Abiotic and biotic influence on biomineralization
Models indicated differential variations of shell traits with SST and SSS regimes (Figure 4; Table S7). Shell deposition was correlated with changing SST (SST-PC1: F 2.9,250 = 5.32, p = .0019; SST-PC2:   Figure S10). This suggests that temporal trends in shell deposition were generally more correlated with changes in environmental variability than altered mean water conditions. We identified variations in the thickness of each shell layer under different predation types (durophagous vs. drillers) and relative/recorded abundances ( Figure 5; Table S8). We observed a 30% increase in deposition of calcareous layers (prismatic: t 342 = 3.94, p = .0001; nacre: t 342 = 6.68, p = .0061) and formation of 45% thinner periostraca (t 342 = 6.68, p < .0001) with higher decapod abundances. Mussels were characterized by 73% thicker prismatic layers  (1914-1918 and 1940-1945) for which we extrapolated our predictions. (Right panels) Perspective plots of the prediction plane resulting from interactions between long-term trend and seasonal components of the data. Note the different scales of the ordinate (response, marks as in the colour scale) and abscissa (year axis), highlighting differences in the magnitudes of change and temporal availability for each parameter.  us to identify unexpected natural responses of shell deposition and to evaluate the differential weighting of marked environmental (physicochemical conditions) and ecological (predators abundance) changes on temporal biomineralization patterns. Ecological impacts via rapid changes in predator communities had a stronger influence than longterm environmental variations on observed compensatory responses in adult mussels that would be expected to follow climate changes according to experimental predictive models.

| Long-term environmental patterns
Our models indicated formation of more variable surface regimes along the Belgian coast over the last century, with temporal patterns depending on the descriptor considered (Figure 3a; Figure S6; Table 1; Table S5). Long-term increases of SST with TA B L E 1 GAMMs summary statistics of environmental descriptors and predator abundance. F-values and the significance (n.s., p > .01; *, p < .01; **, p < .001; ***, p < . warmer spring and summer temperatures are in line with global ocean trends (Desmit et al., 2020;Kirtman et al., 2013). No change in long-term and seasonal SSS patterns was detected, reflecting expected long-term stability in local water cycles for areas with intermediate salinity (mid-latitudes; Kirtman et al., 2013). Increasing Chl-a concentration along the Belgian coast has been widely correlated with growing human activity and land use over the last 60 years, enhancing inorganic nutrients and organic carbon river loads (Desmit et al., 2020;Gypens et al., 2009). These processes have led to increased coastal eutrophication that boosted primary production (Borges & Gypens, 2010;Gypens et al., 2009). The observed changes in timing and onset of Chl-a peaks, with formation of an earlier annual peak in the last two decades is consistent with Desmit et al. (2020). Lower dissolved oxygen concentrations over time are in line with higher SST (decreased oxygen solubility) and primary production, overall increasing biological activity ( Figure 3a).
Our results suggest that formation of more variable sea surface regimes of SST over time has produced energetically more demanding conditions (e.g. from increased basal metabolic costs at higher temperatures) for intertidal mussels (Bayne, 1976;Thomsen et al., 2013;Watson et al., 2017). Local-scale alterations may, however, benefit indirectly calcification via increased food supply (higher Chl-a concentration) and buffered ocean acidification (Connell et al., 2017;Thomsen et al., 2013; Figure 6). From an historical perspective, observed shell changes, with no evident alteration of shell growth, suggest current mussel populations may have acclimated or adapted to altered environmental conditions along the Belgian coast.

| Historical change of predation regime
Our results indicate increasing predation pressures on M. edulis with a rapid shift to shell-breakers (durophagy) dominated regimes (i.e. decapods and seagulls) from the early 1990s on that followed a decline in predation from the local extinction of a keystone predatory driller (dog whelks) at the end of the 1970s (Figures 5 and 6; Table 1).

Long-term increases in abundance and range of decapods in the
North Sea  have been explained by a temperature-driven, abrupt ecosystem shift during the 1980s . These changes reflect increasing decapod predation on fish recruits and benthic bivalves Lord et al., 2017), that have both decreased in numbers in the 1990s and 2000s ; is significantly controlled by prey abundance and composition (Schwemmer & Garthe, 2005), with decapods, pelagic fishes, and fishing discard being dominant components (Votier et al., 2004).
A greater abundance of decapods  may, therefore, benefit seagulls (Luczak et al., 2012) and especially chick development during the breeding season (May-June; Schwemmer & Garthe, 2005), which coincides with the decapod peak abundance (Figure 3). The North Sea supports important fisheries, generating substantial bycatch and discard for seabirds. However, food depletion due to declining volumes of discard (Votier et al., 2004), reduced fish recruitment , and overfishing , has been reported to drive prey-switching tendencies in generalist scavenging birds, such as seagulls, increasing predation on coastal taxa including mytilids (Votier et al., 2004).
Seagulls have the potential to remove up to 70% of natural and cultured mussel stocks (Nehls et al., 1997). Hence, the population increase in Belgian seagulls (Luczak et al., 2012), coupled with enhanced prey-switching behaviour (Votier et al., 2004), are key factors boosting predation on coastal resources, among which are the dominant, easily accessible intertidal mussel beds. This illustrates how climate and anthropogenic-driven marine changes can extend from the avian fauna and propagate to intertidal and terrestrial food webs around seabird colonies (Luczak et al., 2012;Votier et al., 2004; edulis under future ocean acidification scenarios (Sadler et al., 2018).
Thus, their rapid disappearance would have changed significantly predator-prey coastal dynamics (Figure 6), alleviating predation pressure on mussels.

| Triggers of shell compensation
Over the last century, human impact was a major driver of Belgian coastal dynamics with an end-to-end ecosystem effect, acting through direct and indirect pathways on the physical-chemical environment, as well as across trophic levels and species interactions governing the blue mussels' ecological network (Figure 6). A limited influence of environmental water regimes on shell deposition ( Figure 4) suggest ecological alterations of predator communities in a highly productive system were the principal triggers of the observed responses in biomineralization.
Water temperature is a key factor in controlling growth and physiological activity in marine calcifiers (Bayne, 1976;Gazeau et al., 2013). Previous research demonstrated rising temperatures within a species ecological range can increase calcification in F I G U R E 5 Temporal responses of blue mussel shell biomineralization to changes in predator pressure and strategy. Deposition of individual shell layers (prismatic, blue; nacreous, red; periostracum, yellow) under changing abundance of durophagous and drilling predators. (a) Variations under low (1970-1988) and high (1989-2016) decapod abundances (n = 117 × layer) indicating deposition of thicker calcareous layers, and thinner periostraca under high decapod abundances. (b) Shell deposition before (1980-1990) and after (1991-2016) the exponential increase in the number of gull breeding pairs (n = 70 × layer) showing thicker prismatic layers and thinner periostraca under increased gull foraging. (c) Composition under presence and absence (before and after 1981, respectively) of Nucella lapillus (n = 117 × layer) showing decreased thickness of prismatic layers and periostraca after the disappearance of dog whelks. Note the dualy-axis highlighting scale differences between shell layers. Error bars indicate 95% CIs (*p < .01; **p < .001; ***p < .0001)  showing optimal growth between 10°C and 20°C (Bayne, 1976).

Studies on M. galloprovincialis grown under acidified (variable seawa-
ter CO 2 partial pressure, or pCO 2 ) and warmer waters have reported an average +2.63% and +2% increase in shell deposition per unit size for a 1°C increase under low (~400 μatm) and high (~1,200 μatm) pCO 2 conditions, respectively (Kroeker et al., 2014). Moreover, a large-scale study (Telesca et al., 2019) predicts an increase of +3.3% in both prismatic and nacreous layer thicknesses for a 1°C increase under SSS and Chl-a conditions recorded along the Belgian coastline.
Although increases in average annual (+1°C) and summer (+2.7°C) SSTs along the Belgian coast may have contributed positively to shell calcification, temperature patterns alone explain only a small proportion of historical increases in biomineralization (a total of +57% between 1904 and 2016).
Food supply is potentially the most important factor supporting growth in mussels (Bayne, 1976;Thomsen et al., 2013), but this strongly covaries with water temperature. Increased food supply (Chl-a) supports Mytilus calcification (Thomsen et al., 2013), and increased organic content in shells (Telesca et al., 2019) and allows sustained growth over a wide range of environmental temperatures (5°C-20°C; Bayne, 1976). Although increased Chl-a concentrations in response to eutrophication processes along the Belgian coast (Figure 3a) could be beneficial to mussels (Thomsen F I G U R E 6 Inferred historical evolution of the blue mussel ecological network. Conceptual representation of end-to-end ecosystem interactions (considered in the current study) governing blue mussels beds along the Belgian breakwater system, including primary producers, and primary, secondary, and tertiary consumers from both marine (sublittoral and eulittoral) and terrestrial systems. The graph represents ecological top-down and bottom-up control pathways (blue arrows) through which abiotic and biotic interactions within the Mytilus edulis ecological network can stimulate compensatory responses of shell biomineralization. Red arrows indicate external environmental (ENV) and anthropogenic (AI) impacts. Signs indicate positive (+) or negative (−) responses. Solid and dashed arrows indicate direct and indirect effects, respectively. Thinner lines indicate likely but speculative interactions. Compensatory adjustments of blue mussel shells represent the shifting balance between responses to propagating climatic and anthropogenic impacts through the durophagous linkages (increasing decapod larvae → decapod adults → seagulls), the disappearance of drillers (dog whelks), and energetically more demanding environmental regimes (increasing sea surface temperature [SST] and fluctuating sea surface salinity [SSS]) in an increasingly productive coastal system (boosted primary productivity and chlorophyll-a); and their chain of direct and indirect feedbacks straightened or weakened by climatic-and anthropogenic-driven stressors.    , 2013), the positive correlation between food supply and shell growth decreases with increasing mussel size, and has no effect on adults larger than 40-50 mm (Page & Hubbard, 1987) that represent the lower tail of the shell size classes included in our study. This correlation is also supported by Telesca et al. (2019) who reported no relationship between Chl-a levels and the deposition of calcareous shell components in adult mussels (>35 mm) across latitudes.
Bivalve calcification is particularly sensitive to changes of marine carbonate chemistry resulting from the increased absorption of atmospheric CO2 by the ocean. CO2 uptake increases the dissolved inorganic carbon (DIC), in particular seawater pCO2 and HCO − 3 , controlling both pH and the thermodynamic mineral solubility (saturation state). Increasing pCO2 decreases seawater CO 2− 3 and pH, causing a corresponding decrease in CaCO3 saturation state and acidification of the oceans, that interfere with calcification in marine organisms (Gazeau et al., 2013;Thomsen et al., 2015;Waldbusser et al., 2015). Changes in CO 2− 3 and, therefore, saturation state impact bivalve calcification directly, through dissolution or by preventing early shell formation (Thomsen et al., 2015;Waldbusser et al., 2015). Differently, reductions in seawater pH and increased pCO2 impact calcification indirectly via disturbance of physiological processes and reductions in scope for growth (Thomsen et al., 2015).
Under future climate scenarios changes in marine carbonate chemistry are predicted to have deleterious effects on calcification in bivalves and their ecological functions (Gazeau et al., 2013;Waldbusser et al., 2015).
In the Belgian coastal zone, Desmit et al. (2020) reported a peak of pH values between 1975 and 1988, followed by an overall longterm stability from 1989 to 2015 and progressively more marked seasonal changes over time, that corroborates the observed shift of seasonal Chl-a concentration peaks (Figure 3a; Figure S8). Reported pH trends are in line with the strong impact of altered phytoplankton communities on coastal carbon cycling and sequestration (Gypens et al., 2009), as well as eutrophication and increased riverine nutrient inputs (increased alkalinity) buffering ocean acidification in surface coastal waters (Borges & Gypens, 2010). This pattern supports how the coastal zones will not experience acidification gradually, as in the open oceans, but rather as increased frequency, duration, and magnitude of unfavourable events for coastal ecologically and economically valuable species (Waldbusser et al., 2015). However, reported pH changes may reflect nonlinear alterations in carbonate parameters (i.e. pCO 2 , saturation state, DIC) that cannot be easily reconstructed for our study location. This covariation of carbonate system variables indicates the importance of complete carbonate chemistry monitoring and the potential limitations of historical datasets, which are mainly focused on pH records, in linking calcifiers' biology to water chemistry observations. Historical alterations in other local processes due to the rapid urban development of the Belgian coast (Speybroeck et al., 2008), which could not be directly included in our analysis, such as changes in riverine input and sediment transport, might have influenced small-scale coastal gradients with potential (in)direct effects on mussel beds. Most freshwater inputs are characterized by lower pH than oceanic waters, mainly due to low alkalinity and high pCO 2 originating from higher loads of organic matter (Desmit et al., 2020).
Historical trends of riverine inputs in the southern North Sea regions indicate decreasing volumes of freshwater discharge after the 1980s (van Beusekom et al., 2019). Reductions of river runoffs might lead to alterations of cross-shore SSS gradients across short distances (higher salinity), increasing the saturation state of calcium carbonate (Thomsen et al., 2015), with the potential formation of localized areas along the Belgian coast with more suitable environment for mussel calcification (Telesca et al., 2019). Moreover, long-term increases in dune volume have been reported along most of the Belgian coastline after the 1980s (Strypsteen et al., 2019). However, changes in associated sediment volume and transport might have had only indirect implications for shell biomineralization through physical influences on mussel beds dynamics.
Mytilus edulis fitness is tightly linked to shell integrity, deposition, and composition (Telesca et al., 2019), and to its capacity to change shell production in response to the presence or absence of predators (Freeman, 2007;Lowen et al., 2013;Sherker et al., 2017). The organic periostracum provides protection against predatory and endolithic borers as well as corrosive waters (Harper, 1997;Harper & Skelton, 1993). Increased deposition of calcareous layers improves defences in thicker-shelled individuals against durophagous and drilling predators (Freeman, 2007;Sherker et al., 2017). Given the differential defensive value of shell layers, differences in energetic costs of CaCO 3 (low) versus organics (high) production (Watson et al., 2017), as well as marked environmental and ecological alterations, mussels will be under strong selection pressure favouring the capacity to modify shell production to cope with a range of disturbances.
Although reported environmental changes could have a positive but relatively minor contribution on Mytilus calcification, temporal changes in predation pressure and strategy were the key drivers of the observed historical increase in shell deposition (Figures 1c and 5).
Shell thickening and increased calcification are in line with increasing decapod abundances and the inducible defence responses to crab predation (Freeman, 2007;Lowen et al., 2013;Figure 5). Prismatic layer thickness has been documented as the factor largely determining M. edulis shell vulnerability to seabird predation (Le Rossignol et al., 2011). This suggests the rapid prismatic layer thickening after 1995 ( Figure 5) likely reflects a response to the exponential growth of seagull colonies (Figures 3b and 6). Similarly, the rapid periostracal thinning (−29%) after 1984 and less organic matrix (−13%) in modern shells (Figures 2a and 5) suggest compensatory responses to the disappearance of N. lapillus (1981) and to buffered ocean acidification effects (Borges & Gypens, 2010), both selecting for shells with higher resistance to dissolution (i.e. thick periostracum and abundant organic matrix; Harper & Skelton, 1993;Telesca et al., 2019;Thomsen et al., 2017). Although past shell patterns cannot be fully explained due to the limited historical records, our findings support the potential of mussels for adjustments in the deposition of calcareous and organic shell components as a compensatory response to predation strategy-specific pressures.
Our results highlight how complex, local-scale changes in ecological conditions can trigger compensatory mechanisms and lead to unanticipated biological outcomes that run contrary to global-scale predictions. In this regional scenario (Figure 6), the disappearance of a keystone predator, selecting for energetically 'expensive' responses (periostracum-and organic-enriched shells), increased food supply via eutrophication, and buffered ocean acidification, likely led to compositional trade-offs with formation of 'cheaper', calcareous-dominated (decreased organic content) shells as a response to increased durophagous predation.

| Archival collections as key tools for environmental change research
The use of archival specimens collected regularly from geographically well-delimited locations over the last century  allowed the characterization of overall temporal responses, incorporating combined impacts of long-term acclimation and adaption, in a calcifying foundation species within a system that closely mimics intertidal hard rock environments ( Figure 6).
Our findings demonstrate a surprising capacity of Mytilus to withstand a wide range of environmental perturbations and altered predation impacts through unexpected levels of qualitative and quantitative compensatory adjustments in shell biomineralization.
Increasing evidence supports the potential for resistance in calcifiers to environmental changes (Cross et al., 2018(Cross et al., , 2019Telesca et al., 2019;Thomsen et al., 2017), suggesting compensatory biomineralization to disturbance may act in a range of organisms, being driven by species-specific skeletal structures, physiological plasticity, and genetic adaption. This highlights its potential value in calcifiers to maintain their ecological functions in rapidly changing environments via formation of skeletal structures reflecting the interactions between climatic, ecological, and anthropogenic selective forces (Figure 6), and the cost of compensation.
Our results also highlight the potential role of local ecological changes of predator communities in stimulating unexpected compensatory responses of biomineralization in adult mussels that we often assume will be less accentuated than in early life-history stages and follow environmental change according to experimental predictive models. Altered species interactions represent a key leverage point through which compensatory mechanisms could indirectly attenuate the direct impacts of rapidly changing climates on natural populations. As climate change accelerates, local anthropogenic impacts may determine future shifts in the balance of environmental and ecological forces on calcifiers and lead to unanticipated biological patterns and ecosystem alterations.
Compensatory biomineralization represents a potentially important, but still poorly explored mechanism for resistance in calcifiers to ecological and environmental change complementing other physiological, population, and community-level adjustments (Kroeker et al., 2017;Peck et al., 2015). A better understanding of the mechanisms driving biological variation is needed to inform more realistic projections of organismal responses to climate and anthropogenic impacts (Vargas et al., 2017), and to guide more effective environmental management policies (Urban et al., 2016). In this regard, historical studies of archival records represent a rarely used but powerful approach that, as here, can identify unexpected responses to change and also complement experimental approaches to provide a more integrative assessment of the direct and indirect pathways of ecological change or stability.

ACK N OWLED G EM ENTS
We thank the curator Yves Samyn (RBINS) for making available archival collections. We are grateful to Guy Hillyard (British Antarctic Survey) for assistance with microscopy. We also thank the Management Unit of