Genome-wide association study of resistance to M. tuberculosis infection in Vietnam

First, we conducted a GWAS in the Vietnamese sample using 185 uninfected and 353 infected subjects, consisting of 201 individuals with positive TST/QuantiFERON-TB Gold In-Tube test (QFT-GIT) results and 152 pulmonary TB (PTB) patients. A total of 5,591,951 high quality variants were tested with a genomic inflation factor (λ) at 0.997, suggesting that effects from the familial study design were well controlled (S1 Fig). The corresponding Manhattan plot is shown in Fig 1A. We observed a genome-wide significant association on chromosome 10q26.2, corresponding to a cluster of 12 variants and 6 additional variants in high linkage disequilibrium (LD) with P < 5 × 10⁻⁷ in the intronic regions and upstream of C10orf90 (or FATS, HGNC: 26563) (Fig 1B). The top-associated variant was the imputed rs11245088 (odds ratio (OR) = 0.42, 95% confidence interval (CI) 0.39–0.45, P = 1.58 × 10⁻⁸) while the top-associated genotyped variant was rs17155120 (P = 3.71 × 10⁻⁸) (Table 1). Each copy of the minor allele T of rs17155120 conferred protection against M. tuberculosis infection with an OR of being infected for CT vs. CC or TT vs. CT at 0.42 (95%CI 0.35–0.49) (Fig 1C). The intensity cluster plot for rs17155120 showed that the genotype calling was of high quality and separated clearly into 3 genotype groups (S2 Fig). Since all 18 variants in the locus were in high LD (S3 Fig), the imputed variants were likely to have a high imputation quality as suggested by their info score (S1 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Genome-wide association study of resistance to M. tuberculosis infection in Vietnam.

A) Manhattan plot showing results from a genome-wide association study between 185 uninfected subjects (negative for both tuberculin skin test and QuantiFERON-TB Gold In-Tube test) and 353 infected subjects (201 infected individuals positive for both tests and 152 patients with a history of pulmonary tuberculosis) for 5,591,951 variants (minor allele frequency > 5% and info > 0.8) with an unadjusted additive genetic model. The -log₁₀(P value) for each variant (y-axis) is presented according to its chromosomal position (x-axis, build hg19). The dashed line indicates the genome-wide significant threshold at P = 5 × 10⁻⁸. B) Locus zoom plot showing association for the 10q26.2 locus, in a 500 kb window surrounding the top imputed variant rs11245088 (purple diamond). Colors represent pairwise linkage disequilibrium (r²) with rs11245088 as calculated for the Vietnamese Kinh population of 1000 Genomes phase 3. C) Proportion of Vietnamese individuals resistant to M. tuberculosis infection by genotype for the variant rs17155120. Each bar represents the proportion of uninfected subjects among CC individuals (n = 93/344), CT individuals (n = 76/173) and TT individuals (n = 16/21) for the variant rs17155120 in Vietnam.

https://doi.org/10.1371/journal.pgen.1009392.g001

We also performed a GWAS between the 185 uninfected and the 201 infected subjects, excluding the 152 PTB patients (S4 Fig). Despite a smaller sample size, all the 18 variants of the locus were still associated with protection against infection with P < 5.0 × 10⁻⁶_, with similar ORs (for rs17155120, OR = 0.40, 95%CI 0.32–0.49, P = 2.55 × 10⁻⁷) (S2 Table). Similar findings were also observed when considering only the 152 PTB patients as infection reference, with an OR for rs17155120 estimated at 0.50 (95%CI 0.41–0.59, P = 2.10 × 10⁻⁴). These results indicate that PTB patients are an appropriate infection reference group in this analysis.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Association between an additive genetic effect of variants on chromosome region 10q26.2 and resistance to M. tuberculosis infection in Vietnam, France and South Africa.

https://doi.org/10.1371/journal.pgen.1009392.t001

Replication of variants associated with resistance to M. tuberculosis infection in France and South Africa

