Assembly of 500,000 inter-specific catfish expressed sequence tags and large scale gene-associated marker development for whole genome association studies

Twelve cDNA libraries from two species of catfish have been sequenced, resulting in the generation of nearly 500,000 ESTs.


Background
Catfish is the major aquaculture species in the United States, accounting for over 60% of all US aquaculture production. While channel catfish (Ictalurus punctatus) accounts for the majority of commercial aquaculture production, the closely related blue catfish (Ictalurus furcatus) possesses several economically important traits that led to the production of an inter-specific hybrid (channel catfish female × blue catfish male) available for commercial use [1]. This specific hybrid shows strong heterosis and superior performance traits in disease resistance, growth rate, feed conversion efficiency, processing yields, and seinability. Channel catfish is also an important model species for the study of comparative immunology, reproductive physiology, and toxicology. The channel catfish immune system is among the best characterized of any fish species, with decades of research leading to identification and characterization of catfish immune genes [2,3], establishment of clonal functionally distinct lymphocyte cell lines [4], characterization of much of the machinery of catfish innate [5,6] and adaptive immunity and production of panels of specific monoclonal antibodies for detection of catfish immunocytes [7][8][9].
Genome research requires the development of a number of resources that facilitate both structural and functional analysis of the genome. Many of the required resources have been developed in catfish, including a large number of polymorphic markers [10,11], linkage maps [12][13][14], bacterial artificial chromosome (BAC) libraries [15,16], physical maps [17,18], and BAC end sequences (BES) [19,20]. However, expressed sequence tag (EST) resources were low from catfish [21][22][23][24][25][26], hindering both functional and comparative genome analysis in catfish. Large numbers of ESTs have been produced for most model species as well as a number of agriculturally important species [27][28][29][30][31][32], including cattle (1.5 million), swine (1.4 million), chicken (600,000), Atlantic salmon (471,000), and rainbow trout (281,000). The availability of such EST resources has allowed efficient gene discovery and gene identification in these species, and rapid progress has been made through comparative genome analysis in understanding structural, organizational, and functional properties of the genomes of these species.
A whole genome sequence is not available for most aquaculture species. In the absence of the whole genome sequence of catfish, we initiated this large-scale EST project to provide transcriptomic resources in channel catfish and blue catfish. These ESTs will serve as resources for gene discovery and gene identification, supply the framework for high-density microarray platforms, provide a foundation for the analysis of full-length cDNAs, and assist in the identification of genetic markers such as microsatellites and single nucleotide polymorphisms (SNPs). In this study, we have taken a unique inter-specific approach. The inter-specific approach will help develop markers that are inter-specific and species specific. These resources will also be of great use for comparative genome analysis. The interspecific EST approach to produce parallel EST resources from two closely related Ictalurid species will allow the resolution of some of the most difficult issues in teleost genome research, such as paralog confusions involving duplicated genomes [33][34][35]. Here we report the generation and analysis of nearly 500,000 ESTs from catfish, including 354,377 ESTs from channel catfish and 139,475 ESTs from blue catfish.

