Genomic analysis of Isometamidium Chloride resistance in Trypanosoma congolense.

Isometamidium Chloride (ISM) is one of the principal drugs used to counteract Trypanosoma congolense infection in livestock, both as a prophylactic as well as a curative treatment. However, numerous cases of ISM resistance have been reported in different African regions, representing a significant constraint in the battle against Animal African Trypanosomiasis. In order to identify genetic signatures associated with ISM resistance in T. congolense, the sensitive strain MSOROM7 was selected for induction of ISM resistance in a murine host. Administered ISM concentrations in immune-suppressed mice were gradually increased from 0.001 mg/kg to 1 mg/kg, the maximal dose used in livestock. As a result, three independent MSOROM7 lines acquired full resistance to this concentration after five months of induction, and retained this full resistant phenotype following a six months period without drug pressure. In contrast, parasites did not acquire ISM resistance in immune-competent animals, even after more than two years under ISM pressure, suggesting that the development of full ISM resistance is strongly enhanced when the host immune response is compromised. Genomic analyses comparing the ISM resistant lines with the parental sensitive line identified shifts in read depth at heterozygous loci in genes coding for different transporters and transmembrane products, and several of these shifts were also found within natural ISM resistant isolates. These findings suggested that the transport and accumulation of ISM inside the resistant parasites may be modified, which was confirmed by flow cytometry and ex vivo ISM uptake assays that showed a decrease in the accumulation of ISM in the resistant parasites.

Introduction

African trypanosomes are protozoan parasites transmitted by tsetse flies and responsible for African Trypanosomiasis in both humans (HAT) and animals (AAT). T. congolense, T. vivax and T. brucei brucei are the major pathogens responsible for AAT, a disease also known as Nagana and considered as one of the principal causes of malnutrition and poverty in 37 countries in sub-Saharan Africa (Mattioli et al., 2004). Nagana puts 55 million cattle at risk and leads to the death of three million animals every year, inflicting a direct annual loss of US$ 1.0–1.2 billion in cattle production (Hursey and Slingenbergh, 1995, Mattioli et al., 2004, Cecchi and Mattioli, 2009, Cecchi et al., 2014).

Isometamidium Chloride (ISM) and Diminazene Aceturate (DA) are the main drugs currently available to cure AAT (Giordani et al., 2016). ISM belongs to the phenanthridine drug group of which many members are documented to be mutagenic (Uilenberg, 1998), and DA is an aromatic diamidine. About 35 million doses of trypanocidal drugs are administered every year in an attempt to control the disease (Geerts and Holmes, 1998, Mattioli et al., 2004). However, those compounds have been in use for over 50 years (Lourie and Yorke, 1939, Watkins and Woolfe, 1952, Brown et al., 1961, Sahin et al., 2014), and numerous cases of ISM and DA resistance and treatment failures have been reported in 17 African countries, representing a serious problem in the battle against AAT (Sinyangwe et al., 2004, Delespaux et al., 2008a, Delespaux et al., 2008b, Moti et al., 2012, Mungube et al., 2012, Sow et al., 2012, Dagnachew et al., 2015, Moti et al., 2015).

ISM and DA accumulate inside the cytoplasm and mitochondrion of trypanosomes and have a common target: the mitochondrial DNA (Shapiro and Englund, 1990). ISM also enters the nucleus where it binds to the nuclear DNA (Kaminsky et al., 1997). The mitochondrial DNA, known as kinetoplast DNA (kDNA), represents about 10–20% of the genome of the trypanosomes (Daniels et al., 2010) and consists of a giant network of interlocked DNA rings (Liu et al., 2005). ISM selectively inhibits the kDNA type II topoisomerase, a protein involved in the replication of the mitochondrial DNA, which leads to the shrinking and disappearance of the kDNA network, and ultimately to cell death (Shapiro and Englund, 1990, Wang and Englund, 2001, Roy Chowdhury et al., 2010). DA has high affinity for A-T base pairs and binds to circular double strand DNA to induce complete and irreversible loss of the kDNA network (Brack et al., 1972, Mahler and Perlman, 1973, Gonzalez et al., 1997).

While the mechanisms of ISM and DA resistance have been described in few trypanosome species (Brun and Lun, 1994, Neupert, 1997, Schnaufer et al., 2005, Dean et al., 2013, Gould and Schnaufer, 2014, Baker et al., 2015), little is known for T. congolense. Here, the identification of a GAA insertion in a gene coding for an ATP-binding cassette (ABC) transporter (i.e. energy-dependent transport of various substrates across biological membranes) in a range of ISM resistant parasites led to the hypothesis of a potential involvement of ABC transporters in ISM resistance (Delespaux et al., 2005). However, this GAA insertion was not shared by all ISM resistant T. congolense strains, suggesting that other mechanisms of drug resistance acquisition exist. The mitochondrial membrane potential (ΔΨm) was also suggested to be an important determinant of the rate and accumulation of ISM inside the mitochondrion, with a diminished ΔΨm resulting in a decreased ISM accumulation inside the mitochondrion and drug resistance (Wilkes et al., 1997). Active extrusion by plasma membrane transporters has also been proposed to contribute to the intracellular ISM accumulation level. Here, significant differences in drug transport were suggested to be correlated with the resistant phenotype (Sutherland et al., 1992).

To facilitate the identification of genomic determinants and mechanisms associated with drug resistance, Whole Genome Sequencing (WGS) technologies have been used in numerous studies, including some on protozoan parasites such as Leishmania spp (Downing et al., 2011, Mondelaers et al., 2016) and Plasmodium falciparum (Miotto et al., 2015) allowing the successful identification of the genomic signatures of drug resistance in these pathogenic organisms. Therefore, we have used this WGS methodology to investigate the genomic changes that occur during the experimental induction of ISM resistance in the initially sensitive MSOROM7 T. congolense strain. In three independent experiments, immune-suppressed mice were treated with increasing doses of ISM until the parasites acquired resistance to 1 mg/kg ISM, the maximal dose used in livestock. Flow cytometry and ex vivo uptake/extrusion assays revealed a reduced accumulation of ISM in the three independently induced resistant parasite lines. WGS analysis comparing these ISM resistant lines with the sensitive parental line identified allele frequency shifts in a variety of genes of which some were coding for different transporters and transmembrane proteins. In complement, we performed a similar WGS comparison of 50 T. congolense field isolates with a documented ISM susceptibility profile.

