Does an Invasive Bivalve Outperform Its Native Congener in a Heat Wave Scenario? A Laboratory Study Case with Ruditapes decussatus and R. philippinarum

Simple Summary Global climate change is responsible for more frequent heat waves, which offers opportunities for invasive species to expand their range. Two congener bivalves, the native Ruditapes decussatus and the invasive R. philippinarum, were exposed to a heat wave aquaria simulation and analysed for ecological and subcellular biomarkers responses. Despite reduced responses on the ecological level (bioturbation and nutrient concentration), there were differential responses to the heat wave at the subcellular level, where the invasive species seems to be less impacted than the native by the heat wave. This reinforces the common notion that climate change events may provide opportunities for biological invasions. Abstract Global warming and the subsequent increase in the frequency of temperature anomalies are expected to affect marine and estuarine species’ population dynamics, latitudinal distribution, and fitness, allowing non-native opportunistic species to invade and thrive in new geographical areas. Bivalves represent a significant percentage of the benthic biomass in marine ecosystems worldwide, often with commercial interest, while mediating fundamental ecological processes. To understand how these temperature anomalies contribute to the success (or not) of biological invasions, two closely related species, the native Ruditapes decussatus and the introduced R. philippinarum, were exposed to a simulated heat wave. Organisms of both species were exposed to mean summer temperature (~18 °C) for 6 days, followed by 6 days of simulated heat wave conditions (~22 °C). Both species were analysed for key ecological processes such as bioturbation and nutrient generation—which are significant proxies for benthic function and habitat quality—and subcellular biomarkers—oxidative stress and damage, and energetic metabolism. Results showed subcellular responses to heat waves. However, such responses were not expressed at the addressed ecological levels. The subcellular responses to the heat wave in the invasive R. philippinarum pinpoint less damage and higher cellular energy allocation to cope with thermal stress, which may further improve its fitness and thus invasiveness behaviour.


Introduction
Invasive species (IS) are recognized as a worldwide nuisance [1,2] and are often reported as able to occupy niches that were emptied due to abiotic pressures, such as pollutants [3,4], floods, heat waves or droughts [5,6], warming, changes in hydric regimen systematically applied in ecotoxicology. By addressing both lower and higher levels of biological organization responses, researchers can spot some of the present impacts or later trade-offs that may occur in both native and non-native species under different environmental scenarios [23]. Although this question is fundamental to understand the success or failure of non-native invasive species, there is a consistent lack of literature in which both native and non-native closely related species are compared under similar scenario exposures [30].
The Manila clam Ruditapes philippinarum (A. Adams & Reeve, 1850) has its origins in the Indo-Pacific coasts, and due to its worldwide expansion, high fecundity, and growth rates, is featured among the most relevant coastal IS. It lives under a wide range of environmental conditions, between 14 and 42 PSU [23], and 6 to 30 • C [31], as well as in different types of habitat (coastal mid-tidal levels, coastal lagoons, or estuaries) [32]. The Manilla clam was introduced in Portugal for aquaculture and harvesting purposes (first reported in 1984 [33], probably as secondary introduction from Spanish populations [34]), after the decline in native populations of grooved carpet shell R. decussatus (Linnaeus, 1758) and pullet carpet shell Venerupis corrugata (Gmelin, 1791), and cultured oyster stocks [32,33]. Worldwide, in 2014, R. philippinarum accounted for 25% of total mollusc production [31]. While not licensed for aquaculture production in Portugal, the species has been actively dispersed by harvesters throughout several estuaries [33], and is often a dominant species in invasion sites. It is a highly valued species: in 2018, the R. philippinarum official production in Portugal accounted for 862.3 tons/€ 1,501,350.53 [35]. The native congener R. decussatus shares its ecological niche with the invasive R. philippinarum, which overlaps with 40% of the former's trophic niche [32]. Ruditapes decussatus is widely used in bivalve aquaculture and is distributed along the European Atlantic and Mediterranean coasts [36]. It shows a wider salinity tolerance than R. philippinarum, between 7 and 42 PSU [23], but seems to be less sensitive to high temperatures [37]. Both species have been successfully used as biological monitors for coastal contamination and environmental quality [23,38].
Biomarkers are sub-individual endpoints that are recognized for giving potential responses on effects of a given stress, earlier-in-time than the other endpoints measured on the population or community [21]. These may include oxidative stress and damage endpoints, as stress induces the increase of mitochondrial respiration and of reactive oxygen species (ROS), that have the ability to damage lipids, DNA, and proteins [39]. This may lead to irreversible damage and even death. Superoxide dismutase (SOD) and catalase (CAT) are antioxidant enzymes responsible to reduce ROS levels [40]. A more active metabolism will also have energetic needs, and isocitrate and lactate dehydrogenases (IDH and LDH, respectively) are important biomarkers to pinpoint if the organism is favoring aerobic or anaerobic metabolism for energy acquisition [41]. Moreover, at the energetic level, the measurement of the electron transport system (ETS) activity will inform the increase of cellular metabolism and energetic expenditure (energy consumption; Ec), while, by measuring lipid, carbohydrate, and protein contents the energy available is known (Ea) [42]. The ratio between the Ea and the Ec will give the cellular energy allocation (CEA) at a certain period of time, indicative if the stress is conditioning the organisms' energy [42]. These biomarkers, in addition to the mechanistic information given, will also pinpoint potential trade-offs and ecological consequences [43].
The present study aims to predict the effects of a summer heat wave on competing native and non-native invasive congener bivalve species, by comparing their intrinsic fitness and overall activity, based on multilevel responses (bioturbation, nutrient generation, and subcellular biomarkers, including oxidative stress and energetic biomarkers). This will contribute to understanding and comparing the impacts of a heat wave scenario on a native vs invasive bivalve species, addressing: (a) if a summer heat wave will have more severe impacts in the native species ecological engineering behavioral endpoint compared to the invasive species, and/or (b) the heat wave will induce more oxidative stress and energetic expenditure in the native species that may result in trade-offs that will potentially make the invasive outperform the native species.

