Phenotypic and transcriptional profiling in Entamoeba histolytica reveal costs to fitness and adaptive responses associated with metronidazole resistance

Antimicrobial chemotherapy is critical in the fight against infectious diseases caused by Entamoeba histolytica. Among the drugs available for the treatment of amebiasis, metronidazole (MTZ) is considered the drug of choice. Recently, in vitro studies have described MTZ resistance and the potential mechanisms involved. Costs to fitness and adaptive responses associated with resistance, however, have not been investigated. In this study we generated an HM-1 derived strain resistant to 12 μM MTZ (MTZR). We examined its phenotypic and transcriptional profile to determine the consequences and mRNA level changes associated with MTZ resistance. Our results indicated increased cell size and granularity, and decreased rates in cell division, adhesion, phagocytosis, cytopathogenicity, and glucose consumption. Transcriptome analysis revealed 142 differentially expressed genes in MTZR. In contrast to other MTZ resistant parasites, MTZR did not down-regulate pyruvate:ferredoxin oxidoreductase, but showed increased expression of genes for a hypothetical protein (HP1) and several iron-sulfur flavoproteins, and downregulation of genes for leucine-rich proteins. Fisher's exact test showed 24 significantly enriched GO terms in MTZR, and a 3-way comparison of modulated genes in MTZR against those of MTZR cultured without MTZ and HM-1 cultured with MTZ, showed that 88 genes were specific to MTZR. Overall, our findings suggested that MTZ resistance is associated with specific transcriptional changes and decreased parasite virulence.


Introduction
Currently, 500 million cases of amebiasis are reported each year (World Health Organization Amoebiasis, 1997). Most of these cases are children from developing nations, who are particularly vulnerable to infectious diseases (Marie and Petri, 2014). Among the drugs available for the treatment of amebic dysentery and amebic liver abscess, metronidazole (MTZ) is one of the most widely used, and is often considered the drug of choice (Löfmark et al., 2010). Treatment with MTZ is usually very effective, i.e., no clinical isolates with high levels of resistance have been observed.
Unlike aerobic cells, anaerobic organisms are sensitive to MTZ because their electron transport proteins have sufficient negative redox potential that can activate MTZ (Samuelson, 1999). After entering the cell, the drug is reduced to a cytotoxic nitro radical anion by electron donors like thioredoxin reductase and ferredoxin (Fdx), which in turn serve as electron acceptors for pyruvate:ferredoxin oxidoreductase (PFOR) (Samuelson, 1999;Leitsch et al., 2007). The nitro radical anion is reduced further to a nitrosoimidazole, and further reduction leads to the formation of hydroxylamine and eventually to an amine (West et al., 1982;Moreno and Docampo, 1985;Ludlum et al., 1988). In E. histolytica, the actual mechanism of action is not fully understood, but evidences from other organisms suggest inhibition of DNA synthesis, and damage to DNA, protein, and other cell components by oxidation and adduct formation (Ludlum et al., 1988;Sisson et al., 2000;Leitsch et al., 2007).
A few mechanisms have been proposed to explain how pathogens survive MTZ treatment. The process varies among organisms, but the general principle includes at least one of the following: altered reduction efficiency (Leiros et al., 2004), drug inactivation (Ralph and Clarke, 1978), reduced drug uptake (Lacey et al., 1993), active efflux (Pumbwe et al., 2007), and increased DNA damage repair (Land and Johnson, 1999). In E. histolytica, Samarawickrema et al. (1997) reported the increased expression of iron-containing superoxide dismutase (Fe-SOD) in a strain resistant to 10 µM MTZ, while Wassmann et al. (1999) showed that in addition to Fe-SOD upregulation, resistance to 40 µM MTZ was associated with increased peroxiredoxin (Prx) expression and decreased expression of Fdx. On the contrary, Tazreiter et al. (2008) showed no substantial increase in Fe-SOD expression and less than 3-fold upregulation of Fdx and Prx in cells exposed to 50 µM MTZ. In addition, all three studies did not find any significant changes in PFOR expression, indicating that its downregulation may not be mandatory for low to modest drug resistance in this parasite. Studies on MTZ resistance in T. vaginalis and G. lamblia also showed variable results regarding the role of PFOR (Quon et al., 1992;Rasoloson et al., 2002;Müller et al., 2008;Leitsch et al., 2011). Strains of T. vaginalis deficient in PFOR activity, for example, showed only low levels of anaerobic resistance to MTZ (Rasoloson et al., 2002). In addition, some clinical strains exhibited resistance only under aerobic conditions, and were completely sensitive to MTZ in the absence of oxygen (Rasoloson et al., 2002). These strains also showed no decrease in PFOR activity (Rasoloson et al., 2002). Other reports, however, have shown decreased Fdx levels in MTZ resistant trichomonads (Ralph et al., 1974;Yarlett et al., 1986) and anaerobic resistance that was correlated with decreased PFOR activity (Kulda et al., 1989). From these results, it is conceivable that different pathways are involved in drug activation, and that various mechanisms for MTZ resistance exist in protozoa.
Currently, most of the work done in E. histolytica focused either on pathways known to activate, i.e., chemically reduce, MTZ (Leitsch et al., 2007) or antioxidants that counteract oxidative stress, such as pyruvate:ferredoxin oxidoreductase, superoxide dismutase, and peroxiredoxin (Samarawickrema et al., 1997;Wassmann et al., 1999). Fitness costs and adaptive responses have not yet been examined. In this study we generated an E. histolytica HM-1-derived isogenic strain resistant to 12 µM MTZ (MTZR). By comparing its phenotype and transcriptional profile against the parental strain, our goal was to identify phenotypic and transcriptional changes related to or involved in MTZ resistance. We compared the transcriptome of a single MTZR line cultured with MTZ [MTZR (+)] against MTZR cultured without the drug [MTZR (−)], and HM-1 exposed to MTZ [HM-1 (+)] to determine the significance of the genes we identified. Finally, to verify their involvement in resistance, some genes were overexpressed by episomal transfection, followed by drug challenge.

