Analysis of Theileria orientalis draft genome sequences reveals potential species-level divergence of the Ikeda, Chitose and Buffeli genotypes

Reference genome read mapping and de novo assembly

SNV validation and analysis

Variation between the three isolates examined in this study was relatively high, reflecting the diversity of this parasite. Moreover, many contigs from the Fish Creek and Goon Nure genome assemblies did not align to the Shintoku reference, or other isolates from this study (Additional file 1). BLAST searches of these contigs using the NCBI non-redundant database produce no significant hits. Novel sequence and differences in gene content are discussed in more detail below.

When comparing variation between genotypes, the Robertson isolate showed expectedly high similarity to the Shintoku reference sequence, while Fish Creek and Goon Nure strains showed similar levels of divergence (Table 2). Total numbers of variants were similar in the Fish Creek and Goon Nure strains, however, these variants were largely found in differing positions reflecting the divergence between all three isolates. Coding sequences showed a higher variant density than non-coding sequences (Table 2), but this is likely a consequence of greater coding sequence coverage (Table 1). In validation experiments with simulated reads and Sanger sequencing, we noticed that in sequence regions with a high number of short indel mutations read mapping coverage was reduced. These high indel regions seem to be more common in the non-coding genome and result in reduced coverage. Excepting chromosomal ends where sequencing coverage was low, variants were evenly distributed in Fish Creek and Goon Nure sequences (Additional file 4). Some high variation regions were identified in the Robertson sequence and these areas largely correspond to hypothetical genes (Additional files 4 & 5).

Ortholog clustering and novel gene prescence

To explore differences in gene content between the four genomes, we examined predicted proteins using OrthoFinder [21]. Using an evidence-based annotation methodology for eukaryotic genomes, we identified 3675, 3624 and 4789 genes in the Robertson, Fish Creek and Goon Nure strains, compared with 4002 genes identified in the Shintoku genome sequence. Ortholog clustering identified that the higher gene content detected in the Goon Nure sequence is due to a higher number of genes unassigned to orthologous groups shared by the four genomes (515 vs 59–114) and a higher amount of gene duplication, which is demonstrated by increased number of orthogroups with multiple genes from Goon Nure compared to Robertson/Chitose (Additional file 6). This increase is present when OrthoFinder is run with an increased similarity threshold (e-value = 1e-100), indicating higher gene copy number could be produced by haplotype diversity as outlined above or by gene duplication within the Goon Nure genome. To examine this further we mapped sequencing reads to predicted coding sequences for each isolate with bowtie2, using the -k 1 option so that each read was only aligned once. We then examined if there was a stoichiometric-like relationship between average gene coverage and the per-isolate number of genes in an orthogroup. If duplication of genes was caused by the generation of highly similar quasi-species contigs during the assembly process, then orthogroups with higher numbers of genes per-isolate should show lower coverage in stoichiometric-like ratios. Robertson and Fish Creek showed relatively stable coverage regardless of per-isolate orthogroup gene number (Additional file 7). For Goon Nure, average gene coverage was higher in orthogroups with lower numbers of genes but differences were less than expected for a stoichiometric-like ratio, indicating that increased genome duplication or other reasons may contribute to the fragmentation observed.

To identify unique genes present within Robertson, Fish Creek and Goon Nure strains we first examined the similarity of predicted proteins to those identified from the Shintoku whole genome sequence using blastp, as this method would identify proteins that were both present in each genome and also predicted by each annotation pipeline. A total of 29, 32 and 69 predicted coding sequences from Robertson, Fish Creek and Goon Nure strains respectively did not show similarity (e-value cutoff 1e-5) to any predicted protein sequence from the T. orientalis Shintoku genome.

