Two loci control tuberculin skin test reactivity in an area hyperendemic for tuberculosis

We studied 128 informative families (including 186 parents and 350 children) comprising two to six children with available TST induration size readings (Table I). TST reactivity of the 350 children is shown in Fig. 1 A. Reactivity extended from 0 to 43 mm with a median size of 11.15 mm. Mean age at the time of TST was 11.5 yr (SD, 4.9 yr) and the sex-ratio was 1. As expected as a result of the cumulative exposure to M. tuberculosis in this high incidence area, age had an important impact on TST positivity (Fig. 1 B). Closer inspection of the extent of TST reactivity showed a clear bimodal distribution with 140 individuals (40%) having a value of zero and 210 (60%) having values that approximately followed a normal distribution centered around 16 mm. Only two children had a TST between 1 and 5 mm. Given the age distribution of study subjects (i.e., mean age of 11.5 yr), such a high proportion of TST⁺ individuals is consistent with an annual risk of TB infection that is >4% in the study area (Kritzinger et al., 2009). This TST distribution is strongly suggestive of a gene or group of genes having an impact on TST positivity by itself and a different gene or set of genes controlling the extent of TST reactivity. To test this hypothesis, we genotyped ∼6,000 single nucleotide polymorphism (SNPs) in all 536 individuals of the 128 nuclear families and we performed two complementary linkage analyses. The first analysis focused on the phenotype of positivity by itself, i.e., TST = 0 versus TST > 0 (TST-BINa). The second analysis focused on the quantitative phenotype of extent of TST reactivity (TST quantitative trait locus [QTL]), with a particular emphasis on individuals with a TST > 0. Before linkage analysis, the two phenotypes were adjusted on age, gender, and previous clinical TB to remove as many environmental or nonspecific genetic sources of TST variation as possible. Subsequent linkage analyses were therefore performed on the resulting adjusted residuals shown in Fig. 1 C (TST-BINa) and Fig. 1 D (TST-QTL) that are, by definition, quantitative traits.

We first searched for genetic factors controlling TST-positive response by itself. Results of the TST-BINa linkage analysis are shown in Fig. 2. Information content (IC) was very high across all autosomes with mean genome-wide information of 94.4% (from 80.6 to 98.7%). This IC is excellent and corresponds to the level only accomplished at the fine-mapping stage in microsatellite-based genome scans. Linkage analysis of TST-BINa identified a significant linkage signal on chromosomal region 11p14 (LOD score = 3.81) at chromosomal position 26.37 Mb (IC = 91.15%, P = 1.4 × 10⁻⁵; Fig. 2, bottom). This level of statistical support exceeds the stringent threshold of significance (LOD score = 3.6) for GWL scans (Lander and Kruglyak, 1995). The 1-LOD confidence interval for the location of the major locus (corresponding to the 90% confidence interval for the location of the QTL underlying the linkage peak) spanned from 22.35 to 28.82 Mb (Fig. 2, bottom). A suggestive linkage signal was also observed on chromosomal region 5p15 (LOD score = 2.39, P = 0.0005, IC = 90.5%) at the same position as the TST-QTL locus (see following paragraph). In addition, five weaker linkage peaks with P < 0.01 (i.e., LOD score > 1.17) were observed on chromosomal regions 3p24 (LOD score = 1.40, P = 0.0056, IC = 96.8%), 4q28 (LOD score = 1.57, P = 0.0036, IC = 94.7%), 15q26 (LOD score = 1.20, P = 0.0094, IC = 93.4%), 19q13 (LOD score = 1.83, P = 0.0018, IC = 92.4%), and 20p13 (LOD score = 1.23, P = 0.0087, IC = 93.9%; Fig. 2, top). As a substantial proportion of TST = 0 persons are most likely resistant to infection with M. tuberculosis (see Concluding remarks), these data therefore point toward the identification of one major locus (TST1) controlling human resistance to M. tuberculosis. A list of the known genes located in the 1-LOD confidence interval is given in Table S1.

