Sadd et al. BMC Genomics 2010, 11:110 http://www.biomedcentral.com/1471-2164/11/110 (Aesch, Basel), and were kept under red-light illumina- tion at 27 ± 1°C with sugar water (ApiInvert®) and pol- len provided ad libitum. Worker bees used for the EST library originated from four colonies (two workers each) that had been confirmed as Bombus terrestris by visual checks and based on a species-specific length poly- morphism marker in the region between mitochondrial CO1 and CO2 genes. The colonies were also confirmed to be free of common parasites by microscopic investi- gation of faeces. Worker bees were age controlled, and one bee from each colony was 7 and one bee 14 days post adult eclosion at the time of sacrifice. Bees were sacrificed by snap freezing in liquid nitrogen, and were subsequently stored at -80°C. Tissue from both the thorax and abdomen was used as material for the EST library. RNA isolation and cDNA library construction Total RNA was extracted with Solution D [58]. cDNA was synthesized using the SMART approach (Mint-Uni- versal cDNA synthesis kit, Evrogen, Russia), subse- quently normalized using duplex-specific nuclease (Trimmer kit, Evrogen, Russia) according to manufac- turer’s instructions, and directionally cloned into pAL32 plasmid vector (Evrogen, Russia). Plasmids were trans- ferred via electroporation into E. coli DH10B (Invitro- gen, U.S.A.). EST sequencing, quality control and assembly Plasmids from the normalized libraries were 5’ end sequenced using the pALforward primer (5’- CTCGGGAAGCGCGCCATT-3’) and Big Dye Termina- tor chemistry (ABI). Collected reads from 3730XL capil- lary sequencers (ABI) were base-called using Phred http://www.phrap.org and subsequently quality and vec- tor clipped using Lucy http://compbio.dfci.harvard.edu/ tgi/software with standard parameters. For clustering and assembly, the TGI tools developed at TIGR http:// compbio.dfci.harvard.edu/tgi/software were used. Annotation and comparisons Where explicit packages and scripts are not mentioned, results were obtained by a mixture of custom Perl scripts implementing BioPerl modules [59] and custom R functions in R2.9.2 for Mac [60]. Once clusters were created, the resulting unique sequences were used for various comparisons. BLASTX searches with a cut-off of 1e-05 were run against each of the nr database in Genbank (version dated 19th October 2009), and databases comprising of all pre- dicted proteins in the genomes of A. mellifera, N. vitri- pennis, D. melanogaster and T. castaneum. Redundancy in the unique sequences was determined by comparison with A. mellifera. This was calculated as Page 7 of 9 the number of unique A. mellifera proteins that were “best-hits” in the BLASTX search divided by the number of unique B. terrestris sequences with hits. GC-content was calculated using a custom Perl script that amalga- mated all unique sequences and calculated the propor- tion of G and C bases in all unambiguous bases. Annotation of GO terms [61] to the current B. terres- tris EST derived unique sequences were carried out using the Perl based annot8r [38]. A database containing SwissProt and Trembl (Uniprot Knowledgebase Release 15.9) with manually curated GO terms was used. Terms that had been derived electronic annotation were omitted to improve reliability of the annotations. BLASTP against this database was carried out with a cut-off of 1e-05. In addition, predicted proteins from the entire genomes of A. mellifera [32] and T. castaneum [39] were compared with the same database for compar- isons of the distribution of GO terms. For statistical comparisons, GO terms were matched to higher-level GO Slim categories. Pairwise comparisons between B. terrestris and both A. mellifera and T. castaneum were then carried out for the distribution of sequences in each GO Slim category using Fisher’s exact tests in R2.92 for Mac [60]. The EST unique sequences were compared with the A. mellifera genome [32] using BLASTN, cut-off 1e-08. A custom Perl script was then used to extract positional information of all the “best-hits”. Linkage group and position on the linkage group was then compared with the position of predicted protein coding regions within the genome to obtain a relative position of the “hits” (inside a protein-coding region; spanning a protein-cod- ing region boundary; outside a protein coding region). For those outside a predicted protein-coding region, the distance to the closest was calculated. For those B. ter- restris sequences with hits greater than 4’375 bases (the largest A. mellifera UTR, obtained from the UCSC gen- ome browser [42]) away from a protein coding region, a custom Perl script was used to interface with ORF Fin- der http://www.ncbi.nlm.nih.gov/gorf to obtain ORF information. Of these sequences, the 20 with the largest ORFs were investigated manually. BLASTP searches were carried out against the nr database, and if the “best-hit” was to another insect protein (cut-off = 1e- 20), the “best-hit” was taken and blasted (TBLASTN) against the A. mellifera genome. It was confirmed that the “hit” location was the same as for the EST derived unique sequence and the Blast alignments were manu- ally examined for potential features inhibiting protein production. SSR and SNP marker identification To identify potential genetic markers we used the Perl scripts MISA [50] and SNPidentifier [57] to locate