Cultivation and Induction of MTZ Resistance
E. histolytica strain HM-1:IMSS cl6 (HM-1) (Diamond et al., 1972) was cultured under axenic conditions in 6 mL BI-S-33 medium (BIS) at 35.5 • C in 13 × 100 mm Pyrex screw cap culture tubes (Corning, Corning, NY) (Diamond et al., 1978). Induction of MTZ resistance was initiated when cells in mid-logarithmic phase of growth were exposed to 1 µM of the drug until late logarithmic phase. In this study, cultures with cell concentrations between 4 × 10 4 and 2 × 10 5 /mL were considered to be in mid-logarithmic phase, and corresponds to cell densities after 1 and 2 days, respectively, after the cultures were initiated with at 2-5 ×10 4 trophozoites per mL. The drug concentration was increased by 1 µM until it reached 4 µM. Cells were maintained with the same dose for 4 weeks before further increase in drug concentration was made. A strain growing at 12 µM MTZ [MTZR (+)] was obtained approximately after 28 weeks, and was used in all experiments. MTZR was also cultured without the drug [MTZR (−)] for 24 h up to 1 week when we tested for reversibility of resistance and when samples were prepared for microarray.

Growth Kinetics
After trophozoite cultures in mid-log phase of HM-1 and MTZR were placed on ice for 5 min to detach cells from the glass surface, cells were collected by centrifugation 500 × g for 5 min at room temperature. After centrifugation, spent medium was discarded and the pellet was resuspended in 1-2 mL of BIS medium. The cell number was estimated on a haemocytometer. Approximately 6×10 4 trophozoites were inoculated in 6 mL BIS medium with or without 12 µM MTZ, and the cultures were examined every 24 h for 5 days. Trypan blue exclusion assay was used to determine the number of viable cells in duplicate cultures (Penuliar et al., 2011).

Half Maximal Inhibitory Concentration (IC 50 ) and Cross-Resistance
Trophozoites in mid-log phase of HM-1 and MTZR were harvested, resuspended in BIS medium as described above, and the concentration was adjusted to 1×10 5 cells/mL. About 1×10 4 cells in 100 µL BIS were seeded per well of a 96-well microtiter plate (Iwaki, Tokyo, Japan) and incubated at 35.5 • C for 1 h under anaerobic conditions using Anaerocult A (Merck, Darmstadt, Germany). One hundred microliters of BIS containing 2-fold increases in concentration of the drugs listed in Table 1 was added, and the plates incubated for 24 and 48 h. The medium was removed and replaced with 200 µL 10% WST-1 reagent (Roche, Indianapolis, IN, USA) in phosphate-buffered saline (PBS) and the plates incubated for 20-40 min. Optical density at 450 nm was measured using a DTX 880 Multimode Detector (Beckman Coulter, Fullerton, CA, USA). The percentage of viable cells was calculated after subtraction of background absorbance as % viability = (absorbance of treated cells/absorbance of untreated cells) × 100%. IC 50 values were calculated using the sigmoidal dose-response equation in GraphPad Prism 5.0 (GraphPad Software, La Jolla, CA, USA). All experiments were repeated three times with two replicates per experiment. To assay the hydrogen peroxide sensitivity of the transformants harboring one of the following plasmids: HP1-HA, ISF1-HA, ISF2-HA, ISF4-HA, and pEhEx-HA, about 1 × 10 4 cells were seeded per well of a 96-well microtiter plate and incubated at 35.5 • C for 12-16 h under anaerobic conditions. The cells were then exposed to 25-1600 µM of hydrogen peroxide for 4 h. Then the cell viability was measured by using WST-1 reagent as described above.

CHO Monolayer Destruction Assay
Rate of CHO monolayer destruction was measured as described previously with minor modifications (Penuliar et al., 2011). Briefly, CHO cells in logarithmic phase of growth were prepared as follows. After spent medium was discarded, CHO cells were incubated with trypsin-EDTA at 37 • C for 5 min to detach cells. Cells were transferred to a 50 mL tube containing 5 mL of Ham's F-12 medium and cell density was estimated on a haemocytometer. Cells were collected by centrifugation at 500 × g for 10 min at room temperature. After the supernatant was discarded, the cell pellet was re-suspended in Ham's F-12 medium at 1 × 10 5 CHO cells per mL. Approximately 100 µL of the CHO suspension was dispensed into each well on a 96-well plate and the plate was incubated at 37 • C as described (Penuliar et al., 2011) for 24-48 h until CHO cells form a monolayer. The medium was removed and the plates washed twice with warm PBS. Approximately 1 × 10 4 cells of HM-1 and MTZR were resuspended in 200 µL Opti-MEM medium (GIBCO, Invitrogen Co., Auckland, New Zealand) and added to each well. The plates were incubated under anaerobic conditions at 35.5 • C for 15 min intervals up to 2 h. The plates were placed on ice for 15 min to release adhered trophozoites and washed twice with cold PBS. The percentage of CHO monolayer destroyed was determined relative to wells without the ameba using 10% WST-1.

Substrate Gel Electrophoresis
Proteinase activity was detected by substrate gel electrophoresis as described previously (Hellberg et al., 2000).

Phagocytosis Assay
Trophozoites of HM-1 and MTZR in mid-logarithmic phase of growth were harvested, resuspended in BIS medium as described above, and the concentration was adjusted to 1 × 10 6 cells/mL. Approximately 1.0 × 10 6 of HM-1 and MTZR trophozoites were seeded on each well of 12-well plates. Spent medium was removed and replaced with 500 µL warmed BIS with 1.5 × 10 7 FluoSpheres carboxylate-modified (2.0 µm) Nile Red fluorescent beads (Invitrogen, Eugene, OR, USA). The plates were sealed and incubated at 35.5 • C for up to 80 min. After incubation, 500 µL of PBS containing 2% galactose was added and the plates incubated on ice for at least 20 min. Formaldehyde was added to a final concentration of 4% and plates incubated on ice for 1 h. The cells were collected, washed three times with cold PBS and analyzed by flow cytometry as previously described (Nakada-Tsukui et al., 2005). Cells that ingested the beads were gated based on increased fluorescence compared to control non-phagocytic cells.

