Genetic analysis of the orthologous crt and mdr1 genes in Plasmodium malariae from Thailand and Myanmar

Background Plasmodium malariae is a widely spread but neglected human malaria parasite, which causes chronic infections. Studies on genetic polymorphisms of anti-malarial drug target genes in P. malariae are limited. Previous reports have shown polymorphisms in the P. malariae dihydrofolate reductase gene associated with pyrimethamine resistance and linked to pyrimethamine drug pressure. This study investigated polymorphisms of the P. malariae homologous genes, chloroquine resistant transporter and multidrug resistant 1, associated with chloroquine and mefloquine resistance in Plasmodium falciparum. Methods The orthologous P. malariae crt and mdr1 genes were studied in 95 patients with P. malariae infection between 2002 and 2016 from Thailand (N = 51) and Myanmar (N = 44). Gene sequences were analysed using BioEdit, MEGA7, and DnaSP programs. Mutations and gene amplifications were compared with P. falciparum and Plasmodium vivax orthologous genes. Protein topology models derived from the observed pmcrt and pmmdr1 haplotypes were constructed and analysed using Phyre2, SWISS MODEL and Discovery Studio Visualization V 17.2. Results Two non-synonymous mutations were observed in exon 2 (H53P, 40%) and exon 8 (E278D, 44%) of pmcrt. The topology model indicated that H53P and E278D were located outside of the transmembrane domain and were unlikely to affect protein function. Pmmdr1 was more diverse than pmcrt, with 10 non-synonymous and 3 synonymous mutations observed. Non-synonymous mutations were located in the parasite cytoplasmic site, transmembrane 11 and nucleotide binding domains 1 and 2. Polymorphisms conferring amino acid changes in the transmembrane and nucleotide binding domains were predicted to have some effect on PmMDR1 conformation, but were unlikely to affect protein function. All P. malariae parasites in this study contained a single copy of the mdr1 gene. Conclusions The observed polymorphisms in pmcrt and pmmdr1 genes are unlikely to affect protein function and unlikely related to chloroquine drug pressure. Similarly, the absence of pmmdr1 copy number variation suggests limited mefloquine drug pressure on the P. malariae parasite population, despite its long time use in Thailand for the treatment of falciparum malaria.

is higher than previously assumed [3]. There is limited information on the biology and molecular genetics of P. malariae as well as on anti-malarial drug resistance in this species for which long-term in vitro culture methods are lacking. Small scale in vitro drug susceptibility testing and clinical efficacy [4] showed that P. malariae in Thailand responds well to pyronaridine [5] and artesunate [6] but that fever clearance time was delayed after treatment with chloroquine [7]. However, P. malariae exhibits a 72-hour erythrocytic stage, which suggests that it may require a longer time to clear parasites after treatment. Therefore, the extended clearance time may not be an indicator for chloroquine resistance in P. malariae [8].
In vitro drug testing in P. malariae is further complicated by the low parasitaemias present in patients and the high proportion of mixed infections with other Plasmodium species. An alternative approach for monitoring antimalarial drug resistance in P. malariae is the assessment of polymorphisms in anti-malarial drug target genes. Since the molecular markers for drug resistance are not well characterized in P. malariae, the orthologous gene markers associated with drug resistance in Plasmodium falciparum and Plasmodium vivax were considered in P. malariae. Earlier studies have reported polymorphisms in the P. malariae dhps, dhfr and kelch orthologues, which relate to sulfadoxine, pyrimethamine and artemisinin resistance, respectively.
Plasmodium malariae infection rate is generally underestimated and overlooked because it is often asymptomatic and mostly mixed with other Plasmodium species. Drug-resistant P. falciparum and P. vivax have been reported in many regions, including Thailand [9,10] and Myanmar [11,12]. Chloroquine-resistant P. falciparum was first documented in 1957 in Thailand which is now widely spread around the globe [13]. Polymorphisms of pfcrt and pfmdr1 genes have been linked to chloroquine resistance with the key determinant of K76T mutation in pfcrt gene [10]. Several studies of pfmdr1 supported that amplification and polymorphism were useful for prediction of chloroquine and mefloquine resistance [14][15][16]. The pvmdr1 Y976F mutation has been found to correlate with reduced susceptibility to chloroquine [17,18]. An increased copy number of pvmdr1 has been found in Thailand [17,18], suggesting that there is a relationship between pvmdr1 copy number and mefloquine pressure. In the Greater Mekong Sub-region, previous drug pressure from chloroquine and mefloquine on the P. malariae parasite population is expected to be considerable, since these drugs have been used widely for the treatment of vivax and falciparum malaria, respectively. The current study focuses on the P. malariae crt and mdr1 orthologue genes, which in P. falciparum are involved in chloroquine and mefloquine resistance. These genes were evaluated in P. malariae isolated of patients from Thailand and Myanmar between 2002 and 2016.