Methods

Ethics statement

The animal work was approved by the Institute of Tropical Medicine (ITM) animal ethics committee, clearance numbers DG020 and VP2014-2.

In vivo drug susceptibility testing

To evaluate the susceptibility of a T. congolense strain to ISM, six mice were injected with the strain of interest at a concentration of 105 parasites (bloodstream forms)/mouse. Twenty-four hours later, each mouse was treated with a single ISM dose of 1 mg/kg. The trypanosome presence in these mice was monitored for two months by microscopical examination and PCR analysis using the primers and protocol described in (Cox et al., 2005). At the end of the experiment a strain was considered sensitive to 1 mg/kg ISM when all mice (6/6) were cured, and was considered resistant when five or six relapses (reappearance of parasitemia) were observed. A strain showed partial resistance (intermediate profile) when one to four mice relapsed following drug administration (Eisler et al., 2001).

Effect of ISM pressure on the genome of T. congolense

To obtain insights into the genomic consequences of an ISM drug pressure on a parasite cell population, we maintained the ISM resistant isolate TRT15 in OF-1 mice for 10 passages under two different conditions: without drug pressure versus with a continuous ISM pressure. TRT15 is a T. congolense isolate that was originally collected from cattle in 1996 in Zambia (Delespaux et al., 2008a) (Table S1). To start the experiment, mice were infected by intraperitoneal injection with 105 parasites/mouse. When the level of parasitemia in the blood reached 7.9 × 106 parasites/ml, the mouse was either treated with 1 mg/kg ISM or left untreated. When parasitemia reached 1.3 × 108 parasites/ml, blood was collected and purified using an anion exchanger Diethylaminoethyl (DEAE)-cellulose resin suspended in Phosphate Saline Glucose (PSG) buffer, pH 8 (Lanham and Godfrey, 1970) to obtain parasites for a new infection (passage). Parasites were collected at the start of the experiment and after passages 5 and 10 for DNA extraction and sequencing (see below).

ISM resistance induction experiment

The sensitive T. congolense strain KTT/MSOROM7C1 (further referred as MSOROM7) was chosen to induce ISM resistance in immune-suppressed OF-1 mice (Fig. 1). MSOROM7 is a clone derived from the Zambian isolate KTT/MSOROM7 that was obtained from cattle in 2003 (Masumu et al., 2009). This clone was generated by infecting an OF-1 mouse with a single parasite using the micro-drop method (Masumu et al., 2009). Then, the parasites were passaged only once before the start of the ISM resistance induction experiment. Therefore, we can assume that the parasites were isogenic at the start of this experiment.

Fig. 1

Induction of ISM resistance, experimental procedure. One day prior to infection with the sensitive strain MSOROM7, mice were injected with Endoxan® (D0). When parasitemia reached the concentration of 7.9 × 106/ml (D5), mice were treated with 0.001 mg/kg ISM. Mice became microscopically non-parasitemic but parasitemia reappeared within a few days (Dx), and mice were treated with an increased dose of ISM (0.005 mg/kg ISM). When parasitemia reappeared and reached 1.3–5x108/ml (Dy), mice were euthanized and parasites were collected to inject a new group of naïve mice that were first treated with 0.005 mg/kg ISM, and then with 0.01 mg/kg ISM. This cycle was repeated using increasing doses of ISM (0.02; 0.04 and 1 mg/kg) until reaching resistance to 1 mg/kg ISM.

One day prior to infection (D0), three naïve mice were injected intraperitoneally with 200 mg/kg cyclophosphamide (Endoxan®, Baxter) in order to obtain immune-suppression. The next day (D1), these mice were infected intraperitoneally with a starting dose of 105 bloodstream form trypanosomes. When the level of parasitemia in the blood reached 7.9 × 106 trypanosomes/ml, mice were treated with 0.001 mg/kg ISM (i.e. first ISM treatment) and became microscopically non-parasitemic within a few days. When parasitemia in blood reappeared and reached 7.9 × 106 trypanosomes/ml (DX), mice were injected with an increased dose of ISM (0.005 mg/kg). As soon as the parasitemia reappeared and reached a level between 1.3 × 108 and 5 × 108 trypanosomes/ml (i.e. peak of parasitemia; DY), blood from each individual mouse was collected on a DEAE resin for parasite isolation and DNA extraction. A part of the collected blood was also used to infect a new group of three Endoxan® treated mice with a starting dose of 105 trypanosomes. These newly infected mice received two ISM injections: a first injection with 0.005 mg/kg ISM, and after parasites reappearance, a second injection with 0.01 mg/kg ISM (Ndoutamia et al., 1993). This infection-ISM treatment cycle was repeated several times with increasing doses (0.02; 0.04 and 1 mg/kg ISM) until full resistance to 1 mg/kg ISM was acquired. Using this methodology, three independent ISM resistant lines were obtained: MSOROM7-G; MSOROM7-H and MSOROM7-I. After each treatment cycle (annotated as passages one, two, three, four and five; or referred to as MSOROM7-G2, -G3, -G4, -G; MSOROM7-H2; -H3, -H4, -H; and MSOROM7-I2, -I3, -I4, -I; Table S1) parasites were collected for DNA extraction and sequencing (see below).

A similar resistance induction experiment in immune-competent OF-1 mice (i.e. no Endoxan® treatment) failed to give rise to fully resistant parasites.

Bioinformatic and statistical analyses

