Skip to main content

Transcriptome analysis provides insights into copper toxicology in piebald naked carp (Gymnocypris eckloni)

Abstract

Background

Copper was used for many years in aquaculture operations as an effective algaecide or a parasite treatment of fish. It is an essential nutrient with numerous functions in organisms, but is toxic at high concentrations. However, the toxicity of copper to fish remains unclear. In this study, we used the piebald naked carp, Gymnocypris eckloni, as a model. RNA-seq data from different tissues, including gills, kidney, and liver, were used to investigate the underlying mechanism of copper toxicology in G. eckloni.

Results

We compared the transcriptomes from different tissues with different time durations of copper ion treatment. After 72 h copper ion treatment, the number of genes with different expression in gills and liver changed dramatically, but not in kidneys. In KEGG functional enrichment, the pattern of differentially expressed genes (DEGs) was also similar in the gills and liver. The most enriched pathway of DEGs was “Ribosome” in both tissues. Furthermore, we analyzed the expression levels of genes involved in oxidative stress response and protein synthesis using qPCR and RNA-seq data. Our results showed that several genes involved in oxidative stress response were up-regulated both in gills and liver. Up-regulation of these genes indicated that copper treatment caused oxidative stress, which is likely to result in ribosome damage. In addition, our results showed that the expression of Eef1b2, a transcription elongation factor, was decreased in the liver under oxidative stress, and the expression of translation initiation factors Eif4ebp1 and eIF2α, and elongation factor eEF2 was up-regulated. These results supported the idea that oxidative stress inhibits protein synthesis in cells.

Conclusions

Our results indicate that copper exposure caused different responses in different tissues, since the gene expression patterns changed substantially either in the gills or liver, while the effect on the kidney was relatively weak. Furthermore, our results indicated that the expression pattern of the genes involved in the ribosome, which is a complex molecular machine orchestrating protein synthesis in the cell, together with translation initiation factor and elongation factors, were affected by copper exposure both in the gills and liver of piebald naked carp. This result leads us to speculate that the downregulation of global protein synthesis is an acute response strategy of fish to metal-induced oxidative stress. Moreover, we speculate that this strategy not only exists in the selective translation of proteins but also exists in the specific translation of functional proteins in tissues and cells.

Background

Increasing global contamination of aquatic systems is a critical environmental problem. In addition to cadmium, copper is one of the most common pollutants in water. Copper is a basic nutrient with multiple functions in living organisms [1], but it is toxic if its concentration is high [2,3,4]. It has been used for many years as an effective algaecide or a parasite treatment of fish in aquaculture operations. However, the ambient copper can be absorbed by the fish through the gills, skin, and gastrointestinal tract and transferred through blood to the internal organs [5, 6]. The accumulation of copper causes damage to the homeostasis of essential metals in different organs. Oxidative stress is one of the known mechanisms of copper toxicity to fish [5, 6].

Oxidative stress can be caused by metal ions in model organisms, and biochemical ways to prevent oxidative damage are particularly significant in toxicology. Some fish species, such as zebrafish (Danio rerio), goldfish (Carassius auratus), and brown trout (Salmo trutta), have been well studied as models of the effects of environmental stresses and metal pollutants on oxidative stress and antioxidant resistance in freshwater fish [7,8,9]. In previous studies, the effects of metal contamination on fish were investigated by measuring the activity of superoxide dismutase (SOD), which is involved in the antioxidant response, and glucose-6-phosphate dehydrogenase (G6PDH) activity which participates in energy accumulation but also in anti-oxidant responses [10]. RNA-seq is an effective method for comprehensively identifying the molecular pathways affected by metal ion exposure, but so far, few studies have been conducted using the technique to study metal toxicology in fish [7, 11], especially many aspects of fish copper poisoning, which remain unclear.

In this study, we use the piebald naked carp (Gymnocypris eckloni) as a model. It is a representative species of the subfamily Schizothoracinae, and is the dominant fish in the upper reaches of the Yellow River at an elevation of 3000 m [12, 13]. G. eckloni is economically valuable due to its high production capacity. As mentioned above, copper is added to commercial aquaculture or exists in metal-polluted water. However, the effect of copper on the physiological function of different tissues in G. eckloni has not been evaluated, but is an important issue in food safety and aquatic toxicology.

In the present study, we used RNA-seq to investigate the response of different G. eckloni tissues to a copper concentration of 0.01 mg/L. This concentration is ecologically relevant since it is comparable to those found in water or used in commercial aquaculture [14,15,16]. Transcriptome data from three different tissues – gills, liver, and kidney – were used to investigate the response to copper exposure. We compared transcriptome changes in these tissues to assess the effects of copper ions on different tissues. In addition, by analyzing the changes in gene expression levels related to oxidative stress and the protein synthesis pathway, we proposed the potential molecular mechanism of copper toxicity in the piebald naked carp.

Results

RNA-seq data and de novo transcriptome assembly

To systematically investigate the transcriptome dynamics of the gills, liver, and kidney at different time points after copper ion treatment, we obtained 36 transcriptome samples. After removing low-quality sequences and trimming adapter sequences, 7.18–9.36 Gb of 125-bp paired-end clean reads were generated from each library. The total counts of clean reads were between 47,873,998 and 62,386,764.

