Blue mussel shell shape plasticity and natural environments: a quantitative approach

Shape variability represents an important direct response of organisms to selective environments. Here, we use a combination of geometric morphometrics and generalised additive mixed models (GAMMs) to identify spatial patterns of natural shell shape variation in the North Atlantic and Arctic blue mussels, Mytilus edulis and M. trossulus, with environmental gradients of temperature, salinity and food availability across 3980 km of coastlines. New statistical methods and multiple study systems at various geographical scales allowed the uncoupling of the developmental and genetic contributions to shell shape and made it possible to identify general relationships between blue mussel shape variation and environment that are independent of age and species influences. We find salinity had the strongest effect on the latitudinal patterns of Mytilus shape, producing shells that were more elongated, narrower and with more parallel dorsoventral margins at lower salinities. Temperature and food supply, however, were the main drivers of mussel shape heterogeneity. Our findings revealed similar shell shape responses in Mytilus to less favourable environmental conditions across the different geographical scales analysed. Our results show how shell shape plasticity represents a powerful indicator to understand the alterations of blue mussel communities in rapidly changing environments.

Step 8: calculation of harmonic coefficients (4 × harmonic × outline) Step 9a: multivariate analysis (PCA) Step 4: outline normalisation (smoothing, centreing, scaling) Step 6: Procrustes superimpostionon and normalisation of starting points Step 1: digitization of orthogonal shell views Step 2: alignment and conversion to black mask Step 5: pseudo-landmark configurations (1000) X Y Step 7: EFA and calibration Step 10a: multivariate analysis (MANOVA, LDA) Step 10b: modelling (GAMMs) Step 9b: mean shape and TPS analysis Figure S1. Blue mussel shell morphology and shape analysis.  Figure S2. Variation in outlines and shape features (Atlantic system). (a) Scatterplots of the first two PCs from a PCA performed on the elliptic Fourier coefficients, of lateral and ventral shell views, showing a clear separation and marked shape variation among specimens from the pooled locations. Confidence intervals for each group and the reconstructed morphospace (background) were represented. (b) Mean shape differences of lateral (VL) and ventral (VV) views between populations at the extremes of the morphospace. Population A showed rounder and narrower shells, with bigger height and more convex ventral sides than population 9 with elongated, curved and wide shells.
. Mean shape differences at the extremes of the morphospace. Deformation grids (left), depicting the bindings required to pass from an extreme of the morphospace to another, and iso-deformation lines (right), representing the outline regions subjected to different degrees of change (blue: low deformation; red: strong deformation), for (a) System 1, populations 8 and 3, (b) System 2, populations A and E, (c) System 3, batches I and VII, and (d) Atlantic system, populations A and 9. Figure S4. PCs contribution to shape reconstruction (System 1). Contribution of the first five shape variables (PCs) to shape variation (large-scale study system). The average shell shapes, for both lateral and ventral views, were represented for increasing values along each PC (−3 SD, Mean, +3 SD) and extreme shapes were compared (±3 SD). Figure S5. PCs contribution to shape reconstruction (System 2). Contribution of the first five shape variables (PCs) to shape variation (medium-scale study system). The average shell shapes, for both lateral and ventral views, were represented for increasing values along each PC (−3 SD, Mean, +3 SD) and extreme shapes were compared (±3 SD).  Figure S6. PCs contribution to shape reconstruction (System 3). Contribution of the first five shape variables (PCs) to shape variation (small-scale study system). The average shell shapes, for both lateral and ventral views, were represented for increasing values along each PC (−3 SD, Mean, +3 SD) and extreme shapes were compared (±3 SD).  Figure S7. PCs contribution to shape reconstruction (Atlantic system). Contribution of the first five shape variables (PCs) to shape variation (pooled locations). The average shell shapes, for both lateral and ventral views, were represented for increasing values along each PC (−3 SD, Mean, +3 SD) and extreme shapes were compared (±3 SD). (Significance, n.s. p > 0.01, * p < 0.01, * * p < 0.001, * * * p < 0.0001)

10/21
System Location ID n Longitude Latitude Status of mussels

11/21
System Selected final model  Table S3. MANOVA output. Summary of MANOVA models on the first 10 calculated shape variables (PCs) for each study system, showing the effects of location of origin and shell length (size) on shape variance.