Predicted protein coding sequences that did not match Shintoku predicted proteins were examined further with tblastn to determine if they were present within the Shintoku genome sequence and identify if differences in gene content were due to mutations, annotation pipline differences, pseudogenes or inserted/deleted coding sequences. For strains Robertson, Fish Creek and Goon Nure a respective 28, 19 and 22 zero-hit blastp putative proteins matched regions of the Shintoku genome with tblastn. This resulted in a respective total of 1, 13 and 39 coding sequences from Robertson, Fish Creek and Goon Nure strains did not show similarity (e-value cutoff 1e-5) to the translated Shintoku genome sequence with either blastp or tblastn. Of the zero-hit blastp putative proteins that matched with tblastn 20 (Robertson), 8 (Fish Creek) and 11 (Goon Nure) matched intergenic regions of the Shintoku genome. A further 8, 10 and 10 predicted proteins of the respective Robertson, Fish Creek and Goon Nure genome sequences matched annotations in regions of hypothetical proteins in the Shintoku sequence, but predicted a protein in a different translational frame. Finally, one predicted hypothetical protein from both the Fish Creek and Goon Nure isolates matched the opposite strand of a proteasome component (XP_009692024.1) gene. Further inspection of this proteasome component revealed it is truncated by four C-terminal exons in Fish Creek and Goon Nure when compared to Shintoku and Robertson. However, this truncation is consistant with other piroplasmid species (T. parva, T. annulate, B. microti and B. bigemina) which are of comparable length to the predicted Fish Creek and Goon Nure proteasome component.

Zero-hit blastp putative proteins were also compared to the non-redundant protein database with BLAST (accessed 29–1-2017) to identify if any were homologous to proteins from other species. A total of 8, 7 and 10 of these putative proteins from the Robertson, Fish Creek and Goon Nure isolates respectively showed similarity to those found in the non-redundant database. These included a conserved apicomplexan specific protein (all isolates), a putative exonuclease (all isolates) and a RNA methyl-transferase (Robertson, Fish Creek). Results of blastp searches of the nr database are shown in Additional file 8.

Recombination analysis

Potential for recombination between T. orientalis genotypes has previously been postulated [1, 22]. To explore this, we examined the T. orientalis SNV datasets for recombination events as previously done with Theileria parva [23]. The frequencies of all 7 possible genotype combinations are shown as a Venn diagram (Fig. 1) and show high values for Fish Creek and Goon Nure containing combinations and low numbers for Robertson combinations, which reflect the genetic distance between these genotypes and the Shintoku reference. When these allele combinations were graphed against their genomic loci, little evidence of SNV clustering associated with recombination events could be observed (data not shown), potentially due to the high number of variant positions present. To further examine potential recombination, six recombination detection tests, namely Geneconv, MaxChi, Recombination Detection Program (RDP), BootScan, 3Seq and SiScan were used to analyse SNV-derived alignments. A total of 83 potential recombination events were detected in the concatenated SNV dataset (Fig. 2; detail in Additional file 9). All detected events were relatively small, with genomic sizes ranging 3–13,821 bp in size and a median of 989 bp. Furthermore, no evidence of recombination was detected between the Ikeda genotype samples (Robertson/Shintoku) and Chitose (Fish Creek) or Buffeli (Goon Nure) genotypes.

Mitochondrial genomes of Ikeda, Chitose and Buffeli isolates

Theileria mitochondrial genomes have been previously observed to be linear and relatively small at approximately 6 kbp in size [24]. In T. orientalis isolates that were examined in this study (Fig. 3), we were able to confirm a linear structure with inverted PCR utilizing outward facing primers at the edge of each mitochondrial genome. Furthermore, some Babesia mitochondrial genomes have been described to undergo multiple inversions [25]. To determine if this occurs in T. orientalis mitochondrial genomes, we examined sequencing reads mapped to assembled Robertson, Fish Creek and Goon Nure mitochondrial genomes for evidence of split-reads, i.e. reads partially mapped to the contig and soft-clipped. We found no evidence of soft-clipped reads indicative of an inversion event in any of the three T. orientalis mitochondrial sequences.

dN/dS analysis

MLST phylogeny of Ikeda, Chitose and Buffeli and their place within the Piroplasmida

To assess if whole genome sequences could be used to further refine the taxonomy of T. orientalis, we explored two commonly used methods of species definition, phylogeny and average nucleotide identity (ANI), using T. orientalis whole genome assemblies and representative whole genome sequences of the Order Piroplasmida (Fig. 4). Support values based on bootstrap analysis show very high support for all branches. Both whole genome phylogeny and ANI reveal a very close relationship between Shintoku and Robertson genomes. In contrast, ANIs between strains Robertson/Shintoku and Fish Creek/Goon Nure are low for organisms considered to be of the same species and compare well with those observed between T. parva and T. annulata (Fig. 4b). ANI clustering was also used to compare species of the Plasmodium genus (Additional file 11). For comparison, murine species Plasmodium berghei, Plasmodium yoelii, Plasmodium vinckei, and Plasmodium chabaudi have pairwise ANI that range from 83.5–90.1%.

