Generation and analysis of ESTs from the eastern oyster, Crassostrea virginica Gmelin and identification of microsatellite and SNP markers

Generation of ESTs and contig assembly

The ESTs generated from this study were searched against ESTs of C. virginica in dbEST using BLASTN. Of the 4,688 unique sequences generated in this study, 901 (19.2%) had significant hits (e ≤ 10^-5) with existing ESTs from C. virginica, while 3,787 (80.8%) were found to be novel C. virginica ESTs. Therefore, this EST collection represented a significant addition to the existing oyster EST resources.

Putative identities of the ESTs

In order to make an assessment for the putative identities of the sequenced ESTs, all ESTs were subjected to BLASTX similarity searches. Of the 4,688 unique sequences, 2,174 (46.4%) had significant matches (E ≤ 10^-5) to the non-redundant (nr) protein database using BLASTX (Additional file 1). Of the unique sequences with matched homologies in the nr protein database, seven (0.32 %) had an E-value of 1e-100 or less (Figure 1) and are therefore considered highly significant homology [28]. About 58% had significant homology (E-values between 1e-20 and 1e-99) and thus, are considered moderately similar [28, 29]. About 42% of the ESTs were assigned weak homology (E-values between 1e-05 and 1e-19) and are therefore weakly similar. Among the various organisms that have protein sequences in GenBank, the C. virginica ESTs generated from this study had the highest number of BLASTX hits to the purple sea urchin Strongylocentrotus purpuratus (16%), followed by zebrafish Danio rerio (10%), honey bee Apis melifera (6%), and chicken Gallus gallus (6%) (Figure 2). Among the organisms listed in "other" (43% of total), 13 ESTs (0.6 %) had significant hits to C. virginica proteins in the database and 35 ESTs (1.6%) had hits with the Pacific oyster, Crassostrea gigas.

Of the 2,514 unique sequences without significant matches to the nr protein database, 6 (includes 4 contigs consisting of 57 ESTs and 2 singletons) had matches using the BLASTN algorithm (E-value ≤ e-5). All of these correspond to mtDNA sequences.

Gene ontology annotation

Gene Ontology (GO) categories were assigned to 1,104 unique ESTs with BLASTX hits using Blast2GO. Figure 3 shows the percentage distributions of gene ontology terms (3^rd level GO terms) according to the GO consortium. Cellular physiological process (80 %) was the most dominant 3^rd level term out of the 761 unique sequences which were annotated to the Biological Process GO category. This was followed by metabolism at 62%. Only 4% were assigned to the Biological Process subcategory response to stress. Protein binding (26%) was the most dominant out of 952 ESTs with significant protein hits which were assigned to Molecular Function category at 3^rd level. This was followed by hydrolase activity and nucleotide binding at 17% each. GO terms at higher and lower levels for each of the three main GO categories as well as the unique ESTs which fall on each term are given in additional file 2.

Microsatellites-containing ESTs

Single nucleotide polymorphism

Tissue source and RNA isolation

Oysters from Rutgers University (NEHY-53) were used in this study. Samples of tissues (50 mg from 4 oysters for each tissue type) were excised from gill, mantle, digestive tract and gonad, and adductor muscle. In addition, about 100 mg of tissue was obtained from one whole juvenile oyster and hemocytes were extracted from 12 oysters. Total RNA was prepared from each sample by using RNeasy Mini kits (Qiagen, Valencia, CA). Total RNA quality and quantity were checked by agarose gel electrophoresis containing formaldehyde and by using a spectrophotometer. Equal amounts of total RNA from each sample listed above were mixed for cDNA library construction.

Construction of normalized and subtracted cDNA library