DNA extraction, sequencing and bioinformatic analyses were performed as described earlier (Tihon et al., 2017). In short, reads were mapped against the reference genome T. congolense IL3000 (Jackson et al., 2012) with SMALT v7.4 (http://sourceforge.net/projects/smalt/), and duplicate reads were filtered with SAMtools V1.0 (Li et al., 2009), resulting in an average read depth of 30.9× per TRT15 sample and 33.4× per MSOROM7 sample. Single-Nucleotide-Polymorphisms (SNPs) and insertions/deletions smaller than 15 bp in length (indels) were called using the Genome Analysis Toolkit v3.4 (GATK) (McKenna et al., 2010, DePristo et al., 2011) and annotated using SnpEff (Cingolani et al., 2012) based on the gene annotation file of T. congolense IL3000 (http://tritrypdb.org/common/downloads/release-6.0/TcongolenseIL3000/gff/data/). Gene deletions were measured by both the normalized read depth and the fraction of the gene covered by reads (length). We considered a gene to be totally deleted (i.e. gene deletion) when 75% or more of its depth coverage was missing. In the depth analysis, we used MSOROM7-G4 (i.e. mean depth coverage of 20×) instead of MSOROM7-G (i.e. mean depth coverage of 13×) to estimate long deletions in samples with sufficient depth coverage, and this choice was made because MSOROM7-G4 was as ISM resistant as MSOROM7-G (see results). Gene duplications were measured solely by the normalized read depth. To provide statistical significance for the gene duplication analysis, we calculated the differences between the sensitive and resistant lines and calculated the standard z scores. We converted the standard z scores to p-values using python scipy survival function. See (Tihon et al., 2017) for details on all computational parameters and filtering criteria for each of the above mentioned computational analyses.

We checked whether novel SNPs and small indels appeared in TRT15 and MSOROM7 over the course of the experiments. As this was not the case (see results), we checked whether there was a change in allele frequencies at heterozygous loci between the different passages of TRT15 and MSOROM7. To this end, we calculated the proportion of reads that support a given allele at a given SNP site, providing an estimate of the read depth allele frequency. For TRT15 we compared the read depth allele frequencies between passage 0 (start of the experiment) and passages 5 (1 mg/kg ISM) and 10 (1 mg/kg ISM). For MSOROM7 we compared the read depth allele frequencies between passage 0 and passages 2 (0.01 mg/kg ISM), 3 (0.02 mg/kg ISM), 4 (0.04 mg/kg ISM) and 5 (1 mg/kg ISM). If a sample at a given passage did not have sufficient read depth coverage (i.e. less than 10× coverage) at a given SNP site, then this site was removed from further analysis. In order to measure statistical significance we calculated standard deviations σ and a corresponding z-score of the differences in read depth allele frequencies between the passages as mentioned above. Finally, we only retained those sites for which the read depth allele frequency shifted of at least 0.25 and for which the z-score was at least 3. Shifts in read depth allele frequencies higher than 0.8 were considered homozygous, which was similar to the definition of (Nielsen et al., 2011) where a homozygous genotype was described as a proportion of an allele of less than 0.2 or larger than 0.8.

A similar approach was used to identify changes in allele frequencies between drug resistant and drug sensitive isolates from natural populations. To this end, we used published sequence data from 50 T. congolense isolates from 10 different countries in sub-Saharan Africa (Tihon et al., 2017). Population genomics analyses revealed two clusters: the TC1 cluster contained 21 parasites from The Gambia, Mali, Burkina Faso, Togo, Cameroon, Ethiopia and Uganda and the TC2 cluster contained three parasites from Congo, Tanzania and Zambia (Tihon et al., 2017). The remaining 25 parasites from Zambia and one parasite from Burkina Faso were of uncertain ancestry, showing signatures of introgressive hybridization (Tihon et al., 2017). The ISM sensitivity of all these isolates were determined in previous studies (Geigy and Kauffmann, 1973, Pinder and Authie, 1984, Clausen et al., 1992, Afewerk et al., 2000, Knoppe, 2002, Delespaux et al., 2005, Mamoudou et al., 2006, Mamoudou et al., 2008, Masumu et al., 2009, Moti et al., 2012): 30 strains were resistant to 1 mg/kg ISM, 14 were sensitive to 1 mg/kg ISM and six showed an intermediate profile with partial resistance to 1 mg/kg ISM (Table S1). Within each group (i.e. TC1, TC2 and Zambia) we compared the averaged read depth allele frequencies between resistant and sensitive isolates. Differences in the averaged allele frequencies ≥0.8 between the resistant and sensitive strains were classified as homozygous. Differences in the allele frequencies greater than 0.25 (but lower than 0.8) were classified as heterozygous. We selected genes with the averaged allele frequency shift higher than 0.25 between resistance and sensitive strains with a p value ≤ 0.05 (Mann-Whitney-Wilcoxon test). For a more targeted analysis, we focused on a set of 185 genes including some playing critical roles in drug resistance in various pathogens (Table S6) to further investigate allele frequency shifts between sensitive and resistant field isolates. This analysis identified 217 alleles in a total of 54 genes whose averaged allele frequency shifted ≥0.25 between the sensitive and resistant parasites of each population (Table S7).

The dataset generated and analyzed during the current study is available in EBI under the accession number PRJEB21066, http://www.ebi.ac.uk/ena/data/view/PRJEB21066.

ISM accumulation and extrusion

The accumulation of ISM in the sensitive and resistant T. congolense MSOROM7 lines as well as in the ISM resistant field isolate TRT1 was first monitored by flow cytometry, using the fluorescent properties of the drug (ISM: λex: 374 nm; λem: 590 nm). For an unknown reason, the field isolate TRT15 was very sensitive to the DEAE purification and many cells were lost during the isolation procedure, precluding the collection of a sufficient amount of parasites to perform the experiment. To overcome this issue, we performed the flow cytometric analysis on the ISM resistant field isolate TRT1 that is genetically nearly identical to the TRT15 strain (Tihon et al., 2017).

DEAE-column purified parasites from infected mouse blood were maintained in a cold PSG buffer pH 8. For each ISM condition, 106 trypanosomes were suspended in ice-cold PSG containing 0–10 μg/ml ISM supplemented with the Fixable Viability Stain (FVS, BD Horizon - λex: 649 nm; λem: 660 nm) and incubated for 30 min at 4 °C in the dark. Cells were washed twice in PSG and were analyzed with a BD FACSVerse flow cytometer and the BD FACSuite software v1.0.3.

The accumulation and extrusion of ISM in the MSOROM7 sensitive and resistant lines (except MSOROM7-H for which no sufficient amount of parasites could be harvested) and in the field isolate TRT1 were also measured following the protocol described in (Eze et al., 2016). The ISM fluorescence was measured using a VICTOR Multilabel Plate Reader (PerkinElmer) (ISM: λex: 374 nm; λem: 590 nm) and quantified using a standard curve of known ISM concentrations (i.e. the fluorescence of known concentrations of ISM (i.e. ISM in assay buffer) were determined using the VICTOR Multilabel Plate Reader to make a standard curve of known concentrations).

Real time qPCR

We identified allele frequency shifts in the gene TcIL3000_8_1200 coding for the drug/metabolite transporter (DMT) in the resistant MSOROM7-H line compared with the sensitive parental line (see results). In addition, the DMT gene had two copies in the sensitive parasites and one of these copies disappeared in the resistant lines (see results). To identify a potential link between allele frequency shifts, copy number variations and gene expression level, the expression level of the DMT gene TcIL3000_8_1200 was evaluated in MSOROM7 and in the resistant MSOROM7 lines challenged or not with 1 mg/kg ISM. PCR primers were designed with Primer-BLAST (www.ncbi.nlm.nih.gov/tools/primer-blast/) and evaluated with OligoAnalyzer (https://eu.idtdna.com/calc/analyzer). Each primer was blasted against trypanosome genomes (http://blast.ncbi.nlm.nih.gov/Blast.cgi) to evaluate their specificity. RNA extractions were conducted with TRIzol Reagent (Invitrogen), and the cDNA synthesis was performed using the SensiMix™ SYBR® No-ROX Kit (Bioline) according to manufacturer's instructions. The qPCR reactions were conducted on a LightCycler 480 (Roche) with SYBRGreen I master mix (Qiagen) and 500 nM of each primer. The cycling conditions were as followed: 95 °C for 10 s; 40 cycles at 95 °C for 10 s; 52 °C for 10 s; and 72 °C for 30s.

We selected TERT and GAPDH as reference genes after assessing the stability of the gene expression level of eight putative housekeeping genes using BioGazelle qBase + v1.5 software (Hellemans et al., 2007). The primers used were:

DMT (gene of interest) F: 5′-ACATCGGGCATTTGGACCTT-3′; R: 5′-CAGTGCACCAACGAAAGCAA-3′

TERT (housekeeping gene) F: 5′-TTTCGCCCTCGTTTTCCTCA-3′; R: 5′-AGAAATCACGACCACACGCT-3′

GAPDH (housekeeping gene) F: 5′-CGTTGAGACGGAGTGCTTGA-3′; R: 5′-GATGGACTTCTCGGCACTGA-3′

Statistical analyses of laboratory experiments

For the ISM accumulation/extrusion and gene expression level experiments, two-tailed unpaired Mann-Whitney-t test were carried out in GraphPad Prism 5. Data were represented as means ± standard error of the mean. P values ≤ 0.05 were considered to be statistically significant.

Results

Genomic changes in a resistant parasite under drug pressure

To evaluate the impact of a continuous ISM pressure on the T. congolense genome we maintained the drug resistant field isolate TRT15 for 10 passages in mice with 1 mg/kg ISM (TRT15p10-ISM). Although infected mice became microscopically non-parasitemic after each round of ISM treatment, parasitemia in mice reappeared within less than a week.

Bioinformatic analyses revealed a total of 157,797 SNPs and 3098 indels in TRT15 in comparison to the T. congolense IL3000 reference genome when SNP calling was done with the combined sequence data of all TRT15 samples obtained during the experiment. Following drug pressure, we could not detect any novel SNPs nor indels appearing in TRT15. However, we found that the read depth frequency of alleles at heterozygous loci changed over time (Table 1). A shift in read depth allele frequency of at least 0.5 was observed at 11 loci after 5 passages and at 62 loci after 10 passages under drug pressure, while no such shifts were observed when TRT15 was not subjected to drug pressure (Table 1). Of these 62 loci, 17 shifts resulted in read depth allele frequencies larger than 0.8 following drug pressure, but these were mostly found within genes coding for hypothetical proteins (results not shown).

Table 1

Number of allele frequency shifts identified in each experiment. In the first experiment, we maintained the drug resistant field isolate TRT15 for 5 and 10 passages in mice either with 1 mg/kg ISM (i.e. TRT15p5-ISM and TRT15p10-ISM, respectively) or without the drug pressure (i.e. TRT15p5 and TRT15p10, respectively). In the second experiment, three independent ISM resistant lines were generated starting from the same initial MSOROM7 clone, namely: MSOROM7-G, MSOROM7-H and MSOROM7-I. In the last experiment, we compared the resistant and sensitive field isolates of each natural populations identified in (Tihon et al., 2017).

TRT15 compared to	Number of shifts
	TRT15p5	TRT15p10	TRT15p5-ISM	TRT15p10-ISM
Shift ≥0.50	0	0	11	62
Shift ≥0.333	2	2	93	499
Shift ≥0.25	25	27	248	941
MSOROM7 compared to	MSOROM7-G	MSOROM7-H	MSOROM7-I
Shift ≥0.50	163	4	26
Shift ≥0.333	1496	190	627
Shift ≥0.25	1784	972	1955
Resistant compared to Sensitive	TC1	TC2	Zambia
Shift ≥0.50	2	12	3
Shift ≥0.333	10	32	42
Shift ≥0.25	30	70	132

While the continuous ISM pressure did not induce novel SNPs nor indels it caused copy number variations. Compared to the TRT15 genome before drug pressure we identified a total of 36 and 116 gene duplications and 7 and 38 gene deletions in TRT15 after 5 and 10 passages under drug pressure, respectively (Table S2) (Fig. S1 A and B). In contrast, only 19 and 13 gene duplications and 3 and 13 gene deletions were identified after 5 and 10 passages without drug pressure, respectively (Table S2) (Fig. S1 A and B).

Induction of drug resistance in a sensitive parasite

Three independent ISM resistant lines were generated in the murine host (i.e. immune-suppressed) starting from the sensitive clone MSOROM7 by a gradual ISM induction pressure. These three lines were referred to as MSOROM7-G, MSOROM7-H and MSOROM7-I. Full resistance to 1 mg/kg ISM was acquired after three to five months under increasing drug pressure (i.e. five passages in the murine host) (Fig. 2). MSOROM7-H acquired full resistance more rapidly, already after 88 days under ISM challenge, while MSOROM7-G and MSOROM7-I acquired full resistance after 108 and 145 days, respectively. Notably, once ISM resistance reached a threshold value of 0.04 mg/kg (i.e. passage four), these parasites were also fully resistant to the 1 mg/kg ISM dose (i.e. passage five) (Fig. 2). All three cell lines maintained their acquired ISM resistant phenotype when grown in immune-competent mice for six months without any drug pressure.

Fig. 2

Induction of ISM resistance in the immune-suppressed murine host (Endoxantreated). The T. congolense sensitive strain MSOROM7 acquired ISM resistance between three to five months in three independent experiments. In total, five ISM-dose passages in the murine hosts were sufficient to acquire full ISM resistance: passage one (P1): resistant to 0.005 mg/kg ISM; passage two (P2): resistant to 0.01 mg/kg ISM; passage three (P3): resistant to 0.02 mg/kg ISM; passage four (P4): resistant to 0.04 mg/kg ISM; passage five (P5): resistant to 1 mg/kg ISM. The three independently ISM resistant lines that acquired full resistance after passage five (1 mg/kg ISM) are referred to as MSOROM7-G, MSOROM7-H and MSOROM7-I. In immune-competent mice, the acquisition of fully resistant parasites failed: only one parasite strain (MSOROM7-C) acquired a partial resistance to 0.04 mg/kg ISM after a period of >600 days.

In contrast, the same induction set-up in immune-competent mice resulted in only one line (MSOROM7-C) acquiring partial ISM resistance to 0.04 mg/kg after 20 months of drug pressure with the reappearance of parasitemia only observed in 1/3 of the ISM-treated animals. However, this line did not acquire full resistance to 1 mg/kg ISM after more than two years of drug pressure (Fig. 2). All other attempts to acquire resistant parasites in the immune-competent mice failed because of parasite clearance following drug administration.

Genomic differences between sensitive and resistant parasites

Bioinformatic analyses revealed a total of 93,850 SNPs in comparison to the T. congolense IL3000 reference genome when SNP calling was done with a combined sequence dataset of all MSOROM7 samples that were obtained during the experiment. We also identified 1786 small indels in the MSOROM7 lines compared to the reference genome IL3000, of which 86 in coding regions, including 51 frameshift mutations. Following drug resistance induction, we examined genomic changes but did not identify any novel SNPs nor indels appearing in any of the MSOROM7 ISM resistant lines. However, we did find many genomic changes in the read depth frequency of allele at heterozygous loci (Table 1), while similar shifts were not observed for indel sites (results not shown).

A total of 1784; 972 and 1955 loci showed an allele frequency shift larger than 0.25 (z score ≥ 3) in MSOROM7-G, MSOROM7-H and MSOROM7-I respectively (Table 1). Of these, 55, 25 and 23 resulted in a read depth frequency larger than 0.8 in MSOROM7-G, MSOROM7-H and MSOROM7-I lines, respectively (Table S3), of which four were common among the three resistant MSOROM7 lines but none were common with the resistant TRT15p10-ISM line (Fig. 3A and C; Table S3).

Fig. 3

Number of frequency shifts in the alternate alleles identified in: A) the MSOROM7 resistant lines compared with the original MSOROM7 sensitive line. Four homozygous changes in nucleotide sequences were common to all resistant strains and included one synonymous change in the gene TcIL3000.11.2040 encoding a hypothetical protein, and three non-coding changes (see Table S3: Red); B) Between the resistant and sensitive parasites of each T. congolense populations (i.e. TC1, TC2 and Zambia) identified in (Tihon et al., 2017); C) Between the resistant and sensitive field isolates (i.e. combined TC1, TC2 and Zambian strains), between the MSOROM7 resistant lines compared with the sensitive MSOROM7, and between TRT15p10-ISM compared with TRT15. In total, 11 shifts were common between the resistant field isolates and the MSOROM7 resistant lines, but no common shifts were identified between the field isolates, the MSOROM7 resistant lines and TRT15p10-ISM D) The allele frequency shifts identified in these experiments were found in a total of 53 genes in the field-isolated strains, in 30 genes in the MSOROM7 resistant lines, and in 314 genes in the TRT15p10-ISM isolate.

