Association of HLA class I and II genes with cutaneous leishmaniasis: a case control study from Sri Lanka and a systematic review

Background The outcome of leishmaniasis is an interplay between Leishamania and the host. Identifying contributory host genetic factors is complicated by the variability in phenotype, ethnicity and parasite species. Leishmaniasis is caused exclusively by L. donovani in Sri Lanka with localized cutaneous leishmaniasis (LCL) being the predominant form. We report here an association study of human leucocyte antigen (HLA) class I and II genes with LCL in Sri Lanka, the first on HLA associations in cutaneous leishmaniasis in a South Asian population. Methods An existing DNA repository of 200 each of patients and controls was typed for HLA-DQ by PCR-SSP. Next generation sequencing-based typing for HLA-A, HLA-B and HLA-DRB1 alleles was done in a subset of 280 samples. Association tests were performed on 28,489 genotyped and imputed SNPs spanning a region of 1.4 Mb across the HLA genes. To compare our results with similar studies, we carried out a systematic review to document all HLA associations reported to-date for cutaneous and muco-cutaneous leishmaniasis. Results DRB1*04 DQB1*02 (P = 0.03; Pc = 0.09), DRB1*07 DQB1*02 (P = 0.03; Pc = 0.09) haplotypes were absent in patients. B*07 (P = 0.007; Pc = 0.13; OR = 0.36; 95 % CI = 0.17–0.77) allele and DRB1*15 DQB1*06 (P = 0.00; Pc < 0.01; OR = 0.3; 95 % CI = 0.2–.0.6) haplotype were over represented in controls and DRB1*15 (P = 0.002; Pc = 0.01) allele was over represented in patients. Two SNPs (rs281864595/rs1050517) in the antigen recognition region of HLA-B, comprised a haplotype more frequent in controls (P = 0.04). The alleles identified by the systematic review to predispose or to protect from cutaneous/mucocutaneous leishmaniasis remained highly heterogeneous in different populations studied. Conclusions Our preliminary findings suggest a role for some class I and class II HLA genes in determining predisposition to LCL in this population which should be corroborated with further studies. The systematic review reiterates this need, as the purported susceptibility or protection gained by certain HLA alleles or haplotypes has rarely been independently verified. Electronic supplementary material The online version of this article (doi:10.1186/s12879-016-1626-8) contains supplementary material, which is available to authorized users.


