Pharmacogenomic Profile and Adverse Drug Reactions in a Prospective Therapeutic Cohort of Chagas Disease Patients Treated with Benznidazole

Chagas disease remains a major social and public health problem in Latin America. Benznidazole (BZN) is the main drug with activity against Trypanosoma cruzi. Due to the high number of adverse drug reactions (ADRs), BZN is underprescribed. The goal of this study was to evaluate the genetic and transcriptional basis of BZN adverse reactions. Methods: A prospective cohort with 102 Chagas disease patients who underwent BZN treatment was established to identify ADRs and understand their genetic basis. The patients were classified into two groups: those with at least one ADR (n = 73), and those without ADRs (n = 29). Genomic analyses were performed comparing single nucleotide polymorphisms between groups. Transcriptome data were obtained comparing groups before and after treatment, and signaling pathways related to the main ADRs were evaluated. Results: A total of 73 subjects (71.5%) experienced ADRs. Dermatological symptoms were most frequent (45.1%). One region of chromosome 16, at the gene LOC102724084 (rs1518601, rs11861761, and rs34091595), was associated with ADRs (p = 5.652 × 10−8). Transcriptomic data revealed three significantly enriched signaling pathways related to BZN ADRs. Conclusions: These data suggest that part of adverse BZN reactions might be genetically determined and may facilitate patient risk stratification prior to starting BZN treatment.


Introduction
The protozoan parasite Trypanosoma cruzi is the causal agent of Chagas disease (CD), a neglected tropical disease. Although the global prevalence of CD has fallen from 18 million in 1991 to 5.7 million in 2010 due to effective vector control programs [1][2][3], it still represents a major concern in endemic settings. Most people infected with T. cruzi live in Mexico, Central America, and South America [3]. For example, in Brazil, one million people are currently estimated to be infected with the parasite. It is notable then, given the substantial disease burden, that therapeutic options are limited, and benznidazole (BZN) is the only drug available for treatment in Brazil [4].
BZN is poorly tolerated due to frequent adverse drug reactions (ADRs), with 20-25% of adult patients unable to complete the standard 60-day treatment course [5][6][7][8]. Different Int. J. Mol. Sci. 2021, 22,1960 2 of 13 BZD schemes have been proposed by observational studies, including shorter courses and even intermittent treatment. These may lead to comparable clinical and serological outcomes and, especially, to fewer ADRs and better compliance [9,10].
Despite this major limitation, few studies have explored the biological basis for increased ADR susceptibility. It was previously shown that the development of BZN-related ADRs is independent of dose [11]. In addition, BZN ADRs are negatively associated with Black race and positively associated with female sex and schooling [12]. Notably, however, only one study investigated the genetic basis of BNZ-associated ADRs, in which patients carrying HLA-B*3505 were more likely to abandon treatment [13].
Therefore, a better understanding of the pharmacogenetic aspects of BZN could significantly impact the management of patients. The purpose of this study was to conduct a multi-omics evaluation of BZN ADRs and describe genetic and transcriptional features associated with ADRs among individuals undergoing BZN treatment in a therapeutic cohort of patients chronically infected with T. cruzi.

Results
A total of 102 subjects were enrolled to receive BZN. Of these, 73 (71.5%) experienced at least one ADR and 22 (21.6%) did not complete the full treatment scheme due to severe ADRs. Table 1 summarizes their clinical and epidemiological characteristics. None of the variables evaluated were associated with ADRs, except for the Chagas disease form, with patients presenting with the cardiac form having a lower rate of ADRs.
Of the 102 patients, dermatological, gastrointestinal, and neurological ADRs were reported by 46 (45.1%), 30 (29.4%), and 28 (27.4%) patients, respectively. Treatment interruption was more frequent among patients who reported an ADR during the first 15 days of treatment with relative risk (RR) of 14.1, 95% confidence interval (CI) 2.0-100. Dermatological ADRs were strongly associated with treatment cessation compared to individuals who did not develop this type of ADR (RR 5.1, 95% CI 2.2-11.9). Other less frequent ADRs included arthralgia (12 subjects, 11.8%), reduced renal function (5 subjects, 4.9%), myelotoxicity (4 subjects, 3.9%), and hepatotoxicity (4 subjects, 3.9%). No ADRs required hospital admission, and most were controlled with a short course of low-dose corticosteroids, antihistamines, and analgesics. Sixteen patients (15.6%) reported ADRs on all study visits during treatment. Figure 1 summarizes the frequency of ADRs over time. Dermatological ADRs were more common at the beginning of treatment.

