1 INTRODUCTION

The scallop Chlamys farreri is an important commercial species that is cultivated widely by the rope-caging method in the northern coastal provinces of China. In recent years, anthropogenic activities, such as high summer temperatures and environmental pollution, have damaged scallop aquaculture (Nagashima et al., 2005). Many chemical contaminants and heavy metals are now considered to be harmful to aquaculture environments (Haynes and Johnson, 2000; Torres et al., 2008). Chlamys farreri are filter feeders that have been proved to be effective concentrators of heavy metals. Previous related studies on mollusk have focused mainly on the effects on antioxidant enzyme activity (Shan et al., 2011), lipid peroxidation (Lukyanova and Khotimchenko, 1995), testing bioaccumulation in tissues (Chelomin et al., 1995), DNA strand breakage (Slobodskova et al., 2012), and the discovery of genetic markers (Sato et al., 2005; Wang et al., 2012). However, these studies were limited due to the lack of genomic resources such as genome and transcriptomic sequences. Understanding the mechanism of detoxification in C. farreri is a high priority for the aquaculture industry. In order to achieve this goal, it is necessary to obtain a global understanding of the transcriptomic profile of C. farreri.

Cadmium (Cd) is a heavy metal that can be accumulated easily and it is one of the top 17 types of poison in contaminated food identified by the World Health Organization. Cd pollution can be accumulated through the food chain, and be enriched and amplified. According to an ecological assessment conducted from 2001 to 2005, Cd pollution is a serious problem in Bohai Bay, northeast China (Peng et al., 2009). Cd pollution has become a major problem that affects the quality and safety of aquatic products, especially scallops. Therefore, it is necessary to understand the toxicological mechanism in C. farreri.

In recent years, next-generation sequencing technologies have been used successfully for gene expression analysis, identification of differentially expressed genes, and discovery of novel transcripts (Garber et al., 2011) of genome reference-free species and non-model organisms. For example, RNA-seq has been applied in various invertebrates, such as Eriocheir sinensis (He et al., 2012; Ou et al., 2012), Bactrocera dorsalis (Hsu et al., 2012), and Pinctada martensii (Zhao et al., 2012) and other bivalves (Gavery and Roberts, 2012; Meng et al., 2014).

In the present study, we analyzed the digestive gland of C. farreri exposed to Cd based on Illumina sequencing and bioinformatics methods. The main objective of this study was to annotate functional genes related to Cd accumulation and identify potential molecular mechanisms involved in metabolic and immune pathways in C. farreri.

2 MATERIAL AND METHOD 2.1 Experimental scallop

Adult C. farreri scallops were collected from a commercial farm in Qingdao, China. Housing and care of the scallops and the collection of tissue samples were conducted in accordance with the International Guiding Principles for Biomedical Research Involving Animals. Only healthy scallops of a homogeneous size were used. The scallops were kept in running aerated sea water (salinity 30) at 20±1℃ that was pumped from Aoshanwei seaside of Qingdao and filtered prior to use. All of the seawater was changed daily and the scallops were fed with dried powder of Spirulina platensis (3 g/cm³ water per day).

2.2 Cadmium exposure and RNA preparation

The scallops were divided into two groups: one group was the control, and the other group was exposed to CdCl₂.2.5 H₂O. The final Cd²⁺ concentration reached 0.2 mg/L, which is 40 times of the standard formulated by 'Water Quality Standard for Fisheries of China' (Cd²⁺≤0.005 mg/L). No mortality was observed during the experimental period. After 15 days exposure, tissues of the digestive glands from five scallops in the control and five scallops in the Cd-treated groups were collected for RNA-seq. Two independent biological replicates (Cd-treated groups: XC1 and XC2; control groups: XD1 and XD2) were prepared for each group. Total RNA was isolated from tissues using TRIzol (TaKaRa, Dalian, China) according to the manufacturer's protocol. The quality and quantity of RNA was tested on an Agilent 2100 Bioanalyzer (Agilent Technologies) and a NanoDrop spectrophotometer (Thermo Scientific, USA), respectively. The total RNA was then treated with DNase Ⅰ (Ambion, USA) for 1 h at 37℃ to remove residual DNA from the total RNA.