We tested the effects of the variants of the 10q26.2 locus in a French multi-ethnic HHC cohort (30 uninfected vs. 157 infected subjects) and an admixed familial sample from South Africa (118 uninfected vs. 136 infected subjects). In the French cohort, 17 variants of the cluster could be genotyped or imputed (including 11 genome-wide significant variants), and 3 were replicated at the P < 0.025 level with effect sizes in the same direction as in the Vietnamese cohort. The most significant associated variant was rs56106518 (OR = 0.40, 95%CI 0.30–0.51, P = 2.98 × 10⁻³) (Table 1). In the South African population, 12 variants of the locus were genotyped or successfully imputed (including 10 genome-wide significant variants) and 3 of them showed evidence for replication. The most significant associated variant was rs118037357 (OR = 0.55, 95%CI 0.43–0.67, P = 7.38 × 10⁻³) (Table 1). The top genotyped variant rs17155120 was also replicated in the two cohorts (P_France = 1.51 × 10⁻² and P_SouthAfrica = 1.74 × 10⁻²) (Table 1 and Fig 2A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Genome-wide association study of resistance to M. tuberculosis infection in 3 cohorts from Vietnam, France and South Africa.

A) Forest plot of the association between an additive genetic effect of rs17155120 on chromosome region 10q26.2 and resistance to M. tuberculosis infection. Odds ratios and 95% confidence intervals derived from a linear mixed model, P values, sample sizes and frequency of the effect allele (EAF) are reported by individual cohort and for the random effects meta-analysis. B) Manhattan plot showing results from a genome-wide association study between 333 uninfected subjects and 616 infected subjects for 3,967,482 variants (minor allele frequency > 5% and info > 0.8) with an unadjusted additive genetic model. The -log₁₀(P value) for each variant (y-axis) is presented according to its chromosomal position (x-axis, build hg19). The dashed line indicates the genome-wide significant threshold at P = 5 × 10⁻⁸.

https://doi.org/10.1371/journal.pgen.1009392.g002

Interestingly, the variants across the two replication cohorts were in high LD with each other presenting similar LD patterns as in Vietnam (S5 Fig). The frequencies of the rs17155120 T allele for all 3 cohorts were also similar, ranging from 0.16 to 0.20. The rs17155120 T allele frequency of 0.20 in the Vietnamese cohort was consistent with the Kinh allele frequency of 0.18 in 1000 Genomes (1000G) phase 3 (S6 Fig [18], and S3 Table). The rs17155120 T allele frequency of 0.18 in the multi-ethnic French cohort was also close to the global frequency of 0.16 in all 1000G populations. However, the frequency of 0.16 for rs17155120 T allele in the Souh African cohort was higher than in any 1000G African population that ranged from 0.05 to 0.11. Discrepancy of allele frequencies between Souh African subjects and African populations of 1000G was confirmed by different patterns of LD in a 30 kb region around rs17155120 (S7 Fig), which could be explained by the specific ethnic origin of the Souh African subjects.

Trans-ethnic meta-analysis

Next, we performed a meta-analysis of resistance to M. tuberculosis infection using the GWAS data of the 3 cohorts. The combined analysis was carried out in 333 uninfected and 616 infected subjects. The only genome-wide significant result was observed with the variants of the 10q26.2 region (Figs 2B and S8). The most significant signal was observed at the genotyped variant rs17155120 (OR = 0.50, 95%CI = 0.45–0.55, P = 1.26 × 10⁻⁹) and no heterogeneity was observed across the 3 studies (P_het = 0.31) (Fig 2A and S4 Table).

Functional annotation