Glucose Consumption
Trophozoites of HM-1 and MTZR in mid-logarithmic phase of growth were harvested, resuspended in BIS medium as described above, and the concentration was adjusted to 1 × 10 7 cells/mL. Approximately 1 × 10 7 cells were seeded in 6-well dishes and incubated for 12 and 24 h. Spent medium was decanted and centrifuged at 500 × g for 5 min. Glucose consumption was determined from the clarified medium using a Glucose (GO) Assay kit (Sigma Aldrich, Saint Louis, MO, USA) according to manufacturer's protocol.

Estimation of Total Protein Amount in a Single Cell
Trophozoites of HM-1 and MTZR in mid-logarithmic phase of growth were harvested as described above. Cell pellet was resuspended in PBS at 1 × 10 6 /mL. Approximately 500 µL of the cell suspension were dispensed to a 1.5 mL tube and the cells were collected by centrifugation at 500 × g for 3 min at 4 • C. After the supernatant was carefully removed, the cell pellet was lysed with 50 µL of lysis buffer [50 mM Tris-HCl, pH 7.5, 150 mM NaCl, 1% Triton-X 100, complete mini EDTA-free protease inhibitor cocktail (Roche Molecular Biochemicals, Mannheim, Germany) and 200 mM E-64]. Protein amount of the lysate were quantified by DC protein assay kit (Bio-Rad Laboratories, Inc., Hercules, CA, USA).