Samples

Samples were sourced from cattle testing PCR positive for a single T. orientalis genotype based on the major piroplasm surface protein sequence (ie. Ikeda, Chitose or Buffeli) and were collected by private or district veterinarians as part of routine clinical monitoring. The isolates were sourced from cattle in Robertson, New South Wales (Collected 2009, breed: Angus), Fish Creek in Southern Victoria (Collected 2014, breed: Friesian) and Goon Nure in East Gippsland, Victoria (Collected 2012, breed: Freisian) respectively and named according to their location of isolation. Confirmation that each isolate was of a single genotype was performed with two separate PCR assays [10, 40]. Two to three days following the initial bleed, approximately 80 mL of blood was collected into anticoagulant (heparin or citrate) and shipped to a laboratory cold for propagation (1 passage) and extraction. On arrival samples were mixed with cryopreservative and stored as stabilates at − 80 °C until required.

Propagation and purification of Theileria piroplasms

This research was carried out in accordance with the Australian Code of Practice for the Care and Use of Animals for Scientific Purposes at the Tick Fever Centre, Wacol, Queensland. To propagate Theileria piroplasms, individual splenectomised calves were inoculated with stabilate generated from Robertson, Fish Creek and Goon Nure strains. Blood samples were drawn from each calf at regular intervals to monitor the infection and Giemsa-stained blood smears were used to estimate the level of parasitaemia. The packed cell volume (PCV) was also monitored. When parasitaemia had reached an appropriate level (6–20%), approximately 3.5 L of calf blood was collected into anticoagulant (CPDA1 or heparin). The blood was transferred into 300–500 mL blood bags using a dialysis pump and subsequently passed through Terumo leukocyte filters (300 mL blood per filter) under gravity. Leukocyte-depleted blood was then transferred to centrifuge tubes, centrifuged at 2500×g for 20 min and the serum and any remaining buffy coat removed. The remaining erythrocytes were washed 3 × with Dulbecco’s phosphate-buffered saline (D-PBS) followed by centrifugation as described above and diluted 5 × with D-PBS and loaded into a cell disruption vessel (“nitrogen bomb”). The erythrocytes were lysed and the parasite harvested by differential centrifugation as described previously [55]. Briefly, the vessel was infused with nitrogen gas to a pressure of 1000 psi for 1 min and then the pressure was released. The lysed erythrocytes were collected into a clean vessel and then transferred to centrifuge tubes. The cell lysate was centrifuged at 670×g for 10 min to pellet the red blood cell debris. The supernatant was then harvested and centrifuged at 2700×g for 10 min to pellet the piroplasms. The piroplasms were resuspended in D-PBS and a smear prepared from the suspension, which was stained with Giemsa stain to check purity. The piroplasms were then aliquoted into microfuge tubes and centrifuged at 2700×g for 10 min. The supernatant was removed and the piroplasm pellets stored at − 80 °C. Piroplasm preparations were transferred on dry ice and maintained at − 80 °C thereafter.

Nucleic acid extraction and whole genome sequencing

Piroplasm pellets were resuspended in 200 μL of phosphate-buffered saline and DNA was extracted using the DNeasy Blood and Tissue kit (Qiagen) according to the manufacturer’s instructions. For each sample, tagmentation of genomic DNA, and PCR amplification of tagged DNA were performed in triplicate using the Nextera system (Illumina). Sequencing libraries were pooled and normalized using bead size selection (SPRI beads, Beckman Coulter). The Agilent 2100 Bioanalyzer, with the High Sensitivity DNA kit was used to quantitate the pooled library before loading onto the Illumina platforms (Miseq or HiSeq). Paired-end 250 nt reads were generated using MiSeq V2 chemistry and paired-end 150 nt reads were generated using the HiSeq2500 system.

Protein extraction, electrophoresis and mass spectrometry