Since there is no published genome of this species or a closely related species, we assembled the transcripts de novo and used this assembly as a reference. De novo assembly using all clean reads identified 704,405 transcript fragments, ranging in length from 201 to 23,517 bp, with an average length of 774 bp. The N50 (the length of the contig at 50% of the assembly arranged in descending order according to contig length) and N90 (the length of the contig at 90% of the assembly arranged in descending order according to contig length) values of the obtained transcripts were 1265 and 298 bp, respectively. In further assembly analysis, 486,221 gene fragments were identified across all transcripts, with an average length of 994 bp. The N50 and N90 values of the obtained gene fragments were 1480 and 438 bp, respectively.

The number of transcripts obtained in this work is similar to many published transcriptome assemblies [13, 17,18,19,20,21], including de novo transcriptome assemblies in several teleost species [13, 20, 21]. This indicates that our assembly quality is reliable to be used for subsequent transcriptome data analysis.

Gene functional annotation

All 486,221 gene fragments were annotated against seven public databases: NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences (Nt), Protein family (Pfam), Clusters of Orthologous Groups of proteins/eukaryotic Ortholog Groups (KOG/COG), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Ontology (GO). A total of 152,017 gene fragments were annotated against the Nr protein database (31.26%). The most common species with high sequence similarity in the Nr database was Danio rerio (64.8%), followed by Astyanax mexicanus (5.7%), Clupea harengus (2.8%), Ichthyophthirius multifiliis (2.1%), and Oncorhynchus mykiss (2.0%). Most of these species are teleost fish, suggesting the gene fragments used in this study were correctly assembled and annotated.

After GO annotation, the predicted G. eckloni genes were further cataloged in the next subdirectory under the three categories (Biology Process, Cellular Component, and Molecular Function). A total of 56 functional groups are displayed in Fig. 1.

Fig. 1
figure 1

A total of 56 functional groups were identified by GO annotation under the three categories (Biology Process, Cellular Component, and Molecular Function)

KOG annotation assigned the 55,116 genes into 26 groups, which are shown in Fig. 2. Among these categories, the most sequences were assigned to “Signal transduction mechanisms” (10,376 genes), followed by “General function prediction only” (8966 genes), “Posttranslational modification, protein turnover, chaperones” (5602 genes), “Intracellular trafficking, secretion, and vesicular transport” (4198 genes), and “Transcription” (3644 genes).

Fig. 2
figure 2

Using the KOG database, a total of 55,116 genes were annotated, which were categorized into 26 groups. The vertical axis represents the proportion of the number of gene fragments in each functional category

A total of 50,669 genes were annotated using the KEGG database (Fig. 3). These genes were further divided into five groups: A, Cellular Processes (9905 genes); B, Environmental Information Processing (11,646 genes); C, Genetic Information Processing (6162 genes); D, Metabolism (11,077 genes); and E, Organismal Systems (20,366 genes).

Fig. 3
figure 3

Results of KEGG annotation. a, Cellular Processes (9905 genes); b, Environmental Information Processing (11,646 genes); c, Genetic Information Processing (6162 genes); d, Metabolism (11,077 genes); and e, Organismal Systems (20,366 genes)

Differentially expressed genes and functional enrichment

The responses of different tissues of the piebald naked carp to the copper ion treatment were different. In gills, after 6 h of metal ion treatment, 410 and 361 genes were respectively up- and down-regulated. After 36 h of treatment, the numbers of up and down-regulated genes were 327 and 324, respectively. The number of differentially expressed genes (DEGs) was abruptly increased after 72 h of treatment; 2238 and 609 genes were up- and down-regulated, respectively (Fig. 4a). The trend of gene expression patterns in the liver was similar to the gills. The period with the most differences in gene expression level occurred after 72 h of metal ion treatment, with 401 and 555 genes up- and down-regulated, respectively. The respective numbers were 285 and 246 after 6 h, and 325 and 553 after 36 h (Fig. 5a). The numbers of DEGs were relatively constant in different periods in the kidney. After 6 h of copper ion treatment, 422 and 326 genes were up- and down-regulated, respectively. The respective numbers were 332 and 357 after 36 h, and 295 and 255 after 72 h.

Fig. 4
figure 4

The differentially expressed genes (a), functional enrichment (b), and length distribution (c) in the gills after 72 h of copper ion treatment compared with control

Fig. 5
figure 5

The differentially expressed genes (a), functional enrichment (b), and length distribution (c) in the liver after 72 h of copper ion treatment compared with control

To further explore the mechanism of damage of metal ions to organs and the pathways of responses of the organs, we analyzed the functional enrichment of the DEGs. In GO functional enrichment in gills and liver, the most DEGs were enriched in “metabolic process” (Fig. 4b and 5b) under the Biological Process category. In the gills, the DEGs were enriched in the terms “organic substance metabolic process”, “primary metabolic process” under the Biological Process category; “cell” and “cell part” under the Cellular Component category; and “catalytic activity” under the Molecular Function category. In the liver, the DEGs were enriched in the terms “organonitrogen compound metabolic process”, “cellular component biogenesis” under the Biological Process category; “intracellular” and “intracellular part” under the Cellular Component category; and “heterocyclic compound binding” and “organic cyclic compound binding” under the Molecular Function category.

De novo transcriptome assembly often results in a large number of short contigs. To evaluate the potential impact of short contigs on the identification of DEGs, we plotted the length distribution of DEGs (Fig. 4c and 5c). The results indicated that the proportion of short contigs (< 500 bp) in the DEGs was minor (6.81% in gills and 6.49% in the liver).