Genome-Wide Association Study Results
Of the initial 96 genotyped samples, 6 were excluded because they showed high kinship with other samples. A total of 90 patients were included in the Genome-Wide Association study (GWAS) analysis: 61 with at least one ADR (group gADR) and 29 without any ADRs (group gNADR). No single nucleotide polymorphism (SNP) reached our predefined genome-wide significance threshold. Nonetheless, a number of markers showed a strong association with ADR occurrence. Table 2 presents the 30 SNPs with the lowest p-values.  On chromosome 16, rs1518601, rs11861761, and rs34091595 SNPs (p-value = 5.652 × 10 −8 ) were associated with the occurrence of ADRs. Figure 2 describes the region in detail, with many adjacent SNPs in linkage disequilibrium, categorized by r 2 values between 0.6 and 0.8. These SNPs are located at the gene LOC102724084. On chromosome 3, the SNP rs3968927 (p-value = 2.945 × 10 −7 ) was in strong linkage disequilibrium with other six SNPs, with a p-value below 1.95 × 10 −6 , and several other dispersed across the genes LMOD3, ARL6IP5, FRMD4B, UBA3, TMF1, MIR3136, EOGT, and FAM19A4. On chromosome 15, there were five highly associated SNPs in strong linkage disequilibrium (p-value = 6.9 × 10 −7 ) and in proximity to the genes GABRG3 and OCA2. However, none of these genes have a clear biological relationship with ADRs, drug metabolism, CD, or direct interactions with the immune system or inflammatory response. For a better representation of GWAS results and to assess the quality of the data, a Manhatan plot and a Quantile-quantile plot are provided in the Supplementary Material ( Figures S1 and S2).

Transcriptomic Results
To understand the impact of BZN on gene expression, we compared samples prior to and at the end of BZN treatment (T0d × T60d). A total of 39 patients, 26 with at least one ADR (gADR) and 13 without any ADR (gNADR), were included in this analysis. There were 11 differentially expressed genes (DEGs) in the gNADR and 173 in the gADR group (comparing before and after BZD use). There were no DEGs in common (Table S1). The gNADR group presented two activated signaling pathways that were completely independent of the 25 pathways in the gADR group (Table S2). The Ingenuity Pathway Analysis (IPA) score value classification identified three main pathways in the gADR group, different from those found in the gNADR group.
Patients with dermatological ADRs had the most significantly enriched pathway (28 upregulated genes of 34 genes comprising the pathway, 82%). This was clearly associated with the biology of the skin (Figure 3). Patients with dermatological ADRs had the most significantly enriched pathway (28 upregulated genes of 34 genes comprising the pathway, 82%). This was clearly associated with the biology of the skin (Figure 3). Gastrointestinal ADRs, the second most common in this cohort, were related to a pathway linked with gastrointestinal functioning and disease. This pathway had a score of 39 and consisted of 34 genes, 79.4% of which were upregulated, indicating a strong relationship with gastrointestinal clinical manifestations developed in patients with ADRs ( Figure 4).
Neurological manifestations were associated with a short pathway consisting of two genes. However, the GABRB1 gene had a log Fold change (logFC) above three, the highest gene upregulation observed, and a 200% increase compared to the pretreatment period ( Figure 5). Gastrointestinal ADRs, the second most common in this cohort, were related to a pathway linked with gastrointestinal functioning and disease. This pathway had a score of 39 and consisted of 34 genes, 79.4% of which were upregulated, indicating a strong relationship with gastrointestinal clinical manifestations developed in patients with ADRs ( Figure 4).
Neurological manifestations were associated with a short pathway consisting of two genes. However, the GABRB1 gene had a log Fold change (logFC) above three, the highest gene upregulation observed, and a 200% increase compared to the pretreatment period ( Figure 5).

