Finding the missing honey bee genes: lessons learned from a genome upgrade

Testing combined gene set approaches

We tested two approaches, MAKER2 [28] and GLEAN [17], for combining outputs from multiple gene prediction tools to create an Official Gene Set (OGS). Each approach has advantages and disadvantages. An advantage of MAKER2 is that given sufficient transcript data, it can predict multiple alternative splice forms per locus, while GLEAN computes only one gene model per locus. MAKER2 provides both a conservative biological sequence (transcript and homolog) alignment-based set of gene models that best agree with transcript and protein homolog alignments, and a set of ab initio gene predictions that do not overlap the biological sequence alignment-based set.

The final set of MAKER2 biological sequence alignment-based gene set contained 12,268 genes, encoding 12,575 transcripts, with an additional set of 31,047 ab initio gene predictions that did not have support from transcript or protein alignments. Without filtering the ab initio set for likely false positives, the choices of final gene sets generated with MAKER2 would have been 1) a conservative set of only 12,268 genes, still much lower than gene numbers found in other sequenced insect genomes or 2) a set of 43,315 genes that included ab initio gene models. Other hymenoptera genome projects using MAKER have selected ab initio gene models with matches to InterPro domains [29] to include in the final predicted gene set [20–22]. This approach would allow us to add 1,215 MAKER2 ab initio genes to the final predicted gene set, for a total of 13,303 genes. However, since only 8,602 of the 12,268 MAKER2 biological sequence alignment-based genes had InterPro matches, we were concerned that the InterPro approach to boost ab initio genes to the final gene set would miss some rapidly diverging genes.

An important difference between GLEAN and MAKER2 is the way ab initio predictions without biological sequence evidence are handled. GLEAN collects evidence for genes (ab initio gene predictions and biological sequence alignment-based), generates maximum likelihood estimates of accuracy and error rates for these signals for each gene prediction set, constructs consensus gene models that maximize the scores for the signal sites in each gene model, and labels each model with a probabilistic confidence score reflecting the underlying support for that gene model [17]. In addition to identifying consensus models supported by transcript and protein alignments, GLEAN identifies ab initio predictions with high computed probabilities, filtering out likely false positive ab initio predictions. This approach would not require filtering ab initio predictions using an InterPro search, and would potentially include high confidence ab initio predictions that are true, rapidly diverging genes without significant InterPro domain matches.

Since GLEAN weighs different sources of gene prediction evidence based on agreement among sources, without a priori information on the quality of individual input gene prediction sets, the choice of input datasets impacts the results. Evaluating GLEAN runs with different combinations of input sets allowed us to determine the optimum selection of datasets. We ran GLEAN 32 times with different combinations of input sets, and resulting GLEAN sets ranged in size from in 10,403 (with three input gene sets) to 17,724 genes (with seven input gene sets), (Additional file 2). In some cases, GLEAN predicted a larger number of genes than the MAKER2 set supported by either biological sequence alignment or InterPro. Although GLEAN generates only one transcript per gene locus, we decided to use GLEAN to create the OGS, because we prioritized GLEAN’s potential identification of a more comprehensive set of gene loci over MAKER2’s identification of 307 alternative transcripts. Since with GLEAN we did not need to filter ab initio genes based on InterPro, it might allow the genome-conservation-based N-SCAN dataset to provide support for other ab initio predictions even without detectable conservation with known InterPro domains. Another consideration was that assembly methods for Illumina RNASeq were not well established. We were concerned about the possibility of merging genes or missing introns due to potential genomic contamination in RNASeq experiments. GLEAN uses transcript alignments to support splice predictions, but does not merge or split genes based on transcript alignments, while MAKER2 creates gene models that best agree with transcript alignments, and works best with highly reliable transcript assemblies.

Selecting an official gene set

We evaluated the 32 GLEAN sets based on several criteria, including overlap with a conservative evidence-based set (RefSeq), transcript sequences, peptides and the CEGMA [30] conserved core set (Additional file 2). No single gene set was optimal with respect to all criteria. We chose to rank sets based on number of peptide matches, which would prioritize completeness of a protein-coding gene set rather than correctness of gene structure.

Assessing the new official gene set

Many split and merged gene models were noted when comparing the 32 GLEAN sets, including OGSv3.2, to the conservative RefSeq gene set (Additional file 2). Since it is difficult to computationally resolve disagreements among gene sets, the OGSv3.2 genes in question are provided as a separate track on the BeeBase Amel_4.5 genome browser [31]. Web Apollo [32] annotation tools allow registered BeeBase users to modify gene models, and can be used to manually correct split or merged genes [33].

