Reproduction Under Stress: Acute Effect of Low Salinities and Heat Waves on Reproductive Cycle of Four Ecologically and Commercially Important Bivalves

The impacts of climate change on the structure and functioning of estuaries is a major focus of concern, even more when the affected species support important fisheries as the bivalves Ruditapes decussatus, Ruditapes philippinarum, Venerupis corrugata, and Cerastoderma edule in Europe. Their reproductive performance, in the context of climate stressors, had not been investigated so far. Our objective was to experimentally evaluate acute stress effects over gonad development after 6 days of low salinity stress in autumn, winter and spring as well as 4 days of heatwave stress during emersion in summer. These are the most probable extreme events that bivalves should face in our latitudes. Four different salinity ramps (5–20, 10–25, 15–30, 30–30) were created during simulated tidal cycles in mesocosms for the low salinity experiments. Also four sediment heatwaves at emersion (20–20, 20–27, 20–32, 20–37°C) were done during simulated tidal cycle. Both low salinity and heatwave stresses over such short periods compromised reproduction; the acute response was species-specific and varied with the time of the year, and therefore, with the stage of the gametogenic cycle. In December, during sexual resting and the beginning of gametogenesis, a delay in gametogenesis at lower salinities was recorded in the four species. However, at the peak of the reproductive period (March and May), different responses were observed: abnormal oocytes in R. decussatus and resorption of gametes with haemocytic infiltration in R. philippinarum and V. corrugata. Likewise sediment temperatures higher than 32°C provoked gonadal resorption and severe haemocytic invasion in V. corrugata, R. decussatus, and C. edule but had no effect in R. philippinarum. These responses to both environmental stressors might be related to the allocation of energy from reproduction toward defense and repair mechanisms to ensure survival. Contrastingly, low salinities triggered massive spawning in C. edule that could lead to a mismatch between the presence of larvae and phytoplankton, causing potentially starvation and thus reducing recruitment success. Reproduction of theses bivalves would be compromised if low salinity episodes in winter and spring, even for short periods of time such as those in these experiments, are followed by a heatwave in summer. Furthermore the impact would be magnified if this situation happens during consecutive years preventing replenishment of the shellfish beds.

