In total, 1950 unrelated Chinese Han participants, consisting of 636 TB patients and 608 healthy controls in the discovery study as well as 301 TB patients, 201 LTBI subjects and 204 uninfected healthy controls (HC) in the replication study, were consecutively recruited from the West China Hospital between July 2013 and August 2017. For the discovery study, the diagnosis of TB was based upon the following criteria: histopathological evidence of TB disease and/or culture positivity for MTB and/or smear positivity for MTB in at least two separate specimens and/or radiological and clinical findings consistent with TB, with positive clinical response to anti-TB therapy (18). TB cases were divided into two subgroups: (1) pulmonary TB patients (PTB, pathological changes limited to the lung) and (2) extra-pulmonary patients (EPTB, pathological changes involving other tissues or organs merely or in combination with the lungs). The healthy controls who did not have a history or evidence of TB on the basis of their syndromes and radiographic examination results were enrolled during their routine health examination in the West China Hospital. For the replication study, both the uninfected HC group and LTBI participants were close contacts of PTB patients. LTBI was defined as adults (age > 18) with a documented positive result of IGRAs and negative results of radiological and clinical manifestations. HC individuals were defined as adults with negative results of IGRAs, radiological and clinical findings. The inclusion criteria of TB patients enrolled in the replication study were the same as the discovery study. Another 60 healthy volunteers were recruited for comparison of mRNA levels of neutrophils. All participants diagnosed with HIV infection, diabetes mellitus, autoimmune disorders, tumors or treated with immunosuppressive drugs were excluded.

NFE2L2 is located on chromosome 2 (2q31) and has five exons and four introns. TagSNPs were selected from the region 3,000 base pairs upstream to 300 base pairs downstream of the NFE2L2 gene based on the Chinese Han Beijing data of the HapMap database (http://hapmap.ncbi.nlm.nih.gov, HapMap Data Rel 27 Phase II + III, Feb09, on NCBI B36 assembly, dbSNP b126) by Haploview software 4.2. Using criteria of minor allele frequency (MAF) ≥ 5% and the Tagger pairwise method (r² ≥0.8), five tagSNPs of NFE2L2 (rs10497511, rs2364723, rs13005431, rs6726395 and rs1962142), representing eleven SNPs with MAF ≥ 0.05 in the covered gene region, were selected for genotyping.

A total of 5 ml peripheral venous blood was collected from each subject and genomic DNA was extracted using the AxyPrep genomic DNA Mini kits (Axygen, USA) following the manufacturer’s protocol. SNPs were genotyped on the MassArray Analyzer 4 system (Sequenom,USA). All probes and primers are available upon request. 5% samples were genotyped twice for quality control. The genotyping related primers were listed in Supplementary Table 1 .

Polymerase chain reaction (PCR) was performed to amplify a 651bp sequence surrounding rs13005431 from genotyped genomic DNA with the TT genotype of rs13005431. Primers were designed using Primer-BLAST (www.ncbi.nlm.nih.gov/tools/primer-blast/) and sequences for restriction sites NheI and XhoI were introduced ( Supplementary Table 2 ). Both gel purified PCR product and pGL3-promoter reporter vector (Promega, Madison, WI, USA) were digested by NheI and XhoI restriction enzymes (Takara Bio Inc., Kusatsu, Japan) at 37°C for 4 hours and then ligated by T4 DNA ligase (TaKara Bio Inc.) at 16°C for 8 hours. After transforming into competent E. coli DH5α cells and extraction by plasmid miniprep kit (TianGen, Beijing, China), the recombinant plasmid vector was verified by direct sequencing and named PGL3-rs13005431T-promoter ( Supplementary Figure 1 ). A site-directed mutagenesis strategy was applied to acquire another recombinant plasmid vector of PGL3-rs13005431C-promoter (see Supplementary Table 2 for primer sequences) which was confirmed by DNA sequencing.

HEK293T cells were plated into 96-well culture plates at a density of 3.2×10⁴ cells/well and cultured for 48 hours to acquire 90% confluence at the time of transfection. The PGL3-promoter, PGL3-rs13005431T-promoter and PGL3-rs13005431C-promoter were transfected into the HEK293T cells together with PRL-CMV plasmid vectors which acted as an internal reference using Lipofectamine^® 2000 Reagent (Invitrogen (Thermo Fisher Scientific), Waltham, MA USA). Dual-luciferase reporter assays (Promega, USA) were carried out after 30h transfection on a GloMax™ 96 microplate luminometer. The transfection experiment was performed in triplicate. Results are presented as relative luciferase activity by dividing the firefly luciferase activity of each well by the Renilla luciferase activity. The normalized luciferase activity of recombinant plasmids was further normalized to the PGL3-promoter group.

10 ml of peripheral venous blood in a heparin anticoagulant tube was obtained from 60 healthy volunteers. Neutrophils were isolated using Ficoll gradient density centrifugation and confirmed by Wright staining and flow cytometry (CD11b, CD16). The cells were then cultured in 24-well culture plates at a density of 2×10⁶ cells/well for 3 hours, activated or not by GM-CSF (in a final concentration of 1 ng/ml). The Trizol method was then used to extract the total RNA of the neutrophils. Reverse transcription was conducted using the PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara Bio Inc.). Real-time PCR for NFE2L2 was then performed on an ABI7300 Sequence Detection System (Applied Biosystems (Thermo Fisher Scientific), Waltham, MA USA) using SYBR Premix Ex Taq II (Takara Bio Inc.), with B2M as the house-keeping gene (19). Primers for real-time PCR are shown in the Supplementary Table 2 .

