UBC Faculty Research and Publications

Finding the missing honey bee genes: lessons learned from a genome upgrade Elsik, Christine G; Worley, Kim C; Bennett, Anna K; Beye, Martin; Camara, Francisco; Childers, Christopher P; de Graaf, Dirk C; Debyser, Griet; Deng, Jixin; Devreese, Bart; Elhaik, Eran; Evans, Jay D; Foster, Leonard J; Graur, Dan; Guigo, Roderic; Hoff, Katharina J; Holder, Michael E; Hudson, Matthew E; Hunt, Greg J; Jiang, Huaiyang; Joshi, Vandita; Khetani, Radhika S; Kosarev, Peter; Kovar, Christie L; Ma, Jian; Maleszka, Ryszard; Moritz, Robin F A; Munoz-Torres, Monica C; Murphy, Terence D; Muzny, Donna M; Newsham, Irene F; Reese, Justin T; Robertson, Hugh M; Robinson, Gene E; Rueppell, Olav; Solovyev, Victor; Stanke, Mario; Stolle, Eckart; Tsuruda, Jennifer M; Vaerenbergh, Matthias V; Waterhouse, Robert M; Weaver, Daniel B; Whitfield, Charles W; Wu, Yuanqing; Zdobnov, Evgeny M; Zhang, Lan; Zhu, Dianhui; Gibbs, Richard A Jan 30, 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


52383-12864_2013_Article_5733.pdf [ 1007.67kB ]
JSON: 52383-1.0307140.json
JSON-LD: 52383-1.0307140-ld.json
RDF/XML (Pretty): 52383-1.0307140-rdf.xml
RDF/JSON: 52383-1.0307140-rdf.json
Turtle: 52383-1.0307140-turtle.txt
N-Triples: 52383-1.0307140-rdf-ntriples.txt
Original Record: 52383-1.0307140-source.json
Full Text

Full Text

