Single-Nucleotide Polymorphisms in Genes Predisposing to Leprosy in Leprosy Household Contacts in Zhejiang Province, China.

PURPOSE: Genome-wide association studies (GWAS) have identified multiple genetic variants associated with leprosy. To investigate the single and combined associations between single-nucleotide polymorphisms (SNPs) and the development of leprosy, we therefore performed generalized multi-analytical (GMDR) analysis in Chinese leprosy household contacts and constructed a risk prediction model. PATIENTS AND METHODS: This case-control study included 229 leprosy cases and 233 healthy household contacts in Zhejiang province, China. Participants were genotyped for 17 polymorphisms selected from GWAS. The Pearson χ2 test, logistic regression and GMDR analysis were performed to investigate gene-gene interactions and construct a risk prediction model for leprosy.
RESULTS: The genotype and the allele distributions of rs142179458, rs2275606, rs663743 and rs73058713 were significantly different between patients and controls. rs2275606, rs6478108, rs663743 and rs73058713 showed an association after adjusting for sex and age in the logistic regression. A five-way interaction model consisting of rs2058660, rs2275606, rs4720118, rs6478108 and rs780668 was chosen as the optimal model for determining leprosy susceptibility. The model classified 237 (51.3%) into the low-risk group and 225 (48.7%) individuals into the high-risk group. The area under the curve (AUC) of this model was 0.757 (95% CI: 0.712-0.803), and the odds ratio for leprosy between the high- and low-risk groups was 9.733 (95% CI: 6.384-14.960; P<0.001). The sensitivity and specificity of the model were observed to be 74.7% and 76.8%, respectively.
CONCLUSION: Our results suggest that rs2058660, rs2275606, rs4720118, rs6478108 and rs780668, five SNPs with a significant sole effect on leprosy, interact to confer a higher risk for the disease in leprosy household contacts (HHCs).

Introduction

Leprosy, which is caused mainly by Mycobacterium leprae, remains a significant health problem, especially in low- and middle-income countries. China is still one of the countries where the disease remains endemic, mainly in southwestern provinces, including Sichuan, Hunan, Yunnan and Guizhou.1 A total of 521 new cases were detected in China in 2018, with a case detection rate of 0.036 per 100,000 people. The registered cases were 970 by the end of 2018, accounting for a prevalence rate of 0.068 per 100,000 population.2

As an infection disease, an individual’s exposure to M. leprae serves as essential in leprosy developing, so a gradually higher proportion of new cases is found amongst the contacts of known cases.3 However, only a small group of contacts will develop leprosy.4 Genetic factors play important roles in the host’s immunologic status, which can determine the clinical manifestations of leprosy; a prospective cohort study demonstrated that genetic relationship was a relevant and independent risk factor.5 At the population level, GWAS have been performed for better understanding of the human genetic background involved in the eventual outcomes of M. leprae infection, and identified a number of SNPs associated with leprosy genetic predisposition.6–12 Among the studies, SNPs located in RAB32, HIF1A, BATF3, LACC1, CTSB, TNFSF15, CDH18, SLC29A3, DEC1, FLG, NOD2, IL18RAP/IL18R1, NCKIPSD and CARD9 were significantly associated with leprosy. Several association studies have replicated these genes in other independent populations.13–17 Despite GWAS, some other risk genes have been reported or validated in Han Chinese with leprosy.18–22

Generalized multifactor dimensionality reduction (GMDR) was applied to explore the best prediction model of higher order gene–gene interactions. GMDR reduces high-dimensional genetic data to a single dimension, explores an interaction model through cross-validation and calculates score-based statistics of each subject using maximum likelihood estimates to classify subjects into two different groups (either high-risk or low-risk group), and was regarded as an objective analytical tool for evaluating multifactorial impact on an association, allowing adjustment for confounding factors.23,24

In this study, we utilized published genetic risk variants in GWAS with a case–control study and tried to evaluate the role of higher order gene–gene interaction to predict leprosy risk in leprosy HHCs in Chinese Zhejiang province.

Patients and Methods