In MSOROM7-G, four out of 55 allele frequency shifts were detected at non-synonymous positions in the ABC transporter genes TcIL3000_7_4990 and TcIL3000_10_9180, while one shift was detected at a non-synonymous position in the gene TcIL3000.11.12160, an orthologue of the T. b. brucei gene Tb11.01.3390 encoding the mitochondrial Type II topoisomerase (Table S3: Bold).

In MSOROM7-H, we discovered that 10 out of 25 shifts were located in the gene TcIL3000_8_1200 that encodes for the permease of the drug/metabolite transporter (DMT). Eight of these shifts were detected at non-synonymous positions and reached a read depth frequency larger than 0.8 at passage three and remained there until the end of the experiment (Table 2, Table S3: Bold). We closely examined all SNPs including SNP clusters (i.e. which were excluded from the original analysis, see (Tihon et al., 2017)) in the chromosomal region spanning the DMT gene. In the sensitive MSOROM7, we observed that the alternate allele frequencies of the genes TcIL3000_8_1190 (block ‘d’ in Fig. 4A and D) and TcIL3000_8_1200 (block ‘e’ in Fig. 4A and D) were of 0.25, 0.5, 0.75 or 1 (Fig. 4B). Read depth analysis of this genomic region revealed the presence of two copies of both genes in the sensitive strain (Fig. 4D), possibly explaining this peculiar tetrasomic-like alternate allele frequency pattern. In addition, we identified that almost all alleles within the DMT gene shifted to higher frequencies in MSOROM7-H compared with MSOROM7, while this was not the case for alleles found upstream of the gene (Fig. 4C). We further observed that in MSOROM7-H the normalized read depth started to decrease from passage two until full resistance, revealing the presence of only one copy of both genes in this resistant line (Fig. 4D). Similar observations were made for MSOROM7-G and MSOROM7-I (Fig. 4D) and for the resistant TRT15 (data not shown) that also had only one copy of the DMT gene. We examined the expression level of the DMT gene and observed a significant downregulation in all induced resistant lines but MSOROM7-G (p > 0.05) compared with MSOROM7 (Fig. 4E). Interestingly, MSOROM7-H, which had more definitive allele frequency shifts at non-synonymous positions, showed the lowest expression of the DMT gene (Fig. 4E).