In KEGG functional enrichment, the pattern of DEGs was also similar in the gills and liver. The most enriched pathway of DEGs was “Ribosome” in both tissues. The “Protein processing in endoplasmic reticulum”, “Oxidative phosphorylation”, “Lysosome”, and “Glycolysis/Gluconeogenesis” were also enriched in both tissues (Supplementary file 1 and 2) (Fig. 6).

Fig. 6
figure 6

KEGG enrichment of DEGs in gills (a) and livers (b). The DEGs are from the 72 h copper treatment group vs. the control group

Genes involved in oxidative stress and protein synthesis after copper exposure

The DEGs of both gills and liver were mostly enriched in “Ribosome” and oxidative stress response in the copper treated group as compared with control. Therefore, the expression levels of genes involved in oxidative stress response and protein synthesis were measured using qPCR. Our results showed that the expression levels of these genes displayed different patterns after copper treatment. Among the genes involved in oxidative stress, the expression of Hmox1 (Cluster-47,180.157989), Nqo1 (Cluster-47,180.117564), Gpx1 (Cluster-47,774.0), Cu/Zn-SOD (Cluster-47,180.196300), Prdx1 (Cluster-47,180.233418), Gadd45b1 (Cluster-47,180.146985), Hspa9 (Cluster-47,180.10385), Hspa14 (Cluster-47,180.185583), and Mn-SOD (Cluster-47,180.180330) was greatly up-regulated in the gills (Fig. 7), while the expression of Nqo1 (Cluster-47,180.117564), Cu/Zn-SOD (Cluster-47,180.117564), Prdx1 (Cluster-47,180.233418), Gadd45b1 (Cluster-47,180.233418), Hspa14 (Cluster-47,180.185583), and Mn-SOD (Cluster-47,180.180330) was greatly up-regulated in the liver (Fig. 8). Among the genes involved in protein synthesis, the expression levels of translation initiation factors Eif4ebp1 (Cluster-47,180.374620) and eIF2α (Cluster-47,180.297183), and elongation factors eEF2 (Cluster-47,180.277438) and Eef1b2 (Cluster-47,453.0), were up-regulated (Figs. 7 and 8). In addition, the expression of Hmox1 (Cluster-47,180.157989), Gpx1 (Cluster-47,774.0), Eef1b2 (Cluster-47,453.0), and Hspa8b (Cluster-47,180.176206) was down-regulated in the liver (Fig. 8). In summary, it can be seen that some transcriptome data showed large differences among duplicate samples, but there was general agreement between transcriptome results for whole heavy metal stress and qPCR validation tests, and transcripts basically represented true expression differences.

Fig. 7
figure 7

Expression levels of genes involved in oxidative stress and protein synthesis in gills after copper ion stress. The vertical axis represents the relative expression of these genes, the horizontal axis is 0 h, 6 h, 36 h, and 72 h of copper treatment

Fig. 8
figure 8

Expression levels of genes involved in oxidative stress and protein synthesis in liver after copper ion stress. The vertical axis represents the relative expression of these genes, the horizontal axis is 0 h, 6 h, 36 h, and 72 h of copper treatment

Discussion

Different tissues in fish respond differently to copper exposure

Among the major pollutants in the environment, metal ions pose a serious threat to aquatic organisms due to their unique bioaccumulation and toxicity. When these harmful substances are absorbed by fish, they damage many organs. Cadmium (Cd), for example, is a toxic heavy metal that adversely affects the survival and growth of aquatic organisms [22, 23]. Thus, it is a model oxidative stressor to study the mechanism of chemical toxicity. A recent study found that this toxic substance can lead to a liver malfunction in zebrafish as indicated by RNA-seq data [7]. Also, this heavy metal could elicit oxidative stress during zebrafish larval development [24]. However, in a previous study using yellow perch (Perca flavescens) as an experimental model, there was no significant difference in the accumulation of renal metals in Cd-exposed fish compared to controls [10].

Therefore, in this study, we investigated the response of different tissues to exposure to copper. To this end, RNA-seq data were analyzed in three different tissues: gills, liver and kidney. After 72 h of copper ion treatment, there were acute changes in the numbers of DEGs in the gills and liver, but not the kidney. This result indicates that copper exposure caused different responses, since the gene expression patterns changed substantially in the gills and liver, while the effect on the kidney was relatively weak, which is consistent with a previous study (Supplementary Table S2) [10].

Pathway of metal-induced response in gills and liver

Fish gills are particularly vulnerable to many water-borne pollutants, including heavy metals and chemically toxic substances. These pollutants often cause considerable ultrastructural damage, including epithelial cell proliferation, chloride cell necrosis, lamellar cell swelling and fusion, and mucous cell hyperplasia [25, 26]. In addition, the liver tissue of fish is also vulnerable to environmental stress, such as chemical pollutants and heavy metal ions [27, 28]. Therefore, we explored the response of gills and liver of fish to copper exposure.

From the KEGG enrichment, the differentially expressed genes (DEGs) in both gills and liver were enriched in “Ribosome” (Fig. 6). The ribosome is a complex ribonucleoprotein-based molecular apparatus that that is responsible for the synthesis of proteins within cells. Ribosome proteins and RNA of ribosomes are affected by reactive oxygen radicals, which may impair the homeostasis and stability of cells and even lead to the complete loss of ribosomal function [29]. The binding of divalent metal ions can cause chemical changes in biomacromolecules, either through metal ions directly reacting with the bound molecules or by promoting the local generation of ROS [29]. For example, combining Pb2+ with purified E. coli ribosomes can results in site-specific cleavages in rRNA [30]. ROS can perform essential biochemical and signaling functions or cause damage to cell components, including ribosomes [31]. In order to minimize damage caused by ROS, cells have developed various defense enzymes that modulate the stability of the system in an optimal or tolerable range. However, once the load of reactive oxygen species exceeds the capacity of the cellular antioxidant system, the organelles in the cells may be damaged and many cell functions may suffer.