Study Participants

Two independent sets were enrolled in this study. A total of 245 leprosy patients and 245 healthy leprosy HHCs were enrolled in 2015. People affected by leprosy including those in the course of treatment and cured patients in Zhejiang province, China, were enrolled as the case group. Subjects with any kind of autoimmune-related disease were excluded. The diagnosis of leprosy was confirmed by initial clinical evaluation based on clinical manifestations, slit skin smears and histopathological examinations. Healthy leprosy HHCs from the same geographic region were included in the control group. HHCs were defined as people living under the same roof and sharing food with the patient for at least six months among the past six years and excluded those who refused to provide informed consent, and any person who received treatment for tuberculosis or leprosy within one year. HHCs in our study included relatives (first-, second- and third-degree family members) and genetically unrelated contact individuals (spouses and neighbors). The follow-up set consisted of 939 health contact individuals recruited in 2015 and followed up from 2015 to 2019. Demographic characteristics of the study individuals are listed in Table 1.

Table 1

General Characteristics of the Study Groups

	Leprosy	Controls	P-value	Follow-Up Set
Number	229	233		939
Age (years, mean±SD)	54.98±17.255	61.67±16.868	<0.001	55.17±17.04
Male (%)	166 (52.4)	122 (72.5)	<0.001	461 (49.1)

SNP Selection and Genotyping

We selected 17 SNPs from previously published GWAS and one study combined whole-exome sequencing and targeted next-generation sequencing within the GWAS loci,6–12 wherein was a genome-wide significant association (p<5×10−8) between the SNPs in the genes RAB32, HIF1A, BATF3, LACC1, CTSB, TNFSF15, CDH18, SLC29A3, DEC1, FLG, NOD2, IL18RAP/IL18R1, NCKIPSD and CARD9 and leprosy. This was based on multiplex polymerase chain reaction (PCR) to precisely genotype SNPs with next generation sequencing. After the PCR amplification, the products were genotyped according to the manufacturer's protocol using the Illumina Hiseq X-10 platform. Genotype call rates >95% for both SNPs and individuals, and minor allele frequency (MAF) >1% were found. Ultimately, two variants with MAF <1% were eliminated (rs149308743 and rs145562243). A total of 15 variants and 462 subjects (229 participants with leprosy, 233 controls free of leprosy) were included in the case–control analysis.

Statistical Analysis

The Pearson χ2 test was used to analyze the differences in gender between the patient group and the control group, and the two-sample t-test was used to compare age by SPSS 23.0 software. The P value was considered to be statistically significant at the 0.05 level.

The Pearson χ2 test or Fisher test using PLINK v 1.07 were used to analyze the differences in distributions of genotypes and alleles between the leprosy cases and the control group. We also tested the associations between phenotypes and SNPs based on a logistic regression model adjusted for age and gender, and we constructed a weighted genetic risk score (wGRS) for each individual by summing the number of risk alleles (0/1/2) for each of the four SNPs with P <0.05 weighted by the β coefficient of each allele obtained from the logistic regression. Two-tailed P values of 0.0033 or less after Bonferroni multiple testing correction were considered statistically significant.

GMDR analysis with 10-fold cross-validation was done for assessing gene–gene interactions as an extension of the multifactor dimensionality reduction method.25 It is applicable to both dichotomous and quantitative phenotypes that allow adjustment for covariates and is more precise and accurate than multifactor dimensionality reduction. The GMDR models are evaluated on the testing balanced accuracy (TBA), the cross-validation consistency (CVC) and the statistical significance. The results of GMDR analysis were confirmed by χ2 tests. The AUC and the 95% confidence intervals (CI) were calculated for the risk model with SPSS (version 23; IBM Corporation).

Results

Baseline Characteristics of Participations