Table 2

Alternate allele frequencies (including these of SNP clusters) in the gene The alternate allele frequencies ranged from 0.5 (i.e. heterozygous SNP) to 1 (i.e. homozygous SNP). Alternate allele frequencies of 0.8 or higher are shown in bold and illustrate the homozygous SNPs identified in the intermediate and resistant MSOROM7-H parasites compared with the MSOROM7 sensitive line. Ref: Reference allele. NYS: Non-synonymous SNP; SYN: Synonymous SNP. P2: Resistant to 0.01 mg/kg ISM; P3: Resistant to 0.02 mg/kg ISM; P4: Resistant to 0.04 mg/kg ISM; P5: Resistant to 1 mg/kg ISM.

	Impact	position (bp)	Ref.	SNP	Alternate allele frequencies
					MSOROM7	MSOROM7-H
						P2	P3	P4	P5
	NSY	77	G	A	0.5	0.3	1.0	1.0	1.0
SNP cluster	NSY	122	G	A	0.5	0.4	0.8	0.9	1.0
SNP cluster	NSY	124	T	G	0.5	0.4	0.9	0.9	1.0
SNP cluster	NSY	128	G	T	0.5	0.4	0.9	1.0	1.0
	NSY	224	G	A	0.5	0.4	0.8	1.0	1.0
	NSY	401	T	C	0.5	0.6	0.8	0.9	1.0
