# Blood transcriptional profiles distinguish different clinical stages of cutaneous leishmaniasis in humans

Fariborz Bahrami, Nasrin Masoudzadeh, Suzanne van Veen, Josefine Persson, Arezou Lari, Hamzeh Sarvnaz, Yasaman Taslimi, Malin Östensson, Björn Andersson, Iraj Sharifi, et al.

► To cite this version:

Fariborz Bahrami, Nasrin Masoudzadeh, Suzanne van Veen, Josefine Persson, Arezou Lari, et al. Blood transcriptional profiles distinguish different clinical stages of cutaneous leishmaniasis in humans. Molecular Immunology, 2022, 149, pp.165-173. 10.1016/j.molimm.2022.07.008. hal-04263733

HAL Id: hal-04263733

https://hal.science/hal-04263733v1

Submitted on 29 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Blood transcriptional profiles distinguish different clinical stages of cutaneous leishmaniasis in humans

Fariborz Bahrami<sup>a</sup>, Nasrin Masoudzadeh<sup>b</sup>, Suzanne Van Veen<sup>c</sup>, Josefine Persson<sup>d</sup>, Arezou Lari<sup>e</sup>, Hamzeh Sarvnaz<sup>b</sup>, Yasaman Taslimi<sup>b</sup>, Malin Östensson<sup>f</sup>, Björn Andersson<sup>f</sup>, Iraj Sharifi<sup>g</sup>, Vahid Mashayekhi Goyonlo<sup>h</sup>, Tom HM Ottenhoff<sup>c</sup>, Mariëlle C. Haks<sup>c</sup>, Ali M. Harandi<sup>d,i,1,*</sup>, Sima Rafati<sup>b,1,**</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

Cutaneous leishmaniasis

Leishmania tropica

dcRT-MLPA

IFI16

CCL11

## ABSTRACT

Cutaneous leishmaniasis (CL) is a neglected tropical disease with severe morbidity and socioeconomic sequelae. A better understanding of underlying immune mechanisms that lead to different clinical outcomes of CL could inform the rational design of intervention measures. While transcriptomic analyses of CL lesions were recently reported by us and others, there is a dearth of information on the expression of immune-related genes in the blood of CL patients. Herein, we investigated immune-related gene expression in whole blood samples collected from individuals with different clinical stages of CL along with healthy volunteers in an endemic CL region where *Leishmania* (*L.*) *tropica* is prevalent. Study participants were categorized into asymptomatic (LST+) and healthy uninfected (LST-) groups based on their leishmanin skin test (LST). Whole blood PAXgene samples were collected from volunteers, who had healed CL lesions, and patients with active *L. tropica* cutaneous lesions. Quality RNA extracted from 57 blood samples were subjected to Dual-color reverse-transcription multiplex ligation-dependent probe amplification (dcRT-MLPA) assay for profiling 144 immune-related genes. Results show significant changes in the expression of genes involved in interferon signaling pathway in the blood of active CL patients, asymptomatics and healed individuals. Nonetheless, distinct profiles for several immune-related genes were identified in the healed, the asymptomatic, and the CL patients compared to the healthy controls. Among others, *IFI16* and *CCL11* were found as immune transcript signatures for the healed and the asymptomatic individuals, respectively. These results warrant further exploration to pinpoint novel blood biomarkers for different clinical stages of CL.

## 1. Introduction

Cutaneous leishmaniasis (CL) caused by *Leishmania* (*L.*) species is a neglected tropical disease in many parts of the developing world (Alvar

et al., 2012). Although CL is considered a non-fatal localized skin infection, it has a major impact on public health systems of the affected nations, causing social and economic problems due to disfigurement, morbidity, loss of productivity and medical costs for the public sector.

<footnote><sup>a</sup> Department of Immunology, Pasteur Institute of Iran, Tehran, Iran</footnote>

<footnote><sup>b</sup> Department of Immunotherapy and Leishmania Vaccine Research, Pasteur Institute of Iran, Tehran, Iran</footnote>

<footnote><sup>c</sup> Department of Infectious Diseases, Leiden University Medical Center, Leiden, the Netherlands</footnote>

<footnote><sup>d</sup> Department of Microbiology and Immunology, Institute of Biomedicine, Sahlgrenska Academy, University of Gothenburg, Gothenburg, Sweden</footnote>

<footnote><sup>e</sup> Systems Biomedicine Unit, Pasteur Institute of Iran, Tehran, Iran</footnote>

<footnote><sup>f</sup> Bioinformatics Core Facility, University of Gothenburg, Gothenburg, Sweden</footnote>

<footnote><sup>g</sup> Leishmaniasis Research Center, Kerman University of Medical Sciences, Kerman, Iran</footnote>

<footnote><sup>h</sup> Cutaneous Leishmaniasis Research Center, Mashhad University of Medical Sciences, Mashhad, Iran</footnote>

<footnote><sup>i</sup> Vaccine Evaluation Center, BC Children's Hospital Research Institute, University of British Columbia, Vancouver, BC, Canada</footnote>

<footnote>* Corresponding author at: Department of Microbiology and Immunology, Institute of Biomedicine, Sahlgrenska Academy, University of Gothenburg, Gothenburg, Sweden.</footnote>

<footnote>** Corresponding author.</footnote>

<footnote>E-mail addresses: ali.harandi@microbio.gu.se (A.M. Harandi), s_rafati@yahoo.com (S. Rafati).</footnote>

<footnote n="1">¹ These authors contributed equally to the work</footnote>
The complex and interactive factors involved in the transmission cycle of the protozoan parasite in the human host, the vectors and the reservoirs, along with social and environmental factors, have imposed significant challenges to control the disease (Banuls et al., 2011). The addition factor to this complexity is the dual role of interferon (IFN-γ) which determines the pathologic outcome of the infection. Activation of phagocytic cells such as macrophages due to production of IFN-γ by CD4+ Th1 cells was shown to limit the parasite replication (Scott and Novais, 2016). Conversely, IFN-γ is able to activate the downstream proinflammatory mechanisms by inflammasomes, leading to cytolysis and exacerbated pathology. It is reported that IFN-induced GTPases named as guanylate binding proteins (GBPs) are able to direct and accelerate the inflammasome machinery (Fisch et al., 2019; Kim et al., 2016). Moreover, treatment failure in *L. braziliensis* infected patients was found to be associated with a high lesional expression of IFN-γ and activation of cytolytic pathways during the course of disease (Amorim et al., 2019).