Background
Leishmaniasis is an anthropo zoonotic disease with a global incidence of approximately two million cases [1]. Caused by a protozoan parasite of the genus Leishmania and transmitted by the sandfly, this is a disease with a spectrum of phenotypes ranging from self-healing localized skin lesions to potentially fatal visceral disease. Classically, different species of Leishmania are implicated in the distinct forms of the disease.
Cutaneous disease, even if not fatal can be chronic, disfiguring and stigmatizing and bears the brunt of the disease burden with 1.5 million new cases reported every year [2]. The cutaneous form itself has a range of clinical presentations with localised cutaneous leishmaniasis (LCL) being the commonest. The severity of cutaneous disease depends on both the extent of parasite replication and the relative induction of immunopathologic responses by the host. The genetic background is considered the single most important host factor which modulates these immune responses [3][4][5][6].
The human major histocompatibility complex (MHC) gene cluster on the short arm of chromosome 6 has been central to many studies investigating host genetic factors as determinants of outcome of leishmaniasis. The classical human leukocyte antigen (HLA) genes of class I and II located within the MHC, encode molecules identified for their role in presentation of antigens to CD8 + and CD4 + T cells respectively. CD4 + T cells are the major source of interferon (IFN)-γ which controls parasite multiplication during the early phases of Leishmania infection [7,8] whereas CD8 + T cells have been described to contribute to IFN-γ production and the differentiation of Th1 responses [9,10]. Furthermore, localised cutaneous lesions due to L. (V.) braziliensis or L. (L.) amazonensis which demonstrated a balanced immune response and healing, have been shown to have a high prevalence of MHC class I restricted CD8+ T cells by immunocytochemistry analysis [11,12]. Hence, variations at HLA class I and II loci, with resultant variability in interactions with parasite antigen, are likely to affect susceptibility to leishmaniasis [4].
Knowledge on the role of HLA genes as determinants of cutaneous leishmaniasis has mostly originated from studies in South American populations, with the predominant causative parasite species implicated in these infections being L. braziliensis [6,[13][14][15][16]. It has been shown that several class I (HLA-B) and class II (HLA-DQ) alleles are associated with susceptibility or protection to LCL in Venezuelan patients [13]. In a study undertaken on Mexican Mestizos, DRB1*0407, DPA1*0401 and DPB1*0101 were reported to be risk alleles while DPB1*0401 and DR2 were protective towards LCL [15]. More recently, a study conducted in southern Brazil reported a trend towards susceptibility to cutaneous lesions for several class 1 (HLA-B) and class II (HLA-DRB1) alleles while another HLA-B allele tended to be protective [16]. A single study on patients of South East Asian origin with cutaneous leishmaniasis due to Leishmania (V.) guyanensis reported a decrease of HLA-Cw7 antigen [17].
In contrast, the focus on South Asian populations has mainly been on the visceral form of the disease [18]. Leishmaniasis in Sri Lanka, until recently has been exclusively of the localised cutaneous type with over 2500 cases being diagnosed during the last 10 years [19] with numbers reported to the health system considered to be largely under-represented [20]. The causative parasite of both cutaneous and visceral disease has been identified as L. donovani MON 37; a zymodeme very closely related to L. donovani MON 2 which causes visceral disease in neighbouring India [21,22]. However, LCL continues to be the predominant phenotype in the island in both known and new foci [19,23,24].
Identifying human genetic factors which operate in diverse backgrounds of parasite species and host ethnicity would contribute to understanding the complex pathogenic mechanisms underlying leishmaniasis and help identify common molecular pathways for treatment and control. We have earlier reported the observation that LCL in Sri Lanka appears to more commonly affect those of Sinhalese ethnicity [25]. A previous study by us which investigated selected polymorphisms of TNF, LTA and SLC11A1 genes did not show any association with LCL in this population [26]. The aim of this case control study was to investigate the association of HLA class I/II loci with LCL in the local Sri Lankan population. We also carried out a systematic review of published literature to document all similar associations reported in other geographical regions and compared it with our findings.

Study population and DNA samples
The study was carried out using an already established DNA resource for conducting host genetic studies on leishmaniasis in Sri Lanka. Briefly, the study population comprised 200 unrelated patients with a laboratory confirmed diagnosis of LCL and 200 healthy, unrelated controls matched for age, gender and ethnicity. The patients were recruited at the referral clinic conducted by the Centre for Training, Research and Diagnosis of Leishmaniasis of University of Colombo and were from 11 districts of Sri Lanka; with the highest numbers from Matara and Hambantota districts in the Southern province followed by Anuradhapura district in the North Central province. Controls were recruited from the same endemic regions during the same time period as the patients. Controls were excluded if they had a previous history or diagnostic evaluation suggestive of LCL. All subjects were Sinhalese. Genomic DNA was extracted from peripheral blood lymphocytes using QIAamp DNA Mini Kit (Qiagen, USA) according to the manufacturer's instructions. Recruitment of participants and establishment of the DNA resource has been described before [25].
HLA typing HLA-DQB1 alleles were determined by a PCR sequencespecific priming (SSP) technique using AllSet Gold SSP low resolution kit according to the specifications of the manufacturer (Invitrogen, USA). PCR products were analyzed on a 2.0 % agarose gel stained with ethidium bromide (0.5 g/ml) for specific amplification patterns and alleles were assigned with SSP UniMatch (v5.1) software. Genotyping was performed at the Human Genetics Unit, Faculty of Medicine, University of Colombo. Exons 2 to 3 of HLA-A and HLA-B genes and Exons 2 to 4 of HLA-DRB1 gene were sequenced in a subset of 140 each of patients and controls, on an Ion Proton next generation sequencing platform at a commercial sequencing laboratory [27]. The HLA alleles were analyzed and assigned at the same facility, by aligning against reference sequences of IMGT/HLA database and using an in house protocol adapted from published methods [28].

Construction of a SNP panel of the HLA region
The variant calls derived from sequencing reads were used to construct a discovery SNP panel for cases and controls. All bi allelic SNPs identified in any of the subjects were extracted and included in the panel. Each individual sequence was analysed for the presence of the SNPs in the panel and zygosity at each SNP position was determined accordingly using custom built scripts. Further to this, variants within the HLA region on chromosome 6 from 6_29910603 to 6_31324183 were imputed in our discovery panel using IMPUTE2 (version2.3.1) [29] and an integrated dataset from 1000 Genomes Phase 1 [30,31] as the reference panel.

Statistical analysis
HLA allele and haplotype frequencies were estimated by the maximum likelihood method, using the Arlequin (version 3.5) software [32]. Exact tests of Hardy-Weinberg equilibrium for each locus were performed using the same. The frequencies of multi locus haplotypes were estimated under an unknown gametic phase.
The distribution of these alleles and haplotypes was compared with Pearson χ 2 test with Yates correction and Fisher exact test with SPSS (version 15.0). The resulting P values were corrected for multiple comparisons (Pc) by the false discovery rate method [33] and significance was determined at P < 0.05. Both adjusted and unadjusted P values are presented. All statistical analyses were performed at the two digit-allele group level.
Comparison of SNPs between cases and controls was carried out using PLINK [34,35]. SNPs were tested at allelic level and genotype level under dominant, additive and recessive models of inheritance. P value of 0.05 adjusted for multiple testing by FDR method was considered as significant. Haplotypes were constructed for selected SNPs and also compared between cases and controls using the same software. Functional prediction of the SNPs was carried out using standard bioinformatic tools including SNPNexus [36], Mutation Taster [37] and SIFT [38] and further analysed with PDB, Structure and InterPro databases [39][40][41]. All genomic positions reported correspond to NCBI SNP build 137.

Systematic review
The aim of the systematic review was to document all HLA associations reported in literature that increased susceptibility to cuteanous/mucocutaneous leishmaniasis and to correlate these findings with our results.
We searched PUBMED, Web of Science (SCI), EMBASE and CINAHL for research on HLA associations with different clinical types of leishmaniasis. All databases were searched with keywords "Leishmania*" AND "HLA", "Leishmania*" AND "Histocompatibility" in title, keywords and abstracts. All relevant studies on human populations (or case control/family studies) were selected for further review without time or language  Figure S1).

