The Combined Effects of Toxic Microcystis aeruginosa and Thermal Stress on the Edible Clam (Corbicula fluminea): Insights into Oxidative Stress Responses and Molecular Networks

Cyanobacterial blooms (CYBs) have become a global environmental issue, posing risks to edible bivalves. Toxic cyanobacteria and thermal stress represent the two key co-occurring stressors to bivalves experiencing CYBs. To investigate the combined effects of these stressors on the edible bivalve Corbicula fluminea, the responses to oxidative stress and the molecular mechanisms of physiological adaptations in C. fluminea were examined under co-exposure to toxic Microcystis aeruginosa and thermal stress. The activity of antioxidant enzymes, including GST, SOD, CAT, GPx and GR, was significantly influenced by the interaction between temperature and M. aeruginosa (p < 0.05). A positive correlation was observed between toxic M. aeruginosa exposure and elevated SOD and GPx activities at 30 °C, demonstrating that SOD and GPx may help C. fluminea defend effectively against MCs under thermal stress. Furthermore, significant interactive effects between toxic M. aeruginosa and temperature were also observed in ROS and MDA (p < 0.05). The results of the PCA and IBR index also evidenced the apparent influence of toxic M. aeruginosa and thermal stress on oxidative stress responses of C. fluminea. The eggNOG and GO annotations confirmed that a substantial portion of differentially expressed genes (DEGs) exhibited associations with responses to oxidative stress and transporter activity. Additionally, KEGG analysis revealed that abundant DEGs were involved in pathways related to inflammatory responses, immune functions and metabolic functions. These findings improve our understanding of the mechanism of the physiological adaptation in bivalves in response to cyanotoxins under thermal conditions, potentially enabling the evaluation of the viability of using bivalves as a bioremediation tool to manage CYBs in eutrophic waters.


