The Plasmodium berghei RC strain is highly diverged and harbors putatively novel drug resistance variants

Four-day suppressive test for in vivo drug sensitivity

BALB/c mice were injected intravenously with 1 × 10⁷ P. berghei infected red blood cells. For PbRC, five to six mice were used per group. For PbANKA, four to eight mice were used per group. Oral doses of artesunate were given at 4, 24, 48 and 72 h post infection. Four days post infection, parasitemia was determined by manually counting infected and uninfected red blood cells in Giemsa-stained thin blood smears. Percent inhibition was calculated using the following formula: $$\text{Percent inhibition}=100-((100\times \text{parasitemia of each dose})/\text{parasitemia of untreated control})$$

Data of percent inhibition at different doses of artesunate were fitted to the two-parameter sigmoidal dose response equation using the drc package in R (Ritz & Streibig, 2005). This study was carried out in accordance with the guidelines in the Guide for the Care and Use of Laboratory Animals of the National Research Council, Thailand. All animal experiments were performed with the approval of BIOTEC’s Institutional Animal Care and Use Committee (Permit number BT-Animal 02/2557). At the end of the experiments, mice were euthanized by CO₂ asphyxiation. All efforts were made to alleviate pain and suffering.

Sanger dideoxy sequencing of PCR products from selected regions

Selected genomic regions with candidate variants in the P. berghei RC strain were PCR-amplified using high-fidelity Phusion DNA polymerase (New England Biolabs, Ipswich, MA, USA) using primers listed in Table S1. The resulting PCR products were sent for Sanger dideoxy sequencing (1st BASE, Selangor, Malaysia).

Whole genome DNA sequencing

Parasitized blood obtained from a single mouse infected with the PbRC parasite was passed through a CF11 (Whatman/Sigma Aldrich, St. Louis, MO, USA) column to remove white blood cells. Cells were harvested by centrifugation and parasites were liberated from red cells by lysis of the red blood cell membrane with 0.2% saponin. Parasites were washed twice with phosphate buffered saline (137 mM NaCl, 2.7 mM KCl, 4.3 mM Na₂HPO₄, 1.47 mM KH₂PO₄, pH 7.4), resuspended in lysis buffer (8.5 mM Tris–HCl, 342 mM NaCl, 2 mM Na₂-EDTA, 0.14 mg/ml proteinase K, 0.14% SDS, pH 8.2) and incubated at 37 °C overnight. Proteins were precipitated by the addition of NaCl to 1.5 M, centrifuged and discarded. Isopropanol was added to the supernatant to precipitate nucleic acids. The nucleic acid pellet was dried and resuspended in TE buffer (10 mM Tris, 1 mM EDTA, pH 8.0). RNA was removed by digestion with RNAseA. Phenol chloroform extraction and ethanol precipitation of DNA were performed. The integrity of genomic DNA was checked by agarose gel electrophoresis. Genomic DNA was submitted to the Chulalongkorn Medical Research Center, Bangkok, Thailand for genome sequencing. Genomic DNA was sheared by sonication and DNA fragments <1 kb were gel-purified. Sequencing libraries were constructed from sheared genomic DNA using a TruSeq kit (Illumina, San Diego, USA). About 2 × 150 bp reads were obtained using a MiSeq instrument (Illumina). The raw data are deposited in the NCBI Sequence Read Archive, accession number PRJNA277169.

Sequence data analysis

Raw sequencing data were obtained in FASTQ format. Genome sequence data for strains K173, NK65NY, NK65E, SP11_RLL, SP11_A (Otto et al., 2014) were downloaded from the European Nucleotide Archive. Filtering of raw data to remove poor quality reads (average Q-score <30) and removal of adapter sequences were performed using FASTQC (Andrews, 2010). Preprocessed read data were aligned to the P. berghei ANKA version 3 reference genome (downloaded from GeneDB; Logan-Klumpler et al., 2012) using BOWTIE version 2.2.2.6 software (Langmead et al., 2009) under the setting of length of seed substring 22 bases without clipping. GATK software version 3.3.0 (McKenna et al., 2010) was used to identify potential base-substitution single nucleotide variants (SNV) and small insertion/deletions (INDEL) between all strains different from the reference strain by the HaplotypeCaller method combined with the GenotypeGVCFs method, under the settings: --minReadPerAlignmentStart (10), --min_base_quality_score (10), --sample_ploidy(2), --heterozygosity(0.001). SnpEff software version 4.3g (Cingolani et al., 2012) was used to annotate the variants using the genome annotation file downloaded from GeneDB. The .vcf file containing all raw variants called by GATK is provided in Data S1. The raw variants were filtered to retain only high-confidence SNVs for genome analyses. INDEL variants were not included, since many occur within simple sequence repeats that could be prone to PCR and sequencing artifacts. Potential false positive SNVs resulting from sequencing artifacts were removed by excluding variants with heterozygous or missing genotypic calls in any strain. Potential false variants resulting from read alignment error in repetitive regions were removed using the DustMasker program (Morgulis et al., 2006). Variants present among multigene families listed in Additional File 4 of Otto et al. (2014) were also treated as potential false positives and were removed. After filtering, 8,681 SNV markers remained.

