Giardia duodenalis trophozoites were generously provided by Professor Jacqueline Upcroft and Professor Peter Upcroft. Metronidazole (Mtz)-sensitive trophozoites (BRIS/83/HEPU/106, BRIS/87/HEPU/713, and WB1B; referred to hereafter as 106-s, 713-s, and WB-s) were maintained in flat-sided 10 mL tubes (Nunclon delta) filled with 10 mL of complete TYI-S33 medium (Davids and Gillin, 2011) and sub-cultured twice weekly. The standard medium contains 55 mM glucose, which is well in excess of the concentration required to sustain trophozoites throughout a four-day culture passage. As such, with a view to correlating transcriptomic results with future metabolomics results (in which glucose can saturate spectral profiles), the TYI-S33 medium was modified to contain 6 mM glucose (designated “TYI medium”). Mtz-resistant lines (BRIS/83/HEPU/106-2ID₁₀, BRIS/83/HEPU/713-M3 and WB1B-M3; hereafter 106-r, 713-r, and WB-r) were maintained under the same conditions with Mtz (Sigma Aldrich; 100 mM stock dissolved in dimethyl sulfoxide) added to a final concentration of 30 μM. References describing the original axenization of these lines and resistance selection methods are provided in Table 1.

Trophozoites from confluent flasks were chilled on wet ice (5 min), pelleted (680 × g, 5 min, 4°C), diluted in fresh TYI medium to 5 × 10⁴ cells/mL (susceptible lines) or 1.5 × 10⁵ cells/mL (resistant lines; cf Tejman-Yarden et al., 2011), and added in 40 μL volumes to wells of a 96-well clear-bottom plate (Corning #3610). 200 μL of sterile water was added to the peripheral wells to limit evaporation. Serial dilutions of Mtz were created in TYI medium and added to wells containing trophozoites, in duplicate. The wells of the resulting plates contained ~2,000 (susceptible) or 6,000 (resistant) trophozoites in a 50 μL volume, and either Mtz (2.5–200 μM) or 1% DMSO in TYI medium (negative control). Plates were incubated at 37°C within GasPak EZ anaerobe pouches, with fresh anaerobic sachets for 48 h. Reconstituted CellTiter-Glo reagent (Promega) was stored at −20°C and thawed to room temperature (22–24°C; RT) 30 min prior to use; 50 μL of reagent were added to wells containing trophozoites, and plates were incubated for 15 min, shaking, at RT. ATP-based luminescence (corresponding to live cells/well) was measured using a luminometer (BioTek), and transformed to a proportion of the negative control value. Prism software (GraphPad) was used to fit Hill plots to dose-response data and to calculate IC₅₀ values, via the “log(inhibitor) vs. normalized response—variable slope” module. In order to estimate generation times, standard curves relating ATP luminescence and trophozoite numbers were created for Mtz-susceptible and -resistant lines, and found to be log-linear within 5 × 10⁴–10⁶ cells per well. Trophozoite numbers in the negative control wells of dose-response experiments were interpolated from the standard curve, and generation times were calculated according to Boreham et al. (1984).

For Mtz-susceptible lines, samples for mRNA sequencing were generated in triplicate. Falcon t25 flasks were seeded with 10⁵ trophozoites in 56 mL of fresh TYI medium and incubated at 37°C for ~60 h, at which time trophozoites formed a confluent monolayer on the flask walls. For Mtz-resistant lines, up to 8 flasks were seeded with 1.5–3 × 10⁶ trophozoites (according to growth-rate), and Mtz was added to a final concentration of 30 μM. After incubation, flasks were inverted, and supernatant and suspended trophozoites (likely to contain dead/dying and dividing cells, which could confound the mRNA integrity and comparative transcription analyses, respectively), were decanted to waste; 50 mL of ice-cold complete TYI medium were added and flasks were incubated on ice-water for 10 min to detach trophozoites. The suspended cells were transferred to 50 mL falcon tubes and pelleted (680 × g, 5 mins, 4°C). Supernatant was removed, and pellets were re-suspended, transferred to UV-irradiated 1.5 mL microfuge tubes, and pelleted (720 × g, 2 mins, RT). Samples from Mtz-resistant lines were consolidated into three microfuge tubes. Supernatant was removed and pellets were dissolved in 1 mL of TriPure reagent (Roche), incubated at RT for 5 min and stored at −80°C.

RNA was extracted according to the TriPure manufacturer's protocol within 4 weeks of sample preparation. The dried RNA pellet was re-suspended in reverse-osmosis deionized water (H₂O) and treated with Turbo DNAse (Ambion) according to the manufacturer's instructions. RNA concentration was estimated by fluorometry (Qubit), and quality control was performed using a BioAnalyzer (Agilent). Polyadenylated RNA was purified from 10 μg of total RNA using Sera-mag oligo(dT) beads, fragmented to a length of 100–500 bases, and reverse transcribed using random hexamers. Strands were labeled using the dUTP second-strand synthesis method, end-repaired and adaptor-ligated, and then treated with uracil-specific excision reagent (USER; NEB) before amplification by PCR. Products were purified over a MinElute column (Qiagen) and single-end sequencing was performed using an Illumina HiSeq 2500 (YourGene Biosciences, Taiwan). Raw RNA sequencing reads are available through the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/SRP075868).

