Neospora caninum genome sequence

We sequenced the genome of N. caninum Liverpool strain using Sanger sequencing to ∼8-fold depth. It was assembled into 585 supercontigs with an N50 of 359 kb totaling 61 Mb (Table 1). We constructed a set of N. caninum pseudochromosomes by aligning 242 supercontigs to the fourteen publicly available T. gondii Me49 chromosomes [25] based on predicted protein sequence similarity (Figure 1A). It has been shown previously using our partially assembled sequencing data that N. caninum and T. gondii genomes are largely syntenic [26]. Here we show that for almost all regions where conservation of gene order (synteny) is interrupted, corresponding orthologous regions are found elsewhere in the genome. This suggests that while there may have been a small number of chromosomal rearrangements, there has been very little net gain or loss of genetic content (Figure 1B). The N. caninum Liverpool genome sequence has been added to the European Nucleotide Archive as project CADU00000000.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Chromosomal alignment of N. caninum Nc-Liv and T. gondii Me49 highlighting surface antigen gene families.

(A) Aligned chromosomes of N. caninum (above) and T. gondii (below) showing conservation of synteny and distribution of SRS and SUSA surface antigen gene families. Tandemly repeated genes are shown clustered together. Uncoloured genes had less than 20% unique sequence and expression levels could not be accurately determined. 49 additional NcSRSs were found in UnAssigned Contigs (UACs), while three further TgSRSs were not assigned to chromosomes. (B) Shows putative rearrangements between N. caninum and T. gondii chromosomes. Large (>30 kb) insertions in one genome relative to the other are numbered on the chromosomes of N. caninum (orange) and T. gondii (blue). Red ribbons show regions of protein sequence similarity between these regions. The plot shows that most insertions have a pairwise relationship, e.g. region 13 from T. gondii chromosome VIIa is putatively orthologous to region 24 in N. caninum chromosome IX. Thus these regions are shared and not specific to one organism. The arrow symbol refers to sequence similarity with parts of the comparator genome not currently assigned to chromosomes (UACs). A capital ‘T’ identifies a region with no similarity in the comparator genome. These regions include genes belonging to novel families (TSF and KRUF).

https://doi.org/10.1371/journal.ppat.1002567.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Composition of N. caninum genome in comparison with T. gondii and P. falciparum.

https://doi.org/10.1371/journal.ppat.1002567.t001

Transcriptome sequencing of N. caninum and T. gondii

To determine gene expression differences between species and to improve genome annotation we sequenced the transcriptome of the invasive stage (tachyzoite) of N. caninum Liverpool and T. gondii VEG using mRNA sequencing (mRNAseq) on an Illumina GAIIx machine (Tables S5 & S6). The parasites were grown asynchronously for a period of six days in cell culture. We took samples of RNA at days three, four and six. We found that days three and four showed fairly similar expression profiles within species and so we have pooled this data for most analyses. We found however that day six N. caninum parasites were showing expression of bradyzoite (quiescent stage) marker genes (Text S1). These parasites had not fully converted into bradyzoites, but may be preparing to do so. We did not observe expression of these markers at day six in T. gondii, so we did not seek to compare transcriptomes of the two species at this timepoint. Transcriptome sequencing data has been submitted to ArrayExpress with accessions E-MTAB-549 for N. caninum sequences and E-MTAB-550 for T. gondii sequences.

Gene content is largely conserved between species

Combining de novo gene predictors and mRNAseq evidence we identified 7121 protein-coding genes in N. caninum and produced a revised T. gondii ME49 gene count of 7286, a reduction of 9% from previous predictions (Table 1). This was achieved predominantly by merging adjacent genes based on mRNAseq evidence. In N. caninum we detected the expression of 74% of genes during the tachyzoite stage. In T. gondii 80% were expressed, significantly more than the 49% recently reported using microarrays suggesting greatly improved sensitivity [27].

Using a combination of automated orthologue identification and manual curation we identified a small number of unique (i.e. organism-specific) genes in both genomes that might underlie their phenotypic differences (Figure 2A). Excluding surface antigen families (discussed later), we found 231 genes unique to T. gondii and 113 to N. caninum, i.e. with no orthologue or paralogue based on our orthologue analysis. Of these, 72 from T. gondii and 43 from N. caninum had Pfam domains or proteomics-based evidence (Table S1). These genes represent good candidates for understanding organism-specific differences and are enriched for those involved in host-parasite interactions. The remainder had no detectable homologues or proteomics-based evidence, although most had good transcriptome evidence.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Protein-coding gene content and metabolic activity are largely conserved between the two species.