The impacts of climate change on the structure and functioning of estuaries is a major focus of concern, even more when the affected species support important fisheries as the bivalves Ruditapes decussatus, Ruditapes philippinarum, Venerupis corrugata, and Cerastoderma edule in Europe. Their reproductive performance, in the context of climate stressors, had not been investigated so far. Our objective was to experimentally evaluate acute stress effects over gonad development after 6 days of low salinity stress in autumn, winter and spring as well as 4 days of heatwave stress during emersion in summer. These are the most probable extreme events that bivalves should face in our latitudes. Four different salinity ramps (5-20, 10-25, 15-30, 30-30) were created during simulated tidal cycles in mesocosms for the low salinity experiments. Also four sediment heatwaves at emersion (20-20, 20-27, 20-32, 20-37 • C) were done during simulated tidal cycle. Both low salinity and heatwave stresses over such short periods compromised reproduction; the acute response was species-specific and varied with the time of the year, and therefore, with the stage of the gametogenic cycle. In December, during sexual resting and the beginning of gametogenesis, a delay in gametogenesis at lower salinities was recorded in the four species. However, at the peak of the reproductive period (March and May), different responses were observed: abnormal oocytes in R. decussatus and resorption of gametes with haemocytic infiltration in R. philippinarum and V. corrugata. Likewise sediment temperatures higher than 32 • C provoked gonadal resorption and severe haemocytic invasion in V. corrugata, R. decussatus, and C. edule but had no effect in R. philippinarum. These responses to both environmental stressors might be related to the allocation of energy from reproduction toward defense and repair mechanisms to ensure survival. Contrastingly, low salinities triggered massive spawning in C. edule that could lead to a mismatch between the presence of larvae and phytoplankton, causing potentially starvation and thus reducing recruitment success. Reproduction of theses INTRODUCTION Current predictions show that on the Atlantic coast of Europe there will be "very likely" more frequent heat waves, longer in duration and greater in intensity, as well as increases in the intensity and frequency of extreme precipitations that will modify coastal salinity (Christensen et al., 2007;Trenberth, 2012;IPCC, 2014;Jacob et al., 2014;Hoegh-Guldberg et al., 2018;Viceto et al., 2019;Carvalho et al., 2021; but see Cardoso Pereira et al., 2020;Lorenzo and Álvarez, 2020). In fact, from 1982 to 2010 almost half of the world's coastlines had already experienced a significant decrease in the frequency of extremely cold days and more than a third had faced an increase in extremely hot days (Lima and Wethey, 2012).
These acute events are important in shaping communities and their effects can be widespread and long lasting provoking important biogeographical changes (e.g., Daufresne et al., 2007;Mulholland et al., 2009;Byrnes et al., 2011;Wethey et al., 2011;Frölicher and Laufkötter, 2018;Shanks et al., 2020). Temperature affects primarily physiology and metabolism of ectotherms resulting in changes in energy allocation among metabolism, maintenance, growth, reproduction, and protective responses such as heat shock proteins or antioxidant pathways, which require energy (Pörtner and Farrell, 2008;Eymann et al., 2020). Likewise, changes in salinity produce energetic demands of osmoregulation with allocation of energy in growth and cellular maintenance at the expense of reproduction (Brett, 1979).
The reproductive cycle of bivalves is controlled by the interaction between endogenous and environmental factors (Normand et al., 2008). Among those environmental factors, temperature is a forcing variable for the gametogenic cycle in bivalves (e.g., Sastry, 1966;Xie and Burnell, 1994;Albentosa et al., 2007) defining both the starting point and the rate of gonadal development, whereas food influences the duration of the reproductive cycle (Lubet, 1959;Ruíz et al., 1992) and the percentage of mature individuals (Delgado and Pérez-Camacho, 2007). It has been demonstrated that under increased stresses, such as those predicted from climate change scenarios, bivalves may allocate energy away from growth and reproduction toward costly physiological defenses in order to survive (Petes et al., 2008). And the potential for fecundity loss is one of the most likely consequences of that scenario. For instance, elevated temperature reduces gamete development in Mytilus galloprovincialis (Múgica et al., 2015) and fertilization success in Tridacna maxima (Armstrong et al., 2020); abrupt changes in salinity affect spawning of Ruditapes philippinarum (Dang et al., 2010). Thus, the very energetically costly reproductive processes (Gosling, 2003) may be impaired or suppressed (Wingfield and Sapolsky, 2003) with enormous negative consequences for the structure and dynamics of populations of coastal organisms. Failures in recruitment of juveniles will prevent replenishment of shellfish beds (Beukema and Dekker, 2005;Petes et al., 2007;Talmage and Gobler, 2011), and, ultimately, the decline or local extinction of populations (Parmesan, 2006;Berg et al., 2010).
Small-scale fisheries (SSFs, Small-Scale and Spatially Structured Fisheries targeting sedentary resources with artisanal gears, sensu Orensanz et al., 2005) are of enormous importance to both harvest and employment (Mills et al., 2011). This is particularly true in Galicia (NW of Spain) where fisheries in the intertidal and shallow subtidal of the native clams Ruditapes decussatus (Linnaeus, 1758) (grooved carpet shell) and Venerupis corrugata (Gmelin, 1791) (pullet carpet shell), the introduced species R. philippinarum (Adams and Reeve, 1850) (Manila clam), and the cockle Cerastoderma edule (Linnaeus, 1758) contribute an annual value of ∼74 million € and involve ∼7,100 fishers. 1 These four bivalves are dioecious broadcast spawners without sexual dimorphism, each displaying a different reproductive cycle. In Galicia (Spain) the reproductive cycle of R. decussatus starts in January with the initiation of gametogenesis although the proliferation of gametes occurs during spring, followed by spawning in summer (June-August) (Rodríguez-Moscoso, 2000;Ojea et al., 2004). This species has a great capacity for gonadal regeneration thus allowing multiple spawnings per season (Matias et al., 2013). Gametogenic reabsorption after spawning seems to be unimportant (Rodríguez-Moscoso, 2000) with only some degenerate oocytes (cytolysis) and no major hemocyte infiltration. Resting phase happens in November-December with the presence of reserve tissue formed by interfollicular muscle and vesicular cells (polygonal cells with an eccentric nucleus) inside the follicles, which disappear as the gonadal development advances (Delgado and Pérez-Camacho, 2007).
Ruditapes philippinarum, native to the Pacific, was introduced for culture in Europe in the 1970s because of its high adaptability (Latrouite and Claude, 1976;Pérez-Camacho and Cuña, 1985;Drummond et al., 2006). Nowadays it is one of the most commercially exploited molluscs in the world. It has a long reproductive period and its gonadal development is quite asynchronous within the gonad of each individual and between individuals. It is in sexual resting from October until January, when the gametogenesis starts, reaching advanced gametogenesis in March. In April the gametes are mature and a period of successive spawning interspersed with gonad recovery lasts from April to September (Rodríguez-Moscoso et al., 1996;Delgado and Pérez-Camacho, 2007).
In Galicia V. corrugata shows a long reproductive cycle (Villalba et al., 1993) with ripe and spawning stages observed throughout the year, although with larger number of ripe follicles between February and June (Cerviño-Otero et al., 2007). During autumn the gonad develops asynchronously, showing simultaneously in the same gonad mature and early gametogenic follicles. Gametogenesis occurs between February and March, reaching maturity in March although spawning starts in May continuing until July. In August and September a new maturation period begins with sporadic spawning in autumn. Sexual rest is very occasional since as soon as residual gametes appear, vesicle cells are formed and new germinal lines develop (Cerviño-Otero, 2011).
The onset of gametogenesis in the cockle C. edule in Galicia takes place at the end of the summer (September to October), progresses throughout the winter, and reaches maturation in the spring. The first spawning occurs in April and May and, after gonad restoration, another spawning episode takes place in May and June. During July and August, most of the population shows signs of gonad resorption, although a less-intensive spawning event can occur at the end of summer/beginning of autumn (Martínez-Castro and Vázquez, 2012).
The physiological responses to temperature and salinity fluctuations of these four bivalves have been recently investigated by our research team in mesocosm experiments in which field conditions were simulated (Macho et al., 2016;Peteiro et al., 2018;Domínguez et al., 2020). Early recruits of C. edule showed a lethal limit at salinity of ∼15 and a marked reduction of activity (i.e., continuous valve closure and inhibition of respiration and ammonium excretion rates) at salinities between 15 and 30 (Peteiro et al., 2018). After 6 days of exposure to different salinity fluctuations, adults of C. edule, R. decussatus, R. philippinarum, and V. corrugata showed a dramatic reduction of pumping activity, scope for growth (SFG) and burrowing activity at the lowest salinities (i.e., 5, 10, and 15), but responses under higher salinity (20, 25, and 30) varied among species and across seasonal periods . While C. edule was the most affected species in autumn showing no fast recovery after the end of the stress episode, R. decussatus was more resistant in all seasons .
Similarly, temperature stress experiment performed for these species, agree to identify acute responses to short heatwaves (Domínguez et al., 2021). After 2 days of exposure to sediment heatwave during low tide, C. edule and V. corrugata suffered significant mortalities and also a dramatic decrease in SFG as well as reduction in burrowing activity compared to R. decussatus and R. philipinnarum (Domínguez et al., 2021).
While the effects of acute environmental stress on the SFG and behavior (see above), or early-life history stages (e.g., Talmage and Gobler, 2011) in bivalves have been documented, the effect on the reproductive performance in the context of climate stressors has not yet been investigated (but see Xu et al., 2016) even though the reproductive response of clams and cockles under future climate change scenarios is crucial for evaluating population dynamics (Morgan et al., 2013;Xu et al., 2016).
Thus, the aim of the present study was to experimentally evaluate the potential effect of short episodes of low salinity as well as sediment heatwaves on the reproductive cycle of these bivalves. Both low salinity and temperature stress treatments consisted of four ramps that followed similar profiles to those experienced by bivalves in the field during tidal cycles. Consequently, a low salinity stress experiment was repeated in autumn, winter and spring when episodes of heavy rains occur in the Atlantic coast of Europe while a heatwave experiment was performed during summer. Considering the available physiological data from previous studies, it is reasonable to hypothesize that each of these four species may respond differently to increasing environmental stress, potentially impairing or suppressing reproductive success. Also, we expect a more acute effect during gametogenesis and spawning than during the resting period of the gametogenic cycle.