Introduction
Due to elevated eutrophication from intensified human activities (e.g., agricultural activities, urbanization) and global warming, cyanobacterial blooms (CYBs) have been increasing in frequency, intensity and duration worldwide [1].One of the primary risks associated with CYBs is the synthesis of various cyanotoxins.Microcystins (MCs), known as one of the most toxic cyanotoxins produced by Microcystis strains, have been the subject of numerous studies and have attracted extensive attention due to their high toxicity and widespread occurrence [2].MCs are generally recognized as tumor promoters, with chronic exposure to MCs associated with the appearance of liver cancer, hence posing a serious risk to human health [3].As typical hepatotoxins, MCs can directly affect the liver or hepatopancreas of aquatic organisms, causing oxidative stress, liver damage and even reducing growth and reproductive capacity [4].
As global warming increases, many eutrophic waters are experiencing significant fluctuations in temperature.Numerous studies have indicated that elevated temperatures have the potential to provide advantages to toxic bloom-forming cyanobacteria by increasing their growth rate and expanding their ecological niche [5][6][7].Therefore, the observed rise in temperatures has been associated with an increased proliferation and a prolonged duration of toxic CYBs [8].For aquatic animals, elevated temperatures could enhance their physiological metabolism's performance and increase oxygen consumption, thus potentially leading to an increase in reactive oxygen species (ROS) levels and the induction of oxidative stress [9].Although the individual effects of toxic CYBs and thermal stress on aquatic animals have been well characterized, there is a lack of information on the combined effects of toxic CYBs and elevated temperatures.Previous studies have assessed the toxicity of cyanobacteria at different temperatures in aquatic animals, including zebrafish (Danio rerio), mussels (Hyriopsis cumingii), freshwater daphnids (Moina macrocopa) and rotifers (Brachionus calyciflorus), indicating that the combined effects of temperature and cyanotoxins are complicated and species-specific [9][10][11][12][13].Consequently, a comprehensive understanding of the combined impacts of toxic cyanobacteria and thermal stress on various aquatic organisms will provide a broader perspective on the sustainability of aquatic ecosystems.
Bivalves generally inhabit the benthic zone and play an important role in the trophic structure of freshwater ecosystems, which can effectively facilitate energy flows and nutrient cycling within the ecosystem [14].Compared to other aquatic organisms, bivalves exhibit more susceptibility to the effects of CYBs due to their ability to filter feed on cyanobacteria as well as their restricted migration capacity [15].The edible bivalve Corbicula fluminea, formerly endemic to Eastern Asia and Southeast Asia, has achieved global distribution as a result of its high fertility and adaptability to various habitats [16].C. fluminea has significant economic value, not only for direct human consumption, but also as a feed option for livestock and poultry, as well as a natural bait for carnivorous fish and waterfowl [17].As a model bivalve, C. fluminea has been widely used for the evaluation of the toxic effects of various pollutants due to its low mobility, ease of cultivation, economic practicality and wide distribution [18][19][20].Furthermore, C. fluminea has also been considered a valuable tool for the restoration of eutrophic lakes where toxic CYBs occur due to their efficient filtration of bloom-forming cyanobacteria and their ability to tolerate cyanotoxins [21,22].Previous studies have demonstrated that toxic M. aeruginosa or thermal stress could affect the growth performance, biochemical parameters and physiological functions of C. fluminea [23][24][25].However, little is known about the combined effects of toxic M. aeruginosa and thermal stress on this species.Furthermore, there is still a lack of information on the molecular mechanisms of tolerance to MCs in C. fluminea at high temperatures.
Most of the evidence indicates that the adverse effects of MC exposure and thermal stress on aquatic animals could be attributed to ROS production [26].The presence of an excessive amount of ROS has the potential to cause lipid peroxidation in cell membranes, protein denaturation and DNA damage [9].Cell antioxidant defense systems, including a series of antioxidant enzymes such as superoxide dismutase (SOD), catalase (CAT), glutathione peroxidase (GPx) and glutathione reductase (GR), possess the ability to maintain redox equilibrium and protect cells against oxidative damage [25].Oxidative stress arises as a result of an imbalance between ROS-induced oxidative effects and ROS elimination mediated by antioxidant defense systems.Therefore, these biochemical parameters are commonly employed as biomarkers in the evaluation of the biological toxicity responses of MCs in bivalves.Additionally, the deployment of high-throughput RNA sequencing technology (RNA-seq) has enabled the comprehensive analysis of transcriptomes, resulting in their increasing use in the study of the molecular regulatory mechanisms of aquatic animals in response to MC exposure [27,28].
This study investigated the combined effects of toxic M. aeruginosa and thermal stress on oxidative stress responses, including GST, SOD, CAT, GPx and GR activity, ROS level and malondialdehyde (MDA) content of the edible clam C. fluminea.Moreover, a principal component analysis (PCA) and an integrated biomarker response (IBR) analysis were performed to highlight the integrated variations in biochemical parameters between different treatment groups (three different concentrations of toxic M. aeruginosa × three temperature gradients).Furthermore, the RNA-seq was conducted in the hepatopancreas to determine the molecular response in C. fluminea exposed to toxic M. aeruginosa under thermal stress.The present study could contribute to further understanding of the physiological adaptation mechanisms of edible bivalves to co-exposure to toxic cyanobacteria and thermal stress, potentially enabling the evaluation of the feasibility of bivalves as a bioremediation tool to control CYBs in eutrophic waters.

Experimental Clam Acclimation and Algae Culture
Healthy C. fluminea was obtained from Hongze Lake (China, 118.8 • E, 33.3 • N).Clams (length 2.16 ± 0.23 cm) were transported to the laboratory and acclimated for at least 14 days under controlled conditions, with a water temperature of 20 ± 1 • C, dissolved oxygen content greater than 6.0 mg/L, a pH of 7.4 ± 0.3 and a 12 h light-dark photoperiod.The clams were fed Chlorella vulgaris (1 × 10 5 cells/mL) twice daily during the acclimatization period.The algae species, namely toxic M. aeruginosa (FACHB-905) and C. vulgaris (FACHB-8), were obtained from the Institute of Hydrobiology, CAS (Wuhan, China).The microalga was cultured in BG11 medium and harvested in stationary phase with a density of 0.5-1 × 10 7 cells/mL.

Experimental Design
The acclimated clams were exposed to nine treatments (3 × 3 factorial design), including three temperature levels (20, 25 and 30 • C) and three toxic M. aeruginosa concentrations (0, 1 and 10 × 10 5 cells/mL).The treatment groups without M. aeruginosa exposure (0 cells/mL) were fed with C. vulgaris (1 × 10 5 cells/mL) instead.Forty clams were randomly assigned to each treatment group and cultured in 20 L tanks.Each treatment was performed in triplicate.During the experiment, full water renewal was performed once daily with the addition of microalgae.
Five clams from each tank were randomly collected for biochemical parameter determinations on days 1, 3, 5 and 7.After 7 days of exposure, ten clams from the control group (0 cells/mL toxic M. aeruginosa × 20 • C) and the treatment group (10 × 10 5 cells/mL M. aeruginosa × 30 • C) were collected for transcriptome analysis.