(A) Most protein-coding genes in N. caninum have a one-to-one orthologous relationship (yellow) with a gene of T. gondii. A larger proportion of the T. gondii genome consists of genes with no N. caninum homologue than vice versa (organism-specific genes in red). The increase in shared multi-gene families (blue) in N. caninum reflects the expansion of SRS genes in this organism. The increase in organism-specific multigene families (red) in T. gondii reflects, for instance, the TSF gene family identified by us in this work. (B) Of the one-to-one orthologues shared by T. gondii and N. caninum we identified those which have orthologues in three or more non-apicomplexan eukaryotes (yellow), are not present in three or more apicomplexans but in all apicomplexan groups sequenced to date (grey), are in at least one other apicomplexan group (blue) or are specific to T. gondii and N. caninum (red). (C) Pooled day three and four RNAseq experiments were used to determine orthologous genes differentially expressed between T. gondii and N. caninum. Differentially expressed genes were examined for enrichment with enzymes from different KEGG pathways as described in methods. No pathways were identified with a p-value less than 0.01, although putative differences were identified at a p-value cutoff of 0.05. The width of bars in the chart relates to the number of genes in the pathway which were differentially expressed, e.g. Porphyrin metabolism had three enzymes differentially expressed, while Nitrogen metabolism had two and Tyrosine metabolism one. Only pathways with a p-value (adjusted for multiple hypothesis testing) below 0.1 are shown. (D) Gene Ontology terms over represented amongst genes upregulated in N. caninum and T. gondii. All terms shown are significantly upregulated (P<0.05). The terms ‘membrane’, ‘regulation of transcription, DNA-dependent’ and ‘ATP binding’ are found more often than expected in genes upregulated in T. gondii. SRS surface antigens, rhoptry kinases and AP2 transcription factors respectively are associated with these terms, suggesting that SRSs, ROPs and AP2s are amongst the most highly upregulated groups of genes in T. gondii relative to N. caninum. The term ‘protein amino acid phosphorylation’ is overrepresented amongst genes upregulated in N. caninum relative to T. gondii. Many of the genes in this group are rhoptry kinases suggesting that while some are upregulated in T. gondii, others are upregulated in N. caninum. These findings are explored in more detail in Figure 3.

https://doi.org/10.1371/journal.ppat.1002567.g002

Only one organism-specific multigene family, with no homologues in the other species was identified: a family that we have named Toxoplasma-specific family (TSF; Figure S1). This family is located largely in chromosomal regions with no similarity to N. caninum (regions 3, 5 and 17 in Figure 1B) and varies in size between T. gondii strains. We found that all ten members of TSF from T. gondii Me49 were expressed during the tachyzoite stage. No significant domains, motifs or signal peptides were identified; however a putative transmembrane helix was predicted between amino acids 195 and 217 on TgTSF1.

Another previously unidentified family was present in N. caninum, but was expanded in T. gondii (Figure S2). This family comprises three genes from N. caninum and seven from T. gondii Me49. Sequences were scanned using InterProScan [28] but no significant domains or motifs were identified. Lysine-arginine rich motifs are however present in the sequences suggesting possible nuclear localization signals. We have therefore named this family Lysine-Arginine rich Unidentified Function (KRUF). KRUF genes appear to be highly expanded in the GT1 strain of T. gondii, with up to twenty members [25]. Two of the three N. caninum members are expressed in the tachyzoite and early bradyzoite stages (NCLIV_002020 and NCLIV_002030). Most of the T. gondii members are expressed, some at very high levels (esp. TGME49_051170).

While 32% of the genes shared by T. gondii and N. caninum have orthologues in a range of eukaryotes, we found that ∼39% of the shared genes do not have orthologues in other apicomplexans sequenced to date (Figure 2B). Furthermore, while ∼29% of the shared genes not found outside apicomplexans have orthologues in at least one apicomplexan, only 0.3% are shared between all apicomplexans with completed genome sequences. Due to the assumptions behind this analysis we have likely underestimated the similarity between Apicomplexa and more detailed manual analysis will no doubt reveal more divergent orthologues. However our results suggest that the genome content of apicomplexans is very diverse and that many novel and divergent genes are found within the Coccidia.

Metabolic capacity is conserved between species

