Exploring the Sublethal Impacts of Cu and Zn on Daphnia magna: a transcriptomic perspective

Abstract

Metal contamination of aquatic environments remains a major concern due to their persistence. The water flea Daphnia magna is an important model species for metal toxicity studies and water quality assessment. However, most research has focused on physiological endpoints such as mortality, growth, and reproduction in laboratory settings, as well as neglected toxicogenomic responses. Copper (Cu) and zinc (Zn) are essential trace elements that play crucial roles in many biological processes, including iron metabolism, connective tissue formation, neurotransmitter synthesis, DNA synthesis, and immune function. Excess amounts of these metals result in deviations from homeostasis and may induce toxic responses. In this study, we analyzed Daphnia magna transcriptomic responses to IC₅ levels of Cu (120 µg/L) and Zn (300 µg/L) in environmental water obtained from a pristine lake with adjusted water hardness (150 mg/L CaCO₃). The study was carried out to gain insights into the Cu and Zn regulated stress response mechanisms in Daphnia magna at transcriptome level. A total of 2,688 and 3,080 genes were found to be differentially expressed (DEG) between the control and Cu and the control and Zn, respectively. There were 1,793 differentially expressed genes in common for both Cu and Zn, whereas the number of unique DEGs for Cu and Zn were 895 and 1,287, respectively. Gene ontology and KEGG pathways enrichment were carried out to identify the molecular functions and biological processes affected by metal exposures. In addition to well-known biomarkers, novel targets for metal toxicity screening at the genomic level were identified.

Introduction

Daphnia magna is a freshwater crustacean that has been widely used as a model organism in ecotoxicology due to its sensitivity to environmental stressors and ease of handling [1,2,3,4]. Copper (Cu) and zinc (Zn) are two heavy metals commonly found in aquatic environments as a result of anthropogenic activity. Cu is involved in numerous biological processes including antioxidant defense, muscle development, and immunological functions [5]. Cu also serves as a cofactor for oxidoreductases, oxygenases, hydroxylases, and transferases [6]. Zn is an integral component of several enzymes and plays a key role in growth, development, and immune functions [7]. Deficiencies in Cu and Zn have been linked to a range of health problems. Cu deficiency may lead to impaired growth, anemia, skeletal abnormalities, and compromised immune function [8]. Insufficient Zn intake has similarly been associated with impaired growth, as well as immune dysfunction, delayed sexual maturation, and impaired cognitive functions [9].

To ensure optimal health it is crucial to maintain a sufficient intake of Cu and Zn. However, excessive intake of these trace elements may also cause adverse effects such as metal toxicity and oxidative stress [10]. Daphnia magna has been extensively used to study the effects of these heavy metals in aquatic environments [3]. Exposure to Cu and Zn can induce a range of toxic effects in Daphnia magna, including impaired growth, reduced reproductive capacity, and premature death [11]. It is important to note that much of the research related to this topic has been performed using exposure mediums prepared in laboratory settings and has focused on physiological endpoints such as mortality and reproduction while neglecting the toxicogenomic endpoints [4, 12]. The need for transcriptomics in ecotoxicology is becoming increasingly important as traditional methods often fail to provide sufficient information about the mechanisms of toxicity and the effects of environmental stressors on living organisms [3, 13]. Transcriptomics provides a comprehensive view of gene expression changes and can identify key genes and pathways involved in stress responses [14]. Traditional methods of toxicity testing may thus provide information on the physiological and behavioral responses to the stressor, while transcriptomics provides a deeper understanding of the molecular mechanisms underlying these responses [11, 15, 16]. Furthermore, transcriptomics can aid in the identification of novel biomarkers of exposure and effect in aquatic organisms [16, 17]. Biomarkers are measurable responses at the sub-organism level that indicate exposure to stressors [18]. They can be used to monitor the health status of organisms in their natural environment [16, 19, 20]. Novel transcriptomic biomarkers identified from the exposure of Daphnia magna to heavy metals are promising tools for monitoring the effects of metal exposure [17, 21]. With this aim, we investigated the transcriptome of Daphnia magna exposed to environmental water samples with Cu and Zn in concentrations resulting in 5% immobilization (IC₅). The use of IC₅, rather than higher values like IC₁₀ or IC₂₀, or environmental concentrations of Cu or Zn, was employed to ensure that the effects being studied are subtle and do not cause overt harm to the organisms. This allows for a more nuanced understanding of the organism’s response to these stressors [22]. Studying sublethal exposure transcriptomics offers a better understanding of the molecular mechanisms underlying the effects of environmental stressors, which in turn aids in development of effective strategies for environmental monitoring and risk assessment and identification of safer and more sustainable alternatives to harmful chemicals and pollutants.

Materials and methods

Daphnia magna culture and exposure