Determination of Biochemical Parameters
The hepatopancreas of five sampled clams, weighing about 0.2 g, were homogenized via sonication in the ice bath.The supernatant was produced via centrifugation (2500 rpm for 10 min at 4 • C) and then used for the determination of biochemical parameters.Biochemical parameters were evaluated using assay kits (Nanjing Jiancheng Institute of Bioengineering, Nanjing, China) as per the manufacturer's protocols.In summary, the determination of GST activity was conducted using the colorimetric approach, employing the GST Assay Kit.The measurement of SOD activity was conducted via the SOD Assay Kit using the WST-1 method.The CAT Assay Kit was utilized to evaluate the CAT activity by measuring the rate of the hydrogen peroxide (H 2 O 2 ) decomposition reaction.The activity of GPx was evaluated using the GSH-PX Assay Kit, which quantifies the rate of the catalytic reaction involving glutathione.The assessment of GR activity was performed by quantifying the reduction in absorbance at 340 nm, which is indicative of NADPH consumption.Furthermore, the determination of MDA content was performed via the thiobarbituric acid method, employing the MDA Assay Kit.The ROS Assay Kit was employed to quantify ROS levels using the DCFH-DA probe method.The protein content was evaluated via the Protein Assay Kit, applying the Coomassie Brilliant Blue method.

The IBR Index Analysis
The integrated biomarker response index was calculated from the biochemical variables investigated.Calculation of IBR values was carried out using the formulae proposed by Sanchez et al. [29]: The variables X i and X 0 represent the individual biomarker data and the mean reference data, respectively.Y i denotes the logistic transformation of the individual biomarker's data.µ and σ denote the average and standard deviation of Y i .Furthermore, Z i represents the average of the standardized biomarker response, while Z 0 represents the average of the reference biomarker data.Finally, A denotes the deviation in reference values for each biomarker.

Transcriptome Analysis
The hepatopancreas of sampled clams were dissected for RNA extraction.Libraries were prepared using the VAHTS ® Universal V6 RNA-Seq Library Prep Kit and then subjected to sequencing on the Illumina NovaSeq 6000 platform.The differential expression analysis was conducted using the DESeq2 package developed by Anders and Huber [30].To identify genes that exhibit differential expression, a fold change (FC) > 2 and a false discovery rate (FDR) < 0.05 were used.Functional annotations and enrichment analyses of differentially expressed genes (DEGs) were conducted using the evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG), the Gene Ontology (GO) and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) databases.

Statistical Analyses
Statistical analysis was performed using SPSS 25 software (IBM, New York, NY, USA).Data for the biochemical parameters determined were expressed as mean ± SD.Statistical tests of normality (Shapiro-Wilk's test, 1% risk) and equal variance (Levene's test, 5% risk) were employed for the results prior to analysis.The pair correlation coefficients across biochemical parameters calculated by Pearson correlation were shown in Figure S1.A three-way analysis of variance (ANOVA) was used to analyze the effects of temperature, toxic M. aeruginosa and time course on variations of biochemical parameters, followed by Tukey's test to identify pairwise differences.The level of significance was defined as p < 0.05.Principal component analysis of these parameters was performed using GraphPad Prism 9.3 software (GraphPad Software, Boston, MA, USA).GST is responsible for facilitating the conjugation of electrophilic substrates with glutathione (GSH), thus serving a crucial function in detoxification and defending aquatic organisms from the harmful effects of oxidative metabolic by-products [31].At low and medium temperatures (20 and 25 • C), the treated clams exhibited a gradual increase in GST activity as the levels of toxic M. aeruginosa increased, regardless of the duration of exposure (Figure 1).No statistically significant differences were observed between different diet treatments at 20 • C (p > 0.05).However, under high temperature conditions (30 • C), GST activity in the treatment groups that were fed low concentrations (1 × 10 5 cells/mL) of toxic M. aeruginosa was significantly increased compared to those fed C. vulgaris at sampling times on days 5 and 7, indicating that C. fluminea could exhibit a higher metabolic rate towards low concentrations of M. aeruginosa at high temperatures (30 • C).The three-way ANOVA results showed that temperature and M. aeruginosa significantly affected the GST activity of C. fluminea (p < 0.05).Furthermore, significant interactions between temperature and M. aeruginosa for GST activity were observed in treated clams (p < 0.05).This may indicate that GST-mediated metabolism enables C. fluminea to endure cyanobacterial blooms but may be affected by thermal stress.In H. cumingii, similar elevated GST activity was observed and GST activity was obviously affected by temperature and toxin concentration [9].The comparable phenomenon of increased GST activity under MC exposure can also be observed in marine bivalves, such as M. edulis and C. gigas [32].Burmester et al. [33] found that the GST of invasive bivalves could respond more strongly to cyanotoxins than that of indigenous species.Consequently, the physiological adaptations to CYBs observed in C. fluminea could be attributed to the increased activities of biotransformation enzymes in response to cyanotoxins under thermal stress.

