Main Drivers of Fecundity Variability of Mussels along a Latitudinal Gradient: Lessons to Apply for Future Climate Change Scenarios

: Bivalve relevance for ecosystem functioning and human food security emphasize the importance of predictions of mussel performance under different climate stressors. Here, we ad-dress 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 inﬂuence 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.


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][2][3][4][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][3][4][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][8][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 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][17][18][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][21][22][23][24]. Temperature is usually regarded as a trigger for spawning, but also regulates the rate of gonadal development with a cumulative effect [20][21][22][23][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][28][29][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][32][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 information might help to provide a better parameterization of reproduction, and to develop more accurate models for mussel populations forecasting.

Sampling Design
Mediterranean mussels, Mytilus galloprovincialis, were collected at 7 sites along the Portuguese coastline in 2014 and 2017 ( Figure 1)  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.

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.   , 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].

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: Vd Vs 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.

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, 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].

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).
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 parame ters (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 in clusion in the model. Model validation included the verification of homogeneity, normal ity and independence assumptions [39].

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 esti mated coefficients revealed warmer temperatures and lower salinities in 2017 (Table 1).   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).

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). A positive linear relationship between AF and female size was observed for both sampling years, in all locations but Sines in 2014 (Tables 3 and 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).  A positive linear relationship between AF and female size was observed for both sampling years, in all locations but Sines in 2014 (Tables 3 and 4). In this case, the lack of  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). 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 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 p-value; 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).