cDNA libraries and EST sequencing
Twelve cDNA libraries were constructed from various tissues, organs, and cell lines, including four blue catfish libraries and eight channel catfish libraries (Table 1). More than 600,000 sequencing reactions were attempted to sequence a total of 307,296 cDNA clones from both ends. A total of 438,321 ESTs were generated from this project, of which 128,711 sequences were from blue catfish and 309,610 were from channel catfish (Table 1). Of these EST sequences, 219,831 were sequenced from the 5' end of the transcripts, and 218,490 were sequenced from the 3' end. A total of 194,136 clones had paired reads from both 5' and 3' ends of the same transcripts. The lengths of the ESTs range from 100 to 877 bp, with an average length of 576 bp and a median length of 655 bp ( Figure 1). Addition of these sequences to the 10 Table  2). Details of the catfish EST assembly are available online [36].
In order to identify inter-specific SNPs, we also conducted the assembly of all available 493,852 ESTs from blue catfish and channel catfish. This assembly allowed the formation of 45,306 contigs (66,272 singletons), from which potential inter-specific SNPs could be identified. The number of inter-specific contigs was significantly larger than that from either species, potentially due to the formation of new contigs of related transcripts that were singletons in either species (see Discussion). A majority of contigs contained only two (43%) or three (13%) sequences (Figure 2), and average contig depth was nine sequences. With the ESTs being sequenced mostly from normalized libraries, the vast majority of contigs had 50 or fewer sequences. However, some extremely large contigs were found. The largest contig, containing 7,208 sequences, putatively identified as apolipoprotein, was repeatedly sequenced from all libraries, and was prevalent in the pre-existing non-normalized libraries in GenBank. As previously reported [37], contig depth is one of the two most important factors affecting EST-derived SNP qualities. Therefore, the information on contig depth is highly useful.
To assess the similarity between blue catfish and channel catfish sequences, the consensus sequences from each species were compared to each other using BLASTN at a stringency of 1E-10. Of blue catfish and channel catfish sequences with a minimal 200-bp matching region, sequence similarity ranged from 77% to 100%, with an average similarity of 95%. Over 50% of  blue catfish and channel catfish homologous sequences have similarity levels over 95% ( Figure 3). , and the remaining 28,380 sequences (25%) contained no ORFs (Figure 5a). These ORF-less ESTs were likely ESTs sequenced within the untranslated regions (UTRs) or within intron-retaining cDNAs.

Gene identification and annotation
There was a positive correlation between the length of ORF and BLASTX match. Of the identified ORFs, 91% had a length of more than 100 bp. Within these ORFs, 53% had significant (1E-10) BLASTX matches (Figures 5b, c). However, only 6% of the ORFs with less than 100 bp had significant BLASTX matches (Figure 5d).
A total of 41,311 (37%) unique sequences had significant BLASTX matches within the nr database, and 34,860 (31%) had significant BLASTX matches within the Uniprot database (Table 3). Over 98% of unique sequences with significant matches were identified with ORFs, which indicated the reliability of ORF searching. After examination of putative protein identities from the BLASTX searches, homologous sequences were identified from the catfish ESTs. Of the 41,311 sequences with BLASTX hits, 22,642 (approximately 55%) and 17,948 (approximately 43%) unique proteins were identified through searches against the nr and the Uniprot protein databases after removing the redundant protein hits, respectively.

Assessment of the sequenced catfish transcriptome
In order to assess the level to which the catfish transcriptome has been captured, the unique catfish sequences (111,578) were also searched against the NCBI Refseq and Ensembl databases. A number of significant hits were identified within zebrafish, medaka, Tetraodon, human, mouse, and chicken reference protein databases (Table 3). After removal of the redundant protein hits, 14,988 -11,059 unique reference proteins were identified within zebrafish, medaka, Tetraodon, human, mouse, and chicken databases respectively ( Table 3). The unique catfish sequences had hits to 54% to 57% of the unique proteins of zebrafish, medaka, and green-spotted pufferfish. To allow comparison of catfish unique protein coverage with that expected between species with complete genome sequences, all Tetraodon Ensembl proteins were searched against the medaka Ensembl protein database. A total of 22,150 Tetraodon proteins have significant hits to 15,054 (61% of total unique) medaka Ensembl proteins with a cutoff E-value of 1E-10. Similarly, zebrafish Refseq proteins were searched against the human Refseq protein database. In this case, 24,971 zebrafish proteins have significant hits  to 13,789 (36% of total unique) human proteins with a cutoff E-value of 1E-10. Taken together, these numbers provide strong evidence that this project has captured a large majority of the catfish transcriptome. A total of 14,776 cumulative unique genes were identified from catfish based on BLASTX searches against the Refseq/Ensembl database (Table 3), including 8,075 genes identified from both blue catfish and channel catfish, 1,881 unique genes from blue catfish and 4,820 unique genes from channel catfish ( Figure 6). As expected based on sequencing coverage depths, significantly more unique genes were identified from channel catfish than blue catfish.
To assess the evolutionary conservation of the identified unique genes, the number of hits to unique genes in each species of zebrafish, medaka, Tetraodon, human,   mouse, and chicken were compared. A total of 8,592 (58% of total number of unique catfish genes) putative known unique genes were found in all six species: 11,303 (76%) were found in all three fish species; and 14,515 (98%) were found in at least one of the three fish species (Figure 7), indicating a high level of conservation of gene content among catfish and other teleost fish species.

