Gene discovery in the horned beetle Onthophagus taurus

Abstract

Background: Horned beetles, in particular in the genus Onthophagus, are important models for studies on sexual selection, biological radiations, the origin of novel traits, developmental plasticity, biocontrol, conservation, and forensic biology. Despite their growing prominence as models for studying both basic and applied questions in biology, little genomic or transcriptomic data are available for this genus. We used massively parallel pyrosequencing (Roche 454-FLX platform) to produce a comprehensive EST dataset for the horned beetle Onthophagus taurus. To maximize sequence diversity, we pooled RNA extracted from a normalized library encompassing diverse developmental stages and both sexes.

Results: We used 454 pyrosequencing to sequence ESTs from all post-embryonic stages of O. taurus. Approximately 1.36 million reads assembled into 50,080 non-redundant sequences encompassing a total of 26.5 Mbp. The non-redundant sequences match over half of the genes in Tribolium castaneum, the most closely related species with a sequenced genome. Analyses of Gene Ontology annotations and biochemical pathways indicate that the O. taurus sequences reflect a wide and representative sampling of biological functions and biochemical processes. An analysis of sequence polymorphisms revealed that SNP frequency was negatively related to overall expression level and the number of tissue types in which a given gene is expressed. The most variable genes were enriched for a limited number of GO annotations whereas the least variable genes were enriched for a wide range of GO terms directly related to fitness.

Conclusions: This study provides the first large-scale EST database for horned beetles, a much-needed resource for advancing the study of these organisms. Furthermore, we identified instances of gene duplications and alternative splicing, useful for future study of gene regulation, and a large number of SNP markers that could be used in population-genetic studies of O. taurus and possibly other horned beetles.

Figure 1

Sequence matches to NR protein database and Tribolium genome and proteins. Venn diagram showing the number of Onthophagus contigs and singletons (in parenthesis) with sequence matches against the NCBI NR database [52], Tribolium genome sequence [92] and Tribolium annotated proteins [92]. The number of sequence matches at E-value cut-offs of 1 × 10^-5, 1 × 10 ^-20 and 1 × 10^-50 are shown in black, red and blue, respectively.

Figure 2

Taxonomic distribution of sequence matches. Phylogenetic tree showing the number of O. taurus non-redundant sequences assigned to branches. The MEGAN algorithm used in this analysis assigns each sequence to the lowest common ancestor of the set of taxa with corresponding sequence matches. The total numbers of O. taurus sequences assigned to each branch are indicated in decimals and by the pie chart area (Log scale). Pie graph colors indicate the proportion of contigs (red) and singeltons (blue) assigned to each branch.

Figure 3

Bi-connected components and alternative splicing. A. An example of a bi-connected component structure (BCCs). A BBC composed of three contigs 28477, 25928 and 04341 that share three independent sets of 31, 28, and 8 broken reads, respectively (indicated by dashed line) in relation to the homologous T. castaneum transcript (Tc XM963744). Our analysis suggests that this pattern is reflective of two alternative splice variants present in the Onthophagus transcriptome. B. The three conceptual polypeptide sequences from these contigs align to a contiguous region of the Disabled protein from Tribolium, supporting this hypothesis. Shown are, from top to bottom, two alternative Drosophila Disabled transcripts (dark blue lines; thin light blue lines indicate first methionine (M) and stop codon (*)), the homologous Tribolium sequence (green; no alternative transcripts are known from Tribolium) and the relative positions of contigs 28477, 25928 and 04341. Note that the contig 28477 (light orange), which based on our analysis is a putatively alternatively spliced exon, does not share similarity with the exon that is alternatively spliced in Drosophila.

Figure 4

GO categories. GO annotations associated with 8,504 HomoloGene sequence groups. The distributions of the second and third levels of the GO term annotations of the sampled O. taurus sequences (left charts), were remarkably similar to those on the complete T. castaneum proteome (right charts).

Figure 5

Gene expression patterns are correlated with patterns of genetic variation. Residual SNP frequency (number of SNPs in a contig controlling for contig length and read number) was negatively related to overall expression levels ("A") and the number of tissues (head horn epidermis, thoracic horn epidermis, legs and central brain) in which differential expression was detected in a previous microarray study (N = 48 arrays, reported in Snell-Rood et al. 2010). Statistics are presented in Table 6.

Figure 6

Flow diagram of sequence assembly and annotation. Flow diagram illustration the steps involved in the sequence analysis. Computational steps are indicated by purple boxes, sequences are indicated by blue boxes, and analysis with respect to sequence or annotation databases are indicated by green cylinders.