To determine whether their divergent lifestyles are associated with differences in metabolism we compared the predicted repertoires of metabolic enzymes and pathways of N. caninum to those of T. gondii [25]. The pathways identified in N. caninum appeared identical to those in T. gondii and we found no single metabolic gene specific to either species suggesting that changes in metabolism do not play a large role in host restriction and zoonotic compatibility in these species. Although a small number of metabolic genes were differentially expressed between species, we found little evidence that these were clustered in any particular pathway, although there is some evidence that nitrogen metabolism may be upregulated in N. caninum and porphyrin metabolism may be upregulated in T. gondii (Figure 2C).

T. gondii and N. caninum diverged from their common ancestor around 28 million years ago

Previous estimates from rRNA analysis have suggested that N. caninum and T. gondii diverged between 12 and 80 million years ago (mya) [6], [29]. To gain a more accurate estimate we examined a large number of orthologue alignments, determining synonymous substitution rates between N. caninum and T. gondii and between malarial parasites of human and non-human homonidae: Plasmodium falciparum and P. reichenowi, respectively. We assumed constant evolutionary rates between the Plasmodium spp. and Coccidia, excluding genes which were found to have evolved in a non-clock-like manner. We used a previously determined estimate of 2.49 mya for the split between P. falciparum and P. reichenowi [30]. This allowed us to date the speciation of N. caninum and T. gondii to 28.0 mya or between 21.7 and 42.7 mya using the confidence intervals of the P. falciparum and P. reichenowi divergence time. This suggests that speciation of N. caninum and T. gondii occurred after the speciation of their definitive hosts (estimated at 54–67 mya) [31]. The ability to reject non-clock-like genes is dependent on gene length and so we also calculated the divergence time using only the longest 25% of the orthologous groups. This led to a divergence time of 26.9 mya, very close to that calculated using all groups, suggesting that a tendency to exclude longer genes using the clock test has not biased our results.

Surface antigen gene families are greatly expanded in N. caninum

Examination of gene gain and loss and differential expression implicated two host-interaction gene families: SAG1-Related Sequence (SRS) and ROPK, as among the most divergent features of the two species (Figure 2D, Figure 3). SAG1 was the first SRS protein identified and is the major surface antigen of Toxoplasma. SRS proteins localize to the cell surface of both T. gondii and N. caninum. They are thought to play a key role in attachment to host cells, modulation of host immunity and regulation of parasite virulence [32]. Wasmuth et al. (submitted) found 109 functional genes and 35 pseudogenes in T. gondii Me49 with similar numbers across several different strains. They are present sometimes in single copies, often in tandem arrays. They are dispersed across all chromosomes rather than showing a preference for subtelomeric regions as is found for some large gene families in Plasmodium, Babesia and Theileria (Figure 1A). It has been suggested that the large number of SRS genes is present in T. gondii to accommodate the wide spectrum of potential host-cell molecular interactions presented by its exceptionally large host range [33], [34]. However, our data refute this; we found the SRS gene family to be substantially expanded in N. caninum compared to T. gondii with a total of 227 N. caninum SRS genes (NcSRSs) and 52 NcSRS pseudogenes (Figure 1A). Expression data suggest however that T. gondii expresses a greater number of its SRS repertoire (55 vs. 25 in N. caninum) during the tachyzoite stage (Figure 1A). In N. caninum we found in most cases that only a single SRS gene was expressed at a multigene locus, whereas in T. gondii we often found several. Extending our gene expression studies beyond the rapidly growing and invasive tachyzoite stage, we noticed that N. caninum cultures maintained until day six showed expression of known bradyzoite-specific genes (e.g. BAG1, SRS13, SAG4), suggesting they were beginning to convert into the slow-growing stage (Text S1). We observed a greater number of NcSRS genes (36 vs. 25) expressed at day six than at earlier points in the culture. Despite this it remains unclear whether most members of this expanded family in N. caninum are expressed and further expression data are required from all life-stages before the role of these genes can be better understood.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Repertoires and differential expression of known and predicted host-interaction genes and AP2 transcription factors.