MATERIALS AND METHODS
Four experiments were performed in a mesocosm system at Estación de Ciencias Mariñas de Toralla (ECIMAT) 2 of the Universidade de Vigo (Galicia, NW Spain). Low salinity stress experiments were independently run in December 2015, March and May 2016; the heatwave stress experiment was done in July 2016.

Experimental Setup
Four small plastic tanks, one per species, (16 L, 21 cm tall × 30 cm width × 40 cm length) with four 2 cm bottom orifices covered with 80 µm mesh to allow water flux and filled to the top with sediment to prevent any water remaining at low tide, were placed inside each of eight big tanks (480 L, 50 cm tall × 80 cm width × 120 cm length) (Figures 1, 2). The sediment was collected from an intertidal shellfish bed (42 • 11.68 N; 8 • 47.81 W) (median grain size of 0.19 mm); old sediment was discarded and new sediment was collected for each experiment. These big tanks, which were supplied with running 50 µm-filtered seawater that entered via inlets in the bottom of the large tanks and exited via ∼30 cm tall standpipes, were randomly placed in a controlled temperature room at 18 • C. There were two big tanks for each experimental treatment (Figure 2). During daylight hours the ebb tides were simulated every day in each tank by drainage followed by an emersion period. Water salinity and temperature were recorded by mini-CTDs (Star Oddi) and sediment temperature by iButtons at the surface of the sediment and at 5 cm depth within the sediment.
The day before starting each experiment, adult clams and cockles of standardized length (clams: ca. 40 mm; cockles: 30 mm) to minimize bias from size, were manually collected in the shellfish beds of Cambados (42 • 30 55 N; 08 • 48 53 W) and General setting for heatwave stress experiment with the ceramic lamps in the position during low tide and control tanks. bp, dual bellows pumps; cl, ceramic lamps; et, experimental tanks; ht, head tanks with fresh and salt water; i, inlet; mt, mixing tanks for salinity treatments (S5, S10, S15, S30); to, outflow standpipe; st, 16 L species tanks (one per species). Plastic beakers (b) were used to study behavioral responses of the bivalves using pressure sensors; those results were published in Woodin et al. (2020).
Noia (42 • 47 0 N; 8 • 53 0 W) respectively, and transported to the laboratory in refrigerated coolers. A total of 28 individuals of each species, 30 in the December experiment, were marked and placed in each 16 L tank at approximately the densities found in shellfish beds (220 ind. m −2 ) and allowed to burrow. Those that did not burrow within 8 h were discarded and replaced by new individuals. Animals were fed in the evenings during the day before the experiment and during the experiment with a microalgae mixture of Isochrysis galbana, Tetraselmis suecica, Chaetoceros gracilis, and Rhodomonas lens, constituting a 1% diet based on dry weight of the clams (a dry weight of 0.88 g was assumed for each individual based on the mean size).

