Transcriptome analysis of head kidney in grass carp and discovery of immune-related genes

Background Grass carp (Ctenopharyngodon idella) is one of the most economically important freshwater fish, but its production is often affected by diseases that cause serious economic losses. To date, no good breeding varieties have been obtained using the oriented cultivation technique. The ability to identify disease resistance genes in grass carp is important to cultivate disease-resistant varieties of grass carp. Results In this study, we constructed a non-normalized cDNA library of head kidney in grass carp, and, after clustering and assembly, we obtained 3,027 high-quality unigenes. Solexa sequencing was used to generate sequence tags from the transcriptomes of the head kidney in grass carp before and after grass carp reovirus (GCRV) infection. After processing, we obtained 22,144 tags that were differentially expressed by more than 2-fold between the uninfected and infected groups. 679 of the differentially expressed tags (3.1%) mapped to 483 of the unigenes (16.0%). The up-regulated and down-regulated unigenes were annotated using gene ontology terms; 16 were annotated as immune-related and 42 were of unknown function having no matches to any of the sequences in the databases that were used in the similarity searches. Semi-quantitative RT-PCR revealed four unknown unigenes that showed significant responses to the viral infection. Based on domain structure predictions, one of these sequences was found to encode a protein that contained two transmembrane domains and, therefore, may be a transmembrane protein. Here, we proposed that this novel unigene may encode a virus receptor or a protein that mediates the immune signalling pathway at the cell surface. Conclusion This study enriches the molecular basis data of grass carp and further confirms that, based on fish tissue-specific EST databases, transcriptome analysis is an effective route to discover novel functional genes.


Background
Grass carp (Ctenopharyngodon idella) is one of the most important freshwater fish, with fast growth, low cost of breeding, and delicious meat. It is widely distributed in China's major river systems. Grass carp is a farmed species that is easily affected by diseases induced by viruses and bacteria; this can cause tremendous economic losses. To date, no excellent breeding varieties have been obtained by the oriented cultivation technique. Because of the long breeding cycle (4-5 years), a hybrid breeding strategy is not feasible. Further, because of the lack of understanding of the genetic background of grass carp, no molecular breeding technology has been applied. The discovery of economically important trait-related genes and their functional study may help to establish a molecular breeding technology system in the fish.
ESTs (expressed sequence tags) are partial cDNA sequences obtained after sequencing the ends of random cDNA clones. ESTs were first used in 1991 as an effective new method to discover human genes. Using EST sequences, unknown genomes could be explored at a relatively low cost [1]. With the development of DNA sequencing technology, the cost of sequencing is becoming lower, and the application of large-scale EST sequencing combined with bioinformatics tools for analyzing data is being widely used in different species to find novel genes, for genome annotation, for the identification of gene structure and expression, and in the development of type I molecular markers [2]. In fish, large scale EST sequencing was used in channel catfish (Ictalurus punctatus) [3], common carp (Cyprinus carpio) [4], and zebrafish (Danio rerio) [5].
In recent years, high-throughput data analysis methods have gradually improved and the genomes of many kinds of fishes have been studied. The fishes that have been studied include zebrafish [6] and fugu [7], as model organisms, and the commercial fishes such as Atlantic salmon [8], sea bass [9,10], rainbow trout [11], Atlantic halibut [12], bluefin tuna [13], turbot [14,15], and Senegal sole fish [16]. In contrast, the molecular biology of grass carp is relatively unknown; currently, there are only 6,915 grass carp ESTs in NCBI's dbEST database. Most functional genomic research on economically important fish is focused mainly on the development of molecular markers, genetic map construction and gene interval mapping, and other basic data accumulation. Research into gene function and its application to breeding is still in the initial stages.
Head kidney is an important immune organ in teleost fish; its role is equivalent to mammalian bone marrow [17]. Head kidney contains a large number of T and B lymphocytes, macrophages and granulocytes that are the basis upon which specific and non-specific immunity is acquired.
In this study, we constructed a non-normalized cDNA library for the head kidney of grass carp and obtained 3,027 unigenes including 221 genes of unknown function. We compared the head kidney expression profiles of grass carp infected with grass carp reovirus (GCRV) with normal controls and obtained 22,144 differential expressed tags. Based on a comparison of the differential expressed tags and potential genes with unknown function in the cDNA library, and by identifying gene expression response to GCRV and predicting protein structure, we discovered a novel immune-related gene. This study provides a method for the discovery of novel genes, and reveals the function and the network regulation mechanism of immune-related genes. The results provide a theoretical foundation for molecular design breeding in grass carp.