To address the question of whether the increased gene number was due to splitting genes, we looked at the total number of coding nucleotides (nt) and the distribution of coding sequence lengths in OGSv1.0 and OGSv3.2. The ranges in coding sequence lengths were 24-53,649 and 75-70,263 nt for OGSv1.0 and OGSv3.2, respectively. Note that the gene prediction parameters did not allow coding sequences less than 75 nt for OGSv3.2, but a minimum coding sequence parameter was not applied in the generation of OGSv1.0; however there were only 6 OGSv1.0 genes less than 75 nt. The OGSv3.2 gene with the largest coding sequence overlapped a single RefSeq gene prediction (XM_623650.3), and was found to be homologous to genes predicted to code for titin. The total number of coding nucleotides increased from 16,484,776 in OGSv1.0 to 19,342,383 in OGSv3.2. The increase in coding nucleotides by only 2,857,607 (17.3%) compared to a 50.8% increase in the number of genes suggested that new genes were shorter and possibly a result of splitting OGSv1.0 gene models. However, the distribution of coding sequence length in the two gene sets suggested that splitting genes was not the major source of the increased gene number (Figure 3). Although OGSv3.2 possessed a large number of short coding sequences (< ~1.5 kb) compared to OGSv1.0, the larger 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 overlap with genes in OGSv1.0, and thus were not likely to be a result of splitting OGSv1.0 genes.

Biological evidence for the new official gene set

In assessing the gene set, we define biological evidence to include biological sequence alignments used in gene prediction pipelines (transcript, peptide, protein homolog) as well as matches to InterPro domains and alignment to other bee genomes. We define biological gene evidence as all of the datasets considered as biological evidence, except alignment to other bee genomes, since we did not use information about genes or expression for the other bee genomes. Some form of biological gene evidence supported the majority (14,084, 92%) of the OGSv3.2 genes. An additional 754 (4.9%) OGSv3.2 genes aligned to either the Apis florea or Bombus terrestris genome for a total of 14,836 (96.9%) of the OGSv3.2 genes with biological support. Requiring that transcripts be spliced reduces the number of supported genes to 13,264 (86.8%) or 14,661 (95.7%) if alignment of the gene to another bee genome is considered as support.

Identifying new genes due to improved assembly or improved gene prediction

In order to compare the old and new OGS (OGSv1.0 and OGSv3.2), we mapped OGSv3.2 coding sequences to the Amel_2.0 assembly because it was the assembly used to predict OGSv1.0 genes [1], which were later mapped to Amel_4.0. We chose to use Amel_2.0 in this step, because generation of Amel_4.0 and intermediate assemblies had involved rearranging and splitting scaffolds and incorporating unscaffolded contigs using an updated genetic map [1]. We wished to minimize mapping errors in our comparison of gene sets by using Amel_2.0 to determine which OGSv3.2 genes correspond to OGSv1.0 genes.

Of the 5,157 additional genes in OGSv3.2 compared to OGSv1.0, 4,735 genes did not overlap OGSv1.0 genes. The other 422 additional genes in OGSv3.2 were due to either splitting OGSv1.0 genes (discussed above) or to the annotation of additional recent paralogs in OGSv3.2, and are included in the total count of 10579 Previously Known genes. The 4735 genes that did not overlap with OGSv1.0 genes could be classified as 782 genes discovered due to the additional sequencing and reassembly of the bee genome for the Amel_4.5 assembly (Type I New genes; Table 5) and 3,953 genes detected by improved gene prediction, either through the use of new expressed sequence and protein data or improved gene prediction algorithms (Type II New genes; Table 5). We could map 405 additional Type I New genes to the Amel_2.0 assembly if we required only 50% of the gene to be covered (Table S2 in Additional file 1), rather than the more stringent 80% gene coverage reported in Table 5, consistent with the Amel_2.0 assembly being less continuous than Amel_4.5. This lack of contiguity and resulting fragmentation of genes in the Amel_2.0 assembly likely impaired the initial gene prediction efforts. While the frequency of genes with some form of biological support (transcript, peptide, protein homolog, InterPro domain, and/or alignment to another bee genome) was highest for Previously Known genes (99.7%), most of the Type I New genes (93.8%) and Type II New genes (89.9%) were also supported (Table 5).

Characteristics of new genes

We analyzed the data that went into the annotation of OGSv3.2 and evaluated which pieces of evidence contributed to the prediction of the genes to understand why they were missed in OGSv1.0. To determine whether genes that were not detected in OGSv1.0 have common characteristics that make them more challenging to predict, we compared them to the Previously Known genes. We evaluated features such as tissue expression specificity, coding feature length, GC content, overlap of protein homolog alignments on Amel_4.5 and the proportion of non-canonical splice sites. Additional file 3 provides sources of evidence for each OGSv3.2 gene.

The mean coding sequence lengths of Type I New (1,172 bp) and Type II New genes (331 bp) were shorter than that of Previously Known genes (1,623 bp) (P = 2.3 × 10^-12 and P < 2.2 × 10^-16 respectively) (Table 5). This difference may be due to a higher fraction of single coding exon genes among new genes. Thirteen percent of Type I New genes and 31% of Type II New genes contained one coding exon, while only 6.8% of Previously Known genes contained 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 sites was not significantly different between the Previously Known and Type II New genes (Table S3 in Additional file 1).

Type I New genes were found to reside in GC compositional domains with mean 26.4% GC, significantly lower than that of Previously Known genes (28.9%) (P = 2.188 × 10^-13), supporting improvement in the assembly of the high-AT-content regions. On the other hand, Type II New genes were found in GC compositional 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 mean GC content of the genome (32.7%)