Low-Salinity Stress Experiments
Treatments consisted of four salinity ramps in which salinity varied from 5 to 20, 10 to 25, 15 to 30, and the control treatment 30 to 30 (hereinafter S5, S10, S15, and S30), salinity profiles similar to those experienced by bivalves in the field ( Figure 3A). Salinity ramps were created automatically through the use of timers controlling dual bellows pumps (Iwaki DP80-30) that mixed dechlorinated fresh water and 50 µm-filtered sea water at ambient temperature in different proportions. After the daylight ebb tide produced by drainage, an automatic change in water source and thus salinity was performed (∼1.5 h). During night tides the automatic change in water source occurred but was not preceded by drainage so salinity change occurred more slowly. Each salinity ramp was desynchronized 1.5 h, the time to process the bivalves, so all animals stayed the same time at low tide ( Figure 3B). The lower salinity occurred during the ebb tide period and the higher salinity during the flood tide period as it happens in the field.
The salinity stress was carried on during six consecutive days (D6) and, just in the December 2015 experiment, individuals were left three more days at control salinity (S30) to measure their ability to recover (D9).

Heatwave Stress Experiment
Temperature treatments, similar to those experienced by bivalves in the field during heat waves in previous years (Macho et al., 2016;Domínguez et al., 2021;Figure 3C), were set at 20, 27, 32, and 37 • C (hereafter T20, T27, T32, and T37, respectively). To produce experimental temperatures we used infrared ceramic lamps (FTE 150-watt heaters, Ceramicx) over each tank controlled by digital controllers (Omega CN7853) with feedback thermocouples placed at 2 cm depth in the sediment in the middle of each tank. All big tanks experienced simulation of daily tidal cycles at ambient temperature (18-22 • C) at continuous flow and the heat waves were simulated during the emersion period of the diurnal ebb tide (∼4.5 h) ( Figure 3D). The temperature treatments were applied on four consecutive days (stress period) (D4).