RNA extraction and construction of the cDNA library
Total RNA was extracted from the head kidney of healthy adult grass carp using Trizol reagent (Invitrogen, Carlsbad, CA, USA). The mRNA was isolated using the Oligotex mRNA Kit (QIAGEN, Hilden, Germany). Full length cDNA was synthesized by the Creator TM SMART TM cDNA Library Construction Kit (Clontech, CA, USA) following the method described previously [18]. cDNA segments longer than 1 kb were isolated by electrophoresis, then ligated into pDNR-LIB vector (Clontech) and used to transform competent E. coli DH5α cells. After growing the colony for 12 hours on an LB plate containing chloramphenicol, the cDNA library was constructed by selecting mono-clones from the 96well plate. Ethical approval for the work was obtained from Expert Committee of Biomedical Ethics, Institute of Hydrobiology of the Chinese Academy of Sciences. The Reference number obtained was Y12202-1-303.
An optimal peak chart was obtained by processing the raw sequence data with basecalling. Next, FASTA format sequences (raw ESTs) were obtained by processing the optimal peak chart using the Phrap program [19] with the Q20 standard. We used crossmatch (Smith and Green, unpublished observations) to remove the pDNR-LIB vector sequences and after excluding EST sequences that were less than 100 bp long, we obtained a cleaned EST data set. Clustering of the cleaned ESTs was performed using UIcluster [20]. The UIcluster sequences were assembled using the Phrap program to build a unigene data set for the ESTs from the head kidney of grass carp.

BLAST searches, GO functional classification and KEGG pathway analysis
We used the NCBI BLAST server [21] to identify sequences that were similar to the sequences in the NCBI nucleotide sequence database (Nt), the protein sequence database (Nr) [22] and the Swissprot database [23] using BLASTN, and BLASTX [24]. Using the EST sequence with the highest homology as a guide, we set the threshold E-value to E < 1e-6.
We used the BLASTX search results from the Swissprot database and the Blast2GO tool [25] to assign GO functional classification to the unigene sequences. Blas-t2GO parameters were set as follows: E-Value-Hit-Filter < 1e-6; annotation cutOff = 55; other parameters remained at the default values.
KAAS [26] was used to assign the unigene ESTs to pathways based on KEGG Orthology (KO) [27]. Unigenes were mapped to the corresponding KEGG pathways using the comparison method of bi-directional best hit.

GCRV infection of grass carp and preparation of RNA sample
The GCRV-873 strain was provided by the Gaobo biotechnology company (Wuhan, China). One-year-old grass carp with an average weight of 180-210 g were intraperitoneally injected with 150-200 μL GCRV, a dosage of approximately 10 6 TCID 50 kg -1 body weight. The injected grass carp were raised in clean tanks at 28°C . Three infected grass carp with typical hemorrhage symptoms (infected group, n = 3) and three uninfected grass carp (healthy control group, n = 3) were selected at 5d after infection for further study. Total RNA was extracted from the head kidney of both groups using Trizol reagent. cDNA was obtained after reverse transcription and used for Solexa sequencing.
Three-month-old grass carp with an average weight of 30-60 g were intraperitoneally injected with 50-80 μL GCRV, a dosage of approximately 10 6 TCID 50 kg -1 body weight; fish in the control group were injected with same amount of saline. The grass carp were raised in clean tanks at 28°C. At 1d, 2d, 3d, 4d, 5d after infection ten GCRV-infected carp were selected for further study (n = 10). Ten uninfected fish were selected from the control group at 0d (n = 10). The whole fish was immediately used for RNA isolation. cDNA was obtained after reverse transcription and used for the detection of gene expression.

