Comparative transcriptome analysis of the gills and hepatopancreas from Macrobrachium rosenbergii exposed to the heavy metal Cadmium (Cd2+)

Heavy metal Cadmium (Cd2+) pollution has become a severe environmental problem for aquatic organisms. In crustaceans, gills (Gi) and hepatopancreas (Hp) play a vital role in the toxicology. However, in Macrobrachium rosenbergill, there are few researches about gill and hepatopancreases responding to Cd2+ stress at a molecular level. In this study, transcriptomic analysis was applied to characterize gene expression profiles of gills and hepatopancreas of M. rosenbergill after Cd2+ exposure for 0 h, 3 h and 3 d. Six cDNA libraries (Gi 0 h, Gi 3 h, Gi 3 d, Hp 0 h, Hp 3 h, and Hp 3 d) were constructed and a total of 66,676 transcripts and 48,991 unigenes were annotated. Furthermore, differentially expressed genes (DEGs) were isolated by comparing the Cd2+ treated time-point libraries (3 h and 3 d group) with the control library (0 h group). The results showed that most of the DEGs were down-regulated after Cd2+ exposure and the number of DEGs among gill groups were significantly higher than those among hepatopancreas groups. GO functional and KEGG pathway analysis suggested many key DEGs in response to the Cd2+ stress, such as metallothionein and Hemocyanin. Additionally, a total of six DEGs were randomly selected to further identify their expressional profile by qPCR. The results indicated that these DEGs were involved in the response to Cd2+. This comparative transcriptome provides valuable molecular information on the mechanisms of responding to Cd2+ stress in M. rosenbergii, which lays the foundation for further understanding of heavy metal stress.

results in various physiological damages to animal tissues and organs 13,14 . Furthermore, Cd 2+ causes impairment of reproductive activity and disrupts endocrine function in fish 15 . Cd 2+ induced cell apoptosis has been confirmed to be attributed to caspase-dependent and independent pathways of the mitochondria or endoplasmic reticulum (ER) 16,17 . Additionally, Cd 2+ , a non-essential and potentially toxic metal, can be accumulated in humans via food chain 18 , which may result in morphological deformities, physiological dysfunctions and even death 19 . Previous researches showed that Cd 2+ is known to accumulate in marine organisms and induced rapid genetic changes in many crustaceans, such as Sinopotamon henanense 19 and Eriocheir sinensis 20 . Hence, it is essential to focus on the potential response mechanism caused by Cd 2+ stress in crustaceans.
As an important respiratory organ, the gill (Gi) is involved in ion transport, acid-base balance and osmoregulation in crustaceans 21,22 . Due to the crustacean gills being exposed to the water in which they live, the gills play a vital role in the toxicological interactions, such as with heavy metals 23 . Furthermore, the hepatopancreas (Hp), a sensitive organ similar to the liver of higher organisms, is susceptible to be damaged by waterborne pollutants in crustaceans [24][25][26] . Therefore, gills and hepatopancreas are model organs for studying the response to heavy metal stress in crustaceans.
The giant freshwater prawn, Macrobrachium rosenbergii, is an important commercial prawn and widely cultured in China and other Pacific Rim countries 27 . As a freshwater cultured species, the prawn is susceptible to metal accumulation. Previous studies showed that structural changes of gills and hepatopancreas of M. rosenbergii could be caused by the Cu 2+ accumulation, and the degree of damage observed was related to the elevated waterborne copper concentration 24 . Additionally, transcriptomic analysis of gills of M. rosenbergii showed that 19,417 and 8,989 differentially expressed genes (DEGs) were identified at 3 h and 48 h after Cu 2+ exposure, respectively 28 , revealing that a large number of genes were involved in response to Cd 2+ stress. Further research showed that the accumulation of Cd 2+ also manifested histopathological changes in the gills and hepatopancreas of M. rosenbergii under Cd 2+ exposure, and Cd 2+ levels in tissues followed the order of: gills > hepatopancreas 29 . To date, however, limited researches were focused on the Cd 2+ -related stress response and regulatory gene in M. rosenbergii.
In this study, transcriptome sequencing of gills and hepatopancreas in M. rosenbergii was performed to analyze transcriptional responses under Cd 2+ pollution. Many vital genes in response to the Cd 2+ were identified. The study provided valuable and reliable data for aquaculture and environmental monitoring management, and elucidated the potential toxicological mechanism in M. rosenbergii.

Results
Transcriptome sequencing and functional gene annotation. Six