Ephippia from the Daphtoxkit (MicroBioTests Inc., Belgium) were activated by rinsing with tap water and incubation for 72 h in standard freshwater (prepared using 67.75 mg/L NaHCO₃, 294 mg/L CaCl₂, 123.25 mg/L MgSO₄ and 5.75 mg/L KCl) under continuous illumination at a temperature of 22 ± 1 °C. Newly hatched juveniles (< 24 h) were transferred to the test waters directly after hatching and kept in 16-h light and 8-h dark cycle. Daphnia magna was fed a mixture of microalgae (Raphidocelis subcapitata) and yeast with adjusted ratios during metal exposure to ensure proper feeding.

Test waters used in the present study were prepared by adjusting the water hardness to 150 mg CaCO₃/L and the pH to 7 of the water obtained from the pristine Kiimatievanjärvi Lake (67°44′33.2"N, 22°10′49.5"E), Sweden. The concentration of potassium in the exposure waters was set to 3 mg/L with KHCO₃ to ensure proper maintenance for daphnids. The natural concentration of Cu and Zn in the lake water were 1 μg/L and 0.9 μg/L, respectively (Table 1). For the immobilization assay, Cu concentrations were set to 16, 32, 64, 128, 144, 160, 176, and 256 μg/L using CuCl₂, whereas Zn concentrations were set to 25, 50, 100, 200, 400, 500, 600, 700 and 800 using ZnCl₂. Control water without any addition of Cu or Zn was also used. Newly hatched Daphnia magna juveniles (< 24 h) were transferred and kept in 8 ml of exposure water in 6 well plates (Sarstedt, Germany). Each well contained 20 animals and four replicates were used for each exposure (total animals = 80). Immobility was screened daily using a light microscope. Daphnids not able to swim after gentle agitation were considered immobile. The percentage of immobilized daphnids at 96 h was plotted against the test concentrations. Slopes and curves were calculated to determine the IC₅ value with 95% confidence limits (p < 0.05) using concentrations transformed to logarithmic values with a base of 10.

Table 1 Composition of the water used in this study

For the transcriptomics assay daphnids were exposed to approximate IC₅ concentrations of 120 μg Cu/L and 300 μg Zn/L, along with control water. For each treatment, eight replicates containing a total of 35 daphnids were used for RNA extraction and subsequent transcriptomics analysis. In our earlier research [11], we noted that the onset of gene expression related to reproduction occurs around 96-h post-hatching (hph). This time point signifies the transition from the juvenile to the adult phase. Given the significance of reproduction in terms of physiology and its crucial role in ensuring the survival of a population, we aimed to gain insight into this vital physiological process. Following 96 h exposure, daphnids were collected and snap-frozen using liquid nitrogen and stored at -80 °C until RNA extraction.

RNA isolation for transcriptomics and qRT-PCR

Daphnia magna was lysed using 350 μL of TRI Reagent (Sigma) using tissue homogenizer (Precellys Evolution, Bertin Technologies, USA) and RNA extraction was performed using Direct-zol Kit (Zymo Research, USA) according to the manufacturer's instructions. A DeNovix DS-11 spectrophotometer (Wilmington DE, USA) was used to measure RNA concentration and purity (Table S1). cDNA libraries were generated from 100 ng of RNA. Library preparation was carried out by purification of poly-A containing mRNA, mRNA fragmentation, random primed cDNA synthesis, adapter ligation, and adapter specific PCR amplification by INVIEW Transcriptome Discover (Eurofins Genomics, Germany).

RNA sequencing and transcriptomics analysis