SNP cluster	NSY	419	A	C	0.5	0.5	0.8	0.9	1.0
SNP cluster	NSY	422	A	G	0.5	0.5	0.8	0.8	1.0
SNP cluster	SYN	423	G	A	0.5	0.5	0.7	0.8	0.9
	SYN	483	A	G	0.5	0.5	0.9	0.9	0.9
	SYN	696	C	T	0.6	0.5	0.9	0.8	0.9
	NSY	697	G	C	0.6	0.5	0.9	0.8	0.9
	NSY	709	G	T	0.6	0.5	0.9	0.9	0.9
	NSY	731	C	G	0.6	0.5	0.9	0.9	0.9
	NSY	733	A	G	0.6	0.5	1.0	0.9	0.9
	NSY	791	C	T	0.6	0.5	0.9	0.9	1.0
SNP cluster	SYN	855	T	C	0.5	0.3	1.0	1.0	1.0
SNP cluster	NSY	856	A	G	0.4	0.3	1.0	1.0	1.0
SNP cluster	SYN	862	C	T	0.5	0.3	1.0	1.0	1.0

Fig. 4

Alternate allele frequencies (including these of SNP clusters) in the Note that the TcIL3000_8_1200 gene is longer in the MSOROM7 lines compared with the reference TcIL3000 genome because of a frameshift indel at the gene position 823. A) Genomic region analyzed. a. TcIL3000_8_1150 b. TcIL3000_8_1160 c. TcIL3000_8_1170 d. TcIL3000_8_1190 e. TcIL3000_8_1200. This figure is based on the number of SNPs found in this genomic region to highlight the distribution of SNPs. The scale does not reflect the physical size of these genes. A gap in the reference genome after the gene TcIL3000_8_1200 prevented further analysis of the downstream region of the gene. B) Alternate allele frequencies in the sensitive MSOROM7; C) Alternate allele frequencies in the resistant MSOROM7-H. Each dot/triangle represents one allele position. Only the frequencies of the alternate alleles located in the TcIL3000_8_1200 gene shifted and became fixed in the resistant population (black box) D) Normalized read depth. In the MSOROM7 sensitive strain, the normalized read depth indicates the presence of two copies of the genes TcIL3000_8_1190 and TcIL3000_8_1200. These gene duplications were confirmed by sequencing the MSOROM7 sensitive strain with TruSeq DNA libraries using the methodology described in (Dumetz et al., 2017) (data not shown). The normalized read depth decreased in the resistant lines, suggesting the loss of one of these copies in response to drug pressure. E) Expression level of the DMT gene TcIL3000_8_1200 (i.e. evaluation of the mRNA expression level) in MSOROM7 sensitive and resistant lines challenged or not with 1 mg/kg ISM by Real-Time qPCR. Average ± standard error of the mean; *: P < 0.05; **: P < 0.01 (n = 6).

In MSOROM7-I, notable allele frequency shifts included three non-synonymous SNPs in the genes TcIL3000_8_410 and TcIL3000_10_13700 encoding transmembrane proteins (Table S3: Bold).

Six of the genes identified in these ISM resistant lines (i.e. TcIL3000_7_4990, TcIL3000_10_9180, TcIL3000.11.12160, TcIL3000_8_1200, TcIL3000_8_410 and TcIL3000_10_13700) also showed shifts larger than 0.25 in T. congolense natural populations (i.e. TC1, TC2 and Zambia; see (Tihon et al., 2017)) when the read depth allele frequencies were compared between ISM sensitive and ISM resistant parasites (Table 3: Red). Shifts in the allele frequency were also found in additional common genes in drug resistant MSOROM7, the natural resistant isolates and/or TRT15 maintained for 10 passages under drug pressure (Table 3).

Table 3

List of genes in which allele frequency shifts were detected in two or more independent experiments (i.e. selection experiment; TRT15 experiment; comparison between resistant and sensitive field isolates). Allele frequency shifts in the gene TcIL3000.11.4410 were identified in all three experiments (bold). Two genes showed allele frequency shifts in the MSOROM7 resistant lines and TRT15p10-ISM isolate (underlined). Allele frequency shifts were detected in 10 common genes between the MSOROM7 resistant lines and the field isolates (italic; red). An additional four genes showed allele frequency shifts in both TRT15p10-ISM and the field isolates.

While we did not identify novel SNPs nor indels in the MSOROM7 resistant lines, we observed new structural variations of the genome. Gene copy number analysis revealed 31 gene deletions in MSOROM7-G4, 17 gene deletions in MSOROM7-H and 39 gene deletions in MSOROM7-I compared with the MSOROM7 sensitive line (Table S4). Overall, 22 gene deletions were shared between two or more resistant lines, which were found in the gene TcIL3000_0_00780 encoding a transposase and in genes coding for unspecified products (Fig. 5A; Table S4: Bold). Interestingly, four of the deletions identified in the TRT15p10-ISM were common to gene deletions detected in one or more MSOROM7 resistant lines (Fig. 5A; Table S4: Red), but these affected genes are coding for unspecified products. We also identified 58 gene deletions in the resistant TC2 field isolates compared to their sensitive counterparts, three of which were also observed in the resistant MSOROM7 lines (i.e. non coding RNA genes; Table S4: Green, Fig. S1 C).