12/21
System PC %Variance Contribution to the shell shape System 1 PC1 38.7 Variation in shell height, ligament angle, ventral margin shape and shell width: high values corresponded to elongated and narrow shells, with acute ligament angles and convex ventral margins; low values to round and wide shells, with big ligament angles and concave ventral margins.
PC2 30.0 Variation in the shape of ventral margin, ligament angle and shell height, with small contribution to the symmetry of ventral profile: low values corresponded to curved shells, with concave ventral margins, small ligament angles and symmetric ventral profiles; high values corresponded to round shells, with big ligament angles, convex ventral margins and less symmetric ventral profiles.
PC3 10.7 Variation in the shape of umbo, ligament and ventral margin: low values corresponded to "curved" shells with concave ventral margins and ligaments, and an umbo oriented towards the ventral side; high values described elliptical shells with convex ventral margins and ligaments, and an umbo oriented toward the anterior side.
PC4 5.6 Variation in the shape of ventral margin and differences between elliptical and "curved" shells, with a concave ventral margin.    Table S4. PCs contribution to the shell shape. Proportion of shape variance captured by individual shape variables and description of their contributions to the shell features and mean shape reconstruction for each study system (Supplementary Figure S5, S6, S7).

13/21
Supplementary Document S1 Information on the environmental parameters investigated (water temperature, salinity and chl-a concentration) for System 1, 2, and 3. Details on datasets used for measurement/calculation of mean annual values are provided.