Results
The present study investigated a potential association between LCL and class I and II HLA loci in a South Asian population by comparing the distribution of HLA-A, HLA-B, HLA-DRB1 and HLA-DQB1 alleles in patients and healthy controls. The analysis was complemented by comparing the distribution of SNPs in the above HLA gene regions between the two groups.

Clinical characteristics of patients
A majority of patients presented with painless, single and dry ulcers which had been present for less than 6 months. While 20 % of patients had multiple lesions, 5.5 % of patients had other affected family members. The characteristics of the patients have been described in detail previously [25] and demographic details are summarized in Table 1.

Distribution of HLA alleles
The frequencies of all identified HLA (A, B, DRB1 and DQB1) alleles in patients with LCL and the control group are listed in Table 2. B*07 followed by B*40 were the two most common HLA-B alleles in both patients and controls out of 10 alleles which were common to both groups. At the DRB1 locus DRB1*15 was the only allele seen in patients while 7 alleles were seen in controls with DRB1*15 being present in the majority. The   interpretations cannot be made due to limited sample size = 0.01) that was over represented in patients. Overall, successful allele identification could not be performed for many samples in both control and patient groups due to some regions of DNA in the samples being refractory to amplification. Especially the number of alleles identified for HLA-A was limited and therefore it is not possible to draw conclusions about HLA-A allele distribution in cases and controls.

