Molecular Biomarkers of the Mitochondrial Quality Control Are Differently Affected by Hypoxia-Reoxygenation Stress in Marine Bivalves Crassostrea gigas and Mytilus edulis

Coastal environments commonly experience strong oxygen fluctuations. Resulting hypoxia/reoxygenation stress can negatively affect mitochondrial functions, since oxygen deficiency impairs ATP generation, whereas a surge of oxygen causes mitochondrial damage by oxidative stress. Marine intertidal bivalves are adapted to fluctuating oxygen conditions, yet the underlying molecular mechanisms that sustain mitochondrial integrity and function during oxygen fluctuations are not yet well understood. We used targeted mRNA expression analysis to determine the potential involvement of the mitochondrial quality control mechanisms in responses to short-term hypoxia (24 h at <0.01% O2) and subsequent reoxygenation (1.5 h at 21% O2) in two hypoxia-tolerant marine bivalves, the Pacific oysters Crassostrea gigas and the blue mussels Mytilus edulis. We hypothesized that the genes involved in the mitochondrial quality control will be upregulated during hypoxia, and the less hypoxia-tolerant of the two studied species (M. edulis) will show a stronger dependence on transcriptional upregulation of these pathways than C. gigas. To test these hypotheses, mRNA expression of 17 (C. gigas) and 11 (M. edulis) marker genes involved in mitochondrial fusion, fission, proteolysis and mitophagy was analyzed in the digestive gland of M. edulis and C. gigas in normoxia and during hypoxia-reoxygenation (H/R) stress. In the mussels, the mRNA expression of the transcripts related to mitochondrial dynamics and quality control was strongly altered during H/R stress showing a shift toward fission, suppression of fusion, an increase in mitochondrial proteolysis and onset of mitophagy. These changes indicate that H/R stress induces mitochondrial injury in M. edulis requiring upregulation of the protective mechanisms to segregate the dysfunctional mitochondria by fission and degrade the oxidative damaged proteins and/or organelles. Unlike mussels, the transcript levels of all studied genes in the oysters remained at the baseline (normoxic) levels during H/R stress. This muted transcriptional response of C. gigas is in agreement with earlier findings showing better ability to maintain cellular homeostasis and higher resistance to apoptosis during H/R stress in the oysters compared with the mussels. The revealed species-specific differences in the expression of the mitochondrial quality control pathways shed light on the potentially important mechanisms of mitochondrial protection against H/R-induced damage that might contribute to hypoxia tolerance in marine bivalves.