From approximately 0.1 g of T. orientalis Ikeda (Robertson strain) pelleted piroplasms, proteins were extracted, enriched, excised from 1-D sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) gels and prepared for mass spectrometry as previously described [56]. Briefly, proteins were separated by SDS-PAGE, stained, excised from the gel lane and separated into 16 slices. Each gel slice was then diced into ~ 1 mm cubes, stripped of dye, washed and digested in-gel with trypsin for analysis. Protein identification was performed by liquid chromatography with tandem mass spectrometry (LC-MS/MS) using a QSTAR Elite hybrid quadrupole time-of-flight mass spectrometer (AB Sciex). MS/MS data files were analysed using Mascot (provided by the Australian Proteomics Computational Facility, hosted by the Walter and Eliza Hall Institute for Medical Research Systems Biology Mascot Server) against the non-redundant MSPnr100 database compiled from all known reference protein sequences including NCBI, Refseq, UniProt, EuPathDB and Ensembl. Peptides matches at p < 0.05 (Ion score > 60) were classified as hits.

Assembly and annotation

Genome assembly was achieved with raw reads using the A5-miseq (v2015–05-22) pipeline [57]. Levels of remaining host DNA were assessed and filtered using BioBloomTools [58] and the B. taurus genome sequence [59]. Whole genome alignments were performed using progressiveMauve from Mauve v2.4.0 and nucmer from MUMmer v3.23 [60, 61]. For analysis of gene presence/absence, genes were annotated using Maker v3.00 beta [62] to combine the annotation tracks derived from expressed sequence tag (EST) and protein alignments and ab initio gene predictors. Firstly, repeat regions of the Robertson, Fish Creek and Goon Nure strain genome assemblies were masked using RepeatMasker v4.0.7. With repeat masked sequences, ab initio annotation tracks were derived from gene predictors Augustus v3.2.3 [63], Snap v2013–11-29 [64], and GeneMark ES Suite v4.32 [65] utilizing both GeneMark ET (Robertson and Fish Creek), and GeneMark ES (Goon Nure). A training set based on EST alignments (RNA data sourced from the Shintoku WGS project [15]) and transcoding regions was generated by PASA v2.1.0 [66] and used to train all ab initio predictors except GeneMark ES. Alignment of EST and protein evidence was performed with Exonerate v2.2.0 [67] using the coding2genome and protein2genome models respectively. Maker was run using the trained gene predictors and performed EST and protein alignments. Following annotation each independent gene track was assessed using Evidence Modeler v1.1.1 [66] to generate the final annotation with the following weights applied (EST: coding2genome = 7, blastn = 2, tblastx = 2. Protein: protein2genome = 10, blastx = 2. Ab initio: augustus = 2, snap = 1, genemark = 1). Additionally, annotations were manually curated by visual inspection and comparison to those transferred by RATT [68] for further quality checking. All parameter files used by Maker and Evidence Modeler to generate each annotation are included as Additional files 12, 13 and 14. Functional annotation was performed with Blast2GO v4.1.9 [69]. Predicted protein sequences were examined with blastp-fast (e-value cutoff 1e-5) against a local version of the non-redundant protein database (downloaded 2017–07-01). From blastp results, a gene ontology (GO) annotation was generated and extended through merging with results from InterProScan [70]. Finally, GO annotations were extended with ANNEX [71].

Ortholog clustering and gene presence analysis

Ortholog clustering of translated proteins from annotated draft genomes was performed with Orthofinder v2.1.2 [21] using 10 whole genome sequences from species representing the Piroplasmida order and one outgroup species Plasmodium falciparum strain 3D7. The 10 genomes include the three T. orientalis isolates examined in this study and previously published and annotated whole genome sequences of T. orientalis¸ Theileria parva, Theileria annulata, Theileria (formerly Babesia) equi, Babesia bovis, Babesia bigemina and Babesia microti [15–18, 27, 72–74]. For characterisation of the presence of T. orientalis genes in sequences generated by this study and absence in the Shintoku genome sequence, translated Robertson, Fish Creek and Goon Nure proteins were scanned against translated T. orientalis Shintoku proteins using blastp (e-value cutoff 1e-5). Robertson, Fish Creek and Goon Nure proteins that did not match any T. orientalis Shintoku proteins were further examined using tblastn (e-value cutoff 1e-5) of the T. orientalis Shintoku genome sequence. All proteins that did not match Shintoku sequences were further examined using blastp and the non-redundant protein database (2018–01–30).