Distribution of HLA haplotypes
The distribution of specific 2 locus haplotypes is shown in Table 3. The most common haplotype was DRB1*15-DQB1*06 in both patients and controls. Haplotypes

Association analysis of SNPs in HLA region
The discovery SNP panel consisted of 28,489 SNPs spanning a region of 1.4 Mb on chromosome 6 extending from 6:29910603 to 6:31324183. After quality control 14 SNPs in 92 cases and 106 controls remained for analysis. None of the SNPs were associated with LCL at allelic level or genotypic level under different models of inheritance. At haplotype level, a haplotype consisting of alleles of two SNPs at 6:31324638 (rs281864595) and 6:31324643 (rs1050517) was seen to be higher in the controls (84.0 % vs; 93.3 %; P = 0.04). Annotation details of the SNPs which constituted the haplotype which tended to confer protection are given in Table 4.
Structural analysis in relation to above, showed the amino acids implicated by above SNPs to be situated in the alpha 1 -alpha 2 domain of HLA*B which comprises its peptide binding region (Additional file 2: Figure S2).

Systematic review
The initial search with specified search methods yielded 601 abstracts. After removing duplicates and selecting studies that had reported on associations of HLA allele frequency in case control studies, 11 articles remained. Further four articles were excluded as they discussed visceral leishmaniasis. We have summarized the findings of these studies in Table 5. Overall, the reported alleles that purportedly increased or reduced susceptibility to cutaneous/mucocutaneous leishmaniasis remained highly heterogeneous in different populations studied.