Solexa sequencing and expression profile analysis
The NlaIII and MmeI digestion method [28] was used to build a 21-bp cDNA tag library of the two groups (oneyear-old), the control group and the GCRV-infected group. The tags in the two libraries end with different Illumina adapter sequences. The raw sequencing read length was 35 bp. The Solexa sequencing was performed by the Beijing Genomics Institute (BGI, Shenzhen, China).
The raw sequence data was processed through basecalling, the adapter and low quality sequences were removed, and cleaned 21-bp tags were obtained. We converted the cleaned tag number into the standard (relative) number of transcripts per million (TPM), and calculated the logarithm of TPM for each of the cleaned tags from the control and GCRV-infected groups. We used a dual limit of P <0.01 and FPR (false positive rate) <0.01, to find cleaned tags with log2Ratio ≥ 1 or log2Ratio ≤ −1 [29]. The selected tags have differential expression levels of more than 2-fold in both groups. We then compared the differential expressed tags with the    unigenes from the cDNA library using SeqMap [30]; mismatch was set to 0, and sense and antisense strands were considered in the mapping.

Semi-quantitative RT-PCR and RACE cloning
Total RNA was used to synthesize the first strand cDNA. Upstream and downstream primers (Table 1) were designed based on the unigene sequences. β-actin (primers, β-actin-F and β-actin-R) was used as the internal reference. PCR and electrophoresis was used to detect the change of expression level. 3' and 5' RACE was performed using the BD SMART RACE cDNA Amplification Kit (Clontech) according to the manufacturer's instructions. Upstream and downstream primers used in the 3' and 5' RACE were designed based on the EST sequences (Table 1). Full length cDNA sequences of each gene were assembled using the 3' and 5' terminal sequences.

Head kidney cDNA library of grass carp
The storage capacity of the original library was 6 × 10 5 , in the form of the E. coli DH5α cells that were stored on the 532 96-well plates in a total of 51,072 clones. One hundred randomly selected clones were used for further study. The PCR test results showed that the size of inserts was between 1-3 kilobases, the library reorganization was 97.85% and the no-load rate was 2.15%.

EST sequence analysis
10,464 EST clones were sequenced, and 10,282 FASTA sequences (raw ESTs) with an average read length of 470 bp were obtained. After removing the vector and sequences less than 100 bp long, 7,918 cleaned ESTs (accession no. JK847435-JK855352) were obtained. After clustering and assembly, we obtained 3,027 unigene EST sequences, 802 (26.5%) of which were contigs and 2,225 (73.5%) of which were singletons; the library redundancy was 61.78%. Most genes in the library exhibited lowlevel expression, only a small number of genes exhibited high-abundance expression. The number of low expression unigenes, the singletons, was 2,225 (73.5%); the number of medium expression unigenes, those containing 2-5 ESTs was 641 (21.2%); and the number of high expression unigenes, those that contained six or more ESTs, was 161 (5.3%). Only 23 unigenes contained more than 20 ESTs. The average length of the unigenes was 431 bp and 77.33% of the unigenes were 300-500 bp long (Figure 1).