Identification of differentially expressed genes (DEGs). Pearson correlation analysis showed good
correlation among different replicates of the same sample, whereas significant differences were observed between the gill and hepatopancreas groups (Fig. 2). To identify genes displaying significant changes in expression level in the face of Cd 2+ stress, we analyzed the expression level of each unigene by TPM method and found many DEGs by comparing the Cd 2+ treated time-points libraries (3 h and 3 d group) with the control library (0 h group) (Supplementary Table S1). Compared to gill control group (Gi 0 h), a total of 6264 (2,010 upregulated and 4254 downregulated) and 5175 (2,186 upregulated and 2989 downregulated) DEGs were identified in the Gi 3 h group and Gi 3 d group (Fig. 3), respectively. Long duration of Cd 2+ exposure (Gi 3 d group) caused 4222 genes to be differentially expressed compared with short duration (Gi 3 h group) (Fig. 3). Furthermore, Venn analysis showed that 3375 genes were differentially expressed at both time-points, while 2889 DEGs were regulated just at Gi 3 h group and 1,800 genes were altered just at Gi 3 d group (Fig. 4) 5). Among the categories of biological process, cellular component, and molecular function, the top 2 enriched GO terms for each category were "cellular process and metabolic process", "membrane part and cell part", and "catalytic activity and binding", respectively. A total of 59 genes were annotated to "response to oxidative stress" GO terms (Supplementary Table S2)  GH) involved in metal transportation and stress response was in accordance with the results from RNA sequencing (Fig. 7). Obviously, the change range of expression level of many DEGs in gill treated groups is greater than those in hepatopancreas treated groups compared to control group. For example, the expression level of GH gene significantly increased by eight and seven times in gill at 3 h and 3 d group, respectively. In addition, the DEGs

Discussion
In crustaceans, the gill epithelium is generally regarded as a major organ of respiration and osmoregulation, and the first site to be exposed to environmental pollutants 30 . During waterborne exposure to heavy metals, gills act as a protective barrier between the internal and external environment 23 . Waterborne heavy metals are initially absorbed into epithelium cells of gill and transported into hemolymph, and finally infiltrated into internal organs 23 . Hepatopancreas is usually considered as a vital target organ for heavy metal toxicity and other environmental stresses in crustaceans and plays a major role in metal storage and in the detoxification process 31 . Additionally, crustaceans increase metabolic efficiency by promoting the digestive enzyme activities in hepatopancreas in response to heavy metal 32 . Therefore, the gill and hepatopancreas are considered as a good indicator of water quality, and a suitable model for studies of heavy metal pollution. To better understand the molecular mechanisms of Cd 2+ toxicity in M. rosenbergii., RNA-Seq was used to investigate gene expression differences of gill and hepatopancreas in response to Cd 2+ exposure (0 h, 3 h, and 3 d). Six cDNA libraries were constructed and a total of 48,991 unigenes were functionally annotated. GO term enrichment and KEGG pathway enrichment were performed to find important genes and pathways during Cd 2+ exposure in gill and hepatopancreas of M. rosenbergii. We analyzed DEGs by comparing the Cd 2+ treated timepoint libraries with the control library. The results showed that the number of down-regulated DEGs is larger than up-regulated DEGs (Fig. 3), indicating that gene expressions were mainly inhibited by Cd 2+ , which leads to impairments in M. rosenbergii. The results were similar with Sinopotamon henanense and Danio rerio under Cd 2+ stress 19,33 . The number of DEGs among Gi groups were significantly higher than those among Hp groups (Fig. 3), suggesting that the gill has a stronger stress response than hepatopancreas in short time. Additionally, the number of DEGs in Gi groups decreased with the increment of exposure time, while in Hp groups, the number of DEGs increased with the increment of exposure time (Fig. 3). The above results might be attributed to the reason that the gill acts as the entry site and transient store organ of the heavy metal for a short period of exposure time, and Cd 2+ is gradually transferred from the gills to hepatopancreas via the haemolymph with the prolongation of exposure time 34 .
Many genes related to oxidative stress were found in response to the Cd 2+ stress (Supplementary Table S2), and present various expression patterns, as identified by qPCR (Fig. 7). For example, the expression level of metallothionein (MT) was significantly increased at 3 h, then decreased at 3d, which may be related to the accumulation of Cd 2+ . Many studies have shown that MT is critical to heavy metal detoxification 35,36 in addition to storage of essential elements that are necessary for metalloenzymes 37,38 . Some studies have proven that the accumulation of heavy metal has significant time effects. For instance, in Oncorhynchus mykiss, Cu 2+ uptake increased during the 1-2 h under radiolabelled copper exposure, and after 2 h, Cu 2+ level significantly decreased in the gill 39 . A similar tendency was found in Acrossocheilus fasciatus 40 , in which the expression level of zinc-finger BED domaincontaining protein (Zbed) was significantly decreased after exposure to Cd 2+ , which is also consistent with what has been observed in Mytilus galloprovincialis exposed to Cu 2+41 . In contrast, hemocyanin-like protein, a crucial immune protein in arthropods [42][43][44] , had significantly increased expression after exposure to Cd 2+ . Heavy metals   45 . In addition, the expression level of heat shock proteins (Hsps), common stress-inducible proteins, has been known to increase under various stressors, such as oxidative stress, heavy metals, and viral infections [46][47][48] . For instance, Hsp70, Hsp40, and Hsp105, were significantly up-regulated in Eubalaena glacialis exposed to Cd 2+ . Interestingly, in M. rosenbergii, Hsp67B2 was consistently decreased in the hepatopancreas for three days under Cd 2+ exposure, suggesting that Hsp67B2 may be suppressed by Cd 2+ in this prawn. On the other hand, the expression level of IFRD1 was consistently increased in hepatopancreas under Cd 2+ stress, which was consistent with the high upregulation of this gene in hepatopancreas of M. rosenbergii after virus infection 49 . IFRD1 protein has been proven to be involved in the regulation of inflammatory responses 50 , indicating that the increased expression of IFRD1 is intended to cure inflammation caused by Cd 2+ . Nevertheless, further study is required to illustrate the regulatory mechanism of M. rosenbergii after exposure to Cd 2+ . The degree of histological damage of the gills and hepatopancreas under different concentrations and exposure days of Cd 2+ is worth exploring in future research. Additionally, the effects of Cd 2+ on the mitochondrion structure in the gill and on superoxide dismutase (SOD) activity still need to be investigated.