The variants of the 10q26.2 region map to the intronic regions and upstream of the tumor suppressor gene C10orf90 (or FATS, for Fragile-site Associated Tumor Suppressor, HGNC: 26563) (Fig 3). According to ENCODE, two associated intronic variants, rs28703703 and rs77513326, are located in a regulatory genomic region characterized by H3K4me1 and H3K27ac histone marks and an active enhancer signature in lymphocytes T helper 17 (Th17). The variant rs77513326 also overlaps ATAC peaks in Th17 cells, memory T cells, natural killer cells and CD8⁺ T cells. We further explored the variants in various expression quantitative trait loci (eQTLs) databases of relevant tissues for the phenotype and found an association between them and expression of the nearby gene ADAM12. In particular, rs28703703 and the genotyped variant rs17155120 displayed decreasing expression of ADAM12 with each minor allele (having a protective effect against M. tuberculosis infection) in monocytes (P = 2.10 × 10⁻³ and P = 4.70 × 10⁻³ respectively) (from https://immunpop.com/kim/eQTL [19]). No association was observed in other immune cell types.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Genomic annotation of the locus on chromosome region 10q26.2.

The upper panel is adapted from Integrative Genomics Viewer (http://igv.org/app/) and the lower panel is adapted from UCSC Genome Browser (http://genome.ucsc.edu/). The 3 vertical grey lines represent the associated variants rs1715120, rs28703703 and rs77513326 (from left to right) that overlap regulatory regions. From top to bottom: H3K4me1 and H3K27ac histone marks from ENCODE, active enhancer in Th17 cells (chromHMM annotation from ROADMAP), chromatin accessibility as represented by ATAC peaks in Th17 cells, memory T cells, NK cells and CD8⁺ T cells [20].

https://doi.org/10.1371/journal.pgen.1009392.g003

Ethics statement

Signed informed consent was obtained from all the participants, and from the parents of enrolled minors. The study was approved by the regulatory authorities in Binh Duong, Vietnam (1366/UBND-VX), the French Consultative Committee for Protecting Persons in Biomedical Research of Henri Mondor Hospital (Créteil, France) and the Stellenbosch University Health Research Ethics Committee (Tygerberg, South Africa).

Study settings and populations

Genotyping, quality control and imputation

A total of 724 individuals from Vietnam (S9 Fig) and 573 from France were genotyped using the Illumina Infinium OmniExpressExome-8-v1 chip (960,212 single nucleotide polymorphisms, SNPs). For the South African cohort (n = 374), the Illumina HumanOmni2.5–8 BeadChip (~2 million of SNPs) was used. All quality control steps were done in each cohort with PLINK v1.9 [50]. Autosomal SNPs with a minor allele frequency (MAF) > 0.01, a genotype call rate > 0.99 and a Hardy-Weinberg (HWE) equilibrium P > 1.00 × 10⁻⁵ were retained. Individuals with a call rate < 95% were excluded (n = 2). Identity-by-descent (IBD) analysis was done to detect duplicated individuals and the members of the pairs with the lower call rate were excluded (n = 1). After the quality control, imputation was performed on 720 individuals and 598,090 variants from Vietnam and 573 individuals and 886,471 variants from France using the Michigan Imputation Server [51] with Eagle2 [52] for the pre-phasing and the 1000G Project as reference panel [53]. For the South African cohort (374 individuals and 1,347,846 variants), the imputation was done on the Sanger Imputation Server with Eagle2 for the pre-phasing and the African Genome Resources as reference panel which includes ~2000 African samples in addition to the individuals from 1000G [54]. Imputed SNPs with an imputation quality info score > 0.8 and MAF > 0.05 were retained for further association analyses (5,591,951 variants in Vietnam, 7,737,070 variants in France and 6,922,541 variants in South Africa).

For each cohort, principal component analysis (PCA) was conducted to evaluate population structure. Genotypes of the individuals from 1000G were used to calculate principal components and data for subjects from the cohorts were projected onto the eigenvectors. The Vietnamese cohort, which was sampled from the Vietnamese Kinh group, was homogenous and clustered with the 1000G East Asian populations (S10 Fig). By contrast, the families in the Val-de-Marne sample showed genetic diversity at the population level with a majority of individuals of European, North African and Sub-Saharan ancestries (S11 Fig). The admixed Souh African subjects, who exhibited genetic diversity at the individual level, were forming a distinct cluster close to the African populations of 1000G (S12 Fig).

Definition of the M. tuberculosis infection phenotype

The definition of the M. tuberculosis infection phenotype relied on both TST and IGRA results. In particular, we used a 5 mm cut-off to determine TST status, based on previous studies in similar settings [7,17] and published guidelines [55]. We explored covariates associated with our infection definition in the entire cohort of enrolled individuals and the subset of those with genotype information were retained for the GWAS (S9 Fig).

For the Vietnamese study, resistance to M. tuberculosis infection was defined by the presence of a negative TST < 5mm and a negative QFT-GIT test result following a protocol suggested by the manufacturer (S1 Text). Infected individuals were defined as subjects presenting both a positive TST ≥ 5mm and a positive QFT-GIT test result. A total of 188 subjects were classified as double negative and 512 as double positive (S13 Fig), among which 185 and 201 subjects were genotyped, respectively. In order to increase the sample size, we also added 152 genotyped PTB patients, consisting of 146 index cases and 6 subjects with a history of PTB, to the infected group (S5 Table). We investigated covariates associated with our infection definition and no significant association was found (S6 Table). Therefore, we conducted an unadjusted genetic association analysis of our binary infection phenotype.

For the French study, which included only HHCs and none of their PTB index cases, contacts could have had one or two screening visits (V1 and V2) with a TST measurement (S14 Fig). A TST was considered negative when the skin induration was (i) < 5 mm at both V1 and V2, (ii) < 5 mm at V1, when only one visit was done. A TST was considered positive when the skin induration was (i) ≥ 5 mm at both V1 and V2, (ii) < 5 mm at V1 and ≥ 10 mm at V2, which reflected true conversions. In-house IGRA was used in this study and provided quantitative levels of IFN-γ production upon early secretory antigenic target 6 (ESAT6) stimulation (S1 Text). A negative IGRA result was defined by a null production of IFN-γ. To determine the optimal positivity cut-off, we built a receiver operating characteristic (ROC) curve with TST status as the observed outcome and the corrected IFN-γ levels (ESAT6 response minus non-stimulated control value) as the predicted outcome among all the contacts enrolled. We selected as positivity threshold the highest sum of sensitivity plus specificity, which was equal to 175 pg/mL (S14 Fig). Then, we defined uninfected subjects as HHCs with a negative TST and a null IFN-γ production (n = 33) and infected subjects as HHCs who presented both a positive TST and a positive IGRA result (IFN-γ production > 175 pg/mL) (n = 147) (S14 Fig). We also looked for covariates associated with our infection definition in this sample of 180 individuals (S7 Table). Age was the only factor significantly associated with our infection definition and was included as covariate in the genetic association analysis that finally included 30 genotyped uninfected and 127 infected subjects (S8 Table).

For the South African study, a 5 mm cut-off was used to identify negative and positive TSTs (S15 Fig). In-house IGRA was also based on the production of IFN-γ upon ESAT6 stimulation (S1 Text). We determined the optimal IGRA positivity cut-off similarly to the French study by building a ROC curve, leading to a threshold of 20.9 pg/mL (S15 Fig). Then, we defined uninfected subjects as HHCs with a negative TST (< 5 mm) and a null IFN-γ production (n = 128), and infected subjects as those with both positive TST and IGRA result (IFN-γ production > 20.9 pg/mL) (n = 152) (S15 Fig). Age was the only factor significantly associated with our infection definition in the whole cohort (S9 Table) and was included as covariate in the genetic association analysis that finally included 118 genotyped uninfected and 136 infected subjects (S10 Table).

Covariates associated with our infection definition were investigated using mixed-effects logistic regression with a random effect per family in each cohort. All the analyses were carried out using R software (version 3.5.2) and related packages “pROC” and “lme4”[56–58].

Genetic association analyses

We conducted genetic association analyses of uninfected vs. infected subjects in the 3 cohorts using a linear mixed-model (LMM) assuming an additive genetic model as implemented in GEMMA v0.98 [59]. To account for the familial relationships, a genetic relatedness matrix (GRM) was used as random effects. In each cohort, the GRM was estimated using centered genotypes after the quality control described above. P values from the likelihood ratio test were reported. For better interpretability, we reported odds ratios (OR) and their 95% confidence intervals (95%CI) after transforming the regression coefficients of the LMM [60]. Manhattan plots of the -log₁₀(P values) and quantile-quantile (QQ) plots were generated using “CMplot” package in R [61]. Regional plots were generated using LocusZoom Standalone v1.4 [62]. Haplotype plots were generated using Haploview [63]. Replication of genome-wide associated variants (P < 5 × 10⁻⁸) in the primary cohort from Vietnam was assessed in the two cohorts from France and South Africa. The observed LD in France and South Africa was slightly weaker as compared with Vietnam, and two LD blocks were inferred by Haploview (S5 Fig). We therefore considered variants at a nominal one sided P value < 0.025 and with a consistent direction of the effect size as replicated.

We also conducted a trans-ethnic meta-analysis by using summary statistics (i.e. beta estimates and their standard errors) from the Vietnamese discovery cohort and the two replication datasets. We used the random-effects model of Han and Eskin implemented in METASOFT [64]. This model assumes effect sizes of exactly zero in all the studies (i.e. no heterogeneity) under the null hypothesis of no associations and allows the effect sizes to vary among studies (i.e. heterogeneity) under the alternative hypothesis. The effect size consistency across studies were determined using the Cochran’s Q statistic. Allelic effect estimates were also derived on the log-odds scale.

Functional annotation

We used the UCSC Genome Browser [65] to identify associated variants which may overlap with known regulatory regions: 1) histone marks from the ENCODE project [66], 2) chromatin state annotated by ChromHMM on the basis of ROADMAP [67,68] and 3) chromatin accessibility determined by assay for transposase-accessible chromatin using sequencing (ATAC-seq) from immune cell-types [20]. We also looked at the associated variants in eQTL databases which focused on gene expression in monocytes [19,69], T cells [70], macrophages [71], and various types of immune cells [72].