The effective number of codons is an estimate of the deviation from equality of synonymous codon usage of all codons of a gene, and ranges from 20 (extreme bias where 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 effective 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 extreme the deviation from equal proportions of G + C and A + T nucleotides in coding sequences, the lower the potential diversity of synonymous codons. Consistent with their locations in less extreme GC compositional domains, Type II New genes had a significantly higher mean ENC than the Previously Known genes, 45.69 vs. 44.90 respectively (P = 3.923 × 10^-05) (Table 5).

Expression evidence for new genes

The majority of OGSv3.2 genes were supported by transcript 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 coverage by transcript dataset (Table S4 in Additional file 1) showed that the fraction of genes represented in a transcript 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 transcript 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 I New and Type II New genes were more likely to be narrowly expressed or overlap transcript alignments from only one tissue (P = 2.136 × 10^-14 and P < 2.2 × 10^-16,, respectively) 18.7% of Previously Known, 29.9% of Type I New 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 overlap transcript alignments from all tissues than either Type I New or Type II New genes (P < 2.2 × 10^-16 for both tests). 20.5% of Previously Known, 7.7% of Type I New and 2.4% of Type II New genes overlapped transcript alignments from all tissues. However, the fractions of narrowly expressed (18.7%) and broadly expressed (20.5%) genes were similar to each other for Previously Known genes.

Homolog alignment evidence for new genes

Results suggested that the use of dipteran proteins as the main source of protein homolog evidence for OGSv1.0 may have been a contributing factor to an incomplete gene set; 44.3% of the OGSv3.2 genes overlapped at least one protein homolog alignment, but only 7.9% overlapped an alignment of a D. melanogaster protein (Table 5). Differences among gene sets further supports the limited value of Dipteran proteins as a primary source of homolog evidence for A. mellifera gene prediction. Both Type I New and Type II New genes were less likely than Previously Known genes to overlap a D. melanogaster homolog alignment (P = 1.394 × 10^-07 and P < 2.2 × 10^-16, respectively). The effect was more drastic for Type II New genes; only 0.3% of Type II New genes overlapped a D. melanogaster homolog 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 less likely than Previously Known genes to overlap any homolog alignment from another sequenced genome (P < 2.2 × 10^-16 for both tests), but the difference was more extreme for Type II New genes, potentially implying that a greater number of the Type II New genes are specific to either A. mellifera, bees, or Hymenoptera than the Previously Known genes (Table 5). 59.8% of Previously Known genes overlapped a homolog alignment, while only 34.5% and 4.7% of Type I and Type II new genes, respectively, overlapped a homolog alignment (Table 5).

Genome alignment evidence for new genes

Over 80% of the OGSv3.2 genes aligned to both the Aflo_1.0 and Bter_1.0 genome assemblies (Table 5). A notable difference between Type II New genes and the other gene sets was the proportion of genes that were supported by genome conservation, but not other sources of evidence. Of the 3,555 Type II New genes supported by any evidence, 17.5% were not supported by other sources. On the other hand, only 2.9% of supported Type I New genes and 1% of supported Previously Known genes were supported by only genome conservation.

The remaining 20% of the OGSv3.2 genes have the potential to be Apis-specific (8% of genes that align to Aflo_1.0, but not Bter_1.0) or A. mellifera specific (12% of OGSv3.2 genes that do not align to Aflo_1.0) (Table 5).

Use of peptide evidence

Peptide data were used in four ways. First, peptide data were used in gene prediction by AUGUSTUS, the only program used in this study that accepts this source of gene evidence. Second, peptides were compared to all consensus gene prediction sets to identify the set with the highest representation of peptides (described above; Additional file 2). Third, peptide support was used in characterizing Previously Known, Type I and Type II New genes. Fourth, venom peptides were used to manually annotate venom genes associated with the sting of the honey bee.

Gene set comparison with all peptides

Peptide sequences aligned to a greater number of Previously 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 Type I and 2% of Type II New genes (Table 5 and Additional file 3).

Venom peptide analysis

Despite the lower proportion of new genes compared to Previously Known genes with peptide matches, analysis of venom peptides (a subset of the peptide data) showed that the Amel_4.5 assembly provides a significant contribution to venom proteome research. 705 unique venom peptides provided biological evidence for 102 genes (described in [35]). Searching the venom mass spectra against OGSv1.0 and OGSv3.2 showed that the improved assembly allowed detection of 21 additional peptides supporting 9 new venom protein identifications in OGSv3.2 (Additional file 4). Additional tryptic peptides were discovered for 7 venom proteins as a result of improved gene predictions (Additional file 4).

Manual annotation of genes supported by the venom peptides (Additional file 5) using Apollo [36] showed that most honey bee venom genes are fully (76.5%) or partially (19.6%) covered by transcriptome evidence. The putative biological functions of these genes are described elsewhere [35]. We discovered that one of the annotated genes (GB40695) codes for tertiapin. While the tertiapin peptide was already found in bee venom [37], no genomic or transcriptomic evidence had been described, an issue which is now solved as the genome improvement project supplies both a gene prediction (GB40695, NCBI Gene ID 100576769) and EST evidence (Genbank:HP466647.1). The gene is positioned on linkage group 12, next to the apamin and mast cell degranulating peptide venom genes. The three genes are tandemly arranged which suggest their origin by gene duplication via unequal crossing over, and may also point to a joint control of transcription [38].

Gene ontology analysis

Gene Ontology analysis (Additional file 6) showed enrichment of some specific functions in new genes. Type I New genes were enriched for the biological processes “apoptosis” (corrected e-score = .0016), “neurotransmitter release” (corrected e-score = .0025), and “secondary active organic cation transmembrane transporter” (corrected e-score < .0001). The detection of these Type I New genes due to new assembly data is likely due to their location in low GC content regions, which were underrepresented in the older assembly; this is in agreement with the lower mean GC content of GC compositional domains containing Type I New genes. Type II new genes were enriched for the molecular function “nuclease activity” (corrected e-score = .011). Rapidly evolving and single exon genes are more difficult to annotate automatically, so it is perhaps not surprising that we found significant differences in some activities.

InterPro analysis

Comparison of InterPro protein domains found in OGSv3.2 and OGSv1.0 proteins showed that 269 of the InterPro domains present in OGSv3.2 were not present in OGSv1.0 (Additional file 7). There were also significant differences in the fraction of genes with the domains IPR004117 (Olfactory receptor, Drosophila, p < .02) and IPR001888 (Transposase, type 1; p < .01). There were 141 and 92 genes with olfactory receptor domains in OGSv3.2 and OGSv1.0, respectively. Olfactory receptor genes contain a single olfactory receptor domain, so the domain count corresponds to the gene count. The additional olfactory receptor genes were found among both the Type I and Type II New genes. The family of olfactory receptor genes is expanded in hymenopteran insects, with rapidly diverging family members that are often arranged in tandem arrays [20, 21, 39]. Many of the 166 known A. mellifera olfactory receptor genes were identified by manual annotation of the previous assembly, because they were not found in OGSv1.0 [20, 21, 39]. Improved computational prediction of olfactory receptors in the new assembly was likely due to increased assembly continuity which would improve identification of tandemly repeated genes, and new sources of biological sequence alignment evidence (e.g. transcript and protein homolog sequences) which improve identification of rapidly diverging genes.

There were 18 and 3 genes with transposase type 1 domains in OGSv3.2 and OGSv1.0, respectively. Thirteen of the additional transposase type 1domains were in Type 1 New genes. Transposases are associated with transposable elements, which are a class of interspersed repeats. Since the new assembly information filled in gaps, and gaps in de novo assemblies are found where the assembly process fails to find an unambiguous path (e.g. repetitive loci), we were not surprised to identify InterPro domains associated with repetitive sequences among the Type I New genes. Predicted genes with matches to transposable element domains are often removed from gene sets. We further investigated OGSv3.2 for presence of interspersed repeats, and decided not to remove the small number of genes associated with transposable element domains, because evidence suggested that some were real host genes (see Genome-wide repeat analysis).

RNAseq data

We sequenced samples from a number of tissues using the 454 Titanium technology. Tissues from ovary, testes, mixed antennae (worker, drones, and queens), larvae, mixed embryos, abdomen, and a combined library from brain and ovary samples were sequenced.

The testes cDNA library was prepared with the Clontech “full-length” amplification protocol. A gel cut of 400 to 800 bp fragments was combined with nebulized products from the larger cDNA fragments to generate fragments of an optimum sized sequencing library for the 454 Titanium platform. Other libraries were prepared by isolating total RNA using Trizol and RNeasy columns followed by mRNA isolation using the Qiagen kit and cDNA genration using the Invitrogen Superscript kit #11917-010 with random primers. A gel cut of 400 to 800 bp fragments was combined with nebulized products from the larger cDNA fragments to generate fragments of an optimum sized sequencing library for the 454 Titanium platform.

Peptide data

Summary of gene prediction sets

Genes were predicted using a variety of methods including the NCBI RefSeq and Gnomon pipelines, AUGUSTUS, SGP2, GeneID, Fgenesh++ and N-SCAN. Unless stated otherwise below, gene prediction pipelines used the masked assembly available from NCBI [63], which was generated using RepeatMasker [64]. New transcriptome data was used either directly as evidence or in the generation of training sets for the predictors. The Augustus analysis also incorporated available peptide data, and the N-SCAN analysis leveraged nucleotide sequence conservation between the A. mellifera genome and the other bee genomes. SGP2 leveraged conservation between the A. mellifera genome and the previously published Nasonia genomes [23] based on translated alignments.

NCBI RefSeq and Gnomon

The final RefSeq annotation included models that were completely or partially supported by alignment evidence. The pipeline had some provisions to predict models that were disrupted by frameshifts or stop codons in the assembly, some of which were assigned protein accessions (XP_ prefixes) and named with the prefix ‘LOW QUALITY PROTEIN’, whereas others were conservatively classified as pseudogenes and not assigned protein accessions. The RefSeq annotation is available from NCBI's genomes FTP site as NCBI build 5.1 [63].

AUGUSTUS

AUGUSTUS can be used as an ab initio gene prediction tool, but can also integrate extrinsic evidence from various sources [67].

Generation of training gene structures and training AUGUSTUS

We used Scipio [68] to generate gene structures on the Amel_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 AUGUSTUS parameters for A. mellifera, and constructed UTR models from 78,274 A. mellifera ESTs from GenBank, which were used to train AUGUSTUS UTR parameters. We then predicted genes in the Amel_4.5 assembly. Individual RNAseq reads from 454 sequencing (described above) were mapped against predicted transcripts, and fully covered transcripts (2,151) were selected as training genes for optimizing a final AUGUSTUS parameter set.

Gene prediction with AUGUSTUS

We generated three gene sets using AUGUSTUS: a gene set with extrinsic evidence from ESTs and RNAseq data (AU9), a particularly inclusive gene set that contained many alternative transcripts for peptide identification (AU11), and a gene set with extrinsic evidence from ESTs, RNAseq data and peptide data (AU12) (Table S6 in Additional file 1).

We created extrinsic evidence “hints” for protein coding genes and transcripts from the A. mellifera ESTs and from 454 transcriptome libraries. Genes were predicted with AUGUSTUS, allowing the prediction of alternative transcripts and allowing the splice site AT-AC (in addition to GT-AG and GC-AG) in case of supporting extrinsic 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.fa

For the purpose of peptide identification, alternative transcripts were extensively sampled with AUGUSTUS, resulting in the gene set AU11. AU11 genes were predicted issuing the following command: augustus –UTR = on –print_utr = on –hintsfile = all_but_no_peptide.hints –extrinsicCfgFile = extrinsic.M.RM.E.W.cfg –exonnames = on –codingseq = on –species = honeybee1 –alternatives-from-evidence = true –alternatives-from-sampling = true –sample = 100 –minexonintronprob = 0.1 –minmeanexonintronprob = 0.4 –maxtracks = -1 –allow_hinted_splicesites = atac genome.fa

We generated hints from peptide sequences by mapping the peptides against the protein set of AU11 and against a six-frame translation of the genome using BLAT [69] in a non-redundant way (i.e. parts of the six-frame translation that were included in the protein set AU11 and redundant parts of the AU11 protein set that were removed). A final AUGUSTUS gene set AU12 was generated using all available extrinsic evidence and the following 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 for A. mellifera. We used Illumina transcriptome data from A. mellifera forager and nurse brains available from the SRA (SRP003528), which were a total of 181.8 M spots, 18.3G bases of 100 bp single end reads generated on an Illumina Genome Analyzer II in 2010. We mapped the Illumina RNASeq using the ReadsMap program [72], which provided intron position information to Fgenesh for building sample-specific gene models. Predictions were made based on individuals used in Illumina RNASeq libraries (five nurses and five foragers), and then were combined into forager and nurse group predictions, and redundant genes (those having coinciding coding sequences) were removed from each set. Finally, the forager and nurse predictions were combined into a final set of predictions, and redundant genes (those having coinciding coding sequences) were removed from this combined set to produce the final Fgenesh++ prediction set.

GeneID

GeneID is an ab initio gene prediction program used to find potential protein-coding genes in anonymous genomic sequences. The training of GeneID to obtain a parameter file for A. mellifera was based on the method described to obtain a Drosophila melanogaster GeneID parameter file [73]. Training was performed in a “semi-automated” manner by employing a recently developed GeneID training tool that computes position weight matrices (PWMs) or Markov models of order 1 for splice sites and start codons, and derives a model of coding DNA, which, in this case, is a Markov model of order 5. Furthermore, once a preliminary species-specific matrix is obtained it is further optimized by adjusting two internal matrix parameters: -the cutoff of the scores of the predicted exons (eWF) and the ratio of signal to coding statistics information to be used (oWF).

The initial A. mellifera training set was comprised of the 2,151 gene models used to train AUGUSTUS (described above). Of these gene models 80% (1,720) were used to train GeneID while the remaining 20% (431) were set aside to test the accuracy of the newly developed matrix. The 1,720 A. mellifera protein-coding gene models included 7,913 canonical donor splice sites/7,939 canonical acceptor sites and 1,720 start codons. The start codons were used to compute PWMs while the donor and acceptors were used to derive Markov matrices of order 1. Given the large number of sequences we also had enough coding (2,338,800) and non-coding (5,561,154) nucleotides to derive a Markov of order 5 for the coding potential. We tested accuracy of the GeneID A. mellifera parameter file on an artificial contig consisting of the 431 evaluation-set concatenated gene models with 800 nucleotides of intervening sequence between each of the genes (Table S7 in Additional file 1). We then used the GeneID parameter files to predict genes on an assembly consisting of Amel_4.5 sixteen chromosomes and 5,304 "unplaced" scaffolds files that had repeat sequences masked using Repeatmasker. GeneID predicted 24,554 protein-coding genes.

SGP2

SGP2 is a syntenic gene prediction tool that combines ab initio gene prediction (GeneID) with TBLASTX searches between two or more genome sequences to provide both sensitive and specific gene predictions, and it tends to improve GeneID’s performance, especially by reducing the number of false-positive predictions. SGP2 requires one or more reference genomes to which the target genome (in this case A. mellifera) is compared. We decided to use the genomes of Nasonia vitripennis, N. giraulti and N. longicornis[23] as references to develop our A. mellifera parameter file for SGP2 because the genus Nasonia is at an appropriate evolutionary distance from Apis such that mostly the coding regions of the genes, not the introns or intergenic regions, are significantly conserved between these two genomes.

Obtaining the SGP2 A. mellifera-specific parameter file was based on the methodology described by Parra et al[74] used to obtain a human SGP2 parameter file using mouse homology evidence. The starting point to obtain a parameter file for SGP2 was the previously described GeneID A. mellifera matrix. The GeneID-derived SGP2 parameter file was optimized on an artificial contig comprising the same concatenated 1,720 sequences used to train GeneID, with 800 nucleotides between each of the gene models. We optimized the SGP2 matrix by modifying not only the eWF internal parameter (as previously for the GeneID parameter file) but also two SGP2-specific internal parameters (“NO_SCORE” and “HSP_factor”). The “NO_SCORE” parameter provides a penalty for no overlap between TBLASTX-derived HSPs (High-Scoring Pairs) and GeneID ab initio predictions in the same region whereas the “HSP_factor” parameter reduces the score assigned to the HSPs in order to maximize the prediction accuracy. We evaluated the newly developed SGP2 parameter file on the same artificial contig consisting of the 431 concatenated gene models with 800 nucleotides of intervening sequence between each of the genes used to evaluate the GeneID bee matrix (Table 1). We then used the SGP2 parameter file2 to predict genes on the same assembly file used for GeneID, and generated 20,179 predictions.

N-SCAN

We used the N-SCAN package [75] to leverage conservation between the A. mellifera genome and genomes of two other bee species, A. florea (Aflo_1.0 ) and Bombus terrestris (Bter_1.0). We first masked Amel_4.5 for simple sequence repeats using RepeatMasker [64]. We ran LASTZ [76] using default parameters with Amel_4.5 as the target genome, and either Apis florea (Aflo_1.0) or Bombus terrestris (Bter_1.0) as the informant genome. We then used iParameterEstimation to generate both an Amel_4.5-Aflo_1.0 specific parameter set as well as an Amel_4.5-Bter_1.0 parameter set using the training set described above for AUGUSTUS gene prediction, including UTR features. Finally, we ran N-SCAN using each of the A. mellilfera specific parameter sets with the respective LASTZ informant alignments to produce two N-SCAN prediction sets, one set based on Aflo_1.0 as the informant genome and the other set based on Bter_1.0 as the informant genome.

Input data for combined gene sets

We used MAKER2 and GLEAN to generate combined gene sets. The MAKER2 and GLEAN analyses used the same set of input data. Both analyses combined the gene predictions described above (NCBI, AUGUSTUS, Fgenesh++ with RNAseq, N-SCAN using A. florea as an informant genome, GeneID and SGP2) with transcript and protein homolog alignments. Transcript data included the new 454 transcriptome data described above, A. mellifera ESTs from GenBank and Illumina nurse and forager reads downloaded 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 the option "--butterfly-search" for more sensitive splice junction detection, and then generated predicted transcripts for each set of pooled data using Cufflinks version 1.0.3 with default parameters. The 454 reads were assembled into contigs de novo using Newbler (2.3-PreRelease-9/14/2009) with the cDNA option. We aligned 454 contigs and ESTs to Amel_4.5 using MAKER2 v2.15, which uses WU-BLAST [77] and Exonerate est2genome [78], with minimum 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], Harpegnathos saltator (OGSv3.3) [18], Linepithema humile (OGSv1.1) [20], Pogonomyrmex barbatus (OGSv1.1) [21]). Proteins in the SwissProt dataset annotated as transposable elements were removed prior to alignment. We aligned protein sequences to Amel_4.5 using Exonerate protein2genome with a minimum 60% percent identity and 60% alignment coverage.