Species and Aquaria Setup
Both native Ruditapes decussatus and invasive R. philippinarum were supplied by fishermen who harvested them in shellfish banks in Ria de Aveiro (40 • 37 00 N; 8 • 44 27 W), a shallow, mesotidal, coastal lagoon in the north-western coast of Portugal-R. decussatus size was 41.52 ± 0.63 mm, and R. philippinarum 32.41 ± 0.65 mm (individual mean width ± SE), corresponding to 15.70 ± 0.45 g and 8.86 ± 0.16 g (wet weight ± SE), respectively. Biomass in each experimental aquarium (12 × 12 × 35 cm, internal dimensions) was adjusted to be as close as possible by manipulating densities; biomass was calculated after Brey [44], which are close to those found in natural environment [32]. Before the experimental procedures, specimens were acclimated in the laboratory for 7 days, at 17.5 • C and salinity 28, in 6 acclimation tanks (35 × 50 × 28 cm, with~8 cm of sand, overlain with~26 L of seawater diluted with ultrapure water). The acclimation tanks were aerated, and 1/3 of the seawater was replaced every other day to avoid the accumulation of toxic nitrogen compounds in the water. During the acclimation period, specimens were fed with a solution of dried mixed green microalgae (Phytobloom Reef Feed, Necton: Nannochloropsis sp., 45%, Tetraselmis sp., 45% and Isochrysis sp. 10%,~1 µg per g of bivalve wet weight) every other day, after the water replacement. Each experimental aquarium was filled up to 10 cm with sand (medium sand (82.7%) [45]; organic matter loss on ignition (450 • C, 8 h) = 0.18 ± 0.008% mean ± SE), previously air dried for 3 weeks to avoid the introduction of unwanted living macrofauna, and overlain to a total of 30 cm with sand-filtered, UV-treated, seawater (~3 L; adjusted with ultrapure water at salinity 28). Before the beginning of the experimental period, the water was replaced to avoid nutrient pulses associated with assembly [22].
The experiments run for 12 days, within a temperature-controlled room under a 9:15 light:dark regime. Both species were exposed to 2 different temperature treatments: (a) constant temperature (CT; 17.5 • C, which is close to the average water temperature in Ria de Aveiro (17-18 • C) [46][47][48]) and (b) heat wave simulation (HW; 6 days at 17.5 • C, followed by 6 days at 22 • C). Present simulated MHW fits the definition provided by Hobday et al. [16] as a period of at least 5 consecutive days when the temperature is above the 90th percentile. This includes several climatic events that occurred in the Iberian Peninsula in recent decades [13]. The increase in temperature in HW aquaria was achieved with submersible heaters (Jäger 3612, EHEIM, Deizisau, Germany). The HW temperature was achieved within the day. To avoid block effects, the experimental aquaria were randomised using the block.random function of the package psych [49] in R statistical and programming environment [50]. Pumped air for each aquarium ensured aeration and water circulation. Physicochemical parameters (salinity, temperature, pH, and dissolved oxygen), as well as mortality, were recorded every other day. Therefore, the experimental design consisted of 28 aquaria for the following approach: at day 6 (D6), before the temperature elevation, bivalves from 4 aquaria per species (R. decussatus: n = 12; R. philippinarum: n = 16; 8 aquaria in total) were recovered and immediately frozen at −80 • C for subsequent determination of biochemical parameters (8 aquaria in total); at day 12 (D12), at the end of the HW, and integrating the subcellular effects of the previous 6d, individuals from 4 aquaria per species and per temperature treatment (R. decussatus: n = 24; R. philippinarum: n = 32; 16 aquaria in total) were recovered as explained above. Four control aquaria, without bivalves, were also included and kept at CT until the end of the experiment, for procedural control. Individuals were measured and weighed at the beginning and at the end of the experiment. Bivalves were not fed during the experimental period.