It was only possible to reliably identify orthologues for 22 SRS genes due to the way they have expanded, often being in large tandem arrays subject to gene conversion. While the AP2 transcription factors are not directly involved in host-parasite interaction they may be important in regulating expression of invasion genes. Each report card details the comparative repertoires of a particular group of genes in these species, the names of the genes specific to each organism and those which are differentially expressed between organisms. Further details of these relationships, including reference numbers, are included in Table S4. Arrows show the fold change in expression (RPKM; Reads Per Kilobase per Million mapped reads) between the two species on a log₂ scale. The ticks are 2, 6 and 8 on this scale. Green arrows highlight increased expression in N. caninum tachyzoites. Red arrows highlight genes with increased expression in T. gondii tachyzoites vs. N. caninum tachyzoites. A fold change is infinite where the gene is not expressed at all in one organism.

https://doi.org/10.1371/journal.ppat.1002567.g003

SRS genes consist of one or more copies of the SAG domain family, which has been classified into eight subfamilies (Figure S3; Wasmuth et al., submitted). The doubling of SRS gene numbers in N. caninum compared to T. gondii is largely accounted for by expansion of a particular subfamily with a 7–8 domain architecture. No novel SAG domain subfamily has evolved in either lineage, however several domain combinations are found in low copy numbers in only one or the other species (Table S2). Since a particular SRS locus tends to contain genes with the same domain architecture in both species, expansion has likely occurred by tandem duplication. We found evidence that gene conversion may have occurred at, at least, one locus (SRS19; Figure S4A), whereas one of the most highly expressed loci in both organisms (SRS29, containing the SAG1 gene) showed no evidence of gene conversion (Figure S4B), perhaps due to functional constraints.

SUSA genes (SAG-Unrelated Surface Antigen genes) are a superfamily of surface antigens unrelated to SRSs but which are also postulated to interact with the host immune system [35]. In common with the SRS superfamily we found that N. caninum had an expanded number of SUSA genes but that a greater number were expressed in T. gondii (Figure 1A). In fact none of the NcSUSA genes were expressed in the tachyzoite stage. Two NcSUSA genes (NCLIV_067570 and NCLIV_067920) were however expressed at day six of culture.

Variation in rhoptry kinase genes highlights evolution of host-parasite interaction

The apical complex is the defining characteristic of the Apicomplexa. It includes the rhoptry, microneme and dense granule secretory organelles, which are essential for cell invasion. Figure 3 shows how the repertoires and expression of gene products known or predicted to be localized to these organelles differs between T. gondii and N. caninum. Several T. gondii rhoptry genes (ROP18, ROP16 and ROP5) have been implicated in virulence based on a genetic cross between the type II and III [3] and type I and III [4] lineages of T. gondii. N. caninum differs from T. gondii at each of these loci, but shares some similarities with low virulence strains.

Pseudogenisation of ROP18 in N. caninum prevents ROP18-mediated inactivation of immunity-related GTPases

Expression of TgROP18 is associated with virulence in mice [4] and in some hosts high ROP18 expression may reduce parasite fitness by causing rapid host death [36]. It is involved in preventing the host interferon-gamma (IFN-γ) response, during which the host loads immunity-related GTPases (IRGs) onto the parasitophorous vacuole (PV) leading to its disruption and parasite cell death in avirulent strains [37]. Virulent T. gondii strains express high levels of ROP18, which phosphorylates and inactivates IRGs to safeguard the PV [24], [37]. We found that in Neospora Nc-Liv ROP18 is a pseudogene due to several interrupting stop codons in the sequence syntenic with the Toxoplasma gene. We confirmed the presence of these stop codons in a further four strains of the parasite isolated from different geographic locations and hosts (Table S3). To determine whether N. caninum is able to phosphorylate IRGs without a functional copy of ROP18 we examined the loading of Irga6 (a member of the host IRG GTPase family) onto the PV by immunofluorescence studies. We observed that, in both N. caninum and T. gondii infections, host cells responded by loading Irga6 onto the PV but only T. gondii was able to phosphorylate Irga6 and thereby presumably inactivate the IRG protein (Figure 4A). This suggests that N. caninum is unable to prevent its host from using IRGs to attack the PV.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. ROP18-mediated IRG inactivation is absent in N. caninum and ROP5 has a lower copy number in N. caninum than T. gondii.