Microarray Hybridization
Samples from two independent RNA extractions were processed using Ambion MessageAmpTM Premier RNA Amplification Kit (Applied Biosystems, Foster City, CA, USA) according to manufacturer's instructions. Fragmented cRNA was then hybridized onto a probe array chip (Eh_Eia520620F) that was custom-made by Affymetrix (Santa Clara, CA, USA) (Husain et al., 2011;Penuliar et al., 2012) using kits and protocols specified in the Affymetrix GeneChip Expression Analysis Technical Manual (Affymetrix, Santa Clara, CA, USA). Following hybridization, arrays were washed, stained with streptavidinphycoerythrin (Molecular Probes, Eugene, OR, USA) using Affymetrix GeneChip Fluidics Station 450, and scanned with an Affymetrix GeneChip Scanner 3000 at 570 nm. Each array image was visually screened to check for scratches, signal artifacts, and debris.
The microarray used in this study was made based on E. histolytica and E. invadens sequences stored at TIGR and Pathema databases (Loftus et al., 2005; http://amoebadb.org/ amoeba/). It contained 9230 probe sets for E. histolytica and an additional 25 and 81 control probe sets for Entamoeba and Affymetrix, respectively (Eh_Eia520620F_Eh). Nomenclature for the IDs was based on whether the probe set is unique to either Pathema (e.g., EHI_123456) or TIGR (e.g., 12.m00345), or is found in both databases (e.g., 98.m00765_234567). Probe sets labeled with "_at" represent a single gene, while those labeled with "_s_at" represent probe sets that share all probes identically with at least two sequences. The "_s_at" probe sets represent highly similar transcripts, shorter forms of alternatively polyadenylated transcripts, or common regions in the 3 ′ ends of multiple alternative splice forms. Probe sets labeled with "_x_at, " on the other hand, represent probe sets where it was not possible to select either a unique probe set or a probe set with identical probes among multiple transcripts. These probe sets could cross-hybridize with other genes in an unpredictable manner (Affymetrix, Santa Clara, CA, USA).
The microarray data reported in this paper has been deposited in the Gene Expression Omnibus (GEO) database (http://www. ncbi.nlm.nih.gov/geo/) with accession number GSE35990.

Generation of Transgenic Ameba Overexpressing HP1 and ISFs
Full-length gene sequences were amplified from MTZR cDNA with sense-and anti-sense primers containing appropriate restriction sites (Table S2). PCR primers were designed as follows: sense primers contained 26∼30 nucleotides corresponding to a region of each gene encoding the amino-terminal portion of each target proteins, fused at the 5 ′ end with three additional nucleotides and SmaI restriction enzyme recognition site (tccCCCGGG); reverse primers contained 22∼33 nucleotides corresponding to a region of each gene encoding the carboxylterminal portion of each target proteins, fused at the 5 ′ end with three additional nucleotides and XhoI restriction enzyme recognition site (ccgCTCGAG). PCR was performed with the cDNA as a template using a DNA Engine Peltier Thermal Cycler (Bio-Rad, Hercules, CA, USA). The PCR cycling conditions consisted of an initial step of denaturation at 94 • C for 1 min, followed by 30 cycles of denaturation at 98 • C for 10 s, annealing at 55 • C for 30 s, and extension at 72 • C for 60 s with Phusion DNA polymerase (New England BioLabs, Tokyo, Japan). Resultant PCR fragments were precipitated with ethanol, digested with appropriate restriction enzymes, and then purified from an exised piece of agarose gel by using GENECLEAN kit (Funakoshi, Tokyo, Japan). After purification, the gene was inserted into an expression plasmid, pEhEx, as previously described (Penuliar et al., 2012). The plasmid was introduced into E. histolytica G3 strain by lipofection, with minor modifications as previously described (Penuliar et al., 2012). Briefly, E. histolytica G3 trophozoites in the mid-logarithmic growth phase were harvested as described above, and re-suspended with Opti-MEM medium supplemented with 5 mg/mL L-cysteine and 1 mg/mL ascorbic acid at 1 × 10 5 cells/mL. Approximately 5 mL of the suspension was dispensed into each well on a 12-well plate and the plate was incubated under anaerobic condition at 35.5 • C for 30 min. Following incubation, 4.5 mL medium from each well was removed and 500 µL liposome-plasmid-mixture (5 µg plasmid, 10 µL PLUS reagent, 20 µL Lipofectamine in Opti-MEM medium) was added. After 5 h of transfection, cells were harvested by placing the plate on ice for 15 min, then added to culture tubes with 5.5 mL cold BIS, and incubated at 35.5 • C for 24 h. Transformants were initially selected in the presence of 1 µg/mL of G418, and drug concentration was gradually increased to 10 µg/mL during the following 6 weeks before the transformants were analyzed.

Reverse Transcriptase (RT)-PCR and Quantitative Real-Time (qRT)-PCR
Total RNA from semi-confluent cultures, about 2 × 10 5 cells/mL, of HM-1, MTZR, and transformants harboring one of the following plasmids: pEhEx-HA (mock plasmid), HP1-HA, ISF1-HA, ISF2-HA, and ISF4-HA, cultured in 25 cm 2 flasks, was extracted and quantified as described previously (Penuliar et al., 2012). cDNA was reverse transcribed from 5 µg RNA with the SuperScript III First-Strand Synthesis System, and RT-PCR was performed with the resulting cDNA as template and primers listed in Table S2 on a DNA Engine Peltier Thermal Cycle (Bio-Rad Laboratories, Inc., Hercules, CA, USA). The PCR parameters used were: an initial denaturation step at 95 • C for 2 min followed by 30 cycles of denaturation at 94 • C for 15 s, annealing at 55 • C for 30 s, and extension at 72 • C for 30 or 60 s. The products were resolved on 1.5% [w/v] agarose gel with 0.5 µg/mL ethidium bromide, visualized, and photographed under ultraviolet illumination.
The Fast SYBR Green Master Mix (Applied Biosystems, Foster City, CA, USA) was used for qRT-PCR in accordance with the manufacturer's instructions. The list of primers for genes whose expression was quantified by qRT-PCR can be found in Table  S1, and includes a housekeeping gene, RNA polymerase II gene (rpoII), as control. Each PCR reaction contained 5 µL (1:50 dilution) of cDNA and 15 µL primer mix, composed of 10 µL 2X Fast SYBR Green Master Mix, sense and antisense primers, and nuclease-free water, to bring the volume to 20 µL. qRT-PCR was performed using StepOne Plus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) with the following cycling conditions: enzyme activation at 95 • C for 20 s, followed by 40 cycles of denaturation at 95 • C for 3 s and annealing/extension at 60 • C for 30 s. All test samples were run in triplicate including an RT-negative control for each sample set along with a blank control consisting of nuclease-free water in place of cDNA. Quantification for each target gene was determined by the Ct method with rpoII as reference gene (Livak and Schmittgen, 2001).

Statistical Analysis
Correlation coefficients were calculated using the Student's ttest function of Microsoft Excel statistical package (Microsoft Corp., Redmond, WA, USA). Probability levels (P) < 0.05 were considered significant.

Generation and Phenotypic Characterization of MTZR
Exposure of HM-1 to MTZ concentrations up to 3 µM for 24 h did not result in any significant changes in cell morphology. Cells continued to grow and cultures formed confluent lawns after 72 h. Starting at 6 µM MTZ, however, cells became rounder than usual and were slightly bigger in size. Cells also took longer to attach to surfaces and for cultures to form confluent lawns. Cell viability, however, was verified by cell adhesion and trypan blue exclusion assay. After 29 weeks, cells growing at 12 µM MTZ (MTZR) were obtained. We tried to increase the concentration of MTZ further, but this resulted in further increase in doubling time.
The population doubling time of MTZR was significantly longer compared to HM-1 (26.8 ± 0.9 and 18.2 ± 0.5 h, respectively; P-value < 0.001, Figure 1). MTZR also appeared rounder, bigger, and more granular than HM-1. In contrast, parental cells exposed to 12 µM MTZ showed signs of morphological changes as early as 2 h, rounding in 5 h, cell shrinkage, and detachment within 24 h. The slight enlargement of MTZR (+) compared to HM-1 was suggested when we measured the total protein content of MTZR and HM-1, 4.81 ± 0.04 and 4.24 ± 0.03 pg/cell, respectively. Flow cytometric analysis also confirmed that MTZR was slightly larger than HM-1 (data not shown).
The IC 50 of MTZ in MTZR and HM-1 was 12.9 ± 0.3 and 6.5 ± 0.3 µM, respectively, and its resistance index was 2.0 ( Table 1). When MTZR was cultured without the drug for 1 week followed by MTZ challenge for 24 h, percent survival was only 28%, indicating that resistance to the drug was reversible (data not shown). To determine if MTZR was cross-resistant with other amebicides, the IC 50 of the drugs listed in Table 1, against MTZR and HM-1 was determined. Results showed that MTZR had reduced sensitivity to the drugs.
MTZR had a slower rate of CHO monolayer destruction, especially at early time points, indicating a slight decrease in virulence (P-value < 0.05 for 30, 45, 60, and 75 min; Figure 2A). This was consistent with the observed decrease in intensity of bands corresponding to cysteine protease 1 and 2 (CP1 and CP2, synonymous to EhCP-A1 and EhCP-A2, respectively) in the zymogram (Figure 2B), and by an ImageJ (http://rsbweb. nih.gov/ij/) analysis of the same bands that indicated an approximately 33% decrease in band areas. The band for CP5 (EhCP-A5) and control bands (arrow head) were comparable between HM-1 and MTZR. The decrease in cytopathy in MTZR was also indicated by the lower the number of engulfed beads ( Figure 2C) and percentage of cells with engulfed beads in the phagocytosis assay ( Figure 2D). MTZR also consumed glucose at a reduced rate compared to HM-1 (P-value < 0.05, Figure 2E). In addition, MTZR adhered to plastic surfaces, and to fibronectinand collagen-coated plates at relatively comparable levels, but not as efficiently as HM-1 (Supplementary Figure S2).

Transcriptomic Analysis of MTZR
Our expression data showed variations in the number of differentially expressed genes between MTZR and HM-1 cultured with and without MTZ (GEO ID: GSE35990). To identify genes whose expression was most highly affected, a general filtering of the expression data, as indicated in Materials and Methods, was used. This resulted in the identification of 142 probe sets for MTZR, 58 probe sets for MTZR (−), and 37 probe sets for HM-1 (+) that showed at least 3-fold changes compared to HM-1 (−) (Figure 3A, Table S3). In MTZR, 95 (66.9%) of the genes were upregulated, while 47 (33.1%) genes were downregulated (Tables 2, 3). Only 55 (38.7%) genes were annotated after conducting BLAST search, while 87 (61.3%) genes encoded for proteins of unknown function, i.e., hypothetical proteins (Table  S4). The most upregulated gene was a hypothetical protein (HP) (EHI_006850) with some similarity to a zinc finger protein in E. histolytica (E-value = 1.0e −46 ), designated as HP1 in this study. Its expression in MTZR was increased by 36.9 folds. Genes for DNA polymerase (EHI_164190), iron-sulfur flavoprotein (ISF) (EHI_025710), AIG1 family protein [XP_648192 (522.m00018)], and four other HPs were also upregulated by more than 10 folds. The most downregulated genes were three leucine-rich proteins [EHI_033560, EHI_077280, EHI_124070 (371.m00031)] which had fold changes of more than 160, and two AIG family proteins [XP_648014 (628.m00011), EHI_180390] and two CPs (EHI_160330, EHI_121160) genes that were downregulated by at least 5 folds. The two CPs identified had high similarity to CP7 (synonymous to EhCP-B1) (81%) based on amino acid sequence (Tillack et al., 2007). The differential expression of 16 genes was validated by qRT-PCR, using rpoII as an endogenous control ( Table 4). Pearson's correlation coefficient (0.81) between microarray and qRT-PCR results indicated a reasonable degree of concordance between the fold-changes of the genes selected for this comparison.

Go Classification, Features, and Clustering on the Genome
To have a general view of the cellular functions regulated as a result of MTZ resistance, we performed a GO enrichment analysis. A total of 132 GO terms were mined from Amoeba DB representing 63 unique GO IDs (Table S5). However, based on Fisher's exact test, only 24 terms were significantly enriched in MTZR ( Table 5). Genes involved in GTP binding and proteolysis were the most highly enriched, followed by genes associated with cysteine-type peptidase activity, small GTPase-mediated signal transduction, oxidoreductase activity, protein kinase activity, and electron carrier activity. Roughly, the enriched terms could be classified into four functional categories, namely (1) nucleotide binding, (2) metabolism, (3) oxidative stress response, and (4) signal transduction.
Only a fraction of the differentially expressed genes in MTZR (about 5%) encoded for membrane proteins, based on positive results for transmembrane and signal peptide predictions (Table S6). Only 21 genes (14.8%) contained at least one transmembrane domain. Eleven of these were upregulated, while 10 were downregulated. Proteins encoded by 25 genes (17.6%) also had predicted signal peptides. Only seven genes have both transmembrane domains and signal peptides. Only one, however, was annotated as tyrosine kinase (EHI_118410), while the rest were all genes for HPs.
Most of the modulated genes were encoded on different scaffolds. Thirty-four genes (30.6% of EHI-probe sets), however, formed 15 sets of dyads or triads that occupied the same contig (Table S7). Included were genes for four AIG1 family proteins (EHI_126550, EHI_126560, EHI_176700, EHI_176580) and two Ras family proteins (EHI_045450, EHI_045600). In four sets of genes, the members were found to be located on opposite strands, but the regulation of their expression was the same and their foldchanges were comparable. For some genes found on the same strand, the distance between the genes was less than 1 kb. It was therefore possible for these genes to be co-regulated by the same transcription factor.

Functional Analysis
To study further the relationship between the expression of some of the modulated genes in MTZR with drug resistance, generation of stable transformants overexpressing HA-tagged HP1, ISF1, ISF2, and ISF4 was performed ( Figure 4A). The mRNA level of these proteins was estimated by qRT-PCR. The levels of mRNA of HP1 and ISFs were about 5.5 or 2.5 to 3 fold, respectively, higher in MTZR compared to their expression in HM-1 transfected with mock vector (pEhEx-HA). Transformants cultured with 10 µM G418 were challenged with MTZ for 48 h (Figure 4B). The highest concentration shown was based on the IC 50 of MTZ against HM-1, performed under the same conditions (6.5 µM). No significant difference in percent survival was observed at all concentrations tested between the overexpressors and control ( Figure 4B). Neither was significant difference in survival against hydrogen peroxide, observed, among the transformants ( Figure 4C).

Differences in MTZR Generation in This Study and Previous Works
Currently, only 4 studies have been published regarding E. histolytica and its response to or resistance to MTZ (Samarawickrema et al., 1997;Wassmann et al., 1999;Leitsch et al., 2007;Tazreiter et al., 2008). These studies either exposed wild-type cells to serum level concentrations of the drug (70-100 µM) or generated drug resistant trophozoites, similar to what we did in this study. However, the MTZ concentrations used in these studies and ours differ, which might explain some of the differences observed. The concentration used here is slightly higher compared to those used by Samarawickrema et al. 10 µM (Samarawickrema et al., 1997), but significantly lower than those used by Wassmann et al. 40 µM (Wassmann et al., 1999), Tazreiter et al. 50 µM (Tazreiter et al., 2008), and Leitsch et al. 50 µM (Leitsch et al., 2007). Our strain adapted more quickly to 1 µM MTZ compared to the strain used by Samarawickrema (7 vs. 38 days) but overall the length of time used to generate MTZR was comparable. For reasons we do not understand, our attempts to further increase the level of resistance failed, and no resistance above the concentration used here could be established. We could  not culture MTZR consistently at drug concentrations higher than 12 µM. In contrast, strains of Trichomonas and Giardia can be readily adapted to grow at 584 µM and 115 µM MTZ, respectively (Kulda et al., 1984;Townson et al., 1992). Our results are therefore consistent with previous reports indicating that MTZ resistant strains are more difficult to generate in E. histolytica (Samarawickrema et al., 1997). While the concentration we used is several times lower than serum levels after MTZ treatment (Van Oosten et al., 1986), the changes we observed may still provide insights as to what occurs in the parasite in vivo. It was previously reported that in an abscess, parasites encounter significantly lower levels of drug than those found in the serum as drug penetration is limited by poor perfusion and mechanical barriers such as fibrin clots and the abscess wall (Sirinek, 2000).

Decrease in Growth and Reversibility of Resistance in MTZR
Formation of drug resistance is often accompanied by fitness costs (Andersson and Hughes, 2010). One particular cost, reduction in growth rate, is the main parameter that determines the rate by which resistance develops and level of resistance reached. Similar to resistant strains previously generated (Samarawickrema et al., 1997;Wassmann et al., 1999), the doubling time of MTZR was longer compared to HM-1. Fitness cost in growth associated with drug resistance are well documented in mammalian cells, e.g., lymphoma and cancer cells (Lee, 1993;Bishop et al., 2001). It should be noted, however, that the growth kinetics of MTZR in the absence of drug pressure improved over time. This indicates that some of the drugselection-induced changes developed by MTZR were temporal and reversible. This also indicates that cells in the culture may have varying degrees of resistance. In addition, MTZR (−) that was re-exposed to MTZ showed significantly lower growth compared to MTZR, but slightly higher compared to HM-1, which may indicate that some drug-selection-induced changes in MTZR were still operative even after the removal of drug pressure.

Cross-Resistance in MTZR
The resistance level displayed by MTZR is significantly lower compared to previous reports: the IC 50 value of 12 µM in this study compared to 40 µM in the previous study (Samuelson, 1999;Wassmann et al., 1999). However, its survival in otherwise lethal concentration of MTZ for long periods of incubation, does indicate reduced sensitivity and resistance. The cross-resistance of MTZR to ornidazole and tinidazole is not surprising, because both are 5-nitroimidazoles like MTZ. It is therefore likely that their mode of action and mechanism of resistance were similar (Pasupuleti et al., 2014). T. vaginalis isolates that were characterized as having very high resistance to MTZ were similarly insensitive to tinidazole (Narcisi and Secor, 1996). In the case of paromomycin, emetine, chloroquine, and hydrogen peroxide, it is possible that the reduced sensitivity was nonspecific and could be due to an increased capacity of MTZR to deal with stress. Studies have shown that cell lines resistant to one drug are often cross-resistant to closely related drugs or to unrelated compounds when their mode of action is the same or when they carry membrane alterations (Upcroft and Upcroft, 1993).

Decreased Adhesion, Phagocytosis, Cytolysis, and Virulence in MTZR
In E. histolytica, establishment and outcome of infection depend heavily on adhesive and invasive capacity (Flores-Romo et al., 1997). Recently, drug resistance has been shown to modulate these factors and vice versa (Giha et al., 2006). The decreased rate of adhesion observed in MTZR indicated that MTZ selection negatively affected the synthesis or expression of surface adhesion molecules. While other studies have shown that cell adhesion is a key determinant in drug resistance (Shain and Dalton, 2001), this apparently was not the case in MTZ resistance. Inhibition of adherence in this parasite has previously been shown to decrease host cell cytotoxicity (Ravdin et al., 1985) and experimental evidence strongly suggests the involvement of surface adhesins in phagocytosis (Heron et al., 2011). As shown in the results, MTZR also had decreased cytopathy and phagocytosis (Figures 2A,C). The decreased band intensity for CP1 and CP2 in the zymogram may partially explain this phenotype. CP1 has been shown to be upregulated in a mouse model of amebic colitis following invasion (Gilchrist et al., 2006), while CP2 has been reported to contribute to intestinal damage and liver abscess formation (Hellberg et al., 2001). Our microarray data, however, did not show any significant repression of these genes, except for two uncharacterized CPs that were repressed by more than 7 folds. These CPs had high similarity (81%) to CP7 (also called EhCP-B1) based on amino acid sequence (Tillack et al., 2007). No significant change was observed in the expression of the genes known to be involved in the inactivation and intracellular trafficking of CPs such as intrinsic inhibitors of CPs (Sato et al., 2006) and cysteine protease binding protein family 1  and other lysosomal hydrolase carriers (Furukawa et al., , 2013Marumo et al., 2014). It is also possible, however, that resistance to MTZ may not necessarily result to decreased virulence in vivo.

Implications from Transcriptomic Profiling of MTZR
Our transcriptomic profiling suggests that the transcriptomic changes we observed in MTZR strain cultured with MTZ were a combination of the changes responsible for MTZ resistance and those associated with the adaptation to stresses caused by the exposure to MTZ. However, we believe that among the observed changes, the changes that were specific to MTZR strain with MTZ and absent in HM-1 with MTZ were responsible, at least in part, for resistance, rather than adaptation. This is based on the following observations: the transcriptomic profile of E. histolytica HM-1 cultured with MTZ was remarkably different from that of MTZR strain cultured with MTZ, while MTZR strain cultured with or without MTZ showed similar transcriptomic profiles.
In other words, a larger number of common changes in gene expression were identified between MTZR strain cultured with MTZ and those without MTZ when compared between MTZR strain cultured with MTZ and HM-1 cultured with MTZ. GO enrichment of MTZR-associated transcriptomic changes revealed regulation of genes encoding for proteins potentially involved in nucleotide binding, metabolism, oxidative stress response, and signal transduction. These are apparently the processes most influenced by MTZ resistance. However, since most of the genes are without GO annotations, it is likely that additional processes are involved. When we examined the genomic locations of these genes, we found that 30% of them formed dyads or triads that occupied the same contig. While some genes forming dyads and triads are apart by more than 10 kb, the distance between some of the adjacent genes in dyads and triads was less than 1 kb. Thus, co-expression of some genes in dyads and triads in MTZR could be attributed to a shared regulatory system. Consistent with previous findings, no significant difference in PFOR mRNA levels was observed between MTZR and HM-1 (Samarawickrema et al., 1997;Wassmann et al., 1999;Tazreiter et al., 2008). This observation could mean 3 things. First, it's likely that the level of resistance achieved in this study was not sufficient to downregulate PFOR, although a similar finding was reported in strains resistant to 40 µM MTZ (Wassmann et al., 1999). Second, since E. histolytica lacks pyruvate dehydrogenase (Clark et al., 2007;Tazreiter et al., 2008) or pyruvate decarboxylase (Lo and Reeves, 1978), the significance of PFOR cannot be overstated. The lack of an alternative enzyme to compensate for its downregulation makes this enzyme difficult to repress without severe consequences to the ameba. Third, in contrast to Giardia and Trichomonas, lack of PFOR downregulation may indicate the existence of an alternative mechanism for MTZ resistance in this parasite. It should be noted, however, that clinically resistant isolates of Trichomonas have been reported with no decrease in PFOR activity (Müller and Gorrell, 1983).
Two previous studies have reported the upregulation of Fe-SOD and downregulation of Fdx in the MTZ resistant ameba (Samarawickrema et al., 1997;Wassmann et al., 1999). In the previous work, fold changes were determined based on enzyme activity assay or northern blot (Samarawickrema et al., 1997;Wassmann et al., 1999). Disparity, however, between changes in mRNA abundance and enzyme activity have been reported, and in some cases, increases in enzyme activity can exceed those of mRNA abundance (Glanemann et al., 2003). Interestingly, the lack of Fe-SOD modulation in cells exposed to MTZ was also reported by Tazreiter et al. (2008). On the contrary, in this work, neither gene was modulated significantly. It is possible that similar to PFOR, differences in resistance levels could explain the lack of modulation in the levels of Fe-SOD and Fdx. Other studies also indicated that Fe-SOD and Fdx are not tightly associated with MTZ resistance. It was shown that parasites in which Fdx was repressed did not display full resistance to MTZ (Rasoloson et al., 2002) and in bacteria, increased SOD activity was not always associated with MTZ resistance (Smith and Edwards, 1995).
In this study, genes involved in nucleic acid synthesis and binding, stress response, metabolism, and vesicular trafficking were mostly affected in MTZR. The significance of these genes relative to phenotypic changes and drug resistance will be discussed here through their known roles in E. histolytica or in other organisms, in regulating key cell functions. One should also note that since only a single MTZ resistant line was analyzed in this study, it is possible that MTZ resistance can be also conferred by the mechanisms that are independent of those described in this study.

Nucleic Acid Synthesis and Binding
In prokaryotes, upregulation of DNA polymerases imparts plasticity to their genome, allowing them to adapt in unfavorable environmental conditions (Joseph et al., 2006). It has been suggested that pathogenic bacteria utilize this adaptation to develop resistance against therapeutic agents (Karpinets et al., 2006). Since it has been shown that primary effect of MTZ is rapid inhibition of DNA replication (Ludlum et al., 1988;Sisson et al., 2000;Leitsch et al., 2007), the upregulation of a DNA polymerase (EHI_164190) in this study may indicate that this mode of action is operative in MTZ resistance. One hypothetical protein identified in this study (HP1, EHI_006850) has some similarity to zinc finger proteins, which are known to be involved in DNA recognition, RNA packaging, and transcriptional activation (Laity et al., 2001). It has an NAD+ binding pocket and a tetrachlorodibenzo-p-dioxin (TCDD)inducible poly(ADP-ribose) polymerase (PARP)-like domain that is involved in the attachment of ADP-ribose units to DNA-binding proteins (Ma et al., 2001). It's also known to be a regulatory component induced by DNA damage and an enzyme involved in DNA repair (Nguewa et al., 2006;Marchler-Bauer et al., 2011). The dUTP nucleotidohydrolase (EHI_072960) upregulated in MTZR may be involved in the removal of dUTP from the dNTP pool, de novo biosynthesis of dTTP, and DNA replication as reported (Richards et al., 1984;Gadsden et al., 1993). In cancer cells, high level of dUTPase is implicated as a mechanism of resistance to chemotherapy (Ladner et al., 2000).

Stress Response
ISFs that were modulated in MTZR are part of family of redoxactive proteins mainly found in anaerobic prokaryotes (Andrade et al., 2005). In NCBI, 18 ISF genes are annotated in the E. histolytica genome, although 50% of these are identical proteins (http://www.ncbi.nlm.nih.gov/). In this study 7 ISF probe sets were upregulated by at least 4 folds in MTZR. The overexpression of ISF1 (EHI_138480) was also observed in HM-1 (+), while ISF2 was upregulated in both MTZR (−) and HM-1 (+). These 2 ISFs were upregulated in HM-1 (+) by least 20 folds, compared to 6 to 11 folds in MTZR, hinting that these genes are involved in stress response. These two genes were also induced in E. histolytica cultured under L-cysteine deprived and oxidative stress conditions (Vicente et al., 2009;Husain et al., 2011), but were downregulated in a mouse model of intestinal amebiasis (Gilchrist et al., 2006). Our findings together with these reports suggest that induction of these ISFs, particularly ISF1, is an adaptation that may enhance survival when cells experience environmental stress.

Metabolism
Several metabolic genes were upregulated in MTZR. NADspecific glutamate dehydrogenase (EHI_075150) catalyzes the oxidative deamination of glutamate to α-ketoglutarate using NAD as cofactor (Plaitakis and Zaganas, 2001;Girinathan et al., 2014). Most likely, its upregulation is a part of stress response. In humans, its activity is regulated through ADP-ribosylation, and interestingly, Arf1 (EHI_189960) was also upregulated in this study (Herrero-Yraola et al., 2001). Chitinase (EHI_092100) (Escueta-De Cadiz et al., 2013) was also upregulated in MTZR and MTZR (−) by at least 3 folds, although its precise role in drug resistance is not known. In plants, its upregulation plays a defensive role against toxins and antibiotics (Gooday, 1999). Aldose reductase (EHI_029620), on the other hand, can reduce a broad spectrum of substrates including cytotoxic aldehydes (Ramana, 2011). Its upregulation in human liver cancer cells was associated with resistance to daunorubicin and was confirmed by sensitizing cells to the drug by the addition of aldose reductase inhibitors .
The decrease in cell growth observed in MTZR may be partially explained by the repression of phosphoserine aminotransferase (EHI_026360) (Ali and Nozaki, 2006;Mishra et al., 2010), which catalyzes the conversion of 3-phosphohydroxypyruvate to L-phosphoserine, the second step of phosphorylated serine biosynthetic pathway (Ali and Nozaki, 2006). L-serine has a well-recognized role in cell proliferation, providing precursors for amino acids, protein synthesis, and nucleotide synthesis Tabatabaie et al., 2010). The downregulation of this enzyme, therefore, may contribute to the decreased growth rate. The repression of lecithin: cholesterol acyltransferase (EHI_020250) may also be associated with decreased growth rate, because the enzyme is crucial in the maturation, remodeling, and metabolism of high density lipoprotein, HDL (Calabresi et al., 2011).

BspA1 Gene Family
Several genes encoding for leucine-rich repeat (LRR) proteins of the BspA family were also modulated in MTZR, and were the most downregulated genes with fold changes of more than 150. LRR is known to provide a structural framework for the formation of protein-protein interactions and serve as recognition motifs for surface proteins (Kobe and Kajava, 2001). BspA family protein was originally identified in Bacteroides forsythus and the LRR motifs of the BspA protein have been linked to cell binding to components of the extracellular matrix (Sharma et al., 1998). It is therefore possible that their downregulation may partially explain the reduced adhesion seen in MTZR. It was previously shown that an LRR BspA family protein in E. histolytica was located primarily on the plasma membrane, although it did not show whether the protein was involved in surface interactions (Davis et al., 2006).

AIG1 Gene Family
The E. histolytica AIG1 gene family consists of 47 members (Biller et al., 2010), and in this study 10 genes were found to be differentially expressed by 3 to 11 folds. Members of this protein family contain a domain found in Arabidopsis protein AIG1, which is likely involved in plant recognition and resistance to bacterial pathogens (Reuber and Ausubel, 1996). In E. histolytica, however, some AIG proteins appear to be linked with virulence. Comparative genomic analysis of Japanese clinical isolates indicated that an AIG family protein (EHI_176590) and the region around it were uniquely present in KU50, a strain isolated from a diarrheic patient, but absent in a non-pathogenic isolate called KU27 (Nakada-Tsukui et al., unpublished data). Interestingly, the same gene was downregulated by 2.7 folds in MTZR, and two adjacent genes encoding for AIG family proteins (EHI_176580 and EHI_176700) were also repressed by 3 to 4 folds. A previous study also reported the upregulation of an AIG family protein (EHI_180390) in a pathogenic strain of E. histolytica (Biller et al., 2010), although here the same gene was downregulated by 5 folds. It is therefore possible that the downregulation of these AIG family proteins might be related to the decrease in MTZR's cytopathogenicity.
Overall, we believe that these transcriptional changes participate in a global regulatory response resulting to the fitness costs and adaptive responses we observed in MTZR. The fact that more genes were upregulated than downregulated also suggests that increasing the amount of cellular components that play roles in stress response, DNA repair and others, is more important in developing low levels of resistance to MTZ than repression of proteins involved in activating the drug. Indeed, the most profound difference we observed between MTZR and HM-1 was the upregulation of genes for several ISF proteins. It should be noted, however, that more than half of the genes reported in this study are hypothetical. It is also likely that the effects of the repressed genes, although fewer in number, are more important in explaining the growth defect.

Functional Analysis of Potential Resistance Genes
In this study, we examined possible role of genes encoding HP1 and several ISF proteins in MTZ resistance. However, their episomal transfection and overexpression did not confer MTZ resistance to the extent observed in MTZR. Thus, at present, direct causal connection between HP1 and ISF overexpression and MTZ resistance has not been demonstrated. It is likely that in E. histolytica, resistance to the drug is multi-factorial in nature and requires the global regulation of genes involved in several key cell processes and not just DNA repair and stress response.
In conclusion, here we described the biological state of E. histolytica that was selected for MTZ resistance in vitro. Phenotypic and transcriptional profiling revealed fitness costs and adaptive responses that we could associate with drug resistance. Some of the consequences of resistance included decreased cell growth and virulence, which have significant implications on host-pathogen interaction and infection outcome. Comparative transcriptional analysis, on the other hand, revealed genes not previously associated with MTZ resistance in this parasite. The list of differentially transcribed genes also did not indicate in MTZR a reduced capacity to activate MTZ, but hinted on adaptations to tolerate the lethal effects of MTZ. Identifying the signaling mediators for the proteins that these genes encode, however, requires further investigation. Finally, while our data showed that high levels of MTZ resistance cannot be readily induced in vitro, the data presented here may still help researchers in better understanding the mechanisms involved should clinical cases of high levels of resistance are reported in the future.