System 1 and System 2 -Modelled datasets
List of the datasets used for the calculation of mean annual values of environmental descriptors. This study has been conducted using the Copernicus Marine Service Products: COPERNICUS -Marine Environment Monitoring System (http://marine.copernicus.eu/).

DATASET #1
Product name BALTIC SEA PHYSICS ANALYSIS AND FORECAST

Short description
This Baltic Sea physical model product provides forecasts for the physical conditions in the Baltic Sea. The Baltic forecast is updated twice daily providing a new two days forecast with hourly data for sea level variations, ice concentration and thickness at the surface, and temperature, salinity and horizontal velocities for the 3D field. The product is based on the 3D ocean model code HBM developed within the Baltic ocean community.

Short description
This Baltic Sea biogeochemical model product provides forecasts for the biogeochemical conditions in the Baltic Sea. The Baltic forecast is updated twice daily providing a new two days forecast with hourly data for the parameters dissolved oxygen, nitrate, phosphate, chl-a. The product is produced by the biogeochemical model ERGOM one way coupled to the Baltic 3D ocean model HBM, which provides the CMEMS Baltic physical ocean forecast product.

Spatial resolution 2 km
Vertical coverage from -400 m to 0 m

Update frequency Daily
Production unit BAL-DMI-COPENHAGEN-DK

Short description
The Operational Mercator global Ocean analysis and forecast system at 1/12 degree is providing 7 days of 3D global ocean forecasts updated daily and ocean analysis updated weekly. The time series is aggregated in time in order to reach a two full year's time series sliding window. This product includes daily mean files of temperature, salinity, currents, sea level and ice parameters from the top to the bottom of the Ocean over the Global Ocean. It also includes 2-hourly mean surface fields for temperature, currents and sea level.

Update frequency Daily
Production unit NWS-METOFFICE-EXETER-UK

System 3 -Habitat monitoring
Variation in key environmental parameters was measured fortnightly (May 2015 -June 2016) at a traditional longline mussel farm "Glencoe Shellfish" in Loch Leven (56.7097 N, -5.0292 W; Scotland, UK). Temperature ( • C) and salinity (psu), from vertical CastAway CTD profiles, and chl-a concentration (mg m −3 ), using fluorometry, were measured at one, three, five, and nine meters depth along the rope (Michalek manuscript in preparation).

Supplementary Methods
Mytilus shell shape was analysed through an elliptic Fourier analysis (EFA) of outlines [15][16][17] . This geometric morphometrics approach 18,19 was performed in order to examine shell shape variation both within and between groups of individuals. As other modern morphometrics approaches, it considers outlines as a whole, taking into account all the geometrical relationships of the input data. EFA is a powerful method to extract this geometric information and has been implemented on the concept of Fourier series: to decompose a periodic function into a sum of more simple trigonometric functions, such as sine and cosine 17,20 . These simple functions have frequencies that are integer multiples, therefore they are harmonics of one another.
This approach fits Fourier series separately on the x and y coordinates of an outline, projected on the Cartesian plane, as a function of the curvilinear abscissa 16,17,20 . EFA is then used to extract the geometrical information from outlines, described as periodic functions 16,21 , through their decomposition into the harmonic sum of trigonometric functions, called harmonics. Low-frequency harmonics approximate coarse-scale trends in the original outlines, while high-frequency harmonics fit their fine-scale variations 17 . These can be normalized to remove homothetic, translational or rotational differences between shapes and smoothed in order to remove outline noise during the digitization process 22 . Harmonic coefficients are then extracted and used as shape variables capturing shape information. The geometrical information contained in the outlines is thus quantified and can be analysed with classical multivariate tools (i.e. multivariate analysis of variance, principal component analysis and linear discriminant analysis).
EFA of outlines allows shape reconstruction from the numerical signature and this improved method has great advantages compared to more traditional approaches 17,19,22 : complex shapes can be fitted, outlines smoothed, starting points and coefficients can be normalised to remove homothetic, translational and rotational differences between outlines 17,19,[21][22][23] . Shape analyses were carried out using the Momocs ("Morphometrics using R") 17 package with the R v3.3.0 software 24 .
EFA of outlines: acquisition, processing and analysis (Supplementary Figure S1) • Digital images of orthogonal lateral and ventral shell views (left valves) were acquired with a high-resolution digital camera (Nikon D3300 camera, fitted with Sigma 105mm f/28 EX DG Macro lens); • Photographs were processes with an image analysis software ( c Adobe Photoshop), centred and consistently aligned; • Photographs were converted into black masks on a white background (greyscale, 8-bit) and only the shapes of intact shells were retained; • Outlines were isolated, converted into a list of (x; y) pixel coordinates and used as input data.
Lateral and ventral views of each shell were processed independently and later combined for analysis: • Prior to calculation of elliptic Fourier transforms, an outline alignment through geometric operations was directly performed on the list of coordinates. This a priori normalisation was required to avoid potential bias introduced by the numerical adjustment of shapes prone to bad alignment (usually circular or with bilateral symmetry). Indeed, for mussel outlines, "consuming" their first harmonic in order to normalise higher rank harmonics 20 , determined a poor numerical alignment, resulting in not homologues elliptic Fourier descriptors 21,22 ; • Outlines were first smoothed to remove any noise introduced during the digitization process, centred and outline coordinates were rescaled by their centroid size; • Equal number of points were sampled along each outline (1000 pseudo-landmarks); • Point configurations were aligned through a Procrustes superimposition 20, 25 and starting points normalised; • An EFA was then computed on the resulting coordinates from shapes invariant to outline size, rotation and position; • After preliminary calibration, through inspection of i) the outline reconstruction efficiency, ii) the deviation from the optimal fit and iii) the spectrum of harmonic Fourier power, seven harmonics were chosen to encompass 98% of the total harmonic power 17, 23 ; • Four coefficients per harmonic (28 descriptors) were extracted for each outline and used as variables quantifying the geometrical information 20, 21 .

18/21
The shape information contained in the outlines was then quantified and analysed with classical multivariate tools: • Principal component analysis (PCA), with a singular value decomposition method, was performed on the matrix of coefficients, without rescaling, in order to define axes capturing the most of shape variation among individuals 26,27 ; • The first 10 PCs accounting for the 97% of outline variability were used as new shape variables; • Multivariate Analysis of Variance (MANOVA) was performed on the new shape variables to test for a significant effect of location of origin and shell size on shape variances 20,26,27 ; • Linear discriminant analysis (LDA) was performed on shape variables, with a leave-one-out cross-validation procedure, in order to evaluate whether the linear combination of shape features (discriminant function) was able to discriminate between Mytilus species. A priori classification probabilities were set to be proportional to group sizes and Wilks' λ was calculated to test for significant discriminations. Discriminant coefficients were estimated in order to identify shell shape features that optimized the between-species differences "relative" to the within-species variation 20, 26 ; • Shape differences at the extremes of the morphospace were visualised with representation of mean shapes, deformation grids 28 and iso-deformation lines, which were obtained though mathematical formalisation of thin plate splines (TPS) analysis 25 ; • The first five shape variables (PCs), describing distinguishable shell features along the mussel outlines, were then analysed with generalized additive mixed models (GAMMs) in order to identify relationships between mussel shape variation and environmental gradients.