Since no effective human vaccine against CL exists, it is crucial to understand better factors that may affect the disease's outcome, including the pathogen's impact on the host immune system. In this regard, comprehensive and holistic approaches, such as transcriptomic analyses in different human populations within endemic regions may be helpful to muster an effective strategy to counter CL. Recently, transcriptional profiles, molecular pathways, and gene networks of human responses during distinct states of visceral leishmaniasis (VL) caused by *L. infantum* have been reported (Gardinassi et al., 2016). Our group has recently pinpointed the molecular signatures of *L. tropica*-infected human skin lesions with different clinical manifestations. This report has identified unique molecular pathways and inflammatory proteins in the skin lesion of CL-infected patients (Masoudzadeh et al., 2020b). We have recently reviewed the state-of-the-art on transcriptomic profiling of CL patients (Masoudzadeh et al., 2020a). Similar approaches to identify blood signatures of populations at different clinical stages of the CL infection in the endemic areas may inform the epidemiological studies and the development of preventive measures against CL. A recent transcriptomic study on the peripheral blood of CL patients due to *L. braziliensis* infection has revealed transcriptional signatures such as IFN-stimulated genes (Amorim et al., 2021).

Among people living in CL endemic areas, there is a large group of asymptomatically-infected individuals who, despite their exposure to *Leishmania* species (as detectable by their positive delayed-type hypersensitivity (DTH) reactions to leishmanial antigens used in leishmanin skin test (LST), have no history of CL lesion. Both individuals with healed CL lesions and the asymptomatic infected residents in the endemic areas are protected against the development of CL lesions during their lifetime (Banuls et al., 2011; Bahrami et al., 2018). Therefore, further analyses of the immune responses at the transcriptomic level of these two groups along with the patients with active CL lesions may inform rational design of diagnostics and intervention measures against CL. Herein, we analyzed whole blood samples retrieved from individuals of the groups above in two *L. tropica* endemic regions in Iran by dual-color reverse-transcription multiplex ligation-dependent probe amplification (dcRT-MLPA). Results showed that despite causing local damage to the skin, CL leads to significant alterations in the expression of immune-and inflammation-related genes in the blood, depending on the clinical status of the disease.

## 2. Materials and methods

### 2.1. Ethics statement

The study's protocols, questionnaire, and consent forms were approved in advance by the Ethics Committee of Biomedical Researches at Pasteur Institute of Iran (authorization # IR.PII.REC.1396.31). The subjects for this study were chosen in accordance with the following process. The healthy volunteers were the residents of CL endemic

neighborhoods whom were invited by the local health personnel to assemble in local clinics or were the staff of the affiliated universities in Kerman and Mashhad. The healed volunteers were invited after consulting the medical records of patients who had one or two CL scars during the previous CL season (September-December 2017) at a local clinic. Patients with active cutaneous lesions referred to a special CL clinic were invited to participate. All volunteers were informed in plain language, orally as well as in writing about the main aims of the project, LST, and the blood sampling. Signed written informed consent forms were obtained from all volunteers. The exclusion criteria from the study were: pregnancy, lactating, major illnesses such as immunosuppression, HIV infection, cancer, renal failure, heart diseases, tuberculosis, jaundice, diabetes mellitus, allergies, and recent significant surgeries.

### 2.2. Study design, demographics, and blood sampling

Whole blood samples (2.5 ml) were collected in PAXgene tubes (Becton Dickinson, USA) according to the manufacturer's instructions from two *L. tropica* endemic regions of CL, namely the cities of Kerman and Mashhad in Iran. The residents of the endemic areas of Kerman (capital of Kerman province, ~1000 km SE of Tehran) were invited, through the Leishmaniasis Research Center of Kerman University of Medical Sciences, to participate in the study. The study's participants are described below. Twenty-eight volunteers who appeared healthy and did not meet the exclusion criteria were enrolled in the study. Following the blood collection, these volunteers were subjected to LST, and 48 h later, the LST+ and the LST- groups were assigned as the asymptomatic and the healthy groups, respectively. Moreover, fourteen healed volunteers from Kerman and eleven patients with active cutaneous lesions were referred to a special CL clinic at Mashhad Medical School along with four healthy LST- residents from *L. tropica* endemic regions of Mashhad (capital of Khorasan Razavi province in Iran) were enrolled in this study. All patients had clinically-verified cutaneous lesions due to *L. tropica* for less than one year without prior leishmanial treatment and a confirmed positive result by diagnostic *L. tropica*-specific PCR as described below. The filled PAXgene tubes were transferred to Pasteur Institute of Iran in Tehran and kept at -80°C till RNA extraction. The demographics of the study participants are shown in Table 1.

### 2.3. PCR diagnostics of patients

As described previously, *L. tropica* patients were diagnosed by an ITS1 PCR-RFLP assay (Masoudzadeh et al., 2020b). In short, samples were obtained using a non-invasive tape disc-based sampling method, and genomic DNA was used for a PCR, targeting ITS1 region of *Leishmania*, followed by RFLP for species-specific identification, as described elsewhere (Taslimi et al., 2017).

### 2.4. LST

According to the procedures described previously, the TDR/WHO reference leishmanin produced under good manufacturing practices

<table><caption>**Table 1**
Demographics of the studied groups.</caption><thead><tr><th>Characteristics</th><th>Healed</th><th>Asymptomatic</th><th>Patients</th><th>Healthy</th></tr></thead><tbody><tr><td>Females<br>participated<br>(mean age ± SD)</td><td>6<br>(31.7 ± 5.6)</td><td>9<br>(30.9 ± 4.7)</td><td>7<br>(50.3 ± 9.1)</td><td>4<br>(35.5 ± 9.8)</td></tr><tr><td>Males participated<br>(mean age ± SD<br>age)</td><td>8<br>(30.6 ± 6.2)</td><td>8<br>(32.6 ± 7.2)</td><td>4<br>(45.8 ± 23.0)</td><td>11<br>(31.1 ± 10.4)</td></tr><tr><td><strong>Total participants</strong><br>(mean age ± SD)</td><td><strong>14</strong><br>(31.1 ± 6.0)</td><td><strong>17</strong><br>(31.7 ± 6.0)</td><td><strong>11</strong><br>(48.6 ± 15.8)</td><td><strong>15</strong><br>(32.3 ± 10.4)</td></tr></tbody></table>
guidelines at the Department of Immunology of Pasteur Institute of Iran was used for LST (Alimohammadian et al., 1993). Briefly, 0.1 ml of leishmanin was injected intradermally into the inner surface of the left forearm. The diameter of induration at the injection site was measured by a millimeter graduated ruler 48 h later using the ball-point pen technique (Sokal, 1975). Indurations more than or equal to 5 mm were considered as positive LST.

## 2.5. RNA extraction

RNA from each 2.5-ml PAXgene blood sample was extracted by PAXgene Blood RNA kit (PreAnalytiX, UK), according to the manufacturer's instruction, and kept at -80°C until use. The average RNA integrity number (RIN) of the extracted RNA was assayed by an Agilent 2200 TapeStation (Agilent Technologies, USA). The average RIN ± SD of the total RNA samples obtained from PAXgene tubes was 7.9 ± 0.5.

## 2.6. dcRT-MLPA

Targeted gene expression profiling was performed as described in detail previously (Joosten et al., 2012; Haks et al., 2015). Briefly, a specific reverse transcription (RT) primer was designed for each target-specific sequence, located immediately downstream of the left- and right-hand half-probe target sequence. RNA was reverse-transcribed, and the left- and right-hand half-probes were hybridized to the cDNA at 60 °C overnight. Annealed half-probes were ligated by ligase and were subsequently amplified by PCR (33 cycles of 30 s/95 °C, 30 s/58 °C, 60 s/72 °C, followed by 1 cycle of 20 min/72 °C). Primers and probes were from Sigma-Aldrich Chemie (Zwijndrecht, The Netherlands) and MLPA reagents from MRC-Holland (Amsterdam, The Netherlands). Reverse-transcription primers and half-probes were designed by Leiden University Medical Center (Department of Infectious Diseases, Leiden, The Netherlands) (Joosten et al., 2012) and comprised sequences for 4 housekeeping genes along with 144 selected genes to profile innate, adaptive, and inflammatory immune responses (Supplementary Table S1). Of the 4 housekeeping genes included, Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH) was most stably expressed and used as reference gene in this study. The raw and analyzed data were submitted to the Gene Expression Omnibus (GEO) repository (accession number GSE197222).