Prediction of full-length cDNAs
The catfish EST sequences provide a platform for the identification and characterization of full-length cDNA clones without having to use expensive and labor-intensive primer walking sequencing. In the context of this work, full-length cDNA inserts were defined as a cDNA from a single clone with the start codon and poly (A) tail contained within the clone. A total of 10,037 channel catfish and 7,382 blue catfish putative full-length cDNAs were identified from the assembly with a cutoff E-value of 1E-5. A well characterized full-length cDNA set from catfish will be crucial in ongoing studies of teleost gene duplication and gene family structure, as well as aiding in annotation of the catfish whole genome sequence. Current efforts are focused, therefore, on characterization and re-sequencing of these full-length cDNAs.

Microsatellite and SNP marker identification
A total of 20,757 microsatellites were initially identified from 15,082 unique sequences, including di-, tri-, tetra-, penta-and hexa-nucleotide repeats (Table 4). After removing the microsatellites without enough flanking sequence for primer design, 13,375 unique sequences with microsatellites had sufficient flanking sequences (50 bp) on both sides of the microsatellites to design primers for genotyping. Our previous research indicated that over 72% of EST-derived microsatellites were polymorphic in one resource family [12], suggesting the potential utility of these microsatellites. A total of 48,702 putative SNPs and 14,803 putative insertions/deletions (indels) were identified from the blue catfish EST dataset assembly; 102,252 putative SNPs and 41,660 putative indels were identified from channel catfish EST dataset assembly (Table 5). These putative SNPs indicated an SNP rate of 3.2 SNPs per kilobase of transcribed sequences in blue catfish, and 4.1 SNPs per kilobase of transcribed sequences in channel catfish. Obviously, such SNP rates were calculated from the total consensus sequence length and, therefore, the deeper the EST sequencing was, the greater the possibility for the identification of a SNP within the consensus sequences.
Over 303,000 putative SNPs and 100,000 putative indels were identified from the all catfish EST assembly results (Table 5). EST-derived SNPs are often prone to sequencing errors. Therefore, the putative SNPs were subjected to filtering using only those with contig sizes of at least four sequences and the minor allele presence of at least two sequences in the contigs, and indels were not used for further analysis [37]. These parameters were previously shown to select markers with 70% success rate for genotyping. After filtering, 3,891 and 15,827 SNPs were identified from the blue catfish and channel catfish EST dataset assemblies, respectively. A subset of 48,594 filtered SNPs were obtained from the all catfish EST assembly; these SNPs included 32,235 transitions and 16,359 transversions ( Table 5). The filtered SNP frequency in the transcribed sequences was 0.25 SNP in blue catfish, 0.64 SNP in channel catfish, and 1.6 SNP in the all catfish assembly per kilobase. A total of 19,398 filtered insertions and deletions (Indels) were discovered, that is, 0.64 indels per kilobase of the transcribed sequences. Of the 48,594 SNPs, over 90% were identified from the contigs containing 5 or more sequences ( Table 6).
The assessment of the rates of inter-specific SNPs and intra-specific SNPs may have practical applications. We therefore assessed these SNP rates using the EST data. First, SNPs were identified from contigs containing at least four sequences with at least two sequences from either channel catfish or blue catfish in the all catfish EST assembly. Inter-specific SNPs were defined as those that have sequence variations between blue catfish and channel catfish, but no sequence variations within blue catfish or within channel catfish; similarly, SNPs were identified within blue catfish but not within channel catfish or vice versa; and SNPs were identified within both channel catfish and blue catfish at the same SNP positions ( Figure 8). Of the 48,594 filtered SNPs, 42,080 were identified from contigs comprising both channel catfish and blue catfish ESTs, and 6,514 were identified from contigs composed of ESTs from either channel catfish or blue catfish, including 5,396 from channel catfish contigs and 1,118 were identified from blue catfish contigs. Of the contigs containing ESTs from both blue catfish and channel catfish, the estimation of percentage of inter-and intra-specific SNPs was conducted based on the identification of SNPs from 1,000 randomly  selected contigs. Of the 48,594 filtered SNPs identified from the all catfish assembly, over 18,000 (39%) were inter-specific SNPs; with 523 (1%) intra-specific SNPs for blue catfish, 2,352(5%) intra-specific SNPs for channel catfish, and 3,790 (8%) intra-specific SNPs for both channel catfish and blue catfish. However, approximately 17,000 SNPs could not be determined because overall the SNPs qualified as SNPs with at least two minor allele sequences, but only one of the minor allele sequences was from one of the two species of blue catfish or channel catfish (Table 7). Additionally, the number of inter-specific SNPs may be overestimated, due to failure to capture minor allele sequences from one or both species in the current EST data. However, the sequence differences between species should be greater than those within species. Although a large number of filtered inter-specific SNPs were identified (18,000 out of 48,000 total filtered SNPs), they were identified from a relatively small number of contigs. The 18,000 filtered inter-specific SNPs were identified from approximately 2,800 contigs, with an average of 6.6 SNPs per contig.