RESEARCH ARTICLE Open AccessFinding the missing honey bee genes: lessonslearned from a genome upgradeChristine G Elsik1,2*†, Kim C Worley3*†, Anna K Bennett2†, Martin Beye4, Francisco Camara5, Christopher P Childers2,6,Dirk C de Graaf7, Griet Debyser8, Jixin Deng3, Bart Devreese8, Eran Elhaik9, Jay D Evans10, Leonard J Foster11,Dan Graur12, Roderic Guigo5, HGSC production teams3, Katharina Jasmin Hoff13, Michael E Holder3,Matthew E Hudson14, Greg J Hunt15, Huaiyang Jiang16, Vandita Joshi3, Radhika S Khetani17, Peter Kosarev18,Christie L Kovar3, Jian Ma19, Ryszard Maleszka20, Robin F A Moritz21, Monica C Munoz-Torres2,22, Terence D Murphy23,Donna M Muzny3, Irene F Newsham3, Justin T Reese2,6, Hugh M Robertson24, Gene E Robinson25, Olav Rueppell26,Victor Solovyev27, Mario Stanke13, Eckart Stolle21, Jennifer M Tsuruda28, Matthias Van Vaerenbergh7,Robert M Waterhouse29, Daniel B Weaver30, Charles W Whitfield31, Yuanqing Wu3, Evgeny M Zdobnov29, Lan Zhang3,Dianhui Zhu3, Richard A Gibbs3, on behalf of Honey Bee Genome Sequencing ConsortiumAbstractBackground: The first generation of genome sequence assemblies and annotations have had a significant impactupon our understanding of the biology of the sequenced species, the phylogenetic relationships among species,the study of populations within and across species, and have informed the biology of humans. As only a fewMetazoan genomes are approaching finished quality (human, mouse, fly and worm), there is room forimprovement of most genome assemblies. The honey bee (Apis mellifera) genome, published in 2006, was notedfor its bimodal GC content distribution that affected the quality of the assembly in some regions and for fewergenes in the initial gene set (OGSv1.0) compared to what would be expected based on other sequenced insectgenomes.Results: Here, we report an improved honey bee genome assembly (Amel_4.5) with a new gene annotation set(OGSv3.2), and show that the honey bee genome contains a number of genes similar to that of other insectgenomes, contrary to what was suggested in OGSv1.0. The new genome assembly is more contiguous andcomplete and the new gene set includes ~5000 more protein-coding genes, 50% more than previously reported.About 1/6 of the additional genes were due to improvements to the assembly, and the remaining were inferredbased on new RNAseq and protein data.Conclusions: Lessons learned from this genome upgrade have important implications for future genomesequencing projects. Furthermore, the improvements significantly enhance genomic resources for the honey bee,a key model for social behavior and essential to global ecology through pollination.Keywords: Apis mellifera, GC content, Gene annotation, Gene prediction, Genome assembly, Genome improvement,Genome sequencing, Repetitive DNA, Transcriptome* Correspondence: elsikc@missouri.edu; kworley@bcm.edu†Equal contributors1Division of Animal Sciences, Division of Plant Sciences, and MU InformaticsInstitute, University of Missouri, Columbia, MO 65211, USA3Human Genome Sequencing Center, Department of Molecular and HumanGenetics, Baylor College of Medicine, MS BCM226, One Baylor Plaza, Houston,TX 77030, USAFull list of author information is available at the end of the article© 2014 Elsik et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the CreativeCommons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, andreproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedicationwaiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwisestated.Elsik et al. BMC Genomics 2014, 15:86http://www.biomedcentral.com/1471-2164/15/86BackgroundHoney bee, Apis mellifera, provides a key model for un-derstanding the diverse and widespread group of socialinsects. Thanks, in part, to resources generated duringthe initial honey bee genome sequencing effort [1],honey bees are now widely used for elucidating the mo-lecular basis of behavior [2], development [3], diseasetransmission [4], epigenomics [5] and gene regulation[6]. Honey bees are also a vital pollinator of many of theworld’s crops [7,8], and genomic tools are being used toaddress recent serious and unexplained declines in somehoneybee populations [9,10]. In light of the currentinterest in honey bee genetics, improved genomic toolsfor this species are required.Because fewer than 11,000 genes were predicted, oneof the major questions coming out of the initial honeybee genome sequencing effort was, “Where are the miss-ing genes?” Were there issues with the genome sequencethat meant that the genes were filled with gaps andtherefore not annotated? Or were the genes located inislands of repetitive sequence and therefore not identi-fied? Were the gene prediction methods poorly adaptedto a genome with a bimodal distribution in GC content?Or was the tuning of gene prediction algorithms to anaverage GC content failing to represent genes in regionsof more extreme GC content? Was there insufficientmRNA evidence? Were the genes different enough fromknown genes that the prediction tools failed to find thegenes due to lack of protein evidence? Or might thehoney bee have thousands fewer genes than the fewinsects with sequenced genomes at the time?All early genome sequencing projects used Sanger dataand either a BAC-based hierarchical-sequencing model[11] or a whole-genome-shotgun model [12]. In eithercase, a completed draft genome has contiguous sequencepieces (contigs) spaced by paired clone end sequences,forming scaffolds. Missing sequences between consecu-tive contigs are represented by segments of ambiguousbases, denoted as Ns, the lengths of estimated gap sizes.Few Metazoan genomes have been improved to finishedquality by filling in the missing sequence [13-16]. Thegenome scaffolds, with islands of sequence informationand lack of information (in a draft genome), are anno-tated with gene features using gene prediction methods.When gene signals such as splice sites, start codons andtermination codons are missing from the assembly, acomputational gene prediction tool may miss an exon oreven an entire gene. Thus, gaps within genic regions ofan assembly can hinder annotation and lead to an in-complete gene list.The first A. mellifera genome sequencing project [1] re-vealed genome characteristics with potential missing as-sembly information that could impact the gene list. Thegenome assembly had the lowest mean GC content(percent of G +C nucleotides) and the most heterogeneousGC content of any sequenced metazoan genome at thattime, with GC content ranging from 10% to 70% across thepublished genome assembly (Amel_4.0) [1]. Analyses of ini-tial assemblies showed that regions with low GC contentwere under-represented in libraries, so additional shotgunlibraries were generated after fractionating DNA with CsCl-bisbenzimide density gradient centrifugation, and appro-ximately 800,000 reads of <30% GC content were added tothe data that lead to Amel_4.0 [1].In addition to potentially missing assembly informa-tion, other factors including limited EST data and thelarge evolutionary distance between honey bee and otherinsects with sequenced genomes were thought to havecontributed to difficulty in gene prediction. The 78,001A. mellifera ESTs that were available at the time pro-vided only 9,408 unique consensus sequence alignmentsto the genome [17]. The closest organisms with availableprotein sets were the Dipterans, Drosophila melanogasterand Anopheles gambiae, with an estimated divergencetime from honey bee (Hymenopteran) of 300 millionyears.Honey bee researchers suspected that the first OfficialGene Set (OGSv1.0), generated by creating a consensusgene set with GLEAN [1,17], was incomplete, as it com-prised only 10,157 genes, far fewer than expected basedon the estimate of ~13,600 genes in D. melanogaster atthe time. Recently sequenced Hymenoptera genomes(parasitoid wasp and seven ant species), predicted tocontain ~17,000-18,500 protein-coding genes [18-24],further support the suspicion that a large number of A.mellifera genes had not yet been detected. Whole ge-nome tiling arrays provided experimental evidence formissing genes, with the detection of signals “intergenic”to OGSv1.0 genes [1]. Finally, evaluation of the genome-wide distribution of OGSv1.0 genes revealed another po-tential factor in computational gene prediction, whichrelied on algorithms optimized for gene and genomecharacteristics known at that time. Unlike any previouslysequenced metazoan genome, OGSv1.0 genes occurredin gene regions that were more GC-poor (29% GC) thanthe mean genome GC content (33%), and were locatedin regions with GC content as low as 11% [1].Genome sequencing and assembly methods have chan-ged with the advent of next-generation sequencing tech-nologies. Newer technologies typically produce shorterreads at greatly reduced cost and allow the completion ofprojects with much deeper sequence coverage using manymore reads than a Sanger-based project. While contiguoussequence may suffer from these shorter reads that cannotspan repetitive sequences, scaffolding can benefit from theincreased mate-pair information. Furthermore, different se-quencing technologies have different biases related to nu-cleotide composition, so combining technologies providesElsik et al. BMC Genomics 2014, 15:86 Page 2 of 29http://www.biomedcentral.com/1471-2164/15/86the potential for one technology to fill in gaps produced byanother.A bigger impact of the reduced sequencing cost ofsecond-generation sequencing methods is the ability togenerate much more transcript sequence than ever be-fore. These transcript data are very useful as evidencesupporting gene model prediction. As will be discussedbelow, the alignment of new transcriptome data pro-vided evidence for most of the previously un-annotatedgenes in the honey bee genome.This paper presents a de novo annotation of the A.mellifera genome based on an upgraded genome assem-bly and new transcript data. We generated additionalgenome sequence data using ABI SOLiD and Roche 454paired-end sequencing technologies to superscaffold andfill gaps in the Amel_4.0 assembly using the Atlas-linksoftware, creating the improved assembly (Amel_4.5).We also deep-sequenced the transcriptome of seven A.mellifera tissues and sequenced the genomes of twoclosely related species, the dwarf honey bee (A. florea)and buff-tailed bumble bee (Bombus terrestris) (both ofwhich will be published elsewhere).To avoid perpetuating errors in OGSv1.0 and avert dif-ficulties in tracking genes between assemblies (e.g. re-ported in [25]), we chose not to track and update theprevious annotation on the Amel_4.5 assembly. Insteadwe decided to annotate the Amel_4.5 assembly de novoand generate a new official gene set (OGSv3.2) and anew list of repetitive elements based on the most currentevidence and methodologies. To learn which factorswere most important in identifying previously unknowngenes, and to inform potential strategies for future ge-nome projects, we compared both gene sets and charac-terized genes that were missing from OGSv1.0.ResultsImprovements to the genome assemblyAdditional sequence coverage generated using the ABISOLiD [26] technology (~20x) and Roche 454 [27] tech-nology (~4x) from small insert fragments (~2 k) was in-corporated into Amel_4.5. The sequence data aresummarized in Table 1 and Table S1 in Additional file 1and are available from the NCBI Sequence Read Archive(SRA). We compared assembly statistics for Amel_4.5with the previous assembly, Amel_4.0, which is de-scribed in [1] and is the assembly used by the researchcommunity since 2006. The genome assembly improve-ments increased the contig N50 from 40 kb to 46 kband the scaffold N50 from 359 kb to 997 kb, with anadditional 5.5% of the sequence anchored to linkagegroups (Table 2).The completeness (extent of coverage of the genome) ofthe new Amel_4.5 and the earlier, Amel_4.0, assemblieswere measured by BLASTN comparison to available denovo assembled 454 contigs from seven libraries. Amel_4.5showed slightly more complete coverage than Amel_4.0(88.7% vs. 88.2%; Table 3).As hoped, the new sequence data had a large impact onthe regions of the genome that were of low starting quality.While the new sequence did not specifically target the GC-poor regions of the genome that were underrepresented inthe initial Sanger libraries, the sequence reads were moreevenly distributed across the regions of different GC con-tents (Figure 1). A segmentation analysis of Amel_4.5 toidentify GC compositional domains that are homogenousin GC content within domains, but different in GC contentacross domain boundaries, showed that the genomic pro-portion of compositional domains low in GC content in-creased in Amel_4.5 (Figure 2). As a consequence, manynewly discovered genes were annotated in GC poor regions(Figure 2).Generating and assessing a New official gene SetTesting combined gene set approachesWe tested two approaches, MAKER2 [28] and GLEAN[17], for combining outputs from multiple gene predictiontools to create an Official Gene Set (OGS). Each approachhas advantages and disadvantages. An advantage ofMAKER2 is that given sufficient transcript data, it can pre-dict multiple alternative splice forms per locus, whileGLEAN computes only one gene model per locus.MAKER2 provides both a conservative biological sequence(transcript and homolog) alignment-based set of genemodels that best agree with transcript and protein homologTable 1 Additional sequence data for the improvedhoney bee genomeMate-pairdistanceRoche 454fragmentRoche 4542.75 kbSOLiDNumber of reads 2.9 M 5.24 M 90 MRead length 290 bp 56 bp 50 bpPairs No Yes YesSequence coverage 3.6× 1.3× 20×Table 2 Assembly statistics for the improved honeybee genomeAssembly Number N50(kb)aBases + Gaps(kb)Bases(kb)Amel_4.5 Anchored 340 1,209 203,000 200,000Scaffolds 5,644 997 250,271 229,734Contigs 16,501 46 229,734 229,734Amel_4.0 Anchored 626 621/135 217,195 183,323Scaffolds 10,742 359 315,719 231,029Contigs 18,944 40 231,029 231,029aFor Amel_4.0 Anchored scaffolds, N50 was calculated separately for 320oriented and 306 non-oriented scaffolds.Elsik et al. BMC Genomics 2014, 15:86 Page 3 of 29http://www.biomedcentral.com/1471-2164/15/86alignments, and a set of ab initio gene predictions that donot overlap the biological sequence alignment-based set.The final set of MAKER2 biological sequence alignment-based gene set contained 12,268 genes, encoding 12,575transcripts, with an additional set of 31,047 ab initio genepredictions that did not have support from transcript orprotein alignments. Without filtering the ab initio set forlikely false positives, the choices of final gene sets generatedwith MAKER2 would have been 1) a conservative set ofonly 12,268 genes, still much lower than gene numbersfound in other sequenced insect genomes or 2) a set of43,315 genes that included ab initio gene models. Otherhymenoptera genome projects using MAKER have selectedab initio gene models with matches to InterPro domains[29] to include in the final predicted gene set [20-22]. Thisapproach would allow us to add 1,215 MAKER2 ab initiogenes to the final predicted gene set, for a total of 13,303genes. However, since only 8,602 of the 12,268 MAKER2biological sequence alignment-based genes had InterPromatches, we were concerned that the InterPro approach toboost ab initio genes to the final gene set would miss somerapidly diverging genes.An important difference between GLEAN and MAKER2is the way ab initio predictions without biological sequenceevidence are handled. GLEAN collects evidence for genes (abinitio gene predictions and biological sequence alignment-based), generates maximum likelihood estimates of accuracyand error rates for these signals for each gene predic-tion set, constructs consensus gene models that maximizethe scores for the signal sites in each gene model, and la-bels each model with a probabilistic confidence scorereflecting the underlying support for that gene model [17].In addition to identifying consensus models supported bytranscript and protein alignments, GLEAN identifies abinitio predictions with high computed probabilities, filter-ing out likely false positive ab initio predictions. This ap-proach would not require filtering ab initio predictionsusing an InterPro search, and would potentially includehigh confidence ab initio predictions that are true, rapidlydiverging genes without significant InterPro domainmatches.Since GLEAN weighs different sources of gene predic-tion evidence based on agreement among sources,Table 3 Assembly comparison to 454 transcriptome dataTissue TotalassembledtranscriptsAmel_4.5 Amel_4.0Number Percent Number PercentAbdomen 14,614 13,980 95.6% 13,987 95.7%Brain + ovary 27,412 26,341 96.0% 26,342 96.0%Embryo 19,616 18,565 94.6% 18,565 94.6%Larvae 18,050 9,061 50.1% 9,041 50.0%Mixed antennae 14,891 13,868 93.1% 13,865 93.1%Ovary 28,451 27,500 96.6% 26,929 94.6%Testes 10,557 9,234 87.4% 9,060 85.8%Total 133,591 118,549 88.7% 117,789 88.2%Note: BLAT alignments of the assembled transcripts to the genome assembliesusing default parameters and counting matches of any length at 95% identity.Figure 1 Distribution of mapped 454 reads with respect to AT content. The genomic reads were mapped to the Amel_4.5 assembly(scaffolds and contigs) using BLAT. With relatively stringent filtering (at least 80% of total length matched and gap size < 30%), 242,284 reads (93%of all reads) were aligned to the assembly. Most reads (236,090, 93%) aligned to fewer than 10 locations, and had unique alignments (210,625,87%). The AT content for each alignment (adding 10% extension on either end) was calculated for reads with≤ 10 match locations.Elsik et al. BMC Genomics 2014, 15:86 Page 4 of 29http://www.biomedcentral.com/1471-2164/15/86without a priori information on the quality of individualinput gene prediction sets, the choice of input datasetsimpacts the results. Evaluating GLEAN runs with differ-ent combinations of input sets allowed us to determinethe optimum selection of datasets. We ran GLEAN 32times with different combinations of input sets, andresulting GLEAN sets ranged in size from in 10,403(with three input gene sets) to 17,724 genes (with seveninput gene sets), (Additional file 2). In some cases,GLEAN predicted a larger number of genes than theMAKER2 set supported by either biological sequencealignment or InterPro. Although GLEAN generates onlyone transcript per gene locus, we decided to use GLEANto create the OGS, because we prioritized GLEAN’spotential identification of a more comprehensive setof gene loci over MAKER2’s identification of 307alternative transcripts. Since with GLEAN we did notneed to filter ab initio genes based on InterPro, it mightallow the genome-conservation-based N-SCAN datasetto provide support for other ab initio predictions evenwithout detectable conservation with known InterProdomains. Another consideration was that assemblymethods for Illumina RNASeq were not well established.We were concerned about the possibility of merginggenes or missing introns due to potential genomic con-tamination in RNASeq experiments. GLEAN uses tran-script alignments to support splice predictions, but doesnot merge or split genes based on transcript alignments,while MAKER2 creates gene models that best agreewith transcript alignments, and works best with highlyreliable transcript assemblies.Selecting an official gene setWe evaluated the 32 GLEAN sets based on several cri-teria, including overlap with a conservative evidence-based set (RefSeq), transcript sequences, peptides and theCEGMA [30] conserved core set (Additional file 2). Nosingle gene set was optimal with respect to all criteria. Wechose to rank sets based on number of peptide matches,which would prioritize completeness of a protein-codinggene set rather than correctness of gene structure.Assessing the new official gene setThe selected GLEAN set, OGSv3.2 (GLEAN31 in Additionalfile 2), represented a significant improvement becauseit included 15,314 protein-coding genes, which is 5,157more genes than the first official gene set, OGSv1.0.The proportion of genes on placed scaffolds as well asthose with expressed sequence coverage also increasedover OGSv1.0 (Table 4). Since GLEAN predicts onlycoding exons, but not untranslated regions (UTRs), weused MAKER2 [28] to add UTR to the final GLEANgene predictions. Out of a total of 15,314 OGSv3.2genes, UTR were added to 7,514 genes (49%).Many split and merged gene models were noted whencomparing the 32 GLEAN sets, including OGSv3.2, to theconservative RefSeq gene set (Additional file 2). Since it isdifficult to computationally resolve disagreements amonggene sets, the OGSv3.2 genes in question are provided as aseparate track on the BeeBase Amel_4.5 genome browser[31]. Web Apollo [32] annotation tools allow registeredBeeBase users to modify gene models, and can be used tomanually correct split or merged genes [33].To address the question of whether the increased genenumber was due to splitting genes, we looked at thetotal number of coding nucleotides (nt) and the distri-bution of coding sequence lengths in OGSv1.0 andOGSv3.2. The ranges in coding sequence lengths were24-53,649 and 75-70,263 nt for OGSv1.0 and OGSv3.2,respectively. Note that the gene prediction parametersFigure 2 GC content of genic regions and overall genomeassemblies. For each gene, the GC content (percent G + Cnucleotides) of genomic regions containing the gene wasdetermined as described in methods. The cumulative distributions ofGC content for overall genome assemblies (thin red line for Amel 4.5and thin black line for Amel_2.0) show that the Amel_4.5 assemblyhas a higher fraction of low GC content regions than does theAmel_2.0 assembly (note the thin red line is to the left of the thinblack line below about 28% GC). The cumulative distributions of GCcontent for the regions containing genes (thick red line for allOGSv3.2, thick green line for Previously known genes, thick blue linefor Type I New genes and thick pink line for Type II new genes)show that regions containing genes are lower in GC content thanthe overall genome. This trend applies for the complete set ofOGSv3.2 genes, as well as the three subsets. The distribution forType I New genes lies to the left of the other distributions, showingthat Type I New genes are located in lower GC content regions thanthe other gene subsets. The distribution for Type II new genes is tothe right of the distributions for the other gene subsets, showingthat the Type II New genes are located in higher GCcontent regions.Elsik et al. BMC Genomics 2014, 15:86 Page 5 of 29http://www.biomedcentral.com/1471-2164/15/86did not allow coding sequences less than 75 nt forOGSv3.2, but a minimum coding sequence parameterwas not applied in the generation of OGSv1.0; howeverthere were only 6 OGSv1.0 genes less than 75 nt. TheOGSv3.2 gene with the largest coding sequence over-lapped a single RefSeq gene prediction (XM_623650.3),and was found to be homologous to genes predicted tocode for titin. The total number of coding nucleotidesincreased from 16,484,776 in OGSv1.0 to 19,342,383 inOGSv3.2. The increase in coding nucleotides by only2,857,607 (17.3%) compared to a 50.8% increase in thenumber of genes suggested that new genes were shorterand possibly a result of splitting OGSv1.0 gene models.However, the distribution of coding sequence length inthe two gene sets suggested that splitting genes was notthe major source of the increased gene number (Figure 3).Although OGSv3.2 possessed a large number of shortcoding sequences (< ~1.5 kb) compared to OGSv1.0, thelarger genes did not appear to be missing. Furthermore,analyzing overlapping gene alignments (described below),showed that 4,735 of the OGSv3.2 genes did not overlapwith genes in OGSv1.0, and thus were not likely to be aresult of splitting OGSv1.0 genes.Biological evidence for the new official gene setIn assessing the gene set, we define biological evidenceto include biological sequence alignments used in geneprediction pipelines (transcript, peptide, protein homo-log) as well as matches to InterPro domains and align-ment to other bee genomes. We define biological geneevidence as all of the datasets considered as biologicalevidence, except alignment to other bee genomes, sincewe did not use information about genes or expressionfor the other bee genomes. Some form of biological geneevidence supported the majority (14,084, 92%) of theOGSv3.2 genes. An additional 754 (4.9%) OGSv3.2 genesaligned to either the Apis florea or Bombus terrestrisgenome for a total of 14,836 (96.9%) of the OGSv3.2genes with biological support. Requiring that transcriptsbe spliced reduces the number of supported genes to13,264 (86.8%) or 14,661 (95.7%) if alignment of thegene to another bee genome is considered as support.Identification and characterization of New genesIdentifying new genes due to improved assembly orimproved gene predictionIn order to compare the old and new OGS (OGSv1.0and OGSv3.2), we mapped OGSv3.2 coding sequencesto the Amel_2.0 assembly because it was the assemblyused to predict OGSv1.0 genes [1], which were latermapped to Amel_4.0. We chose to use Amel_2.0 in thisstep, because generation of Amel_4.0 and intermediateassemblies had involved rearranging and splitting scaf-folds and incorporating unscaffolded contigs using anupdated genetic map [1]. We wished to minimize map-ping errors in our comparison of gene sets by usingAmel_2.0 to determine which OGSv3.2 genes corres-pond to OGSv1.0 genes.We performed the mapping twice using differentalignment criteria, stringent (95% identity and 80% align-ment coverage) and relaxed (95% identity and 50% align-ment coverage). We will report results only for thestringent criteria (Table 5) unless otherwise specified,but have provided the relaxed mapping results in TableS2 in Additional file 1. On the basis of mapping toAmel_2.0 and overlap between the OGSv3.2 coding se-quence alignments with OGSv1.0 gene models, we di-vided the 15,314 OGSv3.2 genes into three sub-sets(Table 5). Any OGSv3.2 gene that did not align to theAmel_2.0 assembly we deemed a “Type I New” gene.Any OGSv3.2 gene that aligned to the Amel_2.0 assem-bly, but whose coordinates did not overlap an OGSv1.0gene by a single coding base pair on the same strand, wedeemed a “Type II New” gene. Finally, any OGSv3.2Table 4 Comparison of OGSv1.0 and OGSv3.2OGSv1.0 OGSv3.2Number of genes 10,157 15,314Number of genes within mapped scaffolds (% of total no. of genes) 5,973 (58.8%) 13,285 (86.8%)Number of genes within un-mapped scaffolds (% of total no. of genes) 4,184 (41.2%) 2,029 (13.2%)Average coding sequence length (bp) 1,623 1,266Average number of coding exons 6.4 5.3Number of single coding exon genes (% of total no. of genes) 795 (7.8%) 2,059 (13.4%)Number of multi-coding exon genes (% of total no. of genes) 9,362 (92.2%) 13,255 (86.6%)Number of genes with spliced EST coverage (% of total no. of genes) 3,039 (29.9%) 12,172 (79.5%)Number of genes with un-spliced EST coverage (% of total no. of genes) 1,734 (17.1%) 11,019 (72%)Number of genes that overlap a protein alignment (% of total no. of genes) 7,940 (78.2%) 6,778 (44.3%)Elsik et al. BMC Genomics 2014, 15:86 Page 6 of 29http://www.biomedcentral.com/1471-2164/15/86gene that both aligned to the Amel_2.0 assembly andoverlapped an OGSv1.0 gene we deemed a “PreviouslyKnown” gene.Of the 5,157 additional genes in OGSv3.2 compared toOGSv1.0, 4,735 genes did not overlap OGSv1.0 genes.The other 422 additional genes in OGSv3.2 were due toeither splitting OGSv1.0 genes (discussed above) or tothe annotation of additional recent paralogs in OGSv3.2,and are included in the total count of 10579 PreviouslyKnown genes. The 4735 genes that did not overlap withOGSv1.0 genes could be classified as 782 genes discov-ered due to the additional sequencing and reassembly ofthe bee genome for the Amel_4.5 assembly (Type I Newgenes; Table 5) and 3,953 genes detected by improvedgene prediction, either through the use of new expressedsequence and protein data or improved gene predictionalgorithms (Type II New genes; Table 5). We could map405 additional Type I New genes to the Amel_2.0 as-sembly if we required only 50% of the gene to be cov-ered (Table S2 in Additional file 1), rather than the morestringent 80% gene coverage reported in Table 5, con-sistent with the Amel_2.0 assembly being less continu-ous than Amel_4.5. This lack of contiguity and resultingfragmentation of genes in the Amel_2.0 assembly likelyimpaired the initial gene prediction efforts. While thefrequency of genes with some form of biological support(transcript, peptide, protein homolog, InterPro domain,and/or alignment to another bee genome) was highestfor Previously Known genes (99.7%), most of the Type INew genes (93.8%) and Type II New genes (89.9%) werealso supported (Table 5).Characteristics of new genesWe analyzed the data that went into the annotation ofOGSv3.2 and evaluated which pieces of evidence con-tributed to the prediction of the genes to understandwhy they were missed in OGSv1.0. To determinewhether genes that were not detected in OGSv1.0 havecommon characteristics that make them more challen-ging to predict, we compared them to the PreviouslyKnown genes. We evaluated features such as tissue ex-pression specificity, coding feature length, GC content,overlap of protein homolog alignments on Amel_4.5 andthe proportion of non-canonical splice sites. Additionalfile 3 provides sources of evidence for each OGSv3.2gene.The mean coding sequence lengths of Type I New(1,172 bp) and Type II New genes (331 bp) were shorterFigure 3 Distribution of coding sequence lengths in OGSv3.2 and OGSv1.0. Histogram plots showing the number of genes having “X”coding sequence length in bins of 20 nt are illustrated using points instead of lines to allow visualization of both distributions. The range incoding sequence length extends to 70,263 and 53,649 in OGSv3.2 (blue) and OGSv1.0 (red), respectively, but this figure zooms in to show lengthsonly up to 5,000 nt. There were 386 and 344 genes with coding sequences longer than 5,000 nt in OGSv3.2 and OGSv1.0, respectively. This figureshows that the increased number of genes in OGSv3.2 is largely due to increased numbers of short genes. The number of larger genes is notdecreased, so gene splitting is not likely a major source of additional genes.Elsik et al. BMC Genomics 2014, 15:86 Page 7 of 29http://www.biomedcentral.com/1471-2164/15/86Table 5 New and previously known OGSv3.2 genesAllOGSv3.2Type I newgenesType II newgenesPreviouslyknown genesNumber of genes (% of total OGSv3.2 genes) 15,314(100%)782 (5.1%) 3,953 (25.8%) 10,579 (69.1%)Scaffold analysis Number of genes within mapped scaffolds (% of no. of gene type) 13,285(86.8%)544 (69.6%) 3,199 (80.9%) 9,542 (90.2%)Number of genes within un-mapped scaffolds (% of no. of gene type) 2,029(13.2%)238 (30.4%) 754 (19.1%) 1,037 (9.8%)CDS analysis Average CDS length 1,266 1,172 330 1,622Average no. CDS Exons 5.3 5.6 2.1 6.5Number of single CDS exon genes (% of no. of gene type) 2,059(13.4%)101 (12.9%) 1,239 (31.3%) 719 (6.8%)Number of multi-CDS exon genes (% of no. of gene type) 13,255(86.6%)681 (87.1%) 2,714 (68.7%) 9,860 (93.2%)Intron analysis Number of introns (% of total OGSv3.2 introns) 66,212(100%)3,585 (5.4%) 4,333 (6.5%) 58,294 (88%)Number of introns validated by EST intron coordinates (% of introns of gene type) 54,514(82.3%)2,573 (71.8%) 1,930 (44.5%) 50,011 (85.8%)Peptide analysis Number of genes with a peptide match (% of no. of gene type) 3,631(23.7%)132 (16.9%) 82 (2.1%) 3,417 (32.3%)Protein analysis No. of genes with overlap to at least one protein alignment (% of no. of gene type) 6,778(44.3%)270 (34.5%) 186 (4.7%) 6,322 (59.8%)No. of genes with overlap to a Dmel protein alignment (% of no. of gene type) 1,205(7.9%)38 (4.9%) 13 (0.3%) 1,154 (10.9%)Total spliced and un-spliced expressedsequence supportNo. of genes with overlap to at least one transcript alignment from any of the ten libraries(% of no. of gene type)13,517(88.3%)704 (90.0%) 2,771 (70.1%) 10,042 (94.9%)Spliced expressed sequence analysis No. of genes with overlap to at least one spliced transcript alignment from each of the tenlibraries (% of no. of gene type)1,062(6.9%)32 (4.1%) 15 (0.4%) 1,015 (9.6%)No. of genes with overlap to at least one spliced transcript alignment from any of the tenlibraries (% of no. of gene type)12,172(79.5%)622 (79.5%) 2,110 (53.4%) 9,440 (89.2%)No. of genes without overlap to any spliced transcript alignments in any of the ten libraries(% of no. of gene type)3,142(20.5%)160 (20.5%) 1,843 (46.6%) 1,139 (10.8%)Genes broadly expressed across four tissues (% of no. of gene type) 2,326(15.2%)60 (7.7%) 95 (2.4%) 2,171 (20.5%)Genes narrowly expressed in only a single tissue (% of no. of gene type) 3,346(21.8%)234 (29.9%) 1,139 (28.8%) 1,973 (18.7%)No. of genes without overlap to any spliced transcript alignments in any of the four tissues(% of no. of gene type)3,632(23.7%)192 (24.6%) 1,985 (50.2%) 1,455 (13.8%)Elsiketal.BMCGenomics2014,15:86Page8of29http://www.biomedcentral.com/1471-2164/15/86Table 5 New and previously known OGSv3.2 genes (Continued)Analysis of alignments to other beegenomesNo. of genes that align to Aflo_1.0 (% of no. of gene type) 13,491(88.1%)566 (72.4%) 2,584 (65.4%) 10,341 (97.8%)No. of genes that align to Bter_1.0 (% of no. of gene type) 12,262(80.1%)527 (67.4%) 1,566 (39.6%) 10,169 (96.1%)Evidence-supported genes No. of genes with overlap to at least one form of biological evidence (% of no. of gene type) 14,084(92.0%)713 (91.2%) 2,930 (74.1%) 10,441 (98.7%)No. of genes that align to Aflo_1.0 and/or Bter_1.0 and/or overlap at least one form ofbiological evidence (% of no. of gene type)14,836(96.9%)734 (93.9%) 3,555 (89.9%) 10,547 (99.7%)GC analysis Number of genes on GC compositional domains >10 kb (% of OGSv3.2 total) 15,224(99.4%)777 (5.1%) 3,923 (25.8%) 10,524 (69.1%)Avg. GC content of compositional domain gene resides in 29.60% 26.40% 32.00% 28.90%ENC analysis Effective number of codons 44.95 41.97 45.69 44.9Genes were mapped to Amel_2.0 assembly with stringent mapping criteria of 80% gene coverage and 95% identity. Biological evidence includes transcript overlap (spliced or un-spliced), peptide hit, protein homologalignment overlap, or InterPro domain presence.Elsiketal.BMCGenomics2014,15:86Page9of29http://www.biomedcentral.com/1471-2164/15/86than that of Previously Known genes (1,623 bp) (P = 2.3 ×10-12 and P < 2.2 × 10-16 respectively) (Table 5). This dif-ference may be due to a higher fraction of single codingexon genes among new genes. Thirteen percent of Type INew genes and 31% of Type II New genes contained onecoding exon, while only 6.8% of Previously Known genescontained one coding exon (P < 2.805 × 10-10 and P < 2.2 ×10-16 for Type I and Type II genes, respectively) (Table 5).The number of canonical versus non-canonical splice siteswas not significantly different between the PreviouslyKnown and Type II New genes (Table S3 in Additionalfile 1).Type I New genes were found to reside in GC com-positional domains with mean 26.4% GC, significantlylower than that of Previously Known genes (28.9%)(P = 2.188 × 10-13), supporting improvement in the as-sembly of the high-AT-content regions. On the otherhand, Type II New genes were found in GC compos-itional domains with a mean GC content of 32.0%,higher than that of than Previously Known genes(P < 2.2 × 10-16), but still slightly lower than the meanGC content of the genome (32.7%)The effective number of codons is an estimate of thedeviation from equality of synonymous codon usage ofall codons of a gene, and ranges from 20 (extreme biaswhere only one codon is used for each amino acid) to 61(no bias, all synonymous codons are used equally) [34].Type I New genes had a significantly lower mean effect-ive number of codons than Previously Known genes,41.97 vs. 44.90 respectively (P = 2.26 × 10-12) (Table 5).This is consistent with the idea that the more extremethe deviation from equal proportions of G + C and A + Tnucleotides in coding sequences, the lower the potentialdiversity of synonymous codons. Consistent with theirlocations in less extreme GC compositional domains,Type II New genes had a significantly higher mean ENCthan the Previously Known genes, 45.69 vs. 44.90 re-spectively (P = 3.923 × 10-05) (Table 5).Expression evidence for new genesThe majority of OGSv3.2 genes were supported by tran-script evidence. When combined, the spliced and un-spliced transcript alignments overlapped 13,517 (88.3%)of OGSv3.2 genes. An analysis of OGSv3.2 gene cover-age by transcript dataset (Table S4 in Additional file 1)showed that the fraction of genes represented in a tran-script dataset ranged from 28.3% (both larvae and testes)to 79.2% (forager brain) (Table S4 in Additional file 1).Both Type I New (79.5%, 622) and Type II New (53.4%,2110) genes were less likely to overlap a spliced tran-script alignment than Previously Known genes (89.2%,9,440) (P = 3.298 × 10-16 and P < 2.2 × 10-16, respectively)(Table 5).Compared to Previously Known genes both Type INew and Type II New genes were more likely to be nar-rowly expressed or overlap transcript alignments fromonly one tissue (P = 2.136 × 10-14 and P < 2.2 × 10-16,, re-spectively) 18.7% of Previously Known, 29.9% of Type INew and 28.8% of Type II New were narrowly expressed(Table 5). Conversely, Previously Known genes (20.5%,2,171) were more likely to be broadly expressed or over-lap transcript alignments from all tissues than eitherType I New or Type II New genes (P < 2.2 × 10-16 forboth tests). 20.5% of Previously Known, 7.7% of Type INew and 2.4% of Type II New genes overlapped tran-script alignments from all tissues. However, the fractionsof narrowly expressed (18.7%) and broadly expressed(20.5%) genes were similar to each other for PreviouslyKnown genes.Homolog alignment evidence for new genesResults suggested that the use of dipteran proteins asthe main source of protein homolog evidence forOGSv1.0 may have been a contributing factor to an in-complete gene set; 44.3% of the OGSv3.2 genes over-lapped at least one protein homolog alignment, but only7.9% overlapped an alignment of a D. melanogaster pro-tein (Table 5). Differences among gene sets further sup-ports the limited value of Dipteran proteins as aprimary source of homolog evidence for A. melliferagene prediction. Both Type I New and Type II Newgenes were less likely than Previously Known genes tooverlap a D. melanogaster homolog alignment (P =1.394 × 10-07 and P < 2.2 × 10-16, respectively). The ef-fect was more drastic for Type II New genes; only 0.3%of Type II New genes overlapped a D. melanogasterhomolog alignment, while the proportions were 4.9%and 10.9% for Type I New and Previously Known genes,respectively.Both Type I New and Type II New genes were lesslikely than Previously Known genes to overlap anyhomolog alignment from another sequenced genome(P < 2.2 × 10-16 for both tests), but the difference wasmore extreme for Type II New genes, potentially imply-ing that a greater number of the Type II New genes arespecific to either A. mellifera, bees, or Hymenopterathan the Previously Known genes (Table 5). 59.8% ofPreviously Known genes overlapped a homolog align-ment, while only 34.5% and 4.7% of Type I and Type IInew genes, respectively, overlapped a homolog alignment(Table 5).Genome alignment evidence for new genesOver 80% of the OGSv3.2 genes aligned to both theAflo_1.0 and Bter_1.0 genome assemblies (Table 5).A notable difference between Type II New genes and theother gene sets was the proportion of genes that wereElsik et al. BMC Genomics 2014, 15:86 Page 10 of 29http://www.biomedcentral.com/1471-2164/15/86supported by genome conservation, but not other sourcesof evidence. Of the 3,555 Type II New genes supportedby any evidence, 17.5% were not supported by othersources. On the other hand, only 2.9% of supported TypeI New genes and 1% of supported Previously Knowngenes were supported by only genome conservation.The remaining 20% of the OGSv3.2 genes have thepotential to be Apis-specific (8% of genes that align toAflo_1.0, but not Bter_1.0) or A. mellifera specific (12% ofOGSv3.2 genes that do not align to Aflo_1.0) (Table 5).Peptide analysisUse of peptide evidencePeptide data were used in four ways. First, peptide datawere used in gene prediction by AUGUSTUS, the onlyprogram used in this study that accepts this source ofgene evidence. Second, peptides were compared to allconsensus gene prediction sets to identify the set withthe highest representation of peptides (described above;Additional file 2). Third, peptide support was used incharacterizing Previously Known, Type I and Type IINew genes. Fourth, venom peptides were used to manu-ally annotate venom genes associated with the sting ofthe honey bee.Gene set comparison with all peptidesPeptide sequences aligned to a greater number of Previ-ously Known genes than to Type I or Type II New genes(P < 2.2 × 10-16 for both tests). Peptides aligned to 32%of the Previously Known genes, but to only 17% of TypeI and 2% of Type II New genes (Table 5 and Additionalfile 3).Venom peptide analysisDespite the lower proportion of new genes compared toPreviously Known genes with peptide matches, analysisof venom peptides (a subset of the peptide data) showedthat the Amel_4.5 assembly provides a significant contri-bution to venom proteome research. 705 unique venompeptides provided biological evidence for 102 genes(described in [35]). Searching the venom mass spectraagainst OGSv1.0 and OGSv3.2 showed that the improvedassembly allowed detection of 21 additional peptides sup-porting 9 new venom protein identifications in OGSv3.2(Additional file 4). Additional tryptic peptides were dis-covered for 7 venom proteins as a result of improvedgene predictions (Additional file 4).Manual annotation of genes supported by the venom pep-tides (Additional file 5) using Apollo [36] showed that mosthoney bee venom genes are fully (76.5%) or partially (19.6%)covered by transcriptome evidence. The putative biologicalfunctions of these genes are described elsewhere [35]. Wediscovered that one of the annotated genes (GB40695) codesfor tertiapin. While the tertiapin peptide was already foundin bee venom [37], no genomic or transcriptomic evidencehad been described, an issue which is now solved as thegenome improvement project supplies both a gene predic-tion (GB40695, NCBI Gene ID 100576769) and EST evi-dence (Genbank:HP466647.1). The gene is positioned onlinkage group 12, next to the apamin and mast cell de-granulating peptide venom genes. The three genes aretandemly arranged which suggest their origin by gene du-plication via unequal crossing over, and may also point toa joint control of transcription [38].Orthology assessmentAnalysis of two different A. mellifera gene sets in com-parison to genes from other insects allowed us to investi-gate the impact of the new genome and gene annotationson the numbers of A. mellifera genes in near-universalorthologous groups. Although true gene losses can anddo occur, these near-universal insect orthologous groupshighlight possible missed gene annotations in eachspecies.Numbers of orthologs were counted for each of the twosets, V2 (similar to OGSv1.0, described in Methods) andOGSv3.2. Figure 4 and Table S5 in Additional file 1 showthe near-universal orthologs missing in each genome. Thepea aphid (Acyrthosiphon pisum) was found to be miss-ing more orthologous groups than the other insects. TheA. mellifera V2 annotation was average for the genomesanalyzed (missing 263 orthologs), but the OGSv3.2 anno-tation was much more complete, missing fewer orthologs(112) than the other genomes. Thus, the new gene setOGSv3.2 reduced the total number of potentially missingorthologs in A. mellifera by 57%, demonstrating the annota-tion improvement that recovered more “universal” insectorthologs than the previous one did.Predicted gene functionsGene ontology analysisGene Ontology analysis (Additional file 6) showed enrich-ment of some specific functions in new genes. Type I Newgenes were enriched for the biological processes “apop-tosis” (corrected e-score = .0016), “neurotransmitter re-lease” (corrected e-score = .0025), and “secondary activeorganic cation transmembrane transporter” (correctede-score < .0001). The detection of these Type I New genesdue to new assembly data is likely due to their location inlow GC content regions, which were underrepresented inthe older assembly; this is in agreement with the lowermean GC content of GC compositional domains contain-ing Type I New genes. Type II new genes were enrichedfor the molecular function “nuclease activity” (correctede-score = .011). Rapidly evolving and single exon genes aremore difficult to annotate automatically, so it is perhapsElsik et al. BMC Genomics 2014, 15:86 Page 11 of 29http://www.biomedcentral.com/1471-2164/15/86not surprising that we found significant differences insome activities.InterPro analysisComparison of InterPro protein domains found inOGSv3.2 and OGSv1.0 proteins showed that 269 of theInterPro domains present in OGSv3.2 were not presentin OGSv1.0 (Additional file 7). There were also signifi-cant differences in the fraction of genes with thedomains IPR004117 (Olfactory receptor, Drosophila,p < .02) and IPR001888 (Transposase, type 1; p < .01).There were 141 and 92 genes with olfactory receptor do-mains in OGSv3.2 and OGSv1.0, respectively. Olfactoryreceptor genes contain a single olfactory receptor do-main, so the domain count corresponds to the genecount. The additional olfactory receptor genes werefound among both the Type I and Type II New genes.The family of olfactory receptor genes is expanded inhymenopteran insects, with rapidly diverging familymembers that are often arranged in tandem arrays[20,21,39]. Many of the 166 known A. mellifera olfactoryreceptor genes were identified by manual annotation ofthe previous assembly, because they were not found inOGSv1.0 [20,21,39]. Improved computational predictionof olfactory receptors in the new assembly was likely dueto increased assembly continuity which would improveidentification of tandemly repeated genes, and newsources of biological sequence alignment evidence (e.g.transcript and protein homolog sequences) which im-prove identification of rapidly diverging genes.There were 18 and 3 genes with transposase type 1domains in OGSv3.2 and OGSv1.0, respectively. Thir-teen of the additional transposase type 1domains were inType 1 New genes. Transposases are associated withtransposable elements, which are a class of interspersedrepeats. Since the new assembly information filled ingaps, and gaps in de novo assemblies are found wherethe assembly process fails to find an unambiguous path(e.g. repetitive loci), we were not surprised to identifyInterPro domains associated with repetitive sequencesamong the Type I New genes. Predicted genes withmatches to transposable element domains are often re-moved from gene sets. We further investigated OGSv3.2for presence of interspersed repeats, and decided not toremove the small number of genes associated with trans-posable element domains, because evidence suggestedthat some were real host genes (see Genome-wide repeatanalysis).PediculushumanusAcyrthosiphonpisumNasoniavitripennisApismelliferaApismelliferaLinepithemahumilePogonomyrmexbarbatusTriboliumcastaneumDanausplexippusAnophelesgambiaeDrosophilamelanogasterV3.2V2 Single-copy in all but 1Single-copy in all but 2Present in all but 1Present in all but 2100 200 300 400 500 600 700Number of Orthologous GroupsFigure 4 Insect orthologs in two A. mellifera gene sets (V2 and OGSv3.2). For each species, counts of near-universal orthologous groups thatare missing an ortholog in that species, or in that species and one other species, are shown. Total counts are divided into groups with onlysingle-copy orthologs and those with gene duplications, further divided into those with only one missing species and those with twomissing species.Elsik et al. BMC Genomics 2014, 15:86 Page 12 of 29http://www.biomedcentral.com/1471-2164/15/86Genome-wide repeat analysisProcessing the Amel_4.5 assembly with the REPET pipe-line yielded 2,401 de novo predicted repetitive elements,of which 1,045 were validated by annotation of at leastone complete copy. In total 9.46% (22.13 Mb) of the gen-ome appears to be repetitive (Table 6). Non-interspersed re-peats (SSR, low complexity, satellite) accounted for 4.05%(9.49 Mb), whereas interspersed repeats represented 5.07%(12.69 Mb) of the Amel_4.5 assembly. The latter estimatefor transposable elements is similar to a previous estimateof 3% [1], whereby only Mariner and R2 elements were re-ported. Thus the honey bee remains a species with an un-usually low amount of repetitive DNA.Most of the groups of retrotransposable elementswere detected in the genome of the honey bee. In com-parison to many other organisms, the most striking dif-ference is the extremely low diversity and abundance ofthese elements. LTR retrotransposons accounted foronly 0.02% of the genome (49.6 kb) and included exam-ples from only one Copia and a putative, unclassifiedelement. Only fragments of elements from the BelPaoand Gypsy superfamily were found. The DIRS groupwas represented by two incomplete elements andaccounted for 0.01% of the genome (12.5 kb). LINEsaccounted for only 0.04% of the genome (83.1 kb), andwere represented by only two R2 elements, one I (nim-bus) element, and a few fragments, potentially belongingto I (R1) and Jockey (CR1). Of the SINE elements de-tected, 14 could not be classified and five had similar-ities to SINEs of the 5S type, all representing 0.03% ofthe genome (70 kb). Together with another unclassifiedelement, all class I retroelements summed up to only224 kb (0.09%) of the genome, among the lowest in theanimal kingdom. TRIMs (terminal repeat retrotranspo-sons in miniature) [40] and LARDs (large retrotrans-poson derivatives) [41] are derivatives of retroelementsand were detected in larger number occupying 9.57 Mb(4.09%) of the genome (Table 6).Class II DNA transposons were more frequent andaccounted for 0.57% of the genome (1.34 Mb). The vastmajority of the elements were TIR (terminal inverted re-peat) transposons of the Mariner superfamily (0.49%,1.15 Mb). Otherwise only two elements of the PiggyBacsuperfamily (0.04%, 88 kb) and five unclassified TIR ele-ments could be detected. Other types of DNA transpo-sons, Crypton, Maverick and Helitron could not befound. The DNA transposon derivatives, MITES (mini-ature inverted transposable elements), were found onlyonce and accounted for less than 0.01% of the genome(3.8 kb) (Table 6).Besides the well-classified sequences, many repetitiveelements could not be assigned to a superfamily or evenclass. The latter includes a larger number of elements(0.65%, 1.52 Mb), which could represent novel types,but need further investigation. We separately annotatedand excluded from the transposable element counts ele-ments that could not be categorized and elements with-out typical transposable element features that containedprofiles from protein coding genes. Both together com-prise 2.96% (6.93 Mb) of the genome (Table 6).The detected elements of the R2 and the majority ofthose of the Mariner type belonged to or were very similarto previously described elements from A. mellifera or otherinsects (RepBase v17.01, [42]). A few other Mariner and thePiggyBac elements showed similarities to elements knownfrom distant animal species. This indicates a potential hori-zontal gene transfer as previously suggested [43].Most of the elements appeared to be fragmented andincomplete. Although some contained sequences of typ-ical transposable element protein domains, they seemedto be inactive due to stop codons and frameshift muta-tions. We detected only four retrotransposons and 27DNA transposons with RT (reverse transcriptase) orTase (transposase) domains, respectively. None of theretrotransposons appeared to possess an active ORFcontaining an entire RT domain, so we classified themas inactive. Among the DNA transposons, six of theMariners appeared to be complete and two were poten-tially active. Five additional Mariner elements possessedan intact ORF spanning at least parts of a Tase domain,so they might have limited activity. The higher abun-dance and higher number of chimeric inserts of Mari-ners (Table 6) suggests that transposons were morerecently active than retrotransposons.Repeats associated with the official gene setSince gene annotation had been performed on an assemblythat was masked for repeats early in the project, prior tothe availability of results from the genome-wide repeat ana-lysis, we expected some OGS genes to overlap newly-detected repeats due to incomplete masking. Coding se-quences of 1,234 genes overlapped interspersed repeats de-tected by REPET. These genes were further investigated forcharacteristics that would support their annotation as hostgenes, including whether the gene was Previously Known,had multiple coding exons, overlapped a spliced transcriptalignment, and had InterPro domain matches (Additionalfile 8). Of these 1,234 REPET-overlapping genes, 739 wereclassified as Previously Known, 185 as Type I New genes,and 310 as Type II New genes. 1,040 of the REPET-overlapping genes had multiple exons and 972 overlappedspliced transcript alignments.Inspecting InterPro results to identify protein domainsknown to be associated with transposable elements, butnot host genes, resulted in the following list of domains:DDE superfamily endonuclease; Integrase, catalytic core;Reverse transcriptase, RNA-dependent DNA polymerase;Retrotransposon, Pao; Ribonuclease H domain; RibonucleaseElsik et al. BMC Genomics 2014, 15:86 Page 13 of 29http://www.biomedcentral.com/1471-2164/15/86Table 6 Repetitive elements in the Apis mellifera genomeType of element No. elements(no. chimeric/nested insert)aNo. Ffagments(no. full length copies)bGenomecoverage (bp)% of genome(234.087 Mb)Ac Bc Cc Dc EcRepetitive DNA 22,134,229 9.46NON-INTERSPERSED REPEATS 94,86,745 4.05SSR 29,697 1,441,651 0.62Low complexity 31,728 8,001,104 3.42Satellite 5 (0) 75 (6) 43,990 0.02 0 na 0 0 0INTERSPERSED REPEATS 881 (65) 28,004 (1102) 12,647,484 5.40 25 7 40 2 6Class I - Retrotransposons 758 (13) 21,244 (903) 9,790,204 4.18 4 1 4 0 0LTR retrotransposons 2 (9) 42 (4) 49,549 0.02 1 1 1 0 0Copia 1 (3) 29 (3) 43,892 0.02 0 1 1 0 0Gypsy 0 (2) 0 (0) 0 0.00 0 0 0 0 0Bel-Pao 0 (4) 0 (0) 0 0.00 0 0 0 0 0Unclassified LTR retrotransposons 1 (0) 13 (1) 5,657 0.00 1 0 0 0 0DIRS retrotransposons 2 (0) 9 (3) 12,472 0.01 0 0 2 0 0LINE (non-LTR) Retrotransposons 3 (4) 140 (3) 83,103 0.04 2 0 1 0 0R2 (NeSL, R2, R4, CRE) 2 (1) 112 (2) 72,107 0.03 2 0 0 0 0Jockey (Rex, Jockey, Cr1, Kiri, L2, crack, Daphne) 0 (1) 0 (0) 0 0.00 0 0 0 0 0I (R1, I, Nimb, outcast, Tad, Loa) 1 (1) 28 (1) 10,996 0.00 0 0 1 0 0Unclassified LINE 0 (1) 0 (0) 0 0.00 0 0 0 0 0SINE 19 (0) 222 (29) 69,938 0.03 0 0 0 0 0SS-Sine 5 (0) 31 (7) 22,660 0.01 0 na 0 0 0Unclassified SINE 14 (0) 191 (22) 47,278 0.02 0 na 0 0 0Unclassified retrotransposons 1 (0) 2 (1) 8,526 0.00 0 na 0 0 0LARD 301 (0) 16,406 (348) 7,256,932 3.10 1 0 0 0 0TRIM 430 (0) 4,423 (515) 2,309,684 0.99 0 0 0 0 0Class II - DNA transposons 51 (52) 3,209 (93) 1,339,131 0.57 7 6 27 2 5TIR 50 (46) 3,200 (89) 1,335,380 0.57 7 6 27 2 5Tc1/Mariner 43 (40) 2,636 (80) 1,147,521 0.49 5 6 25 2 5PiggyBac 2 (6) 184 (2) 87,963 0.04 2 0 2 0 0Unclassified TIR DNA transposons 5 (0) 380 (7) 99,896 0.04 0 na 0 0 0Unclassified DNA-transposons 0 (6) 0 (0) 0 0.00 0 0 0 0 0MITE 1 (0) 9 (4) 3,751 0.00 0 na 0 0 0Unclassified, putative elements 72 (0) 3,551 (106) 1,518,149 0.65 14 na 9 0 1Elsiketal.BMCGenomics2014,15:86Page14of29http://www.biomedcentral.com/1471-2164/15/86Table 6 Repetitive elements in the Apis mellifera genome (Continued)Other DNA elements (not repetitive DNA) 158 (0) 13,760 (250) 6,934,063 2.96 17 0 0 0 0Not categorized 6 (0) 946 (11) 1,233,884 0.53 0 na 0 0 0Potential host gened 152 (0) 12,814 (239) 5,700,179 2.44 17 na 0 0 0For each group, the number of elements (putative families), the number of element fragments or copies in the genome, the cumulative length, the proportion of the genome and other features (elements containingchimeric or nested inserts of other elements (A), elements that appear to be complete with all typical structural and coding parts present even if stop codons or frameshifts are present (B), elements with a RT or Tasedomain detected (C), potentially active elements that contain an intact ORF with all the typical domains although these can lack terminal repeats (D), elements with an intact ORF for the RT domain or parts of theTase domain that could thus be partly active (E) are shown. The elements that could not be categorized or contained features of A. mellifera coding regions are shown at the bottom, these are probably nottransposable elements.aThe numbers of chimeric/nested elements within elements of other categories are not included in the total numbers of elements.bThe software uses alignments to identify the longest fragment, which it deems as full-length. The number of full-length copies is also included in the total number of fragments.cAdditional Columns:A. No. elements containing insertsB. No. complete elementsC. No. elements with RT or Tase domainsD. No. potentially active elementsE. No. potentially partially active elementsdPotential host genes were predicted by software using DNA characteristics, not by overlap analysis with gene predictions. An example of a potential host gene element is a coding sequence for a repeated proteindomain or motif.Elsiketal.BMCGenomics2014,15:86Page15of29http://www.biomedcentral.com/1471-2164/15/86H-like domain; Transposase, Tc1-like; Transposase, type 1;Transposase, Synechocystis PCC 6803. Of the 760 REPET-overlapping genes that had InterPro matches, only 35matched one of these transposable element domains. Someprotein domains, such as zinc fingers, peptidases and heli-cases, are similar between host genes and transposable ele-ments, so cannot be used to classify genes as transposableelements. Among the genes that matched non-transposableelement domains were some known to be members of largegene families, such as olfactory receptors, and genes with re-petitive domains, such as ankyrin repeat-containing domains.We were not surprised that these would be included in a setof de novo detected repeats.Of the 35 genes with transposable element domainmatches, 13 had other domains suggesting they could behost genes, and 28 overlapped spliced transcript align-ments. Only four genes with either transposable elementor unknown/uncharacterized domain family matcheslacked evidence of being a host gene. That is these fourgenes had a single coding exon, lacked spliced transcriptoverlap and had no other InterPro domain match. Wechose not to remove from the OGS genes that overlappedrepeats detected by REPET or genes that matched Inter-Pro transposable element domains, because of evidencesupporting a large number of them as host genes, andthe possibility that transposable elements may havecontributed to the evolution of some host genes (e.g.reviewed in [44]).DiscussionIn the approximately seven years that have elapsed be-tween generation of OGSv1.0 and OGSv3.2, new se-quence was added to the assembly; new sequencingtechnologies, genome assembly methods and gene pre-diction methods were applied; and new sources of geneprediction evidence became available. Several lines of evi-dence indicate that the new assembly is more complete.These include improved continuity of scaffolds, increasedcoverage of low GC content regions of the genome, iden-tification of 782 new genes that do not align to the olderassembly, and increased number of detectable repetitiveelements. The upgraded assembly allowed us to conducta comprehensive analysis of repetitive elements in thegenome, create an improved gene set and confirm previ-ous findings related to GC composition. It also allowedus to confirm that low amounts of transposable elementsand repetitive DNA are bona fide features of the honeybee genome, and not artifacts of incomplete assemblyand annotation.Transposable elements can be a major factor in genomeand gene evolution. Previous analyses found a low numberof transposons and retrotransposons in the assembledgenome compared to other sequenced insect genomes [1],but researchers questioned whether additional elementswere present in unassembled portions of the genome.The upgraded assembly allowed us to better characterizerepetitive elements in the genome. Despite a more com-prehensive repetitive element annotation, the genomiccoverage of transposable elements was extremely low,most striking for the retrotransposons, in agreement withthe previous analyses. The most apparent difference be-tween A. mellifera and other hymenopteran insects is therelative lack of retrotransposable elements; genomes ofother bee, wasp, ant and insect species all contain higherproportions [22,23,45]. Although comparisons acrossstudies are difficult due to methodological differences, ourresults show that the total fraction of repetitive DNA inthe A. mellifera genome (9.46%) is lower than that of mostother sequenced hymenoptera genomes. Some of theant genomes have more than twice as much repetitiveDNA (Atta cephalotes, 25%, Linepithema humile 23.5%,Acromyrmex echinator 27%, Harpegnathos saltator 27%,Solenopsis invicta >23%) and the genome of the parasitoidwasp, Nasonia vitripennis, contains more than three timesthe amount (>30%) [18-20,22-24]. The ants Camponotusfloridanus and Pogonomyrmex barbatus are more similarto A. mellifera, with estimates of 15% and 9%, respectively[18,21]. However, those analyses did not include LARD,TRIM and MITE elements, which make up a considerablefraction of the repetitive elements in A. mellifera (4.1% ofthe genome; 43.3% of the repetitive DNA). Without thesederivative elements, A. mellifera would possess well lessthan 1% transposable elements. This extraordinarily lowproportion of mobile elements suggests that evolutionaryprocesses molding the A. mellifera genome differ fromprocesses working on other hymenopteran genomes, eventhough many of the species listed above are also insectswith eusocial lifestyles.Combined with new biological evidence for gene predic-tions, the upgraded assembly allowed for significant im-provement to the A. mellifera gene set, with >50% moregenes. The identification of 269 additional InterPro do-mains and the 57% reduction in number of missing univer-sal insect orthologs indicate a more comprehensive catalogof protein functions. The presence of nine new venom pro-tein genes, the detection of new tryptic peptides in existingvenom protein genes, and the 53% increase in the numberof computationally identified olfactory receptor domainsare examples of more comprehensive annotation of specificgene families important to bee biology enabled by the im-proved genome.The previously described A. mellifera genome charac-teristics of low and heterogeneous GC content [1] remainafter the addition of new sequence to the assembly. Ex-pansion and improvement of the low GC content regionsin the new Amel_4.5 assembly was supported by theidentification of Type I New genes, which were found inElsik et al. BMC Genomics 2014, 15:86 Page 16 of 29http://www.biomedcentral.com/1471-2164/15/86regions of lower GC content than that of PreviouslyKnown genes. While it is impossible to know whethernew gene prediction evidence (transcript, peptide,homolog, genome conservation) affected the ability topredict the Type I New genes, the identification ofType II New genes in regions with slightly higher GCcontent than Previously Known genes suggested thatthe addition of new gene evidence had a greater impacton gene prediction in higher GC regions. Although anysingle evidence type was less frequent in the Type IINew genes, the total number of Type II New genessupported by the evidence (3,555) was higher than thatof Type I New genes (734). Higher recombination rates[46] and rates of molecular evolution in high GC con-tent regions of the A. mellifera genome [47] may havecontributed to sequence divergence that made Type IINew genes difficult to detect when generatingOGSv1.0. High GC content regions have been shownto be enriched in genes associated with behavioraltraits [47], suggesting that some Type II New genesmay be associated with important bee-specific traits.Despite the higher mean GC content of regions contain-ing Type II New genes, the strong bias for A. melliferagenes to reside in low GC content regions relative to thegenome [1] remains. Among the wide range of insect ge-nomes we have examined thus far, only the hymenop-terans A. mellifera and H. saltator show a strong bias forgenes to occur in compositional domains with low GCcontent, although S. invicta, P. barbatus, L. humile andN. vitripennis show a slight bias [48]. The biological mean-ing of this, and whether this is related to the lifestyles ofthese hymenopterans, is still unclear.Comparisons of the old and new gene sets suggestedthat short and single coding exon genes, with spatially-or temporally-restricted expression patterns and lowprotein homology remain difficult genes to predict. Com-pared with Previously Known genes, Type II New geneswere more rarely and narrowly expressed, had shortercoding sequences, were more likely to be single codingexon genes, and were less likely to have detectable homo-logs. The characteristics of Type II New genes wereconsistent with those of new genes identified in the mod-Encode effort to reannotate the developmental transcrip-tome of Drosophila melanogaster [49]. Their “newtranscribed regions” (NTRs) in Drosophila also had lowexpression levels with temporally restricted patterns. Inaddition more than half of these NTRs were single-exongenes, and the multi-exon NTRs were shorter and lessconserved than previously annotated genes [49]. Whiledifficulty in computational identification of functionalsmall open-reading frames (smORFs; <100 codons) hasled to their low representation in genome annotations,increasing evidence supports functional smORFs in eu-karyotes [50,51].Efficient and effective annotation methods for non-model organisms has become critical in a time when ini-tiatives such as i5K (5000 insect genomes [52]) andG10K (10,000 animal genomes [53]) motivate scientiststo investigate genomes of diverse organisms, many ofwhich will be evolutionarily distant from existing modelorganisms. Our results suggest that investigators wishingto comprehensively annotate protein-coding genes in ge-nomes of non-model organisms should invest in tran-scriptome sequencing that is both deep and broad,similar to findings in re-annotation of the green anolelizard (Anolis carolinensis) genome [54]. Type II Newgenes were more likely to be narrowly expressed, sotranscript evidence from only a single tissue may havemissed a high proportion of the expressed genes. Tran-scriptome sequence from multiple tissues, life stages andconditions may be more useful than protein homologevidence; the number of Type II New genes with tran-script evidence exceeded the number of those with pro-tein homolog evidence, despite the relatively closeevolutionary distance between A. mellifera and some ofthe reference species, which included six species withinthe order Hymenoptera. Even with sampling many tran-scriptomes, rarely expressed genes can be missed, andgenomic sequence from closely related species can aidgene prediction by leveraging nucleotide sequence con-servation; our analysis showed that 752 of the OGSv3.2genes were supported by genome conservation as theironly source of biological evidence. As more genomes aresequenced and annotated, gene prediction should be-come easier. When expression and homolog evidenceare lacking, it is important to train ab initio predictionalgorithms using genes representing the true distributionof gene features such as coding sequence and exonlength, codon usage and GC content. Our analysis ofType I and Type II New genes suggests that the bimodaldistribution of genome GC content in A. mellifera re-sults in at least two classes of genes distinguished by dif-ferent distributions of GC content and codon usage.Ab initio gene predictors may benefit by training withcoding sequences from each class separately.ConclusionsWe have shown that next generation sequencing followedby de novo annotation can substantially improve an unfin-ished first generation genome sequence assembly and an-notation. The upgraded assembly and annotation allowedus to confirm that the honey bee genome contains atypical number of genes relative to other insects. Theimproved honey bee gene set will be invaluable to thehoney bee research community in efforts to elucidate themechanisms behind fundamental biological processes,such as evolution of insect eusociality, as well as agricul-tural issues, such as pollinator health and immunity.Elsik et al. BMC Genomics 2014, 15:86 Page 17 of 29http://www.biomedcentral.com/1471-2164/15/86Furthermore, understanding the reasons genes were notpredicted in OGSv1.0 will lead to more effective gene pre-diction strategies for new genome projects.MethodsGenome sequencing and assemblyWe improved the published genome assembly of 2.7million Sanger reads of the honey bee genome, versionAmel_4.0 [1] by incorporating ABI Solid sequence andRoche 454 paired-end sequence to superscaffold andgap-fill the Amel_4.0 assembly using the Atlas-Linksoftware [55]. The new sequence data as well as theexisting sequence data were used to link the genomecontigs from Amel_4.0 into more contiguous scaffolds.Adjacent contig sequences within these scaffolds wereassessed and overlapping and redundant contigs weremerged.We used the paired-end reads for superscaffoldingand intra-scaffold gap filling. The Amel_4.5 assemblycontains new scaffolds formed from merging existingscaffolds and filling some intra-scaffold gaps with otherscaffolds or contigs.The newly formed scaffolds were anchored to the samelinkage group that their member contigs were anchored toin Amel_4.0. Ten of the new scaffolds could be anchoredto two different linkage groups so a manual break wasinserted to split these scaffolds for consistent anchoring.New data for gene predictionRNAseq dataWe sequenced samples from a number of tissues usingthe 454 Titanium technology. Tissues from ovary, testes,mixed antennae (worker, drones, and queens), larvae,mixed embryos, abdomen, and a combined library frombrain and ovary samples were sequenced.The testes cDNA library was prepared with the Clontech“full-length” amplification protocol. A gel cut of 400 to800 bp fragments was combined with nebulized productsfrom the larger cDNA fragments to generate fragments ofan optimum sized sequencing library for the 454 Titaniumplatform. Other libraries were prepared by isolating totalRNA using Trizol and RNeasy columns followed by mRNAisolation using the Qiagen kit and cDNA genration usingthe Invitrogen Superscript kit #11917-010 with randomprimers. A gel cut of 400 to 800 bp fragments was com-bined with nebulized products from the larger cDNAfragments to generate fragments of an optimum sizedsequencing library for the 454 Titanium platform.Peptide dataHoney Bee peptide atlas We performed liquid chro-matography-tandem mass spectrometric (LC-MS/MS)analysis of bee protein extracts from 253 samples,representing three castes, larvae and virtually all adulthoney bee tissues in both sexes and of different diseasestates. Tissue collection and interpretation of the LC-MS/MS data has largely been described elsewhere [56-61]and is available for public download from the Honey BeePeptide Atlas [62]. The raw data from these studies,amounting to more than 8 × 106 tandem mass spectra,were searched against a six-frame translation (3,596,047forward sequences, with reversed complemented se-quences and common contaminants concatenated) ofAmel_4.5 using Mascot (v2.3, Matrix Science). This exer-cise identified 30,622 unique peptide sequences at a falsediscovery rate of 1%, of which 5,834 peptides mapped toregions of the genome, where genes were previously un-known in Amel_4.0. As these results were based solely ona basic translation of raw genomic sequence, these pep-tides represent only tryptic peptides wholly within a singleexon.Venom peptides We analyzed the honey bee workervenom proteome by integrating a combinatorial peptideligand library (CPLL) with liquid chromatography/Fouriertransform ion cyclotron resonance (LC-FTICR) MS/MS.The collection and interpretation of the MS/MS data aredescribed elsewhere [35]. Gene prediction datasets and asix-frame translation of Amel_4.5 were searched usingMascot (v2.3, Matrix Science). Setting the significancethreshold at p < 0.01 led to a peptide false discovery rateof 5.23% for the search of the AUGUSTUS AU9 gene set,2.51% for the AUGUSTUS AU11 search and 3.77% forthe Amel_4.5 NCBI RefSeq search. MS/MS data generatedfrom the CPLL flow-through fractions, and the elutionfractions separated by Tris-glycine- or Tris-tricine-SDS-PAGE gels were separately searched against the genomesix-frame translation resulting in false discovery rates ofrespectively 1.56%, 3.75% and 3.44%.Annotation methodsSummary of gene prediction setsGenes were predicted using a variety of methods includingthe NCBI RefSeq and Gnomon pipelines, AUGUSTUS,SGP2, GeneID, Fgenesh++ and N-SCAN. Unless statedotherwise below, gene prediction pipelines used the maskedassembly available from NCBI [63], which was generatedusing RepeatMasker [64]. New transcriptome data was usedeither directly as evidence or in the generation of trainingsets for the predictors. The Augustus analysis also incorpo-rated available peptide data, and the N-SCAN analysis lev-eraged nucleotide sequence conservation between the A.mellifera genome and the other bee genomes. SGP2 lever-aged conservation between the A. mellifera genome andthe previously published Nasonia genomes [23] based ontranslated alignments.Elsik et al. BMC Genomics 2014, 15:86 Page 18 of 29http://www.biomedcentral.com/1471-2164/15/86NCBI RefSeq and GnomonThe Amel_4.5 assembly was annotated with NCBI'seukaryotic genome annotation pipeline (version 3.0), as de-scribed at [65]. The NCBI pipeline uses a repeat-masked as-sembly generated by WindowsMasker [66] and transcriptand protein alignment evidence supplemented with abinitio prediction to annotate coding and non-coding tran-scripts and proteins. The evidence used for this annotationrun included alignments of:1) 9 million RNAseq reads from 454 sequencing(described above), which were treated as ESTs in thepipeline.2) All available honey bee ESTs and mRNAs.3) Proteins from the FlyBase annotation of Drosophilamelanogaster (release 5.30).4) Proteins from the RefSeq annotation of human.5) Proteins annotated on insect mRNAs.6) Proteins from the annotations of the ant genomes,Harpegnathos saltator and Camponotus floridanusavailable in GenBank.The final RefSeq annotation included models that werecompletely or partially supported by alignment evidence.The pipeline had some provisions to predict models thatwere disrupted by frameshifts or stop codons in the assem-bly, some of which were assigned protein accessions (XP_prefixes) and named with the prefix ‘LOW QUALITYPROTEIN’, whereas others were conservatively classified aspseudogenes and not assigned protein accessions. TheRefSeq annotation is available from NCBI's genomes FTPsite as NCBI build 5.1 [63].AUGUSTUSAUGUSTUS can be used as an ab initio gene predictiontool, but can also integrate extrinsic evidence from vari-ous sources [67].Generation of training gene structures and trainingAUGUSTUSWe used Scipio [68] to generate gene structures on theAmel_4.5 assembly with a protein set from a previous A.mellifera annotation (Amel_pre_release2_OGS_pep.fa).We used these gene structures to optimize AUGUSTUSparameters for A. mellifera, and constructed UTR modelsfrom 78,274 A. mellifera ESTs from GenBank, which wereused to train AUGUSTUS UTR parameters. We then pre-dicted genes in the Amel_4.5 assembly. Individual RNA-seq reads from 454 sequencing (described above) weremapped against predicted transcripts, and fully coveredtranscripts (2,151) were selected as training genes for opti-mizing a final AUGUSTUS parameter set.Gene prediction with AUGUSTUSWe generated three gene sets using AUGUSTUS: a geneset with extrinsic evidence from ESTs and RNAseq data(AU9), a particularly inclusive gene set that containedmany alternative transcripts for peptide identification(AU11), and a gene set with extrinsic evidence fromESTs, RNAseq data and peptide data (AU12) (Table S6in Additional file 1).We created extrinsic evidence “hints” for protein cod-ing genes and transcripts from the A. mellifera ESTs andfrom 454 transcriptome libraries. Genes were predictedwith AUGUSTUS, allowing the prediction of alternativetranscripts and allowing the splice site AT-AC (inaddition to GT-AG and GC-AG) in case of supporting ex-trinsic evidence. The resulting gene set was named AU9.AU9 genes were predicted using the following options:2augustus –species = honeybee1 –UTR= on –print_utr =on –hintsfile = all_but_no_peptide.hints –extrinsicCfgFile =extrinsic.M.RM.E.W.cfg –exonnames=on –codingseq= on –alternatives-from-evidence= true –allow_hinted_splicesites =atac genome.faFor the purpose of peptide identification, alternativetranscripts were extensively sampled with AUGUSTUS,resulting in the gene set AU11. AU11 genes werepredicted issuing the following command: augustus –UTR= on –print_utr = on –hintsfile = all_but_no_peptide.hints –extrinsicCfgFile = extrinsic.M.RM.E.W.cfg –exon-names = on –codingseq = on –species = honeybee1 –alter-natives-from-evidence = true –alternatives-from-sampling =true –sample = 100 –minexonintronprob = 0.1 –minmeanex-onintronprob = 0.4 –maxtracks = -1 –allow_hinted_splice-sites = atac genome.faWe generated hints from peptide sequences by map-ping the peptides against the protein set of AU11 andagainst a six-frame translation of the genome usingBLAT [69] in a non-redundant way (i.e. parts of the six-frame translation that were included in the protein setAU11 and redundant parts of the AU11 protein set thatwere removed). A final AUGUSTUS gene set AU12 wasgenerated using all available extrinsic evidence and the fol-lowing options: augustus –species = honeybee1 –UTR=on –print_utr = on –hintsfile = all.hints –extrinsicCfgFile =extrinsic.M.RM.E.W.cfg –exonnames = on –codingseq= on –alternatives-from-evidence = true –allow_hinted_splicesites =atac genome.fa.Fgenesh++Predictions were made using FGENESH 3.1.1 [70,71]using the HBEE matrix with parameters specific forA. mellifera. We used Illumina transcriptome datafrom A. mellifera forager and nurse brains availablefrom the SRA (SRP003528), which were a total of181.8 M spots, 18.3G bases of 100 bp single end readsgenerated on an Illumina Genome Analyzer II in 2010.Elsik et al. BMC Genomics 2014, 15:86 Page 19 of 29http://www.biomedcentral.com/1471-2164/15/86We mapped the Illumina RNASeq using the ReadsMapprogram [72], which provided intron position informa-tion to Fgenesh for building sample-specific genemodels. Predictions were made based on individualsused in Illumina RNASeq libraries (five nurses and fiveforagers), and then were combined into forager andnurse group predictions, and redundant genes (thosehaving coinciding coding sequences) were removedfrom each set. Finally, the forager and nurse predic-tions were combined into a final set of predictions, andredundant genes (those having coinciding coding se-quences) were removed from this combined set to pro-duce the final Fgenesh++ prediction set.GeneIDGeneID is an ab initio gene prediction program used tofind potential protein-coding genes in anonymous gen-omic sequences. The training of GeneID to obtain a par-ameter file for A. mellifera was based on the methoddescribed to obtain a Drosophila melanogaster GeneIDparameter file [73]. Training was performed in a “semi-automated” manner by employing a recently developedGeneID training tool that computes position weightmatrices (PWMs) or Markov models of order 1 forsplice sites and start codons, and derives a model of cod-ing DNA, which, in this case, is a Markov model of order5. Furthermore, once a preliminary species-specific matrixis obtained it is further optimized by adjusting two in-ternal matrix parameters: -the cutoff of the scores of thepredicted exons (eWF) and the ratio of signal to codingstatistics information to be used (oWF).The initial A. mellifera training set was comprised of the2,151 gene models used to train AUGUSTUS (describedabove). Of these gene models 80% (1,720) were used totrain GeneID while the remaining 20% (431) were set asideto test the accuracy of the newly developed matrix. The1,720 A. mellifera protein-coding gene models included7,913 canonical donor splice sites/7,939 canonical acceptorsites and 1,720 start codons. The start codons were used tocompute PWMs while the donor and acceptors were usedto derive Markov matrices of order 1. Given the large num-ber of sequences we also had enough coding (2,338,800)and non-coding (5,561,154) nucleotides to derive a Markovof order 5 for the coding potential. We tested accuracy ofthe GeneID A. mellifera parameter file on an artificial con-tig consisting of the 431 evaluation-set concatenated genemodels with 800 nucleotides of intervening sequence be-tween each of the genes (Table S7 in Additional file 1). Wethen used the GeneID parameter files to predict genes onan assembly consisting of Amel_4.5 sixteen chromosomesand 5,304 "unplaced" scaffolds files that had repeat se-quences masked using Repeatmasker. GeneID predicted24,554 protein-coding genes.SGP2SGP2 is a syntenic gene prediction tool that combines abinitio gene prediction (GeneID) with TBLASTX searchesbetween two or more genome sequences to provide bothsensitive and specific gene predictions, and it tends to im-prove GeneID’s performance, especially by reducing thenumber of false-positive predictions. SGP2 requires oneor more reference genomes to which the target genome(in this case A. mellifera) is compared. We decided to usethe genomes of Nasonia vitripennis, N. giraulti and N.longicornis [23] as references to develop our A. melliferaparameter file for SGP2 because the genus Nasonia is atan appropriate evolutionary distance from Apis such thatmostly the coding regions of the genes, not the introns orintergenic regions, are significantly conserved betweenthese two genomes.Obtaining the SGP2 A. mellifera-specific parameterfile was based on the methodology described by Parraet al [74] used to obtain a human SGP2 parameterfile using mouse homology evidence. The startingpoint to obtain a parameter file for SGP2 was thepreviously described GeneID A. mellifera matrix. TheGeneID-derived SGP2 parameter file was optimizedon an artificial contig comprising the same concatenated1,720 sequences used to train GeneID, with 800 nucleo-tides between each of the gene models. We optimized theSGP2 matrix by modifying not only the eWF internal par-ameter (as previously for the GeneID parameter file) butalso two SGP2-specific internal parameters (“NO_SCORE”and “HSP_factor”). The “NO_SCORE” parameter providesa penalty for no overlap between TBLASTX-derived HSPs(High-Scoring Pairs) and GeneID ab initio predictions inthe same region whereas the “HSP_factor” parameter re-duces the score assigned to the HSPs in order to maximizethe prediction accuracy. We evaluated the newly devel-oped SGP2 parameter file on the same artificial contigconsisting of the 431 concatenated gene models with 800nucleotides of intervening sequence between each of thegenes used to evaluate the GeneID bee matrix (Table 1).We then used the SGP2 parameter file2 to predict geneson the same assembly file used for GeneID, and generated20,179 predictions.N-SCANWe used the N-SCAN package [75] to leverage conser-vation between the A. mellifera genome and genomes oftwo other bee species, A. florea (Aflo_1.0 ) and Bombusterrestris (Bter_1.0). We first masked Amel_4.5 for sim-ple sequence repeats using RepeatMasker [64]. We ranLASTZ [76] using default parameters with Amel_4.5 asthe target genome, and either Apis florea (Aflo_1.0) orBombus terrestris (Bter_1.0) as the informant genome.We then used iParameterEstimation to generate both anElsik et al. BMC Genomics 2014, 15:86 Page 20 of 29http://www.biomedcentral.com/1471-2164/15/86Amel_4.5-Aflo_1.0 specific parameter set as well as anAmel_4.5-Bter_1.0 parameter set using the training setdescribed above for AUGUSTUS gene prediction, in-cluding UTR features. Finally, we ran N-SCAN usingeach of the A. mellilfera specific parameter sets with therespective LASTZ informant alignments to produce twoN-SCAN prediction sets, one set based on Aflo_1.0 asthe informant genome and the other set based onBter_1.0 as the informant genome.Combining gene setsInput data for combined gene setsWe used MAKER2 and GLEAN to generate combinedgene sets. The MAKER2 and GLEAN analyses used thesame set of input data. Both analyses combined thegene predictions described above (NCBI, AUGUSTUS,Fgenesh++ with RNAseq, N-SCAN using A. florea as aninformant genome, GeneID and SGP2) with transcript andprotein homolog alignments. Transcript data included thenew 454 transcriptome data described above, A. melliferaESTs from GenBank and Illumina nurse and forager readsdownloaded from the SRA (SRP003528, described above).We aligned Illumina reads to Amel_4.5 in two groups,nurse and forager, using Tophat version 1.3.1 with theoption "–butterfly-search" for more sensitive splice junc-tion detection, and then generated predicted transcriptsfor each set of pooled data using Cufflinks version 1.0.3with default parameters. The 454 reads were assembledinto contigs de novo using Newbler (2.3-PreRelease-9/14/2009) with the cDNA option. We aligned 454 contigsand ESTs to Amel_4.5 using MAKER2 v2.15, which usesWU-BLAST [77] and Exonerate est2genome [78], withminimum 80% alignment coverage and 95% identity.Protein homolog alignments included SwissProt [79]Metazoa homologs, Drosophila melanogaster (fruit fly;r5.31) [80], Nasonia vitripennis (parasitoid wasp;OGSv1.2) [23] and the ants: Atta cephalotes (OGSv1.1)[22], Camponotus floridanus (OGSv3.3) [18], Harpeg-nathos saltator (OGSv3.3) [18], Linepithema humile(OGSv1.1) [20], Pogonomyrmex barbatus (OGSv1.1)[21]). Proteins in the SwissProt dataset annotated astransposable elements were removed prior to alignment.We aligned protein sequences to Amel_4.5 using Exone-rate protein2genome with a minimum 60% percentidentity and 60% alignment coverage.MAKER2To create a combined gene set, we ran MAKER2 v2.15using parameters min_contig = 1000 and pred_gff, whichallowed us to provide as input the gene prediction setsand alignments described above, instead of generatingnew evidence tracks within MAKER2.GLEANWe ran GLEAN [17] 32 times to create consensus genesets using different combinations of the gene predictionsets described above. All of the GLEAN runs includedthe transcript and protein homolog alignments describedabove, and required a minimum coding sequence lengthof 75 nt.Selecting the new official gene SetWe evaluated the 32 GLEAN sets based on several cri-teria including overlap with a conservative evidence-based set (RefSeq), transcript sequences, peptides andthe CEGMA conserved core set [30] (Additional file 2).NCBI’s RefSeq pipeline has been considered reliable andrelatively conservative, so we performed overlap analysesto determine how many of the RefSeq models were cap-tured in the GLEAN consensus sets. We used FASTA[81] to align GLEAN coding sequences with RefSeq cod-ing sequences where a perfect alignment required 99%identity over the entire lengths of both sequences. Weused FASTA to align peptides to GLEAN proteins, withE-value and scoring matrix optimized for short exactmatches (E-value .01 and MD10 scoring matrix). Weparsed peptide alignments to count those with 100%identity and 100% alignment coverage (Additional file 2).In addition to alignments between GLEAN predictionsand RefSeq genes and peptides, we considered numbersof gene models, total numbers of coding nucleotides,average coding sequence lengths, and numbers of splitsand merges compared to RefSeq. We selected theGLEAN31 set (Additional file 2) to be the official geneset, OGSv3.2. It was the set generated using gene predic-tions from NCBI Gnomon, AUGUSTUS, Fgenesh++,N-SCAN using A. florea as an informant genome andGeneID and SGP2 combined as a single set as well as thetranscript and protein homolog alignments.Adding UTRs to OGSv3.2 gene modelsAlthough we did not use MAKER2 to generate the offi-cial gene set, we did use MAKER2 (v2.15) [28] to addUTR to the GLEAN coding exons since GLEAN doesnot produce annotations with UTR features. OGSv3.2coding exons were input as “pred_gff” and transcriptomeevidence (including the 454 transcripts, Illumina contigsand A. mellifera ESTs described above) was input as“est”. MAKER2 aligned the EST evidence to Amel_4.5using WU-BLAST and Exonerate est2genome then usedoverlapping EST evidence to extend OGSv3.2 models toinclude UTR features when possible. Because MAKER2sometimes modified the coding exon coordinates, weprocessed the MAKER2 output and retained UTRs onlywhen the coding exon structure was unchanged. Out ofa total of 15,314 OGSv3.2 genes, UTR were added to7,514 genes (49%).Elsik et al. BMC Genomics 2014, 15:86 Page 21 of 29http://www.biomedcentral.com/1471-2164/15/86Identifying and characterizing of new genesIdentifying of new and previously known OGSv3.2 genesIn order to compare the newest gene set (OGSv3.2) withthe first official gene set (OGSv1.0), we mapped OGSv3.2coding sequences to Amel_2.0, which was the assembly onwhich OGSv1.0 was generated [1]. First, we used Mega-BLAST [82] to identify scaffold/gene matches with 95%identity and E-value < 1 × 10-20. We then aligned coding se-quences to matching scaffolds using GMAP [83], andparsed the output to create two sets of splice-modeledalignments, both requiring 95% identity. One set was basedon a relaxed coverage criterion, requiring that the align-ment cover at least 50% of the coding sequence. The otherset was based on a stringent coverage criterion, requiringthat the alignment cover 80% of the coding sequence. Re-sults of further analyses for the relaxed mapping set areprovided in the supplemental materials, but we discuss onlyresults for the stringent mapping.On the basis of mapping OGSv3.2 to Amel_2.0 andoverlap between OGSv3.2 and OGSv1.0 gene models onthe Amel_2.0 assembly, we divided 15,314 OGSv3.2genes into three sub-sets (Table 6). We deemed anyOGSv3.2 gene that did not align to the Amel_2.0 assem-bly a “Type I New” gene. The additional sequencing andreassembly of the genome for the Amel_4.5 assemblylikely allowed the detection of these genes. “Type IINew” genes were those that did align to the Amel_2.0assembly, but whose coordinates did not overlap anOGSv1.0 gene by a single coding base pair on the samestrand. Additional expressed sequence and proteinhomolog evidence as well as improvements to gene pre-diction algorithms likely contributed to the detection ofthese genes. Finally, any OGSv3.2 gene that both alignedto the Amel_2.0 assembly and overlapped an OGSv1.0gene was deemed a “Previously Known” gene.Coding sequence length analysisThe total length of the coding sequence was calculatedfor each gene and the means for all genes and each genesub-set were calculated (Table 6). We tested the null hy-potheses that the mean coding sequence lengths of TypeI and Type II New genes and Previously Known geneswere equal.Splice site and single versus multiple coding exongene analysisWe assessed the genomic sequence of the two intronic basepairs adjacent to each coding sequence exon-intron splicesite to determine whether they corresponded to the canon-ical …]5’-GT/AG-3’[… splice site sequence. We consideredonly splice sites supported by matching intron coordinatesof spliced transcript alignment evidence. We used chi-square tests with one degree of freedom to compare thefrequencies of single coding exon genes and non-canonicalsplice sites in Type I or Type II New genes with PreviouslyKnown genes.OGSv3.2 gene location relative to GC compositionaldomainsWe used IsoPlotter, a recursive segmentation algorithm[84,85], to partition the A. mellifera genome into GCcompositionally homogeneous domains, contiguous re-gions of the genome with similar GC content (percentG + C nucleotides). We determined the GC content ofthe GC compositional domain for each OGSv3.2 gene.In cases where a gene spanned multiple GC compos-itional domains, we calculated the average GC contentof the GC compositional domains, and weighted it ac-cording to the fraction of the gene's length that occursin each GC compositional domain. We computed theweighted percent GC based only on non-coding nucleo-tides of the compositional domains, because we did notwish to include effects of codon bias.We compared the distribution of Type I and Type IINew genes with respect to GC content to that of Previ-ously Known genes. IsoPlotter cannot segment scaffoldsless than 10 kb into compositional domains, so the 90genes residing in these short scaffolds were not consid-ered. We tested the null hypotheses that the means ofthe weighted GC compositional domain contents forgenes in the Type I and Type II New gene sets wereequal to the mean of the Previously Known genes.Effective number of codons analysisUsing the chips program within the EMBOSS package[86], we calculated the effective number of codons sepa-rately for each OGSv3.2 gene. We tested the null hypo-theses that the mean effective number of codons valuesfor the Type I and Type II New genes were equal to themean of the Previously Known genes.Tissue expression analysisWe determined the number of OGSv3.2 gene overlapsto transcript alignments that were used in creating theGLEAN consensus gene sets (454 reads, Illumina con-tigs and A. mellifera ESTs from GenBank, describedabove). We relied on splice signals to determine the di-rectionality of a transcript read. We tallied spliced andunspliced alignments separately, since we could beconfident when directionality of a spliced alignmentagreed with a gene prediction, but could not be confidentabout unspliced alignments. For spliced transcript align-ments, if the transcript was on the opposite strand fromthe gene then it was discarded from further analysis. Fortranscripts on the same strand or transcripts that wereun-spliced, in which case directionality could not be de-termined, a coordinate overlap of at least one coding basepair was required for a gene to count as overlapping withElsik et al. BMC Genomics 2014, 15:86 Page 22 of 29http://www.biomedcentral.com/1471-2164/15/86a transcript alignment. We counted the number of tran-script data sets in which each OGSv3.2 gene was foundto have an overlapping alignment (Table 6.) We per-formed chi-square tests with one degree of freedom tocompare the frequencies of spliced transcript overlap inthe Type I or Type II New gene sets with the PreviouslyKnown gene set.Of genes that overlapped spliced transcript align-ments, we identified genes that were narrowly expressedand genes that were broadly expressed on the basis ofoverlap to the four single-tissue libraries (brain [com-bined Illumina forager and nurse brain libraries] and454 libraries of mixed antennae, ovary and testes). (For-agers and nurses are worker honey bees that specializeon collecting food and feeding brood, respectively.)Genes were deemed narrowly expressed if they over-lapped at least one transcript alignment in only one ofthe four tissues and broadly expressed if they over-lapped at least one transcript alignment from all fourtissues. We performed chi-square tests with one degreeof freedom to compare the frequencies of narrowlyexpressed genes and broadly expressed genes in theType I or Type II New gene sets with the PreviouslyKnown gene set.Homolog analysisWe determined the number of protein alignments over-lapping OGSv3.2 for protein data sets that were used increating the GLEAN consensus gene sets. We requiredoverlap of at least one coding base pair on the samestrand to deem a gene overlapping with a protein homo-log alignment. We performed chi-square tests with onedegree of freedom to compare the frequency of the exis-tence of overlaps to homolog alignments for Type I orType II New genes with that of Previously Known genes.TBLASTN alignment of OGSv3.2 to A. florea and B. terrestrisgenomesWe aligned OGSv3.2 protein sequences to the A. florea(Aflo_1.0) and B. terrestris (Bter_1.0) genome assembliesusing TBLASTN [87] with an E-value criteria of 1 × 10-06.We did not use information about predicted genes or ex-pression in A. florea or B. terrestris.Statistical methods to test for differences in means betweennew and previously known genesTo test the null hypothesis that the mean of a particularvariable for Type I or Type II New genes was equal tothe mean of the Previously Known genes, we used both anon-parametric Kolmogorov-Smirnov test and a Welcht-test with the correction for non-homogeneity of vari-ances. Testing the hypotheses in this way avoids assum-ing these data were normally distributed or had equalvariances. To be conservative, we report the least signifi-cant P-value for each test.Using peptide data in development and analysis of genesetsWe used the peptide data is described above to evaluatethe GLEAN sets and analyze Type I New, Type II New,and Previously Known genes. We aligned the peptide se-quences to predicted protein sequences with FASTA[81] with a relaxed E-value of 0.1 and the MD10 scoringmatrix. These parameters were found to allow matchesto peptide sequences as short as 6 amino acids with100% identity. Only alignments with 100% identity wereretained. We used chi-square tests with one degree offreedom to compare the frequency of proteins thataligned to peptide sequences in Type I or Type II Newgenes with that of Previously Known genes.The venom peptide data were used separately to evalu-ate gene sets for improved identification of venomgenes. All significant and top ranking venom peptidesfrom the Mascot output, with an ion score ≥30 wereretained in the final peptide lists. Venom mass spectrawere searched against the OGSv3.2 and OGSv1.0 tocompare the contribution of both datasets to the identi-fication of new venom genes. The false discovery rates ofboth searches were set to 1%.Venom peptide data were used to annotate venomgenes using the Apollo annotation tool [36], providedfor Amel_4.5 by the Hymenoptera Genome Database[88].Assessing orthologyThe protein-coding gene annotations from the two A. mel-lifera genome assemblies were compared with orthologsfrom OrthoDB [89] from nine other insects. These includedPediculus humanus (PhumU1.2, 10,772 genes), Acyrthosi-phon pisum (ACYPI v2.1, 36,275), Nasonia vitripennis(nvit2, 24,369), Linepithema humile (OGSv1.2, 16,048),Pogonomyrmex barbatus (OGSv1.2, 17,100), Tribolium cas-taneum (Tcas 3.0, 16,565), Danaus plexippus (OGS2,15,329), Anopheles gambiae (AgamP3.6, 12,670), and Dros-ophila melanagaster (r5.45, 13,927). For annotations on theolder A. mellifera genome assembly, we used the gene setknown to the community as “Amel_prerelease2” (abbrevi-ated as V2), available from BeeBase. It was the gene setresulting from mapping OGSv1.0 from assembly Amel_2.0to assembly Amel_4.0, and included a small number ofmanual annotations, with a total of 10,699 genes. Compar-ing each Apis annotation (V2 and OGSv3.2) to the othernine gene sets identified near-universal orthologous groupswith orthologs in all but one or all but two insects. For eachgene set, Figure 4 and Table S5 in Additional file 1 showcounts of near-universal insect orthologous groups that aremissing orthologs in different species. Total counts wereElsik et al. BMC Genomics 2014, 15:86 Page 23 of 29http://www.biomedcentral.com/1471-2164/15/86partitioned into groups with only single-copy orthologs andthose with gene duplications, further divided into thosewith only one missing species and those with two missingspecies.Predicting protein functionsGO analysisWe used FASTA [81] with an E-value threshold of 1 ×10-6 to compute reciprocal alignments between OGSv3.2proteins and a D. melanogaster protein set consisting ofthe longest protein isoform of each gene (annotationversion r5.42). We identified reciprocal best hits (RBH)and transferred Gene Ontology (GO) [90] annotationsfrom the D. melanogaster protein to the A. melliferaprotein for each RBH pair, using the GO annotation fileavailable at FlyBase [80]. We used GeneMerge [91] totest for enrichment of GO terms in OGSv3.2 proteindatasets, testing for each of the three GO ontologies(Molecular Function, Biological Process and CellularComponent) separately. Several tests were performedwith either the entire set of Gene Ontology terms, thegeneric GO slim set, or the GO slim set developed forA. mellifera by Whitfield et al. [92]. Population and testgene datasets for a particular GO ontology included onlyOGSv3.2 genes with GO annotations for that ontology.Population datasets consisted of all OGSv3.2 genes withGO annotations. Test datasets were 1) all new genesbased on stringent mapping criteria, 2) Type I Newgenes using stringent criteria, 3) Type II New genesbased on stringent mapping criteria.InterPro analysisWe used InterProScan [93] to compare OGSv3.2 andOGSv1.0 proteins with the following InterPro [29] proteindomain and motif databases: PFAM [94], TIGRFAMS [95],SMART [96], PRODOM [97], PROSITE [98], PIRSF [99],GENE3D [100], SUPERFAMILY [101], and PANTHER[102]. The total numbers of proteins annotated with at leastone InterPro domain were 9,479 and 8,552 for OGSv3.2and OGSv1.0, respectively. For each InterPro domain iden-tified in the combined datasets, we determined the numberof proteins containing that domain within OGSv3.2 andOGSv1.0. Then for each InterPro domain, we used 2 × 2chi-square tests with Yates correction and one degree offreedom to determine whether the frequencies of proteinscontaining that domain differed between the OGSv3.2 andOGSv1.0 InterPro-annotated sets (total 9,479 and 8,552proteins, respectively).Detecting genome-wide repetitive elementsWe detected and annotated repetitive elements with theREPET software package ([103], version 2.0) consistingof two pipelines integrating a set of bioinformaticsprograms. First, repeated sequences were detected bysimilarity, using an all-by-all BLAST [104] search viaBLASTER [105]. LTR-retrotransposons were detected bystructural search with LTRharvest [106]. The similaritymatches were clustered with GROUPER [105], RECON[107] and PILER [108], the structural matches were clus-tered with NCBI BLASTclust [109]. From each cluster aconsensus sequence was generated by multiple align-ment with Map. We analyzed the consensus sequencesfor terminal repeats (TRsearch), tandem repeats (TRF),open reading frames (dbORF.py, REPET) and poly-Atails (polyAtail, REPET). We screened the consensusesfor matches to nucleotide and amino acid sequencesfrom known transposable elements (RepBase 17.01, [42])using BLASTER [105], which runs tblastx and blastx[87], and searched them for HMM profiles (Pfam data-base 26.0, [94]) using HMMER3 [110]. Based on the de-tected structural features and homologies, we classifiedthe consensuses using PASTEC according to Wickeret al. [111]. We then removed redundancies identifiedwith BLASTER and MATCHER [105] as well as ele-ments classified as simple sequence repeats (SSRs; >0.75SSR coverage) or unclassified elements built from lessthan 10 fragments.We used the set of de novo detected repetitive ele-ments to mine the genome in a second pipeline withBLASTER (using NCBI BLAST, sensitivity 4, followedby MATCHER), RepeatMasker (using CrossMatch, sen-sitivity q, cutoff at 200) and CENSOR (using NCBIBLAST). We removed false positive matches using anempirical statistical filter. Satellite repeats were detectedwith TRF [112], MREPS [113]and RepeatMasker [64]and were then merged into a single set. We alsoscreened the genomic sequences for matching nucleotideand amino acid sequences from known transposable ele-ments (RepBase 17.01, [42]) with BLASTER. Finally weremoved TE doublons (loci annotated as multiple trans-posable elements) and SSR annotations within trans-posable element annotations, and performed a "long joinprocedure" to connect distant fragments. Sequencesfrom the de novo repetitive element library, which werefound to have at least one perfect match in the genomewere then used to rerun the whole analysis.To ensure compatibility and to avoid introducing a bias,we refrained from a manual curation or clustering of thede novo detected elements before mining the genome.However, post hoc we manually analyzed all elements,which were previously classified into class I retrotransposonor class II DNA transposon elements or unclassified ele-ments with detected coding element features (similarity toknown transposable elements) due to potential chimeric in-sertion. At this stage we excluded derivative elements(LARD, TRIM, MITE) from further detailed inspection un-less carrying a class I or class II element. Some elementsElsik et al. BMC Genomics 2014, 15:86 Page 24 of 29http://www.biomedcentral.com/1471-2164/15/86were classified “potential Hostgene” in the computationalanalysis, based on characteristics of the DNA, not based onoverlap analysis with predicted genes. These “potentialHostgene” elements, as well as unclassified elements(“noCat”), were also excluded from manual analysis. Weperformed manual inspection by checking for open-readingframes (ORF) with the NCBI ORF Finder (NCBI), bysearching the NCBI Conserved Domain Database (CDD)[114], by searching the most up to date online RepBasedatabase (accessed December 2012-February 2013) viaCENSOR [115]. We also performed phylogenetic analysisfor LINE RT domains with RTclass1 [116] in order toachieve a detailed classification for each element, determineits potential relation to a family of known elements, toevaluate the completeness and to detect potential active ele-ments. We defined an element to be complete, if it pos-sessed the relevant coding parts with the element-typicaldomains and the structural features (LTR, TIR). The poten-tial activity was defined according to the region an intactORF, if present, covered. If an intact ORF seemed to covera complete region including the typical domains (e.g. GAGas well as POL, Tase) then the element was considered tobe potentially active. If a Tase domain was covered by atruncated ORF or the Tase itself appeared to be truncatedbut was covered by an intact ORF, or if the RT domain wascovered by an active ORF but not the remaining element-typical domains, then the element was considered to bemaybe potentially active. During the manual classificationto at least superfamily level, novel transposable elementtypes not covered by the system of Wicker et al. [111] werealso considered: Kolobok, Sola, Chapaev, Ginger, Academ,Novosib and ISL2EU class II DNA transposons [117,118].Simple sequence repeats and other low complexity re-gions were extracted from the REPET pipeline databaseand processed with a custom Perl script to calculate thetotal coverage of these types of repetitive DNA by omit-ting overlaps with transposable element or other repe-titive element annotations.Ethical approvalExperimental research followed appropriate guidelines forthe ethical treatment of research subjects. Human subjectand animal protocol approvals were not applicable.Data availabilityThe 454 transcript read data are available in the SequenceRead Archive (SRA) [119] at the NCBI (SRP003261,SRP003260, SRP001899). The assembled 454 transcripts(119,959) are in the Transcriptome Shotgun Assembly(TSA) database, and available through NCBI BioProjectpages for PRJNA51481 and PRJNA51483 [120].The accessions for the assembled 454 transcript dataare: HP542035 to HP552088 (testes), HP527956 toHP542034 (mixed antennae), HP509343 to HP527955(embryo), HP482918 to HP509342 (brain; ovary), HP473811 to HP482917 (larvae), HP459439 to HP473810 (abdomen),and HP552089 to HP579397 (ovary).The SOLID genomic read data are available in theSRA (SRX097020). The 454 genomic read data are alsoavailable (SRX006752, SRX000071).The Amel_4.5 assembly is available from NCBI underthe accession GCA_000002195.1.The following resources are available at BeeBase [121],a division of the Hymenoptera Genome Database [88]:genome browsers with gene annotations and supportingevidence alignments; BLAST databases with the Amel_4.5scaffold assembly, OGSv3.2, all input gene prediction setsand transcript contig assemblies; and a data downloadpage with fasta sequence files and gff for annotations.Peptides and their tandem mass spectra are available fromthe Apis mellifera PeptideAtlas [62].Additional filesAdditional file 1: Is a document containing Tables S1 through S7and Figure S1. Table S1. provides details for genomic sequencing runs,listed by SRA run number. Table S2. provides a comparison of New andPreviously Known OGSv3.2 genes based on relaxed mapping criteria.Table S3. provides a comparison of frequencies of canonical andnon-canonical splice sites in New and Previously Known genes. Table S4.provides the number of genes overlapping expressed sequencealignments for different transcript libraries. Table S5, provides counts ofnear-universal insect orthologous groups that are missing orthologs ineach species. Table S6. provides evidence and sampling options usedfor the three AUGUSTUS gene sets AU9, AU11, and AU12. Table S7.provides gene prediction accuracy of GeneID and SGP2 on an A. melliferaartificial contig. Figure S1. shows the proportions of differenttransposable element groups in the A. mellifera genome.Additional file 2: Is a spreadsheet showing results of evaluation of32 consensus gene sets generated with GLEAN.Additional file 3: Is a spreadsheet showing sources of biologicalevidence and other characteristics for each OGSv3.2 gene, and isthe data used in comparing Type I New, Type II New and PreviouslyKnown genes.Additional file 4: Is a spreadsheet showing new venom proteinsand known venom proteins with newly identified tryptic peptidesin OGSv3.2.Additional file 5: Is a spreadsheet providing results of manuallyannotating venom peptides.Additional file 6: Is a spreadsheet providing reciprocal best hitorthologs between A. mellifera and D. melanogaster and geneontology for D. melanogaster proteins.Additional file 7: Is a spreadsheet providing counts of InterProdomains in OGSv3.2 and OGSv1.0.Additional file 8: Is a spreadsheet providing an analysis of 1234OGSv3.2 genes that overlap interspersed repeats identified byREPET.AbbreviationsBAC: Bacterial artificial chromosome; CPLL: Combinatorial peptide ligandlibrary; DIRS: Dictyostelium intermediate repeat sequence; ENC: Effectivenumber of codons; EST: Expressed sequence tag; FTICR: Fourier transformion cyclotron resonance; GAG: GAG-Protein of retrotransposons; GO: GeneOntology; HMM: Hidden Markov model; HSP: High scoring pair; LARD: Largeretrotransposon derivative; LC: Liquid chromatography; LINE: Longinterspersed element; LTR: Long terminal repeat; MITE: MiniatiureElsik et al. BMC Genomics 2014, 15:86 Page 25 of 29http://www.biomedcentral.com/1471-2164/15/86inverted-repeat transposable element; MS/MS: Tandem mass spectrometry;nt: Nucleotides; NTR: New transcribed region; OGS: Official gene set;ORF: Open reading frame; POL: POL-polyprotein of retrotransposons;PWM: Position weight matrix; RBH: Reciprocal best hit; RT: Reversetranscriptase; SINE: Short interspersed element; SRA: Sequence Read Archive;SSR: Simple sequence repeat; TASE: Transposase; TE: Transposable element;TIR: Terminal inverted repeat; TRIM: Terminal repeat retrotransposon inminiature; UTR: Untranslated region.Competing interestsThe authors declare that they have no competing interests.Authors’ contributionsManuscript: CGE1, 2, KCW3, AKB2. Human Genome Sequencing CenterDirector: RAG3. Honey Bee Genome Phase II White Paper Authors: JDE10, MB4,CGE1, 2, RM20, HMR24, GER25, DBW30, CWW31. Sequence production: YW3(cDNA sequencing libraries at BCM-HGSC), IFN3 (cDNA sequencing libraries atBCM-HGSC), CLK3, MH3, VJ3, DMM3, HGSC production teams3, GER25 (cDNAsequencing at U. Illinois). cDNA: RM20 (Brain & Ovary, Testes, Mixed Antennae,Embryos, Larvae), GER25 (Queen Ovary, Nurse, Forager), MEH14 (Nurse andForager Brain), RSK17 (Nurse and Forager Brain), JDE10 (Abdomen), JD3(transcript analysis and assembly), DZ3 (transcript analysis). GenomeAssembly: KCW3, LZ3, HJ16, JD3. Genetic Map: GJH15, JMT28, JM19, OR26, ES21.Gene Annotation: CGE1, 2 (N-Scan, MAKER, GLEAN), AKB2 (N-Scan, MAKER),TDM23 (RefSeq, Gnomon), FC5 (SGP2, GeneID), RG5 (SGP2, GeneID), KJH13(AUGUSTUS), PK18 (Fgenesh++), VS27 (Fgenesh++), MS13 (AUGUSTUS). OGSEvaluation and Characterization: CGE1, 2, AKB2. GC Composition: CGE1, 2,AKB2, EE9, DG12, JTR2, 6. Orthology Assessment: RMW29, EMZ29. ProteinFunctional Analysis (GO and InterPro): CGE1, 2. Genomic Repetitive Elements:ES21, RFAM21. Peptide data and analysis: MVV7 (venom peptides), LJF11(PeptideAtlas), GD8 (venom peptides), BD8 (venom peptides), DCdG7 (venompeptides). Genome Browser, BLAST Website and Manual AnnotationSoftware: CGE1, 2, AKB2, CPC2, 6, MCM-T2, 22, JTR2, 6. All authors read andapproved the final manuscript.AcknowledgementsFunding for the project was provided by a grant to RG from the NationalHuman Genome Research Institute, National Institutes of Health (NHGRI, NIH)U54 HG003273. Contributions from members of the CGE lab were supportedby Agriculture and Food Research Initiative Competitive grant no. 2010-65205-20407 from the USDA National Institute of Food Agriculture. AKB wassupported by a Clare Luce Booth Fellowship at Georgetown University. Theauthors are grateful for the HGSC sequence production teams (Patil,S., Gub-bala,S., Aqrawi,P., Arias,F., Bess,C., Blankenburg,K.B., Brocchini,M., Buhay,C.,Challis,D., Chang,K., Chen,D., Coleman,P., Drummond,J., English,A., Evani,U.,Francisco,L., Fu,Q., Goodspeed,R., Haessly,T.H., Hale,W., Han,H., Holder,M., Hu,Y.,Jackson,L., Jakkamsetti,A., Jayaseelan,J.C., Kakkar,N., Kalra,D., Kandadi,H., Lee,S.,Li,H., Liu,Y., Macmil,S., Mandapat,C.M., Mata,R., Mathew,T., Matskevitch,T., Muni-dasa,M., Nagaswamy,U., Najjar,R., Nguyen,N., Niu,J., Opheim,D., Palculict,T.,Paul,S., Pellon,M., Perales,L., Pham,C., Pham,P., Pu,L.-L., Qi,S., Qu,J., Ren,Y., Ruth,R.T., Saada,N., Sabo,A., San Lucas,F., Sershen,C., Shafer,J., Shah,N., Shelton,R.,Song,X.-Z., Tabassum,N., Tang,L., Taylor,A., Taylor,M., Velamala,V., Wan,Z.,Wang,L., Wang,Y., Warren,J., Weissenberger,G., Wilczek-Boney,K.B., Yao,J., Yin,B.,Yu,J., Zhang,J., Zhang,L., Zhou,C., Zhu,D., Zhu,Y., and Zou,X.), and the input ofother members of the HGSC genome assembly team.Author details1Division of Animal Sciences, Division of Plant Sciences, and MU InformaticsInstitute, University of Missouri, Columbia, MO 65211, USA. 2Department ofBiology, Georgetown University, Washington, DC 20057, USA. 3HumanGenome Sequencing Center, Department of Molecular and Human Genetics,Baylor College of Medicine, MS BCM226, One Baylor Plaza, Houston, TX77030, USA. 4Institute of Evolutionary Genetics, Heinrich Heine UniversityDuesseldorf, Universitaetsstrasse 1, 40225 Duesseldorf, Germany. 5Center forGenomic Regulation, Universitat Pompeu Fabra, C/Dr. Aiguader 88, E-08003Barcelona, Catalonia, Spain. 6Division of Animal Sciences, University ofMissouri, Columbia, MO 65211, USA. 7Laboratory of Zoophysiology, GhentUniversity, Krijgslaan 281 S2, B-9000 Ghent, Belgium. 8Laboratory of ProteinBiochemistry and Biomolecular Engineering, Ghent University, K.L.Ledeganckstraat 35, B-9000 Ghent, Belgium. 9Department of Mental Health,Johns Hopkins University Bloomberg School of Public Health, Baltimore, MD21205-2103, USA. 10Bee Research Laboratory, BARC-E, USDA-AgriculturalResearch Service, Beltsville, MD 20705, USA. 11Department of Biochemistry &Molecular Biology, Centre for High-Throughput Biology, University of BritishColumbia, 2125 East Mall, Vancouver, BC, Canada. 12Department of Biologyand Biochemistry, University of Houston, Houston, TX 77204-5001, USA.13Ernst Moritz Arndt University Greifswald, Institute for Mathematics andComputer Science, Walther-Rathenau-Str. 47, 17487 Greifswald, Germany.14Department of Crop Sciences and Institute of Genomic Biology, Universityof Illinois at Urbana-Champaign, Urbana, IL 61801, USA. 15Department ofEntomology, Purdue University, 901 West State Street, West Lafayette, IN47907-2089, USA. 16Department of Obstetrics, Gynecology & ReproductiveSciences, University of Pittsburgh, MAGEE 0000, Pittsburgh, PA 15260, USA.17High-Performance Biological Computing (HPCBio), Roy J. CarverBiotechnology Center, University of Illinois at Urbana-Champaign, Urbana, IL61801, USA. 18Softberry Inc., 116 Radio Circle, Suite 400, Mount Kisco, NY10549, USA. 19Institute for Genomic Biology and Department ofBioengineering, University of Illinois at Urbana-Champaign, 1270 DCL,MC-278, 1304 W Springfield Ave, Urbana, IL 61801, USA. 20Research School ofBiology, The Australian National University, Canberra ACT 0200, Australia.21Institut für Zoologie, Molekulare Ökologie, Martin-Luther-UniversitätHalle-Wittenberg, Hoher Weg 4, D-06099 Halle (Saale), Germany. 22GenomicsDivision, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA.23National Center for Biotechnology Information, National Library ofMedicine, National Institutes of Health, Building 45, 8600 Rockville Pike,Bethesda, MD 20894, USA. 24Department of Entomology, University of Illinoisat Urbana-Champaign, Urbana, IL 61801, USA. 25Institute for GenomicBiology, Department of Entomology, Neuroscience Program, University ofIllinois at Urbana-Champaign, 1206 West Gregory Drive, Urbana, IL 61801,USA. 26Department of Biology, University of North Carolina at Greensboro,321 McIver Street, Greensboro, NC 27412, USA. 27Computer, Electrical andMathematical Sciences and Engineering Division, King Abdullah University ofScience and Technology (KAUST), Thuwal 23955-6900, Kingdom of SaudiArabia. 28Extension Field Operations, Clemson University, 120 McGinty Ct,Clemson, SC 29634, USA. 29University of Geneva and Swiss Institute ofBioinformatics, CMU, Michel-Servet 1, Geneva CH-1211, Switzerland.30Genformatic, 6301 Highland Hills Drive, Austin, TX 78731, USA.31Department of Entomology, Neuroscience Program, Program in Ecologyand Evolutionary Biology, University of Illinois at Urbana-Champaign, Urbana,IL 61801, USA.Received: 11 September 2013 Accepted: 27 January 2014Published: 30 January 2014References1. Honey Bee Genome Sequencing Consortium: Insights into social insects fromthe genome of the honeybee Apis mellifera. Nature 2006, 443:931–949.2. Zayed A, Robinson GE: Understanding the relationship between braingene expression and social behavior: lessons from the honey bee.Annu Rev Genet 2012, 46:591–615.3. Begna D, Han B, Feng M, Fang Y, Li J: Differential expressions of nuclearproteomes between honeybee (Apis mellifera L.) queen and worker larvae:a deep insight into caste pathway decisions. J Proteome Res 2012, 11:1317–1329.4. Cornman RS, Tarpy DR, Chen Y, Jeffreys L, Lopez D, Pettis JS, vanEngelsdorpD, Evans JD: Pathogen webs in collapsing honey bee colonies. PLoS ONE2012, 7:e43562.5. Foret S, Kucharski R, Pellegrini M, Feng S, Jacobsen SE, Robinson GE, Maleszka R:DNA methylation dynamics, metabolic fluxes, gene splicing, and alternativephenotypes in honey bees. Proc Natl Acad Sci USA 2012, 109:4968–4973.6. Li-Byarlay H, Li Y, Stroud H, Feng S, Newman TC, Kaneda M, Hou KK, WorleyKC, Elsik CG, Wickline SA, et al: RNA interference knockdown of DNAmethyltransferase 3 affects gene alternative splicing in the honey bee.Proc Nat Aca Sci USA 2013, 110:12750–12755.7. Calderone NW: Insect pollinated crops, insect pollinators and USagriculture: trend analysis of aggregate data for the period 1992-2009.PLos One 2012, 7:e37235.8. Gallai N, Salles JM, Settele J, Vaissière BE: Economic valuation of thevulnerability of world agriculture confronted with pollinator decline.Ecol Econ 2009, 68:810–821.9. Alaux C, Dantec C, Parrinello H, Le Conte Y: Nutrigenomics in honey bees:digital gene expression analysis of pollen's nutritive effects on healthyand varroa-parasitized bees. BMC Genomics 2011, 12:496.Elsik et al. BMC Genomics 2014, 15:86 Page 26 of 29http://www.biomedcentral.com/1471-2164/15/8610. Johnson RM, Mao W, Pollock HS, Niu G, Schuler MA, Berenbaum MR:Ecologically appropriate xenobiotics induce cytochrome P450s in ApisMellifera. PLoS ONE 2012, 7:e31051.11. Lander ES, Linton LM, Birren B, Nusbaum C, Zody MC, Baldwin J, Devon K,Dewar K, Doyle M, FitzHugh W, et al: Initial sequencing and analysis of thehuman genome. Nature 2001, 409:860–921.12. Venter JC, Adams MD, Myers EW, Li PW, Mural RJ, Sutton GG, Smith HO,Yandell M, Evans CA, Holt RA, et al: The sequence of the human genome.Science 2001, 291:1304–1351.13. Church DM, Goodstadt L, Hillier LW, Zody MC, Goldstein S, She X, Bult CJ,Agarwala R, Cherry JL, DiCuccio M, et al: Lineage-specific biology revealedby a finished genome assembly of the mouse. PLoS Biol 2009, 7:e1000112.14. International Human Genome Sequencing Consortium: Finishing theeuchromatic sequence of the human genome. Nature 2004, 431:931–945.15. Celniker SE, Wheeler DA, Kronmiller B, Carlson JW, Halpern A, Patel S, AdamsM, Champe M, Dugan SP, Frise E, et al: Finishing a whole-genomeshotgun: release 3 of the Drosophila melanogaster euchromatic genomesequence. Genome Biol 2002, 3:RESEARCH0079.16. C elegans Sequencing Consortium: Genome sequence of the nematode C.elegans: a platform for investigating biology. Science 1998,282:2012–2018.17. Elsik CG, Mackey AJ, Reese JT, Milshina NV, Roos DS, Weinstock GM:Creating a honey bee consensus gene set. Genome Biol 2007, 8:R13.18. Bonasio R, Zhang G, Ye C, Mutti NS, Fang X, Qin N, Donahue G, Yang P, LiQ, Li C, et al: Genomic comparison of the ants Camponotus floridanusand Harpegnathos saltator. Science 2010, 329:1068–1071.19. Nygaard S, Zhang G, Schiott M, Li C, Wurm Y, Hu H, Zhou J, Ji L, Qiu F,Rasmussen M, et al: The genome of the leaf-cutting ant Acromyrmexechinatior suggests key adaptations to advanced social life and fungusfarming. Genome Res 2011, 21:1339–1348.20. Smith CD, Zimin A, Holt C, Abouheif E, Benton R, Cash E, Croset V, CurrieCR, Elhaik E, Elsik CG, et al: Draft genome of the globally widespread andinvasive Argentine ant (Linepithema humile). Proc Natl Acad Sci USA 2011,108:5673–5678.21. Smith CR, Smith CD, Robertson HM, Helmkampf M, Zimin A, Yandell M, HoltC, Hu H, Abouheif E, Benton R, et al: Draft genome of the red harvesterant Pogonomyrmex barbatus. Proc Natl Acad Sci USA 2011, 108:5667–5672.22. Suen G, Teiling C, Li L, Holt C, Abouheif E, Bornberg-Bauer E, Bouffard P,Caldera EJ, Cash E, Cavanaugh A, et al: The genome sequence of theleaf-cutter ant Atta cephalotes reveals insights into its obligate symbioticlifestyle. PLoS Genet 2011, 7:e1002007.23. Werren JH, Richards S, Desjardins CA, Niehuis O, Gadau J, Colbourne JK,Beukeboom LW, Desplan C, Elsik CG, Grimmelikhuijzen CJ, et al: Functionaland evolutionary insights from the genomes of three parasitoid Nasoniaspecies. Science 2010, 327:343–348.24. Wurm Y, Wang J, Riba-Grognuz O, Corona M, Nygaard S, Hunt BG, IngramKK, Falquet L, Nipitwattanaphon M, Gotzek D, et al: The genome of the fireant Solenopsis invicta. Proc Natl Acad Sci USA 2011, 108:5679–5684.25. Florea L, Souvorov A, Kalbfleisch TS, Salzberg SL: Genome assembly has amajor impact on gene content: a comparison of annotation in two Bostaurus assemblies. PLoS One 2011, 6:e21400.26. Applied Biosystems: Whole genome sequencing of E. coli Bacteria.http://www3.appliedbiosystems.com/cms/groups/mcb_marketing/documents/generaldocuments/cms_057558.pdf.27. Wheeler DA, Srinivasan M, Egholm M, Shen Y, Chen L, McGuire A, He W, Chen YJ,Makhijani V, Roth GT, et al: The complete genome of an individual bymassively parallel DNA sequencing. Nature 2008, 452:872–876.28. Holt C, Yandell M: MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects.BMC Bioinformatics 2011, 12:491.29. Hunter S, Jones P, Mitchell A, Apweiler R, Attwood TK, Bateman A, BernardT, Binns D, Bork P, Burge S, et al: InterPro in 2011: new developments inthe family and domain prediction database. Nucleic Acids Res 2012,40:D306–D312.30. Parra G, Bradnam K, Ning Z, Keane T, Korf I: Assessing the genespace in draft genomes. Nucleic Acids Res 2009, 37:289–297.31. BeeBase Amel_4.5 JBrowse. http://hymenopteragenome.org:8080/Amel_4.5/jbrowse.32. Lee E, Helt GA, Reese JT, Munoz-Torres MC, Childers CP, Buels RM, Stein L,Holmes IH, Elsik CG, Lewis SE: Web Apollo: a web-based genomicannotation editing platform. Genome Biol 2013, 14:R93.33. BeeBase Amel_4.5 Annotation Tools. http://hymenopteragenome.org/beebase/?q=registration_apollo_amel.34. Wright F: The 'effective number of codons' used in a gene. Gene 1990,87:23–29.35. Van Vaerenbergh M, Debyser G, Devreese B, de Graaf DC: Exploringthe hidden honeybee (Apis mellifera) venom proteome byintegrating a combinatorial peptide ligand library approach withFTMS. J Proteomics 2014. In press.36. Lee E, Harris N, Gibson M, Chetty R, Lewis S: Apollo: a community resourcefor genome annotation editing. Bioinformatics 2009, 25:1836–1837.37. Ovchinnikov Y, Miroshnikov AI, Kudelin AB, Kostina MB, Boikov VA,Magazanik LG, Gotgil’F IM: Structure and presynaptic activity oftertiapin, a neurotoxin from honeybee venom. BioorganicheskayaKhimiya 1980, 6:359–365.38. Gmachl M, Kreil G: The precursors of the bee venom constituents apaminand MCD peptide are encoded by two genes in tandem which share thesame 3'-exon. J Biol Chem 1995, 270:12704–12708.39. Robertson HM, Wanner KW: The chemoreceptor superfamily in the honey bee,Apis mellifera: expansion of the odorant, but not gustatory, receptor family.Genome Res 2006, 16:1395–1403.40. Witte CP, Le QH, Bureau T, Kumar A: Terminal-repeat retrotransposons inminiature (TRIM) are involved in restructuring plant genomes. Proc NatlAcad Sci USA 2001, 98:13778–13783.41. Kalendar R, Vicient CM, Peleg O, Anamthawat-Jonsson K, Bolshoy A,Schulman AH: Large retrotransposon derivatives: abundant, conservedbut nonautonomous retroelements of barley and related genomes.Genetics 2004, 166:1437–1450.42. Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J:Repbase Update, a database of eukaryotic repetitive elements.Cytogenet Genome Res 2005, 110:462–467.43. Lampe DJ, Witherspoon DJ, Soto-Adames FN, Robertson HM: Recenthorizontal transfer of mellifera subfamily mariner transposons into insectlineages representing four different orders shows that selection acts onlyduring horizontal transfer. Mol Biol Evol 2003, 20:554–562.44. Feschotte C: Transposable elements and the evolution of regulatorynetworks. Nat Rev Genet 2008, 9:397–405.45. International Aphid Genomics Consortium: Genome sequence of the peaaphid Acyrthosiphon pisum. PLoS Biol 2010, 8:e1000313.46. Beye M, Gattermeier I, Hasselmann M, Gempe T, Schioett M, BainesJF, Schlipalius D, Mougel F, Emore C, Rueppell O, et al: Exceptionallyhigh levels of recombination across the honey bee genome.Genome Res 2006, 16:1339–1344.47. Kent CF, Minaei S, Harpur BA, Zayed A: Recombination is associatedwith the evolution of genome structure and worker behavior inhoney bees. Proc Natl Acad Sci USA 2012, 109:18012–18017.48. Simola DF, Wissler L, Donahue G, Waterhouse RM, Helmkampf M,Roux J, Nygaard S, Glastad K, Hagen DE, Viljakainen L, et al: Socialinsect genomes exhibit dramatic evolution in gene compositionand regulation while preserving regulatory features linked tosociality. Genome Res 2013, 23:1235–1247.49. Graveley BR, Brooks AN, Carlson JW, Duff MO, Landolin JM, Yang L,Artieri CG, van Baren MJ, Boley N, Booth BW, et al: Thedevelopmental transcriptome of Drosophila melanogaster.Nature 2011, 471:473–479.50. Ladoukakis E, Pereira V, Magny EG, Eyre-Walker A, Couso JP: Hundreds ofputatively functional small open reading frames in Drosophila.Genome Biol 2011, 12:R118.51. Kastenmayer JP, Ni L, Chu A, Kitchen LE, Au WC, Yang H, Carter CD, Wheeler D,Davis RW, Boeke JD, et al: Functional genomics of genes with small openreading frames (sORFs) in S. cerevisiae. Genome Res 2006, 16:365–373.52. i5K Consortium: The i5K initiative: advancing arthropod genomics forknowledge, human health, Agriculture, and the Environment. J Hered2013, 104:595–600.53. Genome 10 K Community of Scientists: Genome 10 K: a proposal to obtainwhole-genome sequence for 10,000 vertebrate species. J Hered 2009,100:659–674.54. Eckalbar WL, Hutchins ED, Markov GJ, Allen AN, Corneveaux JJ,Lindblad-Toh K, Di Palma F, Alfoldi J, Huentelman MJ, Kusumi K: Genomereannotation of the lizard Anolis carolinensis based on 14 adult andembryonic deep transcriptomes. BMC Genomics 2013, 14:49.55. ATLAS-Link. https://www.hgsc.bcm.edu/software/Atlas-Link.Elsik et al. BMC Genomics 2014, 15:86 Page 27 of 29http://www.biomedcentral.com/1471-2164/15/8656. Ament SA, Chan QW, Wheeler MM, Nixon SE, Johnson SP, Rodriguez-Zas SL,Foster LJ, Robinson GE: Mechanisms of stable lipid loss in a social insect.J Exp Biol 2011, 214:3808–3821.57. Chan MM, Choi SY, Chan QW, Li P, Guarna MM, Foster LJ: Proteome profileand lentiviral transduction of cultured honey bee (Apis mellifera L.) cells.Insect Mol Biol 2010, 19:653–658.58. Chan QW, Foster LJ: Changes in protein expression during honey beelarval development. Genome Biol 2008, 9:R156.59. Chan QW, Howes CG, Foster LJ: Quantitative comparison of castedifferences in honeybee hemolymph. Mol Cell Proteomics 2006,5:2252–2262.60. Chan QW, Melathopoulos AP, Pernal SF, Foster LJ: The innate immune andsystemic response in honey bees to a bacterial pathogen, Paenibacilluslarvae. BMC Genomics 2009, 10:387.61. Massa AN, Wanjugi H, Deal KR, O'Brien K, You FM, Maiti R, Chan AP, Gu YQ,Luo MC, Anderson OD, et al: Gene space dynamics during the evolutionof Aegilops tauschii, Brachypodium distachyon, Oryza sativa, andSorghum bicolor genomes. Mol Biol Evol 2011, 28:2537–2547.62. Chan QW, Parker R, Sun Z, Deutsch EW, Foster LJ: A honey bee (Apismellifera L.) PeptideAtlas crossing castes and tissues. BMC Genomics 2011,12:290.63. NCBI Apis mellifera Annotation Build 5.1. ftp://ftp.ncbi.nlm.nih.gov/genomes/Apis_mellifera.64. RepeatMasker Open-3.0. http://www.repeatmasker.org65. NCBI Eukaryotic genome annotation process. http://www.ncbi.nlm.nih.gov/genome/annotation_euk/process.66. Morgulis A, Gertz EM, Schaffer AA, Agarwala R: WindowMasker: window-based masker for sequenced genomes. Bioinformatics 2006, 22:134–141.67. Stanke M, Schoffmann O, Morgenstern B, Waack S: Gene prediction ineukaryotes with a generalized hidden Markov model that uses hintsfrom external sources. BMC Bioinformatics 2006, 7:62.68. Keller O, Odronitz F, Stanke M, Kollmar M, Waack S: Scipio: using proteinsequences to determine the precise exon/intron structures of genes andtheir orthologs in closely related species. BMC Bioinformatics 2008, 9:278.69. Kent WJ: BLAT–the BLAST-like alignment tool. Genome Res 2002,12:656–664.70. Solovyev V: Finding genes by computer: probabilistic and discriminativeapproaches. In Current Topics in Computational Biology. Edited by Jiang T, SmithT, Xu Y, Zhang M. Cambridge, MA: MIT Press; 2002:201–248.71. Solovyev V, Kosarev P, Seledsov I, Vorobyev D: Automatic annotation ofeukaryotic genes, pseudogenes and promoters. Genome Biol 2006,7(Suppl 1):S10 11–12.72. ReadMap: Aligning next generation sequences to the genome sequence.http://linux5.softberry.com/cgi-bin/berry/programs/ReadsMap.73. Parra G, Blanco E, Guigo R: GeneID in Drosophila. Genome Res 2000,10:511–515.74. Parra G, Agarwal P, Abril JF, Wiehe T, Fickett JW, Guigo R: Comparativegene prediction in human and mouse. Genome Res 2003, 13:108–117.75. van Baren MJ, Koebbe BC, Brent MR: Using N-SCAN or TWINSCAN topredict gene structures in genomic DNA sequences. Curr ProtocBioinformatics 2007, Chapter 4:Unit 4 8.76. Harris RS: Improved Pairwise Alignment of Genomic DNA. The PennsylvaniaState University; 2007.77. WU-BLAST. http://blast.wustl.edu.78. Slater GS, Birney E: Automated generation of heuristics for biologicalsequence comparison. BMC Bioinformatics 2005, 6:31.79. UniProt Consortium: Update on activities at the Universal ProteinResource (UniProt) in 2013. Nucleic Acids Res 2013, 41:D43–D47.80. Marygold SJ, Leyland PC, Seal RL, Goodman JL, Thurmond J, Strelets VB,Wilson RJ: FlyBase: improvements to the bibliography. Nucleic Acids Res2012, 41:D751–D757.81. Pearson WR: Rapid and sensitive sequence comparison with FASTP andFASTA. Methods Enzymol 1990, 183:63–98.82. Zhang Z, Schwartz S, Wagner L, Miller W: A greedy algorithm for aligningDNA sequences. J Comput Biol 2000, 7:203–214.83. Wu TD, Watanabe CK: GMAP: a genomic mapping and alignmentprogram for mRNA and EST sequences. Bioinformatics 2005, 21:1859–1875.84. Elhaik E, Graur D, Josic K, Landan G: Identifying compositionally homogeneousand nonhomogeneous domains within the human genome using a novelsegmentation algorithm. Nucleic Acids Res 2010, 38:e158.85. Elhaik E, Graur D: IsoPlotter+: a tool for studying the compositionalarchitecture of genomes. ISRN Bioinformatics 2013, 2013:725434.86. Olson SA: EMBOSS opens up sequence analysis. European molecularbiology open software suite. Brief Bioinform 2002, 3:87–91.87. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ:Gapped BLAST and PSI-BLAST: a new generation of protein databasesearch programs. Nucleic Acids Res 1997, 25:3389–3402.88. Munoz-Torres MC, Reese JT, Childers CP, Bennett AK, Sundaram JP, ChildsKL, Anzola JM, Milshina N, Elsik CG: Hymenoptera Genome Database:integrated community resources for insect species of the orderHymenoptera. Nucleic Acids Res 2011, 39:D658–D662.89. Waterhouse RM, Tegenfeldt F, Li J, Zdobnov EM, Kriventseva EV: OrthoDB:a hierarchical catalog of animal, fungal and bacterial orthologs.Nucleic Acids Res 2013, 41:D358–D365.90. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP,Dolinski K, Dwight SS, Eppig JT, et al: Gene ontology: tool for theunification of biology. The Gene ontology consortium. Nat Genet 2000,25:25–29.91. Castillo-Davis CI, Hartl DL: GeneMerge–post-genomic analysis, datamining, and hypothesis testing. Bioinformatics 2003, 19:891–892.92. Whitfield CW, Band MR, Bonaldo MF, Kumar CG, Liu L, Pardinas JR,Robertson HM, Soares MB, Robinson GE: Annotated expressed sequencetags and cDNA microarrays for studies of brain and behavior in thehoney bee. Genome Res 2002, 12:555–566.93. Zdobnov EM, Apweiler R: InterProScan–an integration platform for thesignature-recognition methods in InterPro. Bioinformatics 2001,17:847–848.94. Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, Pang N,Forslund K, Ceric G, Clements J, et al: The Pfam protein families database.Nucleic Acids Res 2012, 40:D290–D301.95. Haft DH, Selengut JD, White O: The TIGRFAMs database of proteinfamilies. Nucleic Acids Res 2003, 31:371–373.96. Letunic I, Doerks T, Bork P: SMART 7: recent updates to the proteindomain annotation resource. Nucleic Acids Res 2012, 40:D302–D305.97. Bru C, Courcelle E, Carrere S, Beausse Y, Dalmar S, Kahn D: The ProDomdatabase of protein domain families: more emphasis on 3D. Nucleic AcidsRes 2005, 33:D212–D215.98. Sigrist CJ, Cerutti L, de Castro E, Langendijk-Genevaux PS, Bulliard V, BairochA, Hulo N: PROSITE, a protein domain database for functionalcharacterization and annotation. Nucleic Acids Res 2010, 38:D161–D166.99. Wu CH, Nikolskaya A, Huang H, Yeh LS, Natale DA, Vinayaka CR, Hu ZZ,Mazumder R, Kumar S, Kourtesis P, et al: PIRSF: family classification systemat the protein information resource. Nucleic Acids Res 2004, 32:D112–D114.100. Lees J, Yeats C, Perkins J, Sillitoe I, Rentzsch R, Dessailly BH, Orengo C:Gene3D: a domain-based resource for comparative genomics, functionalannotation and protein network analysis. Nucleic Acids Res 2012,40:D465–D471.101. Wilson D, Pethica R, Zhou Y, Talbot C, Vogel C, Madera M, Chothia C,Gough J: SUPERFAMILY–sophisticated comparative genomics, datamining, visualization and phylogeny. Nucleic Acids Res 2009,37:D380–D386.102. Mi H, Dong Q, Muruganujan A, Gaudet P, Lewis S, Thomas PD: PANTHERversion 7: improved phylogenetic trees, orthologs and collaborationwith the Gene Ontology Consortium. Nucleic Acids Res 2010,38:D204–D210.103. Flutre T, Duprat E, Feuillet C, Quesneville H: Considering transposableelement diversification in de novo annotation approaches. PLoS One2011, 6:e16526.104. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local alignmentsearch tool. J Mol Biol 1990, 215:403–410.105. Quesneville H, Nouaud D, Anxolabehere D: Detection of new transposableelement families in Drosophila melanogaster and Anopheles gambiaegenomes. J Mol Evol 2003, 57(Suppl 1):S50–S59.106. Ellinghaus D, Kurtz S, Willhoeft U: LTRharvest, an efficient and flexiblesoftware for de novo detection of LTR retrotransposons.BMC Bioinformatics 2008, 9:18.107. Price AL, Jones NC, Pevzner PA: De novo identification of repeat familiesin large genomes. Bioinformatics 2005, 21(Suppl 1):i351–i358.108. Edgar RC, Myers EW: PILER: identification and classification of genomicrepeats. Bioinformatics 2005, 21(Suppl 1):i152–i158.Elsik et al. BMC Genomics 2014, 15:86 Page 28 of 29http://www.biomedcentral.com/1471-2164/15/86109. BLASTclust. http://www.ncbi.nlm.nih.gov/Web/Newsltr/Spring04/blastlab.html.110. HMMER 3.0. http://hmmer.org.111. Wicker T, Sabot F, Hua-Van A, Bennetzen JL, Capy P, Chalhoub B, Flavell A,Leroy P, Morgante M, Panaud O, et al: A unified classification system foreukaryotic transposable elements. Nat Rev Genet 2007, 8:973–982.112. Benson G: Tandem repeats finder: a program to analyze DNA sequences.Nucleic Acids Res 1999, 27:573–580.113. Kolpakov R, Bana G, Kucherov G: mreps: Efficient and flexible detection oftandem repeats in DNA. Nucleic Acids Res 2003, 31:3672–3678.114. Marchler-Bauer A, Anderson JB, Derbyshire MK, DeWeese-Scott C, GonzalesNR, Gwadz M, Hao L, He S, Hurwitz DI, Jackson JD, et al: CDD: a conserveddomain database for interactive domain family analysis. Nucleic Acids Res2007, 35:D237–D240.115. Kohany O, Gentles AJ, Hankus L, Jurka J: Annotation, submission andscreening of repetitive elements in repbase: repbase submitter andcensor. BMC Bioinformatics 2006, 7:474.116. Kapitonov VV, Tempel S, Jurka J: Simple and fast classification of non-LTRretrotransposons based on phylogeny of their RT domain protein sequences.Gene 2009, 448:207–213.117. Yuan YW, Wessler SR: The catalytic domain of all eukaryotic cut-and-pastetransposase superfamilies. Proc Natl Acad Sci USA 2011, 108:7884–7889.118. Kapitonov VV, Jurka J: A universal classification of eukaryotic transposableelements implemented in Repbase. Nat Rev Genet 2008, 9:411–412.119. Kodama Y, Shumway M, Leinonen R: The sequence read archive: explosivegrowth of sequencing data. Nucleic Acids Res 2012, 40:D54–D56.120. Apis mellifera 454 Transcriptome BioProjects. http://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA51481+OR+PRJNA51483.121. BeeBase. http://hymenopteragenome.org/beebase/.doi:10.1186/1471-2164-15-86Cite this article as: Elsik et al.: Finding the missing honey bee genes:lessons learned from a genome upgrade. BMC Genomics 2014 15:86.Submit your next manuscript to BioMed Centraland take full advantage of: • Convenient online submission• Thorough peer review• No space constraints or color figure charges• Immediate publication on acceptance• Inclusion in PubMed, CAS, Scopus and Google Scholar• Research which is freely available for redistributionSubmit your manuscript at www.biomedcentral.com/submitElsik et al. BMC Genomics 2014, 15:86 Page 29 of 29http://www.biomedcentral.com/1471-2164/15/86


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items