Genome sequencing and comparative genomics of honey bee microsporidia, Nosema apis reveal novel insights into host-parasite interactions

Species specificity of N. Apis

Light microscopy revealed that fresh N. apis spores were oval shaped and varied in size, with a minimum of 5.15 μm and maximum of 7.415 μm in length, and a minimum of 2.99 μm and maximum of 4.836 μm in width (average 6.41 ± 0.74 x 3.36 ± 0.39 μm, N = 50) (Figure 1a). A single specific band was produced from the genomic DNA isolated from purified spores after PCR amplification with N. apis specific primers, and the sequence analysis confirmed the species specificity of the PCR fragment. No PCR product was generated when N. ceranae specific primers were added to the PCR mixture (Figure 1b).

General features of N. Apisgenome assembly

Contents of the N. Apisgenome

The N. apis genome is very AT rich, with an overall GC content of 18.78%, 7.2% lower than that of N. ceranae. The protein coding regions have significantly higher GC content (25%) compared to the overall GC content for the entire genome, indicating that GC content is not uniformly distributed across the genome and supporting the notion that GC content is higher in coding regions than in non-coding regions. The GC content also varies by codon position, with the 1st and 2nd positions of the coding regions having greater GC content than the third (synonymous) position.

N. apis displays bias in its codon usage and has a pattern similar to that of N. ceranae. Both N. apis and N. ceranae have similarly preferred or optimal codons that have an A or T ending. However, codons ending with A or T are used more frequently in N. apis while codons ending with C or G are used more frequently in N. ceranae (Figure 2). The Maximum-likelihood codon bias (MCB) analysis showed that N. apis has a relatively lower level of correlation between GC content at the third position and the MCB (R² = 0.379, ρ < 0.05), compared to N. ceranae (R² = 0.425, ρ < 0.05) (Figure 3).

The merged gene sets from the different programs resulted in a mixture of 2771 predicted protein coding regions (CDSs) made up of 2,159 (78%) intact and 612 (22%) intron interrupted protein-coding genes. The introns were inserted in random places within the genes and the genes containing introns fell into different functional groups. The consensus sequences “A(G/T)GT” at the donor site and “A(G/A)” at the acceptor site were identified in most of the intron-containing genes (Additional file 1: Figure S1). The density of genes across the genome is 0.33, which is similar to that for N. ceranae. However, no significant correlation between the length of scaffolds and the level of gene density was identified (Figure 4). Among the predicted CDSs, there are 2,684 complete genes and 87 partial genes that lack either a start or stop codon, or both. The length of CDSs averaged 274.88 amino acid residues per CDS.

Of these 2,771 CDSs, 2,105 (76%) had matches to the nr protein database (Threshold E value < 1e-03). Most of the sequences have an e-value between 1.0e-3 and 1.0e-25 (Additional file 2: Figure S2a). Among the 2,771 CDSs, 545 (20%) were assigned as conserved hypothetical proteins, 978 (35%) were assigned with GO terms, 582 (21%) were unannotated proteins and 665 (24%) had no match to the databases (Figure 5). The sequences without positive BLAST hits may be explained by a low similarity of the genes to those functionally similar genes in the database, or new genes that possess new functions in the parasite. The sequence similarity of N. apis within the database by Blast search ranges from 35% to approx 98% and peaks at 52% (Additional file 2: Figure S2b). The sequences of N. apis genes showed the most significant similarity to the sequences of N. ceranae followed by sequences of E. intestinalis and E. cuniculi (Additional file 2: Figure S2c). The majority of functional prediction of the coding sequences are obtained from the UniProt Knowledge Base (KB), a nonredundant protein database that includes PDS, UniProt, Swiss-Prot, TreMBL, and TAIR (Additional file 3: Figure S3a, 3b, and 3c). The N. apis reference names, gene prediction sources, the gene locations, coding strand, and conceptual translations are shown in a Gene-Finding File (GFF) formatted dataset (Additional file 4: Table S1). Table S1 also lists Blast-hit annotation, Enzyme code (EC), InterPro motifs, signal peptide and transmembrane motifs.