Discussion
This project represents one of the major milestones in catfish research, and brings the catfish EST resources to almost a half million sequences in GenBank [21][22][23][24][25][26]. This EST resource will prove useful for gene discovery, molecular marker development, and genetic linkage and comparative mapping, and it will help facilitate whole genome sequencing and annotation. Parallel EST sequencing in two closely related species, I. punctatus and I. furcatus, may also provide the material basis for the analysis of genome duplication and genome evolution, providing the basis for establishment of orthologies through phylogenetics analysis. The most important outcome of EST sequencing is gene discovery. This project allowed identification of 70,717 unique sequences in channel catfish and 54,815 unique sequences in blue catfish. We also conducted   Total SNP NA 48,594 100 *SNPs were identified from contigs containing at least four sequences with at least two sequences from either channel catfish or blue catfish in the all catfish EST assembly: 1 where there were no intra-specific blue catfish SNPs or intra-specific channel catfish SNPs, but the sequence differed between the two species at the inter-specific SNP position; 2 where there were SNPs within blue catfish, but not within channel catfish; 3 where there were SNPs within channel catfish, but not within blue catfish; 4 where there were SNPs within both blue catfish and channel catfish; 5 undetermined because overall the SNPs qualified as SNPs with at least two minor allele sequences, but only one of the minor allele sequences was from one of the two species of blue catfish or channel catfish; 6 these SNPs were identified from ESTs that have been sequenced from only one of the two species, blue catfish or channel catfish, to date. Examples of SNPs whose categories could not be determined because the minor allele sequence from a given species is fewer than two.
EST assembly using ESTs from both channel catfish and blue catfish. Assembly of all the catfish ESTs resulted in 111,578 unique sequences. Comparison of channel catfish and blue catfish coding regions in this study indicated that the two species share, on average, 95% sequence identity. Therefore, combining genes identified from both species should provide a more complete picture as to what fraction of the catfish transcriptome was captured to date. Such an approach was taken also because of practical considerations. Hybrid catfish produced by inter-specific hybridization of channel catfish × blue catfish is one of the best production lines of catfish used in aquaculture, and many believe that industry-wide application of this hybrid may have a revolutionary impact on the catfish industry. One of the major catfish breeding programs is based on introgression of beneficial genes from blue catfish into channel catfish breeds. Genetic linkage mapping has been conducted in both the intra-specific resource families involving only channel catfish [14] and the inter-specific resource families made from backcrosses of the channel catfish × blue catfish hybrids [12,13]. Given the close phylogenetic relationship of blue catfish and channel catfish, we expected that many of the contigs from the blue catfish and channel catfish EST assembly would merge together in an all catfish EST assembly. However, the all catfish EST assembly generated 45,306 contigs, a much greater number than the contigs generated in either the blue catfish (22,009) or channel catfish (28,941) EST assembly. There could be several reasons for this major increase in contig numbers with the all catfish EST assembly. First, some ESTs belonging to the contigs were only sequenced in blue catfish but not in channel catfish, and vice versa; second, singletons in either blue catfish or channel catfish were brought together to form new contigs; third, sequence variations or splicing differences between the two species may have led to the formation of a larger number of contigs under our assembly parameters; fourth, ESTs derived primarily from transcript untranslated regions of the two species may differ sufficiently to prevent placement in the same contig. Thorough analysis of the dataset has revealed that all four of these factors contributed to the high number of contigs in the all catfish assembly.
Analysis of the all catfish unique sequences showed that a large proportion of the catfish transcriptome has been captured. BLASTX searches identified 37% of total transcripts with significant hits, similar to levels reported in the salmon EST project [30]. The 111,578 unique catfish sequences had hits to 54% to 57% of the unique proteins of zebrafish, medaka, and green-spotted pufferfish using a cutoff value of E-10. This percentage appeared at first glance to be lower than our expectations. We therefore carried out best-hit searches using identical parameters as those used with catfish but comparing protein coverage of species with complete genome (transcriptome) sequences to serve as reference points. We found that the Tetraodon protein set had significant hits to only 61% of medaka proteins. By comparison, our catfish data set had significant hits to 56% and 57% of medaka and Tetraodon proteins, respectively. Similarly, zebrafish Refseq proteins were searched against the human Refseq protein database. In this case, zebrafish proteins had significant hits to 36% of total unique human proteins, compared to 33% in catfishhuman alignments (Table 3). These reference numbers indicate both the high coverage of the catfish transcriptome obtained in this project, and the limitations of simple homology searches given the rapid divergence of many genes following speciation and the complexity of genome-wide and local gene duplication events within teleost species. Given that the identity of only 37% of the unique catfish sequences could be characterized by homology searches, the utility of the dataset should increase significantly with whole-genome sequencing of Ictalurid catfish and additional sequencing in closely related species within the order of Siluriformes. Interestingly, over 40,000 unique catfish sequences containing an ORF did not have a significant hit by homology searches. Further work will be needed to characterize whether low homology rates in these sequences is due to short read length, the rapid evolution of the encoded gene, or 'catfish-specific' gene duplication and divergence.
Large-scale EST sequences provide an enormous resource for molecular marker development. This project allowed identification of over 20,000 microsatellites within ESTs, of which 13,375 were located within unique ESTs and had sufficient flanking sequences for microsatellite primer design for genotyping (Table 6). Therefore, these microsatellites will be a major resource for genetic linkage and comparative mapping [12]. In addition, over 300,000 putative SNP sites were identified, of which over 48,000 were identified from contigs with at least four ESTs and the minor sequence was represented at least twice ( Table 7). The 48,000 filtered SNPs should be highly useful for the development of a SNP panel for whole genome association studies [37].
The parameters of quality SNP assessment may not be applied to the very large contigs. The utilization of a minor allele frequency of six for all the contigs containing 30 sequences or more resulted in higher SNP frequency from these contigs, such as 13.4 SNPs per kilobase in the contigs with 100 to 500 sequences, and 15 SNPs per kilobase in the contigs with 500 sequences or more. Information regarding contigs over 500 sequences can be found in Additional data file 1. High SNP frequency within these large contigs may be caused by the accumulation of sequencing errors or alignment of transcripts from multi-copy loci, so SNPs from large contigs will be avoided in future SNP genotyping.