Measurement of Particle Reworking
The f-SPI (fluorescent sediment profile imaging [22,26] method was used to assess, non-invasively, the amount of particle reworking. Luminophores (30 g; dyed sediment particles that fluoresce under UV light, 125-250 µm diameter, orange colour; Brian Clegg Ltd., Rochdale, UK) were added to each experimental aquarium, before the introduction of the bivalves. At D12, each aquarium was photographed on all sides under UV light with a Canon EOS 7D reflex digital CMOS camera (18.0 megapixels, set for 0.5 s exposure, f/2.8 diaphragm aperture, ISO 800 film speed equivalent), revealing the luminophores' distribution. Images (JPEG, RGB colour) were cropped to the full internal width of the aquaria (12 cm) and all 4 sides' photographs were merged (48 cm = 5396 pixels, effective resolution = 88.9 µm per pixel). Images were analysed with a custom-made plugin that runs within ImageJ (Version 1.51j8, US National Institute of Health, available at http://imagej.nih.gov/ij/, accessed on 10 February 2017). This plugin converts images to a binary data matrix of the distribution and occurrence of luminophore pixels (0 = background sediment, 1 = luminophore), and returns: (a) mean ( f-SPI L mean , timedependent indication of mixing), (b) median ( f-SPI L median , typical short-term depth of mixing), and (c) maximum ( f-SPI L max , maximum extent of mixing integrated over the long term) mixed depths of particle redistribution, as output. The distance between the highest and lowest points of elevation along the sediment-water interface, (d) surface boundary roughness (SBR = highest − lowest), was also assessed, providing a proxy for surficial faunal activity. Surface boundary roughness, for each aquarium, was averaged from all four sides. These are descriptors of faunal mediated sediment particle reworking, that translates, to some extent, different aspects of invertebrate behaviour [52].

Determination of Biochemical Responses and Data Treatment
Subcellular biomarker analyses were performed in whole-body soft tissues of both native (R. decussatus) and non-native (R. philipinarum) clams. Each individual was homogenised (Ystral, Ballrechten-Dottingen, Germany) in cold potassium-phosphate buffer (0.1 M; pH 7.4) following a 1:10 proportion (m:v). All aliquots were preserved at −80 • C until further analysis.
All biomarker spectrophotometric measurements of each sample were analysed in technical triplicates at 25 • C, using a Synergy H1 Hybrid Multi-Mode Microplate Reader (BioTek Instruments, Winuski, VT, USA). Potassium-phosphate homogenization buffer (0.1 M; pH 7.4) was used as blank in all assays.
To allow the direct comparison of responses between species, the biomarker results for each temperature condition at the end of the exposure (D12), for both species, were normalised as % of change in relation to the D6 biomarker values of the respective species [65], as follows: where B n i,j is the normalized biomarker value for species i and temperature treatment j, B D12 i,j is the biomarker value on D12 for species i and temperature treatment j, and B D6 i is the average value for the biomarker on D6 for species i. Biomarkers results were further integrated under the Integrated Biological Response version 2 (IBR_v2) index, after Sanchez et al. [66]. This index represents the stress levels and overall impact of each experimental condition*species, and in this study, the integrated biomarker responses at D12 (CAT and SOD activities, LPO, DNAd, LDH, IDH, and CEA) were compared to its reference value (D6), retrieving an induction or inhibition value, based on the reference deviation concept. For that, each individual biomarker value was firstly log transformed to reduce variation where Y i,j = individual log-value for species i and temperature treatment j. Then, the general mean of each biomarker (µ) and standard deviation (s i,j ) of Y i,j was calculated, for species i and temperature treatment j, and each Y i,j standardized as For the calculation of a biomarker deviation index, the biomarker data of D6 was also standardized as The biomarker deviation index (A i,j ) was then calculated for each biomarker, for species i and temperature treatment j, as A i,j = Z D12 i,j − Z D6 i . Finally, IBR_v2 was assessed as the sum of all absolute values of biomarker deviation indexes, for each biomarker b:

Statistical Analysis
Independent regression models for the dependent ecological variables (particle reworking: SBR, f-SPI L mean , f-SPI L median , and f-SPI L max ; nutrient concentration: NH 3 -N and PO 4 -P) and the normalized biomarkers values (oxidative stress: SOD activity, CAT activity, DNAd, and LPO; energetic metabolism: IDH and LDH activities; energy allocation: ETS activity, lipid content, carbohydrate content, protein content, total energy available, and cellular energy allocation) against the full factorial combination of independent variables (temperature treatment and species as fixed terms) were tested for statistical significance of their terms. As our focus was to establish the effects of different species, rather than presence versus absence effects, data from the procedural control (aquaria without animals) was removed from the statistical analysis of the ecological effects (bioturbation and nutrients) [22,67,68]. The models were extended to include the appropriate variance covariate structure using a generalized least squares (GLS) estimation procedure [69] (minimal adequate model summaries are appended as Supplementary Material). This extension allows unequal variance along explanatory variables and avoids data transformation when heteroscedasticity is verified. The adequate variance covariate structure was established by comparison with the initial regression model without a variance covariate structure, using a restricted maximum-likelihood (REML) estimation, and based on Akaike Information Criteria (AIC) and visual comparisons of residuals plots. The fixed structure of the minimal adequate model was determined by backward selection based on the maximum-likelihood (ML) ratio test, and re-expressed using REML to obtain the numerical output [70]. The data analyses were carried out within the R statistical and programming environment [50] and the package nlme [71]. Plots were made with ggplot2 package [72]. A significance level of 0.05 was considered in all test procedures.

Environmental Variables
The experimental physicochemical conditions are included in Table 1. The recorded values showed minor variations between treatments. Nevertheless, pH showed an overall increase after D6; dissolved oxygen showed a decrease in those aquaria that were actively heated. Water temperature also increased after D6 in CT aquaria, by thermal influence from the HW aquaria. Despite being randomly distributed, the aquaria were close enough for heat to affect the CT, which were dependent on the room temperature. However, the difference between HW and CT, after D6, was still higher than 4 • C (HW: 22.13 ± 0.15 • C; CT: 18.12 ± 0.05 • C (mean ± SE)). Table 1. Physicochemical conditions of the experimental units for each species (native Ruditapes decussatus and invasive R. philippinarum) and temperature treatments, before and after the introduction of the aquarium heaters at day 6 (mean ± SE). CT: constant temperature; HW: heat wave simulation.

Salinity
Temperature Before D6

Sediment Reworking
The statistical significance and minimal adequate models' summaries are appended as Supplementary Material ("Supplemental material: Statistical significance and minimal adequate model summaries"). The heat wave temperature treatment did not affect any sediment reworking measurement. There were no significant differences between levels in each factor, for surface boundary roughness (SBR) and maximum luminophores depth ( f-SPI L max ) (Supplemental Table S1, Figure 1d). There were differences between species for both mean ( f-SPI L mean : df = 2, L-ratio = 11.074, p = 0.004; Supplemental Material Model S1) and median ( f-SPI L median : df = 2, L-ratio = 12.528, p = 0.002; Model S2) luminophore depths, and, for both measurements, values were higher for R. philippinarum (Figure 1b,c; R. philippinarum f-SPI L mean = 2.13 ± 0.27 and f-SPI L median = 2.08 ± 0.77, vs R. decussatus f-SPI L mean = 0.96 ± 0.17 and f-SPI L median = 0.67 ± 0.13 (cm ± SE)).

Nutrients
For both analysed nutrients, NH 3 -N and PO 4 -P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2).

Subcellular and Biochemical Responses
Statistical differences were found in the relative change of almost all tested subcellular and biochemical responses: these were often due to differences between the two species (statistical significance and minimal adequate models' summaries in Supplementary Material). Yet, temperature treatment affected the variation on some of the energy metabolism-related measurements, in interaction with species identity. Nevertheless, IBR responded only to temperature treatment.

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2).

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2). f different temperature treatments and species (Ruditapes decussatus and R. philippinarum) on (a) NH3entrations (mg L −1 ) in the water at day 12. For clarity, jitter has been applied to the × = argument of void overplotting. constant temperature; heat wave; n = 20.