Results and Discussion
medium temperatures (20 and 25 °C), the treated clams exhibited a gradual increase in GST activity as the levels of toxic M. aeruginosa increased, regardless of the duration of exposure (Figure 1).No statistically significant differences were observed between different diet treatments at 20 °C (p > 0.05).However, under high temperature conditions (30 °C), GST activity in the treatment groups that were fed low concentrations (1 × 10 5 cells/mL) of toxic M. aeruginosa was significantly increased compared to those fed C. vulgaris at sampling times on days 5 and 7, indicating that C. fluminea could exhibit a higher metabolic rate towards low concentrations of M. aeruginosa at high temperatures (30 °C).The three-way ANOVA results showed that temperature and M. aeruginosa significantly affected the GST activity of C. fluminea (p < 0.05).Furthermore, significant interactions between temperature and M. aeruginosa for GST activity were observed in treated clams (p < 0.05).This may indicate that GST-mediated metabolism enables C. fluminea to endure cyanobacterial blooms but may be affected by thermal stress.In H. cumingii, similar elevated GST activity was observed and GST activity was obviously affected by temperature and toxin concentration [9].The comparable phenomenon of increased GST activity under MC exposure can also be observed in marine bivalves, such as M. edulis and C. gigas [32].Burmester et al. [33] found that the GST of invasive bivalves could respond more strongly to cyanotoxins than that of indigenous species.Consequently, the physiological adaptations to CYBs observed in C. fluminea could be attributed to the increased activities of biotransformation enzymes in response to cyanotoxins under thermal stress.

Influence on Antioxidant Defense Systems
SOD and CAT generally serve as the primary defense mechanisms in bivalves against oxidative stress induced by exogenous toxic substances [34].SOD can facilitate the production of H2O2 through the catalysis of O2 − , which subsequently involves its degradation by CAT into water and O2 [35].Therefore, the transformation of excess ROS (such as O2 − and H2O2) into harmless metabolites can be linked to elevated SOD and CAT activities.As the second line of defense against oxidative stress, GPx could facilitate the conversion