DNA extraction
Whole blood samples were collected from 95 P. malariae-infected patients from Thailand (51 samples) and Myanmar (44 samples) between 2002 and 2016 (Table 1). DNA was extracted using QIAamp DNA Mini Kit (Qiagen, Germany) and stored at −80 °C until use. Plasmodium species was confirmed by polymerase chain reaction based on 18 small-subunit ribosomal RNA [19,20].

Amplification of orthologous crt and mdr1 genes from Plasmodium malariae
Specific primers for amplification of pmcrt and pmmdr1 were designed based on the reference sequences (accession number LT594622.1 and LT594631.1). The primers and conditions used for amplification of pmcrt and pmmdr1 are listed in Additional files 1 and 2. Seminested and nested PCR were carried out in 20 μl, 2 mM of MgCl 2 , 250 µM of dNTPs, 250 nM of both forward and reverse primers, 0.5U BIOTAQ DNA polymerase (Bioline, UK), and 2.0 μl of DNA template. The approximated concentration of DNA templates were 80 to 150 ng/μl. The primary PCR products were used as a template in the secondary PCR. Cycling conditions were: initial denaturation at 94 °C for 5 min then denaturation at 94 °C for 1 min, annealing for 1 min, and extension at 72 °C for 1 min plus 30 s, using 30 cycles during the first round and 35 cycles during the second round, follow by a final extension of 72 °C for 10 min. The positive PCR products were purified and submitted for DNA sequencing in South Korea (Macrogen Inc., Korea).

Detection of pmmdr1 copy number variation
To assess pmmdr1 copy number variations, relative quantitative PCR was performed on the Applied Biosystems StepOnePlus ™ (Applied Biosystems, USA). The primers for pmmdr1 and pmβ-tubulin were pmmdr1F (5′-CAG ATG TGG GAA ACG ACA ATG-3′), pmmdr1R (5′-TAG AAG CTC CCT CCC CGT TT-3′), pmβ-tubulinF (5′-TGA AGC AAC TGG AGG AAG GT-3′), and pmβ-tubulinR (5′-GGA CCT GCT CGG ACA CTA TC-3′). Sso-Fast ™ EvaGreen ® Supermixes (Biorad, USA) was used as medium, and the thermal cycler profile was prepared according to the manufacturer's instruction. As calibrator, a single copy control was constructed from a plasmid by insertion of pmmdr1 (nucleotides [nt] 1102 to 1993) and pmβ-tubulin (nucleotides [nt] 1132 to 1215) fragments in a ratio of 1:1 into pGEM ® -T Easy Vector (Promega, USA). In addition, the two-copy pmmdr1 control was constructed by insertion of pmmdr1 and pmβtubulin fragments in a ratio of 2:1 into the pGEM ® -T Easy Vector (Promega, USA). Both of the control plasmids were confirmed their insertions by DNA sequencing. The pmβ-tubulin served as an internal control to calculate the relative amount of pmmdr1 gene by comparing the C t readings. Copy numbers were calculated as 2 −ΔΔCt . All reactions were performed in triplicate. Samples with a copy number more than 1.5 were classified as multiple copies. Reactions were repeated at least twice in case the samples contained ΔΔC t spread > 1.5, a C t value > 35, or a pmmdr1 copy number higher than 1.3.