Discussion
To the best of our knowledge, this is the first paper to describe genetics and transcriptomics data in relation to BZN-induced ADRs in CD patients. In our study, most patients (71.6%) experienced ADRs, in line with previous findings [14][15][16]. Dermatological events were the most frequent ADR observed, representing half of the total ADRs. This value is higher than that previously reported [15], likely due to differences in the clinical assessment procedures. One in five patients did not complete treatment due to ADRs. Dermatological events were most strongly associated with incomplete treatment (relative risk of 5.1 compared to those without this type of ADR) and tended to occur in the first

Discussion
To the best of our knowledge, this is the first paper to describe genetics and transcriptomics data in relation to BZN-induced ADRs in CD patients. In our study, most patients (71.6%) experienced ADRs, in line with previous findings [14][15][16]. Dermatological events were the most frequent ADR observed, representing half of the total ADRs. This value is higher than that previously reported [15], likely due to differences in the clinical assessment procedures. One in five patients did not complete treatment due to ADRs. Dermatological events were most strongly associated with incomplete treatment (relative risk of 5.1 compared to those without this type of ADR) and tended to occur in the first

Discussion
To the best of our knowledge, this is the first paper to describe genetics and transcriptomics data in relation to BZN-induced ADRs in CD patients. In our study, most patients (71.6%) experienced ADRs, in line with previous findings [14][15][16]. Dermatological events were the most frequent ADR observed, representing half of the total ADRs. This value is higher than that previously reported [15], likely due to differences in the clinical assessment procedures. One in five patients did not complete treatment due to ADRs. Dermatological events were most strongly associated with incomplete treatment (relative risk of 5.1 compared to those without this type of ADR) and tended to occur in the first two weeks of treatment. Aldasoro et al. (2018) reported that 80% of BZN-related ADRs occur within the first 30 days of treatment. Complaints such as arthralgia and neurotoxicity were more frequent during the last weeks of the standard 60-day treatment regimen, again in agreement with other studies [14][15][16]. We did not find any difference regarding demographic characteristics and ADR frequency, in disagreement with what was reported in other studies [12], which found a higher frequency of ADRs and BZD discontinuation in women and in those who graduated from elementary school. Maybe our study population represents a cosmopolitan region where the population's demographic heterogeneity is less evident.
Although the sample size was small for a GWAS, we detected three SNPs (rs1518601, rs11861761, and rs34091595) that were highly associated with the occurrence of ADRs. The SNPs are intronic variants located at the gene LOC102724084 on chromosome 16. Unfortunately, the function of this gene is unknown, so we cannot verify the biological plausibility of this finding. There were several other adjacent SNPs in strong linkage disequilibrium located in the same region, many with p-values less than 10 −5 , suggesting that this is a true finding. There is only one previous study conducted in Bolivia that attempted to detect genes associated with BZN ADRs [13]. In that cohort, the gene HLA-B 3505 was found to be associated with adverse reactions. However, this human leukocyte antigen (HLA) is not frequent in the Brazilian population, and we did not confirm their findings.
We identified substantially different gene expression profiles between patients with and without adverse reactions. In the ADR group, there were more differentially expressed genes at the end of treatment compared to pretreatment. Dermatological manifestations were correlated with a pathway in which 82.4% of the genes were upregulated. These genes are associated with many dermatological conditions. They have previously been associated with more severe clinical expression of dermatitis [17][18][19][20] and advanced stages of atopic dermatitis [21], increased epidermal hyperplasia [22,23], cyclic itching associated with epidermolytic hyperkeratosis [24], lichen planus [25,26], erythema nodosum [27], and ichthyosis with confetti [28]. Many of the dermatological ADRs reported by patients can be directly associated with this pathway. The pathway linked to gastrointestinal manifestations showed significant upregulation of 79.4% of the genes. These genes have been linked to gastrointestinal functioning and disease, for example, an increased rate of liver injury [29] and biliary cirrhosis [30]. This can be directly (or indirectly) associated with the gastrointestinal ADRs reported by patients. The pathway associated with neurological ADRs consisted of two genes, GABRG1 and GABRG2. These are associated with neuropathic pain [31], hyperalgesia [32], migraines [32,33], and tension headache [31,32], similar to the clinical ADRs reported by our patients.
The limitations of the current study were the small sample size for a GWAS and not enough RNA for transcriptomic assays and gene expression validation by RT-qPCR.