MAKER2

To create a combined gene set, we ran MAKER2 v2.15 using parameters min_contig = 1000 and pred_gff, which allowed us to provide as input the gene prediction sets and alignments described above, instead of generating new evidence tracks within MAKER2.

GLEAN

We ran GLEAN [17] 32 times to create consensus gene sets using different combinations of the gene prediction sets described above. All of the GLEAN runs included the transcript and protein homolog alignments described above, and required a minimum coding sequence length of 75 nt.

Selecting the new official gene Set

We evaluated the 32 GLEAN sets based on several criteria including overlap with a conservative evidence-based set (RefSeq), transcript sequences, peptides and the CEGMA conserved core set [30] (Additional file 2). NCBI’s RefSeq pipeline has been considered reliable and relatively conservative, so we performed overlap analyses to determine how many of the RefSeq models were captured in the GLEAN consensus sets. We used FASTA [81] to align GLEAN coding sequences with RefSeq coding sequences where a perfect alignment required 99% identity over the entire lengths of both sequences. We used FASTA to align peptides to GLEAN proteins, with E-value and scoring matrix optimized for short exact matches (E-value .01 and MD10 scoring matrix). We parsed peptide alignments to count those with 100% identity and 100% alignment coverage (Additional file 2).

In addition to alignments between GLEAN predictions and RefSeq genes and peptides, we considered numbers of gene models, total numbers of coding nucleotides, average coding sequence lengths, and numbers of splits and merges compared to RefSeq. We selected the GLEAN31 set (Additional file 2) to be the official gene set, OGSv3.2. It was the set generated using gene predictions from NCBI Gnomon, AUGUSTUS, Fgenesh++, N-SCAN using A. florea as an informant genome and GeneID and SGP2 combined as a single set as well as the transcript and protein homolog alignments.