Conclusions
In this project, generation and assembly of channel catfish and blue catfish ESTs allowed the identification of 45,306 contigs and 66,272 singletons, and a large majority of the catfish transcriptome was captured. Whole genome sequencing of channel catfish and blue catfish is currently underway, and the comparison between genome and transcriptome sequences will enable better understanding of the gene structure and organization. The analysis of the inter-specific ESTs resulted in the identification of 20,757 gene-associated microsatellites and over 300,000 putative SNPs, of which over 48,000 were filtered SNPs with the presence of the minor allele at least twice. These SNPs have been utilized to design the first generation high-density SNP chips using Illumina iSelect HD SNP genotyping panels for genome association studies. The inter-and intra-specific SNPs identified from the all catfish EST dataset assembly will greatly benefit catfish introgression breeding selection and whole genome association studies.

Materials and methods
cDNA library construction, EST sequencing and processing The cDNA libraries were constructed by consortium investigators using various tissues, organs, and cell lines, including stomach, muscle, olfactory tissue, trunk kidney, head kidney, gill, intestine, spleen, skin, liver, pituitary, ovary and testis (Table 1). Total RNA was isolated from experimental tissues, reverse transcribed using an oligo-dT primer, directionally cloned into either the pSPORT-1 (Invitrogen Corp., Carlsbad, CA, USA) or pDNR-Lib (Clontech Laboratories Inc., Mountain View, CA, USA) plasmid vectors, and electroporated into competent Escherichia coli. One library (CBPN) underwent subtraction for highly expressed clones, ten libraries were normalized, and one library (CBCZ) was processed without normalization. Clone selection, arraying, and sequencing of all 12 libraries were performed at the US Dept. of Energy -Joint Genome Institute. Both ends of the insert were sequenced using Big Dye Terminator (V3.1) chemistry (Applied Biosystems, Foster City, CA, USA), and low quality sequences were trimmed. Contaminant sequences (E. coli, mitochondrial, cloning vector, rRNA, tRNA) were filtered.

