Apoptotic like cell death in Leishmania is well established under various stress conditions, including oxidative stress [1-3]. Roles of various proteases, including caspases, are recognized for cell death mechanism in different organisms and human cell lines [4,5]. However, in Leishmania, studies remains focused on possible role of caspases and metacaspases in the apoptotic processes. Extensive studies about possible role of other proteases in the apoptotic processes of parasites are not reported. Thus, the mechanism of cell death with respect to role of proteases remains elusive. Caspases are the most important regulators of the apoptotic processes in higher eukaryotic organisms. Several groups have reported evidence of caspase like activities associated with apoptotic cell death in Leishmania [1-3,6]. However, caspase that has essential role in apoptosis in higher eukaryotic organism is absent in genome of Leishmania [7,8]. Metacaspases are cysteine proteases distinctly related to caspase are found in Leishmania (but absent in mammal), was initially thought to be responsible for caspase like activity. However, they are recently reported to have trypsin like activity rather than caspase like activity [9]. Metacaspases do not seem to have role in apoptosis mediated cell death of parasite [10]. Thus, it remains crucial to rigorously analyze possible enzymes that may be responsible for caspase like activity in the pathogen. It is possible that some other protease in Leishmania genome has evolved to perform additional function i.e., caspase like activity. Identification of enzyme cleaving caspase substrate may provide fundamental insights into apoptotic pathways in Leishmania. Furthermore, role of other proteases in the apoptotic cell death of Leishmania is also not extensively explored. To get an insight of the mechanism of apoptotic cell death and possible role of peptidases in the process, real time q-PCR analysis of various peptidase genes of L. donovani was done in miltefosine treated cells and compared with control cells (untreated cells) to know the change in mRNA expression level. It is worth mentioning that miltefosine is known to trigger apoptotic cell death in parasite [11].

Parasites, cell lines and chemicals

The Leishmania donovani (MHOM/IN/2010/BHU1081) promastigotes culture was obtained from Prof. Shyam Sundar, Banaras Hindu University. The apoptosis detection kit was procured from Calbiochem. Power syber green PCR master mix was obtained from Life Technologies. AMV first strand cDNA synthesis kit and DNase I was purchased from New England Biolabs. RNeasy Mini Kit was obtained from Qiagen. All the chemicals used in the experiments were of the highest grade procured from Sigma-Aldrich or Merck.

Sequence retrieval and classification of peptidases