(A) Phosphorylation of Irga6 T102 is observed in Mouse Embryonic Fibroblasts (MEFs) infected with T. gondii RH (a virulent strain) but not in N. caninum NcLiv-infected MEFs despite loading of Irga6 onto the PV. (B) Genomic Illumina sequencing reads are shown mapped to a reassembled ROP5 locus in N. caninum. This sequencing was performed using a PCR-free protocol to remove biases commonly introduced by PCR, resulting in more even coverage. The SSAHA mapping algorithm maps reads to multiple locations if they map equally well (mutireads), but gives these reads a mapping quality score of 0. The upper plot shows all mapping reads (including multireads) and slightly elevated read depth is visible over most of each ROP5 gene. Removing reads mapping with a mapping quality of less than ten shows that reads cannot be mapped reliably to most of each gene because they are almost identical in sequence. Therefore, allowing all reads to map, we would expect an average doubling of read depth over these genes if there were really only two copies and an increase in read depth of four times above background if there were four copies. We conclude that there are only two copies of the ROP5 gene, as in the original assembly.

https://doi.org/10.1371/journal.ppat.1002567.g004

ROP16 and other virulence-associated rhoptry kinases

In T. gondii ROP16 directly interferes with host signaling pathways (e.g. Stat3, Stat6) to modulate the proinflammatory host cytokine IL-12 [21], [38]. A single polymorphic residue on TgROP16 determines the strain-specific activation and phosphorylation of Stat3 [39]. We found that ROP16 was highly expressed in T. gondii VEG (type III) tachyzoites and it has been shown elsewhere to be highly expressed in all T. gondii strain types [27]. While its orthologue in N. caninum, NcROP16, possesses the key active-site leucine residue for Stat3 activation, the gene was not expressed in our experiments. Although it is possible that NcROP16 could be expressed in other cell types, our experiments predict that N. caninum infection does not activate Stat3 due to its lack of expression. Several additional T. gondii rhoptry genes are missing from N. caninum (Figure 3), most notably the entirety of the locus which encodes ROP2A, ROP2B and ROP8. The TgROP5 multigene locus accounts for 50% of inherited variation in Toxoplasma virulence [40]. The relationship between ROP5 genotype and virulence in T. gondii is however not clear. The most virulent, type I T. gondii strain (e.g. RH) has six copies, while the less virulent type II T. gondii strain (e.g. Me49) has around ten copies and the least virulent type III strain (e.g. VEG) has four. N. caninum Liverpool encodes only two copies of the ROP5 gene both of which are highly expressed in the tachyzoite stage (Figure 4B).

Divergence among other apical organelllar genes

The secreted proteins of the microneme organelles play a crucial role in host cell attachment and invasion by mediating gliding motility [16]. We identified thirteen previously undescribed genes putatively encoding micronemal proteins by virtue of conserved domain architectures. Of these newly identified genes, MIC26 (a MIC2 paralogue) and MIC19 (a PAN domain-containing gene) are unique to N. caninum. Some differences also exist between the species in dense granule genes which are involved in the modification and function of the parasitophorous vacuole (PV) [41]. Dense granule genes GRA11 and GRA12 were absent from the N. caninum genome sequence.

Serine proteases are important to the maturation of both rhoptry and microneme proteins and their inhibition blocks parasite replication and rhoptry formation [42]. TgSUB2, a subtilisin-like serine protease has been identified as a likely processor of several rhoptry proteins [43] and whilst T. gondii is vulnerable to a variety of protease inhibitors, including serine protease inhibitors, N. caninum invasion is inhibited only by aspartyl protease inhibitors [44]. We found that while all 12 identifiable T. gondii subtilases had orthologues in N. caninum, there was a significant decrease in expression of these proteases (hypergeometric test; p = 0.003) compared with T. gondii. This suggests that subtilisin-like serine protease activity may not be used to the same extent in N. caninum as in T. gondii and may explain why N. caninum is less susceptible to its inhibition.

Evolution of transcriptional regulation

The ApiAP2 family represents the major group of apicomplexan transcription factors. They have been implicated, for example, in control of the intraerythrocytic development cycle (IDC) [45] and sporozoite development [46] of malaria parasites. Twenty-nine such genes have been identified in Plasmodium and 68 in Toxoplasma [47]. We found N. caninum orthologues for all 68 TgAP2 genes but detected significant differences in the expression of eleven of them (Figure 3), which in turn may be responsible for expression differences we have observed in other genes. It has been suggested for instance that rhoptry genes are regulated by AP2 transcription factors in Plasmodium [48]. We found that 54 of 68 NcAP2s and 61 of 68 TgAP2s were expressed during the tachyzoite stage, more than a previous study [49]. This is surprising considering that one would expect the principal family of transcription factors in organisms with a complex life cycle to be highly specific to different life stages.