Influence on Antioxidant Defense Systems
SOD and CAT generally serve as the primary defense mechanisms in bivalves against oxidative stress induced by exogenous toxic substances [34].SOD can facilitate the production of H 2 O 2 through the catalysis of O 2 − , which subsequently involves its degradation by CAT into water and O 2 [35].Therefore, the transformation of excess ROS (such as O 2 − and H 2 O 2 ) into harmless metabolites can be linked to elevated SOD and CAT activities.As the second line of defense against oxidative stress, GPx could facilitate the conversion of reduced GSH to GSH disulfide, thus removing excess H 2 O 2 [32].In addition, GR plays a crucial role in catalyzing the reduction of GSH disulfide to GSH [36].In the present study, no significant differences were observed for SOD activity between different diet treatment groups at 20 • C (Figure 2).However, at conditions of 25 and 30 • C, a significant increase in SOD activity was observed in treated clams fed toxic M. aeruginosa (1 and 10 × 10 5 cells/mL) compared to those fed C. vulgaris at each sampling time.On the contrary, compared to the treatment groups fed C. vulgaris, CAT activities increased significantly in the treatment groups which were fed high concentrations of M. aeruginosa (10 × 10 5 cells/mL) at 20 and 25 • C, while no significant differences were observed between the different diet treatments at 30 • C (Figure 2).GPx activity exhibited similar fluctuations to SOD activity during the co-exposure period to M. aeruginosa and thermal stress (Figure 2).GPx activity gradually increased as the concentration of toxic M. aeruginosa increased at 25 • C and 30 • C, with significantly higher GPx activities observed in treatment groups fed high concentrations of toxic M. aeruginosa compared to those fed C. vulgaris at elevated temperatures and at each sampling time point (p < 0.05).The results demonstrated a positive correlation between toxic M. aeruginosa exposure and elevated SOD and GPx activities at high temperatures (30 • C), indicating that the SOD and GPx enzymes play a crucial role in antioxidant defense systems and may aid C. fluminea in effectively defending against MCs under thermal stress.Additionally, GR activity displayed a significant decrease in treated clams fed toxic M. aeruginosa compared to those fed C. vulgaris at 25 • C and 30 • C (p < 0.05), while this significant inhibition of GR did not appear when exposed to toxic M. aeruginosa at 20 • C (Figure 2).Results of the three-way ANOVA demonstrated that toxic M. aeruginosa and temperature have significant effects on the activities of SOD, CAT, GPX and GR (p < 0.05) (Table 1).Significant interactive effects between toxic M. aeruginosa and temperature were also observed on these tested parameters (p < 0.05) (Table 1).This is in accordance with the findings of Kim et al. [32,36], who showed that microcystin exposure induces the enzymatic activities of antioxidant defense systems including GST, CAT, SOD, GPx and GR in the bivalves Crassostrea gigas and Mytilus edulis.A previous study in C. fluminea also demonstrated that SOD activity was elevated by toxic M. aeruginosa [25].Our results also agree well with the observations of Liu et al. [9], who found that both toxic M. aeruginosa and thermal stress resulted in the enhancement of the enzymatic activities of SOD, CAT and GPx in H. cumingii.In Dreissena polymorpha, elevated GST and SOD activity was observed after exposure to MC, whereas GST, SOD and CAT activity were inhibited or barely affected in Unio tumidus [33].Different trends in the SOD, CAT, GPx and GR activities observed in various bivalve species could be attributed to interspecies differences as well as the significant influence of thermal stress on basal metabolic levels in bivalves [37].
Table 1.Summary of three-way ANOVA results for the effects of temperature (Temp), toxic M. aeruginosa (M) and sampling time points (Time) on biochemical parameters (glutathione S−transferase, GST; superoxide dismutase, SOD; catalase, CAT; glutathione peroxidase, GPx; glutathione reductase, GR; reactive oxygen species, ROS; malondialdehyde, MDA) in C. fluminea as a function of time.The thermal treatments were 20, 25 and 30 • C The concentrations of toxic M. aeruginosa were 0, 1 and 10 × 10 5 cells/mL.The sampling time points were days 1, 3, 5 and 7.

Assessment of ROS Levels and Lipid Oxidation
The MDA level serves as a significant biomarker of lipid peroxidation, providing information on the degree of cellular oxidative damage caused by excessive ROS [9].In this study, ROS production increased after toxic M. aeruginosa exposure, and a significant increase in ROS was observed in treated clams fed toxic M. aeruginosa compared to those fed C. vulgaris at the tested temperatures and at each sampling time point (p < 0.05) (Figure 3).Comparably, a significant enhancement of MDA content was observed in treated clams fed high concentrations of toxic M. aeruginosa compared to those fed C. vulgaris (Figure 3).And the MDA content only significantly increased in treated clams fed medium concentrations of toxic M. aeruginosa compared to those fed C. vulgaris at 30 °C (p < 0.05), while no significant increase was observed when the clams were exposed to low temperatures (20 °C) (Figure 3).These results indicated that thermal stress exacerbated oxidative damage in treated clams from the treatment groups receiving a low concentration of toxic M. aeruginosa.The three-way ANOVA results showed that toxic M. aeruginosa and temperature have significant effects on ROS and MDA (p < 0.05) (Table 1).Additionally, significant interactive effects between toxic M. aeruginosa and temperature were also observed in ROS and MDA (p < 0.05) (Table 1).Similarly, significantly elevated ROS and MDA content was observed in other bivalves, including H. cumingii, C. gigas and M. edulis, after toxic M. aeruginosa or MC exposure [9,32].

Assessment of ROS Levels and Lipid Oxidation
The MDA level serves as a significant biomarker of lipid peroxidation, providing information on the degree of cellular oxidative damage caused by excessive ROS [9].In this study, ROS production increased after toxic M. aeruginosa exposure, and a significant increase in ROS was observed in treated clams fed toxic M. aeruginosa compared to those fed C. vulgaris at the tested temperatures and at each sampling time point (p < 0.05) (Figure 3).Comparably, a significant enhancement of MDA content was observed in treated clams fed high concentrations of toxic M. aeruginosa compared to those fed C. vulgaris (Figure 3).And the MDA content only significantly increased in treated clams fed medium concentrations of toxic M. aeruginosa compared to those fed C. vulgaris at 30 • C (p < 0.05), while no significant increase was observed when the clams were exposed to low temperatures (20 • C) (Figure 3).These results indicated that thermal stress exacerbated oxidative damage in treated clams from the treatment groups receiving a low concentration of toxic M. aeruginosa.The threeway ANOVA results showed that toxic M. aeruginosa and temperature have significant effects on ROS and MDA (p < 0.05) (Table 1).Additionally, significant interactive effects between toxic M. aeruginosa and temperature were also observed in ROS and MDA (p < 0.05) (Table 1).Similarly, significantly elevated ROS and MDA content was observed in other bivalves, including H. cumingii, C. gigas and M. edulis, after toxic M. aeruginosa or MC exposure [9,32].

