Sequencing, Analysis, and Annotation of Expressed Sequence Tags for Camelus dromedarius

Despite its economical, cultural, and biological importance, there has not been a large scale sequencing project to date for Camelus dromedarius. With the goal of sequencing complete DNA of the organism, we first established and sequenced camel EST libraries, generating 70,272 reads. Following trimming, chimera check, repeat masking, cluster and assembly, we obtained 23,602 putative gene sequences, out of which over 4,500 potentially novel or fast evolving gene sequences do not carry any homology to other available genomes. Functional annotation of sequences with similarities in nucleotide and protein databases has been obtained using Gene Ontology classification. Comparison to available full length cDNA sequences and Open Reading Frame (ORF) analysis of camel sequences that exhibit homology to known genes show more than 80% of the contigs with an ORF>300 bp and ∼40% hits extending to the start codons of full length cDNAs suggesting successful characterization of camel genes. Similarity analyses are done separately for different organisms including human, mouse, bovine, and rat. Accompanying web portal, CAGBASE (http://camel.kacst.edu.sa/), hosts a relational database containing annotated EST sequences and analysis tools with possibility to add sequences from public domain. We anticipate our results to provide a home base for genomic studies of camel and other comparative studies enabling a starting point for whole genome sequencing of the organism.


Introduction
Camelus dromedarius, often referred to as the Arabian camel, is probably the most famous member of the camel family. Other members of the camel family include the llama and the alpaca in South America. The Dromedary has one hump on its back, in contrast to the Bactrian camel which has two. The Dromedary is more numerous than the Bactrian camel and represents almost 90% of the genus Camelus. Camel has historically and economically been an important species worldwide and especially in Arab Peninsula. Saudi camels comprise 16% of the animal biomass. Despite their differences, camels do have their shape in common. These differences are among breeds or within breeds. No clear classification for camels exists but generally they can be classified according to habitat, color and function. Camel breeds vary in size, body conformation and color. Color is the most common character used for classification of camel breeds. Some are dark black and others have white or brown colors. Based on their colors, three main breeds of Saudi camels were distinguished, namely black (Magaheem), white (Magateer) and brown (Al Homr and Al Sofr) [1]. Camels are multipurpose animals with females used primarily as milk producers, the males for transport or draught and both sexes providing meat as tertiary product. Camel stores its energy reserves in the form of fat in different depots in their body of which the hump and abdomen depots comprise a considerable amount of the adult body weight; therefore camels can survive long periods without feed. Hump and abdomen fats contain mixtures of fatty acids [2] and most of these are esterified as triglycerides or phospholipids and vary according to their anatomical location in the body [3]. The demand for camel meat appears to be increasing especially in arid regions. Camel meat is healthier as they produce carcasses with less fat as well as having less levels of cholesterol in fat than other meat animals [4]. Camel offers interesting biological perspectives such as heavy-chain antibody based humoral immune system [5], existence of naturally occurring domain antibodies (dAbs) (smallest known antigenbinding fragments) [6], and alternative insights in reproductive biology [7]. Moreover, camel's characteristic ability to adapt its desert lifestyle with remarkable traits such as fluctuating its body temperature from 34 degrees Celsius to 41.7 degrees Celsius throughout the day, tolerating a water loss greater than 30%, and capability of drinking 100 liters of water in as little as ten minutes [8] demands further biological investigation. Camels are usually raised in the harsh conditions of the desert, where water is scarce. These animals are well adapted to dehydration for relatively long periods [9]. Basal plasma glucose levels are significantly higher in monogastrics than in adult ruminants, but not in camels [10]. Camels are exceptional with regard to their carbohydrate metabolism. On the one hand, plasma glucose levels in camels are in the range or of those of monogastrics or even higher [11], although camels ferment carbohydrates to short-chain fatty acids on the same scale as sheep and lactating cows [12] and like ruminants, they also meet their glucose demands by endogenous gluconeogenesis [13]. On the other hand, the whole body insulin sensitivity of camels is even lower than that of adult ruminants [14]. High basal plasma glucose levels and simultaneously low insulin sensitivities also occur in humans suffering from noninsulin-dependent diabetes mellitus (NIDDM) [15], so the insulin sensitivity of camels may be similar to that of NIDDM patients.
Despite its economical, cultural, and biological importance, there has not been a large scale sequencing project targeting camel genome. We initiated Camel Genome Project with the goal of sequencing complete DNA of the organism. Since there is not much available information about the camel genome, we first established and sequenced camel EST libraries, to increase efficiency and accuracy of whole genome sequencing. A successfully completed EST library project can supply wealthy genetic information for a species, often considerably shortening the time-consuming and laborious gene isolation and characterization procedures. Large scale EST projects do not only provide direct information on the transcriptome and indirect information on correlation of genome and phenotype but also most often are very crucial to whole genome research in a given species as we have witnessed in most of mammalian genome projects. EST libraries are also very useful for molecular marker development in linkage mapping, comparative genome analysis, and estimation of gene duplication [16]. This project supplies the building backbone for comparative genomics providing a reference point from which we can compare camel with other organisms. The project helps us begin to understand the molecular basis for differences in phenotypes between camel breeds, correlation of SNPs among different camels in health and disease, disease-susceptibility prediction based on gene sequence variation, mapping and identifying genes involved in complex biological traits and multigene diseases, obtain camel specific genes to develop vaccines and/or treatments for common diseases. It can also help in the production of useful protein products for use in bioremediation and pharmaceutical industries, protein replacement (e.g. factor VIII, TPA, streptokinase, insulin, interferon), and ''Pharm'' animals.
In this study, we have sequenced 70,272 ESTs from three camel cDNA libraries and established a publicly available database for data mining and visualization. In order to obtain a catalogue of gene sequences found in camel irrespective of phenotype or abundance, we used normalized cDNA libraries from different age, breed, and tissue samples. Following data cleanup, clustering and assembly steps, we generated 23,602 non-redundant gene indices. We compared our findings in particular to the following nine species: Homo sapiens, Mus musculus, Rattus norvegicus, Bos Taurus, Sus scrofa, Pan troglodytes, Macaca mulatta, Canis familiaris, and Equus caballus. We chose human, mouse, and rat as model organisms that have been heavily worked on and have well annotated genomic data. Horse, bovine, and boar were chosen as these species are ungulates (like camel) with substantial sequence data and share many physiological properties with camel. Finally, we included monkey, chimpanzee and dog in order to lead way for comparative genomics regarding mammalian species that have found significant attention by research community in a variety of ways including sequence and genomic studies. To our knowledge data presented here represent first EST sequencing in camel and we believe accompanying web portal will become an essential tool for the annotation and assembly of camel's whole genome. Our work may have significant impact on mammalian functional genomic research since we have described over 4,500 sequences, which do not carry any homology to any other genome that is available.

Ethics Statement
International rules on Experimentation on Living Animals have been strictly implemented and consent has been obtained from the Saudi National Commission for Experimentation on Live Animals to carry out the relevant procedures on camels prior to the commencement of experimental work. Certification for this project has been obtained from the Saudi National Committee of Bio and Medical Ethics (Reference No: 229/131).

Tissue Collection
Nine inbred camels of three age groups (young (0-6 months), adult (2-3 years), and aged (4-6 years)) and three different breeds (white, black, and brown coat color; determined by expert breeding techniques and visual inspection) have been used. Camels were housed in the Camel Facility at King Faisal University at Hassa, KSA. Camels were fed on alfa alfa (Medicago sativa) and standard concentrate diet with free access to water. All camels were humanely sacrificed by a high barbiturate dose (thiopentone sodium) for collection of tissues. We collected and pooled 11 tissues (Brain, Liver, Kidney, Heart, Muscle, Lung, Spleen, Pancreas, Stomach, Genitals, and Colon) from each of the nine camels. Pieces of dissected tissues were placed into the RNAlater solution (Ambion, USA) at room temperature. The solution permeates the cells, stabilizing the RNA. The samples were then stored at 4uC for two hours and then stored at 280uC until analyzed. In addition, blood samples were collected into Paxegen (Qiagen) tubes for RNA isolation. Since the goal of this study is not to interrogate gene expression in camel based on age, breed, tissue, or disease state but rather obtain a catalogue of gene sequences found in camel, we chose to pool different age, breed, and tissue groups to generate a collective gene expression data set in camel. This strategy attempts to capture genes expressed based on phenotypic differences. In order to obtain a most representative cDNA library for camel by trying to capture both low and high abundant transcripts expressed in different age, breed, and tissue types, we adapted a normalized cDNA library construction approach. In brief, this approach makes use of biotinylated nucleotides during in vitro transcription off of linearized cDNA, which are subsequently removed using streptavadin beads and phenol-chloroform, resulting in 100-200 fold reduction in highly abundant clones.

RNA Extraction
100 mg of each tissue were homogenized in a vial containing 350 ml lyses buffer and ceramic beads using a Magnalyzer instrument (Roche Diagnostics) at 6,500 rpm. Following homogenization, the lysates were further extracted and purified by QIAamp RNA Mini Kit (Qiagen). Purified RNA samples were then analyzed by Nanodrop 1000 and Agilent 2100 Bioanalyzer RNA Nano Kit for the determination of the quantity and integrity of RNA. RNA pooling system was generated for each library. Total of 3 libraries were generated based on animal's age. Total of 2 mg of RNA for each library were used for establishing and sequencing EST libraries.

cDNA Library Construction
First strand synthesis of the cDNA was generated by using a clamped dT primer, which is designed to attach to the poly-A junction of the mRNA. This ensures a more efficient first-strand synthesis and a shorter poly-A region at the 39 end of the insert. cDNA fragments were size selected by running on an agarose gel. Following size selection, cDNA ends are polished and the cDNAs were digested using a rare cutting enzyme and were directionally cloned into pAGEN-1 vector.
Single stranded DNA was made from a portion of the primary library by phagemid production and to ensure elimination of the double stranded DNA contamination from the single stranded DNA preb; the reactions were then digested with DNAse I. A second portion of the primary library was linearized and transcribed into anti-sense RNA with biotinylated dNTPs. oligo dT and primer extension were used to pre-block the poly-A region prior to hybridization. This prevents hybridization of the poly-A clone and the poly-U of the RNA. The anti-sense RNA and single stranded circular DNA were then hybridized and abundant clones were removed using streptavidin. To prevent any empty vectors, A Not1 oligo and Taq polymerase were used to synthesize double stranded DNA from the single stranded normalized library.

EST Data Analysis
Using raw chromatogram files, we obtained corresponding sequence and quality values, which results in base calls and associated scores regarding the accuracy of the call for each EST sequence, or ''read'' in the libraries. Reads thus obtained were analyzed for vector sequence contamination and removal of low quality bases. The remaining high quality reads were checked for the existence of chimeric sequences. Upon fixing for chimeras, reads were masked for known and estimated repeats. These masked regions were not used for similarity calculation required in the clustering step, which groups reads showing high sequence similarity. Contigs and singletons obtained from alignment of reads in each cluster imply putative gene sequences, which were analyzed for homology searches in existing sequence databases and open reading frame (ORF) discovery. The results were annotated using Gene Ontology functional classification for various species that the putative gene sequences have been matched to. An overall summary of the analysis pipeline is depicted in Figure 1.

Quality Calls and Trimming
We started with 23,424 chromatogram files from each library, equaling 70,272 reads overall. Raw chromatogram files were processed using PHRED, which reads trace data, calls bases, and assigns quality values to bases using Fourier transforms [17,18]. We removed low quality regions, vector contamination and potential polyA regions using Lucy2 [19], determining the longest continuous high quality region for each read. If the remaining sequence after the clean-up procedure was less than 100 base pairs (bp), we eliminated the sequence from downstream analysis.

Identification of Chimeric Sequences
During cloning process, two or more DNA fragments can be ligated together yielding a chimeric EST read. Chimeras are usually PCR artifacts thought to occur when a prematurely terminated amplicon reanneals to a foreign DNA strand and is copied to completion in the following PCR cycles. We identified chimeric sequences by searching for internally inserted adaptors or cloning sites in the reads [20]. When such a sequence is found, we calculated the lengths of the subsequences after the chimeric sequence is divided into two or more pieces. If the length of any such subsequence was less than 50 bp, the subsequence was removed from downstream analysis. We performed chimera analysis on the reads trimmed of low quality regions and vector sequences.

Masking Repeats
Repeat regions could result in artificial sequence similarity by being used as ''seeds'' in alignment process [21]. In order to prevent this, we adapted two strategies to mask repeats. The first is the widely used RepeatMasker algorithm that uses a library of known repeats [22]. For this purpose, we used ''Camelus dromedarius'' as the reference organism in the RepeatMasker algorithm. However, since the whole genome of the Arabian camel is not available, the algorithm uses a library that contains the first common ancestor of Camelus dromedarius for which a repeat library has been constructed.
For this very reason that we lack a whole genome for camel, we adopted a second approach to mask repeats [23]. In this ''libraryless'' repeat masking procedure, RBR, we first identified words (a small stretch of DNA sequence) of length at least 16 bp in all EST sequences. Next, a distribution for the frequency of word repeats is obtained. If a word's frequency was above 5 standard deviation of the expected number of repeats for a word, we concluded that such a repeat could not have been caused by resampling of the same mRNA. Therefore such words with high occurrences are regarded to be represented in the EST data set due to DNA repeats. We masked a base pair if it was called to be in a repeat region by either of the two algorithms.

Clustering and Assembly
We clustered reads using The Gene Indices Clustering Tool (TGICL) [24] and for each cluster performed multiple alignment of its members using CAP3 (Huang and Madan, 1999). The consensus sequence thus obtained for each cluster approximates the target mRNA sequence represented by the reads in the cluster. When performing assembly for a group of reads in a given cluster, reads that do not fit well in the overall alignment are dropped out and not included in the final formation of the ''contigs'', the consensus sequence representing a cluster. Such dropped out reads combined with reads that are placed in self-contained clusters are referred to as ''singletons''.

Homology Search and Functional Annotation
We calculated the longest ORF for each contig and singleton sequence using NCBI's standard genetic code to relate homology searches to the availability of ORFs in the input set. We then compared each contig and singleton to NCBI's nucleotide (nt; containing 6,944,581 sequences) and non-redundant protein (nr; containing 6,655,203 sequences) databases using BLASTN (Evalue cut-off 10 230 ) and BLASTX (E-value cut-off 10 25 ) algorithms, respectively [25]. Resulting ''hits'' were analyzed separately for nine species. We also searched for hits in the NIH Mammalian Gene Collection Project (http://mgc.nci.nih.gov) to identify matches to full-length cDNA sequences for the following species: Homo sapiens, Mus musculus, Rattus norvegicus, and Bos Taurus [26]. This search was conducted using BLASTX and results were analyzed to identify coverage of match regions over start and/or stop codons of full-length cDNA sequences. A hit was identified if the alignment between query and hit sequences was over 30 amino acids (aa) with at least 96% identity.
Functional annotation of reads that are mapped to homologous sequences using BLAST is done by identifying the Gene Ontology (GO) categories for the hits [27]. We identified Biological Process, Cellular Component, and Molecular Function categories seen in our matched sequences using NCBI's GENE database mapping to GO terms. Our sequence and functional annotation strategy was NCBI's gi number centric: for every contig or singleton, the sequence's hits' gi numbers were pooled and matched to unique GeneIDs to remove redundancy. Hence match results presented always refer to unique number of genes found in NCBI's GENE database. These GeneIDs were then mapped to GO terms.
In an attempt to perform comparative genomic analysis, we compared genes mapped from our EST database to human, mouse, rat, and bovine using NCBI's HomoloGene database. We found genes shared by all four organisms and analyzed them with Ingenuity Pathway Analysis (IPA) software (IngenuityH Systems, www.ingenuity.com). IPA generates networks based on the connectivity of genes in the data set using Ingenuity Knowledge Base, which is mostly dependent on scientific literature. IPA performs functional analysis by identifying biological functions and/or diseases that are most significant to the data set using right tailed Fisher's exact test. This functional analysis is similarly applied to molecules in a given network generated by IPA for the data set. IPA also projects molecules in the data set to known biological pathways (referred to as ''canonical pathways'') and finds significantly associated canonical pathways by calculating the ratio of molecules in the data set mapped to canonical pathway divided by total number of molecules in the canonical pathway and by using Fisher's exact test to assign significance to observed association.
Results cDNA library created from 11 tissues of three different inbred camels of different ages produced enough material to sequence around 80,000 clones. 23,424 EST sequences were determined from each library by single pass 59 sequencing yielding a total of 70,272 EST sequences, with an average read length of 14476411 (avg. 6 st.dev.) bp and an average high quality (Phrap score .20) per read of 6146283 bp. After trimming vector contamination, low quality bases, and eliminating trimmed sequences with length less than 100 bp, three libraries resulted in 58,842 high-quality EST sequences with an average read length of 7556171 bp and an average high quality base pairs per read of 6706181. These results are summarized in Table 1. The increase of average high quality base pairs per read from 614 (untrimmed) to 670 (trimmed) indicates successful elimination of reads with a high abundance of low quality bases. Moreover, the ratio of average high quality base pairs per read to average read length increased from 42% to 88% after trimming. Tightness in the variation of read length and number of high quality bases per read following trimming indicates a high quality set of reads for downstream analysis (see Figure 2).
We identified 1,241 chimeric sequences, accounting for about 2% of high-quality reads. After splitting chimeric sequences and removing remaining subsequences that are less than 50 bp, we were left with 59,534 reads, an increase of about 1% in number of reads. RepeatMasker and RBR found repeats in about 10% and 20% of the reads, respectively, though RepeatMasker masked ,1.5610 6 bp, roughly 1.5 times more bp than masked by RBR. When we combined regions masked by either method, about 30% of all the reads and roughly ,2.5610 6 bp were masked suggesting both sequences and regions masked by the two methods were almost exclusive. These results imply that our approach to repeat finding identifies regions that could have potentially been missed if either of the programs were solely used; an approach adapted by most large-scale sequencing projects. High-quality EST sequences masked for repeats were assembled by TGICL resulting in 8,319 contigs and 15,283 singletons yielding 23,602 putative gene sequences. In Figure 3, we show an instance of a cluster and assembled consensus sequence containing thirty reads. In forming this contig, we considered reads with overlaps of at least 98% identity, at least 40bp, and at most 20bp overlap distance of sequence end. Overall, average number of reads in a cluster was 5.2 (see Table 1). Average sequence length for the contigs and singletons were 1,2476459 bp and 6966216 bp (avg. 6 st.dev.), respectively (see Table 1). In Figure 4a, we show sequence length distribution where sequences of length longer than 2,500 are not shown for display purposes. About 2.3% of contig and no singleton sequences have length longer than 2,500 bp. This difference in sequence length is expected as contigs are supported by multiple reads. On the other hand, sequence lengths are long enough to suggest identification of full or partial coding sequences. In Figure 4b, we show length distribution of the longest ORF found in contigs and singletons. As expected, contigs yield longer ORFs (avg. length: 6736415 bp) with ,80% of contigs having an ORF longer than 300 bp suggesting successful characterization of gene coding regions found in camel DNA. ,55% of singletons have an ORF longer than 300 bp (avg. ORF length: 3906220 bp; see Table 1). This difference is expected as the starting average length for contigs are higher than that of singletons. About 20% of the contigs have an ORF longer than 1,000 bp; 2% of which are longer than 1,800 bp, not shown in Figure 4b for display purposes.
In order to identify likely camel genes through sequence homology, BLAST analysis was performed on 23,602 contigs and singletons against NCBI's nucleotide and non-redundant protein databases. Out of 8,319 contigs 7,490 (90%) sequences got a hit and 829 (10%) sequences got no hit. Similarly, for 15,283 singletons, there were 11,480 (75%) and 3,803 (25%) sequences with a hit and no hit, respectively. Overall, 18,970 (80%) sequences were hit with one or more matches and there were 4,632 (20%) sequences that showed no significant similarity to any of the sequences in the databases. 4,632 sequences with ''no hits'' may represent novel camel specific genes or genes that evolve relatively quickly.
We also analyzed length and ORF distributions in contigs and singletons based on their ''hit'' status. These results are shown in Figure 4c,d. The average length 6 st. dev. for contig sequences that got hit and no hit were 1,2686465 bp and 1,0576125 bp, respectively. Similarly, average length 6 st. dev. values were 7136202 bp and 6466247 bp for singleton sequences that got hit and no hit, respectively. These results show about 16.6% and 9.4% decrease in average lengths between sequences with a hit and no hit for contig and singleton sequences, respectively. As the length of reads that got no hit are not dramatically smaller, we believe this indicates potential discovery of novel coding regions in camel. On the other hand, the average length 6 st. dev. of longest ORFs found in contig sequences that got hit and no hit were 7146416 bp and 3066125 bp, respectively. The difference in average ORF length between contigs that got hit and no hit is about 57%, significantly larger than the average length difference for the two groups. We observe a similar behavior for singletons where average length 6 st. dev. values were 4376227 bp and 2486112 bp (,43% difference) for longest ORF lengths found in singleton sequences that got hit and no hit, respectively. Although reads that got no hit show small ORF lengths, there were 362 contig and 1,409 singleton ''no hit'' sequences with an ORF longer than 300 bp suggesting ,40% of potentially novel genes show strong coding capacity.
We further analyzed BLAST results separately for nine species. These results are summarized in Table 2. For example, 85% of all contigs found a hit in human and the average ORF length of these contigs were 740 bp. 84% of the contigs that found a match in human contained an ORF longer than 300 bp. Overall, we generally have higher 70s to 80% of contigs and 50-60% of singletons with a hit in the nine species when analyzed separately. Moreover, the average ORF lengths were 700-800 bp for contigs and 400-500 bp for singletons with about 90% of contigs and 60-70% of singletons having an ORF length greater than 300 bp. These results suggest that when reads with a hit are analyzed in a specific species, we most likely find homologous gene sequences with coding capacity. We also found genes that got most hits by camel sequences for the nine species. A complete list for each analyzed species can be found in CAGBASE and lists of top fifty genes for human, mouse, rat, and bovine are shown in Table 3, Table 4, Table 5 and Table 6. We also compared our EST sequences to all 674 available Camelus Dromedarius mRNA sequences in NCBI GenBank. Out of these 674 sequences only 16 were not matched significantly. The results for all 674 sequences along with top 10 hits from our camel EST database can be found in CAGBASE. Table 7 summarizes functional annotation results for putative camel gene sequences. Overall, regarding species for which GO annotation is available, around 70% of genes that got a hit contain a GO term, corresponding to about 90% of all the GO terms found in the organism. For example, when all camel sequences are considered in combination, compared to Mus musculus, 14,390 GeneIDs (or 65%) out of 22,212 GeneIDs that are hit contain a GO term amounting to 6,098 distinct GO terms, which is 95% of a total of 6,452 distinct GO terms found in the species. These results suggest that identified camel sequences capture most of functional characteristics described in other species. In Figure 5, we show top highly abundant Biological Process GO Terms found in camel among the genes mapped to Homo sapiens. When most frequently matched functional GO categories are studied, we see a broad range of functional terms but interestingly some that could relate to properties more commonly attributable to camel such as ''oxidation reduction'' and ''sensory perception of smell''. A complete list of GO results for all four species can be found in online CAGBASE web portal.
We compared camel ESTs to full-length cDNA sequences of Homo sapiens, Mus musculus, Rattus norvegicus, and Bos Taurus containing 28,133; 23,120; 5,341; and 9,188 sequences, respectively. Summarized in Table 8, we identified 2,089; 1,543; 884; and 1,646 camel ESTs that show high similarity (96% identity over at least 30 aa) to full-length cDNA sequences of aforementioned organisms, respectively. In total, we have 2637 camel ESTs that matched to a full-length cDNA in any of the four organisms, corresponding to about 11% of all camel sequences. A complete  , sequence length distribution for contigs and singletons with hits and no hits (c), and distribution of longest ORF lengths found in contigs and singletons with hits and no hits (d). Average 6 standard dev. values of sequence and ORF lengths are overlaid on corresponding graphs. Sequence lengths up to 2,500 and ORF lengths up to 1,800 bp are shown for display purposes. 2.3% of contig and no singleton sequences have length longer than 2,500 bp (a), 2% of contig and no singleton sequences have an ORF longer than 1,800 bp (b), 1.8% of contigs with a hit and no other sequences in the remaining three groups have length longer than 2,500 bp (c), and 1.6% of contigs with a hit and no other sequences in the remaining three groups have an ORF longer than 1,800 bp (d). doi:10.1371/journal.pone.0010720.g004       list of full length cDNA analysis is hosted at CAGBASE web portal. We found that in 35-50% of matching contigs, the alignment was within 5 aa of the start codon, suggesting that nearcomplete coding sequences have been obtained for these genes. Lack of a similar database for camel unfortunately prevents these results to be processed to better estimate sequences representing full length cDNAs and total number of genes in camel.
EST sequences have mapped to 22,980; 22,212; 14,991; and 14,766 genes in human, mouse, rat, and bovine, respectively. We first matched these genes to HomoloGene Group IDs, resulting in 16,405; 16,638; 11,834; and 11,773 matches in human, mouse, rat, and bovine, respectively. When these lists were compared, we found 8,405 genes shared by all four organisms constituting a large portion of genes identified in each organism. Figure 6 shows Venn diagram depicting comparison of genes found (and then matched in HomoloGene) in four organisms. Regions not shown in the Venn diagram are genes shared by human and bovine only (403 genes) and mouse and rat only (536 genes). Complete lists of genes for each region of Venn diagram can be found in CAGBASE. In Table 9, we list top 50 genes shared by all four organisms based on their collective frequency of occurrence. This list varies from the genes shown in Table 3, Table 4, Table 5, and Table 6 as most of the genes shown in Table 3, Table 4, Table 5, and Table 6 are not mapped by HomoloGene database. IPA analysis of shared genes revealed a network with 35 molecules where most significantly associated biological functions were ''hair and skin development and function'' and ''renal and urological system development'' (p,10 25 ), This network along with sample molecules participating in aforementioned functions are depicted in Figure 7. In Figure 8, we show part of the top canonical pathway, NRF-2 mediated oxidative stress response (p,10 210 ), which is most significantly associated with shared genes based on IPA.

Discussion
We have performed large scale EST sequencing of Camelus Dromedarius as the first phase of whole genome sequencing of the organism. This work is the first large scale sequencing for camel and has provided over 4,500 potentially novel or fast evolving camel gene sequences that do not carry any homology to other available genomes. Our goal in this work was to generate a comprehensive list of coding sequences found in camel genome irrespective of any specific phenotype. For this purpose, we pooled samples from different tissue, age, and breed. We believe inclusion of a more diverse sampling would yield a better coverage of camel transcriptome but current sample set coupled with our cDNA library normalization strategy should suffice for an initial phase study identifying camel genes. Our results based on read quality, homology, ORF, GO based functional analysis, full length cDNA, and comparative genomic analysis suggest that we have successfully found genes in camel matched to known coding sequences in other organisms and novel gene sequences with a protein coding capacity.
When putative camel gene sequences were analyzed for ORF length based on their hit status, we observed significantly shorter ORFs in sequences with no hit. These results suggest that ORF length, not the sequence length, is a better indicator of finding transcripts with protein coding capacity and subsequently getting a hit in sequence databases. On the other hand, still more than one third of the sequences that got no hit contained an ORF greater than 300 bp and we believe that sequences that got no hit with a relatively long ORF represent novel genes with protein coding capacity. Our results also showed a higher no hit percentage in singleton sequences. We believe that the reason for this is possibly due to the fact that singletons represent rare genes in camel genome that is not well described in other organisms.
It is common to see up to 20% redundancy in large scale EST projects [28]. Instead of assessing such a redundancy in our data set, as no reference sequence set for camel exists, we took advantage of this behavior by assessing genes that got most hits by camel sequences for nine analyzed sequences (see Table 3, Table 4,  Table 5, Table 6 and CAGBASE). There is considerable consensus in the lists, which are dominated by immunoglobulin heavy-chain genes (known to be naturally occurring in camel), genes that belong to major histocompatibility classes, and zinc finger proteins. Regarding the IgG isotypes in camels as compared to other species, many authors reported that the numbers of IgG isotypes in different mammalian species vary considerably, depending on the number of functional genes, ranging from one in rabbit [29], three in cattle [30], four in human [31], mouse [32] and rat [33], five in pig [34] and six in horse [35]. It was speculated that increase in the number of IgG isotypes might be accompanied by an increase in functional diversity [36]. Sequences of distinct cDNAs homologous to zinc finger proteins indicate that alternative splicing events adjoin either coding or non coding exons to the finger preceding box sequences [37]. Genes of histocompatibility classes are found in all vertebrates, though the gene composition and genomic arrangement may vary widely due to gene duplication. It was reported that Xanthine oxidoreductase (XOR) activity of camels was markedly lower than that of human, bovine and goat enzymes obtained under the same conditions [38]. This work suggested that the molybdo-form of camel enzyme is totally under desulpho inactive form. It is possible that camel neonates are equipped with an enzymic system that reactivates XOR in their gut and consequently generates antibacterial reactive oxygen species. There are sporadic instances of uniden-tified genes that have been highly matched by camel sequences, which warrant further biological investigation. GO analysis revealed expected categories such as ''transcription'', ''translation'', and ''cell cycle'' (see Figure 5), however, particularly interesting in camel, one of the overrepresented categories was ''oxidation reduction''. Psychosocial stress generally increases oxidative stress promoting a pro-inflammatory environment potentially through facilitation of NF-kappa B between stress and oxidative cellular activation [39,40]. Oxidative phosphorylation and generation of reactive oxidative species diminishes the negative effect of stress [41]. In this context, surfacing of oxidation reduction process in camel among highly abundant biological processes found in our camel EST database, agrees with camel being a low tempered species with relatively slow movements [4]. Moreover, this observation agrees with recent findings that camel milk alleviates oxidative stress [42] and common infections in camel is associated with a state of oxidative process [43]. There were 514 camel ESTs matched to human genes that were associated with oxidation reduction category including cytochrome c oxidase related genes (COX family), ubiquinolcytochrome c reductase binding protein (UQCRB), and genes from proxiredoxin (PRDXx) family. Another GO category revealed to be significant that stand outs from generally expected processes is ''sensory perception of smell'', which is known to be extremely good in camel [9]. There were 195 camel ESTs matched to human genes that were associated with perception of smell including a broad coverage of olfactory receptors, dopamine transporter (SLC6A3) and receptor (DRD2) and glutamate receptors (GRMx).
The list of most frequently matched 50 out of 8,405 genes found in camel and shared by human, mouse, rat, and bovine was dominated by globulin inhibitors (ITIH family) and zinc finger proteins (ZNF family; see Table 9). ZNFx play a role in diverse functions including lipid binding [44], which may be relevant in camel as camel uses extensive reservoirs of fat to adapt to its environment. The top network generated by IPA using shared genes was associated with hair and skin development, which is extremely important in camel as its thick coat and fur are among the most prioritized adaptations in camel. This network, depicted in Figure 7, holds damage specific DNA binding protein 1 (DDB1) and cyclin dependent kinase inhibitor 1 (CDKN1B) as its hubs. DDB1, CDKN1B, along with cullin 4A (CUL4A) and Vpr binding protein (VPRB) participate in hair and skin development by functioning in arrest of epithelial cell lines in G1 phase [45] and G2 phase [46]. The other most significantly associated function with this network was renal and urological system development, which are important functionalities in camel as its kidneys and intestines efficiently retain water as part of camel's adaptation [9]. This function was facilitated by control of aforementioned four genes in size of kidney cell lines [47], arrest of kidney cell lines in G1 phase [45] and G2 phase [46], and arrest of mesangial cells in G1 phase [48]. In Figure 7, associated functions and genes (in light green) are shown along with complete interaction of all 35 molecules participating in this network. In IPA analysis of 8,405 genes, most significantly associated pathway with the data set came out to be ''NRF-2 mediated oxidative stress response pathway'', validating the results of independent GO analysis (see Figure 8). This pathway, which is mostly regulated by mitogen activated protein kinases -all found in the camel ESTs-is anchored by NRF-2. NRF-2 is a key regulator of antioxidant response [49], deficiency of which results in increased quantities of reactive oxygen species [50]. In Figure 8, we show a part of this pathway where genes in the input data set are indicated in red.
In order to share our findings with the scientific community, we developed a web portal hosting the EST database and analysis results at http://camel.kacst.edu.sa. The web portal is made public domain with possibility to add sequences to the existing Figure 8. Canonical Pathway: Most significantly associated pathway (NRF-2 mediated oxidative stress response pathway) by the data set of 8,405 genes found in camel ESTs shared by human, mouse, rat, and bovine using IPA. Molecules that exist in the data set are shown in red. doi:10.1371/journal.pone.0010720.g008 database. A number of tools are also made available such as the ability to BLAST against our EST sequences. The database has been named CAGBASE, after Camel Genome Database. Through this web portal, we present a relational database showing all camel sequences along with corresponding genes and GO terms found for different species. The web portal accepts queries with GeneID, gene symbol, gene name, or GO ID to find camel sequences associated with any one of the query terms. We hope that this web portal provides a home base for genetic studies regarding Camelus Dromedarius.