There were no subsequent cases in the period of 2015 to 2019 in the follow-up set. 229 participants with leprosy and 233 controls free of leprosy who met the inclusion criteria were included in the study. The case group consisted of 229 leprosy patients, with a mean age of 61.67±16.868 years,while 233 healthy leprosy contacts consisted the control group, with a mean age of 54.98±17.255 years. 939 healthy leprosy HHCs consisted the follow-up set, with a mean age of 55.17±17.044 years. The baseline characteristics of participants are listed in Table 1. Significant difference was found in mean age and sex ratio between leprosy patients and controls. The individuals in the case group were more likely to be older and male.

Association Between Leprosy and SNPs

Genotype distributions between leprosy cases and controls were analyzed using the χ2 test or Fisher test. shows the distributions of genotypes and alleles for the SNPs. The genotype and the allele distributions of rs142179458, rs2275606, rs663743 and rs73058713 were significantly different between patients and controls. Significant difference between leprosy and control groups was found in rs663743 (P=0.00299, after Bonferroni multiple testing correction).

Four variants showed an association at P<0.05 after adjusting for sex and age in the logistic regression. These included rs2275606 at the RAB32 locus (P=0.03), rs6478108 at the TNFSF15 locus (P=0.03), rs663743 at the CCDC88B locus (P=0.0002) and rs73058713 at the CDH18 locus (P=0.005). The characteristics and association results of the 15 variants are displayed in Table 2.

Table 2

Association Between SNPs and Leprosya

SNP	Region	Candidate Gene	Risk Allele	OR (95% CI)	P-value
rs2275606	6q24.3	RAB32	A	1.41 (1.03–1.93)	0.03
rs142179458	14q23.2	HIF1A	A	2.63 (0.98–7.07)	0.06
rs2221593	1q32.3	BATF3	A	1.24 (0.94–1.65)	0.13
rs3764147	13q14.11	LACC1	G	1.22 (0.92–1.62)	0.16
rs55894533	8p23.1	CTSB	C	1.30 (0.98–1.71)	0.07
rs6478108	9q32	TNFSF15	C	1.35 (1.03–1.78)	0.03
rs73058713	5p14.3	CDH18	A	1.81 (1.20–2.73)	0.005
rs780668	10q22.1	SLC29A3	A	1.24 (0.95–1.60)	0.11
rs10817758	9q32	DEC1	C	1.09 (0.83–1.43)	0.53
rs146466242	1q21.3	FLG	T	1.11 (0.57–2.16)	0.76
rs9302752	16q12.1	NOD2	C	1.20 (0.89–1.62)	0.23
rs2058660	2q12.1	IL18RAP/IL18R1	T	1.11 (0.84–1.45)	0.47
rs4720118	7p14.3	BBS9	C	1.02 (0.76–1.35)	0.91
rs663743	11q13.1	CCDC88B	A	2.04 (1.40–2.97)	0.0002
rs76418789	1p31.3	IL23R	A	1.42 (0.80–2.53)	0.23

Note: aAssociation tested with logistic regression model adjusted for sex and age.

Abbreviations: OR, odds ratio; CI, confidence interval; SNP, single-nucleotide polymorphism.

Interaction Models by GMDR Analysis

The higher order gene–gene interactions were analyzed using the GMDR software with the adjustment for sex and age. The GMDR analysis yielded the best models for all possible one- to eight-locus models (Table 3). Based on the TBA, CVC and significant P values, the five-way interaction model (rs2058660, rs2275606, rs4720118, rs6478108 and rs780668) was regarded as the optimal model. The testing balanced accuracy was 0.5889, the testing accuracy was 0.5905, the testing sensitivity was 0.6036, the testing specificity was 0.5743, the testing odds ratio was 2.4807 (95% CI: 0.3472–17.7234). For the whole dataset statistics, the training accuracy was 0.7904, the training sensitivity was 0.7879, the training specificity was 0.7930, the training odds ratio was 14.2280 (95% CI: 7.3547–27.5247). There were 76 genotype combinations of the five-way interaction model having high risk of leprosy susceptibility, and the genotype combination (TC, AG, CC, TC, AG) especially showed highest risk for leprosy (Figure 1).

Table 3

Gene–Gene Interaction Models for Leprosy by GMDR Analysis