Discussion
Temperature latitudinal gradients are usually explored to understand the effects of thermal stress on marine invertebrates and to try to predict the subsequent impacts of future climate change scenarios over their distribution and population dynamics [21,[40][41][42]. Our study explores the reproductive output of mussel populations along a latitudinal gradient at the North Atlantic south limit of Mytilus galloprovincialis distribution. The sampled populations cover a temperature range of 1.5 °C from the northernmost to the southernmost locations (Table 1) which encompass the predicted rise of SST by 2100 under certain climate change scenarios (1.2 to 3.6 °C for RCP2.6 and RCP8.5 respectively) [1]. The temperature latitudinal gradient observed in our study area goes along with a trend of rising salinity and decreasing Chl-a availability from North to South (Table 1). Nonetheless, while the temperature differences between extreme locations was stable between years (1.5 °C even when average temperature was much larger in 2017 than in 2014), for salinity and Chl-a, spatial variability was not so consistent (Table 1). Maximum Chl-a and minimum salinity values were detected at the northernmost location in both years, even if lower riverine inputs and a much larger phytoplankton bloom during spring transition occurred in 2014 (Figure 2), but the progressive trend towards the south was not monotonic (Table 1). Other local factors like riverine inputs, upwelling intensity or topography

Discussion
Temperature latitudinal gradients are usually explored to understand the effects of thermal stress on marine invertebrates and to try to predict the subsequent impacts of future climate change scenarios over their distribution and population dynamics [21,[40][41][42]. Our study explores the reproductive output of mussel populations along a latitudinal gradient at the North Atlantic south limit of Mytilus galloprovincialis distribution. The sampled populations cover a temperature range of 1.5 • C from the northernmost to the southernmost locations (Table 1) which encompass the predicted rise of SST by 2100 under certain climate change scenarios (1.2 to 3.6 • C for RCP2.6 and RCP8.5 respectively) [1]. The temperature latitudinal gradient observed in our study area goes along with a trend of rising salinity and decreasing Chl-a availability from North to South (Table 1). Nonetheless, while the temperature differences between extreme locations was stable between years (1.5 • C even when average temperature was much larger in 2017 than in 2014), for salinity and Chl-a, spatial variability was not so consistent (Table 1). Maximum Chl-a and minimum salinity values were detected at the northernmost location in both years, even if lower riverine inputs and a much larger phytoplankton bloom during spring transition occurred in 2014 (Figure 2), but the progressive trend towards the south was not monotonic (Table 1). Other local factors like riverine inputs, upwelling intensity or topography may have a larger influence over the spatio-temporal variability of those variables, preventing in many cases a direct correlation of Chl-a and salinity with latitude [21].
The large spatio-temporal variability on relevant environmental variables (food availability, wave intensity, air exposure, topographical characteristics, etc.) coupled to complex interactions and trade-offs between them and temperature, has been frequently suggested as an explanation for the lack of detection of latitudinal gradients on physiological performance attributable to an increase of thermal stress [21,[40][41][42]. Particularly, for mussel reproductive output, the lack of detection of a thermal stress latitudinal gradient is commonly attributed to the large influence of food availability on the reproductive processes [7,20,23]. Food availability regulates the number of oocytes produced [25], but also modulates stress response capability and metabolic activity in mussels, increasing their ability to cope with different sources of stress [43]. Some authors identify more complex interactions between environmental variables, but associated to oceanographic processes which lead to the persistence of multivariate mosaics of optimal environmental conditions and therefore maximum reproductive investment [21].
Our results showed that 3 locations in the Central-North coast of Portugal (Peniche, Cascais and Costa Nova) registered the largest absolute fecundity values during both years ( Figure 3A). Egg production per individual was much higher in 2017 than in 2014, suggesting a larger influence of temperature over food availability for the reproductive output of those populations, since 2017 was characterized by warmer temperatures and lower food availability (Table 1; Figure 2). Nonetheless, the higher values observed for AF at those locations, both during 2014 and 2017, seems to be driven by the larger size of mussels at these three populations ( Figure 3B). Mussels, as well as many iteroparous invertebrates, show a direct relationship between reproductive effort and body size, which is usually explained as an age-related increase on the proportion of energy allocated to reproduction [24,25,44,45]. Our results highlight the relationship between size and AF for all the populations and years studied (Table 3). Nonetheless, slopes describing the relationship between AF and size varied between locations ( Figure 3C), which might be a consequence of the amount of stress experienced at each population. Previous studies have suggested that environmental stress can push populations of mussels to reach reproductive maturity at smaller sizes [46,47], as well as to attain earlier their maximum rate of reproductive investment, resulting in steeper slopes between reproductive effort and size [44]. Optimal environmental conditions usually lead to less stressed and faster growing populations, which reach reproductive maturity at larger sizes but also attain larger life spans [44][45][46][47]. Several authors have reported a decay on reproductive output on mussels over 50 mm, related to senescence processes [37,44,48], which can also flatten the slope of the relationship between oocyte production and size in populations where average size was larger ( Figure 3). A combination of both processes might determine the relationship of reproductive output with size at each population [44], and agree with the hypothesis of persistent multivariate mosaics of environmental conditions which determine the adaptation of each population to their local habitat [21], and comply with the consistence of the slopes calculated for each location between years in the present study (Table 4; Figure 3C).
The standardization of AF by gonad weight (RF) allows for the integration of size on the reproductive output of each population, and for the detection of further latitudinal patterns which vary between years (Figure 4). In 2014, RF increased almost exponentially towards the south, suggesting almost a linear effect of temperature on reproductive output over certain thresholds (Table 1; Figure 2). Nonetheless, in 2017 maximum RF values were displaced towards the central coast of Portugal, with a sharp drop of RF farther South (Figure 4), which seems to indicate that the warmer temperatures detected that year (Table 1; Figure 2) overpassed reproductive optimums at those locations, even when temperatures recorded were far from the upper thermal limit reported for the species (≈27-35 • C; [7,28]).
Gonadal development of mussels and many other bivalves takes place over a range of temperatures which are usually well below their thermal tolerance limit for survival [34,[49][50][51]. Temperature influences metabolic rates, increasing energetic cost and limiting the energy available for reproduction [24]. For M. galloprovincialis, energy allocation for reproduction seems to decrease fast when temperature rises over the optimum for somatic maintenance (≈17.5 • C; [28]), leading to a reproduction limiting temperature around 19 • C [34]. RF latitudinal gradients observed in the present study ( Figure 4) might reflect the thermal history of mussels during gametogenesis according to those limits, reaching maximum values at those locations and years where temperature (Figure 2) was between the reported optimum (≈14 • C; [34]) for M. galloprovincialis gamete production and its upper limit (19 • C; [34]). Food availability can compensate for the extra demand of energy under certain stress conditions [43], but several studies report decreased rates of gametogenesis under high temperatures for bivalves, even under no limiting food conditions [34,50,51].
Our GAM model also reflects a larger effect of the temperature history over the food availability on gamete production (Table 6), supporting the role of temperature as a switching-on/-off mechanism for gametogenesis ( Figure 5A). The number of days with temperatures over the optimum for gamete production during gametogenesis (4 months) explained the largest amount of the variability observed on RF (p-value < 2.9 × 10 −5 ; Table 6). RF increases slightly with the number of warm days until those overpasses certain limit (≈80 days; Figure 5A) leading to a sharp drop on reproductive output. "Temperature windows" for gametogenesis are well documented for bivalves [25] and their evolution is usually associated with the advantage of matching reproduction with optimal conditions for larval development [52]. For M. galloprovincialis in the Iberian Coast, two peaks of reproduction are commonly reported (spring and autumn) [36,53], when food availability and oceanographic conditions optimize larval survival [52,54]. Therefore, having temperature thresholds for gametogenesis to avoid energy allocation for gamete production during winter and summer, when larval survival is limited by environmental conditions, might be a good strategy for energy optimization.
It is well known that temperature has multiple effects over bivalves' gametogenesis, first acting as an on/off mechanism, and once activated regulating its kinetics [55]. Our results also highlighted a positive effect of temperature over RF when it rises over 16 • C during the last month of gametogenesis ( Figure 5B). The average temperature values recorded during our study did not reach the reproduction limiting temperature (19 • C [34]) at any time ( Figure 5B). Therefore, the positive relationship observed between temperature and RF during the last month of gametogenesis is not expected to be maintained over that temperature threshold.
We also detected a linear and positive effect of Chl-a concentration during the last month of gametogenesis over RF ( Figure 5C), while the maximum Chl-a attained during the spring bloom had no significant effect on the reproductive output ( Figure 5D; Table 6). These results might indicate that mussels are taking an opportunistic reproductive strategy, since the reproductive output is only dependent of the most recent food concentration. Food availability has a direct effect over the amount and quality of gametes produced by bivalves, but according to their reproductive strategy, species can be classified as opportunistic or conservative [24,56]. In opportunistic species the gametogenesis is linked to the food supply and occurs when feeding conditions are favourable, while in the conservative species the nutrients are stored before gametogenesis is initiated [24]. Mussels are, however, quite flexible and can exhibit both types of strategy, adapting it to the prevailing environmental conditions [24]. For some authors, gametogenesis initiation is mostly limited by nutrient availability, either as reserves or recently ingested food [45], while for others it is the temperature window which determines initiation/cessation of gametogenesis, and food availability modulates the amount of gametes produced [24]. Our results support better the second hypothesis, since we detected a much lower influence of food availability than temperature over RF ( Table 6). The patterns observed ( Figure 5) also agree with a major role of accumulative temperature conditions acting as a switching-on/off mechanism for gametogenesis, while more recent environmental conditions (temperature and food availability) regulate gamete production rates.
IPPC predicted a rise of SST by 2100 ranges from 1.2 to 3.6 • C (for RCP2.6 and RCP8.5 respectively; [1]). Our results suggest that under these climate predictions, temperature windows for gametogenesis will limit the number of reproductive events and might also alter the reproductive phenology of mussels to avoid persistent warm temperatures which will start earlier in the year and last for longer. Many studies agree on predicting changes on reproductive phenology of mussels [7,9,28,52], because an earlier activation of gametogenesis but an upper temperature limit controlling the decline of gametogenesis is not commonly included. Therefore, even when mismatch effects of dislodging reproduction towards less favourable conditions for larval survival might negatively affect mussel population dynamics [52,53], many studies based on DEB models consistently predict a maintenance or a slight increase of reproductive output for bivalves, and an increase of the population performance in all the locations where warming future scenarios don't overpass thermal lethal tolerance of the species [7,9,20,23,[27][28][29][30]. As mentioned before, those results might be related to the lack of specific parametrization of non-linear responses of gamete production to environmental conditions. DEB models establish general rules and priorities for energy allocation (somatic maintenance > growth > maturity maintenance > reproduction), and estimations of the cost of different physiological processes [31]. In addition, DEB models establish a series of assumptions for bivalve reproduction which usually include: a threshold of degree-days to activate gametogenesis, a threshold of gonadosomatic index to achieve ripe condition and a temperature threshold to trigger spawning [7,9,20,23,[27][28][29][30]. However, those models do not include an upper limit for gametogenesis [9,20,23,[27][28][29][30]. With the exception of some studies which incorporate a thermal upper limit for spawning [7], the lack of a specific upper limit for gametogenesis and the assumption that reproduction follows the same thermal tolerance rules than other physiological processes [31] forces the predictions towards large reproductive outputs [9,20,23,[27][28][29][30], even at temperatures over the reported limit for reproduction [34]. The lack of correct parameterization of gametogenesis can therefore overestimate mussel population performance under future climate change scenarios.
Our results deepen the understanding of the underlying processes controlling the complex interactions between environmental drivers and gamete production, particularly highlighting the relevance of considering the thermal window for gametogenesis. Our study also points out the relevance of a correct parameterization of those processes, to be able to perform accurate predictions of population performance under future climate change scenarios.

Conclusions
Our study supports the view that gametogenesis responds non-linearly to temperature and chlorophyll concentrations, and that these variables show different effects over fecundity depending on the temporal scale evaluated. We found a major role of long-term temperature, acting as a switching on-off mechanism for gametogenesis, while short-term food availability has a lower influence but also modulates the amount of gametes produced.
Our results highlight the complexity of gametogenesis environmental dependences, and how its inclusion in predictive models is vital to avoid overestimations on the capability of mussel populations to deal with climate change future scenarios. Institutional Review Board Statement: Ethical review and approval were waived for this study because this is a very abundant and non-endangered species. Samples taken were only those strictly necessary for the analyses, and were made under permission of the Agência Portuguesa do Ambiente and the Instituto de Conservação da Biodiversidade e Florestas.

Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets generated and analyzed during this study are available from the corresponding author upon request.