A search for motifs in the 100-bp upstream regions of start codons of N. apis genes showed an enrichment of TATA box motifs in the upstream of genes and also yielded three distinct MEME motifs. Of the three over-represented motifs, the “ACCCTT” motif conserved in N. ceranae and E. cuniculi[52] was particularly prominent (E = 1.14e-242), and was present in 37% of the predicted genes. The motif shows a preference for occurring predominantly at the ~15 bp upstream position of the start codon. The “ACCCTT” sequence logo and representative upstream sequences of start codon of N. apis genes are shown in Figure 6.

The predictions by SignalP and TMHMM showed that 399 proteins contained putative signal peptides with over 70% of the probability and significantly high C, S and Y scores and that 348 proteins contained putative transmembrane segments of some form with a score of equal or more than 0.3. Of 348 proteins that had membrane-spanning domains, 114 were predicted to contain an amino-terminal signal peptide. The signal peptides, and transmembrane motifs are also included in the dataset (Additional file 4: Table S1).

The combination of noncoding RNA (ncRNA) prediction programs and BLAST similarity search algorithms identified tRNA genes, rRNA gene clusters (16S, 23S and 5S), and small nuclear RNA (snRNA) genes. Both ARAGORN and tRNAscan-SE yielded similar results for tRNA prediction and showed a diversity of tRNA genes in the N. apis genome. There were 62 tRNA genes from ARAGORN and 58 genes from tRNAscan-SE as well as 16 tRNA-derived pseudogenes in tRNAscan-SE. By combining predictions from two tRNA prediction programs, there were 60 overlapping tRNAs with 42 distinct anticodons and gene lengths of 57 to 68 nucleotides. The 16S, 23S, and 5S rRNA genes were organized as a typical co-transcribed operon and arranged in order of 16S-23S-5S separated by two internal transcribed spacers (ITSs) in a tandemly repeated manner. The 16S, 23S, and 5S genes were identified in twenty-six scaffolds and four scaffolds contain a complete locus of 16S, 23S, and 5S genes. The sequence comparisons revealed polymorphisms among copies of rRNA genes and ITSs with the majority of polymorphisms found in the IGS region. The degree of intragenomic variability observed in 23S rRNA genes was greater than in 16S and 5S rRNA genes. While both 16S and 5S genes displayed 2% intragenomic variability, the average intragenomic variability observed it the 23S genes was about 10% with one 23S rRNA gene showing extreme divergence and sharing only 82% nucleotide sequence homology with the other 23S rRNA genes. While sequence alignment identified short conserved sequence motifs in ITSs, the intragenomic variability observed in ITSs was demonstrated not only by nucleotide substitutions but also by length variance across the genome.

The comparison of the 2,771 putative N. apis proteins with the 2,614 N. ceranae proteins revealed a broadly similar gene set shared between N. apis and N. ceranae. The reciprocal best BLAST search comparing the N. apis gene set and N. ceranae gene set yielded 1,356 one-to-one orthologs that were each other’s mutual best matches; thus, the orthologs corresponded to 49% of N. apis genes and 52% of N. ceranae gene sets. Of 1,356 ortholog pairs for both N. apis and N. ceranae, 137 orthologs had multiple hits in N. apis, suggesting the presence of paralogous proteins. The list of ortholog pairs in the genome of N. apis and N. ceranae can been found in Additional file 5: Table S2.

Functional annotation of protein coding genes

Among sequences that were assigned with gene ontology (GO) functional classification terms, 103 out of 978 sequences were further annotated with Enzyme Commission (EC) codes. Annotated sequences were further mapped to 50 different biological pathways. The KEGG pathway map revealed the presence of a large number of N. apis sequences involved in nucleotide metabolism including purine metabolism and pyrimidine metabolism. Genes encoded all key enzymes involved in glycolysis, the pentose-phosphate pathway, and the trehalose pathway identified in N. ceranae were also found in N. apis, indicating that carbohydrate metabolism is the main process responsible for energy generation in Nosema (Figure 7). N. apis also contained enzymes participating in amino acid and lipid metabolism. It was not unexpected that N. apis lacked genes encoding enzymes involved in electron transfer chains and the tricarboxylic acid (TCA) cycle, which corroborates previous findings in other microsporidia and indicate a strong host dependence for energy production. The KEGG pathway map also revealed that N. apis harbored genes for aminoacyl-tRNA biosynthesis, a process that is essential for accurate translation of genetic information, and genes required for biosynthesis of secondary metabolites and microbial metabolism that are involved in the parasite’s growth, survival, and fecundity.