Principal components analysis (PCA) was performed using the ipPCA tool (Limpiti et al., 2011), implemented in MATLAB version R2009b. The 8,681 filtered SNV markers were used as input for ipPCA. Tajima’s D scores were calculated using PopGenome, a population genomic analysis tool in R (Pfeifer et al., 2014). Scores were calculated in sliding non-overlapping genomic windows of five consecutive variants (1,730 windows in total), and separately for annotated genes (1,765 genes in total, Table S2). Plots of Tajima’s D scores were made using the ggplot2 package in R (Wickham, 2009). Gene ontology analysis of genes with Tajima D scores greater than 1 was performed using the Gene ontology web service provided in the PlasmoDB website (Aurrecoechea et al., 2009). Terms were considered significant using a Bonferroni-corrected p-value threshold of 0.05.

Copy number variants (CNVs) were identified from the whole-genome sequencing data using the cn.MOPS tool implemented in R (Klambauer et al., 2012). The sequence data for all six strains aligned to the reference genome were used as input for cn.MOPS, which was run using the haploid genome setting and window size of 500 bp used for local modelling of read counts. Mapped reads in genomic regions with putative variants were visualized using the Integrative Genomics Viewer program (Thorvaldsdottir, Robinson & Mesirov, 2013).

PbRC is widely diverged from other P. berghei strains

Illumina sequencing of PbRC genomic DNA was performed. Sequence reads were obtained from 17,801,263 clusters, and 80.5% of the preprocessed reads could be mapped to the PbANKA v3 reference genome. Genomic sequence data of previously sequenced strains K173, NK65NY, NK65E, SP11_RLL and SP11_A (Otto et al., 2014) were aligned to the reference genome using the same analytical procedure. SNV and small INDEL variants were called from the mapped reads for all strains. A total of 27,495 variant markers were identified among the six strains (Data S1), of which 8,681 SNV markers remained after applying stringent filtering criteria. The numbers of variants identified in each strain are shown in Table 1.

The majority of filtered SNVs are private to each strain, and the PbRC strain shows the highest number of variants, indicating that it is markedly different from the other strains. This is confirmed by PCA in which PbRC is clearly separated from the other strains in the first principal component, which captures most of the genotypic variance among these strains (Fig. 2).

Genetic diversity across the P. berghei genome

A previous analysis of P. berghei genotypic diversity revealed a low SNV frequency compared with P. chabaudi, another rodent malaria species (Otto et al., 2014). However, the extreme divergence of PbRC from other strains may provide further insight into P. berghei intraspecific genetic diversity. The Tajima’s D score was calculated for sliding genomic windows along all chromosomes, and also for annotated genes. Tajima’s D is a measure of the observed versus expected genetic diversity, calculated as the difference between the mean number of pairwise differences and the number segregating sites (Tajima, 1989). The mean Tajima’s D is negative for each chromosome, indicating a genome-wide pattern of negative Tajima’s D (Fig. 3A). However, there are several regions with positive Tajima’s D values that are spread throughout the genome (Fig. 3B). At the gene level, 122 genes were identified with a Tajima’s D score greater than one (Table S2). No specific gene ontology terms or biological processes are significantly enriched among these genes with high Tajima’s D scores. However, a few genes homologous to Plasmodium genes with known functions in host cell invasion showed high Tajima’s D scores, including PBANKA_132170 (berghepain 1), PBANKA_041290 (circumsporozoite- and TRAP-related protein, PbCTRP), PBANKA_1115300 (glideosome-associated protein 40, PbGAP40) and PBANKA_1137800 (glideosome-associated connector, PbGAC).

Next, variants private to PbRC located within genes were examined to identify putative causal variants of the drug resistant trait in this strain. To our knowledge, the only prior evidence of PbRC genetic variants compared with the PbANKA reference is an Ile to Lys substitution at residue 413 (I413K) of the gamma glutamylcysteine synthetase (Pbγ-gcs) gene (Pérez-Rosado et al., 2002). We confirmed this mutation in PbRC from whole genome and Sanger sequencing. Comparison of P. berghei strain genomes identified this variant as private to PbRC (Table 2). A previous study reported a possible translocation of the multidrug resistance associated protein (MRP) gene to chromosome 8 in PbRC (Gonzalez-Pons et al., 2009). However, visualization of PbRC mapped reads showed contiguity in chromosome 14 in the vicinity of the MRP gene reference location, and thus no evidence of translocation (Fig. S1).