Integrated Analysis of Biochemical Parameters
The PCA approach is a factor analysis that achieves a significant reduction in the dimensionality of the original data set, revealing a meaningful association among the variables of interest [38].In this study, the PCA of the biochemical parameters showed that 82.13% of the total variance was accounted for by the first two principal components (Figure 4A).PC1 accounted for 65.23% of the overall variation, with GR exhibiting the only negative contribution.PC2 was responsible for 16.90% of the total variance and was driven by CAT, ROS and MDA.The PCA results demonstrated that temperature played an important role when toxic M. aeruginosa was removed, as evidenced by the pronounced clustering of the three temperature treatments without toxic M. aeruginosa exposure (Figure 4B).However, when accompanied by toxic M. aeruginosa exposure, the clustering effects resulting from temperature show a degree of attenuation, indicating the apparent interaction between toxic M. aeruginosa and thermal stress.
The IBR index is a commonly used quantitative measure in toxicological studies that effectively evaluates individual biomarker responses [39].It serves as a comprehensive indicator of the health state of aquatic organisms exposed to environmental stressors and toxic substances, integrating all relevant responses [40].In this study, the IBR values of the treatment groups consistently followed comparable patterns at various time intervals, showing a consistent increase as the concentrations of toxic M. aeruginosa increased at the same temperatures (Figure 5).Our observations were consistent with the results of Li et al. [41], who reported that the elevated IBR values integrated with biochemical parameters in the clam Ruditapes philippinarum and scallop Chlamys farreri were significantly correlated to the chemical contaminants' concentrations.Similarly, a previous study also found that the IBR values in R. philippinarum increased after exposure to heavy metal contaminations [42].In the present study, the results proved that toxic M. aeruginosa exposure has an apparent influence on the oxidative stress responses of C. fluminea at each temperature.

Integrated Analysis of Biochemical Parameters
The PCA approach is a factor analysis that achieves a significant reduction in the dimensionality of the original data set, revealing a meaningful association among the variables of interest [38].In this study, the PCA of the biochemical parameters showed that 82.13% of the total variance was accounted for by the first two principal components (Figure 4A).PC1 accounted for 65.23% of the overall variation, with GR exhibiting the only negative contribution.PC2 was responsible for 16.90% of the total variance and was driven by CAT, ROS and MDA.The PCA results demonstrated that temperature played an important role when toxic M. aeruginosa was removed, as evidenced by the pronounced clustering of the three temperature treatments without toxic M. aeruginosa exposure (Figure 4B).However, when accompanied by toxic M. aeruginosa exposure, the clustering effects resulting from temperature show a degree of attenuation, indicating the apparent interaction between toxic M. aeruginosa and thermal stress.
The IBR index is a commonly used quantitative measure in toxicological studies that effectively evaluates individual biomarker responses [39].It serves as a comprehensive indicator of the health state of aquatic organisms exposed to environmental stressors and toxic substances, integrating all relevant responses [40].In this study, the IBR values of the treatment groups consistently followed comparable patterns at various time intervals, showing a consistent increase as the concentrations of toxic M. aeruginosa increased at the same temperatures (Figure 5).Our observations were consistent with the results of Li et al. [41], who reported that the elevated IBR values integrated with biochemical parameters in the clam Ruditapes philippinarum and scallop Chlamys farreri were significantly correlated to the chemical contaminants' concentrations.Similarly, a previous study also found that the IBR values in R. philippinarum increased after exposure to heavy metal contaminations [42].In the present study, the results proved that toxic M. aeruginosa exposure has an apparent influence on the oxidative stress responses of C. fluminea at each temperature.