As expected, the repertoire of ncRNA genes of known function (e.g. t-RNAs, snoRNAs, snRNAs etc.) is almost identical between Toxoplasma and Neospora. However, we were able to identify an expansion of a previously unidentified candidate structured non-coding RNA family in N. caninum. This suggests that ncRNA repertoires are divergent in these species, although the functions of these RNAs remain to be identified (Text S2, Figures S5 & S6).

Parasite cultivation

Neospora caninum Liverpool strain was originally isolated from the cerebrum of a congenitally infected dog [54]. N. caninum Liverpool and Toxoplasma gondii VEG tachyzoites were maintained as described previously [55].

Genome sequencing and assembly

Paired end reads of N. caninum DNA were generated from random subclone libraries and additional reads were directed to close gaps and improve the data coverage of low quality regions. All sequencing was performed using BigDye terminator chemistry and used AB 3730xl analyzers (Life Technologies). In total, 920k reads were obtained, quality-clipped and screened for contamination. 92% of reads were used in the final assembly. Based on an estimated genome size of 62 Mb for N. caninum the sequencing coverage was ∼8x. Sequence reads were assembled using PHRAP (P. Green, unpublished) into 960 supercontigs with an N50 of 354 kb. To reanalyse the ROP5 locus these reads were reassembled using Arachne [56].

N. caninum pseudochromosomes were generated by aligning supercontigs to T. gondii ME49 chromosomes using PROmer [57]. 242 contigs (90.4% of the sequence) aligned successfully to the 14 chromosomes of T. gondii. Of the remaining 718 contigs, 375 were removed due to contamination, poor quality or if they were <1 kb in length. The remaining 343 contigs were grouped as UnAssigned Contigs (UACs) and used in further analysis alongside the pseudochromosomes. Telomeres were identified by examining chromosome ends for the typical TTTAGGG septameric repeat.

A genome resequencing library for N. caninum were prepared as in [58]. Sequencing was performed on an Illumina GAIIX as for transcriptome libraries. Illumina paired-end reads were mapped using SSAHA2 [59].

The ROP18 region of five N. caninum isolates (Table S3) was amplified in two overlapping sections using the following primer pairs: F1+R3 and F3+R1 (exp. product 1268 bp and 890 bp respectively) (supplied by Eurofins). ROP18_F1 – GAGTGCCACGGTCCTCTAAG, ROP18_R3 – ATTTGTCCGACGCAAAATTC, ROP18_F3 – GGCTTCTGCTCCAGTATTCG, ROP18_R1 – GCCTTATAAACCACCCGTCA. PCR Reagents were supplied by Qiagen.