Coastal environments commonly experience strong oxygen fluctuations. Resulting hypoxia/reoxygenation stress can negatively affect mitochondrial functions, since oxygen deficiency impairs ATP generation, whereas a surge of oxygen causes mitochondrial damage by oxidative stress. Marine intertidal bivalves are adapted to fluctuating oxygen conditions, yet the underlying molecular mechanisms that sustain mitochondrial integrity and function during oxygen fluctuations are not yet well understood. We used targeted mRNA expression analysis to determine the potential involvement of the mitochondrial quality control mechanisms in responses to shortterm hypoxia (24 h at <0.01% O 2 ) and subsequent reoxygenation (1.5 h at 21% O 2 ) in two hypoxia-tolerant marine bivalves, the Pacific oysters Crassostrea gigas and the blue mussels Mytilus edulis. We hypothesized that the genes involved in the mitochondrial quality control will be upregulated during hypoxia, and the less hypoxiatolerant of the two studied species (M. edulis) will show a stronger dependence on transcriptional upregulation of these pathways than C. gigas. To test these hypotheses, mRNA expression of 17 (C. gigas) and 11 (M. edulis) marker genes involved in mitochondrial fusion, fission, proteolysis and mitophagy was analyzed in the digestive gland of M. edulis and C. gigas in normoxia and during hypoxia-reoxygenation (H/R) stress. In the mussels, the mRNA expression of the transcripts related to mitochondrial dynamics and quality control was strongly altered during H/R stress showing a shift toward fission, suppression of fusion, an increase in mitochondrial proteolysis and onset of mitophagy. These changes indicate that H/R stress induces mitochondrial injury in M. edulis requiring upregulation of the protective mechanisms to segregate the dysfunctional mitochondria by fission and degrade the oxidative damaged proteins INTRODUCTION Deoxygenation driven by nutrient pollution and warming is a major stressor in estuarine and coastal oceans worldwide (Breitburg et al., 2018(Breitburg et al., , 2019. The deficiency (hypoxia) or lack (anoxia) of oxygen occurs in the estuarine and coastal habitats where oxygen depletion due to respiration outstrips the oxygen influx from diffusion, water mixing and photosynthesis. As a result, these habitats can experience hypoxic periods from a few hours (such as in a diel cycling hypoxia) to several weeks or months in the case of seasonal hypoxia Rosenberg, 1995, 2008). Most benthic metazoans including important ecosystem engineers and keystone species of bivalves, crustaceans and fish require oxygen to survive and complete their life cycles, and oxygen deficiency can lead to a decrease in growth, performance and survival of marine benthos thereby threatening coastal marine ecosystems (Diaz and Rosenberg, 1995). Low motility of benthic marine invertebrates such as bivalves limits their ability to escape local hypoxia, making them reliant on physiological and biochemical adaptations to survive oxygen deficiency (Grieshaber et al., 1994). Metabolic regulation is at the crux of such adaptations due to the direct dependence of the mitochondria (generating >90% of the cellular ATP) on oxygen availability and the key role of these organelles in the energetic and redox balance of the cell. Investigation of mitochondrial responses to fluctuating oxygen conditions is thus important for understanding the fundamental physiological mechanisms setting limits of hypoxic tolerance in marine organisms.
In benthic invertebrates, severe hypoxia (<4% O 2 ) and anoxia results in transition to a metabolically suppressed state, where ATP production and consumption are downregulated and anaerobic pathways are engaged to compensate for insufficient aerobic ATP generation (Hochachka et al., 1996;Hochachka and Lutz, 2001;Storey, 2002). Under these conditions, the mitochondrial respiration strongly declines or ceases (Zwaan et al., 1991;van den Thillart et al., 1992). In hypoxia-sensitive organisms such as terrestrial mammals, mitochondria become depolarized and eventually damaged during hypoxia due to the ATP depletion, acidosis and Ca 2+ overload (Piper et al., 2003;Chen et al., 2007;Solaini et al., 2010). Post-hypoxic reoxygenation restores the aerobic ATP synthesis in mammalian mitochondria but may come at a cost due to a surge of the mitochondrial production of reactive oxygen species (ROS) that can damage the cellular proteins, lipids and DNA (Cadenas and Davies, 2000;Paradis et al., 2016). Similar patterns of mitochondrial responses during hypoxia-reoxygenation (H/R) stress involving oxidative damage, collapse of the mitochondrial membrane potential and loss of the OXPHOS capacity has been reported in hypoxia-sensitive marine mollusks such as the bay scallops . In hypoxia-tolerant organisms such as freshwater turtles, fish and marine intertidal mollusks, mitochondrial respiration, OXPHOS capacity and membrane potential are preserved and the oxidative damage is mitigated during the H/R stress (Galli and Richards, 2014;Sokolova, 2018;Sokolova et al., 2019). These findings imply that maintenance of the mitochondrial integrity through mitochondrial damage sensing and quality control mechanisms might play an important role in mitochondrial tolerance to H/R stress (Berlett and Stadtman, 1997;Cabiscol et al., 2000;Zorov et al., 2014). However, the mitochondrial quality control mechanisms remain poorly understood in hypoxiatolerant organisms including benthic marine bivalves and require further investigations.
The mitochondrial quality control is well studied in hypoxia-sensitive terrestrial mammals (such as rodents) and encompasses a complex protective network of proteins regulating mitochondrial dynamics (fission and fusion), as well as proteases and mitophagic pathways degrading oxidatively damaged proteins and defective mitochondria (Eisner et al., 2018;Li and Liu, 2018). Under normal conditions, mitochondria tend to form mitochondrial networks by mitochondrial fusion. Fusion of the outer mitochondrial membrane (OMM) is mediated by mitofusins 1 and 2 (Mfn1 and 2), while dynamin-related GTPase OPA1 induces fusion of the inner mitochondrial membrane (IMM) (Chen et al., 2003;Song et al., 2009). Under stress, membrane depolarization activates the metalloendopeptidase OMA1, which degrades OPA1 Xiao et al., 2014) and induces mitochondrial fission through interactions of dynamin-related protein DRP1 and its adaptor proteins (the mitochondrial fission factor Mff and fission protein 1 Fis1) (Otera et al., 2010;Losón et al., 2013). Mitochondrial fission removes the damaged mitochondria and facilitates apoptosis under extreme cellular stress (Youle and van der Bliek, 2012), while the hypoxia upregulated protein (HYOU1) supports the mitochondrial repair by suppressing premature hypoxia-induced cell death (Ozawa et al., 1999).
Mechanisms involved in the degradation of the damaged proteins and organelles play an important role in the cellular responses to hypoxia in mammalian models. Protein quality surveillance in the mitochondria is supported by multiple ATP-dependent and -independent proteases including Lon protease and paraplegin (both selectively degrading oxidatively damaged proteins) (Atorino et al., 2003;Venkatesh et al., 2012;Kuo et al., 2015;Pinti et al., 2015;Shanmughapriya et al., 2015;Sepuri et al., 2017), the IMM protease ATP23, and the caseinolytic matrix peptidase chaperone (Osman et al., 2007;Doyle and Wickner, 2008). On DNA level, the helicase Twinkle is responsible for mtDNA maintenance and replication (Spelbrink et al., 2001). Cells can also degrade complete damaged mitochondria through mitophagy to prevent joining of dysfunctional mitochondria to the network. These pathways involve the PTEN-induced kinase 1 (PINK1) which accumulates under the hypoxia and recruits Parkin to tag damaged mitochondria for degradation (Eiyama and Okamoto, 2015). Mitophagic mediators in the PINK1-Parkin pathway are regulated via reversible phosphorylation by a serine/threonineprotein phosphatase PGAM5 (Chen et al., 2016). Additional to the PINK1-Parkin pathway, an alternative mitophagic pathway degrading unhealthy mitochondria is mediated by the mitochondrial eating protein (MIEAP) (Kitamura et al., 2011;Nakamura and Arakawa, 2017).
The important role of mitochondrial quality control mechanisms (involving the regulation of the mitochondrial fission and fusion, proteolysis and mitophagy) in the responses to H/R stress is well established in the hypoxia-sensitive mammalian models such as rodents (Kulek et al., 2020;Wang and Zhou, 2020). Earlier studies also indicate that high activity and upregulated mRNA expression of mitochondrial proteases correlates with elevated hypoxia tolerance in marine bivalves . However, involvement of the pathways regulating mitochondrial dynamics and quality control in the response to H/R stress have not been extensively studied in hypoxia-tolerant organisms including bivalves. Here, we explored the molecular mechanisms underlying the mitochondrial responses to H/R stress by measuring the mRNA expression of key marker genes in the mitochondrial quality control pathways in the Pacific oyster Crassostrea gigas and the blue mussel Mytilus edulis. The two studied species are commonly exposed to fluctuating oxygen conditions in intertidal, coastal and estuarine habitats and are therefore good models to study the mitochondrial adaptations to hypoxia and reoxygenation stress .
Earlier studies found that the Pacific oyster is more tolerant to abiotic stressors (including hypoxia) than the blue mussel (David et al., 2005;Le Moullac et al., 2007;. At the cellular level, higher hypoxia tolerance of C. gigas was associated with the resistance to apoptosis and muted inflammation response indicating less pronounced mitochondrial and cellular damage during the H/R stress in oysters compared with the mussels (Falfushynska et al., 2020a,b). Therefore, we hypothesized that the less hypoxia-tolerant of the two studied species (M. edulis) will more strongly rely on the regulation of the mitochondrial quality control mechanisms to counteract the H/R-induced cellular damage compared with the more tolerant C. gigas that experiences less cellular injury (Falfushynska et al., 2020b) and thus might require less protection from the inducible mitochondrial quality control mechanisms. To test this hypothesis, mRNA expression of 17 (C. gigas) and 11 (M. edulis) key genes in the mitochondrial fusion, fission, proteolysis and mitophagy pathways were analyzed by quantitative RT-PCR in the digestive gland of the two studied species of bivalves exposed to hypoxia and reoxygenation relative to the normoxic controls. We focused on the transcript levels of the following marker genes: for mitochondrial fission and fusion -mfn2 (encoding mitofusin 2), opa1 (mitochondrial dynamin-like 120 kDa protein), dnm1l (dynamin-1-like protein), mff (mitochondrial fission factor), fis1 (mitochondrial fission protein 1); for protein and DNA quality control -tsfm (encoding mitochondrial translation elongation factor Ts), lonp1 (mitochondrial Lon protease), spg7 (paraplegin), oma1 (mitochondrial metalloendopeptidase OMA1), clpB (mitochondrial caseinolytic matrix peptidase chaperone subunit B), atp23 (mitochondrial inner membrane protease ATP23), twnk (mitochondrial twinkle mtDNA helicase); and for mitophagy -mieap (encoding mitochondrial eating protein), hyou1 (hypoxia upregulated protein 1), prkn (parkin), pink1 (PTEN-induced kinase 1), and pgam5 (mitochondrial serine/threonine protein phosphatase PGAM5).

Animal Maintenance
Adult oysters [mean shell length ± the standard deviation (SD): 99.3 ± 13.1 mm] were obtained from the low intertidal zone of the German Wadden Sea near List/Sylt (55 • 01 42 N 8 • 26 04 E) and transported within 48 h of collection to the University of Rostock, Germany. Adult blue mussels (mean shell length ± SD: 52.5 ± 3.22 mm) were collected at Warnemünde marina "Mittelmole" (54 • 10 49.4 N 12 • 05 19.9 E) and brought to the University of Rostock within an hour of collection. All bivalves were transported in coolers lined with seawater-soaked paper towels, and the shells were cleaned from epibionts upon arrival. Two weeks prior to experiments, oysters and mussels were acclimated in recirculated temperature-controlled aquarium systems (Kunststoff-Spranger GmbH, Plauen, Germany) with aerated artificial seawater (ASW) (Tropic Marin R , Wartenberg, Germany) at salinity 33 ± 1 (C. gigas) and 15 ± 1 (M. edulis), respectively, and temperature 15 ± 0.5 • C. These salinity and temperature were within natural range of the respective oysters' and mussels' habitat conditions. Bivalves were fed ad libitum by continuous addition of a commercial live algal blend (DTs Premium Blend Live Marine Phytoplankton, Coralsands, Mainz Kastel, Germany) according to the manufacturer's instructions (0.01 ml g −1 fresh mass per day) using an automatic aquarium feeder.

Experimental Exposures
Bivalves (three oysters or six mussels in 2 l of ASW) were exposed to 24 h of severe hypoxia (<0.1% O 2 ) by bubbling the ASW with pure nitrogen (Westfalen AG, Münster, Germany) in air-tight chambers at 15 ± 0.5 • C and respective salinity. Oxygen concentration was monitored with an Intellical TM LDO101 Laboratory Luminescent/Optical Dissolved Oxygen (DO) Sensor (HACH, Loveland, CO, United States). During exposure, bivalves were not fed to prevent bacterial growth in the chambers. It is worth noting that oysters and mussels spontaneously close their shells and cease feeding during hypoxic exposures regardless of the presence of algae in surrounding water. After hypoxia exposure, a subset of bivalves was allowed to recover in normoxic ASW (21% O 2 ) for 1.5 h. The H/R incubation times were chosen based on tidal and diurnal cycles of oxygen fluctuations and previous findings that show that strongest mitochondrial response to reoxygenation occurs within the first hours of recovery (Kurochkin et al., 2009;Richards, 2011;Andrienko et al., 2017). The control group was maintained in normoxia (21% O 2 ) in recirculated temperaturecontrolled aquarium systems (10 oysters or 30 mussels per 30 l of seawater) with continuous feeding as described in section "Animal Maintenance." Hepatopancreas tissue was sampled on ice after hypoxia, reoxygenation and in normoxic control, immediately shockfrozen in liquid nitrogen and stored at −80 • C until further analysis. We have chosen hepatopancreas because it is one of the largest metabolically active organs in bivalves involved in digestion, energy storage and maintenance of energy homeostasis (Gosling, 1992;Kennedy, 1996). Furthermore, earlier studies in bivalves showed strong hypoxia-induced shifts in the hepatopancreas metabolome (Haider et al., 2020) and induction of apoptotic and inflammatory responses in this tissue during H/R stress (Falfushynska et al., 2020a). This makes hepatopancreas a useful tissue to investigate the pathways that maintain the mitochondrial integrity during H/R exposure. Samples sizes were 10 and 6 per treatment group for oysters and mussels, respectively.

Quantitative Real-Time PCR (qRT-PCR)
RNA extraction and cDNA synthesis was conducted as described elsewhere (Falfushynska et al., 2020a). Briefly, 20 to 50 mg hepatopancreas tissue were extracted in RNA Extracol (EURX Ltd.-molecular biology products, Gdansk, Poland) according to manufacturer's instructions using an automatic homogenizer for 40 s at 6.0 m s −1 (FastPrep-24 MP Biomedicals, Valiant Co. Ltd., Yantai, Shandong, China). To avoid potential crossreaction with residual DNA, RNA extracts were treated with DNase using the TURBO DNA-free kit (ThemoFisher Scientific, Waltham, MA, United States). RNA with a 280/260 absorbance ratio > 2.0 assessed by NanoVue Plus TM spectrophotometer (GE Healthcare, Chicago, IL, United States) was used for further processing. cDNA was synthesized from 2 µg of total RNA using RevertAid RT Kit (Thermo Fischer Scientific, Waltham, MA, United States). Expression levels were quantified by quantitative PCR (qPCR) using the StepOnePlus Realtime-PCR System instrument (Applied Biosystems, Thermo Fisher Scientific Corp., Foster City, CA, United States) and Biozym Blue S'Green qPCR BlueMix Separate ROX kit (Biozym Scientific GmbH, Hessisch Oldendorf, Germany). Each reaction comprised of 10 µl 2x qPCR master mix with ROX additive, 6.4 µl DNase/RNase-free water, 1.6 µl gene-specific reverse and forward primer mix (each to a final concentration of 0.4 µM) and 2 µl cDNA. For the genes where species-specific annotated transcript sequences were available, gene-specific primers were designed using published gene sequences of M. edulis and C. gigas in the NCBI/GenBank database. This was the case for all C. gigas target genes (except lonp1, opa1, hyou1 and fis1) and for the housekeeping genes from M. edulis. For those Mytilus sequences where no annotated transcripts were found in the NCBI/GenBank, we have used a transcriptome from a closely related species Mytilus galloprovincialis that belongs to the Mytilus edulis species complex (Varvio et al., 1988;Gaitán-Espitia et al., 2016). M. edulis and M. galloprovincialis are closely related, genetically similar and characterized by a high degree of hybridization and introgression in the nature (Gosling and Wilkins, 1981;Varvio et al., 1988;Gardner, 1996;Daguin et al., 2001;Gaitán-Espitia et al., 2016). To identify the candidate regions for primer design in Mytilus, we have aligned M. galloprovincialis sequences with available homologous sequences from multiple molluscan species, including Crassostrea. The regions with the high degree of nucleotide sequence conservation were chosen for primer design, and the qPCR was conducted under stringent conditions to avoid non-specific priming. The relevant fragments of sequence alignments between Crassostrea and M. galloprovincialis are provided in Supplementary Table S1. For lonp1, opa1, hyou1 and fis1 genes, for which no transcript from C. gigas were publicly available at the time of this study, primers were designed using a sequence from a closely related species, Crassostrea virginica (Supplementary Table S2). The identified sequence segments were then used to design gene-specific primers ( Table 1). For primer check, the fragments were amplified using the following cycling parameters: 2 min at 95 • C for polymerase activation, 40 cycles of 15 s at 95 • C and 30 s at 60 • C. Only the primers that produced a single product of expected length were used in the following analyses. All 17 target genes were successfully amplified from C. gigas. Six out of 17 target genes (including tsfm, clpB, mieap, atp23, hyou1, and twnk) could not be amplified from Mytilus due to the low quality or insufficient length of the respective transcript sequences available in the NCBI/GenBank.
To measure mRNA levels of the target and housekeeping genes, the following cycling parameters were used: 2 min at 95 • C for polymerase activation, 40 cycles of 15 s at 95 • C, 15 s at 57 • C (annealing), 20 s at 72 • C (extension) and 30 s at a primer specific reading temperature (for acquiring the product signal). For quality control, the melt curve analysis was conducted at the end of each run. To determine the apparent amplification efficiency of the primers and to correct for the possible plateto-plate variation, a standard dilution series of a single pooled cDNA sample was run on each plate. All measurements were conducted in duplicates. Expression levels of the target samples relative to the reference genes were conducted using the relative standard curve method (Dorak, 2006). Reference genes were chosen as described elsewhere (Falfushynska et al., 2020a). In pilot experiments, we tested five potential housekeeping genes (encoding eukaryotic elongation factor 1, α-tubulin, β-actin, 1 | Primers used for qRT-PCR for genes of mitochondrial quality control in C. gigas and M. edulis.

Gene
Forward primer (5 -3 ) Reverse primer (5 -3 ) NCBI accession No. of targeted sequence ubiquitin and 18S rRNA) as potential reference genes. Only one out of five tested genes (eukaryotic elongation factor 1 in mussels and actin in oysters) showed no significant change in expression (p > 0.05) with ≤1 difference in mean CT values between different experimental treatments. Therefore, target mRNA expression was normalized to the expression of a reference gene (β-actin βAct for C. gigas and the eukaryotic elongation factor 1 eEF1 tsfm for M. edulis) and the mean CT of the normoxic group as described elsewhere (Dorak, 2006)  Significant differences between exposures were tested by a one-way ANOVA within each species using SigmaPlot 13. Oxygen regime was used as a fixed factor with three levels (normoxic control, 1 day hypoxia, reoxygenation). Significant differences between the pairs of means were determined using a Tukey's Honest Significant Differences (HSD) post hoc test. Pearson correlation analysis was carried out by RStudio's packages "mixOmics" and "factoextra" [v. 1.2.5033 and R v. 3.6.3 (R Core Team, 2020)] with cut off R = 0.6. To reduce the dimensionality of the data, the principal component analysis (PCA) was conducted on transformed data with outliers replaced by group means. PCA and graphs were designed in RStudio. All effects were considered significant at p < 0.05.

Mitochondrial Dynamics: Fusion and Fission
In M. edulis, severe hypoxia suppressed mRNA levels of mfn2 encoding mitofusin-2, which were rapidly restored during reoxygenation ( Figure 1A). Transcript levels of opa1 encoding dynamin-like 120 kDa protein did not change during H/R exposure in the mussels ( Figure 1C). In C. gigas, tissue levels of mfn2 and opa1 mRNA did not change in response to H/R stress (Figures 1B,D). Transcript levels of mitochondrial fission mediators dnm1l, mff and fis1 were significantly higher in the hepatopancreas of the mussels exposed to hypoxia than the control group (Figures 2A,C,F). Reoxygenation led to a further significant increase in the mRNA expression levels of these genes. Levels of mff mRNA were ∼4-fold and 6-fold elevated in hypoxia and reoxygenation, respectively, compared to control. Expression of fis1 mRNA was ∼2.5-to 3-fold increased in hypoxia and reoxygenation compared to normoxia. The dnm1l gene was ∼1.5-to 2-fold elevated by hypoxia and reoxygenation in the mussels compared to normoxic control. In C. gigas, transcript levels of the mitochondrial fission factors dnm1l, mff, tsfm and fis1 were not significantly affected by hypoxia and reoxygenation (Figures 2B,D,E,G).

Mitochondrial Protein and DNA Quality Control
Transcript levels of the genes encoding mitochondrial proteases involved in the protein quality control including lonp1 and spg7 were not affected by hypoxia and reoxygenation in M. edulis (Figures 3A,C). In contrast, transcript levels of metalloendopeptidase OMA1 (oma1) increased during hypoxia and reoxygenation by ∼1.5fold and ∼1.7-fold, respectively, compared to the normoxic controls ( Figure 3E).
In C. gigas, the mRNA level of the genes involved in mitochondrial protein quality control including lonp1, spg7, oma1 (Figures 3B,D,F) as well as four additional genes (not measured in M. edulis due to the lack of specific primers) including clpB, mieap, atp23 and hyou1 did not change in response to H/R stress (Figures 4A-D). Transcript levels of the mitochondrial helicase TWINKLE (twnk) responsible for mtDNA maintenance, remained unchanged in C. gigas throughout experimental exposures ( Figure 4E) (no data are available for M. edulis due to the lack of specific primers).

Mitophagy
In M. edulis, hypoxia induced a ∼2-fold higher expression of prkn mRNA (encoding a mitophagy mediator Parkin) compared to normoxic condition ( Figure 5A). The levels of prkn remained elevated during reoxygenation in the mussels (Figure 5A). Expression of its co-mediator pink1 did not change with hypoxia, while reoxygenation led to a ∼2.5-fold increase of pink1 expression compared to normoxic controls ( Figure 5C). Transcript levels of the protein phosphatase PGAM5 showed a decreasing trend with hypoxia, followed by a significant increase during reoxygenation ( Figure 5E). In C. gigas, the expression pattern of the studied mitophagy genes did not show pronounced effects of hypoxia and subsequent reoxygenation (Figures 5B,D,F).

Data Integration
The principal component analysis in M. edulis hepatopancreas revealed two principal components explaining 47.2% (PC1) and 19.1% (PC2) of the data variation ( Figure 6C and Supplementary Table S3). The 1 st component showed high positive loadings for genes of mitochondrial fission (dnm1l, mff, fis1) and fusion (opa1) and mitophagy (prkn, pink1) (Figure 6D and Supplementary Table S4). The 2 nd principal component had high positive loadings of genes of mitochondrial quality control (spg7, lonp1) and a strong negative loading of mfn2, mediator of mitochondrial fusion. In M. edulis, control and reoxygenation group were clearly separated along the 1 st principal component (Figure 6C). Hypoxia and to greater degree reoxygenation induced a shift toward more positive values of the 1 st principal component compared to the normoxic control reflecting upregulation of fission genes mff, fis1 and dnm1l and a gene encoding a mitochondrial protease oma1. All experimental treatment groups showed a high variance along the 2 nd principal component axis.
In C. gigas, the 1 st and 2 nd principal component explained 64.3% and 7.9% of variance, respectively ( Figure 6A and Supplementary Table S5). All genes analyzed in this study had high positive loadings on the 1 st principal component (>0.05) (Figure 6B and Supplementary Table S6).
Gene transcription of C. gigas from all three experimental treatments showed considerable overlap in the plane of the two first principal components consistent with the lack of change in the mRNA expression profiles of the studied genes in response to the H/R stress. In C. gigas, no apparent correlated gene clusters were found with the tight correlation across all studied genes (Figures 7A,B). These data indicate that the observed variability and correlations in the mRNA expression profiles of C. gigas mostly reflects individual differences between bivalves rather than the effects of the oxygen regime.

DISCUSSION
Short-term severe hypoxia and subsequent reoxygenation had different effects on mRNA expression of mitochondrial quality control in C. gigas and M. edulis ( Table 2). In C. gigas, transcript levels were not affected by H/R stress indicating that the mitochondrial injury thresholds needed to induce the protective quality control mechanisms might not have been reached in this hypoxia-tolerant species. In M. edulis, mRNA expression of the

Gene
Crassostrea gigas Mytilus edulis transcripts related to mitochondrial dynamics and quality control was strongly regulated during H/R stress showing a shift toward the predominance of mitochondrial fission over fusion and an onset of mitophagy.

Effects of H/R Stress on Transcriptional Regulation of Mitochondrial Quality Control in the Mussels
Mitochondrial dynamics (in particular fission and fusion) plays a critical role in maintaining functional mitochondria during stress exposures in animals (Youle and van der Bliek, 2012). In mammalian models, ischemia-reperfusion (I/R) stress regulates transcription and post-translational modifications of fusion and fission proteins reflecting the important role of the mitochondrial dynamics in adjustments to interrupted oxygen supply (Li and Liu, 2018). Cellular stress commonly shifts mitochondrial dynamics toward mitochondrial fission to allow degradation of damaged mitochondria and avoid re-joining of non-functional mitochondria into the mitochondrial network (Chen et al., 2016;Eisner et al., 2018;Sokolova et al., 2019). The GTPase Drp1/Dnm1l mediates mitochondrial fission in response to cellular stress signals such as increase of Ca 2+ concentration, ROS production or depolarized mitochondrial membrane as shown during H/R (I/R) stress in mammalian models (Yu et al., 2016;Tian et al., 2017;Yang et al., 2017;Li and Liu, 2018). Drp1/Dnm1l binds to the OMM with the help of the adaptor proteins Mff and Fis1 and initializes division of the membrane (Mears et al., 2011;van der Bliek et al., 2013;Osellame et al., 2016). Drp1/Dnm1l is a transcriptional target of p53 and, in rats, it is regulated at the transcriptional level (Li et al., 2010). The finding of transcriptional upregulation of dnm1l and its two adapter proteins mff and fis1 during H/R stress in M. edulis is consistent with upregulation of the fission pathways possibly as a mechanism to protect mitochondrial integrity by segregation of dysfunctional mitochondria (Twig et al., 2008;Anzell et al., 2018). mRNA expression of the fission-related proteins remained elevated during reoxygenation indicating that 1.5 h reoxygenation might not be sufficient for the full recovery of the mitochondrial fusion-fission balance in M. edulis.
Fusion is an effective strategy to protect against the stressinduced DNA damage, as the damaged mitochondria fuse with the healthy ones and replace the damaged mitochondrial DNA through complementation (Youle and van der Bliek, 2012). In our present study, mRNA levels of a key fusion mediator mfn2 (encoding the mitochondrial Mfn 2 protein) were suppressed during hypoxia in M. edulis indicating downregulation of the mitochondrial fusion pathways. This lack of induction of mfn2 might indicate that mitochondrial DNA damage is not critically involved in the H/R stress responses of the blue mussels. An earlier study in the eastern oysters C. virginica showed an increase in the levels of DNA damage during exposure to 6 days of hypoxia and subsequent reoxygenation in the oysters exposed to a toxic metal cadmium, but not in their counterparts exposed to H/R stress only (Kurochkin et al., 2009). In the mussels M. edulis, hypoxia caused by emersion (72 h) had no effect on the DNA damage assessed by the Comet assay (Singh and Hartl, 2012). In the brown mussels Perna perna, 24 h of emersioninduced hypoxia led to the elevated oxidative DNA damage in the gills but not in the digestive gland (Almeida et al., 2005). In these studies, the oxidative lesions were assessed in the total rather than the mitochondrial DNA. To the best of our knowledge, the effects of H/R stress on the mitochondrial DNA damage have not been studied in marine bivalves. It is worth noting that mitofusin 2 requires energy in the form of GTP to mediate the mitochondrial fusion (Filadi et al., 2018), so that the downregulation of mfn2 in hypoxia-exposed mussels might be a consequence of the metabolic rate depression, a common energysaving strategy to survive the oxygen deficiency in hypoxiatolerant mollusks (Hochachka, 1993;Storey, 1998;Storey and Storey, 2004). During reoxygenation mRNA expression of mfn2 returned to the normoxic (baseline) levels in the mussels, which indicates a rapid re-balancing of mitochondrial fusion during post-hypoxic recovery.
The transcript levels of opa1 (encoding the mitochondrial dynamin-like GTPase OPA1) were not affected by H/R stress in M. edulis. The mitochondrial OPA1 is involved in the regulation of the mitochondrial fusion and the cross-talk between fusion and fission (Cipolat et al., 2004). In mammals, OPA1 activity is regulated transcriptionally (Ryu et al., 2015;Zheng et al., 2019) as well as post-translationally through proteolytic cleavage by the mitochondrial metalloprotease OMA1 (Chen et al., 2016). Mitochondrial stress such as a decrease in the mitochondrial membrane potential or ATP deficiency stimulates the proteolytic processing of OPA1 generating short isoforms unable to support mitochondrial fusion Chen et al., 2016;Ali and McStay, 2018). The posttranslational regulation of OPA1 allows for rapid suppression of mitochondrial fusion and mitochondrial fragmentation but might not be reflected in the changes in mRNA expression (Ali and McStay, 2018). Nevertheless, opa1 was observed to be transcriptionally regulated during formation of immune cells in mice (Ryu et al., 2015), while disorders related to impaired mitochondrial functions were associated with transcriptional regulation of opa1 (Zheng et al., 2019). The observed stability of opa1 mRNA in the mussels might indicate the lack of involvement of OPA1 in the mitochondrial response to the H/R stress in M. edulis. Alternatively, if posttranslational regulation plays a predominant role in the OPA1 activation in the mussels, activation of OPA1 might not be reflected at the transcript level. Further studies are needed to distinguish between these alternative explanations.
In M. edulis, H/R stress led to a strong upregulation of oma1 mRNA encoding the mitochondrial metalloprotease OMA1. As discussed above, OMA1 mediates the proteolytic cleavage of fusion mediator OPA1 during stress thereby suppressing mitochondrial fusion . OMA1dependent degradation of OPA1 is thus essential for the selective removal of defective mitochondrial fragments and prevention of incorporation of those into the mitochondrial network (Head et al., 2009;Baker et al., 2014). Consequently, elevated oma1 expression during hypoxia and subsequent reoxygenation in the mussels might reflect the overall shift in the mitochondrial fusion-fission dynamics toward the predominance of fission (Ali and McStay, 2018) and is consistent with the transcriptional upregulation of the fission regulating proteins dnm1l, mff and fis1 and suppressed expression of a fusion regulator mfn2 in M. edulis during the H/R stress. Furthermore, OMA1 mediates an ATP-independent proteolytic breakdown of misfolded inner membrane proteins and is complementary to the action of m-AAA proteases (Käser et al., 2003). OMA1 activity is strongly upregulated in response to various stresses in metazoans, plants and yeast (Bohovych et al., 2014;Rainbolt et al., 2015;Migdal et al., 2017) and is required for stabilization of the mitochondrial respiratory supercomplexes and maintenance of the mitochondrial respiration (Bohovych et al., 2015;Richter et al., 2015). Therefore, transcriptional upregulation of oma1 during H/R stress in the mussels might assist in the degradation of the damaged mitochondrial proteins and protection of the mitochondrial integrity and function.
Our study shows that unlike oma1, two other studied mitochondrial proteases (including the mitochondrial Lon protease and paraplegin, encoded by lonp1 and spg7, respectively) are not transcriptionally upregulated during H/R stress in M. edulis. Lon protease is highly conserved throughout organisms and essential for mitochondrial proteostasis (Pinti et al., 2015). It is activated by ROS, protein carbonylation and lipid peroxidation (Kuo et al., 2015) and specifically targets the oxidatively damaged proteins for ATP dependent degradation (Bota and Davies, 2002;Bayot et al., 2010). Lon protease activity is regulated transcriptionally and post-translationally suggesting a fine tuning by several mechanisms, although stronger regulatory patterns are commonly seen at the protein level (Pinti et al., 2015). Earlier studies in bivalves showed that Lon protease mRNA expression was upregulated during H/R stress in the hypoxia-sensitive scallops Argopecten irradians but not in the hypoxia-tolerant hard shell clams Mercenaria mercenaria . Upregulation of lonp1 transcripts in the scallops was associated with accumulation of oxidative lesions (i.e., lipid and protein peroxidation products) in the mitochondria, whereas in the hard shell clams the mitochondrial oxidative stress markers remained at the baseline (normoxic) levels . Upregulation of lonp1 mRNA expression was also found in the Eastern oyster C. virginica during simultaneous exposure to high temperature and a toxic metal cadmium (Sanni et al., 2008), a combination of stressors known to cause strong oxidative stress in oysters (Bagwe et al., 2015). Taken together, these findings indicate that transcriptional upregulation of lonp1 in bivalves might occur when the oxidative stress reaches a species-specific threshold in the mitochondria, and the lack of transcriptional upregulation of lonp1 mRNA levels during H/R stress in M. edulis might indicate that this threshold has not been reached under the conditions of our present study. Similar to lonp1, transcript levels of spg7 gene encoding a mitochondrial metalloprotease paraplegin did not change during H/R stress in the mussels. To the best of our knowledge, expression of spg7 has not been studied in other marine bivalves. In model organisms including mammals and yeast, paraplegin is regulated by mitochondrial stress signals such as Ca 2+ accumulation and elevated ROS levels (Ishihara et al., 2006;Shanmughapriya et al., 2015). If similar ROS-dependent mechanisms exist in bivalves, the lack of spg7 (as well as lonp1) transcriptional responses to H/R stress might be explained by the relatively low levels of oxidative stress insufficient to trigger transcriptional upregulation of these proteases. Further studies of the ROS production, oxidative injury and expression of mitochondrial proteases in bivalve mitochondria during exposures to different stressors are required to test this hypothesis and shed light on the thresholds of mitochondrial injury required for induction of these protective mechanisms.
Mitophagy is a selective form of autophagy to degrade damaged mitochondria (Pickles et al., 2018). Based on the studies in the model organisms, the PINK1-Parkin pathways plays a key role in the stress-induced mitophagy (Jin et al., 2010). In healthy mitochondria, the constitutively expressed PINK1 is translocated to mitochondria in a membrane-potential-dependent manner, cleaved by proteases and released to cytosol for further degradation. Membrane depolarization in dysfunctional mitochondria causes stabilization and dimerization of PINK1 on mitochondrial surface (Ding and Yin, 2012). Bound PINK1 recruits and phosphorylates parkin (Dagda et al., 2009;Eiyama and Okamoto, 2015), which tags damaged mitochondria for further degradation by mitophagy (Narendra and Youle, 2011). Activity of the protein phosphatase PGAM5 also contributes to the activation of mitophagy Wu et al., 2014) and promotes the mitochondrial fission by activating the Drp1/Dnm1l protein (Wang et al., 2012). In our present study, hypoxia exposure induced elevated mRNA expression of prkn in M. edulis, while the transcript levels of the associated key genes pink1 and pgam5 rose during reoxygenation mitochondria. This expression pattern indicates an onset of mitophagy during reoxygenation as the PINK1-parkin pathway is increasingly required to tag dysfunctional mitochondria caused by oxidative stress. Transcriptional upregulation of PGAM5 also indicates activation of the cell death pathways (Wang et al., 2012) during reoxygenation in M. edulis, consistent with the earlier findings of upregulated apoptosis in this species during H/R stress (Falfushynska et al., 2020a).
The correlation network analysis showed strong correlation of mRNA expression of all three fission mediators, the PINK1parkin pathway regulators and oma1, underpinning the common regulatory mechanism for these pathways. This is also reflected by the clear separation of mRNA expression profiles of the mussels exposed to hypoxia and (especially) reoxygenation from the normoxic control groups along the 1 st principal component axis in the PCA associated with the high loadings of the fission-mitophagy-related transcripts. An earlier study showed a similar major shift of the mRNA expression profiles of the pro-apoptotic and pro-inflammatory genes in M. edulis exposed to the short-term H/R stress (Falfushynska et al., 2020a). The metabolic profiling in a related study showed onset of anaerobiosis, downregulation of gluconeogenesis and a shifted balance between protein synthesis and breakdown in M. edulis indicating early onset of severe stress response during H/R exposure (Haider et al., 2020).

