Genetic variation in metronidazole metabolism and oxidative stress pathways in clinical Giardia lamblia assemblage A and B isolates

Introduction

Giardia lamblia is a microaerophilic, eukaryotic and tetraploid parasite, annually responsible for up to ~280 million human infections worldwide.¹ Giardia is endemic in low- and middle-income countries (LMIC), and infection with the parasite has been shown to affect growth and cognitive functions in children in resource-poor regions.^2,3 In high-income countries, Giardia is usually transmitted by waterborne outbreaks or found in returning travelers from LMIC.^4,5 Giardiasis is clinically associated with prolonged diarrhea and stomach pain, and it may lead to post-infectious sequela, eg, irritable bowel syndrome and chronic fatigue.^6,7

The antibiotic metronidazole (MTZ) is a prodrug belonging to the 5-nitromidazole chemical group, and a common first-line treatment for giardiasis.⁸ With its activity against microaerophilic and anaerobic pathogens, MTZ has a treatment efficacy of 73–100%.⁹ MTZ enters the parasite through passive diffusion, it is activated through enzymatic reduction and releases toxic intermediates. This reaction can only occur in an anaerobic or microaerophilic environment.^10,11 Effects of activated MTZ on susceptible Giardia include protein malfunction, multiplication disruption, DNA damage and oxidative stress which will lead to cell death.^12–15 MTZ can be activated by four main enzymes in Giardia: pyruvate-flavodoxin oxidoreductases (PFOR-1 and 2), the thiol-cycling associated enzyme thioredoxin reductase (TrxR)¹⁶ and nitroreductase (NR)-1. The iron-sulfur (Fe-S) cluster binding proteins, ferredoxins (Fd) have also been shown to participate in activation of MTZ in Giardia.¹⁶ NR-2, a paralog of NR-1, is hypothesized to detoxify MTZ by fully reducing the drug without making toxic intermediates.^17–19 While Giardia is exposed to MTZ, it also has to constantly detoxify oxygen, nitric oxide (NO) and other harmful chemicals that may affect growth and viability.²⁰

The majority of proposed resistance mechanisms consist of either down-regulation of enzymes that render MTZ cytotoxic (NR-1, PFORs and TrxR), or up-regulation of enzymes that can make MTZ an inert molecule (NR-2).¹⁹ Most work regarding drug resistance mechanisms in Giardia has been done using laboratory induced resistant assemblage A1 isolates, which rarely infects humans, compared to infections with Giardia assemblage A2 and B, which are common. MTZ is probably activated through several pathways and exhibit a pleiotropic mode of action, and resistance to MTZ has generally taken longer time to develop compared to resistance against most other antibiotics.²¹

MTZ resistance in Giardia is a growing problem. Clinical experience and recent studies of returning travelers have shown that MTZ-refractory Giardia cases have a prevalence of ~20%.^22,23 The prevalence has dramatically increased over the last few years, reaching up to 40% in a referral outpatient clinic in London.²⁴ The highest prevalence of MTZ refractory cases have been found in returning travelers from Asia, especially India, or the south and east Mediterranean areas.^22,24–26

No standardized susceptibility assays for Giardia isolates are available,²⁷ and susceptibility testing in the laboratory is challenging due to genomic complexity of the parasite. Moreover, isolates are difficult to culture, and growth rates may differ between isolates and assemblages. The microaerophilic environment in which Giardia thrives is difficult to mimic in vitro for drug susceptibility assays. Experiments may be performed in unfavorable high levels of O₂, affecting growth rates and MTZ metabolism of Giardia. Thus, a reliable evaluation of MTZ’s activation and tolerance can be difficult to obtain.²⁸ Studies of MTZ metabolism pathway genes in Giardia are mostly based on gene expression analysis in isogenic laboratory-derived resistant Giardia of the assemblage AI strains, rarely found in humans.^{13,18,29–34} So far, no specific markers for resistance in Giardia have been discovered, and resistance is likely complex and based on several factors. These factors could include amino acid variations in proteins and altered function, post- or pre-translational variants, changes in expression of activating or protective proteins, modifications of regulatory pathways and epigenetic modifications,^16,35–39 and a recent study of resistant Giardia strains showed that MTZ is capable of inducing significant changes in proteins found in the antioxidant, electron transport, and pyruvate catabolism networks. In addition, a correlation between acetylation of these proteins and MTZ resistance was found.³⁵