Of protein sequences with assigned GO terms, most sequences have between 1 and 6 GO terms annotated. Annotation results and distribution, GO-level distribution, number of GO-terms for N. apis sequences with length (x-axis), annotation score distribution; and percentage of N. apis sequences with length (x-axis) annotated are shown in Additional file 6: Figure S4a, 4b, 4c, 4d, 4e and 4d. Of the proteins that had matches to the nr protein database, the most abundant protein class is the binding proteins. Other highly abundant proteins include intracellular transport proteins, ribosomal proteins, and proteins involved in major functional categories such as RNA processing, transcription and translation initiation and regulation, post-traslational modification and energy production. Among GO terms of molecular function, “Binding” was the most predominant group. Of GO terms identified for biological process, “metabolic”, was the most dominant term with 18 subcategories related to the metabolic process. Of categories enriched in the cellular component, protein complex, intracellular membrane-bounded and non-membrane bounded organelles were the largest groups (Figure 8A, 8B, and 8C).

Virulence factors are molecules expressed and secreted by parasites that are keys for causing disease in the host as the parasites enter into and exit out of host cells, and inhibit certain host functions. Genome analysis of N. apis has identified genes encoding a wide array of potential virulence determinants. These potential virulence factors of N. apis include genes encoding 1) polar tube proteins defining structure and invasion mechanisms; 2) spore wall and anchoring proteins involved in host tissue recognition and the initiation of host cell invasion; 3) endochitinase and chitin synthase involved in spore-host cell attachment; 4) proteins involved in pathogenicity regulations and energy parasitism, 5) components of the spliceosome machinery that plays a major role in the generation of proteomic diversity and proteins regulating splicing activity, and 6) putative secretive proteins that are involved in interactions with the host and are essential for modulation of host immunity towards successful completion of the parasite lifestyle in the host (Table 3).

Functional conservation and divergence of protein sequences between N. Apis and N. Ceranae

While putative proteins of N. apis share significant sequence homology and functional similarity for enriched GO terms with N. ceranae, the comparison of the GO terms based on Fisher’s exact test after a false discovery rate correction (p < 0.05) displayed some differences between the two species. Comparison of biological processes (Figure 8A) showed that genes involved in the developmental process, organismal development, signaling, kinase activity, translational activity were on the same level of abundance for both species. However, genes that are functionally linked to cell amino acid derivative metabolic process and reproduction were significantly more common in N. apis than in N. ceranae. In contrast, genes involved in metabolic processes and stress response were greatly expanded in N. ceranae. The protein transport genes were also found to have greater representation in N. ceranae than in N. apis. N. apis and N. ceranae also displayed different protein profiles responsive to stress and stimuli. Genes responsive to stress were over represented in N. ceranae compared to N. apis, while N. apis contained a slightly greater number of genes responding to abiotic stimuli. Both species had about the same number of genes related to responses to endogenous stimuli. While cell homeostasis and response to external stimulus genes were only identified in N. ceranae, symbiosis and encompassing, mutualism through parasitism genes were unique to N. apis. The genes related to response to external stimuli were present only in N. ceranae and absent in N. apis. These gene families may therefore represent species-specific proteins. For molecular function GO terms, a number of differences were observed between the two Nosema species (Figure 8B). Structure molecule activity proteins that contribute to the structural integrity of a complex or assembly within or outside a cell were found significantly reduced in N. ceranae, compared to N. apis. The calcium binding protein that was not in N. apis was present in N. ceranae. An enrichment test of cellular components showed that genes encoding nucleolus, Golgi apparatus, cytosol, nuclear chromosome, mitochondrion and plasma membrane were at approximately similar levels of abundance for both species. However, intracellular membrane-bounded organelles including mitochondrion and endoplasmic reticulum were overrepresented in N. ceranae relative to N. apis. In contrast, nucleoplasm was significantly more common in N. apis than in N. ceranae. The ribosome was also found to be more represented in N. ceranae in relation to N. apis. The vacuole and vesicle proteins were present in N. ceranae but absent in N. apis, suggesting species-specific genes involved in structure and function of cellular components.