Muted Transcriptional Response of the Mitochondrial Dynamics and Quality
Control Genes During H/R Stress in C. gigas Unlike M. edulis, C. gigas did not change its mRNA expression in any of the 17 studied genes in response to 24 h hypoxia and subsequent 1.5 h reoxygenation. Furthermore, the overlap of gene transcription found by PCA and the apparent tight correlation across all studied genes supports a conclusion of the lack of a transcriptional shift in the studied pathways during H/R stress in the Pacific oyster. A similar transcriptional stability was found for lonp1 expression after 6 days of hypoxia and 1-12 h of reoxygenation in the eastern oyster C. virginica (Ivanina et al., 2010). Previously, metabolic and mRNA expression profiling on apoptotic and pro-inflammatory genes in C. gigas revealed no changes in response to short-term hypoxia (Falfushynska et al., 2020a;Haider et al., 2020). Only 6 days exposure to hypoxia led to a small shift in metabolite profiles indicating metabolic rate depression, maintenance of balanced FAA/metabolite pool and effective detoxification by urea and purine cycles. Furthermore, maintenance of steady activity levels of mitoproteases during the H/R stress was earlier reported in C. gigas . Together with the muted mRNA expression regulation of mitochondrial quality control genes and mitochondrial dynamics, these findings indicate that 24 h of severe hypoxia might not be sufficient to induce strong stress markers in C. gigas. However, shortterm hypoxia (24 h, <1% O 2 ) and subsequent reoxygenation affected a reversible phosphorylation pattern of mitochondrial quality control proteins such as PGAM5 in C. gigas . This indicates that despite the lack of a strong transcriptional induction of the key genes in the mitochondrial quality control pathways, the activity of these pathways might be modulated by posttranslational protein modifications (PTM) in C. gigas. This hypothesis requires further investigation as the role of the PTM in the regulation of the mitochondrial function of marine bivalves has not been extensively studied Falfushynska et al., 2020b). Notably, a muted transcriptional response to short-term hypoxia and subsequent reoxygenation has been earlier shown for other important pathways involved in cellular bioenergetics of the oysters. Thus, hypoxia induced genes (especially respiratory chain genes) were regulated in C. gigas only after 7 days of exposure to 30% O 2 saturation, with no apparent changes at the earlier exposures (David et al., 2005). In the eastern oyster C. virginica, key marker genes in the mitochondrial electron chain and Krebs' cycle, as well as the cytosolic glycolytic genes showed stable expression during 6 days of severe hypoxia ( < 1% O 2 ) and subsequent reoxygenation (Ivanina et al., 2010). Like other hypoxiatolerant mollusks, the Eastern and Pacific oyster rely on metabolic rate depression to cope with oxygen deficiency (Storey and Storey, 2004) and use a reversible global transcription repression to arrest energy consuming transcription (Larade and Storey, 2002). Only specific transcription factors seem to be upregulated to allow transcription of certain hypoxic genes (David et al., 2005). These factors might contribute to the muted transcriptional response during H/R stress in the oysters. Overall, the hypoxia tolerant Pacific oyster seems to maintain its mitochondrial integrity and functions by reversible enzyme regulations sufficient for fast recovery in reoxygenation, which might contribute to the high stress tolerance of this species adapted to the habitat conditions of daily immersion and emersion (Sussarellu et al., 2013). Further investigations in a broader comparative framework (involving multiple species of bivalves with different levels of hypoxia tolerance) are needed to determine whether the hypoxia-tolerant mitochondrial phenotypes generally correlate with the transcriptional stability and stronger dependence on post-translational regulatory mechanisms during H/R stress in marine bivalves.

