Main drivers of fecundity variability of mussels along a latitudinal gradient: Lessons to apply for future climate change scenarios

Journal of Marine Science and Engineering Article Main Drivers of Fecundity Variability of Mussels along a Latitudinal Gradient: Lessons to Apply for Future Climate Change Scenarios Gabriela F. Oliveira 1,* , Hanifah Siregar 2, Henrique Queiroga 1and Laura G. Peteiro 3,*   Citation: Oliveira, G.F.; Siregar, H.; Queiroga, H.; Peteiro, L.G. Main Drivers of Fecundity Variability of Mussels along a Latitudinal Gradient: Lessons to Apply for Future Climate Change Scenarios. J. Mar. Sci. Eng. 2021,9, 759. https://doi.org/ 10.3390/jmse9070759 Academic Editor: Dariusz Kucharczyk Received: 16 June 2021 Accepted: 5 July 2021 Published: 9 July 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1CESAM—Centre for Environmental and Marine Studies, Department of Biology, University of Aveiro, 3810-193 Aveiro, Portugal; [email protected] 2Wildlife Conservation Society, Bogor 16128, Indonesia; [email protected] 3Instituto de Investigaciones Marinas (IIM), Consejo Superior de Investigaciones Científicas (CSIC), 36208 Vigo, Spain *Correspondence: [email protected] (G.F.O.); [email protected] (L.G.P.) Abstract: Bivalve relevance for ecosystem functioning and human food security emphasize the importance of predictions of mussel performance under different climate stressors. Here, we address the effect of a latitudinal gradient of temperature and food availability on the fecundity of the Mediterranean mussel to try to better parameterize environmental forcing over reproductive output. We show that temperature plays a major role, acting as a switching on–off mechanism for gametogenesis, while food availability has a lower influence but also modulates the number of gametes produced. Temperature and food availability also show different effects over fecundity depending on the temporal scale evaluated. Our results support the view that the gametogenesis responds non-linearly with temperature and chlorophyll concentration, an issue that is largely overlooked in growth, production and energy budgets of bivalve populations, leading to predictive models that can overestimate the capability of the mussel’s populations to deal with climate change future scenarios. Keywords: reproductive output; latitudinal gradients; temperature; food availability; Mytilus galloprovincialis; climate change; bivalves 1. Introduction The accumulation of evidence on climate change impact on marine ecosystems increasingly points towards complex cascading effects at different levels of biological organization [ 1 – 5 ]. Rising atmospheric CO 2 has consequences not only for ocean temperature and acidification, but also on circulation patterns, stratification, or oxygen content [ 3 , 4 ]. Physical and chemical changes derived from anthropogenic activity have direct and indirect effects on physiology, behaviour and interactions of marine organisms, with consequences at population and community levels, compromising ecosystem functioning and finally its ability to provide a range of goods and services to humans [2–5] Understanding and foreseeing the response of marine populations to climate change is crucial for human planning and development of mitigation and adaptation measures [3,6,7] . A commonly employed approach is to couple population and climate models to test for the effect of different climate scenarios on species distribution, performance or productivity [3,7–9] . Many of these studies agree to identify reef forming calcified organisms (corals, mussels, oysters, etc.) as one of the most vulnerable, because of their particular sensitivity to ocean acidification and other cascading effects related to the rising temperatures [1,3,10]. Mussels are amongst the most studied reef forming organisms, not just because of their interest as ecological engineers, driving biodiversity patterns [ 11 ] or regulating benthicpelagic coupling process and water quality [ 12 , 13 ], but also because of their relevance for aquaculture and food security [ 14 ]. Recommendations of scientific advisory committees J. Mar. Sci. Eng. 2021,9, 759. https://doi.org/10.3390/jmse9070759 https://www.mdpi.com/journal/jmse J. Mar. Sci. Eng. 2021,9, 759 2 of 17 for a sustainable response to the increasing human demand of animal protein include an increment of bivalve production by 100 Mt before 2030 [ 15 ]. The large vulnerability of bivalves coupled to their relevance for food security emphasizes the importance of accurate predictions of mussel performance under different climate change scenarios, and has prompted the use of this species as a model for the study of the impact of several stressors [16–19]. Reproductive success is a key process controlling population dynamics, and identifying the environmental drivers that control reproductive output is key to be able to forecast the consequences of climate change [ 7 , 20 , 21 ]. The reproductive cycle of bivalves and other invertebrates comprises a sequence of biological events comprising gametogenesis, spawning and a resting period before gonad restoration [ 22 ]. This cycle is controlled by endogenous and environmental factors, and their interactions influence the onset and the duration of the reproductive cycle [ 22 ]. Among these environmental factors, temperature and food availability are considered to be the main factors regulating reproduction in marine invertebrates [ 20 – 24 ]. Temperature is usually regarded as a trigger for spawning, but also regulates the rate of gonadal development with a cumulative effect [ 20 – 24 ]. Food availability has been directly associated to fecundity (number of gametes produced) [ 24 , 25 ] although some authors also pointed to phytoplankton blooms as a spawning cue, which allows matching larval development with optimal environmental conditions [26]. Under unfavourable environmental conditions, bivalves activate costly defence and repair mechanisms, which reallocate energy away from reproduction towards somatic maintenance [ 16 , 20 ]. Therefore, reabsorption of gametes and fecundity loss are expected consequences of stressing scenarios which will become more persistent in the future [1,3] . Nonetheless, studies forecasting reproductive output and performance of bivalve populations, based on the Dynamic Energy Budget (DEB) principle, consistently predict a maintenance or a slight increase of reproductive output for bivalves in all the locations where warming future scenarios don’t overpass thermal tolerance of the species [ 7 , 9 , 20 , 23 , 27 – 30 ]. DEB and other models assume that physiological rates are proportional to size and follow at least a temperature-dependent relationship according to the optimal physiological range for somatic maintenance of the species, but do not usually consider specific environmental limits for gamete development [ 31 – 33 ]. Even when many attempts have been made to increase the complexity of the description of energy allocation to reproduction to make it more realistic [ 23 ], specific non-linear responses of gamete production to environmental conditions are not taken into account. Optimum values of temperature and food availability might differ between somatic and gametogenic processes. For example, while temperature tolerance ranges of Mytilus galloprovincialis range from 5 ◦ C to 35 ◦ C, with an optimum for somatic maintenance around 17.5 ◦ C [ 28 ], some studies have determined a reproduction limiting temperature of 19 ◦ C for the same species, with an optimum for the production of vitellogenic oocytes around 14 ◦C [34]. A correct parameterization of the complex relationships between environmental drivers and gamete production requires a better understanding of the underlying processes. Knowledge about the scales of spatial and temporal variability in reproductive output coupled to environmental variability is essential to gain a better understanding of the reproductive physiology, and to be able to incorporate it properly in the prediction of future scenarios. In this study, we examined the variability in reproductive outputs of a series of populations of Mytilus galloprovincialis along a latitudinal gradient covering the Portuguese coast during the main reproductive season (Spring) for two different years. Consistency of latitudinal patterns of fecundity between years along with environmental variability was explored to identify relationships along an environmental gradient. Our goal was to identify non-linear relationships between gamete production and environmental variables at different temporal scales, to understand their influence when looking at the complete gametogenesis process ( ≈ 4 months) or to the last phases of gamete maturation. That J. Mar. Sci. Eng. 2021,9, 759 3 of 17 information might help to provide a better parameterization of reproduction, and to develop more accurate models for mussel populations forecasting. 2. Materials and Methods 2.1. Sampling Design Mediterranean mussels, Mytilus galloprovincialis, were collected at 7 sites along the Portuguese coastline in 2014 and 2017 (Figure 1) following a latitudinal gradient from North to South: Póvoa do Varzim, Costa Nova, Peniche, Cascais, Galapos and Zavial. Sampling locations ranged from around 37 ◦ N to 41 ◦ N (distance ≈ 560 km, mesoscale), and took place during spring in 2014 and 2017 (14 April to 4 May 2014 and 18 May to 22 June 2017 ). At each site, 150 individuals were collected randomly along the intertidal, for length frequency and fecundity analysis. J. Mar. Sci. Eng. 2021, 9, 759 3 of 18 was explored to identify relationships along an environmental gradient. Our goal was to identify non-linear relationships between gamete production and environmental variables at different temporal scales, to understand their influence when looking at the complete gametogenesis process (≈4 months) or to the last phases of gamete maturation. That information might help to provide a better parameterization of reproduction, and to develop more accurate models for mussel populations forecasting. 2. Materials and Methods 2.1. Sampling Design Mediterranean mussels, Mytilus galloprovincialis, were collected at 7 sites along the Portuguese coastline in 2014 and 2017 (Figure 1) following a latitudinal gradient from North to South: Póvoa do Varzim, Costa Nova, Peniche, Cascais, Galapos and Zavial. Sampling locations ranged from around 37° N to 41° N (distance ≈ 560 km, mesoscale), and took place during spring in 2014 and 2017 (14 April to 4 May 2014 and 18 May to 22 June 2017). At each site, 150 individuals were collected randomly along the intertidal, for length frequency and fecundity analysis. Figure 1. Fecundity sampling locations along the Portuguese coastline. 2.2. Environmental Variables Figure 1. Fecundity sampling locations along the Portuguese coastline. 2.2. Environmental Variables Sea Surface Temperature (SST), Sea Surface Salinity (SSS) and Chlorophyll-a concentration (Chl-a) as environmental variables were obtained from the beginning of the year to the sampling date at each location. Daily averaged values of SST ( ◦ C) and SSS of 2014 were provided by Hycom model (http://www.hycom.org) and values of Chl-a concentration (mg/m 3 ) were obtained from interpolation of satellite data (MODIS, SeaWiFS and J. Mar. Sci. Eng. 2021,9, 759 4 of 17 MERIS) supplied by CERSAT/IFREMER (http://ftp.ifremer.fr/ifremer/cersat/products/ gridded/ocean-color/atlantic/EURL4-CHL-ATL-v01) (both accessed on 4 July 2014). The SST and SSS values of 2017 were obtained from the Mercator Ocean (Daily Global Analysis and Forecast of Ocean Physics, http://www.mercator-ocean.fr, EU Copernicus Marine Environment Monitoring Service, accessed on 7 November 2017), resolution 1/12 ◦ ( ≈ 7 km for medium latitudes). Values of Chl-a concentration (mg/m 3 ) were obtained from CMEMS Products (Copernicus-GlobColour Project, ACRI-ST, EU Copernicus Marine Environment Monitoring Service, accessed on 7 November 2017), which consist of daily data interpolated and reprocessed (multi-year time series) from satellite observations (SeaWIFs, MODIS-Aqua, MERIS and VIIRSN), and have a horizontal resolution of 1 km. Based on these data, a set of variables regarding temperature and food availability were created to identify different effects of environmental variability, at long and short-time scales, on fecundity. Concerning short-term effects, monthly averaged SST (average SST, ◦ C) and Chl-a concentration (average Chl-a, mg/m 3 ) during the month preceding the sampling date were calculated for each location. Data from 2014 for Cascais, Sines and Zavial were unavailable, so these were not included in the analysis. Looking at the long-term effects, we calculated the averaged Chl-a concentration during 10 days around the peak of the spring transition (Max Chl-a) and the number of days with SST over 14 ◦ C (T > 14 ◦ C) during the 4 months previous to the sampling date for each location to match the expected gametogenesis duration [ 35 , 36 ]. The threshold of 14 ◦ C was selected based on the temperature optimum reported in laboratory experiments for the production of vitellogenic oocytes [34]. 2.3. Fecundity Analysis All mussels collected were measured along their anterior–posterior axis to calculate their total length (L, mm), using a digital calliper. Samples were divided into size classes with 10 mm intervals, from 5 mm to 95 mm (no individuals larger than 100 mm were detected in the intertidal). From each size class, 5 females were randomly selected for fecundity analysis. Females were first identified visually based on the orange coloration of their gonads and the presence of oocytes confirmed under a stereomicroscope. From the gonads of each female, a 1 cm diameter sample was taken to assess fertility, according to the methodology presented in Sukhotin and Flyachinskaya [ 37 ]. Each gonad sample was weighed, gently macerated and diluted with 20–25 mL of salted water. From this dilution, 3 sub-samples between 150–250 µ L were taken to count the oocytes under a microscope, with the help of a Bogorov counting chamber. Another gonad sample (1 cm) was taken to estimate the weight of the tissue used for oocytes quantification. The rest of the gonad was separated from the remaining soft tissues to estimate the weight of gonads with regard to the rest of tissues. After dissecting, all samples were dried at 60 ◦ C for 24 h to obtain dry weights. Absolute fecundity (AF, number of oocytes/female) and relative fecundity (RF, number of oocytes/g gonad) were estimated according to the following formulas: AF =NWg Ws Vd Vs 100 95 RF = AF Wg where N is the average number of oocytes observed in 3 subsamples, Wg is the dry weight of the total gonad (g), Ws is the dry weight of the gonad sample extracted for the oocyte count (g), Vd is the volume of water where the oocytes were diluted (mL), and Vs is the volume used in the subsamples for counting oocytes (mL). Since the extraction of oocytes from the samples is only 95% efficient, with 5% of oocytes remaining in the tissues during the extraction process [37], a 100/95 correction factor was applied to the estimate of AF. 2.4. Statistical Analysis In order to evaluate the temporal and spatial variability of the AF and RF we used factorial ANOVAs with location and year as fixed factors. The relationship between size J. Mar. Sci. Eng. 2021,9, 759 5 of 17 and AF was analysed with ANCOVA (Separate Slopes model), separately for each year, including size (L, mm) as a covariable and location (Loc, 7 levels) as a fixed factor. In these analysis AF, RF and Size were log-transformed. The data and residues were evaluated for normality (Q-Q plot and Shapiro-Wilk Normality Test) and homoscedasticity (Levene’s Test). ANOVAs and ANCOVAs were performed using the software Statistica 12 (Stat Soft, Tulsa, OK, USA). Generalized additive models (GAM), as implemented in the mgcv library of R 3.6.2 [ 38 ], were used to investigate the seasonal patterns (day of the year) and local variation (location) of the environmental variables measured (SST, SSS and Chl-a). GAM models were also employed to investigate the relationships between the environmental parameters (days > 14.5, average SST, average SSS, max Chl-a) on RF (log-transformed). The Akaike information criterion (AIC) was used to select the optimal set of variables for inclusion in the model. Model validation included the verification of homogeneity, normality and independence assumptions [39]. 3. Results 3.1. Environmental Variables Seasonality largely explained the variability observed in temperature and salinity (77.7 to 86.5% of deviance explained by the variables location and day of the year; Table 1 , Figure 2 ). Latitudinal patterns were consistent during both sampling years, with a temperature and salinity increase from North to South (Table 1). Nonetheless, the estimated coefficients revealed warmer temperatures and lower salinities in 2017 (Table 1). J. Mar. Sci. Eng. 2021, 9, 759 5 of 18 during the extraction process [37], a 100/95 correction factor was applied to the estimate of AF. 2.4. Statistical Analysis In order to evaluate the temporal and spatial variability of the AF and RF we used factorial ANOVAs with location and year as fixed factors. The relationship between size and AF was analysed with ANCOVA (Separate Slopes model), separately for each year, including size (L, mm) as a covariable and location (Loc, 7 levels) as a fixed factor. In these analysis AF, RF and Size were log-transformed. The data and residues were evaluated for normality (Q-Q plot and Shapiro-Wilk Normality Test) and homoscedasticity (Levene’s Test). ANOVAs and ANCOVAs were performed using the software Statistica 12 (Stat Soft, Tulsa, USA). Generalized additive models (GAM), as implemented in the mgcv library of R 3.6.2 [38], were used to investigate the seasonal patterns (day of the year) and local variation (location) of the environmental variables measured (SST, SSS and Chl-a). GAM models were also employed to investigate the relationships between the environmental parameters (days > 14.5, average SST, average SSS, max Chl-a) on RF (log-transformed). The Akaike information criterion (AIC) was used to select the optimal set of variables for inclusion in the model. Model validation included the verification of homogeneity, normality and independence assumptions [39]. 3. Results 3.1. Environmental Variables Seasonality largely explained the variability observed in temperature and salinity (77.7 to 86.5% of deviance explained by the variables location and day of the year; Table 1, Figure 2). Latitudinal patterns were consistent during both sampling years, with a temperature and salinity increase from North to South (Table 1). Nonetheless, the estimated coefficients revealed warmer temperatures and lower salinities in 2017 (Table 1). Figure 2. Estimated smoothing curves of day of the year for Sea Surface Temperature (SST, °C), Sea Surface Salinity (SSS) and Chlorophyll-a concentration (Chl-a, mg/m3) from 2014 and 2017, expressed as a mean between locations along the Portuguese coastline. Results of the General Additive Model showing the effect of the variable “day of the year” (Julian year). Dashed lines show a 95% Confidence Interval and tick marks along the x-axis represent when observations occurred. Vertical dashed line (2017) is a reference to the last day from 2014 (day 120). Figure 2. Estimated smoothing curves of day of the year for Sea Surface Temperature (SST, ◦ C), Sea Surface Salinity (SSS) and Chlorophyll-a concentration (Chl-a, mg/m 3 ) from 2014 and 2017, expressed as a mean between locations along the Portuguese coastline. Results of the General Additive Model showing the effect of the variable “day of the year” (Julian year). Dashed lines show a 95% Confidence Interval and tick marks along the x-axis represent when observations occurred. Vertical dashed line (2017) is a reference to the last day from 2014 (day 120). J. Mar. Sci. Eng. 2021,9, 759 6 of 17 Table 1. Structure of the General Additive Model describing Sea Surface Temperature ( ◦ C), Sea Surface Salinity (SSS) and Chlorophyll-a concentration (mg/m 3 ) variability along the coastline from January to May 2014 and January to June 2017. S.E.: standard error; e.d.f.: estimated degrees of freedom. SST (◦C) SSS Chl-a (mg/m3) 2014 Parametric Coefficients Parametric Coefficients Parametric Coefficients Location Estimate S.E. t p Estimate S.E. t p Estimate S.E. t p Póvoa do Varzim (Intercept) 13.357 0.031 436.581 p< 0.001 35.294 0.017 2051.909 p< 0.001 2.499 0.138 21.488 p< 0.001 Costa Nova 0.287 0.044 − 6.541 p< 0.001 0.103 0.024 − 4.239 p< 0.001 0.476 0.196 − 6.427 p< 0.050 Peniche 0.486 0.044 11.077 p< 0.001 0.611 0.024 25.197 p< 0.001 −1.258 0.196 − 6.427 p< 0.001 Cascais 0.667 0.052 12.883 p< 0.001 0.418 0.029 14.651 p< 0.001 - - - - Galapos 0.730 0.044 16.614 p< 0.001 0.512 0.024 21.103 p< 0.001 −0.698 0.196 − 3.563 p< 0.001 Sines 1.103 0.052 21.310 p< 0.001 0.611 0.029 24.423 p< 0.001 - - - - Zavial 1.571 0.062 25.430 p< 0.001 0.746 0.034 21.903 p< 0.001 - - - - Smooth terms (non parametric) Smooth terms (non parametric) Smooth terms (non parametric) e.d.f F pe.d.f F pe.d.f F p Day of the Year 8.94 302.6 p< 0.001 6.07 34.04 p< 0.001 8.61 16.3 p< 0.001 R2adjusted: 0.861 % Deviance explained: 86.5% R2adjusted: 0.772 % Deviance explained: 77.7% R2adjusted: 0.278 % Deviance explained: 29.6% 2017 Parametric coefficients Parametric coefficients Parametric coefficients Location Estimate S.E. t pEstimate S.E. t pEstimate S.E. t p Póvoa do Varzim (Intercept) 15.051 0.057 276.091 p< 0.001 33.215 0.032 1059.370 p< 0.001 1.718 0.080 24.265 p< 0.001 Costa Nova 0.623 0.080 − 7.757 p< 0.001 1.112 0.046 − 24.270 p< 0.001 0.224 0.113 − 0.198 0.843 Peniche 0.174 0.080 2.168 p< 0.050 0.667 0.046 14.550 p< 0.001 −0.496 0.113 − 4.384 p< 0.001 Cascais 0.569 0.080 7.084 p< 0.001 0.940 0.046 20.520 p< 0.001 −0.180 0.113 − 1.593 0.111 Galapos 0.921 0.080 11.477 p< 0.001 1.404 0.046 30.640 p< 0.001 −0.495 0.113 − 4.373 p< 0.001 Sines 1.149 0.080 14.306 p< 0.001 1.562 0.046 34.090 p< 0.001 −1.132 0.113 − 10.004 p< 0.050 Zavial 1.489 0.080 18.544 p< 0.001 1.759 0.046 38.380 p< 0.001 −1.241 0.113 − 10.970 p< 0.001 Smooth terms (non parametric) Smooth terms (non parametric) Smooth terms (non parametric) e.d.f F pe.d.f F pe.d.f F p Day of the Year 8.94 767.6 p< 0.001 8.88 119.5 p< 0.001 8.518 12.03 p< 0.001 R2adjusted: 0.862 % Deviance explained: 86.3% R2adjusted: 0.846 % Deviance explained: 84.8% R2adjusted: 0.208 % Deviance explained: 21.7% On the other hand, day of the year had a lower effect on the variability of Chl-a, with around 20%–30% of variance explained, much less than that observed for SSS or SST (Table 1; Figure 2). Both years maintained a similar latitudinal pattern for Chl-a, with a decreasing trend from North to South but with a maximum in Costa Nova (Table 1). Nonetheless, during 2014, the Chl-a peak related to spring transition was much more pronounced than during 2017 (Figure 2), and averaged estimated parameters also indicated larger concentrations of Chl-a during 2014 (Table 1). 3.2. Fecundity Both indicators of fecundity, AF and RF, showed a significant interaction between location and year (Table 2). With regard to the absolute fecundity, larger values were observed in 2017 in comparison with 2014, and in spite of the significant interaction between year and location (Table 2), the latitudinal pattern was quite similar during both sampling years (Figure 3A). With the exception of Zavial, which showed an opposite trend between years, the larger number of eggs produced per individuals were consistently detected at Costa Nova Peniche and Cascais (Figure 3A). Average size of the mussel populations also showed a predominance of larger individuals at those locations (Figure 3B). J. Mar. Sci. Eng. 2021,9, 759 7 of 17 Table 2. Results of Factorial Analysis of Variance (ANOVA) of absolute fecundity (AF, log) and relative fecundity (RF, log) with locations and year as fixed factors (d.f.: estimated degrees of freedom, MS: mean squares). AF (log) Effect df MS F p Location 6 3.918 9.950 <0.001 Year 1 15.472 39.280 <0.001 Loc x Year 6 2.173 5.520 <0.001 Error 255 0.394 - - RF (log) Effect df MS F p Location 6 1.490 8.030 <0.001 Year 1 0.108 0.580 0.446 Loc x Year 6 2.590 13.960 <0.001 Error 255 0.186 - - A positive linear relationship between AF and female size was observed for both sampling years, in all locations but Sines in 2014 (Tables 3and 4). In this case, the lack of relationship between both variables might be related to the narrow range of sizes sampled (25.5–35.9 mm) in that particular location and year. Although AF maintained a positive relationship with size in every location, the significant interaction between location and size indicates differences among locations on the relevance of size on determining the amount of oocytes produced (Table 3) and therefore on the value of the fitted slopes (Figure 3C; Table 4). Nonetheless, similar slopes were fitted for each location at different years (with the exception of Sines), suggesting this relationship to be quite stable at each location (Figure 3C; Table 4). Table 3. Results of Analysis of Covariance (ANCOVA) with Separate Slopes of Absolute Fecundity (AF, log) in different locations with size (mm, log) as covariate (d.f.: estimated degrees of freedom, SS: sum of squares, MS: mean squares). ANCOVA with Separate Slopes: AF (log) vs. Location and Size (log) 2014 Effect Df SS MS F p Intercept 1 1.212 1.212 5.197 p< 0.05 Location ×log size 7 17.152 2.450 10.506 p< 0.05 Location 6 1.590 0.265 1.137 0.346 Error 110 25.653 0.233 - - Total 123 59.150 - - - 2017 Effect Df SS MS F p Intercept 1 2.343 2.343 10.624 p< 0.05 Location ×log size 7 28.736 4.105 18.613 p< 0.05 Location 6 4.015 0.669 3.034 p< 0.05 Error 131 28.892 0.221 - - Total 144 81.505 - - - J. Mar. Sci. Eng. 2021,9, 759 8 of 17 J. Mar. Sci. Eng. 2021, 9, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/jmse Figure 3. Average and standard error of ( A ) absolute fecundity (AF, n ◦ . oocytes per individual), ( B ) Size (mm) and ( C ) Coefficients Slopes of Analysis of Covariance (ANCOVA) between AF (logarithm-transformed) and Size (log), significant values (p< 0.05), except to Sines. Data of the different locations along the Portuguese coastline, from North (left) to South (right), during the years 2014 and 2017. J. Mar. Sci. Eng. 2021,9, 759 9 of 17 Table 4. Coefficients of ANCOVA with Separate Slopes analysis of absolute fecundity (AF, log) in different locations with size (mm, log) as covariate (d.f.: estimated degrees of freedom, SS: sum of squares, MS: mean squares). Separate Slopes Coefficients: AF (log) vs. Location and Size (log) 2014 Location Slope S.E. t p Povoa do Varzim 2.268 0.871 2.603 p< 0.05 Costa Nova 2.996 0.759 3.949 p< 0.05 Peniche 3.565 1.270 2.807 p< 0.05 Cascais 5.503 1.284 4.286 p< 0.05 Galapos 3.862 1.034 3.734 p< 0.05 Sines 0.687 2.737 0.251 0.802 Zavial 5.209 1.576 3.305 p< 0.05 2017 Location Slope S.E. t p Povoa do Varzim 3.237 1.220 2.653 p< 0.05 Costa Nova 1.848 0.829 2.229 p< 0.05 Peniche 2.865 1.039 2.757 p< 0.05 Cascais 5.669 0.797 7.117 p< 0.05 Galapos 2.983 0.526 5.667 p< 0.05 Sines 7.792 1.753 4.445 p< 0.05 Zavial 2.612 0.915 2.855 p< 0.05 When looking at the relative fecundity values, larger RF values were reached in 2014 in comparison to 2014, but also a different spatial pattern was observed between years (Figure 4; Table 2). During 2014, RF values were low at the north and central coast locations and increased almost exponentially towards the southernmost populations (Galapos, Sines and Zavial). Conversely, during 2017, RF increased linearly from the northernmost location to the central coast, reaching a maximum at Cascais which was followed by a sharp drop in RF towards the southernmost locations (Figure 4). J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 10 of 18 Figure 4. Average and standard error of relative fecundity (RF, n°. oocytes per gram of gonad) in the different locations along the Portuguese coastline, from North (left) to South (right), during the years of 2014 and 2017. When looking at the relative fecundity values, larger RF values were reached in 2014 in comparison to 2014, but also a different spatial pattern was observed between years (Figure 4; Table 2). During 2014, RF values were low at the north and central coast locations and increased almost exponentially towards the southernmost populations (Galapos, Sines and Zavial). Conversely, during 2017, RF increased linearly from the northernmost location to the central coast, reaching a maximum at Cascais which was followed by a sharp drop in RF towards the southernmost locations (Figure 4). According to the AIC, the four variables used to assess the effect of long-term (T > 14 °C and Max. Chl-a) and short-term (Average SST and Average Chl-a) environmental variability on RF were included in the GAM, explaining 26.6% of the variability observed across years and locations (Tables 5 and 6; Figure 5). Despite max. Chl-a not showing a significant effect on RF, the variable was maintained in the model because of its traction over AIC (Tables 5 and 6; Figure 5D). The number of days with temperatures over 14 °C during the previous 4 months was the variable with the larger influence on RF (lowest pvalue; Table 6). The highest positive effect of T > 14 °C was around 80 days, with a sharp detrimental effect on RF when warm days overpassed 100 days (Figure 5A). Average SST and Average Chl-a during the previous month had a similar influence on RF (p-values; Table 6), but while Average Chl-a showed a positive and linear relationship with RF (Figure 5C), Average SST only seemed to exert an effect when temperatures rose above 16 °C (Figure 5B). 0×10^5 2×10^5 4×10^5 6×10^5 8×10^5 10×10^5 12×10^5 14×10^5 0×10^5 5×10^5 10×10^5 15×10^5 20×10^5 25×10^5 30×10^5 35×10^5 Povoa do Varzim Costa Nova Peniche Cascais Galapos Sines Zavial RF 2017 (Nº Oocytes/g) RF 2014 (Nº Oocytes/g) 2014 2017 Figure 4. Average and standard error of relative fecundity (RF, n ◦ . oocytes per gram of gonad) in the different locations along the Portuguese coastline, from North (left) to South (right), during the years of 2014 and 2017. According to the AIC, the four variables used to assess the effect of long-term ( T > 14 ◦C and Max. Chl-a) and short-term (Average SST and Average Chl-a) environmental variability on RF were included in the GAM, explaining 26.6% of the variability J. Mar. Sci. Eng. 2021,9, 759 16 of 17 12. Froján, M.; Figueiras, F.G.; Zúñiga, D.; Pérez, F.A.; Arbones, B.; Castro, C.G. Influence of Mussel Culture on the Vertical Export of Phytoplankton Carbon in a Coastal Upwelling Embayment (Ría de Vigo, NW Iberia). Estuaries Coasts 2016 ,39, 1449–1462. [CrossRef] 13. Petersen, J.K.; Saurel, C.; Nielsen, P.; Timmermann, K. The use of shellfish for eutrophication control. Aquac. Int. 2015 ,24, 857–878. [CrossRef] 14. Food and Agriculture Organization (FAO). The State of World Fisheries and Aquaculture 2020: Sustainability in Action; The State of World Fisheries and Aquaculture (SOFIA); Food and Agriculture Organization: Rome, Italy, 2020. 15. Science Advice for Policy by European Academies. Food from the Oceans: How Can More Food and Biomass Be Obtained from the Oceans in a Way That Does Not Deprive Future Generations of Their Benefits; SAPEA: Berlin, Germany, 2017. 16. Farcy, É.; Burgeot, T.; Haberkorn, H.; Auffret, M.; Lagadic, L.; Allenou, J.-P.; Budzinski, H.; Mazzella, N.; Pete, R.; Heydorff, M.; et al. An integrated environmental approach to investigate biomarker fluctuations in the blue mussel Mytilus edulis L. in the Vilaine estuary, France. Environ. Sci. Pollut. Res. 2012,20, 630–650. [CrossRef] [PubMed] 17. Fly, E.K.; Hilbish, T.J.; Wethey, D.; Rognstad, R.L. Physiology and Biogeography: The Response of European Mussels (Mytilus spp.) to Climate Change. Am. Malacol. Bull. 2015,33, 136–149. [CrossRef] 18. Kroeker, K.J.; Gaylord, B.; Hill, T.M.; Hosfelt, J.D.; Miller, S.H.; Sanford, E. The Role of Temperature in Determining Species’ Vulnerability to Ocean Acidification: A Case Study Using Mytilus galloprovincialis.PLoS ONE 2014 ,9, e100353. [CrossRef] [PubMed] 19. Zippay, M.L.; Helmuth, B. Effects of temperature change on mussel, Mytilus.Integr. Zoöl. 2012,7, 312–327. [CrossRef] 20. Gourault, M.; Petton, S.; Thomas, Y.; Pecquerie, L.; Marques, G.; Cassou, C.; Fleury, E.; Paulet, Y.-M.; Pouvreau, S. Modeling reproductive traits of an invasive bivalve species under contrasting climate scenarios from 1960 to 2100. J. Sea Res. 2019 ,143, 128–139. [CrossRef] 21. McCabe, M.; Navarrete, S. Reproductive investment in rocky intertidal mussels: Spatiotemporal variability and environmental determinants. Mar. Ecol. Prog. Ser. 2018,599, 107–124. [CrossRef] 22. Gosling, E. Reproduction, Settlement and Recruitment. In Bivalve Molluscs: Biology, Ecology and Culture; Fishing News Books: Oxford, UK, 2003; pp. 131–168. 23. Bernard, I.; de Kermoysan, G.; Pouvreau, S. Effect of phytoplankton and temperature on the reproduction of the Pacific oyster Crassostrea gigas: Investigation through DEB theory. J. Sea Res. 2011,66, 349–360. [CrossRef] 24. Seed, R.; Suchanek, T. Population and community ecology of Mytilus. In the Mussel Mytilus: Ecology, Physiology, Genetics and Culture; Gosling, E., Ed.; Elsevier: San Diego, CA, USA, 1992; pp. 87–170. 25. Gosling, E. (Ed.) Reproduction, Settlement and Recruitment. In Marine Bivalve Molluscs; John Wiley & Sons: Oxford, UK, 2015; pp. 157–202. 26. Starr, M.; Himmelman, J.H.; Therriault, J.-C. Direct Coupling of Marine Invertebrate Spawning with Phytoplankton Blooms. Science 1990,247, 1071–1074. [CrossRef] 27. Monaco, C.J.; McQuaid, C. Climate warming reduces the reproductive advantage of a globally invasive intertidal mussel. Biol. Invasions 2019,21, 2503–2516. [CrossRef] 28. Montalto, V.; Helmuth, B.; Ruti, P.M.; Dell’Aquila, A.; Rinaldi, A.; Sarà, G. A mechanistic approach reveals non linear effects of climate warming on mussels throughout the Mediterranean sea. Clim. Chang. 2016,139, 293–306. [CrossRef] 29. Sara, G.; Palmeri, V.; Rinaldi, A.C.; Montalto, V.; Helmuth, B. Predicting biological invasions in marine habitats through ecophysiological mechanistic models: A case study with the bivalve Brachidontes pharaonis.Divers. Distrib. 2013 ,19, 1235–1247. [CrossRef] 30. Sara, G.; Palmeri, V.; Montalto, V.; Rinaldi, A.; Widdows, J. Parameterisation of bivalve functional traits for mechanistic eco-physiological dynamic energy budget (DEB) models. Mar. Ecol. Prog. Ser. 2013,480, 99–117. [CrossRef] 31. Kooijman, B. Dynamic Energy Budget Theory for Metabolic Organisation; Cambridge University Press: Cambridge, UK, 2010. 32. Pörtner, H.-O.; Gutt, J. Impacts of Climate Variability and Change on (Marine) Animals: Physiological Underpinnings and Evolutionary Consequences. Integr. Comp. Biol. 2016,56, 31–44. [CrossRef] 33. Christensen, E.A.F.; Norin, T.; Tabak, I.; van Deurs, M.; Behrens, J.W. Effects of temperature on physiological performance and behavioral thermoregulation in an invasive fish, the round goby. J. Exp. Biol. 2021,224, 237669. [CrossRef] 34. Fearman, J.; Moltschaniwskyj, N. Warmer temperatures reduce rates of gametogenesis in temperate mussels, Mytilus galloprovincialis.Aquaculture 2010,305, 20–25. [CrossRef] 35. Cáceres-Martínez, J.; Figueras, A. Long-term survey on wild and cultured mussels (Mytilus galloprovincialis Lmk) reproductive cycles in the Ria de Vigo (NW Spain). Aquaculture 1998,162, 141–156. [CrossRef] 36. Villalba, A. Gametogenic cycle of cultured mussel, Mytilus galloprovincialis, in the bays of Galicia (N.W. Spain). Aquaculture 1995 , 130, 269–277. [CrossRef] 37. Sukhotin, A.A.; Flyachinskaya, L.P. Aging reduces reproductive success in mussels Mytilus edulis.Mech. Ageing Dev. 2009 ,130, 754–761. [CrossRef] [PubMed] 38. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2018. 39. Zuur, A.; Ieno, E.N.; Walker, N.; Saveliev, A.A.; Smith, G.M. Mixed Effects Models and Extensions in Ecology with R; Springer: New York, NY, USA, 2009. J. Mar. Sci. Eng. 2021,9, 759 17 of 17 40. Helmuth, B.; Harley, C.; Halpin, P.M.; O’Donnell, M.; Hofmann, G.E.; Blanchette, C.A. Climate Change and Latitudinal Patterns of Intertidal Thermal Stress. Science 2002,298, 1015–1017. [CrossRef] 41. Sara’, G.; Kearney, M.; Helmuth, B. Combining heat-transfer and energy budget models to predict thermal stress in Mediterranean intertidal mussels. Chem. Ecol. 2011,27, 135–145. [CrossRef] 42. Newell, R.I.E.; Hilbish, T.J.; Koehn, R.K.; Newell, C.J. Temporal Variation in the Reproductive Cycle of Mytilus edulisl. (Bivalvia, Mytilidae) from Localities on the East Coast of the United States. Biol. Bull. 1982,162, 299–310. [CrossRef] 43. Dowd, W.W.; Felton, C.A.; Heymann, H.M.; Kost, L.E.; Somero, G.N. Food availability, more than body temperature, drives correlated shifts in ATP-generating and antioxidant enzyme capacities in a population of intertidal mussels (Mytilus californianus). J. Exp. Mar. Biol. Ecol. 2013,449, 171–185. [CrossRef] 44. Bayne, B.L.; Salkeld, P.N.; Worrall, C.M. Reproductive effort and value in different populations of the marine mussel, Mytilus edulis L. Oecologia 1983,59, 18–26. [CrossRef] 45. Bayne, B.; Worrall, C. Growth and Production of Mussels Mytilus edulis from Two Populations. Mar. Ecol. Prog. Ser. 1980 ,3, 317–328. [CrossRef] 46. Sprung, M. Reproduction and fecundity of the mussel Mytilus edulis at helgoland (North sea). Helgol. Mar. Res. 1983 ,36, 243–255. [CrossRef] 47. Seed, R. The ecology of Mytilus edulis L. (Lamellibranchiata) on exposed rocky shores. Oecologia 1969 ,3, 277–316. [CrossRef] [PubMed] 48. Rodhouse, P.; McDonald, J.H.; Newell, R.I.E.; Koehn, R.K. Gamete production, somatic growth and multiple-locus enzyme heterozygosity in Mytilus edulis.Mar. Biol. 1986,90, 209–214. [CrossRef] 49. Beukema, J.J.; Dekker, R.; Jansen, J.M. Some like it cold: Populations of the tellinid bivalve Macoma balthica (L.) suffer in various ways from a warming climate. Mar. Ecol. Prog. Ser. 2009,384, 135–145. [CrossRef] 50. Heasman, M.; O’Connor, W.; Frazer, A. Temperature and nutrition as factors in conditioning broodstock of the commercial scallop Pecten fumatus Reeve. Aquaculture 1996,143, 75–90. [CrossRef] 51. Martínez, G.; Pérez, H. Effect of different temperature regimes on reproductive conditioning in the scallop Argopecten purpuratus. Aquaculture 2003,228, 153–167. [CrossRef] 52. Philippart, C.J.M.; Van Bleijswijk, J.D.L.; Kromkamp, J.C.; Zuur, A.F.; Herman, P. Reproductive phenology of coastal marine bivalves in a seasonal environment. J. Plankton Res. 2014,36, 1512–1527. [CrossRef] 53. Philippart, C.J.; Amaral, A.; Asmus, R.; Van Bleijswijk, J.; Bremner, J.; Buchholz, F.; Cabanellas-Reboredo, M.; Catarino, D.; Cattrijsse, A.; Charles, F.; et al. Spatial synchronies in the seasonal occurrence of larvae of oysters (Crassostrea gigas) and mussels (Mytilus edulis/galloprovincialis) in European coastal waters. Estuar. Coast. Shelf Sci. 2012,108, 52–63. [CrossRef] 54. Peteiro, L.G.; Labarta, U.; Fernández-Reiriz, M.; Alvarez-Salgado, X.; Filgueira, R.; Piedracoba, S. Influence of intermittentupwelling on Mytilus galloprovincialis settlement patterns in the Ría de Ares-Betanzos. Mar. Ecol. Prog. Ser. 2011 ,443, 111–127. [CrossRef] 55. Navarro, E.; Iglesias, J.I.P. Energetics of reproduction related to environmental variability in bivalve molluscs. Haliotis 1995 ,24, 43–55. 56. Bayne, B.L. Aspects of reproduction in bivalve molluscs. In Estuarine Processes; Wiley, M., Ed.; Academic Press: New York, NY, USA, 1976; pp. 432–448.