BLAST searches and GO functional classification
The 3,027 unigenes were used as queries in BLAST searches of the NCBI nucleotide and protein sequence databases and the Swissprot database. 2,713 unigenes (89.6%) matched sequences in the nucleotide sequence database, 2,162 unigenes (71.4%) matched sequences in protein sequence database and 1,845 unigenes (61.0%) matched sequences in the Swissprot database. In all, 2,806 unigenes (92.7%) matched sequences in at least one of the three databases; the remaining 221 unigenes (7.3%) were not found (E-value <1e-6) in any of the three databases and may be novel gene sequences. Using the gene ontology (GO) classification, we successfully assigned functional annotations to 1,323 of the unigene sequences. In the GO biological process ontology, three terms accounted for the largest proportion of unigenes, they were cellular process, metabolic process and biological regulation; in the GO molecular function ontology, the three most commonly occurring terms were binding, catalytic activity and structural molecule activity; and in the GO cellular component ontology, cell, cell part and organelle were the terms that occurred most frequently (Table 2). Of the 1,323 GO-annotated unigenes, 53 were immune system process-related genes (Table 3), 4 were response to virus, and 9 were response to bacterium process-related genes (Tables 4 and 5). Some unigenes were assigned multiple functions. Not all of the unigenes could be mapped to the lower level GO terms.

KEGG pathway analysis
A total of 989 of the 3,027 were assigned a KEGG ontology (KO) annotation; they were mapped to 201 KEGG pathways. Three most frequently occurring KEGG pathways were ribosome, oxidative phosphorylation, and proteasome. 68 unigenes mapped to immune-related pathways including leukocyte transendothelial migration, antigen processing and presentation, chemokine signalling pathway, and T cell receptor signalling pathway (Table 6). We found that 28 unigenes from head kidney in grass carp have been reported to be involved in the following pathways, Toll-like receptor signalling pathway, RIG-I-like receptor signalling pathway and the NOD-like receptor signalling pathway (Table 7).

Expression profiling analysis
By Solexa sequencing, we obtained 7,696,804 and 6,136,889 raw tags from the transcriptomes of head kidney tissue from grass carp before and after GCRV infection, respectively. After removing low quality sequences, adapter sequences and single copy sequence the cleaned tag numbers were 7,188,005 and 5,724,526, respectively.  The final numbers of non-redundant distinct tags were 152,826 and 105,653 before and after GCRV infection, respectively. All tags were submitted to SRA at NCBI under the accession no. SRA052520.2. Of the distinct tags, 22,144 were differentially expressed by more than 2-fold between the GCRV-infected and uninfected groups. These 22,144 differentially expressed tags mapped to 3,027 unigenes using SeqMap [30]. Of the differentially expressed tags, 679 (3.1%) mapped to 483 differentially expressed unigenes (16.0%); 145 of the unigenes were up-regulated genes, 307 were down-regulated genes. The remaining 31 unigenes mapped to tags that exhibited both up and down regulation, and so these unigenes were not included in the statistics. The up-and downregulated genes were mainly annotated with the GO terms, genetic information processing, metabolism, and cellular processes and 16 unigenes were annotated with the GO term immune-related (Table 8). We found 54 tags that mapped onto 42 of the 221 unknown unigenes. These are potentially infection related novel genes; 15 of them were up-regulated between the GCRV-infected and uninfected groups, and 27 were down-regulated genes (Table 9).

Cloning and expression regulation analysis of the novel genes
Using semi-quantitative RT-PCR, we examined the gene expression changes of the 42 potentially novel unigenes that were detected in the head kidney after viral infection. By comparing the 1, 2, 3, 4, and 5 day post-infection samples and the samples from the control group, we found four unigenes that showed a significant response to the viral infection: cichka_Cluster153 and cichka_Cluster291 were up-regulated in days 1 and 2 post-infection after which their expressions returned to the starting level; cichka_Cluster357 and cichka_Clus-ter788 were up-regulated in days 1 and 2 post-infection, and the increased expression levels were maintained till day 5 (Figure 2).
The full-length cDNA sequences of these four unigenes were 2,057 bp (cichka_Cluster291, JQ412736), 2,288 bp (cichka_Cluster357, JQ412737), 1,044 bp (cich-ka_Cluster788, JQ412738) and 1,387 bp (cichka_Clus-ter153, JQ412739) encoding polypeptides of 586, 322, 142 and 155 amino acids, respectively. BLAST searches revealed that cichka_Cluster291 can encode a protein that is similar to the vertebrate endonuclease domain containing protein, cichka_Cluster357 can encode a protein that is similar to the vertebrate ankyrin repeat domain 10 protein, cichka_Cluster788 can encode a protein that is similar to the CST complex subunit TEN1; for the cichka_Cluster153 encoded protein, no similar sequences were found in the databases that we searched, suggesting that cichka_Cluster153 may represent a novel gene in grass carp. We used the SMART server [31] to predict the domain structure of the 42 novel unigenes and found that 83.02% of them contained the endonuclease domain 1 that is found in proteins that are involved in the apoptosis pathway, and 35.22% contained the ankyrin repeat domain that is present in proteins that are involved in pathways that include the B cell receptor signalling pathway, the T cell receptor