Fig. 5

Structural variations: gene deletions and gene duplications in the resistant MSOROM7 lines and TRT15p10-ISM compared with MSOROM7. A) Number of long deletions detected in the MSOROM7 resistant lines compared with the MSOROM7 sensitive line, and number of long deletions detected in TRT15p10-ISM compared with TRT15. B) Number of genes with at least one additional copy in the MSOROM7 resistant lines compared with the MSOROM7 sensitive line, and in TRT15p10-ISM compared with TRT15. C) Changes in gene copy number of 13 genes coding for VSGs (i.e. 5 genes), Helicases (i.e. 5 genes) and unnamed products (i.e. 3 genes) in the MSOROM7 intermediate and resistant lines compared with the MSOROM7 sensitive line, and in TRT15p10-ISM compared with TRT15, and statistical significance (complete list Table S5).

Furthermore, 272 genes had at least one additional copy in the MSOROM7 resistant lines compared with the sensitive line (Fig. 5B; Table S5). The gene copy number of 79 of these 272 genes also increased in TRT15p10-ISM compared with TRT15 (Fig. 5C for a representative example; Table S5: Bold). Some of the highlights of these 272 genes were 115 genes coding for VSGs or VSG associated proteins; five genes coding for cell surface expressed proteins; six genes coding for ATP-dependent dead h RNA helicases (RNA helicases); seven genes coding for retroposon hot spot proteins and 11 genes coding for the transferrin receptor ESAG6. In addition, the duplication of seven genes was common between TC1 or Zambian resistant field isolates and one or more MSOROM7 resistant lines (Fig. S1 D; Table S5: Underlined).

Reduced ISM accumulation in resistant parasites

A significant decrease in ISM accumulation was observed by flow cytometry in all resistant trypanosomes compared with the sensitive MSOROM7 parental line (Fig. 6A and C). The accumulation of ISM in the sensitive and resistant parasites gradually increased as ISM concentrations augmented, until reaching a plateau at a concentration of ISM of 5 μg/ml in the incubation medium. At this concentration, the sensitive cells displayed a mean ISM fluorescent signal of 599 (standard deviation = 64.6) compared to 327 (standard deviation = 52.6) for the resistant parasites. ISM exposure started to kill parasites at a concentration of 5 μg/ml and higher, as indicated by the increased signal of the Fixable Viability Stain (FVS) (data not shown). No differences in the accumulation of ISM could be detected between the three induced-resistant lines. A similar decrease in ISM accumulation was also observed for the ISM resistant field isolate TRT1 (Fig. 6A).

Fig. 6

Accumulation and uptake of ISM in the sensitive and resistant MSOROM7 lines and in the resistant field isolate TRT1. A) and C) Evaluation of ISM accumulation by flow cytometry (n = 3). A significant decrease in the accumulation of ISM is observed in all resistant strains compared with the sensitive MSOROM7 parental line. B) Measurement of the uptake and extrusion of ISM by an ex vivo fluorescent assay (n = 4, except MSOROM7 where n = 7). The uptake of ISM was significantly reduced in MSOROM7-G, MSOROM7-I and TRT1 (p < 0.05) compared with the sensitive MSOROM7 parental line. After 60 min (dashed line) parasites were placed in a drug free medium for 10 and 20 min to evaluate ISM extrusion. Average ± standard error of the mean; *: P < 0.05, **: P < 0.01.

Then, we assessed the dynamics of the accumulation and extrusion of ISM in the sensitive and resistant parasites and observed a fast accumulation of ISM within the first 10 min of incubation for all cell lines (Fig. 6B). After this initial period, the ISM accumulation remained at the same level over time. The observed ISM accumulation in the resistant MSOROM7-G, MSOROM7-I and TRT1 was significantly lower than for the sensitive MSOROM7 line (Fig. 6B). For both sensitive and resistant parasites, the concentration of ISM accumulating inside trypanosomes exceeded the initial starting concentration of 1 μg/ml, suggesting the involvement of an active transport of ISM. Finally, the extrusion of ISM by the parasites displayed a similar profile for all the strains (Fig. 6B).

Discussion

The drug ISM plays a key role in the treatment of T. congolense infections. However, little information is available about its cellular uptake and accumulation, its mechanisms of action and development of drug resistance in this parasite, impeding a proper management of ISM resistance in the field. Here, we used WGS to study the genetic changes that may underlie the acquisition of ISM resistance within an in vivo model.

We demonstrated the importance of the immune system of the host in preventing the development of ISM resistance. In an immune-competent host, sensitive parasites failed to acquire resistance to 1 mg/kg ISM after more than two years under drug pressure. In immune-suppressed animals, sensitive parasites evolved to fully drug resistant parasites after three to five months of drug pressure. Similar observations were previously described in T. b. evansi parasites (Osman et al., 1992, Mutugi et al., 1994), and it has been hypothesized that trypanocidal drugs, rather than killing parasites, induce restrictions in the replication of the trypanosomes that will then be eliminated by the immune system of the host (Geerts et al., 2001). This mechanism could explain partly why trypanosomes can overcome an ISM pressure in immune-deficient animals. The link between a weakened host immune system and drug resistance acquisition in protozoan parasites is further highlighted by the observation that impaired immunity resulted in treatment failures and drug resistance to visceral leishmaniasis and malaria in HIV co-infected patients (Troya et al., 2008, Mondal et al., 2010, Menendez et al., 2011). In the field, impaired animal immunity could result from injuries, infections with other pathogens such as intestinal worms, malnutrition and stressful situations. Proper livestock health management is therefore highly important to lower the probability of ISM resistance development.