Toxoplasma genome sequences and gene models were downloaded from ToxoDb v5.2 (http://www.toxodb.org); they were generated at the J. Craig Venter Institute and have been kindly provided by the Toxoplasma research community.

Gene finding and annotation

N. caninum gene models were created using several algorithms [60]–[66] trained on T. gondii Me49 (ToxoDB v4.2) and using ESTs from N. caninum Liverpool and NC-1 strains collected by the Gene Index Project [67]. The models were examined using the Artemis Comparison Tool [68] and where possible corrected based on evidence from synteny and sequence conservation with T. gondii and transcriptome sequencing evidence. We found a large number of erroneously unfused gene models. Using our T. gondii transcriptome data a total of 354 pairs of adjacent, same-strand genes were linked by reliably mapped bridging read-pairs. A further 449 genes in T. gondii were found to be parts of an adjacent gene but did not have spanning read pairs, usually being likely UTR segments. This resulted in a large drop in the predicted T. gondii gene count and we incorporated these corrections into our subsequent analysis.

We used orthoMCL [69] to identify an preliminary set of orthologous groups between T. gondii and N. caninum. These results was modified using 679 manually identified orthologue pairs. We identified 6348 one-to-one orthologous gene pairs, which we then used to determine whether genes in these organisms tend to be shared with other apicomplexan species. We performed an orthoMCL with representative genes for the one-to-one (core) T. gondii/N. caninum set as well as predicted protein sequences for the following species: the plant Arabidopsis thaliana, the piroplasmic apicomplexan Babesia bovis, the apicomplexan Cryptosporidium parvum, the slime-mould Dictyostelium discoidium, human, the haemosporian apicomplexan Plasmodium falciparum, the yeast Saccharomyces cerevisiae, the piroplasmic apicomplexan Theileria annulata, the kinetoplastid Trypanosoma brucei and the diatom Thalassiosira pseudonana. Where a core gene was found to have an orthologue in three or more non-apicomplexan eukaryotes, we defined it as eukaryotic. Where it was not eukaryotic, but conserved amongst all apicomplexa, we defined it as conserved apicomplexan, If a gene was not conserved apicomplexan, but found in one or more apicomplexan species, other than T. gondii and N. caninum it was defined as apicomplexan. Remaining genes were considered specific to T. gondii/N. caninum.

SRS genes and SAG domains were identified as in Wasmuth et al. (submitted). We identified pseudogenes as clusters of significant BLAST hits which did not overlap valid gene models. Putative pseudogenes were manually checked to determine whether rational gene models could be made and whether Ilumina resequencing data supported any stop codons. N. caninum SAG domains were clustered as in Wasmuth et al., to identify any novel domain subfamilies.

N. caninum ROP, MIC, GRA and AP2 genes were initially determined by manually identifying orthologues of known T. gondii genes with reference to various studies [70]–[77]. Where homologous families of proteins fell into these groups, e.g. ROPK family for ROP, TRAP and MAR for MIC, novel members were sought using BLAST and HMMer.

Transcriptome sequencing and mapping

Poly A+ mRNA was purified from total RNA using oligo-dT dyna bead selection followed by metal ion hydrolysis fragmentation with the Ambion RNA fragmentation kit. 1st strand cDNA was synthesized using randomly primed oligos followed by 2nd strand synthesis to produce dscDNA. Fragments were selected for 200–250 bp inserts amplified by PCR to enrich for properly ligated template strands. Libraries were sequenced using the Illumina Genome Analyzer IIX in paired end mode for 2×76 cycles using proprietary reagents according to the manufacturer's recommended protocol.

RNA-seq reads were aligned against the reference genomes using SSAHA2 [59]. Reads were included only where one end of the pair aligned uniquely to the genome and the distance between the pairs was within the expected insert size range, plus the expected intron length (80–4000 bp).

Expression analysis

We used Reads Per Kilobase of exon model per Million mapped reads (RPKM) normalised by the unique length of the gene as a measurement of expression level. We excluded positions which were non-unique from the length calculation using a kmer window of 75 bp, 37 bp either side of that position. Non-uniquely mapped reads were excluded by removing reads with a score <10. In order to determine whether or not a gene was expressed we calculated an RPKM threshold (Figure S7).

We used DESeq to determine differentially expressed genes [78]. In each pairwise comparison of two conditions A and B (e.g. N. caninum day 4 tachyzoites with T. gondii day 4 tachyzoites) biological replicates were used for both A and B to gain more accurate estimates of experimental variation. Genes with an adjusted p-value of <1e-5 were considered differentially expressed. When considering differential expression between species rather than between different time points in the same species we considered only genes identified as pairwise orthologues. Orthologues of N. caninum and T. gondii are often different lengths and therefore we normalised the read counts for T. gondii genes based on the gene length of the N. caninum orthologue.

Divergence time between species

We determined orthologous relationships between N. caninum, T. gondii, P. falciparum and P. reichenowi using orthoMCL [69]. Orthologous groups containing a single gene from each species were aligned using muscle [79] and those with less than 50% conserved positions across all four species (including gaps) were excluded, leaving 184 alignments. We further excluded those orthologous groups which we determined have not evolved in a clock-like manner. To do this we used a likelihood ratio test for a constant rate of evolution [80]. Likelihood computations on a fixed species tree under a model where branch lengths are free to vary, and under a model in which branch lengths were constrained to be clock-like were performed in PAUP 4b10 [81], under a general time-reversible model [82] incorporating both a proportion of invariant sites and a gamma distribution of rates across sites [83], [84]. The GTR+I+G model of evolution was applied for each locus independently. We performed this test for single-copy orthologous sequences, determined as above, from T. gondii, N. caninum, P. falciparum, P. reichenowi and P. berghei. Twelve of the 184 alignments were excluded because they were determined to have evolved in a nonclock-like manner. As clock-like evolution is the null hypothesis in this test, failure to reject a molecular clock can be either due to the true process of evolution being clock-like, or very close to clock-like for a particular locus, or because of a lack of power in the data to reject the clock. To test that this was not introducing a significant bias in our results, we looked at the estimate of divergence times from two different sets of loci: all loci that fail to reject the clock model, and only those genes in the top quartile of alignment length (which should have the most statistical power to reject the clock).

Codeml [85] was used to calculate the maximum likelihood value of dS in pairwise runmode with the JTT model allowing 2 or more dN/dS ratios for branches. Using all 172 alignments the median 4-fold coding site synonymous substitution rate (dS) across pairs of N. caninum/T. gondii orthologues was 0.856 substitutions per site. Between P. falciparum and P. reichenowi this was 0.076, similar to that calculated by Neafsey et al. [86] (0.068; 95% CI [0.060–0.077]). We assume that these Plasmodia diverged 2.49 mya (95% CI [1.93–3.79]) [30]. We thus dated the split for N. caninum and T. gondii to 28 mya, or taking into account the confidence intervals for the Plasmodium divergence estimate, between 21.7 and 42.7 mya, after the divergence of the definitive hosts around 52.9 mya [87].

If we calculate the median dS values using only those longest 25% of the 172 orthologous groups, we get a dS of 1.230 for N. caninum/T. gondii and 0.114 for P. falciparum/P. reichenowi. This translates to a divergence time of 26.9 mya, This value is very close to that calculated using all 172 alignments. This suggests that a any tendency to exclude longer genes using the clock test has not biased our results.

Comparison of metabolic capacity

Enzyme Commission (EC) number mappings were extracted from the KEGG database [88] from 23 different species covering prokaryotes, archaea and eukaryotes and were mapped on to the corresponding genes in OrthoMCL database [69]. All N. caninum proteins that shared orthology with these enzymes were transitively assigned one or more EC number.

KEGG pathway mapping/coloring tools were used to map EC numbers to pathways. The final set of N. caninum metabolic pathways was compared to that of T. gondii (EC numbers assigned and used in similar fashion to Neospora). Pathways containing significantly high numbers of gene expression differences were determined as discussed in Statistical analysis.

IRG phosphorylation assay

Cell culture was performed as described in [37]. The following immunoreagents were used (dilutions in parentheses). From J.C. Howard (University of Cologne): mouse anti-Irga6 monoclonal antibody (mAb) 10E7 (1∶500) [89], anti Irga6 phosphopeptide Ab T102-555 (1∶5000) [23] Alexa 350/488/546/555/647-labelled donkey anti-mouse, rabbit and goat sera (Molecular Probes), donkey anti-rabbit- (GE Healthcare), donkey antigoat- (Santa Cruz Biotechnology) and goat anti-mouse- HRP (horseradish peroxidase) (Pierce) antisera (all 1∶1000). From P. Bradley (UCLA): mouse anti-N. caninum ROP2 family member monoclonal antibody 20B5D5 (1∶2000). 4′,6-Diamidine-2′-phenylindole dihydrochloride (DAPI, Invitrogen) was used for nuclear counterstaining at a final concentration of 0.5 mg ml⁻¹. Saponin permeabilization and immunostaining were performed as described in [90], [91], except for slides stained with T102-555 which were permeablized in ice cold methanol for 20 min and blocked in 1% BSA in PBS for 30 min.

1D SDS-PAGE

N. caninum proteins were purified from a tachyzoite pellet and resolved into 127 contiguous bands using acrylamide gel electrophoresis. Bands were excised and digested with trypsin. LC MS/MS was carried out using an LTQ ion trap mass spectrometer (Thermo Fisher Scientific Inc, Waltham, MA, USA) with an electrospray ionization source, Tryptic peptides were eluted using a linear gradient of 0–50% (v/v) acetonitrile/0.1% (v/v) formic acid over 140 minutes followed by 100% (v/v) ACN/0.1% formic acid for 20 minutes and a further 20 minutes of 0% (v/v) acetonitrile/0.1% (v/v) formic acid. Protein identifications were made as in [92], those above 1% false discovery rate were discarded. 1053 proteins were found to have at least one significantly matching peptide.

Statistical analysis

To determine whether certain gene functions were overrepresented in differentially expressed genes we assigned GO terms using InterPro2GO [93]. The hypergeometric test was used to determine overrepresented GO terms in pooled day three and four expression data with a DESeq q-value cutoff of 1e⁻⁵. The Benjamini-Hochberg method was used to correct for multiple hypothesis testing. Values of P<0.05 were considered significant. The hypergeometric test was also used in the same way to determine whether KEGG metabolic pathways were enriched in differentially expressed genes.