Begaydarova et al 2015 investigated clinical isolates and linked PFOR gene expression to MTZ resistance, but PFOR was not found to be a good resistance marker.⁴⁰ Another study by Galeh et al 2016⁴¹ linked non-synonymous single nucleotide variants (nsSNVs) in PFOR and NR genes to MTZ resistance and possibly through reduced translation and thereby reduced MTZ activation. Genetic variation in Giardia genome is not well known, except for a few genotyping genes, ie, triose phosphate isomerase (tpi), glutamate dehydrogenase (gdh) and beta-giardin (bg).⁴² The genetic variation among Giardia isolates has in one study been found to be low for some assemblage A genes, ie, Fd3 (DHA2_154390) with more SNPs in assemblage B.⁴³ However, the background genetic variability of the parasites’ MTZ metabolizing pathways and management of the oxidative and nitrosative stress has not been investigated in detail and constitutes an important basis for further research into potential markers of MTZ resistance.

Laboratory-induced MTZ resistance may be lost during one en/excystation cycle,⁴⁴ indicating that metabolic adaptations are likely to contribute to resistance. Still, the recent emergence of rapidly increasing clinical MTZ resistance raises the suspicion that strains with a genetically inheritable trait for resistance are now circulating. Indeed, one laboratory-induced resistance line has been shown to possess a non-sense mutation in one of its NR-1 alleles.¹⁶ However, this challenge may best be studied by examining current clinical isolates rather than laboratory strains obtained decades ago. The prospect for whole genome sequencing of purified cysts from clinical stool samples makes it possible to perform such studies.⁴⁵ The aim of our study was to investigate the general genetic variation in proteins involved in metabolic pathways associated with MTZ resistance or management of oxidative stress substances.

Material and methods

Genes of interest