Subcellular and Biochemical Responses
heat wave; n = 20.

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2).

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differenc between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2)

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2). ) in the water at day 12. For clarity, jitter has been applied to the × = argument of the plot function to avoid overplotting. constant temperature; heat wave; n = 20.

Subcellular and Biochemical Responses
heat wave; n = 20.

Oxidative Stress and Damage
The variation in oxidative stress and related damage biomarkers, from D6 to D12, was significantly affected by species identity, except for superoxide dismutase (SOD: Figure 3a). Catalase (CAT) activity variation was affected only by species (Supplemental Table S2 and Model S3; Figure 3b) and generally with a decrease in activity for R. decussatus (−8.56 ± 3.70% ± SE) and increase for R. philippinarum (10.5 ± 5.29% ± SE). While non-significant, in R. decussatus this decrease in CAT activity was even more evident when under the heat wave (HW) temperature treatment (CT: −2.220 ± 5.534; HW: −16.161 ± 3.707. (% ± SE)). The variation in DNA damage (DNAd) was affected by species identity (Supplemental Table S2 and Model S4; Figure 3c), with very high values in R. decussatus (206.326 ± 35.575% ± SE) and close to zero in R. philippinarum (0.724 ± 22.585% ± SE). Lipid peroxidation (LPO) was again affected by species (Supplemental Table S2 and Model S5; Figure 3d). Yet, despite the interaction between species and temperature treatment being non-significant (d.f. = 2, L-ratio = 3.045, p = 0.081), LPO variation on R. philippinarum was negative (i.e., there was a reduction in LPO from D6 to D12) when under HW (−19.

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2).

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2).