EST assembly
Three separate assemblies were performed: blue catfish ESTs, channel catfish ESTs, and all catfish ESTs. The new EST sequences and existing EST sequences from channel catfish and blue catfish were clustered and assembled using the Paracel Transcript Assembler, based on the CAP3 assembler [38]. Repeat sequences and poly (A) tails were masked and annotated. Prior to assembly, all EST sequences were compared to 'seed' sequences, which were existing catfish full-length or partial cDNA sequences in GenBank. New sequences sharing 80% similarity to seed sequences were clustered and assembled at 95% identity with at least a 50-bp overlap to generate seed-cluster contigs. The seed cluster assembly reduced the number of sequences for final assembly in order to minimize computational requirements. The remaining EST sequences were then clustered based on local similarity scores of pairwise comparisons with a minimum 88% similarity of at least 100 bp. Clusters containing only one sequence were denoted as singletons. The EST clusters were assembled into contiguous sequences (contigs) by multiplesequence alignment with 95% identity of at least a 50-bp overlap, and a consensus sequence was generated for each cluster. Multiple contigs could be generated from each cluster, since EST clusters may not share enough similarity over their entire length to be assembled as a single contig. Multiple contigs could also be generated when ESTs in the cluster represented splice variant forms or paralogs. Single ESTs remaining in a cluster after the formation of contigs were designated as cluster singletons. The unique sequences for each assembly included the seed-cluster contigs, cluster contigs, cluster singletons, and singletons. All the sequence assemblies are available upon request to the corresponding author.

ORF searching, gene identification and gene ontology annotation
All unique sequences obtained after the assembly were analyzed by ESTScan [39] to search for ORFs, which could be used to distinguish coding and non-coding sequences [39,40]. The putative protein sequences were also generated at the same time by ESTscan, which could be used to analyze splice variation, determine paralogs, and assess gene families. All unique sequences were compared against the nr and Uniprot databases using BLASTX (cutoff E-value of 1E-10) to obtain the putative identity. The NCBI Refseq protein and Ensemble databases (zebrafish, medaka, Tetraodon, human, mouse, and chicken) were also used to annotate the unique catfish genes.

Full length cDNA identification
Putative full-length cDNAs were identified by comparison to full-length genes and start signals in Uniprot databases using TargetIdentifier [34,41] with a cutoff Evalue of 1E-5. Once the start codon (ATG) and poly (A) tail were identified, the cDNA sequence was considered a full-length cDNA.

Microsatellite and SNP marker identification
All the unique sequences were used to search for microsatellite makers using Msatfinder [42] with a repeat threshold of eight di-nucleotide repeats or five tri-, tetrapenta-, or hexa-nucleotide repeats. Clones containing 50bp sequence on both sides of the microsatellite repeat were considered sufficient for primer design [43].
All three assemblies were used for SNP identification using autoSNP [44]. The parameters for minimum minor allele frequency for SNP detection varied with the number of sequences in the contig. A sequence variation was declared a SNP when: a mismatch was identified in contigs with four or fewer sequences; the minor allele sequence existed at least twice within contigs containing 5 to 6 sequences; the minor allele sequence existed at least three times within contigs containing 7 to 8 sequences; the minor allele sequence existed at least four times within contigs containing 9 to 12 sequences, or the minor allele sequence existed at least five times within contigs with 13 or more sequences. One thousand contigs containing filtered SNPs were randomly selected to inspect the inter-and intra-specific SNP calls.