Within the drug resistant parasites, gene deletions and duplications could be detected in contrast to novel SNPs or indels when comparing them with the drug sensitive parasites. Moreover, no significant changes were observed within the genes coding for the F1-subunit γ and other proteins interfering with the mitochondrial F1FO-ATP synthase, suggesting that T. congolense developed different ISM resistance mechanisms than T. b. brucei (i.e. in vitro-induced resistance), T. b. evansi and T. b. equiperdum (Brun and Lun, 1994, Dean et al., 2013, Gould and Schnaufer, 2014, Baker et al., 2015). In addition, the resistant MSOROM7 T. congolense parasites were still fully capable of completing successfully their life cycle in the tsetse fly vector (data not shown), indicating the presence of a functional kDNA. Interestingly, we identified significant shifts in allele frequency among heterozygous loci associated with non-synonymous changes in genes coding for ABC transporters, DMT and transmembrane products. Specifically on the DMT gene, we observed a strong increase in the frequency of 19 alternate alleles at heterozygous sites, resulting in the fixation of 14 non-synonymous changes in nucleotide sequence in MSOROM7-H. This fixation probably resulted from the deletion of one of the two copies of the DMT gene that were present in the original sensitive genome. In order to explain the adjacent alternate allele frequencies of 0.25, 0.5 and 0.75 in the sensitive strain (see Fig. 4B), one must assume that both copies show different genotypes, as explained in the following theoretical example: (i) BB/BA/BB (three sites, B and A standing for alternate and the reference alleles, respectively) and (ii) AA/AA/BA, hereby explaining an alternate allele frequency of 0.5/0.25/0.75 for the respective three sites. The loss of the second gene in the resistant MSOROM7-H would lead to a unique BB/BA/BB, with alternate allele frequencies of 1/0.5/1, resulting in fixation of a series of alternate alleles.

The loss of one DMT gene copy was also observed in MSOROM7-G and MSOROM7-I. However, since we did not observe any fixation of alternate alleles at non-synonymous positions, we concluded that the first copy of the gene was lost in these two resistant lines.

Interestingly, the expression level of the DMT product was significantly downregulated in the resistant parasites, especially in MSOROM7-H where allele frequency shifts at non-synonymous positions were observed. This suggests that the combination of the loss of one copy of the DMT genes and the fixation of non-synonymous mutations had a higher impact on the gene expression level. This may have contributed to the observed decrease in ISM accumulation within the drug resistant parasites.

Members of the DMT superfamily have been involved in antibiotics resistance in Staphylococcus aureus (Grinius and Goldberg, 1994) and Escherichia coli (Yerushalmi et al., 1995), and in chloroquine resistance in P. falciparum (Martin and Kirk, 2004). Similarly, an important role of ABC transporters in drug resistance has previously been proposed in T. congolense (Delespaux et al., 2005) and has been documented in many other organisms including T. brucei, T. cruzi, P. falciparum, Schistosoma mansoni, and Leishmania parasites (Ouellette et al., 2001, Perez-Victoria et al., 2002, Shahi et al., 2002, Alibu et al., 2006, Sanderson et al., 2009, Franco et al., 2015, Garg and Goyal, 2015, Pinto-Almeida et al., 2015, Veiga et al., 2016). Changes in other transport proteins also conferred resistance in Trypanosoma spp. (Barrett et al., 1995, Matovu et al., 2003, de Koning et al., 2004, Witola et al., 2004, Munday et al., 2015). The potential link between the observed allele frequency shifts and ISM transport modifications still has to be confirmed, but the genes listed in this paper represent promising targets for further functional experimental analyses.

Our genomic analyses also revealed that new gene deletions and gene amplifications may possibly play a role in ISM resistance. However, it is difficult to establish a link between ISM resistance and the observed gene deletions/amplifications due to relatively poor annotation of the T. congolense genome. The fact that 79 gene amplifications were common between the experimentally induced MSOROM7 resistant parasites and the TRT15 parasites under drug pressure for 10 passages may suggest that these changes in gene copy numbers reflect the mutagenic properties of ISM. This is possibly the consequence of the high affinity of ISM for double-stranded DNA molecules (Wilkes et al., 1995) that leads to its accumulation in both the mitochondrion and nucleus of the parasites where it induces double strand breaks and subsequently inhibits the kDNA (Shapiro and Englund, 1990) and DNA (Kaminsky et al., 1997) replications. To overcome DNA damages and maintain its integrity, cells possess a DNA repair machinery that consists of a complex of distinct repair pathways (reviewed in (Genois et al., 2014)). This includes a double strand breaks (DSB) repair that involves both homologous and non-homologous recombination (Genois et al., 2014). The potential involvement of the DSB repair in response to an ISM pressure could explain the observed shifts in read depth allele frequency (i.e. homologous recombination where the sister chromatid acts as repair template) and CNVs (i.e. homologous recombination between homologous repeats to amplify or delete DNA sequences) (Genois et al., 2014). This hypothesis is supported by the observation of the deletion of the DMT gene (Fig. 4) and other gene deletions and duplications (Fig. 5), possibly as the result of DSB repair events covering these genes. In addition, allele frequency shifts and CNVs were also observed in the TRT15 isolate under drug pressure but not in the absence of drug pressure, further suggesting that ISM induces strong mutagenic effects.

The potential involvement of the DNA repair machinery in drug resistance mechanisms has been suggested in P. falciparum where the DNA repair machinery was upregulated to respond to exogenous stress (i.e. DNA damaging agents) and promoted cell survival (Gupta et al., 2016), and in Leishmania parasites where the amplification of P-glycoprotein genes through homologous recombination conferred drug resistance (Ouellette et al., 1991, Grondin et al., 1993). Similar mechanisms could potentially exist in T. congolense in response to the DNA damaging drug ISM to promote cell survival and drug resistance.

Other genome modifying mechanisms such as point mutations and sexual recombination could also play a role in some of the genomic modifications observed in the resistant parasites. However, there is currently no evidence about a potential recombination cycle in bloodstream form T. congolense parasites (i.e. meiosis is taking place in the tsetse fly salivary glands in T. brucei (Gibson et al., 2008, Peacock et al., 2011)).

Conclusion

It is clear that the acquisition of ISM resistance in T. congolense is a complex phenomenon and that multiple mechanisms could lead to the acquisition of the drug resistant phenotype known as the ‘many roads to drug resistance’ concept (Hefnawy et al., 2017). While we demonstrated that the immune system of the host plays a key role in hindering the acquisition of ISM resistance, we were able to induce drug resistance in a sensitive isolate in an immune-suppressed host. Subsequent genomic analyses revealed that non-synonymous mutations within genes encoding ABC transporters, DMT and some ‘unidentified’ transmembrane proteins may explain the development of drug resistance. The next step would be to perform functional analyses and confirm the involvement of these genes in ISM resistance.