Histological Analysis
To determine gonadal stage, 20 animals of each species (10 from each of the two tanks at S30) were dissected the day before starting the stress (D0). In the December experiment, 10 animals of each species and treatment per replicate at the end of the Frontiers in Marine Science | www.frontiersin.org stress (D6) and after recovery (D9) were dissected. In the rest of the experiments, due to time constrains of the histological processing, only 6 individuals of each species per treatment and replicated tank were collected at the end of stress (D6 in the low-salinity and D4 in the heatwave experiment). Since in bivalves the gonad is a diffuse tissue, a piece of tissue of ca. 1 cm −2 was dissected from the foot and routinely processed for histology: i.e., fixed in Davidson formaldehyde for 24 h, rinsed with running water for 15 min, dehydrated in an ethanol series (automatic tissue processor Leica TP1020), immediately embedded in paraffin (Paraffin embedding station Leica EG1150H), sectioned at 5 µm (rotary microtome Leica RM2255) and stained with hematoxylin and eosin (multistainer Leica ST5020).
The histologically prepared slides (at least four sections on each slide) were examined under a microscope and each specimen was assigned to a stage which represented the gonadal development. When more than one developmental stage was observed within a single individual, the assignment of a stage criterion decision was based upon the condition of the majority of the follicles of the four sections of the slide. The stages of gonad development were scored according to the scales proposed by Delgado and Pérez-Camacho (2005) for R. decussatus, Holland and Chew (1974) for R. philippinarum, Cerviño-Otero (2011) and Joaquim et al. (2011) for V. corrugata, andIglesias (2006) and Martínez-Castro and Vázquez (2012) for C. edule. To unify the different names of the categories of each species and for clarification and comparative purposes, we defined the following stages: Sexual rest: follicles are scarce, isolated and small. Gonadal follicles are absent and connective tissue occupies the gonadal gland. There is no evidence of gonadal development and sex determination is not possible.
Early gametogenesis: follicles become evident, increasing in size and number. Their walls are full of germ cells: oogonia and previtellogenic oocytes in females and spermatogonia and primary spermatocytes in males.
Late gametogenesis: follicles occupy most of the visceral mass, and germinal cells are present at all stages of gametogenesis. In females, free oocytes in the lumen and immature oocytes attached to the basal membrane are observed at different stages of vitellogenesis. The abundance of free oocytes equals those attached to the wall of the follicle. In males, spermatogonia, spermatocytes, spermatids, and some radially arranged spermatozoa are present.
Ripe: Follicles are almost full of ripe gametes. In females, ripe oocytes are free in the lumen. In males, spermatozoa are distributed radially with the tails pointing toward the center of the follicle.
Spawning: As a result of the release of the gamete, pressure inside the follicle decreases. Depending on the degree of spawning the follicles are more or less empty. In females, empty spaces in the follicular lumen are observed, although some ripe oocytes are present. In males, spermatozoa lose their radial arrangement and follicles are partly empty.
Gonad restoration for R. philippinarum, V. corrugata, and C. edule: after spawning, a new gametogenic cycle begins in the follicles and new germ cells appear in the follicle walls. This stage is similar to advanced gametogenesis, but new germ cells co-exist with ripe gametes. In females, oogonia in division are in the follicle walls, numerous previtellogenic oocytes and scarce free ripe oocytes in the lumen are observed. In males, the number of spermatocytes increases again, in contrast to the small number of spermatozoa.
Resorption is the process of follicular recession process after the last spawning of the reproductive cycle before initiating a new period of sexual rest. Cytolysis begins in the gonad (hemocytes are common) and follicles become very small and practically empty. In females, some oocytes, showing clear signs of atresia, are present and, in males, some spermatozoa remain. In R. philippinarum, vesicle cells, that storage glycogen for the gametogenesis, appear into the follicles.
A mean Gonadal Index (GI) was calculated modifying the method proposed by Seed (1976) and Matias et al. (2013): GI = ( individuals in each stage × stage ranking)/ total individuals in each treatment For each of the gonadal stages a numerical ranking was assigned: sexual rest (0); early gametogenesis (3); late gametogenesis (4); ripe (5); spawning (2); gonadal restoration (4) since this stage is similar to advanced gametogenesis; resorption (1). The GI ranged from 0 (all individuals in sexual rest) to 5 (all individuals ripe).
Variation in the frequency of gonad maturation stages with salinty and temperature was analyzed with multinomial regressions. Analyses were performed separately for each species, day of measurement (D4, D6, or D9 depending on the experiment) and, in the case of salinity experiments, each experimental period (December, March, and May) was also analyzed separately. Multinomial logistic regressions were performed with the multinom function from the nnet package (Fox and Weisberg, 2019) and poshoc pairwise analyses were run using the package emmeans (Estimated Marginal Means, aka Least-Squares Means) (Lent, 2020). Analyses were made in R version 3.6.1 (R Core Team, 2019) with R studio interface (RStudio Team, 2019).
Gonadal Index of each experiment and species was not formally analyzed, but graphically shown, due to the low number of replicates, i.e., 2 tanks (n = 2).
No mortality was recorded during the salinity stress experiments. During the heatwave experiments, mortality of R. decussatus was less than 4% in all treatments. R. philippinarum only suffered mortality of 2% at temperatures <37 • C but 30% at 37 • C treatment. Mortality of V. corrugata was 12% at temperatures <37 • C and 40% at 37 • C treatment. The greatest mortality was experienced by C. edule with more than 75% at 37 • C and 12% at temperatures <37 • C.
At the beginning of the experiment, all individuals of R. decussatus were at sexual rest stage with all their gonadal tissue undifferentiated, thus a GI of 0; after 6 days of stress, 30% of the clams kept at the control salinity (S30) initiated gametogenesis while within all low salinity treatments only one (S5) or two clams (S10, S15) initiated gametogenesis ( Figure 4A). These results could be also observed in the GI with higher values at S30 (Figure 4B). After 3 days of recovery (D9) the gametogenetic stages did not change (Figures 4A,B), and the effect of the previous low salinity stress persisted (χ 2 = 16.086, df 3, p < 0.001).
A similar trend was observed in the other two venerid clams, consistent with more gonadal development at higher salinity. Around 50% of R. philippinarum individuals started the experiment resorbing their gonads with abundant intragonadic vesicular cells and 50% of inactive gonads with abundant interfollicular connective tissue. After 6 days of stress at S5, processes of resorption of the gonad continued since more individuals were at sexual resting whereas at higher salinities, mainly at S15 and S30, gametogenesis started with the appearance of the first follicles (marginal differences in early gametogenesis; post hoc S30-S5 df 8, p = 0.054), which was reflected in the GI, with an increase at S15 and S30 (Figures 4C,D). After recovery, gonads advanced toward sexual resting without significant difference among treatments (χ 2 = 5.872, df 9, p = 0.753), which was also shown by the lower GI value on D9 compared to D6.
At the beginning of the experiment the majority of V. corrugata were resorbing gonads although it was common to observe different stages of gonadal development in the same individuals. Around 20% of clams restarted gametogenesis with development of new gametes ( Figure 4E). After stress, the proportion of sexual resting animals diminished in all treatments although it was more evident at S30; actually 6 clams at S30 underwent spawning, probably the 25% of individuals that were in gametogenesis at D0 (marginal differences in spawning; post hoc S30-S5 df 8, p = 0.071 and S30-S15 df = 8, p = 0.071). On D9 gametogenesis progressed at S5 and S15 without differences between treatments (χ 2 = 14.788, df 12, p = 0.253) as was also shown by higher GI values (Figures 4E,F).
On D0, C. edule did not show follicles (sexual rest) or these were small with connective and muscular tissue occupying the entire zone from the digestive gland to foot. Although statistical differences were not found between salinity treatments (χ 2 = 7.012, df 6, p = 0.310) gametogenesis progressed faster at S10, S15, and S30 (Figures 4G,H). Remarkably, after the recovery period, gametogenesis at the lowest salinity advanced to almost reach the same gonadal stage as in the other salinities (χ 2 = 6.994, df 6, p = 0.321).