Results of the TST-QTL linkage analysis are shown in Fig. 3. A highly significant linkage signal was observed on chromosomal region 5p15 with a multipoint LOD score of 4.00 (P = 9 × 10⁻⁶) at position 2.70 Mb (TST2). Again, this is substantially above the commonly accepted threshold for significance in GWL studies (Lander and Kruglyak, 1995). The 1-LOD confidence interval for location of the QTL was small, spanning ∼2 Mb from 1.07 Mb to 3.23 Mb (Fig. 3, bottom). In addition to this major locus, four weaker linkage peaks with P < 0.01 were observed in chromosomal regions 11p14 (LOD score = 1.47, P = 0.0046, IC = 93.1%), at the same position as the TST-BINa locus 10p15 (LOD score = 1.76, P = 0.002, IC = 93.3%), 13q21 (LOD score = 1.51, P = 0.0042, IC = 96%), and 22q11 (LOD score = 1.34, P = 0.0065, IC = 97%). Therefore, these results support the hypothesis of a second major locus (TST2) on chromosomal region 5p15 controlling the intensity of TST reactivity as a quantitative trait, i.e., the intensity of T cell–mediated DTH to tuberculin. An interesting observation is the mirror effect between the two previous analyses with a suggestive linkage signal on chromosomal region 5p15 in the analysis of TST-BINa and on chromosomal region 11p14 in the analysis of TST-QTL. This is suggestive of some level of redundancy between the two phenotypes with the most straightforward explanation being individuals coded as 0 in both analytical approaches. To further explore this observation, we performed a linkage analysis of TST-QTL among individuals with TST > 0 mm only. Despite a dramatic reduction in family numbers (68 families comprising 164 children), we still found significant evidence of linkage at chromosomal region 5p15 (LOD score = 2.17, P < 8 × 10⁻⁴). In contrast, the linkage signal on 11p14 totally disappeared (LOD score = 0.14, P = 0.21). These results may reflect a sequential genetic control of TST reactivity, with TST1 controlling TST positivity by itself and, presumably, T cell–independent innate resistance to M. tuberculosis, and TST2 controlling the intensity of TST reactivity as a quantitative trait and, presumably, the intensity of T cell–mediated DTH to tuberculin.