Illumina paired-end read sequencing was carried out to generate 2 × 150 bp reads. Reads were trimmed from three prime ends based on quality score. Average Phred score for the reads were determined for all replicates (Table 2). Reads were aligned to the Daphnia magna genome (daphmag2.4) using Hisat2 Aligner followed by quantification to alignment model using Partek genomic software (Partek Inc., St. Louis, USA). Data were normalized by counts per million (CPM) + 1.0E⁻⁴. For the gene expression analysis, we used DESeq2 [23]. This method allows us to identify genes that show statistically significant differences in expression levels between the groups. It uses a model based on the negative binomial distribution to test for differential expression in digital gene expression data. DESeq2 was chosen as it provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. This makes it highly suitable for high-throughput count data from RNA-Seq. FDR step-up value was used for follow-up multiple test correction. Transcripts with p-value ≤ 0.05 and FDR step-up value ≤ 0.05 were considered differentially expressed. A BLAST search was executed on the Daphnia magna genome, focusing on uncharacterized proteins that were among the top ten for fold change, both upregulated and downregulated using the default parameters. The selection of significant hits was based on e-values, with a focus on those with the highest identity and query coverage. After the initial transcriptomics analysis, Partek’s biomarker computing function was used to identify exposure-specific biomarker genes and to design primers for biomarker validation. Bubble plots were created to highlight GO annotations and KEGG pathway analysis for each exposure by employing identified DEGs using online data analysis and visualization platform (https://www.bioinformatics.com.cn/en).

Table 2 Samples were aligned to Daphnia magna reference genome by HISAT 2 algorithm and post alignment quality check for the read pairs were generated

qPCR validation

cDNA was synthesized using the qScript cDNA synthesis kit (Quanta Biosciences, USA) and 1000 ng of RNA obtained from 4 additional biological replicates, according to the manufacturer's instructions. The primer sequences of genes used for validation are shown in Table S2. qRT-PCR was performed to quantify the expression of the genes using qPCRBIO SyGreenMix Lo-ROX (PCR Biosystems, USA) using the CFX384 Real time PCR detection system (Bio-Rad, USA) with thermal cycling profiles of initial denaturation step at 95 °C for 2 min followed by 35 cycles of 95 °C for 5 s and 60 °C for 30 s. Expression ratios were calculated based on the ΔΔCt method [24]. Four different housekeeping genes were tested (actin, tubulin, elong, gapdh) and, ultimately, gapdh was selected based on the stability among the treatment groups.

Statistical analysis

All statistical analyses were performed using GraphPad Prism 8.0.2 software (GraphPad Software, USA) using one-way ANOVA followed by Dunnett’s post-test for multiple group comparison. Statistically significant differences were considered when p-values were < 0.05 (*p < 0.05, **p < 0.01 and ***p < 0.001). Principal component analysis (PCA) was used to analyze multivariate data using the SIMCA software, v13.0.3 (Umetrics, Sweden). A number of significant components were validated using cross validation rules. Tolerance ellipse based on Hotelling's T2 was used to check for outliers (95%). For enrichment analysis, terms and pathways with p-values less than 0.05 are considered significant.

Results

Determination of IC₅ levels

The results of the 96 h immobilization assay are presented in Fig. 1, along with the calculated IC₅ values for Cu (Fig. 1A) and Zn (Fig. 1B). The IC₅ values for Cu and Zn were determined to be 133 μg/L and 267 μg/L, respectively. These values represent the concentrations of Cu and Zn that would result in a 5% immobilization of the test organisms. The slopes and curves used to derive the IC₅ values provide an accurate representation of the dose–response relationship between the concentrations of Cu and Zn and the observed level of immobilization. Following the immobilization test, Daphnia magna neonates were treated with IC₅ concentrations (130 μg/L Cu and 300 μg/L Zn) as well as control water with parameters shown in Table 1 for transcriptomics analysis.

Fig. 1

Acute toxicity determination. Twenty Daphnia magna neonates (< 24 h old) per well were exposed to varying (A) Cu and (B) Zn concentrations for 96 h and the survival rate was recorded at the end of the exposure. Dose–response curves were generated to calculate IC₅ values. Mean ± SEM. N = 4

DEG identification

On average, 33 million reads were generated per sample and over 90 percent of the reads were aligned to Daphnia genome with a read quality of 35.73 Phred quality score (Table 2) reflecting the high accuracy of base calling. Principal component analysis displayed variance between the control and exposure groups. PCA explained 57.38% with three components (Fig. 2A). Control and exposure groups were separated clearly with the first principal component, whereas variance between the Cu and Zn exposed samples were less clear and mostly occurred on the secondary component axis.

Fig. 2

Transcriptomics analysis. A Principal component analysis (PCA) was conducted using transcripts with a Lowest Average Coverage (LAC) greater than one. The first component (PC1) accounted for 25.32% of the variance and distinctly segregated the control and exposure groups. B Venn diagram illustrates the differentially expressed genes (DEGs) between the control and exposure groups, with a significance level of p < 0.05 and a false discovery rate (FDR) less than 0.05. C and D Volcano plots were created to visualize the effects of Cu and Zn exposures

A total of 12,213 transcripts were identified with a lowest average coverage value greater than one and were used for DEG analysis. The Venn diagram was used to compare the DEGs identified between Cu and Zn exposures (Fig. 2B). There were 1,793 differentially expressed transcripts common to both Cu and Zn exposures. Although 1147 and 646 DEGs were identified for upregulation and downregulation, respectively, we did not identify any DEGs for either Cu or Zn that showed opposite expression (i.e., upregulated in one and downregulated in other, or vice versa). The Volcano plots visualize DEG distribution for both exposure groups (Fig. 2C, 2D) highlighting highly regulated annotated genes. Following the DEG analysis, genes were ranked based on their fold change to identify highly affected transcripts between control and exposure groups. The top ten transcripts for both upregulation and downregulation following Cu and Zn exposure are listed in Table 3 and Table 4, respectively. Uncharacterized proteins were blasted against the Daphnia magna genome at default settings and hit with significant e-values and highest identity and query coverage to identify the unknown proteins [25]. Aromatic-L-amino-acid decarboxylase (aadc) was the most upregulated gene for both Cu (57.85-fold) and Zn (93.9-fold) exposure. Additionally, complement factor D (cfd), putative Nose resistant to fluoxetine protein (nrf), leukocyte receptor cluster member 9 (leng9), and coagulation factor IX (fix) were identified in both exposure settings, as well as an unidentified protein with the locus tag of APZ42_015569. There was only one common transcript (APZ42_014255) among the top ten downregulated genes for both exposures (Cu -8.87-fold, Zn -31.21-fold), coding for an unidentified protein product. Identification of uniquely expressed genes following either Cu or Zn exposure revealed 650 and 833 upregulated as well as 245 and 454 downregulated genes for Cu and Zn, respectively.

Table 3 Top ten significantly upregulated and downregulated genes following Cu exposure ranked based on fold change and false discovery rate (FDR) from transcriptomics analysis

Table 4 Top ten significantly upregulated and downregulated genes following Zn exposure ranked based on fold change and false discovery rate (FDR) from transcriptomics analysis

Biomarker selection and qPCR validation

Using the Partek biomarker-computing function, a hierarchical cluster and a heatmap with unique biomarker genes for both Cu and Zn exposures was created (Fig. 3A). In addition to metallothionein (mt) homologs (known for their role in metal sequestration), a subset of the identified potential biomarkers underwent qPCR validation using a different set of biological replicates (Fig. 3B). All potential biomarkers were subsequently validated using qPCR, affirming their status as both unique and significant DEGs for their respective exposure groups. The Daphnia magna mt genes differed in their metal responses, with mt-a being upregulated by Zn exposure, whereas mt-b was upregulated by both Cu and Zn exposure, and mt-c remained unaffected by either metal exposure. The gene expression levels following Zn exposure were higher for mt-a (25-fold) as opposed to mt-b (threefold). The copy numbers obtained from transcriptomics and the relative expression results from qPCR were compared using Pearson's correlation analysis (Fig. 3C). An R-value of 0.96 and significant p-values indicate a strong correlation between the transcriptomics copy numbers and the qPCR relative expression results.

Fig. 3

Biomarker identification. A Heatmap and hierarchical clustering of potential biomarkers was generated using top ranked genes in biomarker analysis using Partek. B qPCR was carried out to to validate potential biomarker genes (Mean ± SEM. N = 8). C Linearity between RNA-Seq and qPCR reseults were analyzed using Pearson correlation coefficient

Gene set enrichment

We conducted a functional enrichment analysis on the differentially expressed genes within their respective groups. The analysis was performed using the STRING database [26]. The aim was to understand the impact of Cu and Zn exposure on different biological processes, molecular functions, and cellular components (Fig. 4A, Fig. 4B, and Table 5).

Fig. 4

Functional annotation of the Gene Ontology (GO). Enrichment was performed using Differentially Expressed Genes (DEGs) for each exposure. To prioritize the GO terms, an enrichment score was applied. The size of the dots in the visualization corresponds to the number of genes associated with the respective terms. The color scheme applied to these dots reflects the percentage of genes upregulated or downregulated within the specified GO term

Table 5 Biological pathways in Daphnia magna enriched by upregulated genes after copper and zinc exposure. Related terms were summarized by the top- and lowest-level pathways/events

The cellular biogenic amine metabolic process emerged as the most enriched Gene Ontology (GO) term in the biological process category for both exposures. A significant proportion of the genes associated with this term exhibited upregulation for both exposures (Fig. 5A). Cu exposure led to significant enrichment in molecular processes such as cellular oxidant detoxification, response to oxidative stress, and lipid metabolic process. Similarly, Zn exposure influenced biological processes such as lipid metabolic process and cuticle development. Both Cu and Zn exposures significantly upregulated ribosomal subunit biogenesis, underscoring their importance in cellular homeostasis. Cellular compartments such as the extracellular space, lysosome, and extracellular region were significantly influenced by Cu exposure, with genes related to the extracellular space showing noticeable downregulation. Zn exposure also impacted the extracellular space associated genes with a similar degree of downregulation. However, cellular compartments such as the lysosome and lytic vacuoles were exclusively enriched with Cu exposure. Moreover, Cu exposure affected various molecular functions such as peroxidase activity, antioxidant activity, and heme binding. Genes associated with both peroxidase and antioxidant activity showed substantial downregulation (Fig. 5B). In contrast, Zn exposure had a more distinct effect, with significant alterations in metallocarboxypeptidase activity and carboxypeptidase activity. All genes associated with metallocarboxypeptidase activity displayed downregulation (Fig. 5C).

Fig. 5

Expression level of genes in highest enriched pathways. A Cellular biogenic amine process was the most enriched biological function for both exposures. B Peroxidase activity and (C) Metallocarboxypeptidase activity were also further identified as the most enriched molecular functions for Cu and Zn, respectively

The KEGG pathway analyses revealed differential gene expression profiles elicited by exposure to Cu and Zn across a variety of biological pathways (Fig. 6). Upregulation of genes in response to Cu exposure was observed predominantly within pathways associated with genetic information processing. These encompassed processes, including RNA transport, spliceosome activity, and protein processing within the endoplasmic reticulum. Notably, certain cellular processes, including lysosome function, were also enriched, indicating a broad cellular response to Cu exposure. Conversely, Zn exposure resulted in a significant shift in gene expression within metabolic pathways. This was particularly evident in pathways related to arachidonic acid metabolism and fatty acid metabolism. Moreover, ribosome biogenesis in the eukaryote’s pathway (categorized under genetic information processing) showed significant upregulation of genes in response to both exposures.

Fig. 6

KEGG analysis and identified affected pathways. Enriched pathways are categorized according to their functional relationships using hierarchical terms. Additionally, the number of upregulated or downregulated genes within each pathway is indicated within the corresponding bars

Using the STRING database, 94 significantly enriched Reactome pathways were identified as a result of Cu exposure. Similarly, 6 significantly enriched Reactome pathways were identified due to Zn exposure (Table 5). Among the 94 pathways enriched by Cu exposure, 93 were predominated by upregulated genes. These pathways are related to oxidative stress and DNA damage, altered cell cycle progression and DNA replication, and impaired gene expression and RNA metabolism. Oxygen-dependent proline hydroxylation of Hypoxia-inducible Factor Alpha and regulation of ornithine decarboxylase suggest that Cu exposure may induce oxidative stress in Daphnia magna. This was also evident by pathways involved in DNA replication and the regulation of the cell cycle checkpoints, which can arrest the cell cycle in response to DNA damage or other stress signals. Cu exposure also enriched the pathways associated with activation and signaling of nuclear factor kappa-B (NF-κB) in B cells and cross-presentation of soluble exogenous antigens (endosomes). These pathways were associated with activation of inflammatory and adaptive immune responses. The interleukin-1 signaling, which is associated with production of ceruloplasmin, a major Cu store was also enriched. On the other hand, genes upregulated with Zn exposure enriched tyrosine metabolism, which involves the biosynthesis and degradation of tyrosine and its derivatives. The major pathway of rRNA processing in the nucleolus and cytosol was the only pathway affected by overexpressed genes following Cu and Zn exposure. Developmental processes and Choline catabolism were overrepresented by down-regulated genes in Cu and Zn exposures, respectively.

Discussion

Daphnia magna is a commonly used model organism for toxicity testing due to its ecological relevance and sensitivity to contaminants such as Cu and Zn [3]. The OECD 202 immobilization assay is used to assess the acute toxicity of contaminants on Daphnia magna [27]. In this study, an extended 96 h immobilization assay was used to evaluate the effects of Cu and Zn on the mobility of Daphnia magna. As both Cu and Zn are essential metals and required for several biological processes, exposure concentrations of less than 100 μg/L did not result in any significant immobilization and was in line with previous reports [28]. Cu concentrations of 133 μg/L and higher resulted in significant immobilization of the test organisms compared to the control groups. The concentration–response curve for Cu showed a clear dose–response relationship, with increasing concentrations of Cu resulting in a proportional increase in immobilization. This was consistent with previous studies that have demonstrated the acute toxicity of Cu on aquatic organisms, including Daphnia magna [4, 11, 29, 30]. The calculated IC₅ value for Cu indicates that Cu contamination at concentrations just above 100 μg/L can have a significant impact on the mobility of aquatic organisms. This suggests that the distinction between the homeostatic boundaries of a healthy environment and a toxic one lies within a very narrow range of Cu concentrations. Therefore, the IC₅ value is an important metric for establishing regulatory limits for contaminants in aquatic ecosystems and for assessing the risk of toxicity to aquatic organisms with established standards by regulating bodies [31].

The concentration–response curve for Zn indicated a lower sensitivity of Daphnia magna to Zn than to Cu. Although Zn is an essential micronutrient for living organisms, excessive accumulation in the environment can cause severe toxicity to aquatic organisms, including the freshwater planktonic crustacean Daphnia magna The IC₅ of Zn was determined to be two times more than the IC₅ of Cu after 96 h exposure. The results show that Zn contamination has a lower impact on the mobility of aquatic organisms as compared to Cu. These results are in line with previous research using genus Daphnia where Zn has been shown to be two to ten times more toxic than Cu [32]. The IC₅ value used in this study also has practical applications in setting regulatory limits in rivers and other streams found downstream of mining operations where effluents increase the concentrations of these metals. We have previously shown the effects of streams contaminated with effluents containing different metals from decommissioned mining sites on several organisms including Daphnia magna [4, 33, 34].

To further understand the underlying mechanisms of this differential toxicity, transcriptomic analysis of Daphnia magna has allowed researchers to investigate the molecular responses to Cu and Zn exposure. Several studies have identified differentially expressed genes (DEGs) in Daphnia magna exposed to Cu and Zn. This includes genes involved in oxidative stress response, metal ion binding, and detoxification. These genes were upregulated in response to Cu and Zn exposure as observed following quantitative real-time PCR (qPCR) validation of the gene expression [35, 36].

In the present study, aadc was identified as the most upregulated gene for both Cu and Zn exposure. It is responsible for the decarboxylation step by converting L-3,4-dihydroxyphenylalanine (L-DOPA) and L-5-hydroxytryptophan (5-HTP) into dopamine and serotonin, respectively in both the catecholamine and the indolamine synthetic pathways, which could result in decreased response to stimuli and decreased mobility [37,38,39,40]. Previous studies have shown neurotoxicity effect of nickel, lead, and mercury in humans and rat pheochromocytoma cell line (PC12) through disruption of dopamine synthesis and transport [41, 42]. The relationship between Zn and aadc has not been reported in previous studies. The aadc expression is upregulated during the innate immune response in Drosophila melanogaster [43]. It is also expressed in a wide range of invertebrates, where it plays a critical role in the biosynthesis of biogenic amines, including dopamine, serotonin, and octopamine [37, 40]. These neurotransmitters are involved in regulating physiological processes in invertebrates, including behavior, development, and reproduction. When metal homeostasis is disrupted by excess Cu or Zn, it can cause neurotoxicity in marine invertebrates by inducing oxidative stress, disrupting calcium signaling, and impairing mitochondrial function through the production of reactive oxygen species (ROS) and interaction with sulfhydryl chemical groups [44, 38, 45,46,47]. The nrf is another common gene that was upregulated significantly by both exposures. NRF has been identified in Caenorhabditis elegans to produce resistance to nose contractions caused by the widely used antidepressant, fluoxetine [48]. Atli and Sevgiler reported increased oxidative stress in Daphnia magna exposed to the antidepressant fluoxetine (20–41 µg/L) and Zn (40–80 µg/L) [49]. Although the regulation of both aadc and nrf in metal induced toxicity are unclear, the present results suggest that they are regulated by metal exposures in Daphnia magna.

Exposure to Zn upregulated the expression of several newly identified biomarker genes such as clca2, enpp4, znt1, and mt-a. These genes may play important roles in maintaining homeostasis within cells exposed to excess Zn, either by excreting Zn ions from cells or by reducing influx through membrane transporters [50] clca2 encodes a calcium-activated chloride channel that allows the transport of chloride and other anions into the cell to maintain the electrochemical balance across the cell membrane [51]. These chloride channels are known to be involved in the secretion of electrolytes and water in epithelial cells, as well as regulating smooth muscle contraction and neuronal excitability in sensory neurons [52, 53]. Bouron and Oberwinkler argued that such channels may also take part in plasmalemmal transport of Zn in Daphnia magna [54]. enpp4 encodes an ectonucleotide pyrophosphatase/phosphodiesterase belonging to the Phospho-/Sulfo-coordinating Metalloenzyme superfamily, which require Zn ions to produce nucleoside-monophosphates from nucleotides and their derivatives [55]. Increased expression of enpp4 may increase bone mineralization and vascular smooth muscle calcification through hydrolysis of inorganic pyrophosphate, as well as generate extracellular nucleosides for cellular uptake, as seen in mouse [50, 56]. This may lead to immobilization of the animals. In Xenopus laevis and rainbow trout, znt1 encodes a Zn transporter that pumps Zn ions out of the cytoplasm to reduce the Zn concentration in cells [57, 58]. Several studies have reported similar findings, showing that exposure to Zn can upregulate the expression of genes involved in Zn transport and detoxification in marine copepods and fish [57,58,59,60]. Taken together, the upregulation of clca2, enpp4 and znt1 could indicate an increased cellular activity induced by excess intracellular Zn which will ultimately be sequestrated by metallothionines (MT), mainly MT-A [4, 61].

Daphnia magna has a sophisticated detoxification system to deal with the toxicity of heavy metals, which includes the synthesis of MT proteins [62]. MTs are a group of low-molecular-weight, cysteine-rich proteins that bind to heavy metals and play a critical role in the detoxification process [63]. MTs not only bind to metals for the purpose of excretion but also serve as a Zn pool that can be used in the event of deficiency [64]. Daphnia magna has three MT isoforms, namely mt-a, mt-b, and mt-c [65]. Although the physiological roles of these MT isoforms during Zn and Cu exposure in Daphnia magna have been studied in several experiments, the specific functions of each isoform are still not fully understood [66]. Upregulation of the mt-a gene indicates that Daphnia magna actively tries to sequester excess Zn. While the expression of mt-a was induced by Zn exposure, mt-b expression was induced by both Zn and Cu exposure, albeit to a lesser extent than mt-a induction by Zn. In contrast, mt-c expression was not altered by either exposure. The higher response of mt-a to Zn than Cu is consistent with previous research on MTs [61]. Furthermore, mt-a appears to be more responsive to Zn exposure than mt-b and mt-c, as its expression is induced at lower Zn concentrations. In contrast to our results where only mt-b, and not mt-c, was upregulated by Cu, others have observed that both mt-b and mt-c expression are induced at higher Cu concentrations [67]. Thus, mt-a appears to be the main isoform that sequesters Zn while mt-b is induced to a lesser extent by both Cu and Zn. Unlike Cu, the redox-inert nature of Zn supports an evolutionary conservation of Zn binding sites in metalloproteins specializing in sequestering Zn only [68]. It has been suggested that this specificity might be a consequence of necessarily strict control of Zn levels since it is required for a vast number of proteins [69]. The differential expression of mt-a and mt-b in response to Zn and Cu exposure suggests that they are regulated differently, which might be a result of having different metal transcription factor binding sites in their enhancer regions [70, 71]

The expression levels of mthfr and tsl genes were found to be upregulated when Daphnia magna was exposed to IC₅ Cu. MTHFR and TSL are important components in DNA synthesis [72]. The mthfr encodes for methylenetetrahydrofolate reductase, which plays a crucial role in folate metabolism [73, 74] This enzyme is involved in the conversion of 5,10-methylenetetrahydrofolate to 5-methyltetrahydrofolate [73], a reaction that is necessary for the synthesis of methionine from homocysteine. On the other hand, tsl encodes for thymidylate synthase-like protein, which is involved in folate metabolism during DNA synthesis [75]. Thymidylate synthase is an essential enzyme for DNA synthesis as it provides the sole intracellular de novo source of dTMP (deoxythymidine monophosphate), which is subsequently phosphorylated to dTTP for DNA replication [75]. The upregulation of both genes could suggest an increased demand for DNA synthesis or repair in Daphnia magna when exposed to IC₅ Cu. Intracellular accumulation of Cu, Zn, and other heavy metals and their toxicity through increased production of ROS (reactive oxygen species) coupled with reduced antioxidant activity, DNA damage and inhibition of relevant repair mechanisms, and protein misfolding disorders in intertidal copepod, Daphnia magna and other eukaryotic organisms are well established [76,77,78,79].

The expression levels of metalloendopeptidase-like protein (mep) and mucin (muc) genes were downregulated in Daphnia magna following exposure to IC₅ Cu. These genes are involved in the degradation of extracellular matrix proteins, which are a major component of the mucus layer in the gastrointestinal tract [80,81,82]. The observed downregulation of mep and muc genes indicates a potential decrease in both the degradation and production of mucus within the gastrointestinal tract of Daphnia magna [83]. This could potentially lead to changes in the organism’s absorption and digestion of nutrients, as well as the ability to protect itself from harmful substances in the environment [84,85,86].

Biogenic amines are low-molecular-weight, organic nitrogen compounds formed either by the decarboxylation of amino acids or by the amination and transamination of aldehydes and ketones during normal metabolic processes in living cells [87]. They participate in signaling as neurotransmitters and hormones, or as components of vitamins, phospholipids, and ribosomes within cells. Biogenic amines also play a number of important roles in physiological processes, from cell proliferation and differentiation to nutrition, immune response, and neurobiology and reproduction [88]. AADCs catalyze the conversion of aromatic L-amino acids into aromatic monoamines [89]. S-adenosylmethionine decarboxylase (AMD) is an essential regulatory component of the polyamine biosynthetic pathway by generating the n-propylamine residue required for the synthesis of spermidine and spermine from putrescin [90]. The upregulation of these central genes indicates a potential need for biogenic amines in the presence of these metals, possibly due to their essential roles in cellular development, metabolism, and ribosomal function.

Peroxidases are enzymes that typically catalyze a reaction where a peroxide (ROOR’) is reduced to water (ROH) or alcohol (R’OH) using an electron donor. This reaction is crucial in lignification, suberization, auxin catabolism, and the response to environmental stresses like wounding, pathogen attack, and oxidative stress [91]. In this study, many genes related to peroxidase activity were downregulated in the presence of both Cu and Zn, possibly as a response to oxidative stress. Cu exposure has been reported to induce changes in the expression of glutathione-related genes, which play a crucial role in protecting cells from oxidative stress [92]. Peroxidases participate in the detoxification of reactive oxygen species (ROS), and their activity is modulated in response to changes in ROS levels.

Cu exposure resulted in the regulation of detoxification processes and several molecular functions, while Zn exposure demonstrated a mix of up- and downregulation in several processes and functions. Cu is a component of enzymes involved in energy production, iron metabolism, neurotransmitter biosynthesis, and connective tissue formation [93]. Cu is also needed for lysosomal, waste removal and recycling within the cell [93]. This role of Cu could explain enrichment in terms related to lysosomes and lytic vacuoles in the presence of Cu. Moreover, signaling and activation of the NF-κB, and cross-presentation of soluble exogenous antigens (endosomes), could be associated with the role of Cu in enhancing immune activity of macrophages, a phenomenon known as immune priming [94, 95]. Cu is also involved in cellular processes both at the cell surface and in cell signaling, which could explain the enrichment of terms related to the plasma membrane [96]. The apical part of a cell is a specialized region involved in secretion and absorption that may be influenced by the presence of Cu. Zn is essential for enzyme function, protein structure, and gene regulation. However, its roles may differ from those of Cu, leading to the enrichment of different cellular components as seen in the functional enrichment results.

KEGG pathway enrichment resulted in distinct pathways for each exposure. Cell toxicity may occur due to high concentrations of Cu, either from environmental intake or abnormal accumulation within cells as a result of genetic mutations [97]. Our results show significant upregulation in genetic information processing pathways, including cell cycle checkpoints, DNA replication, RNA transport, and protein processing in the endoplasmic reticulum, which could be indicative of the cellular response to mitigate Cu induced stress [77, 78]. However, the absence of pathways known to increase the expression of antioxidant and anti-inflammatory enzymes indicates that the oxidative stress may be uncontrolled and potentially lead to cellular damage. This is also supported by enrichment of the interleukin-1 signaling, which induces the major Cu sequestrating protein, ceruloplasmin [98]. On the other hand, Zn is known to influence gene expression through both direct and indirect interactions. Direct interactions involve the binding of Zn to transcription factors, thereby altering the rate of transcription, while indirect interactions may involve signal transduction systems and hormonal or cytokine stimuli [99]. We observed an upregulation in genes associated with metabolic pathways following Zn exposure which indicates the activation of certain regulatory mechanisms. Evidence for this regulation is provided by the enrichment of the tyrosine metabolism pathway. Tyrosine metabolism pathway produces various bioactive molecules from Tyrosine, which serve multiple functions. They regulate antioxidant defense, act as chelators of metal ions (including Zn) and modulate the availability and toxicity of these ions [100, 101]. Moreover, Zn has been shown to affect the activity and expression of antioxidant enzymes, further emphasizing the role of Zn in detoxification processes [102]. Eukaryotic cells have been shown to rapidly modulate ribosome biogenesis in response to oxidative stress, DNA damage, and alterations in amino acid levels [103]. The upregulation of genes within the ‘ribosome biogenesis in eukaryotes’ pathway, observed following both exposures, indicates a stress response. Enrichment of ribosome-related pathways has been reported after exposure to metal stress, carbamates, oxidative stress and microcystis exposure in Daphnia pulex and Daphnia magna [104,105,106]. Under stress conditions the rRNA processing is inhibited and the unprocessed rRNA is stored in the nucleolus until the stress is resolved. This allows the cells to conserve energy and maintain the integrity of the nucleolus, which is a key organelle for ribosome biogenesis and stress response [107]. Activation of the stimuli-sensing channels may activate the Ca²⁺ dependent endonuclease, which can cleave the rRNA and induce nucleolar fragmentation [108]. This may result in reduced ribosome biogenesis and protein synthesis as well as increased apoptosis.

The transcriptomics analysis of Daphnia magna exposed to IC₅ concentrations of Cu and Zn revealed significant pathway enrichment that has direct implications for both individual organisms and population dynamics. Collectively, these molecular disturbances can decrease population resilience, reduce reproductive success, and alter population dynamics, ultimately impacting ecosystem stability [77, 104]. Demonstrating significant pathway enrichment at the IC₅ level is crucial as it highlights the initial molecular disruptions that serve as early indicators of ensuing physiological and phenotypic changes, providing a predictive model for more severe impacts at higher concentrations, and offering valuable insights for regulatory frameworks to preemptively mitigate environmental and ecological risks.

Overall, transcriptomics analysis provided valuable insights into the molecular mechanisms underlying the responses of Daphnia magna to Cu and Zn exposure. In the present study we identify key genes and pathways involved in stress responses, providing a basis for further investigations into the effects of heavy metal exposure on aquatic organisms. In addition, we have identified IC₅ values for Daphnia magna, which increases the resolution of our understanding of the organism’s response to environmental stressors.

In conclusion, our findings provide valuable insights into the distinct cellular responses invoked by Cu and Zn exposure. This approach not only expands our comprehension of the impact of these metals on Daphnia magna, but also establishes a robust foundation for prospective investigations into this key aquatic organism.