CONCLUSION
Our study showed that differential transcriptional regulation of the mitochondrial quality control pathways is involved in the molecular response to H/R stress in the blue mussels M. edulis but not in the oysters C. gigas. In the mussels, the H/R stress induced the pathways related to the mitochondrial fission, mitophagy and ATP-independent mitochondrial proteolysis. Interestingly, the transcripts of the proteins encoding ATPdependent mitochondrial quality control proteins (such as mitofusin 2, Lon protease and paraplegin) remained stable or downregulated during hypoxia, possibly as a consequence of the metabolic rate depression and an associated decrease in the ATP turnover rates. In C. gigas, no clear transcriptional response of the mitochondrial quality control pathways was found during H/R stress, with all studied transcripts maintained at control (normoxic) levels. This might indicate that unlike M. edulis, other mechanisms (such as post-transcriptional and post-translational regulation of the respective pathways) are involved in the regulation of the mitochondrial integrity in oysters. Alternatively, it is possible that the level of the mitochondrial damage in oysters did not reach the threshold required to trigger the respective transcriptional activation of the mitochondrial quality control pathways. The latter explanation appears plausible, given the relative stability of the metabolomics profile and the transcript levels of other stress genes (such as those involved in apoptosis and inflammation) during H/R stress in oysters compared with the mussels (Falfushynska et al., 2020a;Haider et al., 2020). Our present study and earlier published research (Falfushynska et al., 2020a;Haider et al., 2020) also indicate that the mRNA expression of the stress-related genes may not be a good biomarker of hypoxia-induced stress in an extremely tolerant species such as C. gigas. It is worth noting that the natural habitat conditions differ with regard to salinity (ocean vs. brackish, respectively) and the tidal height (intertidal vs. subtidal, respectively) between the two studied populations of the oysters and the mussels. Further studies are needed to determine whether the established links between the mRNA expression profiles of the stress-related mitochondrial genes and hypoxia tolerance reflect the species-specific differences between the oysters and the mussels, or some other aspects of the local adaptation to the divergent habitats in which these organisms were collected.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Zenodo: doi: 10.5281/ zenodo.4266408.

AUTHOR CONTRIBUTIONS
JS and HF conducted the experimental exposures and tissue collections with C. gigas, and M. edulis, respectively. mRNA extraction, cDNA synthesis and quantitative real-time PCR were done by JS and HF, and statistical analysis and data integration was done by JS. HF and JS designed primers according to transcripts provided by bioinformatical analysis from HP. JS, IS, and HP wrote the first draft of the manuscript, and all author contributed to the data interpretation and revisions of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was funded by the Deutsche Forschungsgemeinschaft within the project "MitoBOX: the mitochondrial basis of hypoxia tolerance in marine mollusks" (grant number 415984732). HP was supported by the NIH grant R21AG064479-01 and a Brain Health Research Institute Pilot Award from Kent State University.