article: Localized Cutaneous Leishmaniasis (LCL) and Mucosal Leishmaniasis (ML) are two extreme clinical forms of American Tegumentary Leishmaniasis that usually begin as solitary primary cutaneous lesions. Host and parasite factors that influence the progression of LCL to ML are not completely understood. In this manuscript, we compare the gene expression profiles of primary cutaneous lesions from patients who eventually developed ML to those that did not. Using RNA-seq, we analyzed both the human and Leishmania transcriptomes in primary cutaneous lesions. Limited number of reads mapping to Leishmania transcripts were obtained. For human transcripts, compared to ML patients, lesions from LCL patients displayed a general multi-polarization of the adaptive immune response and showed up-regulation of genes involved in chemoattraction of innate immune cells and in antigen presentation. We also identified a potential transcriptional signature in the primary lesions that may predict long-term disease outcome. We were able to simultaneously sequence both human and Leishmania mRNA transcripts in primary cutaneous leishmaniasis lesions. Our results suggest an intrinsic difference in the immune capacity of LCL and ML patients. The findings correlate the complete cure of L. braziliensis infection with a controlled inflammatory response and a balanced activation of innate and adaptive immunity.
  Leishmaniasis is an important parasitic disease found in 88 tropical and subtropical countries and is still one of the most neglected diseases in the world. Approximately 12 million people are infected with the Leishmania protozoa, 270 thousand new cases occur yearly, and 1. 6 million individuals are considered at risk of infection. American tegumentary leishmaniasis (ATL) is a serious public health problem in Latin America. In Brazil, over 26 thousand new cases are registered annually and most of the infections are caused by Leishmania braziliensis [1], [2]. Localized cutaneous leishmaniasis (LCL) is the most common manifestation of ATL and is characterized by the development of a single lesion at the site of the sandfly' s bite that tends to heal spontaneously and gradually even if left untreated [3]. However, the parasite dissemination from the primary cutaneous lesion to distant skin and mucosal regions, without lesion development, is a common event. L. braziliensis infection can still be detected sometimes after clinical cure [4]. Fortunately, less than 5% of the individuals infected with L. braziliensis may progress to mucosal leishmaniasis (ML), which results in ulceration of oropharynx and nasal mucosa either at the time of primary skin lesion or several months or years later [5]. The natural history of ML consists in the progressive destruction of the nasal septa and soft and hard palates, eventually causing severe facial disfiguration and respiratory disturbances that do not heal spontaneously [6]. Moreover, 22–50% of the individuals with ML are unresponsive to the current antimonial therapy and have an associated increase in morbidity and mortality [7]. Susceptibility to ML progression has been associated with polymorphism of TNF (−308 bp), IL6 (−174G/C) and CCL2/MCP1 (−2518A/G) genes [8]–[10]. This, together with the finding of an increased frequency of HLA-DQw3 allele in ML patients, suggests that host genetic factors may play an important role in individual tendency to ML development [11]. However, these findings alone cannot explain susceptibility patterns entirely. Following the genome sequencing of humans and L. braziliensis, and the development of Next Generation Sequencing (NGS) technologies [12], [13], we employed the RNA-seq method to analyze the transcriptomes of human and L. braziliensis in a high-throughput manner. The resulting data reveal both the identity and quantity of the different transcripts present from both host and pathogen simultaneously in infected tissue samples. A better understanding of the immune environment of the primary infection site may allow us to identify the mechanisms involved in successful host defense, predict the factors involved in the disease outcome and potential pathways for targeted therapeutic intervention. We hypothesized that the quality of the host immune response at the primary cutaneous infection with L. braziliensis is responsible for the ultimate clinical outcome of the disease (cure or development of ML) and can be decoded by specific transcriptional profiles. To test this hypothesis, we characterized and compared the transcriptomes of the primary cutaneous lesions from both LCL and ML patients (before they developed ML). In addition to human transcripts, we were able to detect the L. braziliensis transcripts in situ, as well.
  Skin tissue fragments were obtained from biopsies of active primary cutaneous leishmaniasis lesions from 10 patients before starting treatment with pentavalent antimony. The biopsy procedure was performed under anesthesia and in sterile conditions. The fragments were collected from the border of the ulcer, including epidermis and dermis, immediately embedded in O. C. T. medium (Tissue-Tek, Sakura Finetek, USA) and preserved in liquid nitrogen. All patients lived in Rio de Janeiro state, Brazil, did not present any symptom of mucosal lesions and no other pathologies. After diagnose confirmation for infection with L. braziliensis, those patients were treated with Glucantime (Rhodia Laboratories, Antony, France). All patients completely resolved the primary cutaneous lesions after treatment, with no lesion relapse or need for repeat treatment. Five individuals cured completely after therapy (LCL group) and 5 individuals developed mucosal lesions long after initial healing of the primary lesion (ML group). The age of the patients ranged between 15–67 years in LCL group and 23–63 years in ML group. All five LCL individuals and 3 ML subjects developed only one primary cutaneous lesion. Two ML subjects had more than one cutaneous lesion. By patient recall, the average (± standard deviation) time of lesion evolution was similar in both groups (LCL: 1. 4±0. 5 months; ML: 2. 3±1. 3 months). The criteria for inclusion in this study, therefore, were the time of evolution of the lesion, which should range between 1–4 months; the good response for the treatment, which comprehended the completely re-epithelialization of the cutaneous lesion and no lesion reactivation; and the final outcome of the disease, which should be no disease relapse (LCL) or development of lesions at the nasopharynx mucosa (ML). The research was approved by the Research Ethics Committee and IRB (Institutional Review Board) from both FIOCRUZ and LABIOMED Institutes. Written informed consent was obtained from all the adults and from the parents or legal guardians on behalf of the participants under 18 years of age. The study was conducted in agreement with the principles of the Helsinki Declaration and the Resolution 196/96 of the National Health Council of the Ministry of Health that regulates research involving human subjects in Brazil. The total RNA was isolated using guanidine-chloroform-phenol nucleic acid extraction, treated with DNAse (Qiagen, Valencia, USA) and cleaned with Rneasy Micro-kit (Qiagen). The quality control step was performed using Bioanalyzer (Agilent Technologies, Santa Clara, USA). Afterward, the total RNA was used to prepare the transcriptome libraries with the mRNA-seq sample preparation kit (Illumina, San Diego, USA) according to the manufacturer' s instructions and sequenced using the Illumina Genome Analyzer II platform. The sequence of all cDNA libraries will be freely available on the GEO NCBI database (http: //www. ncbi. nlm. nih. gov/geo), accession number GSE33601. Reads were aligned to the references genomes using the Bowtie aligner and Tophat program. The reference genomes used in this project were from Homo sapiens (UCSC- hg19), Leishmania braziliensis (TritrypDB, version 3. 2), Staphylococcus aureus (Wellcome Trust Sanger Institute - EMRSA15) and Pseudomonas aeruginosa (Wellcome Trust Sanger Institute - LEB58) (Table S8). SAMtools were used for quality control and processing the reads (Table S9). Cuffdiff from the Cufflinks' package was applied to normalize the number of human reads to FPKM (Fragments Per Kilobase per Million) considering gene isoforms and including multi-hit reads (up to 30 regions - Table S2). Student' s t-test was used for dimensional reduction of the data set, considering as threshold P-values lower than 0. 05 and the software SAM (Significance Analysis of Microarray [14]) was selected to perform the multiple testing correction, considering q-values lower than 0. 01. Afterward, the genes with a minimum fold change of gene expression of 1. 5 were considered for biological interpretation. Linear regression analysis was used to evaluate the similarity of gene expression between LCL and ML samples. A MA-plot was used to visualize the gene expression distribution in both groups. A rarefaction curve was calculated for each sample (method described in Text S1) to evaluate if the sequencing depth has provided a sufficient coverage of the transcriptome. Leave-One-Out cross validation analysis was performed as previously described [15], [16] to determine the phenotypic predictive power of the top differentially expressed genes. Functional annotation and canonical pathway identification of the selected genes was performed by Ingenuity Pathway Analysis (IPA) software (Ingenuity System, Redwood City, CA, http: //www. ingenuity. com), and by Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources 6. 7 (NIAID-NIH, http: //david. abcc. ncifcrf. gov/home. jsp). The Euclidean distance method was used for hierarchical clustering of the samples using Gene Cluster 3. 0 software (http: //rana. lbl. gov/EisenSoftware. htm).
  To determine the relative abundance of mRNA of individual genes for both host and pathogen simultaneously and in an unbiased manner, we used the RNA-Seq technique. The samples used for this purpose were skin biopsy fragments obtained from the border of primary cutaneous lesions of patients infected with L. braziliensis from Rio de Janeiro, Brazil. The sequencing generated 13±5 million reads per sample. The reads were aligned to both human (hg19) and L. braziliensis (Tritryp, vs 3. 2) genomes using the Bowtie and Tophat algorithms. This annotation of the human genome has 23,136 genes, while the L. braziliensis (Lbr) genome has 8,504 genes. The results indicated a good efficiency of the RNA-seq assay to identify transcripts from both host and parasite genomes in the cutaneous lesions of human leishmaniasis. On average, 20% of the reads were not assigned to human or L. braziliensis genomes. Almost 80% of the total number of reads aligned to the human genome, while less than 1% of the total reads (range 0. 06–0. 91%) were aligned to L. braziliensis. The ratios of misalignment between human and parasite genomes ranged from 0. 002% to 0. 01% (Table S1). As the absolute number of L. braziliensis reads detected was insufficient to effectively cover the L. braziliensis transcriptome, we did not perform differential gene expression analysis on the parasite transcripts. Bacterial RNA would not be predicted to be detected due to lack of polyadenylation. To verify that bacterial RNA was not present in the samples to a high degree, we aligned the reads to Staphylococcus aureus (EMRSA15) and Pseudomonas aeruginosa (LEB58) genomes, both representing the common bacteria causing secondary infection in cutaneous leishmaniasis lesions in Brazil [17]. Less than 0. 005% of the unaligned reads matched to these potential contaminants. To estimate if the RNA-Seq sequencing depth provided a satisfactory coverage of expressed human genes, we calculated rarefaction curves using reads that were mapped to unique regions of the human genome. The curves indicate that we have achieved good coverage in all samples. Regardless the total number of reads, in all samples we were able to detect a similar number of expressed genes (Figure S1). For gene expression analysis, unique and multi-hit reads that matched up to 30 regions to the human genome were considered and normalized to FPKM. The average percentage of reads that mapped in more than 2 regions was 23% (Table S2). Reads that aligned to more than 30 regions were excluded. The similarity of gene expression patterns of LCL and ML groups was analyzed by linear regression, which indicates that approximately 90% of the genes had comparable levels of expression (r2 = 0. 91; Figure 1A). This similarity can also be observed in the MA-plot (Figure 1B). Additionally, using all genes that were expressed in at least one sample (18,577 genes), we performed hierarchical clustering on the samples using Euclidean distance (Figure 2 left panel). The heat map shows that all of the samples were correctly clustered into the two clinical phenotypic groups, ML and LCL. Before determining the differences in the genes expression and in the biological pathways between the LCL and ML groups, we performed a dimensional reduction of the data set using Student' s t-test, considering as threshold P-values lower than 0. 05. This step resulted in the selection of 1208 genes (Figure 2 right panel). From this gene list, 58. 5% (707) of the genes had increased expression in the LCL samples, while 41. 5% (501) were increased in ML samples. Our choice of multiple testing correction was the software SAM (Significance Analysis of Microarray) which includes a special statistic called d statistic (true null is centered in 0) and includes randomization and permutation. We performed the multiple testing correction using the FPKM values. Considering the resulting q-value lower than 0. 01 for significance, 1130 genes were selected. As an additional filter, we established a minimum fold change in the expression of 1. 5. This step resulted in the selection of 867 genes (434 genes up-regulated in LCL group and 433 genes up-regulated in ML group). Those genes were considered as differentially expressed and potentially biologically relevant, and were then used for further biological interpretation. We hypothesized that the genes that were differentially expressed between the ML (n = 5) and LCL (n = 5) clinical subgroups would represent important biological differences and would thus associate with different functional categories. We used Ingenuity Pathway Analysis (IPA) to determine the functional annotation of the differentially expressed genes (Table S3). As hypothesized, we observed several common biological activity pathways in both ML and LCL comprised of genes with distinct functions (Tables S4, S5, S6, S7). The most significant finding was that many of the genes with increased expression in LCL lesions were involved in immune system activities, such as cell-to-cell signaling and interaction, immune cell trafficking, and inflammatory response. These groupings shared common events and involved 86 genes up-regulated in LCL samples (Tables 1 and 2, Figure 3). Immune response associated genes were related to the recruitment and activation of lymphocytes, granulocytes, natural killer cells and antigen presenting cells. In addition, we detected higher expression of genes related to the control of inflammation and induction of immune tolerance. Conversely, ML samples had only 9 genes with increased expression that were related to inflammation and immunity; all of them were related to recruitment of memory T cells and B lymphocytes (Table 2). The presence and implications of the most important individual genes and pathways is explored in more details in the discussion. The majority of genes with increased expression in ML samples are involved in biological activities that would not have been intuitively predicted for this infectious disease. These biological groupings included cancer, genetic disorders, metabolic diseases, and cell death. Some of these activities were also observed in the LCL group, however, with less significance according to the IPA assay (Table 1, Table S3). We hypothesized that transcriptional differences in the primary infectious lesions of cutaneous leishmaniasis could be used to predict the long term clinical outcome of the infection. To determine if the top differentially expressed genes could be used to predict the clinical form (ML vs LCL), we performed a Leave-One-Out cross validation analysis on the samples [15], [16]. Based on the groups of 5 ML subjects and 5 LCL subjects, the prediction algorithm determined that the top 13 differentially expressed genes showed the best prediction of the clinical phenotype (P = 0. 12; Figure S2). According to the prediction, increased expression of the genes FBXL2, ZDBF2, MORN2, C1orf109, TSPAN13, MRPS31 and STAMBP, and decreased expression of CYB5R3, CREB5, EPB49, STMN3, ADCK1 and TMEM239 in a primary cutaneous lesion would indicate a higher probability of ML development (Table 3). Partly due to the small number of samples, the result was not statistically significant. However, this type of analysis could be effective when using a larger number of clinical samples. Additionally, Ingenuity Pathway Analysis offers a biomarker filter tool based on a database of genes with known biological importance for diagnosis and prognosis for different diseases. Considering both FPKM values and the fold change for biomarker filtering, the IPA assay indicated 15 genes as possible markers of high risk to develop ML (Table 3). Based on the IPA analysis, patients that are likely to develop ML would demonstrate increased expression of the genes SOX9, BRCA2, ARG2, EYA2, PPARG, KLK4, AGR2, ABCG2 and FGF1 and decreased expression of TERT, NTF3, VWA2, TNFRSF6B, ALOX15 and FCER2.
  Although several studies have addressed the immunological differences between LCL and ML patients, most of the observations were performed in situ using mucosal lesions, or systemically using PBMCs from patients with the active disease. The innovation brought by the present work is the comparison of LCL and ML groups starting from the primary cutaneous lesion, which is the beginning of both clinical forms. For this analysis, we used the RNA-seq technology to obtain and evaluate the transcriptome of those lesions. After Glucantime therapy, all patients successfully resolved the initial cutaneous lesions. Despite similar clinical control of the infection, we were able to detect consistent differences in the expression of genes involved in inflammatory and immune responses before the treatment, which could explain the persistence of the infection and ultimately predict the clinical outcome. These findings need to be validated in additional studies with more primary lesions in new patients. Previous studies cannot be directly compared to this analysis mostly because of the different specimens used. However, for discussion, they are reviewed here. Our group had previously shown that mucosal lesions resulting from L. braziliensis infection displayed high levels of transcripts for IL-5, IL-10 and, specially, IL-4 when compared to cutaneous lesions [18]. Other authors observed the presence of increased number of cells producing IL-10 and IFN-γ in these type of lesions [19]. Nevertheless, the systemic approaches revealed that ML patients display higher levels of inflammatory cytokines, such as TNF-α and IFN-γ, and lower IL-10 production as compared with CL patients, besides a weak modulation by IL-10 and TGF-β [20]–[23]. Those studies also concluded that mucosal leishmaniasis is rather an immunopathologic hyperimmune reaction to parasites and their antigens [24]. Our findings revealed that several genes involved in diverse immunological activities, ranging from the innate immune response to adaptive immunity were expressed at lower levels in the primary cutaneous lesion from patients that afterward evolved to ML when compared to those who only presented LCL. This data set suggests that the quality of the immune response against L. braziliensis infection that best correlates with long term control of the disease involves the recruitment and activation of lymphocytes, dendritic cells (DC), macrophages, neutrophils, and natural killer (NK) cells. The specific transcripts detected do not guarantee the presence of specific cells or a predicted role for them. However, the high level of transcripts in LCL samples for genes involved in antigen presentation to both CD4+ and CD8+ T cells suggests this may be an important step for the development of a durable adaptive immune response against Leishmania. The genes TREM-2 and WAS have been associated with the migration of DCs to lymph nodes for antigen presentation [25], [26]. The detection of HLA-DRB1, HLA-DQA1, HLA-DQB1 and CD74 transcripts suggests a significant presence of professional antigen presenting cells (APC) in situ. The participation of HLA has been explored before in the immunopathology of leishmaniasis. The high frequency of some HLA-DR alleles was correlated to the susceptibility to mucosal leishmaniasis [11], while the levels of HLA-DR in cutaneous lesions were down-modulated after drug therapy [27]. Additionally, antigen presentation through MHC class I is supported by the high mRNA levels of endoplasmic reticulum aminopeptidase 1 (ERAP1), Tapasin (TAPBP) and HLA-C [28]–[30]. The transcriptome in LCL primary lesions also indicates a multi-polarized adaptive immune response. Expression of different chemokines, cytokines and receptors, transcription factor and others genes, suggests that the lymphocytes present in the lesions are Th1, Th2 and Th17, CD8+ cells, memory cells CD4+ and CD8+, natural killer T (NKT) cells. Some of these genes had been previously associated with leishmaniasis. IL-17 is detected in high levels in subclinical samples as compared to CL samples, suggesting the participation of Th17 in the control of the infection [31]. IL-21R, CD30 (TNRFSF8) and WAS participate of the Th2 polarization during Leishmania infection [32]–[34]. CXCR3 plays a critical role in the host defense against cutaneous leishmaniasis caused by L. major and helps the Th1 response [35], while CCL17 is associated with diffuse leishmaniasis and late LCL cases [35], [36]. GATA-3 favors the resistance to L. panamensis infection, even though promotes the secretion of Th2 cytokines [37]. TNFRSF18 stimulation improves the therapeutic drug response in visceral leishmaniasis [38] and CD23 (FcεRII) participates of the nitric oxide control of L. braziliensis infection [39]. We also observed high levels of anti-inflammatory genes in LCL lesions, such as ALOX15, BCL3, LRP1 and SOD3. These genes act by controlling migration and proliferation of inflammatory cells, vascular permeability, response to inflammatory factors, and production of pro-inflammatory cytokines [40]–[43]. We also observed genes involved in the induction of immune tolerance, but none of them had been reported in leishmaniasis. Controversially, compared with LCL samples, ML samples showed a small number of genes with higher levels of expression related to inflammation and immunity. However, two potentially important genes are ICOS and PPARG. ICOS regulates both Th1 and Th2 responses and is involved in the clearance of L. mexicana infection [44]. PPARG was previously associated with susceptibility to leishmaniasis. Activation of PPARγ in macrophages not only induces differentiation of macrophages into the non-inflammatory M2 type but also reduces the resistance to L. major infection [45]. Additionally, the presence of increased CXCL13 indicates recruitment of B cells for antibody production [46] and TSLP indicates activation of B lymphocytes, dendritic cells, and Th2 cell differentiation [47]. ABCG2 is a multidrug resistance protein and is involved in the development of Langerhans cells [48]. Taken together, our results using RNA-seq technology and a transcriptomic approach suggest that differences in the quality of the initial immune response of an individual to L. braziliensis infection may influence the long term control of the infection. Our observations of the host transcriptome in primary lesions indicate that individuals prone to develop ML may have an insufficient or maladaptive immune response in the cutaneous stage of infection. Such delay in developing specific immunity could be partially responsible for the exacerbated immune response observed in the patients with mucosal leishmaniasis. However, our findings are strictly observational and do not allow us predict the mechanistic underpinnings of this difference. Do ML patients develop a tolerance to Leishmania antigens after infection or do they have an intrinsic immunological defect that affects the parasite control? The detection of a significant number of genes that segregate consistently with the clinical outcomes suggests that transcriptional fingerprints to predict the clinical development of ML may have a functional prognostic significance. Different approaches can be used to answer these questions. We used leave-one-out cross validation analysis to determine the predictive or prognosticating power of the top differently expressed genes. We also used the database offered by IPA software that enumerates genes with previously reported medical importance for prognosis. This literature bias may provide an important reason to study these genes in more detail, but does not directly impact the interpretation of the leave-one-out cross validation results. While we would not expect the resulting lists of genes to overlap (as they did not), they may all have important predictive characteristics when taken together, but used independently. These results encourage us to increase the number of samples to confirm the biological and medical relevance of these results. A larger prospective study to determine positive or negative predictive values and the sensitivity and specificity of these approaches are also needed before clinical use in patients. In summary, this study provides proof of principle to use next-generation sequencing of the host transcriptome to identify potentially important biological pathways or drug targets in the host-response to L. braziliensis infections. Although pathogen transcripts were detected at a low level in the present study, the ability to simultaneously analyze the host and pathogen mRNAs may provide a valuable tool for fundamental pathogenesis research of leishmaniasis or other organisms in the future. Finally, this study suggests that the use of the transcriptional fingerprint or gene expression signature to predict disease outcome is a viable approach. As the costs and time required for next-generation sequencing continue to decrease, this approach may quickly surpass current laboratory tests for both diagnosis and prognosis of many infectious processes.
summary: In Brazil, American tegumentary leishmaniasis is mainly caused by Leishmania braziliensis infection. Usually, it begins as a solitary skin lesion, which is called Localized Cutaneous Leishmaniasis (LCL). However, after this lesion heals, 5% of the patients may develop destructive lesions of the mucosa of nose and throat, which is called Mucosal Leishmaniasis (ML). Currently, there is no technology to identify individuals at risk for ML, and the factors that control the evolution to ML remain unknown. This work aims to study the human gene expression patterns that may contribute to the clinical manifestation of the disease. We used the RNA-Seq technique to study skin lesions from individuals that had LCL (LCL group) and those who developed ML (ML group). Our results suggest that individuals that progressed to ML expressed low levels of genes involved in the immune and inflammatory responses, which might lead to insufficient control of the infection. We were also able to identify a potential gene expression signature to predict long-term disease outcome.
keywords: immunopathology
  genome expression analysis
  genomics
  gene expression
  molecular genetics
  biology
  computational biology
  immunology
  microbiology
  host-pathogen interaction
  parasitology
title: Transcriptome Patterns from Primary Cutaneous Leishmania braziliensis Infections Associate with Eventual Development of Mucosal Disease in Humans