Interacting SNPs	TBA	CVC	P-value
rs663743	0.5323	8/10	0.1719
rs663743×rs73058713	0.5585	5/10	0.1719
rs2275606×rs4720118×rs663743	0.5768	9/10	0.0547
rs2058660×rs2275606×rs6478108×rs780668	0.4924	2/10	0.6230
rs2058660×rs2275606×rs4720118×rs6478108×rs780668	0.5889	8/10	0.0107
rs10817758×rs2058660×rs2275606×rs4720118×rs6478108×rs780668	0.6286	6/10	0.0010
rs10817758×rs2058660×rs2275606×rs3764147×rs4720118×rs6478108×rs780668	0.6049	4/10	0.0107
rs2058660×rs2221593×rs2275606×rs3764147×rs4720118×rs55894533×rs6478108×rs780668	0.5842	4/10	0.0547
rs10817758×rs2058660×rs2221593×rs2275606×rs3764147×rs4720118×rs55894533×rs6478108×rs9302752	0.6927	4/10	0.0010

Abbreviations: GMDR, generalized multifactor dimensionality reduction; SNP, single-nucleotide polymorphism; TBA, testing balanced accuracy; CVC, cross-validation consistency.

Figure 1

Results of adjusted GMDR analyses for the five-factor model. Dark gray and light gray boxes correspond to the high- and low-risk genotype combinations, respectively. The top number above each bar is the sum of scores for the corresponding group of individuals (cases or controls with particular five-locus genotype). The heights of the bars are proportional to the sum of scores in each group.

Subsequently, the five-locus genotypes were classified into high-risk or low-risk due to the GMDR results and we combined 77 low-risk genotypes and 76 high-risk genotypes to perform the χ2 test to evaluate the five-locus interaction model. The model classified 237 (51.3%) into the low-risk group and 225 (48.7%) individuals into the high-risk group. 171 of 229 leprosy cases were categorized into the high-risk group while 179 of 233 healthy HHCs were categorized into the low-risk group by the model. The estimated AUC of every single risk SNP with P<0.05 showed poor predictive ability. wGRS did improve the predictive ability of the model. However, the GMDR model leads to a much greater AUC increase (Table 4). The AUC of the GMDR model was 0.757 (95% CI: 0.712–0.803), and the odds ratio for leprosy between the high- and low-risk groups was 9.733 (95% CI: 6.384–14.960; P<0.001). The sensitivity and specificity of the model were observed to be 74.7% and 76.8%, respectively. We further evaluated the effectiveness of the prediction by calculating the number of high-risk subjects with this model in the follow-up set. Among the 939 healthy leprosy HHCs in the follow-up group, 375 (39.9%) contact individuals were categorized into the high-risk group and 481 (51.2%) subjects were categorized into the low-risk group by the model, which implies we should pay more attention to the higher risk contacts with longer follow-up.

Table 4

Performance of Risk Models for Leprosy

Risk Models	AUC (95% CI)
rs2275606	0.542 (0.489–0.594)
rs6478108	0.500 (0.447–0.553)
rs663743	0.568 (0.516–0.620)
rs73058713	0.559 (0.507–0.612)
wGRS	0.625 (0.574–0.676)
GMDR	0.757 (0.712–0.803)

Abbreviations: GMDR, generalized multifactor dimensionality reduction; AUC, area under the curve.

Discussion

As a millenary disease, tracing and strengthening active case finding targeted at higher-risk contact individuals of leprosy serves as a key point of interruption in the leprosy transmission. With the advances in molecular biology, numerous studies have begun to assess the contribution of specific genes and SNPs to the development of the leprosy phenotype and supported that human genetics plays an important role in leprosy susceptibility, and different sets of genes modified host predisposition to leprosy per se, its clinical forms and reactional states.26 In view of the diversity of these genes, we chose 17 independent SNP loci associated with leprosy based on GWAS and one candidate-gene study to determine the role of these genes in the development of leprosy in the Chinese population of Zhejiang province.