Raw reads were trimmed using Trimmomatic (Lohse et al., 2012) (sliding window: 4 nt; minimum average PHRED quality: 20; leading and trailing: 3 nt; minimum read length: 40 nt), and mapped as single-ended reads to the accepted G. duodenalis gene models (assemblage A genome release 25; giardiaDB.org; Morrison et al., 2007) using RSEM (Li and Dewey, 2011) with the -forward-prob 0 flag. The locus encoding superoxide reductase (Testa et al., 2011) was added to the reference genome. Reads that did not map under these conditions were likely to include antisense transcripts. To quantify antisense transcription, unmapped reads were submitted to RSEM with the –forward-prob 1 flag, which discards any reads that map with <100% confidence to the coding strand. The RSeQC bam_stats.py module (Wang et al., 2012) was used to calculate read mapping statistics; and infer_experiment.py was used to confirm read mapping orientation. Feature detection was calculated as a function of read mapping depth using the Qualimap counts module (v1.0; García-Alcalde et al., 2012) with the −k 10 flag, denoting a minimum mapped read threshold of 10. Saturation plots were displayed in Excel (Microsoft) (Supplementary Figure 1).

Expected counts from RSEM were submitted to edgeR (Robinson et al., 2009). For preliminary clustering analyses, low abundance genes (fewer than one count per million counts; CPM) were discarded before library re-scaling and TMM normalization. For each susceptible and resistant pair, sample normalization was performed as above and differentially transcribed genes were determined at a false discovery rate of 0.01. CPM-values were multiplied by 1,000 and divided by the reference transcript length to give transcripts per kilobase per million transcripts (TPM), suitable for comparing the transcriptional abundance of different genes across samples. Romer (Wu et al., 2010) was used to test enrichment of KEGG pathways and Gene Ontology (GO) terms among up- and down-regulated genes in resistant lines, and to test enrichment of curated gene sets specific to G. duodenalis (significance threshold: corrected p < 0.05; Supplementary Table 9). Hypothetical G. duodenalis proteins of interest were selected for tertiary structure prediction using a stand-alone implementation of the I-TASSER software suite (v3.0) generously provided by Professor Yang Zhang, as described previously (Ansell et al., 2016). It must be noted that the output from the stand-alone implementation of this software can differ from that available through the online server (zhanglab.ccmb.med.umich.edu/I-TASSER) due to computational resource constraints on the server. As such, the output from the stand-alone software should be more robust. GO terms associated with high-confidence structures (TM score > 0.5; GO confidence score > 0.5) were included in gene set testing. Structural similarity searching was performed using TM-align (Zhang and Skolnick, 2005), and protein structures were visualized using UCSF Chimera (Pettersen et al., 2004). Read coverage data were normalized by total library size (calculated in edgeR), and transcriptional abundance plots and heat maps were generated in R (ggplot2 and gplots packages).

RNA reads mapping in the sense orientation were submitted to mpileup and bcftools (SAMtools suite; Li et al., 2009) to identify single nucleotide polymorphisms (SNPs) using conservative thresholds (≥100 × read coverage for the alternative allele, constituting ≥20% of nucleotides mapped at that position). Selected SNPs were visually confirmed using Integrative Genomics Viewer (Robinson et al., 2011), and amino acid substitutions (products of non-synonymous SNPs), were predicted using getORF (Rice et al., 2000) and custom scripts. For comparative transcriptomics between species, published RNA transcript count data for Mtz-susceptible and –resistant T. vaginalis (Bradic et al., 2016) were submitted to edgeR in a two-sample comparison of transcript abundance. Count data for resistant and susceptible G. duodenalis lines generated in this study were compared using the same design, and orthologous genes (defined in OrthoMCL-DB; Chen, 2006) that were differentially transcribed in the same direction in both species, were retained for further analysis.

All Mtz-susceptible lines exhibited generation times of ~5 h. The generation time for WB-r was 13.85 h, indicating a marked growth defect compared to 106-r and 713-r, for which generation times were 6.36 and 7.46 h, respectively. For susceptible lines, Mtz IC₅₀ concentrations ranged from 7.8 to 9.4 μM. Each resistant line showed significantly greater Mtz tolerance than its susceptible parent (two-sample comparison of fits; p < 0.001; Figure 1A). Specifically, WB-r and 106-r exhibited an approximate 2.6-fold increase in Mtz tolerance relative to their susceptible parents, whereas 713-r was 5.4-fold more tolerant than 713-s, and also more tolerant than the other resistant lines (713-r vs. WB-r: p = 0.027; vs. 106-r: p = 0.004; Table 1).