Material and Methods
We conducted a 12-month prospective treatment study of adults with chronic CD at two large specialized centers, the Hospital das Clínicas and the Institute of Infectology "Emilio Ribas" in Sao Paulo, Brazil. Eligible patients were aged between 18 and 65 years, seropositive, and PCR-positive for T. cruzi. Chagas disease clinical form and treatment indication with BZN were defined according to the Brazilian Ministry of Health Consensus statement [34]. Chronic comorbidities were accessed by self-reporting during study interviews and confirmed by reviewing each patient's medical record. These conditions were diagnosed in agreement with clinical societies' guidelines. Exclusion criteria were previous BZN use; advanced Chagas cardiomyopathy (New York Heat Association gradation III or IV); megaesophagus grade IV (Rezende's classification) [35]; advanced megacolon, hepatic, or renal failure; and inability to attend follow-up appointments. Participants received BZN (5 mg/kg/day) for 60 days, at a maximum dose of 300 mg/day. They underwent regular clinical evaluations that consisted of a medical history and physical examination, a biochemistry panel, and serum analysis. Clinical evaluations were performed 30 days before treatment (-T30d), on the day before BZN initiation (T0d), and on days 15 (T15d), 30 (T30d), and 60 (T60d) after starting treatment. Further study visits were conducted at 6 (T6M) and 12 (T12M) months following treatment initiation. Blood samples were obtained during all follow-up appointments.
ADRs were assessed at T15d, T30d, and T60d by specifically asking about dermatological (pruritus, rash, urticaria), gastrointestinal (nausea, vomiting, abdominal pain, reduced appetite), neurological (headache, paresthesia, or dysesthesia), and joint symptoms, as well as by monitoring for laboratory abnormalities. Hepatotoxicity was defined as an increase in aminotransferases greater than three times the upper normality limit. Nephrotoxicity was defined as an increase in creatinine of more than 1.5 times the upper limit and myelotoxicity as a decrease in hemoglobin, leukocytes, or platelets of more than 10% from pretreatment levels. Data obtained were uploaded on the Research Electronic Data Capture (REDCap) database specifically designed for the present study [36].
Patient characteristics were compared between groups using the chi-square or Fisher's exact test for categorical variables and a Student's t-test for continuous variables. Relative risk and respective 95% confidence intervals were calculated. All statistical tests were two-tailed with an α error of 0.05. Statistical analyses were conducted using Stata (v.13.0).

Genomic and Transcriptomic Sample Extraction
Blood samples were centrifuged, and separate aliquots of plasma and buffy coat/red cells were taken and stored at −20 • C. Genomic DNA was extracted using a QIAamp DNA mini kit (Qiagen, Hilden, Germany), and the concentration was measured using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Samples were accepted for genotyping if they presented a 260/280 ratio between 1.7 and 1.8 or a 260/230 ratio between 2.0 and 2.2. In addition, DNA concentration was required to be above 10 ng/µL. DNA integrity was accessed using an agarose 1.5% gel. A total of six samples did not pass the pre-analytic standards, and genotype testing was not performed. Overall, 96 patients were tested using the Affymetrix assay.
Thirty-nine patients were selected for having transcriptome data assayed using samples collected at T0d and T60d. All 13 patients without ADRs (gNADR) were chosen along with 26 individuals, randomly selected from the 73, who presented with ADRs during treatment (gADR). Whole-blood samples collected using PAXgene tubes were extracted using a QIAcube PAXgene blood RNA kit (Qiagen, Hilden, Germany). The DNase I treatment step was performed during the procedure to remove contaminating genomic DNA. RNA integrity (RIN) and concentration were measured using Agilent Bioanalyzer 2100 with a Nanochip 6000 kit (Agilent Technologies, Santa Clara, CA, USA). All samples had RIN and RNA concentrations higher than 4.5 and 5 ng/µL, respectively.

Transcriptomic Processing and Analyses
Samples were subjected to amplicon RNA sequencing using the Ion Proton Platform. This approach has excellent sensitivity to detect messenger RNA (mRNA) at low abundance. RNA (30 ng) was reverse-transcribed and barcoded. Complementary DNA (cDNA) libraries were generated using an Ion AmpliSeq Transcriptome Human Gene Expression Kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer's protocol for blood samples [37]. Library concentrations were determined by qPCR using an Ion Library TaqMan Quantitation Kit (Thermo Fisher Scientific, Waltham, MA, USA). Template reactions were carried out using an Ion PI Hi-Q OT2 200 Kit and then loaded onto Ion PI chips v3 using an Ion PI Hi-Q Sequencing 200 Kit (Thermo Fisher Scientific, Waltham, MA, USA). All samples were sequenced using the Ion Proton System, and samples with less than 2 × 10 6 reads were resequenced. Ion Torrent Suite 5.0 was used for assessing the quality of the libraries and sequencing runs. Base-calling results for each sample were aligned to a standard reference file (hg19_ampliseq_transcriptome_ercc_v1.fasta) on Ion Torrent Suite (v. 5.04). Files were exported to R (v. 3.6.2). Data were normalized and filtered using the edgeR package. For all individual comparisons, we calculated the log fold-change (logFC) and associated p-value for each gene. Genes that were up-or downregulated by a logFC of 1.25 or more and with a p-value of less than 0.005 were considered significant. Differentially expressed genes were loaded into Ingenuity Pathway Analysis (IPA) (v. 2.4) (Qiagen, Hilden, Germany) to identify over-represented signaling pathways. IPA identifies the key functional and disease associations for given genes using a score value classification.

Genomic Processing and Analyses
Sample DNA (10 ng) was amplified, enzymatically fragmented, dried, and resuspended following the Axiom 2.0 reagent kit protocol (Affymetrix, Santa Clara, CA, USA) [38]. After the hybridization step, samples were processed in GeneTitan equipment (Affymetrix, Santa Clara, CA, USA), where probes identified 842,007 single-nucleotide polymorphisms (SNPs) per individual. In Axiom Suite Analyzes (Affymetrix, Santa Clara, CA, USA), genotype data were normalized, filtered, and exported. SNP cleaning was performed using PLINK (v. 1.9). Based on genetic population study models, a total of 30,316 SNPs did not follow quality control precepts and were removed. SNP imputation was performed on the Michigan Imputation Server platform (https://imputationserver.sph.umich.edu/ (accessed on 13 November 2020)) using the Haplotype Reference Consortium (HRC) reference panel. To protect against association bias due to the population genetic structure, we used the first 4 principal components as covariates in our association models. Polymorphic markers were filtered by Hardy-Weinberg equilibrium using a p-value threshold of 1 × 10 −8 . The generated data were analyzed using PLINK, applying the allelic association model to gADR and gNADR groups and correcting for the population structure using the first two principal components. A p-value of less than 5 × 10 −8 was considered genome-wide significant.

Conclusions
In conclusion, ADRs, especially dermatological reactions occurring within the first 15 days, were the main reason for abandoning BZN treatment. Although the results are preliminary, it was possible to identify a number of SNPs associated with the manifestation of ADRs. In addition, we were able to demonstrate a relationship between the signaling pathways identified and the clinical manifestations in the gADR group, namely dermatological, gastrointestinal, and neurological ADRs. These findings provide new preliminary perspectives for the management of BNZ-related ADRs, since they describe a future research roadmap that could lead to the use of these genetic findings in CD patients' treatment. For example, they may allow better phenotyping of CD patients before initiation of BZN and stratification of the risk of adverse reactions, allowing greater attention to be given to high-risk patients, potentially facilitating treatment adherence.

Institutional Review Board Statement:
The study protocol was approved, in accordance with the Declaration of Helsinki, by the local ethics committee of the Hospital das Clínicas, University of São Paulo, and the Institute of Infectiology "Emilio Ribas" (CEP 042/12).

Informed Consent Statement:
Before recruitment, all patients involved in the study provided written informed consent.

Conflicts of Interest:
All authors declare no conflict.