Oxidized transition metals such as copper ions play a special role in the stability and stability of ribosomes. Copper is thought to be a transitional metal that promotes oxidative damage of nucleic acids and proteins [32]. In this study, we found that DEGs in copper treated fish gills and liver compared with control were enriched in “Ribosome”. In addition, we evaluated expressional levels of genes involved in oxidative stress response by qPCR. Our results indicate that the expression of Hmox1 (Cluster-47,180.157989), Nqo1 (Cluster-47,180.117564), Gpx1 (Cluster-47,774.0), Cu/Zn-SOD (Cluster-47,180.196300), Prdx1 (Cluster-47,180.233418), Gadd45b1 (Cluster-47,180.146985), Hspa9 (Cluster-47,180.10385), Hspa14 (Cluster-47,180.185583), and Mn-SOD (Cluster-47,180.180330) was greatly up-regulated in the gills, while the expression of Nqo1 (Cluster-47,180.117564), Cu/Zn-SOD (Cluster-47,180.117564), Prdx1 (Cluster-47,180.233418), Gadd45b1 (Cluster-47,180.233418), Hspa14 (Cluster-47,180.185583), and Mn-SOD (Cluster-47,180.180330) was greatly up-regulated in the liver (Fig. 7). In previous study, expressional levels of Hmox1, Nqo1, Gpx1, Cu/Zn-SOD, Prdx1, Gadd45b1, Hspa9, Hspa14, and Mn-SOD were used as indicators to detect oxidative stress in zebrafish larvae [26]. Upward regulation of these genes in the gills and liver suggested that copper treatment induced oxidative stress, which may lead to ribosome damage [29,30,31]. Therefore, our results support the notion that divalent metal ions can lead to damage of intracellular organelles [31], since gene expression patterns related to ribosome function were affected by copper exposure.

The biological function of ribosomes is to organize protein synthesis in cells. Protein synthesis is a basic, energy-consuming process for all living things. Cells often respond to adverse stimuli through coordinated changes in gene expression [33]. It is noteworthy that some genes involved in regulating cell growth and proliferation, differentiation, and apoptosis use other strategies for translation. Exposure to certain metals, such as hydrargyrum, copper, and cadmium, may lead to the down-regulation of global protein synthesis in some marine organisms [34, 35]. The effect of oxidative stress on protein synthesis is a basic problem in biomedical research. Many studies have shown that oxidative stress inhibits protein synthesis in cells [36, 37]. But little is known about the mechanism of the process. Numerous studies suggested that the regulatory mechanism of protein synthesis involves transcriptional initiation, elongation, and termination [35,36,37,38]. In the translational responses to oxidative stress, transcription initiation was inhibited by the initiator, eIF2a, and 4E-BPs phosphorylation [35]. Furthermore, elongation factors were also sensitive to oxidative stress. In Saccharomyces cerevisiae, after exposure to H2O2 and cadmium, the elongation factor eEF1Bα, eEF2 and eEF3 were reduced immediately [38]. Similarly, our results relate to the extended stage of translation regulation: Eef1b2 (Cluster-47,453.0) expression level was decreased in the liver under oxidative stress. Moreover, Hmox1 (Cluster-47,180.157989), Gpx1 (Cluster-47,774.0), and Hspa8b (Cluster-47,180.176206) were also decreased in the liver after copper ion stress. Therefore, we speculate that this strategy not only exists in the selective translation of proteins, but is also present in the specific translation of functional proteins in tissues and cells.

Conclusions

RNA-seq from gills, liver, and kidney of piebald naked carp to investigate the mechanism of copper toxicology. Our results show that the response of different tissues to copper exposure varied substantially. The gene expression patterns changed substantially in the gills and liver, indicating that these tissues are vulnerable to waterborne pollutants. In contrast, the kidney is relatively less sensitive. Furthermore, our results support the notion that divalent metal ions can affect the expression of genes involved in ribosome, and thus inflict damage on cellular components, since genes involved in the ribosome displayed different expression patterns. Also, the translation stage of proteins changes on a large scale during oxidative stress. Copper exposure can lead to a decrease in the synthesis of some proteins by down-regulating the expression level of elongation factors, and increase the translation of other proteins to repair damage and adapt to a stressful environment.

Methods

Fish samples and treatment

About 3 years old adult male piebald naked carp (Gymnocypris eckloni), weight 100 ± 2.5 g, Length 19 ± 3.4 cm, were obtained from the Fisheries Environmental Monitoring Station, Qinghai Province, China, and brought back to the laboratory in Qinghai University. During the domestication period-two weeks, fish are raised in fresh water at 15 ± 0.5 °C with a composition almost identical to that of rivers in wild habitats. At the end of this period, fish were placed in several 40 L glass aquariums with circulation systems and kept at 15 ± 0.5 °C during a 14 h:10 h light/dark cycle.

After acclimation, naked carp were exposed to metal ion stressors for different time durations. Exposure levels of oxygen and diet are the same as during adaptation. Twelve fish were randomly sampled and divided into four groups, the control group, heavy metal treated groups for 6, 36, and 72 h, and each group contained 3 fish.