## 2.7. Statistical and bioinformatical analyses

Trace data were analyzed using GeneMapper software 5 package (Applied Biosystems, Bleiswijk, The Netherlands). The areas of each assigned peak (in arbitrary units) were exported for further analysis in Microsoft Excel. Data were normalized to GAPDH and signals below the threshold value for noise cutoff in GeneMapper (Log2-transformed peak area of 7.64) were assigned the threshold value. The normalized data were Log2-transformed for statistical analysis. Differential expression analysis was performed using R (version 3.5.2; http://CRAN.R-project.org/bin/windows/base) statistical computing environment and Bio-Conductor (version 3.9). The Limma package with the empirical Bayes method was used to assess the differentially-expressed genes (DEGs) between the asymptomatic, the healed, and the CL patient groups versus the healthy control group. The *P*-value of less than 0.05 was considered significant, and all *P* values were two-tailed. Heatmap graphical presentation was performed using GraphPad Prism (version 8.0.1; GraphPad Software, La Jolla, CA, USA). Volcano plots were generated in R 4.1.1 using ggplot2 (version 3.3.5) and ggrepel (version 0.9.1) packages. Biological pathways plots were created using ggplot2 (version 3.3.5) package. Principal component analysis (PCA) plot was created using tidyverse package version 1.3.1.9 and ggforce package version 0.3.3. The protein-protein interaction network was analyzed and visualized by STRING (https://string-db.org/). The network data from the STRING database displayed a combination of data, including experiments,

databases, co-expression, neighborhood, gene fusion, and co-occurrence (excluding text mining). The minimum required interaction score was medium confidence (0.400) (Scott and Novais, 2016).

# 3. Results

## 3.1. Transcriptional profiles and the overall magnitude of immune-related genes in healed, asymptomatics and CL patients

The expression of 144 genes related to immune responses was determined using dcRT-MLPA. A total of 27 significant (*P*-values < 0.05) DEGs were identified in whole blood samples of the CL groups compared with the group of healthy individuals (Supplementary Table S2). Some of these DEGs were shared between two groups; however, none was significantly different in all three groups. Several DEGs were involved in IFN signaling (CD274, GBP1, GBP2, GBP5, IFI16, IFIT2, IFIT5, INDO, OAS1, OAS2, SOCS1, STAT1 and TAP2). Other DEGs were related to chemokines (CCL11 & CXCL9), pattern recognition receptors (MRC2 & TLR6), myeloid associated genes (CCL22 & IL12B), T~reg~ associated genes (TGFB1 & TNFRSF18), inflammation (TIMP2), Th1 associated genes (IL2), Th2 associated genes (IL5), T cell subset markers (CD4), cytotoxicity markers (PRF1) and GTPase activating proteins (RAB13). PCA was performed to compare the overall separation of transcriptomes among different cohorts. As depicted in Fig. 1, peripheral blood gene expression distinguished the patients from the healthy group; however, a similar separation could not be attained for either the healed or the asymptomatic group.

## 3.2. The magnitude of fold changes of DEGs in healed, asymptomatics and CL patients

We compared the magnitude of fold changes among 27 DEGs identified as differentially expressed in the whole blood of studied groups to those of the healthy controls. A volcano plot of Log2 fold changes versus the *P* values (-Log10) within the studied groups (Fig. 2), identified genes with large fold changes that were also statistically significant. Compared to the control, out of 144 tested genes, only a few in each group were significantly over-or under-expressed, and as such could potentially be considered as putative signature transcripts.

Moreover, the heatmap shown in Fig. 3 A compares Log2 fold changes of the 27 DEGs mentioned above in each of the main studied groups to the control. The numerical values of Log2 fold changes and the full names of the genes are indicated in Supplementary Table S3. Further, the gene expressions were also compared statistically between the patients and the asymptomatic groups as well as the patients and the healed groups (Supplementary Tables S4 and S5, respectively).

Among the DEGs in the healed group, SOCS1, IFI16, IFIT5, and STAT1 were significantly higher expressed whereas the expression of GBP5, OAS2 and TAP2 was significantly lower, compared to the healthy control group (Fig. 3A). CCL11, RAB13, CXCL9, INDO, IL12B and TNFRSF18 displayed higher expression for the asymptomatic group, while GBP2, GBP5, GBP1, IFIT2 and PRF1 were expressed to a lower extent. The patient group had a larger number of DEGs compared with the other two CL groups; among the patients, CD274, IFIT5, SOCS1, IL2, RAB13, STAT1, OAS1, TGFB1, INDO, TIMP2, and MRC2 were expressed to a greater extent while CCL11, IL5, TAP2, CCL22, and CD4 were expressed to a lower extent. Among the analyzed genes, similar expression patterns were noted for SOCS1, IFIT5, STAT1 and TAP2 in the healed and patient groups. The most profound alterations in fold-change were observed among the patient group members.

As indicated in the Venn diagram (Fig. 3B), the gene expression patterns changed in common between the healed and the patient group (i.e., SOCS1, IFIT5, STAT1 and TAP2) followed the same pattern of over- or under-expression. The healed and the asymptomatic group only had one common DEG (i.e., the under-expressed GBP5). Instead, the healed group shared three DEGs with the same directionality with the patients
<figure><img src="image_4.png" /><figcaption>Fig. 1. Principal Component Analysis (PCA) for dc-RTMLPA transcript data of whole blood obtained from healed (blue circles, n = 14), asymptomatic (orange circles, n = 17) and CL patients (red circles, n = 11), compared to the healthy control (green circles, n = 15).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. A Volcano plot depicting Log2 fold changes of the tested genes in the studied groups versus the *P* values (-Log10). Genes with non-significant (*P* > 0.05) differences are shown as black dots. Using the *P*-value cutoff, the over- and under expressed DEGs are shown as red and blue dots, respectively.</figcaption></figure>

(i.e., the over-expressed RAB13 and INDO). Interestingly, CCL11 was over-expressed in the asymptomatic and under-expressed in the patient groups.

When the number of involved pathways with respect to the DEGs were compared among the studied groups (Fig. 4), pathways related to IFN-signaling was common among all of them. Although the over-expression of SOCS1 and STAT1 were prominent in the healed and pa-tient groups only, the asymptomatic group was unique (especially compared to the patient group) for over-expression of the chemokines pathways due to CCL11 and CXCL9 genes.

### 3.3. Protein network analysis

Interacting protein networks analyzed by the STRING database identified Suppressor of cytokine signaling 1 (SOCS1) and Signal transducer and activator of transcription 1 (STAT1) proteins as major nodes in the CL patient and the healed group. In contrast, Guanylate binding proteins (GBP) formed the main nodes in the asymptomatic group. As shown in Fig. 5, the combination of the predicted networks and information gained by differential expression analysis point to the interactions of an axis formed by over-expressed SOCS1, STAT1, and IFN-induced protein with tetratricopeptide repeats 5 (IFIT5), only for the healed and the patient group, while the impact of the under-
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** (A) A heatmap comparing the 27 DEGs in terms of Log2 fold changes of the studied groups with the controls. Intensities of red and blue colors represent relative increased (>> 0) or decreased (< 0) expression of the genes in terms of Log2 fold changes, respectively. Boxes in grey color indicate no significant change in gene expression relative to the control group. (B) A Venn diagram is indicating common DEGs among the main studied groups. The numbers indicate the total number of common DEGs between the groups. Genes expressed significantly higher or lower than the control are shown in red or blue fonts, respectively, except for *CCL11* which was expressed in opposite directions in the asymptomatic (over-expressed) and the patient (under-expressed) groups.</figcaption></figure>

expressed GBP5 was an integral part of the network recognized for the asymptomatic group. However, the network identified in the asymptomatic group included IFN-signaling genes that were all down-regulated (GBP1/2/5/ and IFIT2), while the patients' network contained IFN-signaling genes that were all up-regulated (SOCS1, OAS1, IFIT5 and STAT1). This finding suggests that transcriptional activation of systemic IFN signaling could contribute to the active state of the disease.

## 4. Discussion

In many parts of North Africa, West Asia, and particularly Iran, CL results from an infection caused by *L. major* and/or *L. tropica* (Mahmoudzadeh-Niknam et al., 2012). In order to determine factors related to long-term immunity and protective responses against CL, aside from patients involved with the active cutaneous lesions, information based on human populations that have resolved or withstood the formation of CL lesions, such as the healed and the asymptomatic individuals in the endemic areas, are deemed to be crucial. This study sought to compare the immune-focused genes expressions in the whole blood of a healthy population residing in CL endemic areas with the three groups mentioned above. We herein used the dcRT-MLPA method to study the expression of 144 genes involved in human immune responses combined with a systems biology approach to study the transcriptomic blood profiles of the CL patients healed and asymptomatic individuals in CL endemic areas. Although the gene expression profiling of 144 genes does not represent a holistic picture of the whole blood inflammatory gene expression profile of the studied groups, our data demonstrated that distinct and overlapping expressions of the inflammatory genes could be detected at the systemic level in different clinical stages of CL.

Considering the documented dual role of IFN in the immunopathogenesis of CL, it was not surprising to observe that 13 out of 27 DEGs identified collectively in the healed, asymptomatic, and CL patients belonged to the IFN signaling pathways. We have recently reported IFN signaling as one of the main five clusters of enriched pathways identified in the skin lesions of ulcerative and non-ulcerative CL patients caused by

*L. tropica* by RNA sequencing (Masoudzadeh et al., 2020b). Moreover, by using RNA sequencing, Amorim et al., have recently reported IFN-stimulated genes signature as the dominant transcriptional response in the blood of *L. braziliensis*-infected CL patients (Amorim et al., 2021). However, we observed for the first time that the components of this pathway are expressed distinctly in different clinical stages of CL due to *L. tropica*. Our *in silico* protein network analysis pointed to an over-expressed SOCS1-STAT1-IFIT5 axis in the healed and the patient groups. The over-expressed SOCS1 in the patient group was shown to be connected to the over-expressed IL-2 node directly as well as the over-expressed STAT1 node, which can potentially favor a Th1 response in the patient group. Both SOCS and STAT proteins are considered as critical mediators of IFN signaling. The expression of the SOCS1 gene was significantly higher in the healed and CL patients but not in the asymptomatic group, compared to that of the healthy control.

SOCS1 belongs to the suppressor of cytokine signaling (SOCS) family that encodes intracellular proteins with regulatory functions in the immune system. These proteins can inhibit inflammatory cascades by inhibiting Janus Kinase/Signal Transducer and Activator of Transcription (JAK/STAT) pathways and TLR signaling (Serezani et al., 2011; Sharma and Larkin, 2019; Scott and Novais, 2016). The family of intracellular tyrosine kinases, such as JAK/STAT, are essential factors with respect to CL, as they play a pivotal role in transducing extracellular cytokine signaling and thereby contribute to major processes such as cellular proliferation, homeostasis, and host defense against the infection (Olivier et al., 2005). It has been demonstrated that STAT1-mediated IFN signaling is crucial in developing Th1 response against CL caused by *L. major* (Rosas et al., 2003). The involvement of SOCS1 in the pathogenesis of leishmaniasis has previously been shown in animal models of the disease due to other *Leishmania* species than *L. tropica*. For instance, mice lacking a single copy of the *SOCS1* gene (*SOCS1 +/−*) were shown to develop larger cutaneous lesions than wild-type mice when given a moderate dose of *L. major* promastigotes. Interestingly, these heterozygous mice sustained a similar parasite burden in their draining lymph nodes and exhibited a normal Th1 type immune response. This suggests that *SOCS1* has a regulatory role for the inflicted pathology; however, in a different manner from the immune clearance of the parasite (Bullen et al., 2003). Another animal model study on VL has shown that microRNA 155 (miR155) promotes Th1 type response with production of IFN-γ, which dampens disease severity by targeting *SOCS1* gene expression (Varikuti et al., 2019). Recently, it has been shown that *L. donovani* induces over-expression of SOCS1 and SOCS3 in macrophages and T-cells, which leads to impairment of macrophage-T cell cross-talk by interfering with macrophage IFN-γ signaling as an immune-evasion strategy (Chandrakar et al., 2020). Moreover, it has been shown that a decrease in the expression of SOCS3 due to a genetic variability might be linked to susceptibility to American CL (Sanchez Luquez et al., 2020). Another gene of interest was *IFIT5*, observed to be over-expressed in both the healed and the patient groups. The gene product has close interaction with both SOCS1 and STAT1 (Fig. 5). The role of *IFIT5* in CL has previously been reported in a study in which *L. major*, but not *L. donovani*, was shown to induce IFIT5 expression, leading to IL-12 production and a strong Th1 response (Favila et al., 2014). These results indicate that SOCS1 acts as a regulatory protein across a wide range of *Leishmania* species, which could have potential as a therapeutic target against different types of leishmaniases.

The under-expressed *GBP5* was the only DEG in common between the healed and the asymptomatic groups and its product was an identifiable node in their expected protein networks. *PRF1* and all tested Guanylate-binding protein genes, namely *GBP1*, *GBP2* and *GBP5*, were significantly under-expressed in the asymptomatic group. The *PRF1* gene encodes perforin, a pore-forming protein that can facilitate the entry of granzymes through the cellular membrane of invading *Leishmania* species, leading to apoptosis (Amorim et al., 2021). The pathologic roles of *PRF1* expression in CD8⁺ cytotoxic T cell machinery are
Differentially-Expressed Genes

<figure><img src="image_7.png" /><figcaption>Asymptomatic</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Patients</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Pathways</figcaption></figure>

Fig. 4. Comparison of the number of involved pathways concerning the immune-related DEGs in whole blood obtained from healed, asymptomatic and CL patients. The Y-axes indicate significant DEGs (*P* < 0.05) and the X-axes indicate the pathways (specified in *Supplementary Table S1*). The colors and sizes of the bubbles represent Log2FC and the *P* values of significant genes mapped to the indicated pathways, respectively. Positive and negative Log2FC denote up- and down-regulated pathways, respectively.
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** *In silico* comparison of multiple proteins interactions for functional protein association networks in the studied groups. The connected nodes computed by the STRING database are marked as shades of red and blue circles based on Log2 FC in each group to indicate up- and down-regulated protein-encoding genes, respectively. Line colors indicate interactions between the proteins: pink, as experimentally determined known interaction; blue and light purple as gene co-occurrence and protein homology predicted interactions, respectively; and black as co-expression.</figcaption></figure>

shown using murine models of CL due to *L. braziliensis* (Novais et al., 2013). Similar models have confirmed the pathological role of cytolytic CD8⁺ T cells in inflammatory responses (Novais et al., 2021) as well as inflammasome activation (Novais et al., 2017). Limited information is currently available regarding the impact of GBP family of IFN-inducible GTPases on CL in humans. It is known that they can regulate inflammasomes, which in turn affect inflammatory caspases that promote pathogen clearance and/or death of the infected cells while offering a defensive role against leishmaniasis (Sohrabi et al., 2018). A microarray genome-wide transcriptional analysis of 25 *L. braziliensis* CL-biopsies has indicated that *GBP5* is among the top 10-induced genes in cutaneous lesions (Novais et al., 2015). Also, in another transcriptomic study of three *L. tropica* CL lesions, *GBP5* has been reported to be among the most up-regulated genes (Masoudzadeh et al., 2017). Since we did not observe the over-expression of *GBP5* in our patient's group, it appears that *GBP5* is regulated differentially at the systemic and the local levels. Altogether, the comparative under-representation of *PRF1* and *GBP* genes in the asymptomatic group as well as the lower expression of *GBP5* in the healed group, could be indications of controlled cytolytic activities due to the infection in these groups which at least in part, contribute to a balance between controlling the infection while avoiding pathological lesion formation. Furthermore, as suggested by Amorim et al., *PRF1* and *GBP* expressions are associated with increased frequencies of cytotoxic T cells and myeloid cells, respectively (Amorim et al., 2021). Here, the decreased expressions of *PRF1* and *GBP* in the blood of the asymptomatic (and *GBP5* in the healed) subjects could possibly reflect a reduced number of those circulating cells which would contribute to reduced systemic inflammation and pathogenesis. This can also suggest that the induced IFN-γ in the blood of CL patients can potentially activate the IFN signaling genes and downstream mechanisms in cells expressing the IFN-γ receptor. This action could potentially prime monocytes before entering the lesion as a protection and parasite control mechanism (Amorim et al., 2021). However, our results suggest that the recruitment of cells that are active for IFN-γ/GBP signaling has a pathogenic role in the active disease.

We also noted a unique transcriptomic profile concerning the asymptomatic group where *CCL11*, *CXCL9*, *IL12B*, and *TNFRSF18* were over-expressed while *GPB1*, *GPB2*, *IFIT2* and *PRF1* were under-expressed. Moreover, distinct profiles comprised of the over-expressed *TLR6* and *IFI16* and the under-expressed *OAS2* were detected in the healed group. The most over-expressed gene in the patient group was *CD274*. Interestingly, none of these profiles had a common gene with the patient group, except for *CCL11* gene expression, which exhibited opposite directions in the asymptomatic and the patient group. This

suggests that *CCL11* can be considered as one of the candidates for further assessments as a putative novel biomarker of asymptomatic CL in the endemic area (Fig. 2 and Fig. 3). Determination of the asymptomatic population in CL endemic areas for epidemiological and vaccine studies is a priority. There is also an increasing demand for the development of minimally invasive alternative methods to replace LST. In this regard, our results could help to inform the development of blood biomarkers of CL asymptomaticity in the endemic area. Additionally, since eotaxin-1 encoded by *CCL11* is involved in recruitment of several types of leukocytes such as eosinophils, basophils, and macrophages during inflammation (Menzies-Gow et al., 2002; Yoshihisa et al., 2011), its over-expression in the asymptomatic group may presumably contribute to pre-activating of the epithelial cells to contain the parasites at their entry port.

Nevertheless, some regulated genes were similarly-expressed among the groups, which cannot potentially be used for the purpose of distinguishing different clinical stages of CL. For instance, the healed and the patient groups had the highest number of common DEGs with *SOCS1*, *IFIT5*, *STAT1* and *TAP2* while displaying the same directionality of (over- or under-) expression. Moreover, NOD-like receptor signaling was a common pathway in all the studied groups, suggesting its involvement in different clinical stages of CL.

However, our analyses revealed that DEGs in the patient group were significantly enriched for several additional pathways, compared with the data from the healed and the asymptomatic groups. This may, at least in part, be explained by the fact that unlike the healed and the asymptomatic individuals, the patients' immune system is actively involved in countering the pathogen and tissue damage repair mechanisms. Moreover, it should be noted that such pathway analyses revealed more significant differences between the groups than when the genes fold-changes were compared. This probably reflects the fact that the identified DEGs are involved in many different biological pathways and not all of the mentioned pathways were simultaneously affected in these groups. Although a comparison of the magnitude of gene expression between RNA-Seq and dcRT-MLPA datasets is not directly possible, we compared the DEGs in whole blood of CL patients studied herein with those of the skin determined by RNA sequencing reported previously by us (Masoudzadeh et al., 2020b). Interestingly, we found 6 common genes (*CXCL9*, *IFI16*, *GBP5*, *GBP2*, *GBP1* and *MRC2*) which among them, *CXCL9* and *IFI16* were similarly over-expressed in both skin and blood of the CL patients, albeit with lower magnitude in the blood. In general, the observed fold changes in the blood were relatively modest, which may be explained by the fact that CL is a localized skin disease.

The gene encoding cluster of differentiation 274 (CD274), also
known as programmed death-ligand 1 (PD-L1), was identified as a unique DEG in the patient group. The gene product is a 40-kDa type 1 transmembrane protein that has a regulatory role in T cell response by reducing the proliferation of antigen-specific T cells in lymph nodes and inducing apoptosis in regulatory T cells (Filippis et al., 2017). Blockade of this protein by specific antibodies has been shown to be helpful for treatment of the chronic forms of leishmaniasis in the animal models of the disease (Amorim et al., 2019; Scott and Novais, 2016). The blockade of PD-L1 leads to increased expression of cytolytic molecules, such as perforin, granzyme B, and granulysin in T cells of tuberculosis and cancer patients (Kim et al., 2016). Our group has recently identified PD-1 signaling among common positively enriched pathways in the skin lesions of ulcerative and non-ulcerative CL patients infected with *L. tropica* (Masoudzadeh et al., 2020b). Interestingly, our data herein indicated that the expression of PD-L1 gene could also be detected in the blood of patients with active CL disease.

## 5. Conclusions

In this study, we analyzed the expression profiles of 144 immune-related genes in the blood of healed and asymptomatic patients in CL endemic areas. Several DEGs and biological pathways were either shared or were specific to the different disease stages. Remarkably, our results showed the involvement of IFN-signaling, and most notably, the SOCS1 pathway in individuals with active or healed CL lesions. This study enhanced our understanding of the systemic immune response, involved during different clinical stages of CL. The components of the identified distinct transcriptomic profiles in each studied groups, particularly *IFI16* for the healed and *CCL11* for the asymptomatic groups may represent putative blood biomarkers. These results warrant further investigation of blood biomarkers of different stages of CL and especially, the asymptomatic population, living in CL endemic areas.

## Funding

This research was funded by research grant No. 992 and contract No 1109 and 1111 from Pasteur Institute of Iran. It has also received funding from the European Union's H2020 LeishShield MATI RISE (Research and Innovation Staff Exchange), under the Marie Skłodowska-Curie grant agreement No 778298.

## CRediT authorship contribution statement

F Bahrami, S Rafati, and AM Harandi conceived the project, wrote and revised the manuscript. MC Haks and THM Ottenhoff supervised dcRT-MLPA and revised the manuscript. F Bahrami, S Rafati, I Sharifi, V Mashayekhi Goyonlo and N Masoudzadeh collected the blood samples and supervised LST. F Bahrami and N Masoudzadeh extracted RNA. N Masoudzadeh and Y Taslimi performed PCR diagnostic of the patients. S Van Veen conducted dcRT-MLPA. I Sharifi and J Persson revised the manuscript. N Masoudzadeh, A Lari, H Sarvnaz, B Andersson and M Östensson analyzed the data statistically. N Masoudzadeh and A Lari performed bioinformatics analyses and N Masoudzadeh plotted the graphs. All authors approved the final version of this manuscript.

## Conflict of Interest

The authors report no declarations of interest.

## Acknowledgments

The authors are grateful to the volunteers in Kerman and Mashhad for their participation in this study, Davoud Iravani for his assistance in performing LST and blood collections, and the staff of Kerman and Mashhad Universities of Medical Sciences for their administrative assistance. The authors also thank Safoora Gharibzadeh (Epidemiology

and Biostatistics Dept., Pasteur Institute of Iran) for her contribution in biostatistical analyses.

## Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.molimm.2022.07.008.

## References

Alimohammadian, M.H., Kivanjah, M., Pak, F., Gaznavia, A., Kharazmi, A., 1993. Evaluation of the efficacy of Iran leishmanin and comparison with leishmanins from Wellcome (UK) and Roma (Italy) in cured cutaneous leishmaniasis patients. Trans. R. Soc. Trop. Med. Hyg. 87, 550-551.

Alvar, J., Velez, I.D., Bern, C., Herrero, M., Desjeux, P., Cano, J., Jannin, J., Den Boer, M., 2012. Leishmaniasis worldwide and global estimates of its incidence. PLoS One 7, e35671.

Amorim, C.F., Novais, F.O., Nguyen, B.T., Misic, A.M., Carvalho, L.P., Carvalho, E.M., Beiting, D.P., Scott, P., 2019. Variable gene expression and parasite load predict treatment outcome in cutaneous leishmaniasis. Sci. Transl. Med. 11.

Amorim, C.F., Novais, F.O., Nguyen, B.T., Nascimento, M.T., Lago, J., Lago, A.S., Carvalho, L.P., Beiting, D.P., Scott, P., 2021. Localized skin inflammation during cutaneous leishmaniasis drives a chronic, systemic IFN-gamma signature. PLoS Negl. Trop. Dis. 15, e0009321.

Bahrami, F., Harandi, A.M., Rafati, S., 2018. Biomarkers of cutaneous leishmaniasis. Front. Cell Infect. Microbiol. 8, 222.

Banuls, A.L., Bastien, P., Pomares, C., Arevalo, J., Fisa, R., Hide, M., 2011. Clinical pleiomorphism in human leishmaniases, with special mention of asymptomatic infection. Clin. Microbiol Infect. 17, 1451-1461.

Bullen, D.V., Baldwin, T.M., Curtis, J.M., Alexander, W.S., Handman, E., 2003. Persistence of lesions in suppressor of cytokine signaling-1-deficient mice infected with *Leishmania major*. J. Immunol. 170, 4267-4272.

Chandrakar, P., Parmar, N., Descoteaux, A., Kar, S., 2020. Differential Induction of SOCS Isoforms by *Leishmania donovani* Impairs Macrophage-T Cell Cross-Talk and Host Defense. J. Immunol. 204, 596-610.

Favila, M.A., Geraci, N.S., Zeng, E., Harker, B., Condon, D., Cotton, R.N., Jayakumar, A., Tripathi, V., Mcdowell, M.A., 2014. Human dendritic cells exhibit a pronounced type I IFN signature following *Leishmania major* infection that is required for IL-12 induction. J. Immunol. 192, 5863-5872.

Filippis, C., Arens, K., Noubissi Nzeteu, G.A., Reichmann, G., Waibler, Z., Crauwels, P., Van Zandbergen, G., 2017. Nivolumab enhances in vitro effector functions of PD-1(+) T-lymphocytes and leishmania-infected human myeloid cells in a host cell-dependent manner. Front. Immunol. 8, 1880.

Fisch, D., Bando, H., Clough, B., Hornung, V., Yamamoto, M., Shenoy, A.R., Frickel, E.M., 2019. Human GBP1 is a microbe-specific gatekeeper of macrophage apoptosis and pyroptosis. EMBO J. 38, e100926.

Gardinassi, L.G., Garcia, G.R., Costa, C.H., Costa Silva, V., De Miranda Santos, I.K., 2016. Blood transcriptional profiling reveals immunological signatures of distinct states of infection of humans with *Leishmania infantum*. PLoS Negl. Trop. Dis. 10, e0005123.

Haks, M.C., Goeman, J.J., Magis-Escurra, C., Ottenhoff, T.H., 2015. Focused human gene expression profiling using dual-color reverse transcriptase multiplex ligation-dependent probe amplification. Vaccine 33, 5282-5288.

Joosten, S.A., Goeman, J.J., Sutherland, J.S., Opmeer, L., De Boer, K.G., Jacobsen, M., Kaufmann, S.H., Finos, L., Magis-Escurra, C., Ota, M.O., Ottenhoff, T.H., Haks, M.C., 2012. Identification of biomarkers for tuberculosis disease using a novel dual-color RT-MLPA assay. Genes Immun. 13, 71-82.

Kim, B.H., Chee, J.D., Bradfield, C.J., Park, E.S., Kumar, P., Macmicking, J.D., 2016. Interferon-induced guanylate-binding proteins in inflammasome activation and host defense. Nat. Immunol. 17, 481-489.

Mahmoudzadeh-Niknam, H., Ajdary, S., Riazi-Rad, F., Mirzadegan, E., Rezaeian, A., Khaze, V., Djadid, N.D., Alimohammadian, M.H., 2012. Molecular epidemiology of cutaneous leishmaniasis and heterogeneity of *Leishmania major* strains in Iran. Trop. Med. Int. Health 17, 1335-1344.

Masoudzadeh, N., Mizbani, A., Taslimi, Y., Mashayekhi, V., Mortazavi, H., Sadeghipour, P., Ardekani, H.M., Rafati, S., 2017. *Leishmania tropica* infected human lesions: Whole genome transcription profiling. Acta Trop. 176, 236-41.

Masoudzadeh, N., Mizbani, A., Rafati, S., 2020a. Transcriptomic profiling in Cutaneous Leishmaniasis patients. Expert Rev. Proteom. 17, 533-41.

Masoudzadeh, N., Östensson, M., Persson, J., Mashayekhi Goyonlo, V., Abgajogu, C., Taslimi, Y., Erfanian Salim, R., Zahedifard, F., Mizbani, A., Malekafzali Ardekani, H., Gunn, B.M., Rafati, S., Harandi, A.M., 2020b. Molecular signatures of anthroponotic cutaneous leishmaniasis in the lesions of patients infected with *Leishmania tropica*. Sci. Rep. 10, 16198.

Menzies-Gow, A., Ying, S., Sabroe, I., Stubbs, V.L., Soler, D., Williams, T.J., Kay, A.B., 2002. Eotaxin (CCL11) and eotaxin-2 (CCL24) induce recruitment of eosinophils, basophils, neutrophils, and macrophages as well as features of early- and late-phase allergic reactions following cutaneous injection in human atopic and nonatopic volunteers. J. Immunol. 169, 2712-2718.

Novais, F.O., Carvalho, L.P., Graff, J.W., Beiting, D.P., Ruthel, G., Roos, D.S., Betts, M.R., Goldschmidt, M.H., Wilson, M.E., De Oliveira, C.I., Scott, P., 2013. Cytotoxic T cells mediate pathology and metastasis in cutaneous leishmaniasis. PLoS Pathog. 9, e1003504.
Novais, F.O., Carvalho, L.P., Passos, S., Roos, D.S., Carvalho, E.M., Scott, P., Beiting, D.P., 2015. Genomic profiling of human *Leishmania braziliensis* lesions identifies transcriptional modules associated with cutaneous immunopathology. J. Investig. Dermatol. 135, 94-101.

Novais, F.O., Carvalho, A.M., Clark, M.L., Carvalho, L.P., Beiting, D.P., Brodsky, I.E., Carvalho, E.M., Scott, P., 2017. CD8+ T cell cytotoxicity mediates pathology in the skin by inflammasome activation and IL-1beta production. PLoS Pathog. 13, e1006196.

Novais, F.O., Nguyen, B.T., Scott, P., 2021. Granzyme B inhibition by tofacitinib blocks the pathology induced by CD8 T cells in cutaneous leishmaniasis. J. Investig. Dermatol. 141, 575-85.

Olivier, M., Gregory, D.J., Forget, G., 2005. Subversion mechanisms by which *Leishmania* parasites can escape the host immune response: a signaling point of view. Clin. Microbiol. Rev. 18, 293-305.

Rosas, L.E., Keiser, T., Pyles, R., Durbin, J., Satoskar, A.R., 2003. Development of protective immunity against cutaneous leishmaniasis is dependent on STAT1-mediated IFN signaling pathway. Eur. J. Immunol. 33, 1799-1805.

Sanchez Luquez, K.Y., Zerpa, O., Paz-Villarraga, C.A., Fernandez-Mestre, M., 2020. Genetic variability of molecules involved in the disease pathogenesis in *Leishmania* infection. Exp. Parasitol. 218, 108007.

Scott, P., Novais, F.O., 2016. Cutaneous leishmaniasis: immune responses in protection and pathogenesis. Nat. Rev. Immunol. 16, 581-92.

Serezani, C.H., Lewis, C., Jancar, S., Peters-Golden, M., 2011. Leukotriene B4 amplifies NF-kappaB activation in mouse macrophages by reducing SOCS1 inhibition of MyD88 expression. J. Clin. Investig. 121, 671-682.

Sharma, J., Larkin 3rd, J., 2019. Therapeutic implication of SOCS1 modulation in the treatment of autoimmunity and cancer. Front Pharmacol. 10, 324.

Sohrabi, Y., Volkova, V., Kobets, T., Havelkova, H., Krayem, I., Slapnickova, M., Demant, P., Lipoldova, M., 2018. Genetic regulation of guanylate-binding proteins 2b and 5 during Leishmaniasis in mice. Front Immunol. 9, 130.

Sokal, J.E., 1975. Editorial: measurement of delayed skin-test responses. New Engl. J. Med. 293, 501-502.

Taslimi, Y., Sadeghipour, P., Habibzadeh, S., Mashayekhi, V., Mortazavi, H., Muller, I., Lane, M.E., Kropf, P., Rafati, S., 2017. A novel non-invasive diagnostic sampling technique for cutaneous leishmaniasis. PLoS Negl. Trop. Dis. 11, e0005750.

Varikuti, S., Natarajan, G., Volpedo, G., Singh, B., Hamza, O., Messick, G.V., Guerau-De-Arellano, M., Papenfuss, T.L., Oghumu, S., Satoskar, A.R., 2019. MicroRNA 155 contributes to host immunity against *Leishmania donovani* but is not essential for resolution of infection. Infect. Immun. 87, e00307-e00319.

Yoshihisa, Y., Makino, T., Matsunaga, K., Honda, A., Norisugi, O., Abe, R., Shimizu, H., Shimizu, T., 2011. Macrophage migration inhibitory factor is essential for eosinophil recruitment in allergen-induced skin inflammation. J. Investig. Dermatol. 131, 925-931.

## Glossary

**CL**: Cutaneous leishmaniasis (CL)

**VL**: Visceral leishmaniasis (VL)

**LST**: Leishmanin skin test

**dcRT-MLPA**: Dual-color reverse-transcription multiplex ligation-dependent probe amplification