March
All clams except R. decussatus were primarily at advanced stages of gametogenesis including individuals with fully mature gametes such as V. corrugata with a GI higher than 4, or even in the spawning period such as C. edule and R. philippinarum. Contrastingly, R. decussatus was at the early stages of gametogenesis (D0, Figure 5).
Gonads of R. decussatus remained at early gametogenesis at all salinities except in the control where 20% of individuals were at advanced gametogenesis ( Figure 5A) (post hoc S15-S5 df 12, p < 0.01; marginally S15-S10 df 12, p = 0.057 and S30-S5 df 12, p = 0.072). Even more, important changes in the morphology of the oocytes, enlarged nucleus with lower salinities, in 60 and 40% of the gonads at S5 and S10, respectively, were observed probably due to the osmotic stress (post hoc S30-S10 and S15-S10 df 12, p < 0.05; S15-S5 and T30-T5 df 12, p < 0.01) (Figures 6A-D). These important deformations did not translate into changes in the GI; it was similar in all treatments ( Figure 5B).
Gonadal index of V. corrugata increased with salinity although the increment at S30 was masked by the high variability between replicates (Figure 5F). A similar pattern was observed in R. philippinarum with the lowest GI value at S5 because of the gamete resorption. In contrast, GI of C. edule decreased as salinity increased due to massive spawning (Figure 5G).

May
The four species were in the peak of the reproductive season: 80% of R. decussatus had fully mature gametes, R. philippinarum and V. corrugata were mostly spawning and C. edule was mostly spawning and reinitiating gametogenesis toward a respawning (D0, Figure 7).  After 6 days of stress, the response was similar to that in March, with an important effect on R. decussatus (χ 2 = 27.921, df 9, p < 0.0001), V. corrugata (χ 2 = 27.421, df 12, p < 0.01), and C. edule (χ 2 = 48.077, df 12, p < 0.0001). Although not statistically significant, individuals of R. philippinarum at S5 and S10 and, to lesser extent at S15 and S30, tended to be in resorption with important presence of hemocytes and vesicle cells (Figures 7C,D).

July
Ruditapes decussatus, R. philippinarum, and V. corrugata were at the end of the reproductive season since spawning was followed by resorption of the gonads instead of a gonadal restoration (D0, Figures 8A-C,E), while cockles continued in the spawning period as they had reinitiated gametogenesis with a new gonadal restoration (D0, Figure 8G,H).
The most obvious consequence of high temperature on gonads of R. decussatus was the presence of more abnormal oocytes in the post-spawning follicles at the higher temperature (post hoc T20-T37 df 12, p < 0.05; marginally T20-T32 df 12, p = 0.067). At T37, resorption, disorganization of the follicles and severe haemocytic invasion were also observed in more individuals the higher the temperature (χ 2 = 26.468, df 9, p < 0.01).
Although not statistically significant (χ 2 = 3.053, df 3, p = 0.383), the gonads of more than 30% of the C. edule at T20 and T27 and more than 80% at T32 entered reabsorption after 4 days of temperature stress and the 3 individuals that survived at T37 were in resorption. GI reflected this situation with lower values than at D0, even at the lowest temperature. Ruditapes philippinarum was not affected by temperature treatments (χ 2 = 9.1768, df 6, p = 0.3197); although there was a tendency to have more gonads in resorption at the highest temperature was observed (Figures 8C,D).