Encouraged by the narrow linkage peak, we decided to fine map the TST2 locus. Based on the 1-LOD confidence interval for the location of the underlying QTL for TST2, we selected a chromosomal region extending from 1.1 to 3.2 Mb on the physical map of chromosome 5 (Fig. 4) for association studies of TST-QTL in a familial sample including the 128 nuclear families used in the linkage study and two additional trios (one child and her/his two parents). This interval contains 13 genes (NKD2→C5orf38; National Center for Biotechnology Information build 36 [http://www.ncbi.nlm.nih.gov]). We genotyped 133 SNPs in the targeted interval, 113 of which were considered suitable for association analysis (Fig. 4 and Table S2). To optimize gene coverage, we used a genocentric strategy for SNPs selection. Tag SNPs were identified for each known gene including 2 kb of their 5′ and 3′ regulatory regions. To capture most genetic variation, a very stringent tag-SNP selection scheme was used (r² = 0.8 with all tagged SNPs and minor allele frequency > 5%). Because no genes have been annotated on chromosomal intervals 2–2.7 Mb and 2.9–3.2 Mb, few SNPs were genotyped in these chromosomal segments and none were significant. In univariate analysis, nine SNPs (rs4975579, rs6554677, rs1801075, rs250682, rs10475030, rs11747565, rs2922061, rs1018120, and rs2232376) were significantly associated with TST-QTL at the 0.05 level (Fig. 4). Analysis of the linkage disequilibrium (LD) patterns between these nine SNPs showed no LD, supporting the hypothesis of independent effects even though it is likely that some of these signals are false positives (unpublished data). Out of these nine SNPs, rs250682 located in SLC6A3 (solute carrier family 6 member 3) displayed the strongest association with TST-QTL (P = 0.001). Allele G of rs250682 (q = 0.36) has a dominant effect and is associated with lower values of the TST-QTL. To ensure that SNP rs250682 was responsible for the observed association, we screened the Yoruba sample of the HapMap database (www.hapmap.org) in a region spanning 5 Mb around rs250682. Only four SNPs (rs250681, rs40358, rs403636, and rs464049) were correlated with an r² > 0.5 with rs250682. All of these four SNPs are located in the SLC6A3 gene. These results identified SLCA63 as a promising target for further study on the intensity of T cell–mediated DTH to tuberculin.

We report in this paper that absence of TST reactivity (zero induration) has a major human genetic component and that the corresponding locus, TST1, maps to chromosome region 11p14. As the risk of developing TB for persons with M. tuberculosis exposure but TST = 0 was previously shown to be extremely small, such TST = 0 persons are most likely not infected with M. tuberculosis (Rose et al., 1995; ATS-CTC-IDSA, 2000). Persons who experience continued high exposure to M. tuberculosis and display persistent lack of T cell sensitization are much more likely to be naturally resistant to infection with M. tuberculosis than intrinsically deficient in mounting a DTH response. An estimate of infection-resistant persons can be obtained from countries where TB is highly endemic. In such conditions, where exposure to the tubercle bacillus is intense and sustained, ∼20% of individuals remain TST negative (Rieder, 1999). In our families, it is not possible to distinguish between individuals who have not been exposed to M. tuberculosis and those who have been exposed but are naturally resistant to infection. However, because the lack of exposure is difficult to reconcile with a genetic component, the most parsimonious explanation of our data is that the TST1 locus reflects T cell–independent resistance to M. tuberculosis infection. Consequently, a substantial fraction of TST = 0 persons in the study families must have been exposed to M. tuberculosis but, as a result of very effective T cell–independent responses, those persons were resistant to M. tuberculosis before T cell sensitization could occur. Our findings make it likely that being spared from infection with M. tuberculosis is not simply a matter of lack of exposure but reflects genetically controlled profound differences of T cell–independent resistance among exposed persons. Although it is possible that this innate resistance can be overcome by extreme exposure pressure, the targeted strengthening of T cell–independent resistance loci in infection susceptible persons offers an attractive avenue of protection from TB disease. The identification of the TST1 locus that belongs to a much larger region than TST2 will follow the strategy successfully applied in leprosy (Mira et al., 2004; Alcaïs et al., 2007).

Our study also provides strong evidence that a major QTL (TST2) mapping to chromosome region 5p15 is involved in the control of TST extent in families living in a hyperendemic region for TB. Consistent with our results, the same region showed some evidence of linkage (P < 0.05) with persistently low TST reactivity in a familial sample from Uganda (Stein et al., 2008). TST2 is the first non-MHC locus that has an impact on intensity of T cell–mediated DTH to tuberculin. This identifies TST2 as important immune modulator and possible confounder in the numerous studies that have aimed at correlating extent of TST with risk for future clinical TB disease. A more direct link of TST2 with risk of disease is perhaps provided by the repeated observation of the 5p15 region as location of a sarcoidosis susceptibility locus (Iannuzzi et al., 2007). Sarcoidosis is characterized by the immune paradox of extensive local inflammation (granuloma and cytokine secretion) associated with suppression of the immune response to tuberculin (Miyara et al., 2006). Multiple studies have investigated the association between mycobacterial antigens and sarcoidosis but no clear consensus has emerged (Iannuzzi et al., 2007). Our results suggest that human genetics could serendipitously connect the mechanisms governing sarcoidosis and T cell–dependent DTH to mycobacteria. The only other example of DTH for which response regulators have been reported in humans is the Montenegro skin test to leishmania antigen (Jeronimo et al., 2007a,b). However, none of these regulators mapped to the 5p15 (or the 11p14) chromosomal region, suggesting that TST2 immune regulation involves an element of specific antigen recognition. Our first genocentric association scan of the TST2 interval identified SLC6A3 as a promising candidate for TST2. Experimental deletion of this gene in mice resulted in their significantly reduced ability to mount a DTH reaction against ovalbumin (Kavelaars et al., 2005). SLC6A3 is related to NRAMP1 (alias SLC11A1), which has previously been shown to have a major impact on the extent of the granulomatous response to Mycobacterium leprae antigen (Alcaïs et al., 2000). Nevertheless, additional studies will be required to fully understand the molecular basis of the TST2 locus.

The gene regions previously reported to play a role in clinical TB (Schurr and Kramnik, 2008), including the only major TB susceptibility locus on chromosome 8p13 (El Baghdadi et al., 2006), do not overlap with the two main chromosomal regions identified in this study. For the TST1 locus on chromosome 11p14, this is not surprising because genetic studies are generally done on TST-positive (i.e., infected) subjects, which precludes identification of infection resistance genes. For the TST2 region on chromosome 5p15, the lack of overlap with TB susceptibility genes is more surprising. However, because TST2 has a strong impact on TST induration size, which is correlated with risk of subsequent TB disease, a further study of TST2 on susceptibility to clinical TB may provide new insights into the mechanism of advancement from infection to TB disease. HLA class II genes are validated TB susceptibility loci, and MHC alleles are known from the mouse model to have a strong impact on DTH (Schurr and Kramnik, 2008). These observations made the chromosome 6 HLA region a prime candidate for loci impacting on TST induration size. Yet we failed to observe any evidence for a role of HLA genes on TST. That the strongest genetic effect on tuberculin reactivity is caused by a non-HLA locus directly leads to the question of whether this or other non-HLA loci will be modulators of vaccine-induced anti-mycobacterial protective immunity. Although the impact of the TST1 and TST2 loci on progression from infection to TB disease is presently unknown, the two loci which control T cell–independent resistance to infection with M. tuberculosis (TST1) and T cell–dependent intensity of tuberculin reactivity (TST2) represent an important expansion in our understanding of TB immunity.

Nuclear families (i.e., parents and offspring) with at least two children were enrolled from Ravensmead and Uitsig. Ravensmead/Uitsig (R/U) is a suburban area of Cape Town, South Africa, with a population of 38,656 living in a 3.5-km² area. HIV prevalence at the time of family enrolment was 5.2% in the overall population and <2% in the pediatric population. The population is stable and there are few homeless people or migratory workers. All individuals of the sample belong to the Cape Colored ethnic group. There was no requirement for subjects to be household contacts of TB cases. However, this area has a notification rate of TB (all cases) of 761/100,000 and of new bacteriologically confirmed cases of 313/100,000 (Verver et al., 2004), which is suggestive of a high level of exposure to M. tuberculosis. The enrolment strategy was to target large households to allow later reconstruction of nuclear families. In addition, if TST reactivity was known at time of enrollment, households that contained both strongly TST-positive and TST-negative subjects were prioritized for enrollment, as this sampling strategy has been shown to be the most powerful for linkage analysis.

Subjects who had had clinical TB disease in the 2 yr preceding the study were excluded. It is of note that the distribution of TST values among siblings of individuals who were excluded because they developed TB in the last two years did not differ from the overall distribution (i.e., 33 vs. 40% of null values and mean TST of 15.8 vs. 16 mm). Individuals who were HIV positive, pregnant, or using immunomodulatory chemotherapy were also excluded at the time of enrolment. BCG vaccination at birth is routine in the study area and was therefore not a confounding factor in our study because several studies have shown that the impact of BCG at birth on the TST vanishes rapidly (Menzies, 2000). Similarly, all individuals belong to the Cape Colored ethnic group, therefore limiting the risk of genetic heterogeneity. Informed consent was obtained from all study participants. Protocols involving human subjects were approved by the Stellenbosch University Health Research Ethics Committee (Tygerberg, South Africa), the University of Cape Town (Cape Town, South Africa), and the Research Ethics Board at the Research Institute of the McGill University Health Centre (Montreal, Canada).

TSTs were performed by specially trained health care providers using the Mantoux method with M. tuberculosis PPD RT23 (2 tuberculin units; Statens Serum Institut, Copenhagen, Denmark). TST reactivity was read between 48 and 72 h after the skin test was performed and the diameter of induration was measured in millimeters using a set of calipers calibrated to the nearest 0.5 mm. Two phenotypic definitions were used (Fig. 1). First, we dichotomized the TST distribution using 0 mm as the threshold to study the TST positivity (TST-BINa). Second, we analyzed the extent of the TST reactivity as a quantitative trait (TST-QTL) by Tobit regression.

Before linkage analysis, TST-BINa and TST-QTL phenotypes were all adjusted on previous clinical TB (at least 2 yr preceding the study), sex (male or female), and age (in years). Pearson and Tobit residuals were used for linkage analysis of TST-BINa and TST-QTL, respectively. TST-BINa was adjusted by means of logistic regression (Fig. 1 C), as implemented in the PROC LOGISTIC of the SAS software v9.1 (SAS institute, Cary, NC). TST-QTL was adjusted by means of the Tobit censored regression as implemented in the PROC QLIM of the SAS software, with the censoring threshold fixed at 0 (Fig. 1 D). Note that because the linkage analysis was performed on the residuals generated by the adjustment procedure, both phenotypes under study (TST-BINa and TST-QTL) were analyzed as quantitative traits. For each phenotype, the best fitting model for the age effect was determined among a set of multivariate fractional polynomials (FP) models as proposed in (Royston et al., 1999). First (FP1) and second (FP2) degree FP models were fitted with power p (i.e., age to the power p) for FP1 and p and q for FP2 chosen from −2, −1, −0.5, 0, 0.5, 1, 2, and 3, with 0 denoting log transformation. Among the 44 possible combinations, the model providing the smaller Akaïke Information Criteria was selected. A simple FP1 model was retained with power p = 0.5 for TST-BINa, whereas a FP2 model was retained for TST-QTL (power p = 1 and q = 2).

High-density genotyping for the linkage study was performed at the Centre National de Génotypage (Paris, France) with the linkage IVb panel (Illumina), containing >6,000 SNPs. 11 nonpolymorphic SNPs and 79 SNPs with a call rate <80% were excluded. None of the remaining SNPs showed departure from Hardy-Weinberg equilibrium among the founders at the 0.001 level. Pairwise LD analysis between adjacent SNPs was performed using Haploview software (Barrett et al., 2005) in the 186 parents of our sample. In our sample, pairwise LD between adjacent SNPs was very weak, with most of SNPs-pairs having an r² < 0.1 (mean pairwise r² = 0.07).

Genotyping for fine mapping of TST2 on chromosome region 5p15 was performed at the McGill University and Genome Quebec Innovation Centre. 133 SNPs spanning the 2.2-Mb targeted interval on chromosome region 5p15 (i.e., the 1-LOD confidence interval for the location of the QTLs underlying the linkage peak) were selected on the basis of their location within known genes in the interval and on information regarding bin structure and allelic frequencies publicly available from the International HapMap project (http://www.hapmap.org/). These 133 SNPs were genotyped on the high-throughput MassARRAY platform (SEQUENOM, Inc.), which uses the iPLEX assay to incorporate mass-modified terminal nucleotides in the SBE step, which are then detected by MALDI-TOF MS (Griffin and Smith, 2000). 10 SNPs were excluded because they were nonpolymorphic and 1 SNP was excluded because of a call rate <80%. SNPs with a minor allele frequency <5% were also excluded. Finally, two SNPs that showed departure from Hardy-Weinberg equilibrium among the founders at the 0.01 level were excluded (Table S1).

As suggested by Thompson et al. (2006), we checked for population substructure to minimize genetic heterogeneity of the sample before linkage analysis. We performed a principal component analysis of the 5,567 autosomal SNPs in the 186 genotyped founders of our sample as implemented in the SMARTPCA software (Patterson et al., 2006). The principle of this method is to determine the major axes of genetic variation in the sample and to output each individual's coordinates along axes of variations, without formally clustering individuals into discrete population. No population substructure was found in our data but we identified five outliers distributed in four families that were excluded from the analysis.

We performed quantitative model-free multipoint linkage analysis on the Pearson residuals of the TST-BINa logistic regression and on the residuals of the TST-QTL censored Tobit regression. As residuals were not normally distributed and pedigrees were extended, we used the maximum-likelihood binomial (MLB) method extended to quantitative trait linkage analysis (MLB-QTL v.3.0 available for download at www.hgid.net; Alcaïs and Abel, 1999) that has been implemented in an extension of the GENEHUNTER program.

The MLB approach considers the sibship as a whole and does not make any assumption about the distribution of the phenotype. The idea of the MLB method is to introduce an individual latent binary variable, which captures the linkage information between the observed quantitative trait and the marker. The probability that the latent variable for an individual takes the value 0 or 1 depends on the quantitative trait value of the individual and on a link function, which can be parametric (like the standard cumulative normal distribution) or empirical (derived from the observed phenotypic distribution).

The test of linkage is a maximum-likelihood ratio test that compares the likelihoods under the null hypothesis of no linkage (H₀) and the hypothesis of linkage. The test statistic is asymptotically distributed as a 50:50 mixture of χ² distribution with 0 and 1 degree of freedom and can be expressed as a classical LOD score. In this study, we used an empirical link function based on the deciles of the residuals distribution as suggested in Alcaïs and Abel (1999). Confidence intervals for the location of QTL underlying the linkage peaks were calculated by use of the support interval method in which a 1-LOD interval corresponds to a 90% confidence interval (Dupuis and Siegmund, 1999). LOD scores 3.6 (P ≤ 2 × 10⁻⁵) and 2.2 (P ≤ 7 × 10⁻⁴) were considered as significant and suggestive of linkage, respectively, as proposed in Lander and Kruglyak (1995).

It has been shown that linkage analysis using a dense SNP panel in familial samples with missing parents could inflate the type I error rate in regions where SNPs display high levels of LD (Huang et al., 2004). Therefore, although the LD was very weak at the genome-wide level and between markers underlying the two different linkage peaks (mean pairwise r² = 0.07), we removed the SNP with the lowest minor allele frequency for each pair of SNPs with an r² > 0.1. This had no impact on the LOD scores or the IC, ruling out the possibility of TST1 and TST2 being false positive linkage signals caused by LD patterns.

We performed family-based association testing on the 113 SNPs located in the 2.2-Mb targeted interval on chromosome region 5p15. The analysis was performed using the extension of the transmission disequilibrium test to quantitative phenotypes as implemented in FBAT software v2.0.3 (Laird et al., 2000).

Table S1 presents a list of the known genes that are located in the 1-LOD interval of TST1. Table S2 shows the SNPs that were genotyped for the fine mapping of TST2.

We sincerely thank all members of the community who participated in this study. We are grateful to Christine Johannisen, Magdalene Kennedy, Ronell Taylor, and Maria Gallant for their help with field work and blood collection. Gian van der Spuy, Neal Rogers, and Mari de Klerk assisted with data capture and database management. We are grateful to the Centre National de Génotypage for providing the genotypes.

This work was supported by grants from the Sequella/Aeras Global Tuberculosis Foundation and the Canadian Institutes of Health Research (CIHR) to E. Schurr and, in part, by a grant from the Grand Challenge exploration of the Bill and Melinda Gates Foundation to A. Alcaïs. C.J. Gallant is the recipient of a fellowship from the CIHR Strategic Training Centre in the Integrative Biology of Infectious Diseases and Autoimmunity. A. Cobat was supported in part by a grant from Fondation de la Recherche Médicale, J.L. Casanova, L. Abel, and A. Alcaïs are supported by the Assistance Publique-Hopitaux de Paris, Program de Recherche Fondamentale en Microbiologie, Maladies Infectieuses et Parasitaires, and the Agence Nationale de la Recherche of the Ministère Francais de l'Education Nationale de la Recherche et de la Technologie, and J.L. Casanova was an International Research Scholar of the Howard Hughes Medical Institute (HHMI; 2005–2008). E. Schurr is a Chercheur National of the Fonds de Recherche en Santé du Québec and an International Research Scholar of the HHMI.

The authors have no conflicting financial interests.