The Creator™ SMART™ cDNA Library Construction Kit (Clontech, Mountain View, CA) and components from the TRIMMER-DIRECT Kit from Evrogen (Moscow, Russia) were used for the construction of the normalized cDNA library. One microgram (1 μg) of mixed total RNA was combined with 1 μl of SMART IV oligonucleotide (Clontech) and 1 μl CDS-3 M adapter (Evrogen) for first-strand cDNA synthesis. The reaction was incubated at 72°C for 2 min followed by immediate cooling on ice for 2 min. Then, 2 μl of 5X first strand buffer, 1 μl of DTT (20 mM), 1 μl of dNTP mix (10 mM), and 1 μl of PowerScript reverse transcriptase were added to the tube and incubated at 42°C for 1 hr in a thermal cycler (PTC-100, Bio-Rad, Hercules, CA) and placed on ice. The SMART cDNA cloning system allows the enrichment of full-length cDNA through the use of a 5'-linker with 3'-GGG tails. Reverse transcriptase has terminal transferase activity that preferentially adds three additional Cs at the end of first strand cDNA. As a result, the first strand cDNA is able to base pair with the 5'-linker with 3'-GGG tails. Once base paired, the reverse transcriptase switches templates and extends into the linker sequences allowing PCR amplification of full-length cDNA using a single primer (the 5'-linker has the same sequences as the linker containing poly-T used for the synthesis of the first strand cDNA). Truncated cDNAs are not able to base pair with the 5'-linker and, therefore, get lost in the PCR amplification of the full-length cDNA. The first strand cDNA was initially amplified by long-distance PCR (LD-PCR) using hot-start amplification. For the reaction, the following were combined in a reaction tube: 1.5 μl of the first strand cDNA, 60 μl of sterile deionized water, 7.5 μl of 10X Advantage 2 PCR buffer, 1.5 μl of 50X dNTP mix, 3 μl of 5' PCR primer and 1.5 μl of 50X Advantage 2 polymerase mix. The tube was mixed and briefly centrifuged and added to a pre-heated (95°C) thermal cycler. Cycle settings were 95°C for 1 min followed by 19 cycles of 95°C for 7 sec, 66°C for 20 sec, and 72°C for 5.5 min. The product was analyzed on a 1.1% agarose gel to determine the sizes and amount of the cDNA products before proceeding to the next step. The LD-PCR reaction was purified and eluted in 30μl of sterile nanopure water using the QIAquick PCR Purification Kit (Qiagen). For the normalization procedure, the TRIMMER-DIRECT Kit from Evrogen (Moscow, Russia) was used. This system is specially developed to normalize cDNA enriched with full length sequences. The cDNA from the LD-PCR was quantified (~200 ng/μl) and 1 μl was mixed with 1 μl of 4X hybridization buffer, 1 μl of sterile water, and 1 μl (~500 ng/μl) of PCR fragments of 135 abundant ESTs (see additional file 3 for the list) from a previous oyster cDNA library developed in our lab [12]. The mix was overlaid with mineral oil and incubated for 3 min at 98°C then at 70°C for 4 hr. Then, 5 μl of 2X DSN buffer (preheated to 70°C) and 0.25 Kunitz units of DSN enzyme were added and incubated at 70°C for 20 min. The DSN enzyme specifically degrades double-stranded molecules. The reaction was inactivated by adding 10 μl of DSN stop solution, and sterile water added to a final volume of 40 μl. Following normalization, two rounds of PCR were performed using 1 μl of the normalization reaction as template. A shorter primer M1 (first 23 bases of the SMART IV oligonucleotide) was used in the first round of PCR with 15 amplification cycles using the same thermal cycling parameters as above; and an even shorter primer M2 (first 20 bases of the SMART IV oligonucleotide) was used in the second round of PCR for 15 amplification cycles of 95°C for 7 sec, 64°C for 20 sec, and 72°C for 5.5 min. Products were checked on a 1.1% agarose gel. The PCR products were quantified and 3 μg were used for treatment with proteinase K. All the subsequent procedures including proteinase K treatment, restriction digestion with Sfi I, size fractionation, and ligation followed the manufacturer's instructions (Clontech). The cDNA was ligated to the pDNR-LIB vector. Electroporation (MicroPulser, Bio-Rad, Hercules, CA) was performed using DH12S electrocompetent cells following supplier's instructions (Invitrogen). A total of about 8 × 10⁵ primary recombinant clones were obtained, the average insert size was 1650 bp, and the library was amplified, titred, and stored in 25% of glycerol stocks in a -80°C freezer. When needed, clones were plated on Luria-Bertani (LB) agar medium containing (50 μg/ml chloramphenicol) and grown overnight at 37°C. Then, colonies were randomly picked, inoculated in 384-well microtitre plate (containing LB medium, 50 μg/ml chloramphenicol, and 10 % glycerol), incubated overnight with shaking at 37°C, and stored in a -80°C freezer until further use.

Plasmid isolation and DNA sequencing

Clones were transferred from 384-well plates to 96-well plates containing LB medium with 50 μg/ml chloramphenicol and grown for about 20 hr prior to plasmid isolation. Plasmid DNA was isolated from 6,528 randomly selected clones using Perfectprep^® Plasmid 96 Vac, Direct Bind Kit (Eppendorf). Sequences of cDNAs were sequenced from their 5' end using Big Dye terminator and M13 primer (5'TGTAAAACGACGGCCAGT3') on an ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA) following manufacturer's protocol.

EST processing, contig assembly and analysis

The chromatogram files were exported to the PHRED program [36] for basecalling and removal of poor quality sequences. Then, vector sequence, adapter region, and poly(A) tails were trimmed from the sequence using Vector NTI Advance™ 10 (Invitrogen Corporation, 2005). Trimmed sequences were further screened by including the vector sequence in the contig assembly using the ContigExpress in Vector NTI Advance™ 10 (Invitrogen Corporation, 2005) and discarding all ESTs that formed a contig with the vector. High quality ESTs (at least 100 bp after removal of vector sequence, adapter, and poly(A) tail) were then assembled into clusters of contiguous sequences (contigs). Three software programs were used for contig assembly, namely, PHRAP, CAP3 [37], and Vector NTI Advance™ 10. After trying the default parameters of the three programs, varying the stringencies of the different parameters, and manually examining the contigs, the results for Vector NTI (using more stringent parameters, i.e, overlap length cutoff of 50 and overlap percent identity of 90) were retained for further analysis. The ContigExpress module in Vector NTI Advance™ 10 also makes use of CAP3 to drive the assembly process (Invitrogen Corporation, 2005). The consensus sequence of each contig and singletons comprised the unique sequences which were compared against the National Center for Biotechnology Information (NCBI) nonredundant protein database using BLASTX. ESTs that did not match sequences in the protein database were further analyzed by searching for similarities at the nucleotide level by using BLASTN. Novel ESTs were also identified by comparison with C. virginica EST sequences in dbEST at NCBI using BLASTN. The E-value cutoff was 1e-5.

Gene ontology annotation

Sequences with BLASTX hits were mapped and annotated according to gene ontology terms (GO) using the program Blast2GO [38]. The distribution of genes in each of the main ontology categories was examined and the percentages of unique sequences in each of the assigned GO terms were computed. In each of the three main categories of GO, namely, biological process, molecular function, and cellular component [39], 100 % was considered as the total number of unique sequences having an assigned GO term. Thus, in each main category, the percentages do not add up to 100% because some deduced proteins have more than one GO category assigned to them [40].

Bioinformatic mining of microsatellites and SNPs

The set of 4,688 unique sequences was searched for microsatellites using the web-based software, Msatfinder v. 2.0.9 [41]. The thresholds or minimum repeat number used for the search was 8 for dinucleotide microsatellites, 6 for trinucleotide microsatellites, and 5 for tetra- and pentanucloetide microsatellites. Microsatellite-containing ESTs were identified as candidates for marker development if they contained sufficient flanking sequences on either side of the repeats for primer design using Msatfinder. All existing C. virginica ESTs in NCBI's dbEST were added to the ESTs generated in this work for detection of putative SNPs using the AutoSNP program with default parameters [42]. Thus, a total of 14,560 C. virginica ESTs were searched for putative SNPs.

Accession numbers

The ESTs generated from this study have been deposited in GenBank and were assigned accession numbers from EH643873 to EH649414.