Discussion
We studied the association of variations at HLA class I (HLA-A, HLA-B) and HLA class II (HLA-DRB1, HLA-DQB1) loci with susceptibility to localised cutaneous leishmaniasis in Sinhalese in Sri Lanka. The results presented in this report are, to our knowledge, the first on HLA associations in cutaneous leishmaniasis in a South Asian population.
The distribution of HLA alleles showed 7 alleles of HLA-A, 18 alleles of HLA-B and 4 alleles of HLA-DQB1 to be present in patients with LCL. Notably, the HLA-B locus was more polymorphic in the patients than in the control population. Eventhough only HLA-DRB1*15 allele was observed in the patients, this is likely to be due to the limited sample number, as evident by the polymorphic nature of the locus in the larger control group. The overall distribution of alleles in the controls conforms to previous reports from the country [42] and adds to the limited data available on HLA polymorphisms in this population.
The case control analysis did not show any statistically significant associations between HLA markers and LCL after correcting for multiple testing. However, comparisons of trends observed in this study with findings from other populations are noteworthy. For instance, the allele frequency of HLA-B*07 (P = 0.007; Pc = 0.13; OR = 0.36; 95 % CI = 0.17-0.77) was higher in controls suggesting a protective tendency by this allele. Interestingly HLA Bw22, which is similar to HLA-B*07 in the amino acid sequence at alpha 1-helix [43], has been shown to be a risk factor for cutaneous leishmaniasis due to L.    braziliensis in Venezuelans [13]. The same study reported a protective effect conferred by HLA-B15, whereas the frequency of HLA-B15 was twice that of controls in our patients. Ribas-Silva et al. [16], investigating genetic susceptibility to different clinical manifestations of American cutaneous leishmaniasis in southern Brazil, reported HLA B*35 and B*44 to confer susceptibility and B*45 to confer protection against cutaneous lesions. Both HLA B*35 and B*44 were over represented in our patient group but did not reach significant proportions. At the HLA-DQ locus, DQw3 (DQB1*03) and DQw8 (DQB1*0302) have been associated with susceptibility to LCL [13] whereas the present study showed an equivalent distribution of DQB1*03 between the two groups. In contrast to the current study where DRB1*15 (P = 0.002; Pc = 0.01) appeared to be a risk allele, DR2 serotype which encompasses DRB1*15 and DRB1*16 was found to be protective in Mexican Mestizos [15]. Further, DRB1*13 and DRB1*04 were considered risk factors in this same population and in Brazilians [16] respectively.
The small sample number is a limitation of this study. This may lead to a bias in estimating allele and haplotype frequencies and thus caution should be exercised in extrapolating them to the study population. This may also explain some of the observations made, for instance the patient group being monomorphic at the HLA-DRB1 locus, and calls for replication of findings in a bigger cohort in validating the significant associations involving this locus.
The apparently contradictory results for the same locus in different global populations as we report here, is a common occurrence in genetic association studies. These differences most probably reflect the LD patterns across the gene regions under focus, with the clinical outcome being the concerted effect of several loci, some of which may have not been considered in the particular study. This variability is further compounded by the extensive polymorphism in the HLA genes [44]. However, notably the only genome wide association study to be conducted on leishmaniasis reported risk variants for the visceral phenotype in the HLA-DRB1-HLA-DQA1 region, which were common to population groups across different geographical areas and affected by different parasite species (Fakiola et al. [18]). Further, patients with other affected family members may indicate shared living environment or a true genetic association. Thus, complimenting population based case control analyses such as that undertaken in the present study, with multi-case family based linkage analyses and transmission disequilibrium testing would also strengthen the findings.
The systematic review showed that there are only a handful of studies that have assessed the HLA polymorphism in cases and controls worldwide. Most of these studies are from Latin America and our study in fact is the first from Asia on cutaneous leishmaniasis. There are several studies that have assessed HLA polymorphism in visceral leishmaniasis including a large scale multicentre study in Brazil and India [45], but not for cutaneous leishmaniasis. It is also noted that while there are two studies that explicitly stated no HLA associations for visceral leishmaniasis (out of 4 studies) when compared with controls [46,47], all studies that assessed the same for cutaneous/muco-cutaneous leishmaniasis report otherwise. There can be reporting bias with nonpublication of negative results but even in the two studies that showed a significant association it was just limited to 3 HLA alleles in both studies combined [45,48]. The synthesis of observations (Table 5) shows that the purported susceptibility or protection gained by certain HLA alleles or haplotypes has rarely been independently verified by subsequent studies except for few exceptions (DR2 serotype). There can be several reasons for this. Firstly, the studies are few in number and small in sample size. Secondly, they are spread across different countries/geographical regions that have populations isolated from each other with plausibly different baseline HLA allele frequencies in the population. Thirdly, the Leishmania species (hence antigens) in each of these regions differ which can lead to heterogeneity in host response. Overall, it raises the importance of screening individual populations locally to identify protective and susceptible characteristics of host genome which would be important to design effective vaccines for local populations.

Conclusions
Overall, several of the alleles in the HLA-B region which have been associated with susceptibility or protection towards cutaneous leishmaniasis in other populations, were also shown to differ between patients and controls of our study group. The two SNPs on HLA-B gene, which result in a haplotype found to be more frequent in controls as shown by the complementary SNP analysis, have not been widely reported to aid definitive conclusions. Nevertheless, the fact that these implicate two adjacent amino acids in the alpha1 -alpha2 region of HLA*B which comprises its antigen recognition region, suggests candidate loci to be considered in expanded work in this study group. Further studies in this population, on larger sample numbers conferring better statistical power will be required to support the findings of this study. The findings should also be complemented with higher resolution HLA typing, high density SNP analyses of the HLA region and also functional studies on interactions between host immune response cells and parasite antigens in a HLA restricted setting.

Additional files
Additional file 1: Figure S1. PRISMA flow chart for systematic review. Abbreviations HLA, human leucocyte antigen; LCL, localized cutaneous leishmaniasis; PRISMA, preferred reporting items for systematic reviews and meta-analyses