CuSO4·5H2O was purchased from Sigma-Aldrich Chemical Co. (St Louis, MO, USA). The National Drinking Water Quality Standards in the People’s Republic of China (GB5749–2006) (copper ion concentration is not greater than 0.01 mg/L) was adopted to detect a series of biomarkers, which can more realistically reflect the response of aquatic organisms to heavy metal pollutants, and have more practical significance for the early warning of heavy metal pollution in water quality and environment. Besides, the metal ion concentration was chosen based on published works [10, 11]. Test solutions were completely changed every 24 h.

RNA extraction and transcriptome sequencing

Animal euthanasia was under the approved QHU Laboratory Animal Care and Use Committee protocols. In short, each fish was placed in a separate 2 L container containing 1.5 L of unbuffered 0.1% (100 mg/L) of clove oil until they lost their sense of direction and the gill cover stopped moving. Generally, adult fish of 19 ± 3.4 cm usually need 2–5 min. After 0.01 mg/L Cu2+ stressed for 6, 36, and 72 h, the gills, liver, and kidney were dissected and immediately frozen in the liquid nitrogen, stored at − 80 °C until RNA extraction.

Total RNA was extracted using RNAiso Plus (TaKaRa) as described by the manufacturer, and then perform DNase I processing. Degradation and contamination of RNA were detected by 1% agarose gel. In order to check the purity of RNA, the A260/A280 ratios of all RNA samples were measured by using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). The RNA integrity was assessed with the RNA integrity number (RIN) values by using RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Thirty-six separate Illumina sequencing libraries were established with 1.0 μg RNA per sample.

Sequencing libraries for each sample (3 replicates per sample) were prepared with NEBNext® Ultra™ RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s recommendations. An index code was added to assign a sequence to each sample. The sequencing was performed on Illumina HiSeq 2500 instrument at Novogene Bioinformatics Institute, Beijing, China.

Quality control and de novo transcriptome assembly

The original data in FASTQ format was filtered and processed with a self-written Perl script. In this step, adapter sequences, read containing poly-N, and low-quality reads in raw data were removed to generate clean data. In addition, levels of Q20, Q30, GC content and sequence duplication of clean data were calculated to monitor data quality. All downstream analyses were based on high-quality clean data.

The high-quality short reading of all samples was performed using Trinity software to produce a single de novo transcriptome assembly a single end-to-end transcriptome component [39], in which min_kmer_cov set to 2 and all other parameters set by default. First, all overlapping k-mers were extracted from the RNA-seq reads to generate transcript contigs using a greedy extension based on (k-1)-mer overlaps. In this step, full-length transcripts for a dominant isoform were generated, but only the unique portions of alternatively spliced transcripts were retained for subsequent analysis. Next, the related contigs were clustered using the raw reads to group transcripts based on shared read support and paired read links, when available. These steps were repeated to connect the contigs until the sequences could no longer be extended. Finally, all plausible transcript sequences, resolving alternatively spliced isoforms and transcripts derived from paralogous genes, were generated by analyzing the paths taken by reads and read pairings in the de Bruijn graph [40]. These sequences were considered to be unigenes.

Gene functional annotation