Molecular Processes Revealed via Transcriptome Analysis
To investigate the molecular processes behind the tolerance of C. fluminea to MCs under thermal stress, the RNA-Seq approach was employed to examine transcriptome alterations in C. fluminea co-exposed to a high concentration of M. aeruginosa (10 × 10 5 cells/mL) and high temperatures (30 • C).After 7 days of co-exposure to M. aeruginosa and thermal stress, a total of 1105 DEGs, including 761 significantly up-regulated genes and 344 significantly down-regulated genes, were observed in the treatment group compared to the control group (fed C. vulgaris at 20 • C) (Figure 6).Subsequently, the DEGs were subjected to eggNOG and GO annotations, as well as KEGG pathway analysis.

Molecular Processes Revealed via Transcriptome Analysis
To investigate the molecular processes behind the tolerance of C. fluminea to MCs under thermal stress, the RNA-Seq approach was employed to examine transcriptome alterations in C. fluminea co-exposed to a high concentration of M. aeruginosa (10 × 10 5 cells/mL) and high temperatures (30 °C).After 7 days of co-exposure to M. aeruginosa and thermal stress, a total of 1105 DEGs, including 761 significantly up-regulated genes and 344 significantly down-regulated genes, were observed in the treatment group compared to the control group (fed C. vulgaris at 20 °C) (Figure 6).Subsequently, the DEGs were subjected to eggNOG and GO annotations, as well as KEGG pathway analysis.eggNOG is a publicly accessible database for the fast functional annotation and orthology prediction of custom genomics or metagenomics data sets [43].This study involved the annotation of 717 DEGs into 21 eggNOG terms.Except for group of function unknown, function class of posttranslational modification, and protein turnover, chaperones have the most DEGs (17.85% of annotated DEGs), followed by intracellular trafficking, secretion, and vesicular transport (7.39%), and carbohydrate transport and metabolism (4.88%) (Figure 7).The top 30 GO terms with the most significant DEGs belong to three domains: molecular function (MF) included 15 subcategories, cellular component (CC) included 3 subcategories, and biological process (BP) included 12 subcategories (Figure 8).In MF, glutathione peroxidase activity, sodium: phosphate symporter activity and sodium-dependent phosphate transmembrane transporter activity constituted the top three clusters.In CC, the top three GO clusters were extracellular region, mitochondrial matrix and cytosol.The three most significant GO groups in BP were responses to oxidative stress, cellular phosphate ion homeostasis and the regulation of translational fidelity.The results revealed by DEG annotations confirmed that a substantial portion of DEGs exhibited associations with responses to oxidative stress and transporter activity which agree well with enzymatic observations in the activation of biotransformation and antioxidant defense systems, suggesting that C. fluminea could initiate integrated metabolic processes in response to oxidative stress caused by microcystins and thermal stress.eggNOG is a publicly accessible database for the fast functional annotation and orthology prediction of custom genomics or metagenomics data sets [43].This study involved the annotation of 717 DEGs into 21 eggNOG terms.Except for group of function unknown, function class of posttranslational modification, and protein turnover, chaperones have the most DEGs (17.85% of annotated DEGs), followed by intracellular trafficking, secretion, and vesicular transport (7.39%), and carbohydrate transport and metabolism (4.88%) (Figure 7).The top 30 GO terms with the most significant DEGs belong to three domains: molecular function (MF) included 15 subcategories, cellular component (CC) included 3 subcategories, and biological process (BP) included 12 subcategories (Figure 8).In MF, glutathione peroxidase activity, sodium: phosphate symporter activity and sodium-dependent phosphate transmembrane transporter activity constituted the top three clusters.In CC, the top three GO clusters were extracellular region, mitochondrial matrix and cytosol.The three most significant GO groups in BP were responses to oxidative stress, cellular phosphate ion homeostasis and the regulation of translational fidelity.The results revealed by DEG annotations confirmed that a substantial portion of DEGs exhibited associations with responses to oxidative stress and transporter activity which agree well with enzymatic observations in the activation of biotransformation and antioxidant defense systems, suggesting that C. fluminea could initiate integrated metabolic processes in response to oxidative stress caused by microcystins and thermal stress.The analysis of KEGG functional classifications indicated that the lysosome pathway was significantly enriched after treatment with toxic M. aeruginosa and thermal stress (FDR < 0.05) (Figure 9).The three most abundant KEGG pathways were lysosome, phagosome and protein processing in the endoplasmic reticulum, which were closely related to the inflammatory response (Figure 9).Furthermore, other abundant KEGG pathways were involved in immune functions, including the ECM-receptor interaction and FoxO signaling pathway, and metabolism functions, including arachidonic acid metabolism and glutathione metabolism (Figure 9).Arachidonic acid, known as the precursor to prostaglandins, has a wide range of physiological functions, including acting as a pro-inflammatory mediator and promoting the oxidation of exogenous compounds via co-oxidation [44].The activation of the pathways related to inflammatory and immune responses could be attributed to the oxidative stress induced by high temperatures and MCs.Additionally, the oxidation and co-oxidation occurring intracellularly play a significant role in the metabolic transformation of cyanotoxins into more polar molecules [45].Subsequently, these polar compounds were subjected to glutathione metabolism-mediated conjugation reactions and eventually discharged [46].Therefore, it is reasonable to assume that the metabolic processes of co-oxidation (phase I reaction) and conjugation (phase II reaction) in combination contribute to the biodetoxification of MCs in C. fluminea under thermal stress.The analysis of KEGG functional classifications indicated that the lysosome pathway was significantly enriched after treatment with toxic M. aeruginosa and thermal stress (FDR < 0.05) (Figure 9).The three most abundant KEGG pathways were lysosome, phagosome and protein processing in the endoplasmic reticulum, which were closely related to the inflammatory response (Figure 9).Furthermore, other abundant KEGG pathways were involved in immune functions, including the ECM-receptor interaction and FoxO signaling pathway, and metabolism functions, including arachidonic acid metabolism and glutathione metabolism (Figure 9).Arachidonic acid, known as the precursor to prostaglandins, has a wide range of physiological functions, including acting as a pro-inflammatory mediator and promoting the oxidation of exogenous compounds via co-oxidation [44].The activation of the pathways related to inflammatory and immune responses could be attributed to the oxidative stress induced by high temperatures and MCs.Additionally, the oxidation and co-oxidation occurring intracellularly play a significant role in the metabolic transformation of cyanotoxins into more polar molecules [45].Subsequently, these polar compounds were subjected to glutathione metabolism-mediated conjugation reactions and eventually discharged [46].Therefore, it is reasonable to assume that the metabolic processes of co-oxidation (phase Ⅰ reaction) and conjugation (phase Ⅱ reaction) in combination contribute to the biodetoxification of MCs in C. fluminea under thermal stress.