Homology modelling of PmCRT and PmMDR1 proteins
The topology structure of PmCRT and PmMDR1 were predicted using Phyre2 [24] and SWISS-MODEL [25]. Homologous sequences of PmCRT and PmMDR1 proteins derived from different haplotypes were searched with PSI-Blast [26]. The secondary structure of PmCRT and PmMDR1 was predicted using Psi-pred 2.5 [27] and Disopred 2.4 softwares [28]. The membrane-spanning domains of PmCRT and PmMDR1 were predicted using Memsat_SVM program [29]. A multi-template approach was selected for prediction of the structural model of PmCRT and PmMDR1 proteins using SWISS-MODEL [25]. The model was evaluated with the VADAR tool and was visualized by Discovery Studio Visualizer V 17.2 [30].
For comparison, PmCRT haplotypes obtained from this study were aligned with the orthologous genes from human Plasmodium spp. The two-point mutations, H53P and E278D, found in P. malariae samples were not corresponded to the point mutations associated with chloroquine resistance described in P. falciparum. To predict the effect of these two-point mutations on PmCRT function, a topology model of PmCRT was constructed and then compared to PfCRT and PvCRT (Fig. 1). This showed that both H53P and E278D mutations are located outside the transmembrane domain (Fig. 1a), and are thus less likely to affect PmCRT function.    Table 4). For comparison, the non-synonymous mutations observed in pmmdr1 gene were aligned to sequences of the orthologous pfmdr1 and pvmdr1 genes. This showed that the observed point mutations in pmmdr1 were not corresponded to the pfmdr1 or pvmdr1 mutations associated with drug resistance (Additional file 3). The PmMDR1 sequences showed 16 haplotypes patterns (Table 5). Samples from Myanmar were classified into 13 haplotypes while the samples collected from Thailand were classified into 6 haplotypes. Fifty-five samples, 30 from Thailand and 25 from Myanmar, were wild type (haplotype 16). Haplotypes 3 (NCLLNTRTAA) and 4 (NYLLNSRTAA) were found in both countries. A total of 10 haplotypes (haplotypes 5-7 and 9-15) were only identified in Myanmar, whereas three other haplotypes (haplotypes 1, 2, and 8), were identified only in Thailand (Table 5).
A topology model of PmMDR1 was constructed and compared to the models of PfMDR1 and PvMDR1 (Fig. 2). PmMDR1 contains two nucleotide-binding domains (NBD1 and NBD2), facing the cytoplasm. Ten PmMDR1 substitution residues were identified: N6I, Y7C, L490I, L1063F, N1248I, T1266S, R1361S, T1406S, A1460S, and A1460T. Among the 10 mutations, N6I, Y7C were located on the cytoplasmic side, L490I located in nucleotide binding domain 1, L1063F located in transmembrane domain 11, N1248I, T1266S, R1361S, T1406S, A1460S, and A1460T located in NBD2. Although the mutations observed in P. malariae were not corresponding to residues associated with drug resistance reported previously in P. falciparum, the residues located in TMD11 and NBD might have an effect on PmMDR1 function. For this, structural models for each haplotype carrying mutations in TMD11 and NBD were constructed and analysed. The structural models representing these haplotypes showed similar characteristics and the mutated residues were predicted to have only a moderate effect on the PmMDR1 structure (Fig. 3).
Predicted changes in the physicochemical properties of each mutation were also addressed. For the L1063F mutant found in TMD11, leucine and phenylalanine share similar physicochemical properties, from which it can be inferred that this mutation has limited effect on the protein three-dimensional structure and function. A similar approach was followed for the four-point mutation identified in NBD, including L490I in NBD1 and T1266S, R1361S and A1460T/S in NBD2. To predict whether these point mutations affect the NBD structure and potentially its function, a structural model was constructed and residues involved in interactions between NBD and ATP were identified. For this, the highest conserved protein PGP1 from human was used as a template. Structural analysis of PmMDR1 model indicated that 20 and 19 residues are involved in ATP1 and ATP2 binding, respectively. Residue L490I was located outside of the ATP binding sites in NBD1. Residue T1266S, R1361S and A1460T/S in NBD2, were also located outside of the ATP binding sites.
A total of 95 P. malariae isolates were assessed for pmmdr1 gene amplification. The pmmdr1 copy number in all samples ranged from 0.75 to 1.25 (Fig. 4), which represent a single copy of the pmmdr1 gene.