Candidate coding regions within transcript sequences were identified using TransDecoder (http://transdecoder.github.io/). The functional annotation of the unigenes was conducted using the BLASTX algorithm of the DIAMOND program against 3 reference databases: the NCBI non-redundant (Nr) protein database, the Clusters of Orthologous Groups of proteins (KOG/COG) database, and the Swiss-Prot database [41]. The functional annotation was also performed using four other databases, including the NCBI non-redundant (Nr) nucleotide database by NCBI blast 2.2.28+, the Protein family (Pfam) database by hmmscan in HMMER 3.0 package, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database by KAAS [42], and the Gene Ontology (GO) database by Blast2GO v2.5 [43].

Differential expression analysis and functional enrichment

The level of gene expression was measured using the criteria of reads per kb per million mapped reads (RPKM) [44]. Differentially expressed genes (DEGs) were detected using the DESeq software [45]. The P values were adjusted using the Benjamini and Hochberg multiple testing for controlling the false discovery rate. The gene with an adjusted P-value < 0.05 was assigned as differentially expressed. GO enrichment analysis was performed by the GOseq R package [46], and the gene length deviation was corrected. KEGG enrichment was performed by the KOBAS software [47]. Significance analysis was performed by Fisher’s Exact Test.

The expression level of genes involved in oxidative stress and protein synthesis

Seven putative Nrf2-dependent antioxidant genes (hmox1, sod1, sod2, gstp1, prdx1, gpx1, nqo1), three stress response gene (Hspa9, Hspa8, Hspa14), an inducible DNA damage repair gene (gadd45bb), two translation initiation factors (eIF2α, Eif4ebp1), and two elongation factors (EeF1b2, eEF2) were identified from the gene prediction results [24, 35]. Expression levels of these genes were measured using qPCR and RNA-seq data [48].

qPCR analysis

To validate the RNA-seq results, qPCR was performed by the following method. The product cDNA was diluted into 50 ng/μL, and 1 μL was used for template. The qRT-PCR reaction mixtures (20 μL) also contained 0.8 μL of each gene specific primer, 10 μL abm®EvaGreen qPCR MasterMix-no dye (abm, Canada) and 7.4 μL RNase-free water. The thermal cycling conditions are as follows: pre-denaturation at 95 °C for 10 min, 40 cycles of denaturation at 95 °C for 30 s, annealing at 60 °C for 1 min, followed by extension at 72 °C for 30 s in Roche LightCycler 480 Real-Time PCR system. The 2−ΔΔCt method was used to calculate the relative gene expression level. The β-actin of Gymnocypris eckloni was chosen as a reference gene to determine the relative level of the target genes. All qPCR experiments were performed in triplicate (Supplementary file 3) and the gene-specific primers were listed in supplementary table (Table S1).

Availability of data and materials

The datasets supporting the results of this article are available in the NCBI BioProject (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA666236) and NCBI Sequence Read Archive (SRA) repository (https://www.ncbi.nlm.nih.gov/sra/) under accession PRJNA666236.

Abbreviations

SOD:

Superoxide dismutase

G6PDH:

Glucose-6-phosphate dehydrogenase

Cd:

Cadmium

ROS:

Reactive oxygen species

QHU:

Qinghai University

RIN:

RNA integrity number

References

  1. Burke J, Handy RD. Sodium-sensitive and -insensitive copper accumulation by isolated intestinal cells of rainbow trout Oncorhynchus mykiss. J Exp Biol. 2005;208(2):391–407. https://doi.org/10.1242/jeb.01379.

    Article  CAS  PubMed  Google Scholar 

  2. Dautremepuits C, Paris-Palacios S, Betoulle S, Vernet G. Modulation in hepatic and head kidney parameters of carp (Cyprinus carpio L.) induced by copper and chitosan. Comp Biochem Physiol C Toxicol Pharmacol. 2004;137:325–33.

    Article  PubMed  Google Scholar 

  3. Sanchez W, Palluel O, Meunier L, Coquery M, Porcher JM, Aït-Aïssa S. Copper-induced oxidative stress in three-spined stickleback: relationship with hepatic metal levels. Environ Toxicol Pharmacol. 2005;19(1):177–83. https://doi.org/10.1016/j.etap.2004.07.003.

    Article  CAS  PubMed  Google Scholar 

  4. Carvalho CS, Fernandes MN. Effect of copper on liver key enzymes of anaerobic glucose metabolism from freshwater tropical fish Prochilodus lineatus. Comp Biochem Physiol A Mol Integr Physiol. 2008;151(3):437–42. https://doi.org/10.1016/j.cbpa.2007.04.016.

    Article  CAS  PubMed  Google Scholar 

  5. Lushchak VI. Environmentally induced oxidative stress in aquatic animals. Aquat Toxicol. 2011;101(1):13–30. https://doi.org/10.1016/j.aquatox.2010.10.006.

    Article  CAS  PubMed  Google Scholar 

  6. Lushchak VI. Contaminant-induced oxidative stress in fish: a mechanistic approach. Fish Physiol Biochem. 2016;42(2):711–47. https://doi.org/10.1007/s10695-015-0171-5.

    Article  CAS  PubMed  Google Scholar 

  7. Gao M, Lv M, Liu Y, Song Z. Transcriptome analysis of the effects of cd and nanomaterial-loaded cd on the liver in zebrafish. Ecotoxicol Environ Saf. 2018;164:530–9. https://doi.org/10.1016/j.ecoenv.2018.08.068.

    Article  CAS  PubMed  Google Scholar 

  8. Beaumont MW, Butler PJ, Taylor EW. Exposure of brown trout Salmo trutta to a sublethal concentration of copper in soft acidic water: effects upon gas exchange and ammonia accumulation. J Exp Biol. 2003;206(1):153–62. https://doi.org/10.1242/jeb.00060.

    Article  CAS  PubMed  Google Scholar 

  9. Husak VV, Mosiichuk NM, Kubrak OI, Matviishyn TM, Storey JM, Storey KB, et al. Acute exposure to copper induces variable intensity of oxidative stress in goldfish tissues. Fish Physiol Biochem. 2018;44(3):841–52. https://doi.org/10.1007/s10695-018-0473-5.

    Article  CAS  PubMed  Google Scholar 

  10. Grasset J, Ollivier É, Bougas B, Yannic G, Campbell PG, Bernatchez L, et al. Combined effects of temperature changes and metal contamination at different levels of biological organization in yellow perch. Aquat Toxicol. 2016;177:324–32. https://doi.org/10.1016/j.aquatox.2016.06.008.

    Article  CAS  PubMed  Google Scholar 

  11. Wang ZJ, Liu XH, Jin L, Pu DY, Huang J, Zhang YG. Transcriptome profiling analysis of rare minnow (Gobiocypris rarus) gills after waterborne cadmium exposure. Comp Biochem Physiol Part D Genomics Proteomics. 2016;19:120–8.

    Article  CAS  PubMed  Google Scholar 

  12. Yang L, Fang J, Peng X, Cui H, He M, Zuo Z, et al. Study on the morphology, histology and enzymatic activity of the digestive tract of Gymnocypris eckloni Herzenstein. Fish Physiol Biochem. 2017;43(4):1175–85. https://doi.org/10.1007/s10695-017-0363-2.

    Article  CAS  PubMed  Google Scholar 

  13. Qi D, Chao Y, Wu R, Xia M, Chen Q, Zheng Z. Transcriptome analysis provides insights into the adaptive responses to hypoxia of a schizothoracine fish (Gymnocypris eckloni). Front Physiol. 2018;9:1326. https://doi.org/10.3389/fphys.2018.01326.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Neal C, Robson AJ. A summary of river water quality data collected within the Land-Ocean interaction study: core data for eastern UK rivers draining to the North Sea. Sci Total Environ. 2000;251-252:585–665. https://doi.org/10.1016/S0048-9697(00)00397-1.

    Article  CAS  PubMed  Google Scholar 

  15. Mansour SA, Sidky MM. Ecotoxicological studies. 3. Heavy metals contaminating water and fish from Fayoum governorate, Egypt. Food Chem. 2002;78(1):15–22. https://doi.org/10.1016/S0308-8146(01)00197-2.

    Article  CAS  Google Scholar 

  16. Morillo J, Usero J, Gracia I. Biomonitoring of trace metals in a mine-polluted estuarine system (Spain). Chemosphere. 2005;58(10):1421–30. https://doi.org/10.1016/j.chemosphere.2004.09.093.

    Article  CAS  PubMed  Google Scholar 

  17. Bentley BP, Haas BJ, Tedeschi JN, Berry O. Loggerhead Sea turtle embryos (Caretta caretta) regulate expression of stress response and developmental genes when exposed to a biologically realistic heat stress. Mol Ecol. 2017;26(11):2978–92. https://doi.org/10.1111/mec.14087.

    Article  CAS  PubMed  Google Scholar 

  18. Faherty SL, Villanueva-Cañas JL, Blanco MB, Albà MM, Yoder AD. Transcriptomics in the wild: hibernation physiology in free-ranging dwarf lemurs. Mol Ecol. 2018;27(3):709–22. https://doi.org/10.1111/mec.14483.

    Article  CAS  PubMed  Google Scholar 

  19. Nespolo RF, Gaitan-Espitia JD, Quintero-Galvis JF, Fernandez FV, Silva AX, Molina C, et al. A functional transcriptomic analysis in the relict marsupial Dromiciops gliroides reveals adaptive regulation of protective functions during hibernation. Mol Ecol. 2018;27(22):4489–500. https://doi.org/10.1111/mec.14876.

    Article  CAS  PubMed  Google Scholar 

  20. Tong GX, Xu W, Zhang YQ, Zhang QY, Yin JS, Kuang YY. De novo assembly and characterization of the Hucho taimen transcriptome. Ecol Evol. 2017;8:1271–85.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Zhao Y, Liu X, Sato H, Zhang Q, Li A, Zhang J. RNA-seq analysis of local tissue of Carassius auratus gibelio with pharyngeal myxobolosis: insights into the pharyngeal mucosal immune response in a fish-parasite dialogue. Fish Shellfish Immunol. 2019;94:99–112. https://doi.org/10.1016/j.fsi.2019.08.076.

    Article  CAS  PubMed  Google Scholar 

  22. Chakraborty PK, Scharner B, Jurasovic J, Messner B, Bernhard D, Thévenod F. Chronic cadmium exposure induces transcriptional activation of the Wnt pathway and upregulation of epithelial-to-mesenchymal transition markers in mouse kidney. Toxicol Lett. 2010;198(1):69–76. https://doi.org/10.1016/j.toxlet.2010.05.007.

    Article  CAS  PubMed  Google Scholar 

  23. Prozialeck WC, Edwards JR. Mechanisms of cadmium-induced proximal tubule injury: new insights with implications for biomonitoring and therapeutic interventions. J Pharmacol Exp Ther. 2012;343(1):2–12. https://doi.org/10.1124/jpet.110.166769.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Mills MG, Gallagher EP. A targeted gene expression platform allows for rapid analysis of chemical-induced antioxidant mRNA expression in zebrafish larvae. PLoS One. 2017;12(2):e0171025. https://doi.org/10.1371/journal.pone.0171025.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Kirk RS, Lewis JW. An evaluation of pollutant induced changes in the gills of rainbow trout using scanning electron microscopy. Environ Technol. 1993;14:577–85.

    Article  CAS  Google Scholar 

  26. Wilson RW, Taylor EW. Transbranchial ammonia gradients and acid-base responses to high external ammonia concentration in rainbow trout (Oncorhynchus mykiss) acclimated to different salinities. J Exp Biol. 1993;166:95–112.

    Article  Google Scholar 

  27. Roda JFB, Lauer MM, Risso WE. Bueno dos Reis Martinez C. microplastics and copper effects on the neotropical teleost Prochilodus lineatus: is there any interaction? Comparative biochemistry and physiology. Part a, molecular and integrative. Physiology. 2020;242:110659.

    CAS  Google Scholar 

  28. Gopi N, Vijayakumar S, Thaya R, Govindarajan M, Alharbi NS, Kadaikunnan S, et al. Chronic exposure of Oreochromis niloticus to sub-lethal copper concentrations: effects on growth, antioxidant, non-enzymatic antioxidant, oxidative stress and non-specific immune responses. J Trace Elem Med Biol. 2019;55:170–9. https://doi.org/10.1016/j.jtemb.2019.06.011.

    Article  CAS  PubMed  Google Scholar 

  29. Shcherbik N, Pestov DG. The Impact of Oxidative Stress on Ribosomes: From Injury to Regulation. Cells. 2019;8:E1379.

    Article  PubMed  Google Scholar 

  30. Winter D, Polacek N, Halama I, Streicher B, Barta A. Lead-catalysed specific cleavage of ribosomal RNAs. Nucleic Acids Res. 1997;25(9):1817–24. https://doi.org/10.1093/nar/25.9.1817.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Schieber M, Chandel NS. ROS function in redox signaling and oxidative stress. Curr Biol. 2014;24(10):R453–62. https://doi.org/10.1016/j.cub.2014.03.034.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Dusek P, Roos PM, Litwin T, Schneider SA, Flaten TP, Aaseth J. The neurotoxicity of iron, copper and manganese in Parkinson's and Wilson's diseases. J Trace Elem Med Biol. 2015;31:193–203. https://doi.org/10.1016/j.jtemb.2014.05.007.

    Article  CAS  PubMed  Google Scholar 

  33. Holcik M, Sonenberg N. Translational control in stress and apoptosis. Nat Rev Mol Cell Biol. 2005;6(4):318–27. https://doi.org/10.1038/nrm1618.

    Article  CAS  PubMed  Google Scholar 

  34. Kalpaxis DL, Theos C, Xaplanteri MA, Dinos GP, Catsiki AV, Leotsinidis M. Biomonitoring of gulf of Patras, N. Peloponnesus, Greece. Application of a biomarker suite including evaluation of translation efficiency in Mytilus galloprovincialis cells. Environ Res. 2004;94(2):211–20. https://doi.org/10.1016/S0013-9351(03)00048-3.

    Article  CAS  PubMed  Google Scholar 

  35. Pytharopoulou S, Grintzalis K, Sazakli E, Leotsinidis M, Georgiou CD, Kalpaxis DL. Translational responses and oxidative stress of mussels experimentally exposed to hg, cu and cd: one pattern does not fit at all. Aquat Toxicol. 2011;105(1-2):157–65. https://doi.org/10.1016/j.aquatox.2011.06.007.

    Article  CAS  PubMed  Google Scholar 

  36. Shenton D, Smirnova JB, Selley JN, Carroll K, Hubbard SJ, Pavitt GD, et al. Global translational responses to oxidative stress impact upon multiple levels of protein synthesis. J Biol Chem. 2006;281(39):29011–21. https://doi.org/10.1074/jbc.M601545200.

    Article  CAS  PubMed  Google Scholar 

  37. Topf U, Suppanz I, Samluk L, Wrobel L, Böser A, Sakowska P, et al. Quantitative proteomics identifies redox switches for global translation modulation by mitochondrially produced reactive oxygen species. Nat Commun. 2018;9(1):324. https://doi.org/10.1038/s41467-017-02694-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Vido K, Spector D, Lagniel G, Lopez S, Toledano MB, Labarre J. A proteome analysis of the cadmium response in Saccharomyces cerevisiae. J Biol Chem. 2001;276(11):8469–74. https://doi.org/10.1074/jbc.M008708200.

    Article  CAS  PubMed  Google Scholar 

  39. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52. https://doi.org/10.1038/nbt.1883.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512. https://doi.org/10.1038/nprot.2013.084.

    Article  CAS  PubMed  Google Scholar 

  41. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12(1):59–60. https://doi.org/10.1038/nmeth.3176.

    Article  CAS  PubMed  Google Scholar 

  42. Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35:182–5.

    Article  Google Scholar 

  43. Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36(10):3420–35. https://doi.org/10.1093/nar/gkn176.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5. https://doi.org/10.1038/nbt.1621.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106. https://doi.org/10.1186/gb-2010-11-10-r106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14. https://doi.org/10.1186/gb-2010-11-2-r14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93. https://doi.org/10.1093/bioinformatics/bti430.

    Article  CAS  PubMed  Google Scholar 

  48. Guðbrandsson J, Franzdóttir SR, Kristjánsson BK, Ahi EP, Maier VH, Kapralova KH, et al. Differential gene expression during early development in recently evolved and sympatric Arctic charr morphs. PeerJ. 2018;6:e4345. https://doi.org/10.7717/peerj.4345.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We are grateful to the reviewers for their constructive comments and suggestions. We also express our gratitude to Mr. Zhenji Wang and Mr. Guojie Wang from Fisheries Environmental Monitoring Station (Xining, China) for their support to conduct the collection and transportation of specimens.

Funding

This work was financially supported by grants from the National Natural Science Foundation of China (NSFC, 31460700) to CZL, and the Major Science and Technology Projects of Qinghai Province (2019-NK-A2) to WJJ. These funding bodies had no role in the design of the study, collection, analysis, and interpretation of data, or in writing the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

WJJ and CZL conceived and designed the experimental plan. ZXL, FXR, SH, KFH, JJL, QSH, GJW, ZJW, SLJ, KML, and BTZ participated in collecting samples and tissues, and performing the experiments. WJJ analyzed the data and interpreted the results. WJJ and CZL drafted and revised this manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Changzhong Li.

Ethics declarations

Ethics approval and consent to participate

Rearing and all experiments were conducted with the approval of the Laboratory Animal Care and Use Committee at Qinghai University.

Consent for publication

Not applicable.

Competing interests

The authors declare they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

G3vsG0.DEG enriched KEGG pathway API.

Additional file 2.

L3vsL0.DEG enriched KEGG pathway API.

Additional file 3.

qPCR data.

Additional file 4: Table S1.

List of primers used in qRT-PCR

Additional file 5: Table S2.

Numbers of transcripts and genes in published papers.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Jin, W., Li, Z., Ran, F. et al. Transcriptome analysis provides insights into copper toxicology in piebald naked carp (Gymnocypris eckloni). BMC Genomics 22, 416 (2021). https://doi.org/10.1186/s12864-021-07673-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-021-07673-4

Keywords