Experimental data were analyzed using the Statistical Package for Social Sciences version 17.0 (SPSS, Chicago, IL, USA). Distributions of clinical characteristics between the control group and TB patients were evaluated by the chi-squared test and Fisher exact probability for categorical variables and the student’s t-test for continuous variables. P values, Odds Ratios (OR) and 95% confidence intervals (95%CI) of association between SNPs and TB susceptibility were calculated by binary or multinomial logistic regression under four genetic models (allelic, additive, dominant and recessive). Linkage disequilibrium (LD) analysis and haplotype blocks were created by Haploview software 4.2 (20), while p values, ORs, 95%CIs and the global test of haplotype analysis were calculated using the SHEsis program (21). Genotype and allele frequencies of SNPs in the control group were assessed to determine whether they conformed to Hardy-Weinberg equilibrium (HWE). Gene-environment interaction was detected using the multifactor dimensionality reduction (MDR) constructive induction algorithm. SNP-by-sex additive interactions by the method of Andersson and multiplicative interactions by logistic regression in the discovery study (22). Normalized luciferase activity and relative mRNA levels were evaluated by the student’s t-test. A statistically significant difference was indicated by a two-sided P-value < 0.05.

Table 1 shows the characteristics of all the study subjects. In the discovery study, there was no significant difference in the distribution of gender and age between the case and control groups. However, smoking status was statistically different (P = 0.003). In the replication study, although gender distribution was similar among the three groups, age significantly differed between groups. Smoking status and TB types of the participants in the replication study were not available.

The locations of the five tagSNPs genotyped are presented in Figure 1 . The genotype call rates varied from 99.7% to 100% and the genotype reproducibility was 100%. All SNPs in the control group of both studies were in HWE.

In the discovery study, we analyzed the data between TB patients and healthy controls. As shown in Table 2 , allele C of rs13005431 and allele G of rs2364723 were associated with increased TB risk (P = 0.010 and P = 0.041). In stratified analyses, the deleterious effect of rs13005431 C was still seen in the non-smoker, EPTB and female subgroups (P = 0.003; P = 0.015 and P<0.001, respectively). In addition, rs2364723 G still demonstrated a relationship with increased susceptibility to TB in the non-smoker and female subgroups (P = 0.022 and P = 0.004). In addition, allele A of rs6726395 was demonstrated to be a risk factor associated with TB risk in the female subgroup (P = 0.003). There was no significant association observed between the other two tagSNPs (rs10497511 and rs1962142) and TB susceptibility ( Table 2 ).

In the replication study, we analyzed the data among TB patients, LTBI subjects and uninfected healthy controls using a multinomial logistic regression analysis model. When comparing the LTBI group with TB patients, only allele A of rs6726395 was associated with increased TB risk (P = 0.035) ( Table 3 ). When the HC group was compared with the LTBI group, none of the tagSNPs was associated with susceptibility to TB infection. In addition, we combined LTBI subjects and uninfected healthy controls as a healthy control group. When comparing the healthy control group with TB patients, we still seen that rs6726395 A was associated with increased TB risk (P = 0.028) ( Table 4 ).