Discussion
In this study, the anti-malarial drug target genes, pmcrt and pmmdr1, were characterized in P. malariae samples collected from Thailand and Myanmar and were compared to the orthologous mutations in pfcrt and pfmdr1 because polymorphisms in P. falciparum are well characterized for their association with chloroquine and mefloquine resistance [14][15][16]31]. Since chloroquine and mefloquine have been used widely in the treatment of vivax and falciparum malaria in Thailand and Myanmar, it can be assumed that the P. malariae parasite population in these countries has also been exposed to these drugs, given the frequent co-infection with other human Plasmodium species in patients with P. malariae infection.
Nine nucleotide polymorphisms were identified in 5 out of the 14 exons of the pmcrt gene, with 2 non-synonymous mutations found in exon 2 (H53P) and exon 8 (E278D), both were not corresponding to the orthologous mutations in P. falciparum involved in chloroquine      [32,33]. The in silico topology model of PmCRT revealed that H53P and E278D are located in cytoplasmic region and food vacuole, which are outside of the transmembrane domain, and are thus considered unlikely to affect PmCRT function. Genetic polymorphisms in the pmcrt gene in both Thailand and Myanmar samples were limited, similar to previous studies describing very limited polymorphisms in the pvcrt gene in P. vivax populations in Thailand [34].
The number of polymorphisms in this study was higher in the pmmdr1 gene. There were 13 SNPs, including 10  non-synonymous and 3 synonymous mutations, and none of these showed an equivalent position in pfmdr1 associated with anti-malarial drug resistance in P. falciparum. Some of the P. malariae studied here were mixed infections with other species, which might have an impact on polymorphisms of pmcrt and pmmdr1. There were 12 and 11 mixed infections found in the samples from Thailand and Myanmar, respectively. The patterns of nucleotide polymorphisms in pmcrt and pmmdr1 between single and mixed infections were compared. The proportion of mixed infected samples carrying mutations in Thailand and Myanmar was accounted for 14.29-22.72% (Additional file 4), suggesting that the mixed infections are unlikely to affect pmcrt and pmmdr1 mutations in this study.
The mutations in pmmdr1 codons L490I and L1063F correspond to the previously reported mutations in P. vivax at pvmdr1 codons L493L and F1076L [17,35]. In P. vivax, these polymorphisms have not been clearly linked to chloroquine resistance [36,37]. The observed polymorphisms in pmmdr1 were translated into a topology model of PmMDR1, showing that the mutations resulted in predicted protein changes located in the parasite cytoplasmic side of the protein, in nucleotide binding domains and one change in TMD11. In P. falciparum, residues in TMD11 were suggested as part of an anti-malarial binding pocket [38][39][40]. An in silico homology model of PmMDR1 was constructed and analysed to predict the effect of amino acid changes in TMD11 and NBD, which showed that the observed mutations were unlikely to affect the tertiary structure of the protein. Additionally, structural analysis of PmMDR1 polymorphisms found in NBD1 and NBD2 suggested that the mutations (L490I and L1063F) had only a moderate effect on the conformation of these domains, unlikely to affect the NBD function. Overall, the described polymorphisms were predicted to have insignificant impact on PmMDR1 protein morphology and function.
The samples used in this study were from Thailand and Myanmar collected at different periods of time, which may potentially affect the pattern of gene mutations [34]. The samples from Myanmar were collected in year 2009 in which artemether-lumefantrine was used as the first line treatment for uncomplicated malaria [41]. For Thailand, P. malariae were collected from two periods of time, year 2002-2008 and year 2012-2016. During those time periods, artesunatemefloquine was used as the first-line treatment before it was changed to dihydroartemisinin-piperaquine in 2015 [41]. The pmcrt in Thailand showed more mutations (haplotypes 1-3, 81.48-91.67%) when compared to Myanmar (31.82%) (Additional file 5). There was no difference in mutation pattern of pmmdr1 in both countries (Additional file 5). Although H53P and E278D mutations found in PmCRT were predicted that they were unlikely to have an impact on PmCRT function, the high proportion of those point mutations in Thailand might refer to geographical characteristics of the parasites. This will need to be confirmed in a larger sample size collected from different areas.
In P. falciparum, amplification of the pfmdr1 gene is strongly associated with mefloquine resistance [14][15][16]. Despite the presumed long-term mefloquine drug pressure on the P. malariae parasite population in the study areas, none of P. malariae samples carried amplification of pmmdr1. Possible explanations include firstly, that pmmdr1 amplification is not involved in mefloquine resistance in P. malariae and thus is not a good marker, and that there might be an alternative mechanism conferring mefloquine resistance in P. malariae other than mdr1 amplification. Secondly, the parasite loads during infection of P. malariae is rarely exceed 1000 parasites per µl of blood. Thus, the number of parasites under selective pressure would be low for P. malariae and the likelihood of selecting resistant is lower. However, some of the potential drug resistance markers that have been studied in P. malariae might be under selective pressure such as pmdhfr [42,43], pmdhps [44], and pmkelch [45]. Thirdly, although clinical P. malariae infection often presents as co-infection with other human Plasmodium species, there might be a large P. malariae reservoir outside of these patients with co-infection or a reservoir in nonhuman primates, so that overall P. malariae population has been limited to the selective pressure. The study of chloroquine-resistant P. vivax revealed that increased expression of pvcrt and pvmdr1 are associated with chloroquine resistance [46,47]. Moreover, the gene copy number of pvcrt was significantly higher in chloroquineresistant P. vivax [48]. In addition to this study, expression level of pmcrt and pmmdr1, and the copy number variation of pmcrt should be evaluated.

Conclusions
Polymorphisms in pmmdr1 were more frequently observed than in pmcrt. The non-synonymous mutations found in both pmcrt and pmmdr1 were unlikely to affect protein function. No amplification of pmmdr1 was observed in this study. If the orthologous resistance genes in P. malariae are indeed associated with anti-malarial drug resistance in this Plasmodium species, the findings suggest limited chloroquine and mefloquine drug pressure on the P. malariae populations in the study regions. Alternatively, anti-malarial drug resistance in P. malariae could differ from that described in P. falciparum and P. vivax, which will require further investigation.