Variant calling and validation

For variant calling, reads were first trimmed with Trimmomatic v0.33 based on leading and trailing base quality (Q < 20) and a 4-base sliding window when average quality scores were less than 20 [75]. Read quality was assessed pre- and post-trimming with FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and sequence files were manipulated using ngs-utils v0.5.7 [76]. For between population variants, trimmed reads were mapped to the T. orientalis Shintoku reference sequence [15] with NextGenMap v0.4.12 [77]. Low quality mapped reads (Q < 10) were filtered using samtools v1.2 [78] and alignments were subsequently sorted, and duplicates removed using picard tools v1.138 (http://broadinstitute.github.io/picard/). Variant calling was performed using VarScan v2.3.8 [20] on mpileup files generated by samtools. Variant calling parameters were based on previous studies into apicomplexan SNV detection and variant calling algorithm comparisons [79, 80]. Specifically, variants were called using a minimum coverage of 14 × [79], minimum average quality of 20, minimum variant supporting reads of 4, minimum variant frequency of 0.01 and minimum p-value threshold of 0.05. A minimum frequency of 0.9 was used to define homozygous variants. Variant quality was assessed by examining subsets of variant sequences with Sanger sequencing and by using simulated alignments. Primers for Sanger sequencing experiments are shown in Additional file 15. Sanger sequencing was performed on PCR products at the Australian Genome Research Facility. For simulated alignments, the T. orientalis Shintoku sequence was modified by the msbar program of EMBOSS [81] to contain approximately 90,000 point, codon and block insertion/deletion events and 900,000 substitution events. Simulated reads were generated from simulated genomes using ART vMountRainier 2016–06-05 [82] and mixed into one fastq file. False discovery rate was calculated by aligning simulated to experimentally determined sequences and comparing true and detected variants.

Recombination and selection analysis

Recombination analysis was performed on concatenated SNV alignments to allow direct comparison to previous attempts at examining recombination in Theileria spp. [23]. These methods are potentially susceptible to error from highly divergent sequences (< 60% identity) and reduced SNV density information. To limit false positive results, we performed recombination analysis using six recombination tests, namely Geneconv, MaxChi, RDP, BootScan, 3Seq and SiScan as previously described with the exception that, due to higher SNV density, SNVs were not discounted if they were within 100 bp of another SNV [23]. For dN/dS analysis, SNV-containing genomes were constructed using homogeneous SNVs from Fish Creek and Goon Nure strain variant analysis with the FastaAlternateReferenceMaker tool from the Genome Analysis Toolkit v3.2.2 [83]. Coding sequence alignments of Reference and SNV-containing genomes were generated for gene regions of 14 × or greater sequencing coverage using customised python scripts and biopython [84]. Coding DNA sequences (CDS) were excluded where sequencing coverage at > 14 × made up less than 50% of the CDS. dN/dS ratios and synonymous/non-synonymous mutations were calculated for coding sequences using KaKs Calculator v2.0 [85]. Predicted surface proteins were identified using SignalP 4.0 [86] and TMHMM [87].

Piroplasmida phylogenomic analysis

Piroplasmida phylogeny was examined using single copy orthogroups identified during ortholog clustering (outlined above). A total of 654 orthologous groups with single copy genes found in all 11 genomes were identified using Orthofinder v2.1.2 [21], with initial search results filtered to retain only groups with an e-value <1e-30. Sequences from each of the 654 orthologous groups were aligned and trimmed using MAFFT v7.310 and trimAl v1.3 [88, 89] as described previously [27] with the following alterations. Sequence alignments were trimmed using the nogaps automated option and two orthologous groups were removed from the final analysis after trimming due to poor alignment and high trimming. A supermatrix alignment (261,036 residues) was generated by concatenating individual gene alignments using a customised python script. Protein evolution model selection was performed for each gene of the supermatrix alignment using Prottest v3.4 [90]. The supermatrix tree was inferred with RAxML v8.2.4 as previously described [27, 91]. Average nucleotide identity comparisons were calculated using pyani v0.2.3 (https://github.com/widdowquinn/pyani) using the ANIb method described in [92].