RNA sequencing of six cell lines in biological triplicate produced a total of 641 million single-ended reads, with an average of 35.6 million reads per sample. The vast majority of reads satisfied quality control (mean = 93.5%; range = 92.4–94.7%), of which 80% were mapped to the accepted G. duodenalis genes in the sense orientation (Supplementary Tables 1, 2; Supplementary Figure 1). Around 5% of reads from each sample mapped to the reverse complement, consistent with antisense transcripts (detailed below). After filtering lowly transcribed genes, 5,792 non-deprecated coding domains were detected in all isotypes. Transcripts from a further 231 coding domains were detected in either one or two isotypes (Supplementary Table 3). Using transcript abundance, susceptible lines clustered together in both principal coordinate analysis and Euclidean clustering; although Pearson correlations were lower on average between susceptible lines compared with resistant lines (Pearson's r = 0.93 vs. 0.96). Omission of transcripts encoding VSPs improved the mean pair-wise correlation for susceptible lines (r = 0.98), with little effect on resistant lines (r = 0.97; Supplementary Figure 2).

Pearson correlations based on log₂(fold-change) values in resistant lines, were strongest between 713-r and WB-r (r = 0.47). These lines also correlated with microarray-based fold-change data for 41 strongly differentially regulated genes in the previously published Mtz-resistant line WB-C4 (0.27 < r > 0.37; Müller et al., 2008). Transcriptional changes in 106-r diverged substantially from 713-r and WB-r (r ≤ 0.2), and correlated negatively with WB-C4 (r = −0.29). Comparisons with transcriptomic data for WB-s trophozoites after exposure to sub-lethal stressors (Ansell et al., 2016) revealed strongest agreement between Mtz stress and 106-r (r = 0.18), and a modest correlation with 713-r (r = 0.03). At a false-discovery threshold of 0.01, 2,010 genes (i.e., non-deprecated coding domains) were up-regulated in at least one resistant line relative to the susceptible parent, and 2,120 genes were down-regulated (Figure 1B; Supplementary Tables 4, 5). In 713-r, 1,178 genes were up-regulated, and 1,169 were down-regulated relative to 713-s. In WB-r, 1,175, and 1,176 genes were up- and down-regulated respectively, of which 917 were differentially transcribed in the same direction in 713-r. Comparatively few DTGs (921 in total) were identified for 106-r, supported by close clustering between the resistant and susceptible lines for this isotype (Supplementary Figure 2B; overlapping DTGs are displayed in Figure 1C). This result may also be affected by greater variance between replicates in 106 (mean biological co-efficient of variance = 0.18) compared to WB and 713 isotypes (mean BCV = 0.1). As such, the minimum change in transcript abundance for genes satisfying the statistical threshold for differential transcription in 106-r was 24%, compared to 17% in WB-r and 713-r. The 30 most strongly up- and down-regulated genes in each resistant line (excluding VSPs) are provided in Supplementary Table 6.

Of 156 DTGs common to all three resistant lines, 103 genes were up-regulated, and 53 were down-regulated (Figure 1C). Common up-regulated genes included 33 VSPs, a high-cysteine membrane protein, the mitotic regulator Mad2 (Vicente and Cande, 2014), zinc and calcium transporters, Nek and DYRK kinases, a sulfatase-like enzyme, and a putative fatty acid α-oxidase. A ParcC-like chaperone (GL_15455), and two putative mTor kinase subunits were also up-regulated (GL_6542 and GL_13469). The 53 common down-regulated genes included nitroreductase-1, putative FMN-dependent chromate and quinone oxidoreductases (GL_9719 and GL_17151), a putative calcium sensor, glutamate dehydrogenase, threonine dehydratase, acetyl co-A acetyltransferase, two α-giardins, and E1 and E2 ubiquitin-related enzymes (Figure 2A; Figure 3; Supplementary Table 5). Proton- and ATP-dependent transporters, a hexose transporter and a putative sugar-nucleotide transporter (GL_9036) were also down-regulated. In all resistant lines, threonine dehydratase and a small hypothetical protein (GL_104062) encoded different non-synonymous SNPs. A further 19 genes encoded identical substitutions in two of three resistant lines (Table 2; Supplementary Table 7; Supplementary Figure 3).

According to previous studies, resistant lines may preferentially express certain VSPs (Chen et al., 1995; Müller et al., 2007a, 2008). As mentioned above, disparities in VSP transcription preferentially decreased the strength of correlation between susceptible lines. When the 20 most abundant VSP-coding transcripts (approximating the top decile) in each line were compared, there was greater convergence in resistant lines. Ten genes among the top decile in at least one resistant line (including GL_137617, which was in the top decile for all three resistant lines), were absent from the equivalent group in all susceptible lines. A further nine genes were among the top decile in all resistant lines, whereas only six genes were consistently highly transcribed in susceptible lines (Supplementary Figure 4). No VSPs were consistently suppressed in resistant- relative to susceptible lines. The VSP TSA417 (GL_113797), which was reported as transcriptionally suppressed in the Mtz-resistant lines WB-C4 and WB-C5 (Müller et al., 2007a, 2008), was among the top decile of VSPs in 713-s, and in both 106-s and 106-r. VSPAS7 (GL_137740), which was suppressed in WB-C4 (Müller et al., 2008), was in the top decile for WB-r only. CRP136 (GL_103142), was previously detected in WB-r but not WB-s (Chen et al., 1995). In this study however it was among the top VSP decile in both lines. Overall, these results show little agreement in VSP transcription between lines, and within lines over time. To test the possibility that particular VSP peptide sequences are under selection in the presence of Mtz, we investigated the proportion of clade I, II, and III VSPs among the top decile in resistant and susceptible lines (Adam et al., 2010). Clade II was over-represented in both resistant and susceptible lines, relative to the background prevalence, with little between-group difference (odds ratio = 1.44 and 1.65, respectively). Similarly, there was little difference in the number of CXXC motifs (average of 18–21 motifs), or cysteine content (9–10% cysteine) in the top decile of vsps between lines (Supplementary Table 8). High cysteine membrane proteins are another prominent Giardia-specific protein family with variable numbers of cysteine motifs (Davids et al., 2006). Transcripts from 60 hcmps were significantly enriched among up-regulated genes in 713-r, but not 106-r, or WB-r. A clade of six hcmps whose transcription tended to increase in all resistant lines contained the only consistently up-regulated HCMP, GL_11309. Sixty-nine CXXC motifs are encoded in this gene, well above the median for the HCMP family (Davids et al., 2006). No HCMPs were consistently down-regulated (Supplementary Figure 5).

We undertook further structural and transcriptional investigation of differentially transcribed FMN-dependent NAD(P)H oxidoreductases in resistant lines. GL_9719 and GL_17151 were down-regulated in all resistant lines. The predicted structure of GL_9719 is highly similar to the crystal structures of FMN-binding chromate reductases from Thermus scotoductus (Opperman et al., 2010) and Bacillis subtilis (Kitzing et al., 2005) (Supplementary Figure 7). The putative quinone reductase GL_17151 has two paralogs in assemblage A G. duodenalis: DT-diaphorase, which is up-regulated in 713-r and WB-r; and GL_17150, which is down-regulated in the same lines. The predicted structures of all three proteins are homologous with the KefF (K⁺ efflux regulator) subunit of the Escherichia coli KefF/KefC complex (Figure 5A; Roosild et al., 2009). In E. coli, KefF associates with a NAD-binding KefC subunit via residues that are not conserved in the predicted G. duodenalis homologs. Nevertheless, we identified two hypothetical proteins (GL_8468 and GL_10682) with predicted NADP-binding pockets and structural similarity to KefC. Further, GL_8468 was down-regulated in 713-r and WB-r, in similarity to GL_17150 (Figure 5A). Interestingly, GL_17150 and GL_17151 are adjacent on chromosome 5, suggesting a gene duplication. Whereas these neighbours share only 43% amino acid sequence identity, GL_17150 is 86% identical to DT-diaphorase, which is encoded at a distant locus. Further, GL_17150 is vastly more highly transcribed than either GL_17151 or DT-diaphorase (~800 vs. 40 TPM). To probe a possible mechanism of transcriptional regulation for these paralogous genes, we quantified antisense transcription from GL_17150 and GL_17151. Transcripts antisense to GL_17151 were 4-fold more abundant than those antisense to GL_17150, however the presence of a potentially confounding pseudogene antisense to GL_17150 (GL_34649), should caution the interpretation of these results (Figure 5B).

We further analysed antisense transcription from 2,268 non-overlapping genes that were detected in all isotypes. A total of 792 antisense transcripts were up-regulated across resistant lines, of which 53 were identified in all lines. Alternatively, 703 down-regulated antisense transcripts were identified across resistant lines (34 in all lines; Supplementary Tables 13, 14). A single antisense transcript encoding a solute transporter (GL_9036) was up-regulated in all lines and corresponded with down-regulation of the coding transcript in the same manner. When the fold change in sense and antisense transcripts from non-overlapping genes in resistant lines was compared, positive correlations were identified for 106-r (r = 0.104) and 713-r (r = 0.155; Figure 5C). This finding corresponds to higher proportions of directionally correlated sense and antisense DTGs in all lines (Supplementary Table 13).

Two-sample comparisons between Mtz-resistant and -susceptible G. duodenalis and T. vaginalis revealed orthologous heat shock protein-coding genes, peroxidases and dullard-like phosphatases that were exclusively up-regulated in resistant lines in both species (Supplementary Table 15; Supplementary Figure 8). Two large orthologous groups comprised genes encoding NEK kinase/ankyrin repeat proteins, and protein transport-related molecules, that were up- and down-regulated in both species. Down-regulated orthologous groups contained threonine dehydratases, 40S and 60S ribosomal proteins, the putative G. duodenalis chromate reductase GL_9719 and a T. vaginalis putative etylmaleimide reductase (TVAG_411220), UTP glucose-phosphate metabolizing enzymes; and the previously mentioned putative quinone reductase (GL_17151) and its ortholog TVAG_311580.

The earliest reports of Mtz-resistant G. duodenalis identified changed intracellular redox conditions, and perturbed enzyme activities (Gillin and Reiner, 1982; Smith et al., 1988; Ellis et al., 1993). Recurring themes have since emerged, in particular, the association of ferredoxin family proteins with resistance. PFORs operate at the interface of glycolytic, amino acid, and fatty acid metabolism in G. duodenalis, and supply electrons to ferredoxins that activate Mtz. Both PFOR paralogs require a thiamine pyrophosphate (TPP) catalytic co-factor, and can decarboxylate pyruvate and oxaloacetate. PFOR-2 can also catabolize α-ketobutyrate (Townson et al., 1996). Transcriptional suppression of genes that act proximal to PFORs in all lines, may impose functionally analogous limits on PFOR activity, and thus Mtz activation. Firstly, PFOR activity may be impaired by lower transcription of serine and threonine dehydratases, which produce pyruvate and α-ketobutyrate, respectively. Transcriptional suppression of thiamine pyrophosphokinase, identified in 713-r, may similarly inhibit PFOR activity by restricting supply of the TPP cofactor. That threonine dehydratase acquired non-synonymous mutations in all lines, and was consistently suppressed, may reflect the lower energetic value of its product (proprionyl-coA), which is a less preferred substrate of the ATP-forming acetyl-coA synthetase (Sanchez and Müller, 1996). Other transcriptional changes suggest that accumulation of pyruvate, which has antioxidant effects in G. duodenalis (Biagini et al., 2001b) may be prioritized in resistant lines. Indeed, the suppression of acetyl-coA acetyltransferase in all resistant lines, lower transcription of histone- and N-terminal acetyltransferases in 713-r, and induction of histone deacetylases in 106-r, are all consistent with acetyl-coA conservation (albeit at the cost of lipid metabolism, protein quality control, and epigenetic regulation). This in turn may promote pyruvate accumulation by inhibiting carbon flux through PFOR. In this context, lower glutamate dehydrogenase transcription in all resistant lines might further bolster the pyruvate pool by limiting its consumption through the “glutamate shunt” (Park et al., 1998; Müller, 2003). Glutamate dehydrogenase suppression would also conserve NADPH, a key source of reducing power for antioxidant enzymes such as TrxR (Figure 3).

Down-regulation of the ferredoxin-nitroreductase chimera, nitroreductase-1, was observed in all three lines investigated here, and is reported in the majority of Mtz-resistant G. duodenalis investigated to date (reviewed in Ansell et al., 2015b). Recombinant nitroreductase-1 reduces Mtz when NADH is supplied as an electron donor (Müller et al., 2015). This NADH-dependent mechanism is reminiscent of canonical nitroreductases that reduce Mtz at the NADH-binding site (Koder and Miller, 1998; Martinez-Julvez et al., 2012), and may therefore be catalytically independent of the ferredoxin domain. We contend that the toxicity of nitroreductase-1 in Mtz-exposed G. duodenalis could be due to an additional mechanism of drug activation that involves electron transport through the ferredoxin domain, via the FMN cofactor, to the NADH-binding site. Indeed, although G. duodenalis is widely claimed to lack a ferredoxin:NADH reductase (FNR), such a function has been theoretically attributed to nitroreductase-1 (Ali and Nozaki, 2007; Figure 3). FNR activity in nitroreductase-1 would represent a parsimonious link between pyruvate catabolism and NAD reduction under drug-free conditions (Ellis et al., 1993). By extension, in resistant lines, functional inhibition of PFOR would directly limit nitroreductase activity. The ability of ferredoxin-nitroreductases to accept electrons from PFOR or reduced ferredoxins (as is likely for Fe-hydrogenase; Mulder et al., 2011) and, thereby, to reduce NAD or Mtz, remains to be tested.

Two other putative FMN-dependent oxidoreductases were down-regulated in all resistant lines, indicating potential roles in Mtz activation. The predicted structure of the putative chromate reductase (GL_9719) is similar to that of enzymes with low substrate specificity that reduce xenobiotic electrophilic compounds such as trinitrotolulene (Blehert et al., 1999). The physiological substrates of these enzymes remain elusive, although a chromate reductase in the human gut commensal Bacilis subtilis, is induced after exposure to hydrogen peroxide, and proposed to reduce the unsaturated products of peroxidized lipids (oxylipins; Kitzing et al., 2005). In this context, the up-regulation of a putative fatty-acid α-oxidase in all resistant lines is particularly interesting. In plants, tissue damage or chemical insults induce oxidation of the haem co-factor bound in fatty acid α-oxidase (Zhu et al., 2013), which drives the conversion of fatty acids into oxylipins. These modified lipids are potent signaling molecules (Blée, 2002) that drive transcriptional changes (Eckardt, 2008) and are implicated in regulating growth, stage conversion and pathogenicity in numerous parasitic protists (reviewed in Noverr et al., 2003). Taken together, it may be that the down-regulation of chromate reductase in resistant lines both limits Mtz activation, and prolongs oxylipin signaling, which could, in turn, promote growth in resistant lines. To test this hypothesis, oxylipins could be added to drug-susceptible trophozoites prior to Mtz exposure, or resistant lines treated with lipoxygenase inhibitors.

A KefF-like quinone reductase encoded by GL_17151, was also down-regulated in all resistant lines investigated here, and was dynamically suppressed in a previous study after exposure of drug-sensitive trophozoites to Mtz (Ansell et al., 2016). These results strongly support a role for GL_17151 in Mtz activation. Putative structural homologs of NADP-binding KefC subunits were also found in our study. In E. coli the Kef system couples cytosolic oxidation with pH regulation, which can minimize oxidative damage to biomolecules (Healy et al., 2014). GL_17151 is a paralog of DT-diaphorase, an enzyme that converts oxygen to hydrogen peroxide (Sánchez et al., 2001; Li and Wang, 2006), which is up-regulated in lines 713-r and WB-r. Interestingly, an orthologous enzyme in T. vaginalis with H₂O₂-forming function, termed flavin reductase-1, is down-regulated in Mtz-resistant lines. Loss of this enzyme in T. vaginalis is proposed to confer Mtz resistance by promoting futile redox cycling (Leitsch et al., 2014). The most parsimonious explanation for the contrasting regulation of these orthologs G. duodenalis and T. vaginalis is that DT-diaphorase can detoxify Mtz, whereas the flavin reductase-1 cannot. Indeed, subtle biochemical differences in paralogous enzymes that markedly affect substrate specificity and reaction stoichiometry are well supported; in the ferredoxin nitroreductases of G. duodenalis (Müller et al., 2013, 2015), the canonical nitroreductases in H. pylori (Olekhnovich et al., 2009; Martinez-Julvez et al., 2012), and the expanded flavin reductase family in T. vaginalis (Leitsch et al., 2014). Effective management of functionally divergent paralogs likely requires concerted transcriptional regulation. Further to this, we found an inverse abundance of sense and antisense transcripts arising from the paralogous chromosomal neighbours GL_17150 and GL_17151, consistent with antisense interference. The role of antisense transcripts could be further probed via episomal overexpression of the pseudogene antisense to GL_17150 (Figure 5B), which is replicated antisense to DT-diaphorase. Transcriptional down-regulation of genes orthologous to GL_9719 and GL_17151 in Mtz-resistant T. vaginalis, and the established link between Mtz-resistant H. pylori and mutations in nitroreductase-coding genes (Sisson et al., 2002; Olekhnovich et al., 2009; Martinez-Julvez et al., 2012), demonstrate the centrality of FMN-dependent enzymes to Mtz metabolism and resistance across pathogenic species. Detailed investigation of the localization, interacting partners, substrate specificity, and reaction kinetics of these enzymes is now clearly warranted.

Expression of certain VSPs, the major surface antigens in G. duodenalis, is associated with the development of drug resistance (Chen et al., 1995; Müller et al., 2007a, 2008). These cysteine-rich molecules are proposed to protect G. duodenalis from environmental stresses such as stomach peptidases, and are constitutively turned over, possibly to evade the host immune system (Nash, 2002). The prominent transcription of 10 VSPs in all resistant lines, suggests positive selection from a more varied VSP population in susceptible trophozoites (Supplementary Figure 5). Whereas these dominant genes seem likely to confer a selective advantage in the presence of Mtz, in our investigations their functional significance has remained elusive. Specifically, we found little difference in the cysteine content and motif number between highly-transcribed VSPs in resistant and susceptible lines, and cysteine substitutions were not highly prevalent. It is possible that resistance-associated VSPs play a role that predominates antioxidant activity, for example, in maximizing the stability and integrity of degraded plasma membranes (Touz et al., 2005). In contrast, the high proportion of CxxC motifs in the consistently up-regulated HCMP GL_11309, supports an important function in protecting membranes from oxidative damage in the presence of Mtz. Functional considerations aside, this molecule and the highly transcribed VSPs in resistant lines, could yet be explored in clinical isolates as potential biomarkers of resistance.

The transcriptional changes discussed thus far are common to all three Mtz-resistant lines, and are therefore most likely to account for increased Mtz tolerance. However, it is clear in this study and others, that cellular phenotypes such as growth rate and infectivity diverge between resistant lines (Tejman-Yarden et al., 2011). Those that remain infectious may be more clinically relevant, and thus concerning, than non-infectious resistant lines. The divergent transcriptomes of 106-r and 713-r belie similar growth rates, although only 106-r can establish experimental infections (Tejman-Yarden et al., 2011). Alternatively, despite the transcriptional similarities between 713-r and WB-r, the latter line grows more slowly. Investigating the transcriptional features that are unique to each line can shed some light on the molecular bases of these different cellular phenotypes, and on transcriptional changes that may augment drug resistance.

106-r remains infective, grows relatively well and exhibits the fewest DTGs, which together suggest a highly efficient resistance phenotype. The non-sense mutation in 29% of transcripts encoding nitroreductase-1, could functionally ablate a substantial proportion of the nitroreductase-1 enzyme pool (assuming for simplicity equimolar transcript and protein concentrations), and may thus contribute to the increased Mtz tolerance of 106-r. This finding constitutes the first report of a functional mutation in a known Mtz-activating gene in G. duodenalis, which, combined with transcriptional suppression of the same enzyme, should substantially decrease nitroreductase-1 function in 106-r. In addition, 106-r is the only line in which the NAD scavenging pathway and the putative Mtz-detoxification enzyme, nitroreductase-2, are up-regulated. If these changes preference the flow of glycolytic electrons to nitroreductase-2, (either indirectly via NADH, or directly via the ferredoxin electron transport chain as proposed above), then 106-r may essentially use Mtz as a terminal electron acceptor in glycolysis. Such an adaptation would allow sustained PFOR activity and possibly ATP generation via acetyl-CoA, whilst negating the detrimental effects of the activated drug. Continued ATP generation might in turn decrease reliance on other sources of ATP, such as the arginine dihydrolase pathway, which is transcriptionally suppressed.

Increased NAD scavenging in 106-r might also support heightened activity of NAD-dependent Sir2 deacetylases, which are redox-responsive histone modification enzymes (Fulco et al., 2003). By contrast, transcriptional suppression of a Sir2-coding gene (GL_16569) was found in WB-r, and is previously reported in WB-C4 (Müller et al., 2008), a line whose transcriptional profile correlates negatively with 106-r but positively with 713-r and WB-r. It would be interesting to test whether chemical or genetic perturbation of Sir2 enzymes could recapitulate aspects of these divergent molecular resistance phenotypes. Further, 106-r is the only line in which a set of V-type ATPases and a MATE-like transporter are up-regulated. Like other protists, G. duodenalis may utilize V-type ATPases at the plasma membrane to regulate cytosolic pH (Biagini et al., 2001a). Increased transcription of dynein and kinesin-related proteins in 106-r is consistent with the recycling of such ATPases via microtubule-dependent transport of ATPase-rich vesicles, evidenced in higher eukaryotes (Breton and Brown, 2013). Active proton efflux, coupled to a sodium:proton exchanger (up-regulated in 106-r and WB-r; Biagini et al., 2001a), may decrease the cytosolic sodium concentration sufficiently to drive Mtz efflux via the MATE transporter, which is a sodium-coupled multidrug efflux pump (Du et al., 2015). Indeed, although MATE transporters are generally associated with extruding cationic compounds, the possibility of Mtz efflux through these channels is supported by experiments in H. pylori demonstrating increased Mtz sensitivity MATE-knockout lines (van Amsterdam et al., 2005; Figure 4).

The transcriptome of 713-r is similar to WB-r, despite a marked difference in growth rate. This may derive from active and passive resistance mechanisms that appear unique to 713-r. In addition to lower transcription of nitroreductase-1, a physiological role for NTR-1 in Mtz activation is supported by lower transcription of the encoding gene in this line. Similarly, marked transcriptional suppression of the putative nitroreductase GL_8377, and GL_17150, provides further evidence that several hitherto understudied FMN-binding enzymes may contribute to Mtz activation. Down-regulation of nitroreductase-2 in 713-r may, however, preclude the energetic efficiencies postulated to occur in 106-r. Induction of iron hydrogenase might compensate for lower nitroreductase-2 abundance, and sustain some PFOR activity by directing the flow of electrons to form hydrogen (Lloyd et al., 2002; Figure 3). Further, the marked transcriptional induction of flavohaemoglobin, might represent an alternative Mtz detoxification enzyme. Flavohaemoglobin homologs in bacteria are shown to bind to imidazole drugs (Helmick et al., 2005); and a homolog of this enzyme in the intestinal parasitic protist Blastocystis hominis, is postulated to contribute to Mtz resistance (Wu et al., 2014)—findings that warrant further testing of the ability of the G. duodenalis flavohaemoglobin to detoxify Mtz.

Despite evidence supporting passive Mtz evasion and active detoxification mechanisms in 713-r, other transcriptional features suggest that this line is under greater oxidative stress than 106-r. Up-regulation of components of the thioredoxin system including TrxR, protein disulfide isomerases-2 and -4 (previously reported in WB-C4; Müller et al., 2008), and peroxiredoxin-1ai, indicate reliance on an antioxidant system that collaterally activates Mtz. It is conceivable that Mtz-induced damage necessitates greater thioredoxin activity, creating a positive feedback loop that constitutes an indirect and metabolically burdensome Mtz detoxification system. Although we did not observe increased transcription of NADPH-forming enzymes that might assist thioredoxin reductase activity in 713-r, the functional significance of the T339A substitution predicted in glucose-6-phosphate dehydrogenase, which catalyzes the first step in the NADPH-forming pentose-phosphate pathway, requires further investigation. Lastly in 713-r, up-regulation of protein chaperones suggests that the thioredoxin system alone is insufficient to manage Mtz-induced damage. Induction of peptidyl-prolyl isomerases, HSPs-90 and TCP-1 chaperonins, may help to maintain the structure and function of damaged proteins, although suppression of proteasomal subunits does not support concomitant recycling of highly degraded proteins. Increased activity of ATP-dependent chaperones may be serviced by induction of glycolytic genes in 713-r, especially the ATP-forming phosphoglycerate kinase. Induction of orthologous peroxidredoxins and ATP-dependent HSPs70 and 90 in Mtz-resistant G. duodenalis and T. vaginalis, may constitute universal active resistance mechanisms that complement passive down-regulation of FMN-dependent enzymes.

The transcriptional changes unique to, and uniquely lacking in WB-r may account for its slower growth. Lower transcription of iron-binding enzymes including ferredoxin-3, both PFOR paralogs, putative lipoxygenases and a Fe-S cluster scaffold protein, is intriguing given that a key oxygen detoxification enzyme, NADH oxidase (Brown et al., 1996), is down-regulated solely in WB-r. We previously reported transcriptomic results suggesting heightened oxidative stress in WB-s trophozoites at log and stationary phases of in vitro culture, when PFOR transcription is comparatively low (Ansell et al., 2015a). As the function of Fe-S-binding enzymes such as PFORs is sensitive to dissolved oxygen (Dan et al., 2000), together these findings suggest that WB-r may endure increased intracellular oxygen, consistent with reliance on futile redox cycling. In this context, concerted induction of genes encoding multiple tRNA biosynthesis enzymes, ribosomal proteins and translation initiation factors, may reflect heightened demand for protein synthesis in WB-r, possibly to replace oxidized proteins. Further, as outlined above, lipoxygenases initiate oxylipin signaling pathways in eukaryotes, and oxylipins, such as prostaglandins are implicated in sustaining the growth and development of other protists (Noverr et al., 2003). If the heightened oxidative stress postulated within WB-r is not compatible with sustained expression of iron-dependent lipoxygenases, then a scarcity of oxylipins may compound the observed growth defect. This phenotype may also relate to insufficient induction of important genes in WB-r, relative to the other resistant lines. For example, inorganic sulfate is liberated from sulfate esters by sulfatases, which require sulfatase maturating enzymes to become catalytically active (Bojarová and Williams, 2008). Although a sulfatase-like enzyme is up-regulated in all lines, the Fe-S-dependent maturating enzyme is only up-regulated in 106-r and 713-r. If sulfur metabolism is less efficient in WB-r as a result, then non-protein antioxidant activity could also be impaired (Brown et al., 1998). More broadly, transcripts that were exclusively detected in 106 and 713 isotypes, such as a highly transcribed putative tRNA methyltransferase which may modulate translation efficiency; or those which were only detected the WB isotype, may also contribute to the phenotypic differences between resistant lines. Future work to detail the genetic differences between these isotypes, and between the resistant and susceptible lines therein, will be necessary to differentiate those transcriptional differences based in genetic architecture, from those based purely in dynamic transcriptional regulation. Although the infectivity of WB-r has not been tested, we expect that this line would lack infectivity due to both the similarity between its transcriptional phenotype and that of the non-infective 713-r, and its retarded growth.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.