A PubMed literature search was performed using the word Giardia in combinations with one or two of the words: MTZ, refractory, resistance and oxidative stress, to identify genes in the MTZ metabolism pathway and oxidative and nitrosative stress management in Giardia. The search aimed to identify genes that were upregulated or downregulated during exposure to MTZ and/or free radicals or associated with the activation or inactivation of MTZ and genes with antioxidant properties. These publications and their references were used to establish a list of 29 candidate genes (see Table 1). In order to identify pairwise shared genes and the correct gene ID for Giardia assemblage A2 and B, supplementary data from Adam et al was used.⁴⁶ Finally, the Giardia DB (http://giardiadb.org) was used to find gene orthologs and paralogs. For the NR genes analyzed in this study, the annotations NR-1 and NR-2 have been used for clarity. NR-2 is also known as GlNR1/Fd-NR2 (gene ID: Gl50803_22677), and NR-1 as GlNR2/Fd-NR1 (gene ID: Gl50803_6175). For the two genes PFOR-1 and 2, the literature annotations were contradicting, and no more specifications could be found in the Giardia DB. The annotations used for the PFORs in this study are according to GenBank. The annotations of Fds used in our study is based on Müller et al 2008 and Ansell et al 2017.^16,47 One hypothetical Fd presented by Nixon et al 2002⁴⁸ and Ansell et al 2017,¹⁶ Fd5, could not be found in the reference genomes of A2 and B, nor by GenBank blast. Another hypothetical Fd (DHA2_153401 and GSB_151614) was named Fd6 in tables and figures. Some of the putative Fds presented in this study may actually be domains in a larger protein but are presented as single Fds until more is known.

Table 1 Genes included in the SNV analysis for metronidazole metabolism, oxidative- and nitrosative-stress management. The ferredoxins have been named according to Ansell et al 2017,16 and another Fd, DHA2_153401/GSB_151614, earlier presented by Nixon et al48, has been given the name Fd6 in this table. The other genes of interest have been named according to annotations from Giardia DB

Giardia isolates

Twelve isolates of Giardia assemblage A2 and eight isolates of assemblage B were cultured from human clinical stool samples between 2011 and 2015 at the Robert Koch Institute, Berlin, Germany. The stool samples were collected by the Institute of Tropical Medicine and International Health, Charité, Berlin, Germany. One of the isolates, P324, was obtained from the Tropenmedizinische Ambulanz, University hospital Düsseldorf, Germany. No clinical data for the samples have been collected. Giardia cysts from the stool samples were purified using a two-step sucrose gradient centrifugation. The cysts were excysted using a standard two-step excystation method,⁴⁹ in order to acquire trophozoites. The trophozoites were grown in Diamond’s TY-S-33 medium supplemented with bile at 37°C according to.⁵⁰ Collection of the trophozoites was done from confluent cultures in 11 mL tubes, and then the DNA was extracted using Maxwell® 16 FFPE Plus LEV DNA Purification Kit (Promega BioTech AB, Sweden, cat AS1135) together with Maxwell® 16 instrument (AS1000). The quality and amount of DNA was validated with NanoDrop and the concentration of the DNA was verified with Qubit 2.0 (Life Technologies, Grand Island, NY).

Genome sequencing

Fragmented Giardia DNA (2 μg, 550 bp insert size) from the trophozoites was first prepared and indexed using Illumina TruSeq DNA PCR-Free Sample Prep Kit (Cat. Nr 20015962 and 20015960S). Library size distribution was validated using the Agilent Technologies 2100 BioAnalyzer, High Sensitivity DNA chip, and quantification was done using the KAPA Library Quantification Kit for Illumina sequencing platforms. Normalized whole genome Giardia DNA-libraries were sequenced using the Illumina MiSeq paired-end technology (2×300 bp). FASTQ-files from the sequencing were assembled and aligned to the corresponding reference genomes of Giardia assemblage A2 and B (v.26) using Bowtie 2-2.2.3.⁵¹ The reference genomes were obtained from Giardia DB in January 2016 versions 2013-11-25.⁵²

Variant analysis

Aligned data were analyzed using Geneious v.10.2.4. Variant calling was done to identify SNVs using the parameters minimum variant frequency 0.1 and a minimum coverage of 10 nucleotides. The coding regions (CDS) of 29 candidate genes (including pre, 150 bp, - and post-coding regions, 50 bp) were extracted from Giardia assemblage A2 and B in Geneious. For genes with introns (Fd3 in both assemblages), genetic analysis was done in the CDS and SNVs in introns are therefore not presented. To obtain all the SNVs in one alignment, the 29 genes from each isolate were concatenated, and then compared against each other assemblage-wise. Comparisons of the 29 genes were also performed in Geneious using a global alignment with free end gaps. Alignments of ferredoxin domains are based on findings from Nixon et al 2002,⁴⁸ and for more detailed information see the Supplementary material. For each gene, the total number of possible SNV positions has been used. The SNVs were divided into nsSNVs or synonymous SNVs (sSNVs) in the coding regions or as SNVs in the non-coding regions. Ratios for sSNV and nsSNVs were calculated for each gene.

Homology modeling of the nitroreductases and ferredoxins

Homology modeling of the NRs and ferredoxins 1–6 were done based on the amino acid sequences used in the gene analysis. Models were created using Protein Homology/analogY Recognition Engine V 2.0 (Phyre2) server with intensive mode.⁵³ Models were visualized and superimposed in PyMol (The PyMOL Molecular Graphics System, Version 2.1.0 Schrödinger, LLC).

Storage of data

All genetic data is kept on a secure server; Norwegian e-infrastructure for Life Sciences.

Results

Whole genome sequencing of 20 cultured Giardia assemblage A2 and B isolates was successfully carried out with an average coverage of 31.1 (110.3) (see Table S1) for the 29 genes of interest. Flavohemoprotein (DHA2_154000 and GSB_151570), was found to have a substantially higher coverage (87.6 (54.2)), where a total of 17 of the 20 Giardia assemblage A2 and B isolates had higher than average coverage.

Based on the criteria in the literature search, 29 candidate genes were selected and are presented in Table 1.

For the 29 candidate genes, the genetic variation was analyzed across the 12 assemblage A2 and eight assemblage B clinical isolates. The average number of nsSNVs per gene length was 0.3 % for assemblage A2 and 1.2% for assemblage B. Averages for SNVs for each gene are presented in Table S2. One of the genes in this study, Fd3, was found to have a 35 bp intron in assemblage A2 and a 36 bp intron in assemblage B, which supports earlier findings.^43,54

Genes involved in the MTZ and oxidative stress and nitrosative stress management

Considerably more nsSNVs could be found in the Giardia assemblage B isolates than in the assemblage A2 isolates. A color-graded heatmap of the nsSNVs causing amino acid changes in the 29 proteins involved in MTZ metabolism, oxidative and nitrosative stress management is therefore only shown for Giardia assemblage B (Figure 1). The genetic variation in the flavohemoprotein (GSB_151570) is not illustrated in Figure 1 due to high frequency of SNVs and variably higher coverage than the average.

Figure 1 Heat map of proteins in the metabolism of metronidazole and oxidative/nitrosative stress management in Giardia assemblage B. Metronidazole (MTZ) passively diffuses through the trophozoite membrane and needs to be activated into toxic intermediates to execute its function as an antibiotic. Several theoretical intermediates exist for MTZ. MTZ can be activated either through the enzymes pyruvate: ferredoxin oxidoreductase (PFOR)-1 and/or 2 with ferredoxin (fd) as a co-factor, by nitroreductase-1 with flavin mononucleotide (FMN) as a co-factor or by thioredoxin reductase with the redox cofactor flavin adenine dinucleotide (FAD). MTZ may also be converted into an inert metabolite through the enzyme nitroreductase-2. The MTZ-NO formed during activation may react with O2− and create the reactive free radical molecule peroxynitrite (ONOO−). In order to remove toxic free radicals many enzymes exhibit protective functions, ie, thioredoxin peroxidase enzymes may convert the peroxynitrite to the harmless molecule nitrate (NO3−). The free radical form of NO· may be converted to nitrate by flavohemoprotein (O2 dependent reaction). In the microaerophilic environment, Giardia is repeatedly exposed to the harmful O2/O2−. The O2 may be metabolized to H2O2 through desulfoferredoxin (SOR) and further H2O2 is converted to H2O by the thioredoxin peroxidase enzymes. O2 may also be converted to H2O by NADH oxidase enzymes or the O2 scavenging enzyme flavoprotein. Other enzymes may cause free radicals in Giardia, ie, NADPH oxidoreductase putative is areactive oxygen species (ROS) contributor due to conversion of O2 to a more reactive free radical form (O2−). Another contributor to ROS is the reduced version of thioredoxin and it may initiate protein misfolding. The color gradation of the proteins represents the number of nsSNVs positions per length of each gene. *Thioredoxin reductase is active in both MTZ metabolism and reduction of the protein thioredoxin. **Flavohemoprotein is a suspected multicopy gene and the number of SNVs has not been defined and is not represented by color in the heatmap.

For two of the genes with highest frequency of genetic variability in assemblage B (GSB_151294 - thioredoxin peroxidase and GSB_153178 – NR-1), the nsSNVs constituted over 4% of the nucleotides in the gene, and only two genes were found to be without any nsSNVs (histone H2A; GSB_151412 and Fd6; GSB_151614). The average percentage of any SNV per CDS length was 5.4% in Giardia assemblage B and 1.1% in assemblage A2.

The majority of the genes was shown to have at least one sSNV, except for the four genes Fd3; DHA2_154390, NR-1; DHA2_153380, Fd1; DHA2_9662 and NADPH oxidoreductase putative; DHA2_17151 in assemblage A2. sSNVs were present in all genes in assemblage B. SNVs in the pre-CDS regions were shown to be higher in assemblage B (average of 10.1 SNVs in 150 bp) than in assemblage A2, where the average was found to be 1.1 SNVs per 150 bps. Post-CDS regions had generally lower variability (0.3 in A2 and 4.2 in B for 50 bp).

Nitroreductases

Three different NRs were investigated for the presence of SNVs in this study, the NR-1, NR-2 and the oxygen-insensitive NTR-1. NR-1, had one of the highest numbers of nsSNVs (4.1%) for Giardia assemblage B, with more nsSNVs than sSNVs in the CDS. However, there were lower than average numbers of SNVs in the non-coding regions of the NR-1 in assemblage B. The NR-1 gene in assemblage A2 was found to be fully conserved without any SNVs in the coding or non-coding regions. The NR-2 gene sequence also had high numbers of nsSNVs in Giardia assemblage B (2.1%), and SNVs in the pre-CDS regions were higher than the average (21 vs 10.1 out of 150 bp). The NR-2 in Giardia assemblage A2 had low numbers of SNVs, where only one was a nsSNV. The NTR-1 gene was fully conserved in assemblage A2, and had lower than average numbers on non-syn SNVs in B.

Ferredoxins

The Fds are proposed to be co-factors for the MTZ activating proteins, PFOR-1 and 2, and changes in these co-factors could potentially affect their function, binding and role. We found that three of the ferredoxins in the Giardia assemblage B isolates, Fd2; GSB_10329, Fd4; GSB_153527 and Fd1; GSB_9662, had higher numbers of nsSNVs than the average of 1.2% (respectively 3.1%, 1.7% and 1.7%). Interestingly the Fd2 was one of the proteins with the most SNVs in the post- and pre-CDS regions of all the genes investigated in assemblage B, and it also had a high number of nsSNVs in the coding region.

For Giardia assemblage A2 two of the ferredoxins had higher numbers of nsSNVs than the average of 0.3%, respectively Fd2; DHA2_10329 (2.1%) and Fd6; DHA2_153401 (1.1%). One fascinating finding was that the Fd6 was fully conserved in Giardia assemblage B while assemblage A2 had some SNVs present.

Antioxidant and metabolic proteins

Two of the free radical protective enzymes, prx (GSB_151294) and NADH oxidase (GSB_9719) and the metabolic protein alcohol dehydrogenase (GSB_13350) was demonstrated to have higher numbers of nsSNVs in assemblage B than the average. Post-CDS analysis revealed that TrxR; GSB_9827 was one of the proteins with the highest number of SNVs, while NADPH oxidoreductase putative; GSB_17150, the two prxs, GSB_153801 and GSB_14521, malate dehydrogenase; GSB_3331, flavoprotein; GSB_10358 and alcohol dehydrogenase; GSB_13350, all had higher numbers of post-CDS SNVs than the average.

Alignment of the ferredoxins and the nitroreductases

Some of the ferredoxins and NRs were found to have high numbers of nsSNVs in many of the Giardia assemblage A and B isolates. The amino acid changing SNVs and their positions in the proteins were therefore analyzed further, and an alignment was made to compare the cysteine domains in the five different ferredoxin proteins included in this analysis in addition to NR-1 and NR-2 and a hypothetical Fd previously described in WB (Fd5)⁴⁸ The proteins NR1, Fd1, and Fd3 in assemblage A2 and Fd3 and Fd6 in assemblage B had no nsSNVs in their ferredoxin domains. The nsSNVs positioned around the cysteine residues in the ferredoxin domains responsible for forming the Fe-S clusters binding sites (electron transport sites) are shown in Figure 2.

Figure 2 Alignment of the ferredoxins (Fd1-Fd6) and nitroreductase 1 and 2. The genetic variability is represented by alternative amino acids under the original protein sequences. Ferredoxin 5 was found to be an uncharacterized protein and is not annotated in the DHA2 and GSB genomes and is thus presented by the hypothetical sequence from Giardia WB. Only the ferredoxin domains in the N-terminus of the proteins are represented, and the nitro-flavin mononucleotide (FMN) reductase domain in the C-terminus of the NRs are not shown. The cysteine residues responsible for forming the iron-sulfur clusters are shaded in either black or magenta to represent how the residues connect. *end of the protein, +the protein has been cut and only the ferredoxin domain is shown.

Homology modeling of the ferredoxin domains

In order to be able to analyze the effect of the amino acid changes, homology modeling was used. The ferredoxin domain sequence alignment presented in Figure 2 and the putative cysteine residues for binding Fe-S clusters are coded with the same color codes as shown in the homology models of the domains in Figure 3. Homology modeling indicates various differences within the ferredoxin domains.

Figure 3 Cartoon tube representation of homology models of Nitroreductase (NR) ferredoxin domains and ferredoxins (Fd) 1–6. Cysteine residues responsible for the binding of [Fe4-S4] iron–sulfur cluster are represented with sticks colored with dark gray and red for clusters 1 and 2, respectively. For the mutations caused by SNVs, the original amino acid residue is shown in stick and possible mutations are marked with one letter code. (A) NR1_A2 (grey), NR2_A2 (cyan), NR1_B (yellow) and NR2_B (brown) (B) Fd1_A2 (grey), Fd2_A2 (cyan), Fd1_B (yellow) and Fd2_B (brown) (C) Fd3_A2 (grey), Fd3_B (cyan). The extended C-terminus has been truncated from the figure and is marked with an asterisk. (D) Fd4_A2 (grey), Fd4_B (cyan). (E) Fd5_WB (grey) (F) Fd6_A2 (grey), Fd6_B (cyan).

Fd1, Fd2 and Fd4 can be superimposed on ferredoxin domains of NR-1 and 2. All of them clearly have two cysteine motifs, with four cysteines each, able to bind Fe-S clusters, probably of the type 4Fe-4S.^55,56 The rest of the ferredoxins differ from the others by shape and ability to bind Fe-S clusters. The Fd6 has two cysteine motifs but one of the motifs only has three coordinating cysteines instead of four, which could indicate that the Fd6 only can bind one 4Fe-4S cluster or potentially one 4Fe-4S and an additional 3Fe-4S cluster (bound to 3-cystein motif). The Fd3 has only one clear motif for binding 4Fe-4S cluster and in addition it has a proposed elongated C-terminus. Fd5 has both an elongated N- and a C-terminus.

Some of the amino acid substitutions were positioned close to the ferredoxin domain with Fe-S cluster binding cysteines and sufficiently different in charge or hydrophobicity to be able to potentially cause a change significant enough to affect protein function. Some of the amino acids in the ferredoxin domains are surface exposed and changes here could potentially affect the proteins’ ability to interact with other proteins in possible intermolecular electron transfers. These changes include changes from a positively charged amino acid such as K to a negatively charged E or D seen in NR-1 assemblage B, or a change from E to K in Fd2 assemblage A2 and in Fd2 from assemblage B where G is changed to R.

In NR-1 from assemblage B, S can be changed to P (pos. 70 Figure 2), and in NR-2 from A2, P can change to S (pos. 41 in Figure 2) which might change the conformation of the protein backbone, as P is structurally more rigid than other amino acids. These SNVs are located next to or between two cysteines coordinating Fe-S cluster binding in NR-1 and NR-2, respectively. SNVs that could potentially have an effect to electrostatic properties of Fe-S clusters are G to R (pos 66) change in NR-2 from assemblage B and G to E (pos 38) change in Fd4 from assemblage B.

Discussion

The genome of Giardia trophozoites is tetraploid, consisting of two diploid nuclei. The variability among and within the assemblages can therefore potentially be high. In our study we show that the nsSNVs are common in the investigated genes in culturable clinical isolates of Giardia assemblage A2 and B. Generally, Giardia assemblage B had higher numbers of SNVs in the 29 genes investigated, and relates well to the earlier finding that the allelic sequence heterozygosity (ASH) in assemblage B is up to 0.53%,⁵⁷ whereas the ASH of two A2 isolates, AS-98 and AS-175, has been determined to be 0.25–0.35%.⁴⁶ In Giardia assemblage B, particularly high numbers of nsSNVs were found in gene sequences encoding MTZ activating enzyme NR-1, the PFOR co-factors Fd2, Fd1 and Fd4, the protective enzymes prx and NADH oxidase, the metabolic enzyme alcohol dehydrogenase, in addition to the MTZ inactivating enzyme NR-2. In assemblage A2, gene sequences of Fd2, Fd6 and the ROS-contributing enzyme NADPH oxidoreductase putative had the greatest number of nsSNVs.

Homology modeling shows that some of the amino acid changing SNVs are in close proximity to the cysteine residues responsible for forming Fe-S clusters and could potentially affect binding affinities of the proteins NR-1, NR-2, Fd2 and Fd4, or how they may exert their roles as electron transferring proteins.

Genetic variation as potential resistance markers

High numbers of nsSNVs found in the Giardia assemblage B isolates may indicate that assemblage B mutates faster than assemblage A2 isolates. The higher variation frequencies found in the B isolates may be associated with greater functional variation, higher prevalence of disease in humans and higher rates of treatment-refractory cases.^58,59

A recent study,⁴¹ using clinical Giardia isolates, suggested that nsSNVs in NR-2 and PFOR genes could be linked to resistant Giardia isolates, but the study was based on only two treatment-refractory isolates. In addition, this study gives potential evidence that different haplotypes of these two genes may exist within the genome. This finding correlates with our observation that the numbers of nsSNVs in the NR-2 protein in Giardia assemblage B was higher than average. The PFORs in our study were found to have lower numbers of nsSNVs than average. However, the co-factors for the PFORs, Fds, were found to be highly variable within three of the five Fds in assemblage B and two of the Fds in Giardia assemblage A2.

The role of the MTZ activating enzyme, PFOR, is thought to be the transfer of electrons, through its bound Fe-S clusters, to soluble Fds.³⁸ Lower levels of PFOR and ferredoxin enzyme activities have been characterized in resistant Giardia cultures and were also connected to MTZ sensitivity in clinical isolates.^16,32,60,61 PFOR gene expression was, however, not found to be of significant difference in patient isolates with different susceptibility to MTZ in vitro.⁴⁰

In other anaerobe organisms, genetic markers of resistance have been identified. The enzyme known as oxygen-insensitive nitroreductase, has been shown to be an MTZ activating enzyme in Helicobacter pylori and Escherichia coli, and bacteria with inactivated NTR or nsSNVs in this gene, may be resistant.^62–64 Additionally, resistance in the parasite Trichomonas vaginalis has been associated to SNVs in the two nitroreductase genes ntr4Tv and ntr6Tv.⁶⁵ SNVs found in intergenic regions of 13 genes in T. vaginalis have been linked to resistance as well, and the study suggests that a panel of genetic markers could be used as a diagnostic tool of resistance.⁶⁶ Intergenic analysis should be done in Giardia as well, in order to predict whether SNVs here may play a role in MTZ resistance.

An enzyme in Giardia known as oxygen-insensitive nitroreductase (NTR-1), is similar to the NTR in bacteria, and it might have been acquired through lateral transfer from bacteria.⁴⁸ The NTR-1 is different from the NR-1 and NR-2 due to missing ferredoxin domains in the protein and may not play a key role in the Fe-S cluster electron transportation. We found the presence of nsSNVs in NTR-1 to be low in our analysis. Whether this enzyme has a role in MTZ activation in Giardia is not known, but it may have a fundamental role in protection against oxidative- and nitrosative stress, as it is O₂ insensitive, and has been shown to be upregulated during MTZ exposure.¹⁵

Giardia could also potentially down-regulate enzymes that detoxify oxygen to protect itself from the activated form of MTZ. This down-regulation is favorable due to the fact that MTZ may only be reduced in microaerophilic or anaerobic environments. Presence of oxygen helps to deactivate MTZ by futile cycling, creating ROS.⁶⁰ The ROS will then require extensive detoxification by enzymes in the antioxidant system to avoid harmful cell damage.⁶⁷ In Giardia assemblage B, high variability was found in the two oxygen detoxifying proteins NADH oxidase (GSB_9719) while the ROS contributing enzyme NADPH oxidoreductase putative (GSB_17150) was found to have higher variability than the average for both assemblages. Oxidoreductases have previously been shown to be upregulated when Giardia is challenged with H₂O_2.³⁷ The functions and characteristics of the NADPH oxidoreductases and their potential paralogs should be investigated in future studies, to understand whether they have a role MTZ resistance, as oxidoreductases with null mutations, causing non-functional proteins, have been associated with MTZ resistance in H. pylori.⁶⁸

Some of the genes investigated in the present study were found to have few, or non-existent, nsSNVs (Fd6 and Histone H2A in assemblage B, in addition to 10 genes in assemblage A2, see Table S2). Changes in these conserved genes could potentially have more serious consequences and affect biological functions such as survival in culture, clinical infectivity and tolerance to free radicals, and could be something to address in future studies of culturable and non-culturable isolates.

Resistance induced in the laboratory may have different manifestations than clinical resistance, and the use of Giardia isolates directly from patients, without undergoing laboratory culturing are likely to supplement studies with isogenic laboratory strains.^69,70

Non-coding regions

SNVs found in non-CDS regions of genes could affect regulatory elements and promoters which may further potentially affect gene transcription, mRNA translation and next protein expression.^35,71 The genes, NR-1 and the PFOR co-factor Fd2, were found to have especially high numbers of pre-CDS SNVs in Giardia assemblage B in our study. The sequence preceding the CDS, harbor transcription promoters, and based on scarce knowledge of these in Giardia, they are short, around 50 bp long, initiator-like AT-rich sequences approximately −65 to −29 base pairs upstream the ATG start codon.^57,72 Low conservation of the promoters in Giardia has previously been reported to be common, and may be related to our findings.⁷² How the presence of SNVs in the promoter regions of genes in Giardia may affect transcription of the genes is however not known yet. The interpretation of the SNV findings in these regions are therefore limited, and in this study, we decided to present them, but not to evaluate their potential regulatory roles.

Ferredoxin domains and SNVs

The alignment of the NRs and Fds in this study, only include the ferredoxin domains, responsible for binding the Fe-S clusters (see Figure 2). This region has earlier been associated with drug activation, as the electrons are thought to be transferred from the ferredoxin domain to the active reduction center in order to reduce MTZ.¹⁸ The active reducing center in the NRs has earlier been characterized as a nitro-flavin mononucleotide (FMN) reductase domain.¹⁸

Most of the nsSNVs present in the ferredoxin domain alignment in Figure 2 will probably not alter the function of the proteins as the changes do not alter the basic properties, eg, polarity or charge, of the amino acid residues.

Some other changes are located close to the cysteine part of the 4Fe-4S binding motif and might, therefore, affect the geometry or the electron transfer properties of the Fe-S cluster. Activated MTZ has previously been speculated to be able to react with nucleotides and cysteines (free thiol groups) and form adducts.³⁸ Whether MTZ can form adducts with the proteins investigated in the present study and whether amino acid substitutions close to the cysteine domains could affect formation must be investigated further using structural- and molecular biology studies.

Other proteins that also have distinct Fd domains include the two large proteins PFOR-1 and PFOR-2. They have been hypothesized to have three motifs to bind 4FeS4 cluster. Even if the PFORs have their own Fd domains, they are still thought to use Fds in the cytosol for electron transfer and MTZ activation.³⁸

Coverage and multicopy genes

Several of the genes in this study were found to have paralogs annotated in the genomes of Giardia assemblage A2 and B, and we aimed to include genes present in both of the Giardia assemblages. A total of 12 paralogs were found for the NADPH oxidoreductase, and an alignment was done to compare their similarity (data not shown). The alignment showed they were only 78.3% similar to one another, and probably not multicopy genes.

We were not able to investigate the true number of sSNVs and nsSNVs of the gene sequence of flavohemoprotein due to its variable higher coverage than the other genes, indicating that it is present in a variable number of copies in different isolates (Table S1). The SNVs could possibly be a combination of SNVs of all the alleles present of the flavohemoprotein gene and not true SNVs of one gene.

Two shorter paralogs of flavohemoprotein, DHA2_153759 and DHA2_152971 were identified, and had high identity to the longer flavohemoprotein DHA2_154000 (data not shown). Functions of these two proteins should be investigated in detail in the future to know if they have similar properties as the DHA2_154000, or if they potentially could be wrongly annotated or fragmented versions of DHA2_154000. The shorter flavohemoproteins could have affected the coverage of the Giardia assemblage A2 isolates, as they each would align to the 5ʹ or 3ʹ end of the DHA2_154000 gene. They do not overlap, leaving approx. 500 bp in the middle of the gene. The coverage was nevertheless even, and the two paralogs in assemblage A2 would also not explain the higher coverage of the assemblage B gene GSB_151570. The potential copy-number variation in this gene, coupled with amino acid changing SNVs could be important for understanding Giardia’s capability of handling oxidative stress, as flavohemoprotein is an important NO and O₂ scavenging enzyme.^{16,20,36,37,73–76}

We also noted that histone H2A had 24 times lower coverage than the average in one of the assemblage A2 isolates. Possibly this gene could be missing, or found in another contig in this isolate, or a sequencing error may have occurred.

WGS limitations

The Illumina MiSeq paired-end technology sequencing used in our study gave a low coverage for several of the isolates and their genes, limiting the SNV analysis. For two assemblage B and one assemblage, A2 isolates the average coverage for the genes was only 11.3–14.0 (Table S1). Low coverage decreases the possibility of finding low prevalence SNVs, as the cutoff value was set to >10 in coverage. In addition, short reads of 2×300 bp prevented identification of full-length gene haplotypes in our study. The relatively crude reference genomes used in this analysis can also potentially hamper the SNV analysis, as the Giardia isolates used for obtaining the reference genomes have been grown in culture for decades and still consist of several hundred contigs, instead of the predicted five chromosomes.

Conclusion

Genetic variation in the form of SNVs is common in both of the human-infecting Giardia assemblages A2 and B. Assemblage B showed higher variation than assemblage A2, and presence of amino acid changing SNVs were especially high in genes associated with MTZ metabolism (ferredoxins, nitroreductase-1 and 2) and for some of the genes important in the detoxification of free radicals (thioredoxin peroxidase and NADH oxidase). Some of the discovered SNVs in this study could potentially affect the electron transfer properties of Fe-S clusters in some of the NRs and Fds investigated. Further studies are needed to address their functional properties, but study design can be guided by the presented findings. Analysis of genetic variation in circulating MTZ susceptible and treatment-refractory clinical Giardia isolates will be important to address the recent increase in MTZ resistant Giardia cases. The present study can be regarded as a basis for further studies into how genes in the metabolizing pathways may differ between resistant and susceptible Giardia isolates.