Conclusions
This study provides the first comprehensive analysis of the oxidative stress responses and molecular mechanisms in C. fluminea under co-exposure to toxic M. aeruginosa and thermal stress.The results indicated that the antioxidant enzyme activity, ROS levels and Figure 9. Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis of the differentially expressed genes (DEGs).Coloring represents the q value.A lower q value represents a more significant enrichment.The point size is proportional to the number of DEGs.

Figure 6 .
Figure 6.Volcano plot of differentially expressed genes (DEGs) in the hepatopancreas of C. fluminea fed with toxic M. aeruginosa at 30 °C compared with the control (fed with C. vulgaris at 20 °C) at 7 days post-treatment.Red dots represent up-regulated genes and green dots represent down-regulated genes.

Figure 6 .
Figure 6.Volcano plot of differentially expressed genes (DEGs) in the hepatopancreas of C. fluminea fed with toxic M. aeruginosa at 30 • C compared with the control (fed with C. vulgaris at 20 • C) at 7 days post-treatment.Red dots represent up-regulated genes and green dots represent downregulated genes.

Figure 7 .
Figure 7. Classification of the differentially expressed genes (DEGs) annotated in eggNOG.The Xaxis represents the names of the annotated groups, and the Y-axis corresponds to the percentage of the number of DEGs in the group accounting for the total number of annotated DEGs.

Figure 7 .
Figure 7. Classification of the differentially expressed genes (DEGs) annotated in eggNOG.The X-axis represents the names of the annotated groups, and the Y-axis corresponds to the percentage of the number of DEGs in the group accounting for the total number of annotated DEGs.

Figure 7 .
Figure 7. Classification of the differentially expressed genes (DEGs) annotated in eggNOG.The Xaxis represents the names of the annotated groups, and the Y-axis corresponds to the percentage of the number of DEGs in the group accounting for the total number of annotated DEGs.

Figure 8 . 16 Figure 8 .
Figure 8. Gene Ontology (GO) categories and the pattern of the differentially expressed genes (DEGs).The distribution of the GO categories was split into three categories: molecular functions (MF), cellular component (CC) and biological process (BP).A lower p value represents more significant enrichment.

Figure 9 .
Figure 9. Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis of the differentially expressed genes (DEGs).Coloring represents the q value.A lower q value represents a more significant enrichment.The point size is proportional to the number of DEGs.