Furthermore, we combined the data of the two studies including the discovery study and the replication study. As shown in Table 4 , allele C of rs13005431, allele G of rs2364723, and allele A of rs6726395 were associated with increased TB risk (P = 0.002; P = 0.013 and P = 0.014, respectively). There was no significant association observed between the other two tagSNPs (rs10497511 and rs1962142) and TB susceptibility.

The LD analysis ( Supplementary Figure 2 ) demonstrated mild-to-moderate levels of LD between the five tagSNPs of NFE2L2. For the discovery study, haplotype analyses of the two blocks are shown in Supplementary Table 3 . The frequency of haplotype rs1049751-rs13005431 AC was significantly higher in the TB group (P = 0.013), while the frequency of haplotype AT was higher in the control group (P = 0.025) and there was a significant global result (P = 0.025). For the replication study, no significant haplotype effect was observed.

The results of gene-environment interactions after MDR analysis are shown in Supplementary Table 4 . In the discovery study, smoking-rs13005431 formed the best interaction model with 55.69% testing balanced accuracy and 10/10 cross-validation consistency. Smoking-sex-rs13005431 formed an interaction model with 55.03% testing balanced accuracy and 8/10 cross-validation consistency. After 1000-fold permutation testing, both models were found to be significant (P = 0.001-0.002 and P = 0.023, respectively). As shown in Supplementary Figure 3 , smoking status formed the high-risk models regardless the genotypes of rs13005431, non-smoking status and rs13005431 TT formed the low-risk model, while non-smoking status and rs13005431 CC/TC formed the high-risk models. As shown in Supplementary Table 5 , the only positive additive interaction was observed between female sex and rs13005431 TC+CC/C genotypes. In addition, there were suggestive positive multiplicative interactions of three SNPs (rs13005431, rs2364723 and rs6726395) with female sex under two genetic models.

We predicted the potential functional significance of SNPs from three aspects including protein coding, splicing regulation and transcription regulation ( Supplementary Table 6 ). The influence factors of transcription regulation can be judged by three terms as following: 1) prediction of transcription factor binding sites altering using TRANSFAC or JASPAR database; 2) other factors affecting transcription regulation, such as high sensitive site of DNase I, histone modification, according to the comprehensive evaluation by GLODEN path; 3) judgment of evolutionary conservatism by PhastCons. According to the predicted results, rs13005431 was identified as a candidate SNP for further functional identification.

rs13005431, located in the first intron of NFE2L2, is outside of the promoter region. However, this SNP might disrupt an enhancer or silencer sequence and thus influence promoter activity indirectly. As previously reported, variants in the first intron might influence gene activity (23). Moreover, rs13005431 was reported to be a cis-eQTL (expression quantitative trait locus) of NFE2L2 in whole blood (P = 1.978×10^-4) (24). Therefore, the dual-luciferase reporter assay was used to investigate the regulatory role of rs13005431. As shown in Figure 2 , the normalized luciferase activity of PGL3-rs13005431C-promoter was significantly lower than the PGL3-rs13005431T-promoter and PGL3-promoter constructs (P <0.001 and P <0.001, respectively), indicating that the mutant allele C of rs13005431 may decrease promoter activity. In addition, no significant difference was observed between the PGL3-rs13005431T-promoter group and the PGL3-promoter group (P =0.077), suggesting that the wild type allele T of rs13005431 had no significant effect on promoter activity.

To further confirm whether allele C of rs13005431 influences the promoter activity of the NFE2L2 gene, 60 healthy subjects were enrolled. The relative expression levels of NFE2L2 in neutrophils were analyzed according to different genotypes under basal and GM-CSF activated conditions. As shown in Figure 3A , the mRNA expression levels of NFE2L2 in GM-CSF activated neutrophils were significantly increased compared with the unstimulated group (P <0.001). As shown in Figure 3B , subjects with genotype TC/CC had lower relative expression level than subjects with genotype TT in the GM-CSF activated group (P =0.017) while no difference was found when the neutrophils were unstimulated (P =0.123).