We searched for candidate drug resistance variants in genes known to modulate CQ and/or ART drug resistance in other Plasmodium spp. The variants in these genes private to PbRC are shown in Table 2. Mutations in the P. chabaudi ubiquitin carboxyl-terminal hydrolase 1 (ubp1) gene modulate CQ and ART resistance in laboratory-selected drug resistant parasites (Hunt et al., 2007; Henriques et al., 2013). Five non-synonymous variants private to PbRC were found in the homologous Pbubp1 gene. Mutations in the CQ resistance transporter (crt) gene modulate CQ resistance in P. falciparum (Martin & Kirk, 2004; Ecker et al., 2012). A non-synonymous V42F variant private to PbRC was found in the P. berghei homologue (Pbcrt). Mutation of the P. falciparum multidrug resistance associated protein 1 (mdr1) gene modulates sensitivity to several antimalarials, including CQ and ART (Sanchez et al., 2010). A non-synonymous V54A variant private to PbRC was found in the homologous Pbmdr1 gene. Among other proteins implicated as modulators of ART sensitivity in P. falciparum, a non-synonymous F320C variant private to PbRC was found in the Pbpi3k gene homologous to phosphatidylinositol-3-kinase (PfPI3K) (Mbengue et al., 2015). PbRC mutations were notably absent from some genes strongly implicated as modulators of ART sensitivity, including the homologues of μ subunit of adaptor protein 2 complex (AP2-μ; PBANKA_1433900), a gene mutated in the P. chabaudi CQ- and ART-resistant AS-ART strain (Henriques et al., 2013) and K13 Kelch propeller (PbK13; PBANKA_1356700), mutations of which confer reduced ART sensitivity in P. falciparum (Ariey et al., 2013; Ghorbal et al., 2014).

Catabolism of hemoglobin is tied to the mechanisms of action of both CQ and ART. CQ prevents crystallization of heme, a product of hemoglobin catabolism (Martin & Kirk, 2004), whereas the antimalarial effect of ART is strongly dependent on catabolism of hemoglobin (Klonis et al., 2011; Xie et al., 2015). Given that the CQ and ART-resistant PbRC is defective for hemoglobin catabolism and hemozoin formation (Peters, 1964; Peters, Fletcher & Staeubli, 1965), we investigated whether PbRC harbored private variants in genes encoding proteins known to function in catabolism of hemoglobin (Ponsuwanna et al., 2016; Lin et al., 2015). Among these genes, a non-synonymous variant private to PbRC was found only in the dipeptidyl aminopeptidase 3 gene (Pbdpap3) (Table 2).

In addition to genes previously associated with CQ and ART-resistance, PbRC private variants in other genes may be responsible for the traits associated with this strain. PbRC does not produce gametocytes (Peters, Fletcher & Staeubli, 1965), and so may harbor variants in genes important for gametocytogenesis. PbRC is also reported as slow growing compared with the drug-sensitive N strain (Peters, 1964), and so may harbor variants in genes important for growth. Furthermore, the ART resistance trait of PbRC may involve mutations in several genes, since ART resistance in P. falciparum is a heritable trait in which the expression patterns of many genes are changed (Mok et al., 2015). PbRC private variants were identified in genes mutated in gametocyte non-producing lines derived from the reference strain PbANKA (Sinha et al., 2014), genes with annotated function with respect to blood-stage growth (Bushell et al., 2017) and P. berghei orthologs of genes with altered expression in ART-resistant P. falciparum (Mok et al., 2015). The summaries of PbRC private variants for each phenotypic category are shown in Table 3, and all variants in each phenotypic category are shown in Table S3. Non-synonymous variants are more prevalent than synonymous for each phenotypic category, including essential genes.

Finally, drug resistance in Plasmodium spp. is often modulated by CNVs; for instance, the P. falciparum gene pfmdr1 is amplified in geographical regions where antimalarial resistance is common (Nair et al., 2006). CNVs among P. berghei strains were identified from whole-genome sequencing data (Table S4). All of the CNVs are present in chromosomal regions near chromosome ends. All of the genes that overlap CNVs are members of multigene families. Amplified regions with copy number greater than one are not private to any strain, including PbRC. Large deletions (>1 kb) are prevalent in all strains except PbSP11_A, and PbRC has the greatest number of large deletions.