We have analyzed the genome of Leishmania donovani available in GeneDB database (http://www.genedb.org/) where 8,021 proteins are found to be annotated [12]. Total 141 proteins (Supplementary Table 1) are identified/predicted as peptidases by BLAST in UniProtKB database (http://www.uniprot.org/help/uniprotkb) [13]; and classified in respective clan and family using MEROPS database (http://merops.sanger.ac.uk/) [14]. The BLAST analysis was performed against the target database UniProtKB to identify similar sequences of known functions. Proteins with low query coverage (< 40%) or low sequence identity (< 25%) were excluded from this study. Threshold E-value for BLAST was 10^-3. Structural homologues were identified using databases like PDB, UniProt. The identified proteases were classified using MEROPS databases in their respective clan and family. Some other family prediction tools namely SVMProt: Protein Functional Family Prediction tool (http://jing.cz3.nus.edu.sg/cgi-bin/svmprot.cgi) [15], ProtoNet tools (http://www.protonet.cs.huji.ac.il/) [16] and SUPERFAMILY [17] (http://supfam.cs.bris.ac.uk/) are used for the peptidases which are not enlisted in MEROPS.

Sub-cellular localization

UniProt database (http://www.uniprot.org/) and BRENDAEnzyme database (http://www.brenda-enzymes.org/) were used to check the experimental information about the sub-cellular localization of proteins. If experimental data of sub cellular localization was not available, various bioinformatics tools were used to predict the localization. Signal P 4.1 server (http://www.cbs.dtu.dk/services/SignalP/) has been used for predicting presence and location of signal peptide sequences. TMHMM tools (http://www.cbs.dtu.dk/services/TMHMM-2.0/) were used for prediction of transmembrane helices. All parameters for TMHMM tools were set as default. LocTree 3 (https://rostlab.org/services/loctree2/) and CELLO v 2.5 (http://cello.life.nctu.edu.tw/) were used for predicting the sub cellular localization based on two level support vector machine (SVM) system, sequence similarity and gene ontology information.

Apoptosis detection

Apoptotic assay was performed as reported earlier to confirm the apoptotic condition after miltefosine treatment [3]. In brief, L. donovani promastigote cells were treated with 50 μM of miltefosine for 24 h. After treatment, cells were harvested by centrifugation at 3000 rpm for 10 min and washed twice by cold PBS. Similarly, untreated cells were also centrifuged and washed by cold PBS. Both treated and untreated cells were suspended in 500 μl of 1X binding buffer, stained with annexin V-FITC and propidium iodide (PI) as per manufacturer’s instructions. After staining, the cells were analyzed using BD FACS Calibur flow cytometer and the fraction of cell population in different quadrants was analyzed by quadrant statistics using CellQuest Pro software.

Expression level analysis of various peptidase genes by real time PCR

Total RNA was extracted from untreated cells (control) and treated cells (50 μM miltefosine for 24 h) using RNeasy Mini Kit – QIAGEN according to the manufacturer’s protocol. In brief, the RNA was treated with DNase I-NEB for 30 min at 37ºC and heat inactivated at 70ºC for 5 min before cDNA preparation to remove contamination of genomic DNA. RNA was quantified by nanodrop at 260 nm absorbance and purity was accessed at A260/A280 ratio (>1.8). Equal amount of RNA was taken for first strand cDNA synthesis using AMV First Strand cDNA Synthesis Kit – NEB. Random hexamer primers were used for cDNA synthesis. Cycling parameters were 70ºC for 5 min, followed by 25ºC for 15 min and final extension at 42ºC for 45 min. The quality of cDNA was assessed by generating expression profiles of L. donovani housekeeping α-tubulin genes with real time PCR. Out of all proteases identified and classified (Figure 1), representative protease(s) from each clan was taken for expression analysis under induced apoptotic condition. Only one clan PB (threonine protease) which primarily includes proteosomes subunits is excluded. mRNA sequences for various L. donovani peptidases were taken from GeneDB (http://www.genedb.org) and primers were designed with primer 3–Biotools [18]. Accession number and primer sequences are shown in Table 1.

Figure 1: All Leishmania donovani peptidases are classified in their respective clan and family. Nomenclature of estimated peptidases is done on the basis of MEROPS database, SVMProt, ProtoNet and UPERFAMILY. Examples of each clan and family are also shown.

Table 1: List of primers used for gene expression analysis.

SyBR green PCR assays were performed on Applied Biosystems 7500 Real-Time PCR System. All quantitative assays were performed with housekeeping α-tubulin gene as endogenous control as reported in literature [19]. Cycling parameters were run at initialization at 50ºC for 2 min, 95ºC for 10 min to activate the DNA polymerase, 40 cycles of denaturation at 95ºC for 30 s, annealing at 60ºC for 30 s and then extension at 72ºC for 30 s. Melt curve was included and analyzed to access the specificity of PCR product. All results were normalized to the expression of α-tubulin gene in control and miltefosine treated cells. In addition to this, the PCR product was run on 1% agarose gel to check the amplicon size and nonspecific amplifications. Data were analyzed using Applied Biosystems SDS v2.0.6 software.

Statistical analysis

All statistical analysis between two groups were performed by Student’s unpaired t-test in SigmaPlot software. Differences between two groups with p-value of less than 0.05 were considered as statistically significant. All results were expressed as mean + SD of at least three independent experiments.

Classification of peptidases

Total 141 proteins are annotated as peptidases in L. donovani. The number of identified proteases are close to the number of proteases reported from another species i.e., L. major (154 peptidases) [20]. Majorly five classes of peptidases are present in parasite namely aspartic peptidase, cysteine peptidase, metallopeptidase, serine and threonine peptidase. We have found two families of aspartic peptidases, ten families of cysteine peptidases, twenty one families of metallopeptidases, seven families of serine peptidases and one family of threonine peptidases in L. donovani. Cysteine peptidase is the largest clan containing 60 peptidases followed by metallopeptidases-52, serine peptidases-15, threonine peptidases-13 and aspartic peptidases containing only two peptidases as shown in Figure 1. All peptidases were classified with high confidence and are reported in Supplementary Table 1.

Sub-cellular localization

Function of proteins is mainly related to its sub-cellular localization. It is important to predict the sub-cellular localization to get the insight of protein. For membrane topology prediction, TMHMM server is run which is based on hidden Markov model (HMM). Statistical reports suggest that accuracy of TMHMM is 84% for prediction of transmembrane helices but it drops when signal peptide sequences are present. Moreover, TMHMM can differentiate between soluble and membrane protein with sensitivity and specificity greater than 99% [21]. To check whether presence of signal peptides in amino acid sequences, SignalP 4.1 server are used which predicts signal peptide/non signal peptide using a combination of several artificial neural networks [22]. Out of total 141 predicted L. donovani peptidases, 24 peptidases having trans-membrane helices and 8 peptidases have signal peptide sequences which counts for 23% of total peptidases. Results of TMHMM and Signal P prediction are shown in Figure 2 and Supplementary Table 1. If experimental data of sub-cellular localization are not available, LocTree3 server and CELLO were used to predict the localization. In eukaryotes, LocTree3 can predict 18 classes of sub-cellular localization with overall accuracy of 80% whereas; its accuracy for extracellular protein and nuclear protein is 88% and 81%, respectively. LocTree3 mainly uses homology based localization annotation using PSI-BLAST, UniProt, PDB and SWISS-PROT databases. In addition to homology based annotation, it also utilizes support vector machine system. The result displayed by LocTree3 prediction shows a score ranging between 0-100 in which 100 is the most reliable prediction. Moreover, the result also shows the expected accuracy in percentage, predicted single localization class, gene ontology terms and type of annotation used [23]. In our studies, the threshold set for predicted accuracy was >80% and prediction score was >80 on 0-100 scale. Based on the result of LocTree3 prediction, the maximum percentage of peptidases present in cytoplasm is 60% followed by nuclear fraction 19%, secreted 15%, mitochondria 8%, mitochondrial membrane 4%, endoplasmic reticulum membrane 6%, endoplasmic reticulum 4%, vacuole and golgi apparatus membrane fraction accounts for 1% of peptidases. The result of LocTree3 prediction is shown in Figure 2B. When the sub-cellular localization of peptidases was checked in experimental database i.e., BRENDA enzyme database, we observed that many of the peptidases are localized at multiple sub-cellular locations. LocTree3 does not predict the multiple sub-cellular localization of signal peptidase; hence it became important to use another sub-cellular localization prediction tool which may predict multiple localization sites with better accuracy. CELLO v 2.5: Subcellular localization predictor server developed by Molecular Bioinformatics Centre, NCTU; uses two level support vector machine systems to predict the subcellular localization in which the first level contain the support vector machine classifiers whereas second level comprises of probability distribution for possible localizations. Overall accuracy for CELLO prediction is 94.9% if sequence identity is >30%. In case of sequence identity <30%, most of the subcellular localization predictor tools shows fairly bad result with low accuracy, whereas CELLO predicts the overall accuracy of 82.6% [24]. The result for CELLO prediction is shown in Figure 2C.

Figure 2: Subcellular localization of 141 peptidases of Leishmania donovani utilizing the prediction result of various bioinformatics tools. The chart shows the (A) transmembrane helices and secretory peptides predicted by TMHMM and Signal P server. (B) subcellular localization prediction by LocTree3 server. (C) subcellular localization prediction by CELLOv2.0 server.

Apoptosis detection

L. donovani promastigotes cells treated with 50 μM of miltefosine for 24 h were stained with annexin V-FITC and PI, and data were analysed by flow cytometry. As also reported in literature [11], the condition provided us significant fraction of parasite in apoptotic condition. The parasite after treatment with miltefosine as mentioned earlier was used for analysis of protease genes in apoptotic conditions.

Expression analysis of L. donovani protease genes by Real time-qPCR

In an initial attempt to identify the involvement of various peptidases in apoptotic cell death (Programmed Cell Death) pathway of parasite, the mRNA level of 24 protease genes of L. donovani were analyzed by Real time-qPCR in control cells and miltefosine treated apoptotic cells. Alpha tubulin was used as an endogenous control. Minimum one peptidase was preferred from each clan and primer was designed for respective peptidase gene of L. donovani to compare the mRNA level in control and treated cells. Our gene expression analysis of Leishmania peptidases in apoptotic condition shows altered expression of several proteases hinting involvement of these peptidases in the process, directly or indirectly. The results of mRNA expression analysis are shown in Figure 3. The previous report has shown that the activity of metacaspases is increased in H₂O₂ induced apoptosis of parasite [9]. In our studies, we have also reported that the mRNA expression level of metacaspases is increased by two fold in treated cells. Moreover, the mRNA expression level of other cysteine peptidases involved in autophagy processes of parasite is increased more than metacaspases. Two cysteine peptidases namely ATG4.1 and ATG4.2 are reported in the autophagic processes of parasite in which ATG4.2 is more important than ATG4.1 [25]. Multiple orthologues of single ATG4 are reported in case of mammals. Orthologues of ATG4 in mammals, mainly ATG4D is involved in apoptotic processes and mitophagy [26,27]. As increased mRNA expression level of ATG4.1 and ATG4.2 in treated cells suggest that there may be the possible role of ATG4 in apoptotic processes of parasite. In addition to this, expression level of several other peptidases genes are increased many folds indicating that the involvement of peptidases directly or indirectly in apoptosis of parasite. Many of the peptidases, whose mRNA expression level is found to be increased in our data, are predicted by bioinformatics tools to be localized at common sub cellular compartment. There is very high possibility that some of the peptidases which may have evolved in such a way to help in PCD processes of parasite, may act as an caspase like activity or act as an effector molecule to activate other peptidases for PCD. As very few experimental data available and limitations of bioinformatics tools, it is difficult to predict the cleavage site position (if present in peptidases) with significant accuracy.

Figure 3: Expression level of different protease genes of Leishmania donovani in apoptotic conditions. Cells were treated with 50 μM of miltefosine for 24 h. Equal amout of cDNA was taken for real time-qPCR anlaysis. Αlpha-tubulin was used as an endogenous control. Peptidase genes belong to clan of (A) Cysteine peptidase (B) Aspartic peptidase and Serine peptidase (C) Metallopeptidase; Results are mean ± SD of three independent experiments. Graph was analyzed using Student’s unpaired t-test in SigmaPlot (*denotes p value ≤ 0.05 and **denotes p value < 0.01). (CCP: Calpain-like cysteine peptidase; CPB: Cysteine peptidase B; CPA: Cysteine peptidase A; CPC: Cysteine peptidase C; UbH: Ubiquitin hydrolase; MCas: Metacspase; ATG4.1-AUT2/APG4/ATG4 cysteine peptidase; ATG4.2 AUT2/ APG4/ATG4 cysteine peptidase; PGP: Pyroglutamyl-peptidase I; PrAP: Presenilin like aminopeptidase; SPP: Signal peptide peptidase; SerP: Serine peptidase; SP1: Signal peptidase type I; PAP: Puromycin sensitive aminopeptidase; ZnMP1: Mitochondrial ATP-dependent zinc metallopeptidase; ZnMP2: ATP-dependent zinc metallopeptidase, putative; ZnCP: Zinc carboxypeptidase; PMP: Pitrilysin-like metalloprotease; MLP: Metalloprotease like protein; APP: Aminopeptidase P; MAP1: Methionine aminopeptidase; MAP2: Methionine aminopeptidase 2; PepT: Peptidase t; AAP: Aspartylaminopeptidase).

Prediction of BH3-like domain and caspase 3 cleavage site

BH3-like domain is known as pro-apoptotic protein which binds to the hydrophobic pocket of multi domain Bcl-2 family, a well-known anti-apoptotic protein and triggers apoptosis via caspase mediated pathway [28,29]. As BH3–like domain having low conservation and shorter peptide length, it is difficult to annotate this domain by common bioinformatics tools. Several studies have reported that human ATG4D (hATG4D) contains Caspase 3 cleavage site on its N-terminal and BH3-domain like protein on its C-terminal whereas BH3–domain like protein contains conserved LXXXXD region [28,30]. The amino acid sequences of ATG4.1 protease of L. donovani was physically verified and run the ESPript 3.o with already known BH3-domain like protein whose results are shown in Figure 4. All protein sequences were retrieved from UniProt database and accession numbers are written in the legend of Figure 4. The highly conserved region LXXXXD was observed in L. donovani ATG4.1 protease using ESPript 3.o alignment. Based on this result, we have proposed the presence of BH3-domain like protein in L. donovani ATG4.1 protease. Moreover, Caspase 3 cleavage site DEVD*T are also identified at C-terminus of L. donovani ATG4.1 protease and shown in Figure 4.

Figure 4: Sequence alignment of Leishmania donovani ATG4.1 protease with known BH3-like domain protein. C-terminal sequence positions are indicated on right. *in red shows highly conserved LXXXXD region. Known BH3–like protein and their UniProt IDs are: BID-55957; TGM2-P21980; NOXA-Q13794; RAD9-Q99638; BIK-Q13323; HRK-O00198; Beclin1-Q14457; ApoL1-O14791. Caspase 3 cleavage sequence of L. donovani ATG4.1 protease (GeneDB ID: LdBPK_300270.1) C-terminal sequence position is also shown.

Mainly five types of peptidases are found in Leishmania–Aspartic, Cysteine, Serine, Threonine and Metallo peptidases. In L. donovani, two aspartic peptidases are present and their predicted sub-cellular localization is in plasma membrane/endoplasmic reticulum membrane. Sequence homology reveals that these proteases are similar to presenilin like aspartic peptidase (PrAP) and signal peptide peptidase. In humans, PrAP are abundantly present in endoplasmic reticulum and helps in the processing of amyloid precursor protein (APP) [31]. Moreover, PrAP are also reported to have role in macroautophagy as PrAP gene knockout mice do not show macroautophagic process [32]. In our mRNA expression analysis, no significant increase in transcription level of aspartic peptidases was found in apoptotic conditions of parasite. As role of PrAP was reported in macroautophagy, hence it is important to establish the role of L. donovani PrAP in autophagic processes and the molecular mechanism behind it.

Classification of peptidases of L. donovani revealed that 61 cysteine peptidases are present which mainly belongs to clan CA. This clan mainly contains cysteine peptidase A, B, C, cathepsins, calpains, ubiquitin hydrolases and ATG4 peptidases. In addition to this, two more clan of cysteine peptidases are present- Clan CD and CF. Metacaspases, an important peptidase which were initially considered to have caspase like activities in parasite, are important member of clan CD and pyroglutamyl peptidases are member of clan CF. In higher organisms, calpain helps in Ca²⁺ regulated signaling pathway, cell differentiation and apoptotic mediated cell death. It has been reported that Ca²⁺ activates the calpain which in turn cleaves the Bcl2- family protein i.e., anti-apoptotic protein eventually increasing the apoptotic cascade [33]. Moreover, activation of calpain is also required in macroautophagic process. Activated calpain cleaves ATG5 protein and in turn truncated ATG5 induces the release of cytochrome c from mitochondria leading to apoptotic cell death [34]. We have observed that in apoptotic condition, transcription of calpain gene has increased more than three folds compared to control. ATG5 gene (LmjF.30.0980) is present in L. major but till now it is not annotated in L. donovani genome. Annotation of ATG5 gene in L. donovani and experimental validation is crucial to understand the increased transcription level of calpain in autophagy and apoptotic conditions. In addition to this, we have observed the increased transcription of ubiquitin hydrolase gene suggesting higher protein turn over in apoptosis.

Interplay of autophagy and apoptosis has been reported in the literature [35]. However, most of these studies are reported in higher organism. Existence of autophagy and apoptosis modes of cell death in Leishmania are already established [36]. Reports suggesting that protozoan parasites evade host cell defense system using autophagy and has important role in infection process [36-39]. However, not much detailed investigation has been done about correlation between apoptosis and autophagy. Over-expression of ATG4.1 and ATG4.2 under apoptotic condition of Leishmania parasite points out toward correlation in there two cell death mechanism in the protozoan parasite. ATG4 proteins, cysteine proteases, reported to have important function in autophagy process as they are involved in formation of autophagosomes and their subsequent targeting to lysosomes [40]. Very significant over-expression of ATG4.1 and ATG4.2 under apoptotic condition of Leishmania parasite suggests some correlation between these two processes in protozoan parasites as well. The link between autophagy and apoptosis is poorly understood and remains matter of controversy [41,42]. Even few proteins are involved in both the processes and the cross talk between autophagy and apoptosis needs further extensive investigation.

Two ATG4 cysteine proteases namely ATG4.1 and ATG4.2 of clan CA, family C54 are present in L. donovani and their role in macroautophagy is well documented in L. major [36]. In case of mammals, four orthologs of ATG4 (ATG4A-D) proteases are present, having various roles in autophagosome formation [35]. Human ATG4D (hATG4D) is reported to have distinct role in mitophagy and apoptosis [29]. It contains Caspase 3 cleavage site DEVD*K at position 63 from N-terminus [29]. After cleavage by Caspase 3, hATG4D localizes to mitochondria where it increases the permeabilization of outer mitochondrial membrane which leads to release of pro apoptotic factor from mitochondria and eventually induction of apoptosis [29,43]. BH3-like domain present at the C-terminus of hATG4D which are exposed after proteolysis by Caspase 3. These are known to bind with Bcl2-family, an anti-apoptotic protein. Binding of BH3- like domains to Bcl2-family inhibits their binding with pro apoptotic proteins and thus induces apoptosis [44]. Moreover, during oxidative stress, ATG4D also involves in mitophagy (selective removal of damaged mitochondria) and decreases the release of pro apoptotic factors eventually limiting the apoptosis [26,27,29]. Interestingly, we have observed that L. donovani ATG4.1 protease also contains Caspase 3 cleavage site DEVD*T at position 386 of C-terminus but till now no Caspase 3 gene is identified in parasite. So it remains elusive to know the importance of Caspase 3 cleavage site in ATG4.1 proteases of parasite. We have also reported the presence of BH3-like domain at N-terminus of ATG4.1 protease. Several reports have suggested that after induction of apoptotic condition, cell lysate of parasite is able to cleave caspase3/7 substrates [2,3]. So, there may be possibility that some of the proteases which are cleaving caspase 3 substrates; may cleave ATG4.1 at DEVD*T site and after cleavage it may expose proposed BH3-like domain to mitochondria and increases the apoptotic demise. Experimental studies are important to validate the significance of BH3-like domain present at the N-terminus of ATG4.1 protease. In our mRNA expression data, we have found increased transcription of ATG4 gene in apoptotic conditions suggesting its role in generation of apoptotic proteins from mitochondria or increasing permeability of outer mitochondrial membrane and escalating PCD.

Eight clans of metallo-peptidases are identified in L. donovani and mainly consist of aminopeptidases, leishmanolysin, carboxypeptidases and dipeptidases. We have reported the increased transcription level of AAP, PepT, MLP, ZnMP2, PAP in apoptotic conditions by >2 times and in MAP2, ZnCP by >3.5 times. Some of the metallo-peptidases like PMP, APP and MAP1 do not show any significant altered expression level in our experiment. Very few of the metallo-peptidases are characterized in Leishmania. So, detailed experimental analysis is required to understand the role of metallo-peptidases in PCD of parasite.

Total 15 Serine peptidases are identified in L. donovani which, belongs to clan SB, SC and SF. Several studies have reported that Serine peptidases are involved in host cell invasion [45,46]. In plasmodium, serine peptidases are involved in proteolytic cleavage of proteins present in cytopasmic membrane of RBC thereby allowing the parasite to infect the host [47]. In trypanosomes, oligopeptidase B, a serine protease belonging to clan SC helps in the invasion of parasite to mammalian host [36]. Real time-qPCR data of our experiment reveals that no significant altered expression level of SP1 gene in apoptotic condition while there is increased transcription level in SerP gene of L. donovani. Various studies have suggested the role of non caspase proteases which includes metacaspases, calpains, cathepsins etc in PCD of parasite. Based on that, it is essential to recognize the non caspase proteases and establish the molecular mechanism of apoptosis in parasite.

We have earlier reported studies on redox metabolism of the parasite [2,3]. We have focused our research on understanding involvement of various Leishmania protease(s) in apoptotic process. The data showed altered expression of various protease mRNA in apoptotic condition of Leishmania. This indicates possible role of various proteases in apoptotic process, directly or indirectly. Further, the data suggests some crosstalk between autophagy and apoptosis modes of death of Leishmania parasite.

Research fellowships to RK by IIT Guwahati are acknowledged. Financial support by Department of Biotechnology, Government of India in the form of research grants (Project no: BT/01/IYBA/2009) to VKD is also acknowledged.

VKD designed the research plan, R.K. and P.M performed the experiments, R.K. and P.M analyzed the data, VKD supervised the research, and R.K. and P.M wrote the manuscript.

No potential conflict of interest was reported by the authors.