DISCUSSION
Successful reproduction is a key factor for the survival of species. Reproduction itself implies physiological stress (Domouhtsidou and Dimitriadis, 2001;Garmendia et al., 2010) so any additional exogenous stress during the reproductive period will be dramatic since, to ensure survival, animals have to allocate a large amount of energy to support maintenance costs (Wingfield and Sapolsky, 2003), thus limiting the energy available for growth, storage, and reproduction (Petes et al., 2008;Kooijman, 2010;Stumpp et al., 2012). The effects of stress on reproduction are very complex, varying with time of the year, food supply or height on the shore in the case of intertidal organisms (Petes et al., 2007(Petes et al., , 2008. Our results support those statements because both low salinity and heatwave stress did have important effects on reproduction of the four studied bivalves. Furthermore, the response was species-specific as well as varied with the time of the year, therefore with the annual gametogenic cycle, being more acute in March and May, at the peak of the reproductive season. During sexual resting or early stages of gametogenesis, in the December experiment, a delay in the gametogenic FIGURE 7 | Gonadal stages (percentage) (A,C,E,G) and gonadal index (mean + SD) (B,D,F,H) in the May low salinity experiment in each salinity treatment at the beginning of the experiment (D0) and after six days of stress (D6). Salinity treatments: S5 corresponds to a salinity ramp from 5 to 20, S10 from 10 to 25, S15 from 15 to 30, and S30 control.  development was observed in the clams R. decussatus, R. philippinarum, and Venerupis corrugata at salinities lower than 15 and at salinities lower than 10 in the cockle Cerastoderma edule. This delay lasted three days after the end of the stress only for R. decussatus. Similarly, delayed gametogenesis has been reported in bivalves subjected to thermal (Fearman and Moltschaniwskyj, 2010) and pollution stress (Gagné et al., 2002;Gauthier-Clerc et al., 2002). It is also known that seawater acidification delays reproductive processes of R. decussatus (Range et al., 2011). Surprisingly SFG of this species did not change with the lowest salinity .
Variations of salinity between 5 to 20 and 10 to 25 in March and May, during the peak of the reproductive season, triggered massive spawning followed by gonadal recovery in C. edule, a similar response to that produced by acidification in R. philippinarum during gonadal maturation (Xu et al., 2016). But early spawning of C. edule under stress could lead to a mismatch between the presence of their planktotrophic larvae and the phytoplankton bloom that serves as a food supply, causing potentially starved larvae (Philippart et al., 2003) and thus reducing recruitment success (Beukema et al., 2009).
Another observed response in these experiments was the change in the morphology of the oocytes of R. decussatus in March and May and of V. corrugata and C. edule in May. Oocytes are quite vulnerable to changes in salinity (Gallo et al., 2020) and changes in the morphology could be attributed to changes in the osmotic pressure. Structural changes in oocytes and breakdown processes related to osmotic stress are common in invertebrates (Ram et al., 2004;Simmonds and Barber, 2016). Marked drops in salinity should be followed by an influx of water, which may harm an oocyte unless it can tolerate such changes in volume or, alternatively, can osmoregulate and, thus, prevent marked changes in volume. Whether or not oocytes of these bivalves can osmoregulate, is a question that remains unanswered.
Resorption of gametes with a severe haemocytic infiltration and formation of vesicle cells was the common response of R. philippinarum and V. corrugata at the lowest salinities in March and May. Despite having the same effect on reproduction, only SFG of V. corrugata drastically diminished by low salinities . Even without stress, some species like R. philippinarum present a considerable number of individuals in gonad resorption during the phase of sexual resting and throughout the process of gonadal maturation that evidences a high capacity for gametogenic regeneration and ability to recover the energy invested in the production (Delgado and Pérez-Camacho, 2007). However, stress also causes resorption. To avoid osmotic stress bivalves frequently close their valves, thereby losing between 1 and 24% of the energy acquired (Gosling, 2015). Specifically, these clams do it at salinity 15 and below (Carregosa et al., 2014;Verdelhos et al., 2015;Domínguez et al., 2020) and the cockle at salinities lower than 10, reducing their SFG (Verdelhos et al., 2015;Domínguez et al., 2020) mainly due to a reduction in the feeding activity . Consequently, their reproductive response to stress, both during early gametogenesis with a delay in the gonadal development and during advanced gametogenesis with resorption of gametes, could be attributed to an energetic limitation due to a cessation of feeding since continuous and adequate energy supply is required for gonad ripening (Beninger and Lucas, 1984;Gosling, 2003;Drummond et al., 2006). For instance, mussels experiencing food shortage are not able to maintain ripe gametes in the follicles (Bayne and Thompson, 1970) and seawater acidification forces them to reabsorb the gametes as an energy saving (Bibby et al., 2008).
Our experimental but realistic atmospheric heat waves also triggered gonadal resorption in V. corrugata, R. decussatus, and C. edule which was almost total at the highest temperature.
This response was expected for V. corrugata and C. edule, since after two days of exposure to the same heatwaves, both species suffered a dramatic decrease in their SFG as well as reduction in burrowing activity. On the contrary no effect was expected for R. decussatus since they live burrowed at 13 cm depth within the sediment and depth in sediment buffers the temperature actually experienced (e.g., Johnson, 1965;Harrison and Phizaclea, 1987;Domínguez et al., 2021). The morphology of their oocytes was also modified although such modification suggests other mechanisms rather than osmotic pressure. Several studies found decreasing quality of oocytes in vertebrates and invertebrates with increasing temperature, affecting mainly processes during maturation (Múgica et al., 2015;Gallo et al., 2020), with different underlying mechanisms, e.g., oxidative stress, DNA methylation (Gallo et al., 2020). Contrastingly, R. philippinarum was not affected by temperature treatments other than a tendency toward a larger proportion of gonads in resorption at the highest temperature. Mostly, higher temperature of sea water results in initiation and acceleration of reproductive processes (Lubet, 1959), hence, raising temperature is a common practice in aquaculture to bring bivalves into spawning condition (e.g., Sastry, 1966;Helm and Bourne, 2004;Delgado and Pérez-Camacho, 2007). The most common response of bivalves to thermal stress is the complete release of gametes (Schreck et al., 2001;Petes et al., 2007Petes et al., , 2008 in order to reallocate the energy from reproduction toward defense and repair mechanisms (Múgica et al., 2015) although these spawned gametes can be of poorer quality (Martínez et al., 2000;Ojea et al., 2008). Even more, gonads can remain spawned-out, not able to restore the gonad after spawning, probably due to a lack of energy to regenerate gametes (Petes et al., 2008) caused by drastic reduction of SFG (Domínguez et al., 2021). Alternatively elevated temperatures can restrain reproduction (Heasman et al., 1996;Jeffs et al., 2002), or slow down gametogenesis by the increased demands of metabolism that limits energy for reproduction as described for Mytilus galloprovincialis (Fearman and Moltschaniwskyj, 2010).
Energy balance is fundamental in environmental stress adaptation and tolerance in marine organisms (Pörtner and Farrell, 2008;Sokolova et al., 2012). As discussed, the gonadal cycle of these species may be affected by changes in SFG (Domínguez et al., , 2021. Gametogenesis is an energy demanding process that results in changes in vulnerability, reducing allocation of energy to grow and burrow . In a previous study we calculated that during high tide in a Galician fishing bed, seawater salinity was ≤15 during 9.5% of the time in autumn and during 21 and 23% of the time in winter and spring, respectively . This yields to a loss of >20% of the original SFG value of these four studied species due to exposures to salinities ≤15, that can be even worse in spring when SFG is already low even without stress. Since the probability of being in energy deficit is very high during spring, that could contribute to a slowdown of reproduction in all species but R. philippinarum in March, and in May with dramatic, not only ecological, but also economic consequences. Bivalve fisheries already experience high spatial and temporal variability in catches that are related to high mortality due to strong fluctuations in environmental conditions such as temperature and salinity (Juanes et al., 2012;Parada et al., 2012;Morgan et al., 2013;Aranguren et al., 2014). This fact is aggravated when reductions in reproductive output could lead to a recruitment failure into adult populations (Shanks et al., 2020) especially given the extremely high mortality (>95%) during pelagic larval development and benthic recruitment (e.g., Rumrill, 1990;Gosselin and Qian, 1997). For short-lived species as these bivalves if low salinity episodes in winter and spring are followed by a heat wave in summer, reproduction may be compromised, and the impact magnified if this situation repeats during consecutive years. Long term consequences in the recruitment of these commercially important species are of crucial importance for the sustainability of the fishery.
Negative impacts of the introduced R. philippinarum on other native bivalves, mainly R. decussatus, such as decrease of densities and changes in the distribution, have been claimed as the reason for the decrease of the native clam populations along the European coast (Pranovi et al., 2006;ICES, 2008). However, interspecific competition with other bivalves (Lee, 1996;Byers, 2005), specifically with R. decussatus (Juanes et al., 2012;Bidegain and Juanes, 2013), for space or resource has not yet been clearly demonstrated. According to these last authors, low salinity events, but not competition, modulate the distribution of R. philippinarum and R. decussatus, with R. decussatus being in more freshwater influenced areas since floods may have a stronger effect on the mortality of R. phillippinarum. However, our results point toward a reproductive failure of R. decussatus over time since under high temperature and low salinity stress, reproductive success of R. philippinarum is higher than R. decussatus, mainly during heat waves in summer when R. decussatus is spawning. This could lead to a failure in the recruitment that, sustained over time, and together with overfishing (Massapina and Arrobas, 1991;Maia et al., 2006) could be the reason for the decrease of the populations of this clam. Consequently, the introduced Manila clam with faster maturation, greater number of spawning events, longer and greater reproductive activity, faster growth (Laruelle et al., 1994;Delgado and Pérez-Camacho, 2007;Moura et al., 2017Moura et al., , 2018 and greater tolerance to salinity and temperature stress (Domínguez et al., , 2021 contribute to the observed larger abundance of the introduced clam. Furthermore, R. philippinarum, but not R. decussatus, is typically seeded in shellfish beds (e.g., Royo et al., 2002;Melià and Gatto, 2005;Mantovani et al., 2006) artificially increasing its abundance. In fact, R. philippinarum has become a natural population and one of the most commercially exploited bivalves along the Atlantic and Mediterranean coast of Europe (FAO, 2020) due to the drastic declines of the native carpet shell clam. Currently in Galicia catches of R. philippinarum are 10 times higher than those of the native R. decussatus and five times higher than 20 years ago 3 (Figure 9). If current predictions for the Atlantic coast of Europe, i.e., more frequent heat waves, longer in duration and greater in intensity, as well as increases in the intensity and frequency of extreme precipitation, hold true, R. philippinarum may replace the native R. decussatus in Galicia as has already happened in Arcachon Bay in France and the lagoon of Venice in Italy (Marin et al., 2003;Caill-Milly et al., 2006). The ecological and socioeconomic consequences need to be analyzed to develop adequate management and conservation strategies. This new information is useful in order to promote government intervention to manage these resources sustainably in a context in which the frequency and intensity of extreme events will increase.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because dataset are still used in ongoing research.
Requests to access the datasets should be directed to EV, eotero@uvigo.es.

AUTHOR CONTRIBUTIONS
EV analyzed the histological slides. EV and CO did the statistical analysis. EV and CO wrote the manuscript which was edited and discussed by DW, SW, and LP. All authors designed and ran the experiments.

FUNDING
This research was supported by grants CTM2014-51935-R from the Spanish Ministerio de Economía y Competitividad to the project MARISCO and the Autonomous government Xunta de Galicia-FEDER (project GRC2013-004) and grants NNX11AP77G and 80NSSC20K0074 from the US National Aeronautics and Space Administration (NASA) and OCE1129401 from the US National Science Foundation to DW and SW.