Our study highlights SNPs of genes encoding proteins involved in the immune response to the M. leprae infection as risk factors for developing leprosy. The genotype and the allele distributions of rs142179458, rs2275606, rs663743 and rs73058713 were significantly different between leprosy cases and HHCs. Bonferroni adjustment method was performed to correct multiple comparisons, and the allele distribution of rs663743 was found to be significant. Meanwhile, we performed logistic regression analysis to determine the role of genetic factors after adjusting the confounding factor. Four variants showed significant association, including rs2275606, rs6478108, rs663743 and rs73058713.

Similar to our findings, a previous study shows that rs142179458 in HIF1A contributes to leprosy risk, and the HIF1A mRNA expression was increased during M. leprae infection.7 rs2275606 located at the RAB32 locus, our research result suggesting that the Rab32-dependent pathway may act as a host defense pathway in humans. In addition, other studies showed that Rab GTPase Rab32 emerged as a key regulator of a host defense pathway that can restrict intracellular bacterial pathogens. Rab32 and its guanine nucleotide exchange factor BLOC-3 are essential to prevent the growth of the human-restricted Salmonella enterica serovar Typhi (S. Typhi) in mice.7 rs663743 at the CCDC88B locus shows significant association with leprosy in our study, and CCDC88B was reported associated with other autoimmunity and inflammatory diseases,29,30 and was required for regulating movement and migration of dendritic cells.31 rs73058713 at the CDH18 locus was found to be associated with leprosy in our study, and a previous study suspects that the CDH18 gene may be involved in neuronal development governing metabolic processes in later life.32 rs6478108 at the TNFSF15 locus was also found to be associated with leprosy in logistic regression analysis, the TNFSF15 gene is involved in mediating the switch from Th1 cells to Th2 cells, and plays an important role in the control of both pro- and anti-inflammatory cytokines.33,34

GMDR was performed as a genetic risk model to determine the influence of the combination of SNPs along with confounding factors.24,35 In our present study, GMDR was used to study the extensive gene–gene interactions of leprosy along with the adjustment of confounding factors. The best gene–gene interaction model was selected across all multi-locus models that maximized CVC, TBA and significant P values for prediction of leprosy. The five-way interaction model (rs2058660, rs2275606, rs4720118, rs6478108 and rs780668) was considered the optimal model. Furthermore, we predicted leprosy risk with this model. This model is considered useful in discriminating between high- and low-risk individuals (AUC=0.757), the AUC is higher than in a previous study,36 single SNP and wGRS models. Compared with the low-risk HHC group, the high-risk group had a 9.733 times higher risk for leprosy. This result demonstrates the considerable value of this model in risk stratification of leprosy developing. We also displayed the clinical effect of this model in the identification of an independent cohort of HHCs at higher risk of developing leprosy, and the higher-risk contact subjects should receive more frequent monitoring and educate them to increase vigilance of developing leprosy.

The strength of the study was that our study subjects were close contacts of leprosy, which leads to a more profound clinical significance compared with the general population. And we investigated 17 SNPs based on GWAS and one candidate-gene study on leprosy risk. There were a few limitations in our study. Firstly, there are some other SNPs in these genes studied in this research that showed a positive association with leprosy (), and it might be necessary to study more polymorphisms associated with leprosy and other epidemiology factors in the future. Moreover, due to the small size of samples in this study, this five-locus model could not cover all genotypes; there were 83 (8.8%) subjects in the follow-up set who could not be categorized because their genotype combinations were not within the scope of the prediction model. And then we were unable to precisely estimate the model’s performance. Thirdly, the HHCs genetically related to index patients were included in the control group, which may have an impact on the results of association analysis, and subjects with or without genetic relationship will need to be studied separately in the future. Finally, further validation studies of the five-locus model in other prospective contact Chinese cohorts are required.

Conclusion

The results presented in our study provide evidence for a five-locus genetic interaction model in identifying higher-risk leprosy contact individuals. This model with good discrimination capacity may be translated into a tool for predicting leprosy developing and risk stratification. Medical workers should pay attention to the higher-risk contact group selected by the model when they make a decision to trace leprosy HHCs.