2.3 cDNA library construction and RNA-seq

One paired-end (PE) cDNA library was generated from the total RNA of the digestive glands. Briefly, oligo (dT) beads were used to isolate poly (A) mRNA (mixture). The first-strand cDNA was synthesized using a random hexamer primer and reverse transcriptase, and the second-strand cDNA using DNA polymerase Ⅰ (TaKaRa) and RNase H (TaKaRa). The gel-purified cDNA product was modified using a TruSeqTM DNA Sample Prep Kit-Set A (Illumina) for sequencing, and then PCR amplified using a TruSeq Paired-End Cluster Kit (Illumina).

2.4 Analysis of Illumina sequencing results

The cDNA library was sequenced on the Illumina platform (GA Ⅱ). The raw data generated were assembled into non-redundant consensus unigenes sequences using two assemblers for very short reads: Velvet (version 0.7.62; http://www.ebi.ac.uk/zerbino/ velvet/) (Zerbino and Birney, 2008) and Oases (version 0.1.8; http://www.ebi.ac.uk/zerbino/oases/) (Schulz et al., 2012). Adapter sequences, low-quality reads, and short sequences ( < 100 bp) were removed using a custom Perl Program (Stajich et al., 2002). The GetORF program in the EMBOSS suite was used to predict open reading frames (ORFs) in the assembled unigenes (Rice et al., 2000). The ORFs were searched using BLASTx against the NCBI nonredundant (NR) and Swiss-Prot protein database with an E-value threshold of le-3. Functional annotation was performed based on the BLASTp results against the Swiss Prot and TrEMBL databases using the GoPipe program based on the gene2go software (http://gopipe.fishgenome.org), to annotate the gene ontology (GO) terms (http://www.geneontology.org). Pathway annotations were assigned using Blastall to search against the Eukaryotic Orthologous Groups (KOG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (Kanehisa et al., 2010).

2.5 Analysis of differentially expressed genes

For analysis of gene expression, the number of expression tags in one library was calculated and normalized to fragments per kilo base of transcript per million mapped reads (FPKM) (Mortazavi et al., 2008). Then, the expression abundance of each gene between the control and Cd-treated groups was counted based on the MARS model (MA-plot-based method with Random Sampling model) (Wang et al., 2010). In our study, we used a stringent FDR (false discovery rate) value of 0.001 and the absolute value of |log2 FoldChange|≤1 as the thresholds to judge a significant difference in gene expression between groups. All differentially expressed genes were mapped to terms in the KEGG database for pathway enrichment analysis.

2.6 Quantitative real-time PCR (qRT-PCR) verification

Six differentially expressed genes were selected for validation by qRT-PCR using a SYBR Premix Ex Taq kit (Takara, Japan) according to the manufacturer's instructions. The qRT-PCRs were performed in triplicate on the Thermo System. Total RNA was isolated from the digestive glands of C. farreri using TRIzol reagent. β-Actin was used as an endogenous control. The specific primers used for the qRT-PCRs are listed in Table S1. The relative gene expression was analyzed using the comparative threshold cycle (CT) method established by Livak and Schmittgen (2001).

3 RESULT AND DISCUSSION 3.1 Illumina paired-end sequencing and reads assembly

A cDNA pool was prepared from the digestive glands and sequenced to obtain a global overview of C. farreri transcriptome and related genes. We obtained approximately 227 673 transcripts. The raw reads produced have been deposited in the NCBI SRA database (accession number: SRR1563103). After adaptors and low-quality reads were removed, the filtered reads were de novo assembled with Trinity software. A total of 105 071 unigenes were obtained with an average size of 829 bp and a N50 length of 1 573 bp (Table 1). Among these unigenes, 61 828 (58.84%) were smaller than 500 bp, 19 393 (18.46%) were between 500 and 1 000 bp, 13 377 (12.73%) were between 1 000 and 2 000 bp, and 10 473 (9.97%) were longer than 2 000 bp. The length distribution of the unigenes is shown in Fig. 1.

3.2 Annotation of predicted proteins

All the predicted protein sequences (unigene sequences were translated by the Transeq program of the EMBOSS suite) were searched against the NR and Swiss-Prot protein sequence databases for gene annotation. A total of 24 576 unigenes matched NR sequences and were assigned biological functions. Then, the functions were classified based on assigned GO terms. A total of 24 493 unigenes were assigned at least one GO term under one of the three categories: biological process, cellular component, and molecular function. Under biological process, cellular process (10.41%) and metabolic process (8.68%) were the most abundant terms; under molecular function, binding (9.81%) and catalytic activity (7.05%) were the most represented terms; and under cellular component, cell (5.41%) and cell part (5.40%) were the most represented terms (Fig. 2). These results are similar to the findings of a previous study on Mizuhopecten yessoensis exposed to Cd (Meng et al., 2013). Additionally, a total of 12 028 unigenes were classified functionally into 26 KOG families, including biochemical metabolism, signal transduction mechanism, defense system, cellular structure, and gene expression (Fig. 3). In total, 7849 unigenes were assigned into different KEGG pathways (Fig. 4). The most represented pathways were signal transduction (1 349 members), carbohydrate metabolism (751 members), and transport and catabolism (563 members) (Table 2). A wellcategorized and annotated resource of the C. farreri transcriptome could become a valuable database for investigating specific bioprocesses and identifying important functional genes in this important commercial species.

3.3 Functional genes involved in environmental stress and metabolic, immune pathways

Based on the GO annotation, we identified 24 493 sequences potentially related to response to a Cd stimulus that could be valuable for future research. Invertebrates have developed unique systems of biological host defense (so-called innate immunity) because they lack immunoglobulins (Iwanaga and Lee, 2005). However, the exact molecular basis of the immune system remains poorly understood. For example, some of the unigenes were annotated as encoding members of the HSP (heat shock protein) family, including HSP22, HSP70, HSP90, and HSP transcription factor, which is known to be involved in the stress response in many species (Gupta et al., 2010). The abundant expression of HSP family transcripts is consistent with a previous report (Hou et al., 2011). In addition, other unigenes that were predicted to encode proteins involved in biological detoxification processes were detected, such as the CYP family, glutathione S-transferases, glutathione peroxidase, catalase, and superoxide dismutase.

Although information of immune-related genes on bivalve has increased, the available knowledge is still fragmentary (Moreira et al., 2012). The KEGG analysis identified unigenes associated with immuneresponse pathways such as Toll-like receptor signaling pathway, chemokine signaling pathway, NOD-like receptor signaling pathway, and complement and coagulation cascades in the immune system. Toll-like receptors are an ancient family of pattern recognition receptors that play major roles in activating the immune system and detecting non-self substances.We also detected unigenes encoding C3 and C1q, which are central components in the complement system that consists of a tightly regulated network of proteins that play crucial roles in host defense and inflammation response. In addition, MAPK (Kaminska, 2005), PI3 K (Wymann and Marone, 2005), and p38 (Kumar et al., 2003) were identified in corresponding pathways that are crucial in cellular responses to external stress signals. Thus, we identified a large set of immune-and metabolicrelated genes in the transcriptome data, which may play important roles in C. farreri responses and may help significantly improve our understanding of the molecular mechanisms involved.

3.4 Function annotation of differentially expressed genes

Genes that were differentially expressed by more than 2-fold between the control and Cd-treated groups, were mapped to terms in the KEGG database and compared with the transcriptome constructed from the pooled cDNA. A total of 3 800 differentially expressed sequences with KEGG annotations were observed between the Cd-treated and control groups (XC vs XD) (Fig. 5). Most of these genes were involved in biological processes such as ATP binding cassette (ABC) transporters, oxidative phosphorylation, ribosome biogenesis in eukaryotes, MAPK signaling pathway, glutathione metabolism, serine and threonine metabolism, and steroid hormone biosynthesis.

Following Cd exposure, several protein families including the ABC transporter, HSP, and CYP families were found to be up-regulated in the digestive gland of C. farreri. Previous work has indicated that ABCBand ABCC-type transporters form an environmenttissue barrier in bivalve gills and confer multixenobiotic resistance (Luckenbach and Epel, 2008). This protective mechanism is supposedly ubiquitous in aquatic organisms. Additionally, another research has suggested that the ABC multidrug transporters play an important role in detoxifying mechanism of various toxic chemicals in aquatic invertebrates (Kingtong et al., 2007). These findings may explain why we found that the ABC families were significantly up-regulated in the digestive gland. CYP3A, a member of the cytochrome P450 family, is known to be involved in the metabolism of a wide range of drugs in many organisms, and it is strongly induced by a variety of structurally unrelated xenobiotics (Hegelund and Celander, 2003). In mammals, erythromycin-N-demethylase (ERND) was found to be partially linked to CYP3A, and is considered as a useful biomarker to evaluate the adverse effect of environmental pollutants (Vaccaro et al., 2003). Previous studies have demonstrated that the expression of CYP3A in some fish species could be induced by certain xenobiotics (Hasselberg et al., 2004; Meucci and Arukwe, 2006). However, in this study, CYP3A was down-regulated in the digestive gland, probably because of the different species, tissues, or exposure times, dosages, and xenobiotics that were used. Metallothionein (MT) is a cysteinerich, low-molecular-weight protein that plays an important role in metal metabolism and detoxification (Olsson et al., 1989). It was reported that MT could be used as a bioindicator of metal contamination in aquatic environments through binding to particular metals. In the digestive gland of the edible mussel Mytilus edulis exposed to 200 μg/L Cd for 4 days, Geret and Cosson (2002) observed a significant increase of MT protein levels compared with a control. As reported by Zorita et al. (2007) an increase on the content of MT in the digestive gland of mediterranean mussels (Mytilus galloprovincialis) was induced by 200 μg/L Cd after 2 and 9 days of exposure. The levels of MT were also induced by 34 μg/L Cd for 5 days and 20 μg/L Cd for 63 days on zebra mussel (Dreissena polymorpha) (Lecoeur et al., 2004; Faria et al., 2009). In this study, the expression of MT was induced by 0.2 mg/L Cd significantly. Thus, the synthesis of MT has been found to increase significantly in response to Cd exposure as a protection mechanism against Cd toxicity. Metal tolerance protein (MTP) belongs to the cation diffusion facilitator protein family, which is widespread in bacteria, fungi, plants, and animals. MTP is a bivalent cation transporter that is necessary for efficient translocation of Zn, Cd, and other heavy metals (Ricachenevsky et al., 2013). In this study, exposure to Cd induced the expression of MTP in the C. farreri digestive gland. The transporter for nonheme iron at the luminal intestinal surface is divalent metal transporter 1 (DMT1), which is an H⁺-coupled and electrogenic membrane transporter that belongs to a large family of integral membrane proteins highly conserved throughout evolution. In our study, the expression of DMT1 decreased in the Cd-treated group, which may indicate that Cd inhibited the transportation of iron, as reported previously (Bannon et al., 2003).

3.5 Experimental validation of differentially expressed genes

Six differentially expressed genes were selected for validation by qRT-PCR. The results showed that five of the six primers produced bands of the expected sizes, and the five PCR products were further confirmed by Sanger sequencing (data not shown). However, all six genes showed the same expression change trends as the corresponding expressions estimated based on the transcriptomics data, which confirmed the effectiveness of our comparative transcriptomics analysis (Fig. 6).

3.6 Biological replicates

The results were replicated in a second experiment. The Pearson correlation was 0.852 for XC1 vs XC2, and 0.86 for XD1 vs XD2 (Fig. 7).

4 CONCLUSION

In this study, we generated a comprehensive transcriptome of adult C. farreri based on the Illumina sequencing platform. A total of 105 071 unigenes were assembled from the transcripts and 24 576 of them were assigned explicit annotation. A significant number of genes involved in the metabolic immune pathway, and response to a Cd stimulus were identified in the transcriptome. These findings substantially supplement existing genomic resources for C. farreri without reference genome information. In addition, differentially expressed genes were identified and functionally annotated based on database. The genetic information provided potential molecular markers in bivalves exposed to marine pollutants for functional studies of genes and better understanding of molecular mechanisms involved in metabolic and immune pathways.

5 ACKNOWLEDGMENT

The authors thank Dr. Yuefeng CAI (Ocean University of China) for useful advice and assistance.