Bee sample collection and spore purification

Adult worker bees infected with N. apis were collected from an apiary infected exclusively with N. apis in Beaverlodge, Alberta, Canada in 2007. The species status of the Nosema infection was determined by a combination of morphological and molecular methods. In order to purify spores of N. apis from the infected bees, alimentary tracts were individually removed by grasping the stinger with forceps and pulling it out with gentle traction. The infected alimentary tracts were homogenized in 3 ml sterile water using a conical ground-glass homogenizer with a rough grinding surface and filtered through a Corning (Lowell, MA) Netwell insert (24 mm diameter, 74 mm mesh size) to remove large-sized tissue debris. The filtered suspension was centrifuged (Rotor: Sorvall RC-5B) at 3,000 g for 5 minutes to pelletize the spores. The resuspended pellet was then purified on a discontinuous Percoll (Sigma-Aldrich, St. Louis, MO) gradient consisting of 10 ml each of 30%, 60%, and 90% Percoll solution by overlaying the spore suspension onto the gradient and centrifuging at 8,000 g for 10 minutes. The pelleted spores were further purified with 30 ml of 100% Percoll solution by centrifuging at 5,000 g for 2–3 minutes to separate the spores and lighter cellular debris into the top and bottom layers, respectively. The spore-containing supernatant was collected and diluted with water in the ratio of 1:1 and centrifuged at 8,000 g for 10 minutes to pellet the spores. The purity and size of spores were examined under a light microscope (Nikon, Eclipse TE 300) and photographed with a Nikon Digital Camera (DXM 1200).

Genomic DNA isolation and N. Apis species confirmation

The genomic DNA was isolated using CTAB (hexadecyltrimethylammonium bromide) extraction buffer following spore disruption. Approximately 10⁶ N. apis spores were suspended in 500 ml CTAB buffer (100 mM Tris–HCl, pH 8.0; 20 mM EDTA, pH 8.0; 1.4 M sodium chloride; 2% cetyltrimethylammonium bromide, w/v; 0.2% 2-mercaptoethanol) and broken by adding 500 mg of glass beads (425–600 mm, Sigma-Aldrich, St. Louis, MO) into the tube and disrupting the mixture at maximum speed for 3–5 minutes using a FastPrep Cell Disrupter (Qbiogene, Carlsbad, CA). The mixture was then incubated with proteinase K (200 mg/ml) at 55°C for five hours and treated with chitinase (0.4 U) with an extended incubation period of 16 hours at 55°C. The genomic DNA was twice extracted in an equal volume of phenol/chloroform/isoamyl alcohol (25:24:1), followed by a single extraction in chloroform. The purified DNA was precipitated with isopropanol, washed in 70% ethanol, and dissolved in 100 μl sterile water. The concentration and purity of the DNA were determined by spectrophotometric absorption at 260 nm, and ratios of absorption at 260 nm and 280 nm.

The Nosema species specificity of extracted genomic DNA was confirmed by a PCR assay prior to 454 pyrosequencing. Primers specific to N. apis and to N. ceranae were used individually in the PCR-mediated amplification of genomic DNA isolated from purified spores. The sequences of N. apis primers that amplified 269-bp fragment of the 18S rRNA were as follows: apisF 5′- CCATTGCCGGATAAGAGAGT -3′ and apisR 5′- CCACCAAAAACTCCCAAGAG -3′. The PCR reaction mixture, thermal conditions, and sequences of N. ceranae primers used were as previously described [41]. The amplified PCR products were purified using the Wizard PCR Prep DNA Purification System (Promega, Madison, WI) and the nucleotide sequence identity of the PCR fragments was determined using a Basic Local Alignment Sequence Tool (BLAST) homology search.

Library construction, sequencing, and assembly