Adding UTRs to OGSv3.2 gene models

Although we did not use MAKER2 to generate the official gene set, we did use MAKER2 (v2.15) [28] to add UTR to the GLEAN coding exons since GLEAN does not produce annotations with UTR features. OGSv3.2 coding exons were input as “pred_gff” and transcriptome evidence (including the 454 transcripts, Illumina contigs and A. mellifera ESTs described above) was input as “est”. MAKER2 aligned the EST evidence to Amel_4.5 using WU-BLAST and Exonerate est2genome then used overlapping EST evidence to extend OGSv3.2 models to include UTR features when possible. Because MAKER2 sometimes modified the coding exon coordinates, we processed the MAKER2 output and retained UTRs only when the coding exon structure was unchanged. Out of a total of 15,314 OGSv3.2 genes, UTR were added to 7,514 genes (49%).

Identifying of new and previously known OGSv3.2 genes

In order to compare the newest gene set (OGSv3.2) with the first official gene set (OGSv1.0), we mapped OGSv3.2 coding sequences to Amel_2.0, which was the assembly on which OGSv1.0 was generated [1]. First, we used MegaBLAST [82] to identify scaffold/gene matches with 95% identity and E-value < 1 × 10^-20. We then aligned coding sequences to matching scaffolds using GMAP [83], and parsed the output to create two sets of splice-modeled alignments, both requiring 95% identity. One set was based on a relaxed coverage criterion, requiring that the alignment cover at least 50% of the coding sequence. The other set was based on a stringent coverage criterion, requiring that the alignment cover 80% of the coding sequence. Results of further analyses for the relaxed mapping set are provided in the supplemental materials, but we discuss only results for the stringent mapping.