Subcellular and Biochemical Responses
heat wave.

Energy Metabolism Related Enzymes
Regarding energy metabolism related biomarkers, only isocitrate dehydrogenase (IDH) presented significant differences (Supplemental Table S2), with species identity as the driving factor. From D6 to D12, R. decussatus had a decrease in IDH levels (−12.761 ± 3.193% ± SE), while it increased in R. philippinarum (28.479 ± 18.292% ± SE) (Figure 4a; Model S6). The invasive species held lower IDH activity levels when under the HW treatment (Figure 4a).

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2). constant temperature; Figure 1. The effects of different temperature treatments on the sediment reworking by Ruditapes decussatus and R. philippinarum (cm, mean ± SE): (a) Surface boundary roughness (SBR), (b) mean mixed depth of luminophores' redistribution ( f-SPI Lmean), (c) median mixed depth of luminophores' redistribution ( f-SPI Lmedian), and (d) maximum mixed depth of luminophores' redistribution ( f-SPI Lmax), at day 12. For clarity, jitter has been applied to the × = argument of the plot function to avoid overplotting. For comparison, measurements in the absence of animals are presented (grey, control). constant temperature; heat wave; n = 20.

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2).

Subcellular and Biochemical Responses
heat wave.

Integrated Biological Response
The integrated biological response (version 2) (IBR_v2, Figure 6) was significantly affected by temperature (Supplemental Table S2 and Model S12), with higher values found under the HW treatment (CT: 8.543 ± 0.369; HW: 10.146 ± 0.417 (% ± SE)). However, the differences between species and temperature treatments are visible in the integrative star plots for all measured biomarkers (Figure 7), where the contribution of each biomarker for IBR is detailed. Again, the HW treatment affected both species, while at CT

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were between any tested levels at D12 (Supplemental Material: Supplem

Nutrients
For both analysed nutrients, NH3-N and PO4-P, there were no significant differences between any tested levels at D12 (Supplemental Material: Supplemental Table S1, Figure 2).

Subcellular and Biochemical Responses
heat wave.

Integrated Biological Response
The integrated biological response (version 2) (IBR_v2, Figure 6) was significantly affected by temperature (Supplemental Table S2 and Model S12), with higher values found under the HW treatment (CT: 8.543 ± 0.369; HW: 10.146 ± 0.417 (% ± SE)). However, the differences between species and temperature treatments are visible in the integrative star plots for all measured biomarkers (Figure 7), where the contribution of each biomarker for IBR is detailed. Again, the HW treatment affected both species, while at CT the normalized values are close to zero. Yet, the HW was responsible for contrasting changes in different biomarkers according to the species. Under the HW, in the native R. decussatus (Figure 7a) a reduction on CEA and on CAT and LDH activities, along with an increase in oxidative damage (LPO and DNAd), are the main responses influencing the IBR index. In opposition, for the invasive R. philippinarum under HW (Figure 7b), the driving responses are a reduction on LPO and DNA damage levels along with an increase on CEA. These observations are in line with the results of the statistical analysis of individual biomarkers. changes in different biomarkers according to the species. Under the HW, in the native R. decussatus (Figure 7a) a reduction on CEA and on CAT and LDH activities, along with an increase in oxidative damage (LPO and DNAd), are the main responses influencing the IBR index. In opposition, for the invasive R. philippinarum under HW (Figure 7b), the driving responses are a reduction on LPO and DNA damage levels along with an increase on CEA. These observations are in line with the results of the statistical analysis of individual biomarkers.  R. philippinarum (n = 32) at constant temperature and after a heat wave simulation: SOD-superoxide dismutase activity; CAT-catalase activity; DNAd-DNA damage; LPO-lipid peroxidation levels; IDH-isocitrate dehydrogenase activity; LDH-lactate dehydrogenase activity; and CEA-cellular energy allocation (which incorporates data on lipid, protein, and carbohydrate contents and electron transport system levels).

Discussion
While both species differed in the intensity and depth of particle reworking, this ecological process did not suffer a differential impact as a consequence of the heat wave (HW). Differences in the luminophore depths between both species may be partially explained by body size differences. Nevertheless, the larger R. decussatus showed smaller

Nutrients
For both an between any teste

Nutrients
For both analysed nutrients, NH3-N and PO between any tested levels at D12 (Supplemental M changes in different biomarkers according to the species. Under the HW, in the native R. decussatus (Figure 7a) a reduction on CEA and on CAT and LDH activities, along with an increase in oxidative damage (LPO and DNAd), are the main responses influencing the IBR index. In opposition, for the invasive R. philippinarum under HW (Figure 7b), the driving responses are a reduction on LPO and DNA damage levels along with an increase on CEA. These observations are in line with the results of the statistical analysis of individual biomarkers.  R. philippinarum (n = 32) at constant temperature and after a heat wave simulation: SOD-superoxide dismutase activity; CAT-catalase activity; DNAd-DNA damage; LPO-lipid peroxidation levels; IDH-isocitrate dehydrogenase activity; LDH-lactate dehydrogenase activity; and CEA-cellular energy allocation (which incorporates data on lipid, protein, and carbohydrate contents and electron transport system levels).

Discussion
While both species differed in the intensity and depth of particle reworking, this ecological process did not suffer a differential impact as a consequence of the heat wave (HW). Differences in the luminophore depths between both species may be partially explained by body size differences. Nevertheless, the larger R. decussatus showed smaller Figure 7. Star plots integrating all biomarkers measured in (a) Ruditapes decussatus (n = 24) and (b) R. philippinarum (n = 32) at constant temperature and after a heat wave simulation: SODsuperoxide dismutase activity; CAT-catalase activity; DNAd-DNA damage; LPO-lipid peroxidation levels; IDH-isocitrate dehydrogenase activity; LDH-lactate dehydrogenase activity; and CEA-cellular energy allocation (which incorporates data on lipid, protein, and carbohydrate contents and electron transport system levels).

Discussion
While both species differed in the intensity and depth of particle reworking, this ecological process did not suffer a differential impact as a consequence of the heat wave (HW). Differences in the luminophore depths between both species may be partially explained by body size differences. Nevertheless, the larger R. decussatus showed smaller mean and median luminophore depths. It should be noted that, however, R. philippinarum held higher densities, which could have contributed to a higher level of intraspecific competition and subsequent avoidance behaviour. With R. philippinarum featuring shorter siphons [73], the number of upward and downward movements within the sediment might also have been higher than those of R. decussatus. Adding to this, smaller individuals may show higher metabolic requirements [74], and therefore an increased active search for food. These concur to increase mean and median sediment reworking depths, which are proxies for the intensity of reworking [52]. A similar tendency was verified in another invasive bivalve, Corbicula fluminea, with smaller individuals responsible for higher bioturbation levels [22]. Surprisingly, maximum luminophore depths were similar on both Ruditapes species, which means that, despite size, density, or siphon length, in the long term both species might reach similar levels of particle reworking, and therefore, differ mostly in the rate of this process. The phylogenetic/taxonomic proximity of both species was expressed clearly in the behavioral response to thermal stress, which is absent in both cases, and therefore not conveyed into ecosystem functioning, under the specific set of conditions under analysis. The results of nutrient content in the water column reinforces this lack of expression at the ecological level, as responses showed no statistical significance, even if densities were slightly biased towards R. philippinarum.
Despite the absence of response on both species to the HW treatment at the ecological and behavioral levels, there were significant responses to thermal stress at the subcellular levels for both species. The analysis of subcellular biomarkers can act as an early warning tool for later-in-time responses in higher level of biological organization, which may occur after the stress exposure stage [21], and not depicted at this stage. Addressing the used ecological endpoints later-in-time, along with others such as growth and reproduction, would provide critical additional information about this potential link of biomarker responses with ecological processes and populational dynamics, and the expected trade-off that give advantage to a certain species [22], and that give relevance to the early-in-time responses. The use of normalized values, with biomarker levels from D6 as baseline (prior to the heat wave simulation), proved to be a useful approach in this context and allowed comparison of the relative changes in response to stress between species. Contrasting biomarker responses were observed, with lower impacts in the invasive species under the HW. The native R. decussatus showed a higher level of oxidative damage, as consequence of a compromised antioxidant and energetic investment given by the reduction on catalase and lactate dehydrogenase activities, while also showing a reduction in their cellular energy allocation. On the other hand, the invasive R. philippinarum seems to benefit from the increase in the temperature by showing a reduction in damage (lower lipid peroxidation and DNA damage), and a reinforcement of cellular energy allocation. The increase in the cellular energy allocation index, which relates energetic reserves to their consumption by cellular metabolism [42], suggests that the invasive species may have a higher energy budget for other functions, such as growth, reproduction, and locomotion/burrowing activity. Due to the defined nature and length of the experiment, growth and reproduction were not assessed but one cannot exclude that these endpoints could be reactive at a later period and translate the trade-offs better than the present studied higher level of organization proxies (behaviour). The increase in cellular energy allocation in R. philippinarum is justified mainly by the increase in lipid content, together with a lesser energetic expenditure, as electron transport system values reflect. One must add that, despite no fresh food being added to the water during the experiment, the sediment was not depleted of organic matter and detritus on which both species may have fed on to alter energy reserves [32]. One may also not disregard that the unexpected increase in lipids may also be related to an eventual r-strategy from the invasive species, with stress triggering gonad maturations and oocyte formation and the biosynthesis of lipids [75][76][77], these lipids being the main reserve of these structures [78,79]-while the native species relied on a k-strategy, as previously reported for other native and invasive bivalves under stress [5]. After integrating all these responses into the integrated biological response index, it becomes clearer that the HW treatment induced alterations in the biochemical stress responses for both species. However, the biomarker responses contributing to the higher index in both species are opposite, indicating a negative impact to the native species (elevated damage levels and lower energetic investments and cellular energy allocation) while those in the IS indicate that the individuals are coping better with the challenge (less damage and higher cellular energy allocation). In sum, these results suggest that, under the same HW conditions, the IS will possess a better condition, which might provide the marginal gains that makes non-native species better contenders in interspecific competition.
There is reduced experimental evidence in the literature that physiological tolerance to multiple environmental conditions is an obligatory characteristic of IS [4,19,80]. However, Zerebecki and Sorte [19] verified a positive link between the geographic range (and the maximum habitat temperature) and the temperature tolerance of epibenthic communities in California (USA), and simultaneously found that invasive species withstood higher lethal temperatures. In the same study, the authors found that the expression of a heatshock protein was higher in the invasive tunicate Diplosoma than in the native Distaplia. For bivalves, a similar tendency was found in two Mytilus species: the invasive M. galloprovincialis showed higher expression of heat-shock proteins, with a higher induction temperature threshold, than the native M. trossulus [81]. Another freshwater bivalve study showed similar results, with the invasive Sinanodonta woodiana dealing better with higher temperatures and contamination levels than the native Anodonta anatina [4]. Sorte et al. [7] also analysed the differential survival between native and non-native species in a marine fouling community, and found that non-native had higher LT 50 values (~2 • C average), which could indicate that these species can cope with higher temperature levels and show sublethal responses later during a thermal anomaly. Simultaneously, the same authors predicted abundance and growth rates increases for most tested IS, under a +4.5 • C temperature increase scenario. All of these are examples of how temperature mediates the survival and success of non-native IS. However, the idea that IS are favored by eury-tolerance, as the capacity to withstand wider ranges of environmental stressors [30], could be challenged. The concept of eury-tolerance includes, for example, tolerance to wide ranges of salinity, temperature, or contaminants [4,30,82], offering competitive advantages towards their native counterparts, particularly in transitional aquatic systems such as estuaries. Nevertheless, the hypothesis of eury-tolerance in IS mostly assumed due to wide geographical ranges as a surrogate measure for extended ranges of temperature or salinity (or other environmental variables). Despite conflicting reports on the tolerance range of IS, global climate change is often reported as an invasion facilitator [1,83], especially because species from warmer regions can expand their geographic ranges to higher latitudes [7], where they can now find more suitable conditions than before. For instance, the native R. decussatus showed a wider tolerance range for salinity than both the invasive R. philippinarum and the native Venerupis corrugata [23]. Other examples are described by McMahon [80] who reviewed the physiological adaptations in aquatic invasive invertebrates and found that the high-profile bivalve invaders Corbicula fluminea and Dreissena polymorpha showed reduced resistance and capacity adaptations (with massive mortality rates under extreme environmental conditions) when compared with endemic native unionoideans. This author sustains that the species invasive prowess is mostly related to r-selected traits (e.g., fast population growth, early maturity, and high fecundity) while native species must rely on their resistance capacity to warrant maintenance under natural disturbance, due to the k-selected traits. Ferreira-Rodríguez and colleagues [5] reached similar conclusions for the invasive C. fluminea vs the native Unio delphinus. A similar trend was found within a group of freshwater pecarid crustaceans [84]. In the present study, results show that R. philippinarum was less affected than the native R. decussatus. This means that, even without proving the extended tolerance range to temperature, i.e., eury-thermality [30], one can hypothesize that the top limit of the thermal range of the IS R. philippinarum is higher than that of its native counterpart R. decussatus.
One of the implications of the results here addressed is that the current climate crisis, not only due to the increase in average temperature, but especially due to the increase on the number of heat waves [6,17], may favor this particular IS, which is already a concern due to human-mediated dispersion as valued species [3,8,33]. It is unclear if the first decline in native species was due to the Manilla clam introduction in the Portuguese territory-which occurred almost simultaneously with the first signs of that decline-or if other factors, such as overfishing, diseases, or habitat modification, might have occurred. In this context, nevertheless, results fit within the narrative that IS are indeed aided by the interactive effects of global changes-in this case, climate change and the human-mediated introduction of species-contributing to the homogenization of the biota [7,9]. This means that the synergy between climate change and human mediated dispersion contributes to R. philippinarum overriding the filters (dispersion and habitat suitability) that often control the establishment of introduced species [30,85,86], with an increase in the pressure on local communities [2,7]. Extreme temperature fluctuation is often accompanied by changes in dissolved oxygen and salinity, which interactions might have disproportionate effects on the bivalve's metabolism and fitness and limit the ability to describe the effects of global changes solely based on thermal characterization. Because changes in ecosystem function due to an introduced species depend not only on the rate and intensity of ecosystem processes, but also on the interactions of the introduced species with the invaded community and the environment [25], there is a high level of unpredictability on the ecological consequences of such introduction. However, Lopes et al. [8], when comparing the same species, verified that the replacement of the native by the IS was not expressed as differences on bioturbation. This pair of species is also known to hybridize [31], so there are hints that the introduction and dispersion of this particular IS may drive small changes in ecosystem functioning, due to their phylogenetic proximity, supported by the concept of functional redundancy [1,68]. However, these genetic outcomes should be monitored for even more severe and unpredictable impacts. If the IS indeed favored by climate change, compositional changes in the invaded communities may also occur due to a disproportionate increase and consequential dominance that is often reported for several IS [2,7,19].

Conclusions
The present work intentionally studies a conservative heatwave, with a relative short period and little temperature range. However, considering that not only the number of HW is increasing, but also the duration of these phenomena, this also triggers the question as to whether longer, repeated, or more intense HW will convey different and pronounced results at the ecological and behavioural level among these species, measurable during the experiment timescale, or even if effects extend after the end of the HW. However, despite the reduced range that was evaluated here, the results from this work contribute to better enlighten the processes that might offer IS higher fitness and competitive advantages towards their native counterparts, supported by biomarker analysis as a set of early responding endpoints that represent potential trade-offs for those advantages later, even before the transposition into ecosystem processes and functioning. Present results indicate that the intensity of responses is not consistent at the different levels under analysis, which means that assumptions based on a single endpoint, and limited timeframe, may not provide enough information on the consequences of environmental stressors on the overall fitness of organisms, as usually the link between different levels of biological organization will emerge over time. Here, lower levels of biological organization, such as biomarkers, provide putative early-warning signs that should add a cautionary tale to the scenario, even when other individual level endpoints fail to depict differential effects. This underlines the relevance of integrative, multilevel studies to disclose the impacts of stressors, and particularly, in the scope of this study, thermal stress, in a context of competing species that share similar ecological niches.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/biology10121284/s1, Table S1: Summary of significant terms found in the generalised least squares models, Table S2: Summary of significant terms found in the generalised least squares models, Model S1: Mean Luminophore Depth ( f-SPI L mean ), Model S2: Median Luminophore Depth Funding: This study had the support of the Fundação para a Ciência e a Tecnologia (FCT) Strategic Projects UID/MAR/04292/2020 and UIDB/04004/2020 granted to MARE and CFE respectively (through national funds; PIDDAC), the grant awarded to Lénia Rato (SFRH/BD/138492/2018) and the contracts attributed to Sara Leston, Filipe Martinho, and Sara Novais in the scope of the Decree-Law 57/2016. Further support was provided by FCT through project MARINE INVADERS-The impact and mechanisms of success of the invasive seaweed Asparagopsis armata on coastal environments (POCI-01-0145-FEDER-031144). The project was also partly funded by the Integrated Programmes of SR&TD "SmartBioR" (reference Centro-01-0145-FEDER-000018) and ReNature (CENTRO-01-0145-FEDER-000007) co-funded by Centro 2020 program, Portugal2020, European Union, through the European Regional Development Fund, and project "Global Invaders-Global trends among valued aquatic invertebrate species: competitive advantages across different latitudes" funded by FCT and DAAD-Deutscher Akademischer Austauschdienst.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.