Although the high% AT vs.% GC ratio of this genome poses challenges for all high-throughput sequencing platforms, a successful sequencing library was constructed, sequenced and assembled through protocol optimizations. The 454 Titanium sequencing platform offers an ideal combination of throughput and long read length to enable de novo sequence assembly of whole genomes. A paired-end library was constructed according to the manufacturer’s protocol, with modifications and optimizations to enable robotic automation of crucial steps and manipulation of limited, high-% AT starting material. An initial sample of 500 ng DNA was sheared using a Covaris E210 high-performance ultra-sonicator. Following adaptor ligation and double AMPure SPRI bead-based size selection, the resulting library had an average fragment size of 2,072 bp with a standard deviation of 208 bp. This library was sequenced on a full two-region 454 Titanium run resulting in 1,070,575 passed-filter sequence reads with an average quality-trimmed length of 304 bp. During data quality control procedures, 126,461 reads, or 11.8% of the total, were eliminated as PCR-induced duplicates. This represents a typical duplicate read percentage for the 454 sequencing platform. An additional 10% of the reads were eliminated as too short for assembly or containing an inconsistent or partial paired-end linker sequence. Of the remaining reads, 360,731 were bisected at the linker to form successful pairs, while 475,944 remained as unpaired shotgun reads.

These data were assembled using the Celera Assembler (http://sourceforge.net/projects/wgs-assembler/) version 6.1 using the CABOG module, an overlap-layout-consensus sequence assembler that has been recently optimized for use with 454 Titanium paired-end sequence data. The Celera Assembler uses a stringent algorithm that heavily weighs paired-end read constraints and favors shorter, higher-confidence contigs over longer, potentially misassembled contigs. This approach ensures a highly accurate result with minimal mis-assembly.

Repetitive element identification

Genome analysis of N. ceranae showed the presence of a large proportion of the repetitive elements which could interfere with the accuracy of de novo gene prediction. We therefore used RepeatMasker (http://www.repeatmasker.org), Tandem Repeat Finder [63], and LTR Finder [64] in an ab initio manner with the default settings to identify the repetitive elements. Repetitive elements and low complexity DNA sequences in the N. apis genome were then masked for elimination. To avoid the incorrect identification of repetitive elements and low complexity from DNA sequences, manual inspection of the list of repeats was carried out. The repeats that had significant Basic Local Alignment Search Tool (BLAST) hits to the NCBI protein database (E-value cutoff of 1.0E-5) were kept for subsequent gene prediction and annotation.

Protein-coding gene prediction and regulatory motif identification

After the exclusion of a few contaminating honey bee host sequences identified through BLASTN searches against the bee genome, a combination of two programs was used for identification of protein-coding sequences (CDSs) in the masked N. apis genome: GeneMark.hmm.ES (version 3.0) using non-supervised training procedures and Augustus trained with E. cuniculi, and S. cerevisiae. The GeneMark.hmm [65] and AUGUSTUS [66] programs for eukaryotic genomes access the coding potential of DNA sequence based on Hidden Markov Models (HMM). Next, outputs of supervised gene prediction based on sequence similarity to reference organisms (extrinsic method) and ab initio gene prediction relying on intrinsic features of the DNA sequence (intrinsic method) were merged and reconciled with other gene models. The CD-HIT-2D [67] was used for identifying overlapping and unique sets of proteins from different gene prediction programs. For overlapping predictions produced by different programs, the prediction with BLAST evidence was kept. If the overlapping predictions by different programs did not have BLAST evidence or if they all had BLAST evidence, the longest prediction was chosen. To verify the correct prediction of protein-coding genes, all inferred genes larger than 50 aa were further compared across an NCBI non-redundant protein database using BLASTP program with default parameters. Frame-shifts and point mutations were detected and corrected where appropriate. Candidates with a protein identity match ≥ 30%, and BLASTP value ≤ 10–5 were retained as protein-encoding gene candidates.

Nucleotide composition and synonymous codon usage of protein coding genes were investigated using the program INCA2.0 by [68]. The maximum-likelihood codon bias (MCB) was used as a measure for codon usage bias where the contribution to the index of the bias of each amino acid is weighted by an estimation of the likelihood of occurrence of bias on each amino acid [69]. The relative synonymous codon usage (RSCU) of individual codons between N. apis and N. ceranae were compared. Putative transcriptive regulators were searched by identifying over-represented motifs in the 100-bp upstream region of the start codon of each protein-coding gene. The MEME 3.0 software program [70] was used for the motif identification.

Signal peptide and transmembrane region prediction

The N-terminal amino acid sequence of predicted protein from this study was analyzed for the presence of putative signal peptides using the SignalP 3.0 program [71], and a neural network and HMM trained on eukaryotes. The transmembrane regions of predicted protein were identified with TMHMM 2.0, a HMM based transmembrane helix prediction program. The architecture of the hidden Markov model and the neural network algorithm was applied for prediction of the signal peptide and transmembrane region of each translated protein. The default parameters were used for both programs.

Non-coding RNAs

tRNA genes were detected using two tRNA identification programs, tRNAScan-SE version 123 (Lowe and Eddy, 1997) and ARAGORN [72] at the default settings. The tRNAScan-SE search was based on tRNA training data from a mixed general tRNA model. The ARAGORN search was based on the homology of recognized tRNA consensus sequences by specifying a yeast nuclear genome (NC_001133-NC_001148) and the ability of input data to form a base-paired cloverleaf to predict tRNA secondary structure. Two programs are also able to distinguish tRNA pseudogenes from true genes with high sensitivity.

Ribosomal genes coding rRNA in a full genome sequence were identified using both RNAmmer 1.2 server [73] with the default parameters and BlastN searches. Matches were defined as those sequences having greater than 90% similarity and greater than 90% of the query length. Small nuclear RNAs (snRNA) were also identified by BlastN searches. rRNA sequences were annotated based on homology to previously published rRNA sequences deposited in GeneBank.

Ortholog and paralog identification

The possible orthologs of N. apis in N. ceranae, E. cuniculi or S. cerevisiae genomes were predicted by searching for the best reciprocal BLASTP matches with InParanoid [74]. In each of the between-species comparisons, we aligned protein sequences from one species to protein sequences from the other, and corresponding best matches were then aligned back to the original species. On the basis of all of the possible between-species sequence pairs, we reported orthologs and in-paralogs that had multiple hits in the same genome with a bit score cutoff of 40 and an aligned region larger than 50% of the longer query sequence in each pair. We also performed bootstrap analysis to obtain confidence for each ortholog identified. Gene duplications were identified with a focus on proteins in the data set that had a better match within the N. apis genome they reside (excluding the self match) when compared against all other genomes used in the study. In order to avoid assigning paralogs to ortholog pairs, the top hit had to have an E-value 10⁵ times lower (more significant) than the next best hit. In addition, the complete genome sequence of N. apis was aligned individually to the complete genome of N. ceranae to find the corresponding orthologous blocks in the genomes and to illustrate homologs between genomes by using the Mauve 2.2. [75].

Protein domain functional annotation

The predicted proteins were analyzed using the BLASTP algorithm to search for homolog sequences and functionally annotated by gene ontology (GO) terms, InterPro terms (Conserved patterns in sequences, InterProScan, EBI), enzyme classification (EC) codes, and metabolic pathways (Kyoto Encyclopedia of Genes and Genomes, KEGG) using the Blast2GO (B2G) software suite v2.4.2 (http://www.BLAST2go.org, [76]. The functional annotation in B2G was performed using the following steps: 1) performing BlastP to final homologous sequences to the query set; 2) mapping to collect gene ontology terms (GO) associated to the hits after a BLASTP search; 3) selecting GO terms from the GO pool collected in the mapping step and assigning them to the query sequence; 4) grouping the descriptions into ‘GOslim” categories for biological process, molecular function, and cellular component, and 5) performing Enzyme code and InterPro annotations support language by GO terms. Finally, an enrichment analysis was performed to assess the differences in functional classes between N. apis and N. ceranae using one-way Fisher’s Exact Tests. The default settings of B2G were used in every annotation step. The false discovery rate (FDR) at level α = 0.05 was applied for controlling procedures. Gene product names are assigned based on sequence homology to known genes in NCBI databases. The putative proteins without GO annotation were classified as “hypothetical protein with no significant database hit”.

Ethics statement

The apiaries for bee sample collection are the property of the USDA-ARS Beltsville Bee Research Laboratory, MD, USA and the Agriculture & Agri-Food Canada, Beaverlodge Research Farm, Beaverlodge, AB, Canada; therefore, no specific permits were required for the described studies. Studies involved the European honey bee (Apis mellifera), which is neither an endangered nor protected species.