On the basis of mapping OGSv3.2 to Amel_2.0 and overlap between OGSv3.2 and OGSv1.0 gene models on the Amel_2.0 assembly, we divided 15,314 OGSv3.2 genes into three sub-sets (Table 6). We deemed any OGSv3.2 gene that did not align to the Amel_2.0 assembly a “Type I New” gene. The additional sequencing and reassembly of the genome for the Amel_4.5 assembly likely allowed the detection of these genes. “Type II New” genes were those that did align to the Amel_2.0 assembly, but whose coordinates did not overlap an OGSv1.0 gene by a single coding base pair on the same strand. Additional expressed sequence and protein homolog evidence as well as improvements to gene prediction algorithms likely contributed to the detection of these genes. Finally, any OGSv3.2 gene that both aligned to the Amel_2.0 assembly and overlapped an OGSv1.0 gene was deemed a “Previously Known” gene.

Coding sequence length analysis

The total length of the coding sequence was calculated for each gene and the means for all genes and each gene sub-set were calculated (Table 6). We tested the null hypotheses that the mean coding sequence lengths of Type I and Type II New genes and Previously Known genes were equal.

Splice site and single versus multiple coding exon gene analysis

We assessed the genomic sequence of the two intronic base pairs adjacent to each coding sequence exon-intron splice site to determine whether they corresponded to the canonical …]5’-GT/AG-3’[… splice site sequence. We considered only splice sites supported by matching intron coordinates of spliced transcript alignment evidence. We used chi-square tests with one degree of freedom to compare the frequencies of single coding exon genes and non-canonical splice sites in Type I or Type II New genes with Previously Known genes.

OGSv3.2 gene location relative to GC compositional domains

We used IsoPlotter, a recursive segmentation algorithm [84, 85], to partition the A. mellifera genome into GC compositionally homogeneous domains, contiguous regions of the genome with similar GC content (percent G + C nucleotides). We determined the GC content of the GC compositional domain for each OGSv3.2 gene. In cases where a gene spanned multiple GC compositional domains, we calculated the average GC content of the GC compositional domains, and weighted it according to the fraction of the gene's length that occurs in each GC compositional domain. We computed the weighted percent GC based only on non-coding nucleotides of the compositional domains, because we did not wish to include effects of codon bias.

We compared the distribution of Type I and Type II New genes with respect to GC content to that of Previously Known genes. IsoPlotter cannot segment scaffolds less than 10 kb into compositional domains, so the 90 genes residing in these short scaffolds were not considered. We tested the null hypotheses that the means of the weighted GC compositional domain contents for genes in the Type I and Type II New gene sets were equal to the mean of the Previously Known genes.

Effective number of codons analysis

Using the chips program within the EMBOSS package [86], we calculated the effective number of codons separately for each OGSv3.2 gene. We tested the null hypotheses that the mean effective number of codons values for the Type I and Type II New genes were equal to the mean of the Previously Known genes.

Tissue expression analysis

We determined the number of OGSv3.2 gene overlaps to transcript alignments that were used in creating the GLEAN consensus gene sets (454 reads, Illumina contigs and A. mellifera ESTs from GenBank, described above). We relied on splice signals to determine the directionality of a transcript read. We tallied spliced and unspliced alignments separately, since we could be confident when directionality of a spliced alignment agreed with a gene prediction, but could not be confident about unspliced alignments. For spliced transcript alignments, if the transcript was on the opposite strand from the gene then it was discarded from further analysis. For transcripts on the same strand or transcripts that were un-spliced, in which case directionality could not be determined, a coordinate overlap of at least one coding base pair was required for a gene to count as overlapping with a transcript alignment. We counted the number of transcript data sets in which each OGSv3.2 gene was found to have an overlapping alignment (Table 6.) We performed chi-square tests with one degree of freedom to compare the frequencies of spliced transcript overlap in the Type I or Type II New gene sets with the Previously Known gene set.

Of genes that overlapped spliced transcript alignments, we identified genes that were narrowly expressed and genes that were broadly expressed on the basis of overlap to the four single-tissue libraries (brain [combined Illumina forager and nurse brain libraries] and 454 libraries of mixed antennae, ovary and testes). (Foragers and nurses are worker honey bees that specialize on collecting food and feeding brood, respectively.) Genes were deemed narrowly expressed if they overlapped at least one transcript alignment in only one of the four tissues and broadly expressed if they overlapped at least one transcript alignment from all four tissues. We performed chi-square tests with one degree of freedom to compare the frequencies of narrowly expressed genes and broadly expressed genes in the Type I or Type II New gene sets with the Previously Known gene set.

Homolog analysis

We determined the number of protein alignments overlapping OGSv3.2 for protein data sets that were used in creating the GLEAN consensus gene sets. We required overlap of at least one coding base pair on the same strand to deem a gene overlapping with a protein homolog alignment. We performed chi-square tests with one degree of freedom to compare the frequency of the existence of overlaps to homolog alignments for Type I or Type II New genes with that of Previously Known genes.

TBLASTN alignment of OGSv3.2 to A. florea and B. terrestris genomes

We aligned OGSv3.2 protein sequences to the A. florea (Aflo_1.0) and B. terrestris (Bter_1.0) genome assemblies using TBLASTN [87] with an E-value criteria of 1 × 10^-06. We did not use information about predicted genes or expression in A. florea or B. terrestris.

Statistical methods to test for differences in means between new and previously known genes

To test the null hypothesis that the mean of a particular variable for Type I or Type II New genes was equal to the mean of the Previously Known genes, we used both a non-parametric Kolmogorov-Smirnov test and a Welch t-test with the correction for non-homogeneity of variances. Testing the hypotheses in this way avoids assuming these data were normally distributed or had equal variances. To be conservative, we report the least significant P-value for each test.

GO analysis

We used FASTA [81] with an E-value threshold of 1 × 10^-6 to compute reciprocal alignments between OGSv3.2 proteins and a D. melanogaster protein set consisting of the longest protein isoform of each gene (annotation version r5.42). We identified reciprocal best hits (RBH) and transferred Gene Ontology (GO) [90] annotations from the D. melanogaster protein to the A. mellifera protein for each RBH pair, using the GO annotation file available at FlyBase [80]. We used GeneMerge [91] to test for enrichment of GO terms in OGSv3.2 protein datasets, testing for each of the three GO ontologies (Molecular Function, Biological Process and Cellular Component) separately. Several tests were performed with either the entire set of Gene Ontology terms, the generic GO slim set, or the GO slim set developed for A. mellifera by Whitfield et al. [92]. Population and test gene datasets for a particular GO ontology included only OGSv3.2 genes with GO annotations for that ontology. Population datasets consisted of all OGSv3.2 genes with GO annotations. Test datasets were 1) all new genes based on stringent mapping criteria, 2) Type I New genes using stringent criteria, 3) Type II New genes based on stringent mapping criteria.

InterPro analysis

We used InterProScan [93] to compare OGSv3.2 and OGSv1.0 proteins with the following InterPro [29] protein domain 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 least one InterPro domain were 9,479 and 8,552 for OGSv3.2 and OGSv1.0, respectively. For each InterPro domain identified in the combined datasets, we determined the number of proteins containing that domain within OGSv3.2 and OGSv1.0. Then for each InterPro domain, we used 2 × 2 chi-square tests with Yates correction and one degree of freedom to determine whether the frequencies of proteins containing that domain differed between the OGSv3.2 and OGSv1.0 InterPro-annotated sets (total 9,479 and 8,552 proteins, respectively).