Discussion
Currently, there are about 6,915 sequences of grass carp in the public databases. This situation does not reflect the extremely important breeding position of grass carp.
In this study, we built a head kidney non-normalized cDNA library of healthy grass carp and obtained 3,027 unigene EST sequences. This library greatly enriches the available genomic data for grass carp and lays an important foundation for the discovery of novel genes and for their functional investigation. GO analysis revealed that the annotated unigenes were mainly related to genes involved in basic biological processes such as cellular process (25.5%), metabolic process (19.1%) and biological regulation (10.8%). This functional distribution is similar to the EST distributions reported earlier in the head kidney of zebrafish [32] and sea bass [10].
Of the unigenes that were similar to immune-related genes, 66 unigenes were annotated as associated with the immune process; 53 were related to the immune system process, 4 were annotated as response to virus, and  9 were related to response to bacteria. Among the 989 unigenes that were assigned KO annotations, 68 were mapped to immune-related pathways that included leukocyte transendothelial migration, antigen processing and presentation, chemokine signalling pathway and T cell receptor signalling pathway. By examining the literature, we found that 28 of the unigenes in grass carp head kidney were related to fish genes that were reported to be involved in the Toll-like receptor signalling pathway, the RIG-I-like receptor signalling pathway and the NOD-like receptor signalling pathway. Clearly, head kidney tissue plays an important role in immune processes in fish. EST databases of head kidney tissue are likely to become important resources in which immune-related genes can be identified. In the 3,027 unigene library of head kidney in grass carp, 7.3% (221) failed to match any of the sequences in the three public databases that were searched. Of the 10 unigenes that were the most highly expressed in grass carp head kidney, 9 were unknown sequences (Table 10). This could be partly because sequence data for fish is still very scarce, and partly because fish head kidney tissue may contain tissue-specific or species-specific genes. EST databases can be important resources for identifying unknown genes in fish [33][34][35]. In recent years, the fish transcriptome has been used to study the regulation of gene expression.  conducted a comparative study of turbot expression profiles in the main immune tissue before and after pathogen infection to find genes that were related to immune response and disease resistance [36]. Chini et al (2008) carried out a comparative study of reproductive development-related tissues in bluefin tuna using transcriptome research methods to explore the molecular mechanism of gonadal development and maturity split [13]. Indeed, comparative transcriptome analysis can be used, not only to investigate the mechanisms of expression and regulation of known genes, but also as an effective means to find important and novel function-related genes.

Conclusion
We carried out a comparative analysis to find differences in the Solexa expression profiles of head kidney in grass carp before and after infection, and identified 42 unigenes of unknown function that showed differential expression in response to the pathogen. After RT-PCR validation of the cDNA and gene structure analysis, we found a potentially novel immune-related gene. Based on its response to viral infection and the prediction that it might encode a membrane protein, we speculate that this novel gene may encode a virus receptor or a protein that mediates the immune signalling pathway at the cell surface. We intend to further investigate the function of this gene in a future study. Our findings confirm that fish tissue-specific EST databases combined with comparative transcriptome analysis are effective tools that can direct the discovery of novel functional genes.