Conclusion
In conclusion, we successfully constructed comparative gill and hepatopancreas transcriptome datasets in Cd 2+ treated groups and control group of M. rosenbergii. Thereafter, 48,991 unigenes were functionally annotated and a series of DEGs were isolated after Cd 2+ exposure. Based on GO functional and KEGG pathway analyses, many DEGs that are potentially relevant to immune responses, antioxidant, and detoxification were identified.

Material and methods
Collection and maintenance of prawns. A total of nine female and nine male M. rosenbergii (23 ± 2.5 g) individuals used in this experiment were collected from Dinghe Aquatic Science and Technology Development Co. LTD (Jiangsu, China) and transported back to our laboratory. The prawns were maintained at 26 ± 2 °C in a 50-L aerated aquarium for three days before treatment. All animals were handled in accordance with guidelines established by the Animal Experiments Ethics Committee of Shanghai Ocean University for the care and use of laboratory animals.
Cadmium exposure experiment. Firstly, Cd 2+ solution (50 mg/L) was prepared by dissolving 102 mg of CdCl 2 ·2.5 H 2 O in 1 L deionized water. After temporary rearing, CdCl 2 solution was added to the culture water and mixed immediately, so as to expose all the prawns to the Cd 2+ (80 μg/L) based on the 96-h LC50 of Cd in M. rosenbergii 51 . Every day, the prawns were fed and the water was renewed by 50% to maintain water quality. Subsequently, the experimental prawns were anesthetized on ice and dissected. The gills (Gi) and hepatopancreas (Hp) were randomly sampled from six individuals (three males and three females) for each of the 3 time points: 0 h, 3 h and 3 d, after Cadmium exposure, and stored at − 80 °C immediately for the following RNA extraction.
Library construction and gene function annotation. Total

Validation of DEGs expression profiles using quantitative real-time RT-PCR (qPCR).
To validate the Illumina sequencing results, the six pooled RNA samples originally used for transcriptome sequencing were analyzed by qPCR. Six randomly selected genes: metallothionein (MT), hemocyanin-like protein (Hemo), interferon-related developmental regulator 1 (IFRD1), heat shock protein 67B2 (Hsp 67B2), zinc finger BED domain-containing protein 4 (Zbed4), and gamma-glutamyl hydrolase (GH), were amplified by specific primers (Table 4). QPCR mixture (20 μL) contained 10 μL of PCR Master with SYRB green, 1 μL Cd 2+ cDNA template (10 ng/ul), 0.25 μL of each primer (10 uM), and 8.5 μL H 2 O. The primers of β-actin were used as the internal control. The relative quantification of the six genes was calculated by the 2-△△ CT method 55 . Analysis of qPCR results was performed in GraphPad Prism 8. All data were presented as means ± SD.
Approval statement. All experimental protocols were approved by the Key Laboratory of Freshwater Aquatic Genetic Resources, Ministry of Agriculture, Shanghai Ocean University in this paper.

ARRIVE guidelines statement.
This study was carried out in compliance with the ARRIVE guidelines.