Introduction
Tuberculosis (TB) caused by Mycobacterium tuberculosis (M. tuberculosis) is a global health epidemic and the leading cause of death due to a single infectious agent (World Health Organization [WHO], 2017). In 2016 1.3 million TB deaths were reported in HIV negative individuals and an additional 374000 deaths related to TB/HIV co-infection were recorded. The majority of these deaths occurred in southeast Asian and African countries (World Health Organization [WHO], 2017). TB is a complex disease, influenced by environmental and behavioral factors such as socio-economic status and smoking, as well as definite human genetic components. The contribution of the host genes to disease has been highlighted by numerous investigations, including animal (Pan et al., 2005), twin (Comstock, 1978; Sorensen et al., 1988; Flynn, 2006), linkage (Bellamy et al., 2000; Greenwood et al., 2000) and candidate gene association studies (Schurz et al., 2015). More recently genome-wide association studies (GWAS) in diverse populations have been done (Thye et al., 2010, 2012; Oki et al., 2011; Mahasirimongkol et al., 2012; Png et al., 2012; Chimusa et al., 2014; Curtis et al., 2015; Grant et al., 2016; Sobota et al., 2016; Qi et al., 2017).
Interestingly another influential factor in TB disease development is an individual’s biological sex, which has been largely ignored in past TB studies and was usually only used as a covariate for adjusting association testing statistics. In 2016, males comprised 65% of the 10.4 million recorded TB cases, indicating that the TB notification rate is nearly twice as high in males as in females (World Health Organization [WHO], 2017). While socio-economic and behavioral factors do influence this ratio, it does not fully explain the observed sex-bias (Jaillon et al., 2017). Another factor that influences sex-bias is the effect that sex hormones (estrogen and testosterone) have on the immune system. Estrogen is an immune activator, upregulating pro-inflammatory cytokines (TNFα), while testosterone is an immune suppressor, upregulating anti-inflammatory cytokines (IL-10) (Cutolo et al., 2006). This could explain why men are more susceptible to infectious diseases compared to females (Jaillon et al., 2017). However, as sex-based differences in immune responses differ even between pre-pubertal boys and girls, as well as between post-menopausal women and elderly men, it shows that sex hormones do not fully explain the sex-bias (Klein et al., 2015). Thus, it has been proposed that the X chromosome and X-linked genes directly contribute to the observed sex-bias.
There are approximately 1,500 genes on the X chromosome, many of which are involved in the adaptive or innate immune system (Brooks, 2010). Since females have two X chromosomes, one requires silencing in order to equalize dosage of gene expression to that of men who only have one X chromosome. This silencing occurs randomly in each cell, making females functional mosaics for X linked genes and giving them a major immunological advantage over males (Jaillon et al., 2017). As males are haploid for X-linked genes any damaging polymorphisms or mutations on the X chromosome will have a more pronounced immunological effect in males than in mosaic females, thereby influencing the sex-bias (Abramowitz et al., 2014).
To date, eleven GWAS investigating susceptibility to clinical TB have been published (Thye et al., 2010, 2012; Oki et al., 2011; Mahasirimongkol et al., 2012; Png et al., 2012; Chimusa et al., 2014; Curtis et al., 2015; Grant et al., 2016; Sobota et al., 2016; Omae et al., 2017; Qi et al., 2017). There has not been significant overlap between the 11 published TB GWAS, but it seems that replication is more likely when populations with similar genetic backgrounds are compared: the WT1 locus was associated with disease in populations from West and South Africa (Thye et al., 2012; Chimusa et al., 2014). Critically, genotyping microarrays that did not fully accommodate African genetic diversity were used in these studies (Thye et al., 2010, 2012; Chimusa et al., 2014; Curtis et al., 2015; Grant et al., 2016). It is therefore possible that unique African-specific susceptibility variants were not tagged by these initial arrays, since LD blocks are shorter in African populations (Campbell and Tishkoff, 2008). Moreover, none of the GWAS included or examined the X chromosome or sex-stratified analysis of the autosomes as was done in an asthma cohort (Mersha et al., 2015). Genetic differences between asthmatic males and females were identified on the autosome, with certain alleles having opposite effects between the sexes. Candidate gene association studies provide independent confirmation of the involvement of the X chromosome in TB susceptibility, through the association of X-linked TLR8 susceptibility variants with active TB. Davila et al. (2008) investigated 4 TLR8 variants (rs3761624, rs3788935, rs3674879, and rs3764880) in an Indonesian cohort and showed that all variants conferred susceptibility to TB in males but not females. The results for males were validated in male Russian individuals (Davila et al., 2008). These results were validated for rs3764880 in Turkish children, but no significant association was found for rs3764879 (Dalgic et al., 2011). Hashemi-Shahri et al. (2014) found no significant TLR8 associations in an Iranian population, while rs3764880 was significantly associated with TB susceptibility in both males and females in a Pakistani cohort (Bukhari et al., 2015). In admixed South African individuals rs3764879 and rs3764880 were significantly associated in both males and females, while rs3761624 was only significantly associated in females (Salie et al., 2015). Interestingly, in this cohort opposite effects were consistently found between the sexes for the same allele in all investigated TLR8 variants (Salie et al., 2015), echoing the asthma findings of Mersha et al. (2015). Finally, in a Chinese cohort rs3764879 was significantly associated with TB disease in males but not females. While many of these variants did not reach genome wide significance they still provide evidence of the involvement of X-linked genes in TB susceptibility.
We report the first TB susceptibility GWAS with a specific focus on sex-stratified autosomal analysis and the X chromosome to elucidate the male sex-bias. Individuals from the unique five-way admixed SAC population, with ancestral contributions from Bantu-speaking African, KhoeSan, European, South and East Asian groups were genotyped in this study (Chimusa et al., 2013; Daya et al., 2013). These genetic contributions are due to both the complex colonization history of South Africa and the country’s importance as a refreshment station on major trade routes during the fifteenth to nineteenth century (de Wit et al., 2010; Uren et al., 2016). This is therefore the first GWAS in the SAC that uses an array (Illumina Multi Ethnic Genotyping Array, see section “Genotyping”) specifically designed to detect variants in the 4 most commonly studied populations, making it the most suitable platform for diverse and admixed populations at the time of genotyping.
Materials and Methods
Study Population
Study participants were recruited from two suburbs in the Cape Town metropole of the Western Cape. These suburbs were chosen for its high TB incidence and low HIV prevalence (2%) at the time of sampling (1995–2005) (Kritzinger et al., 2009). Approximately 98% of the residents in these suburbs self-identify as SAC and have similar socio-economic status, which reduces confounding bias in the association testing (Chimusa et al., 2014). The cohort consists of 420 pulmonary TB (pTB) cases, bacteriologically confirmed to be culture and/or smear positive and 419 healthy controls from the same suburbs. Approximately 80% of individuals over the age of 15 years from these suburbs have a positive tuberculin skin test (TST), indicating exposure to M. tuberculosis (Gallant et al., 2010). All study participants were over 18 years of age and HIV negative.
Approval was obtained from the Health Research Ethics Committee of Stellenbosch University (project registration number S17/01/013 and 95/072) before participant recruitment. Written informed consent was obtained from all study participants prior to blood collection. DNA was extracted from the blood samples using the Nucleon BACC Genomic DNA extraction kit (Illumina, Buckinghamshire, United Kingdom). DNA concentration and purity was checked using the NanoDrop® ND-1000 Spectrophotometer and NanoDrop® v3.0.1 software (Inqaba Biotechnology, Pretoria, South Africa). The study adhered to the ethical guidelines as set out in the “Declaration of Helsinki, 2013 (World Medical Association [WMA], 2018).
Genotyping
Genotyping was done using the Illumina MEGA (Illumina, Miami, United States) which contains 1.7 million markers from various ethnicities making it highly suitable for diverse and admixed populations. The array is based on novel variants identified by the Consortium on Asthma among African ancestry populations in the Americas (CAAPA), the Illumina human core content for European and Asian populations as well as multi-ethnic exome content from African, Asian and European populations. The array also contains ancestry informative markers specific to the SAC population. While the KhoeSan population is not highly represented on the array, which could lead to a certain level of ascertainment bias, at the time of genotyping it was the most suitable platform for this diverse and admixed populations. Genome studio v2.04 (Illumina, Miami, United States) was used for SNP calling to calculate intensity scores and to call common variants (MAF≥5%), followed by analysis with zCall to recall rare genotypes (MAF < 5%) (Goldstein et al., 2012).
Genotyping Quality Control
Quality control (QC) of the genotyping data was done using the XWAS version 2.0 software and QC pipeline to filter out low quality samples and SNPs (Chang et al., 2014; Gao et al., 2015). Data were screened for sex concordance, relatedness (up to third degree of relatedness) and population stratification (as determined by principal component analysis). Genotypes for males and females were filtered separately in order to maintain inherent differences between the sexes. SNPs were removed from the analysis if missingness correlated with phenotype (threshold = 0.01) as well as individual and SNP missingness (greater than 10%), MAF (less than 1%) and Hardy–Weinberg equilibrium (HWE) in controls (threshold = 0.01). Filtering continued iteratively until no additional variants or individuals were removed. Overlapping markers between the sexes were merged into a single dataset. X chromosome genotypes were extracted and variants were removed if the MAF or missingness was significantly different between the sexes (threshold = 0.01). A flow diagram explaining quality control steps and association testing of the data is shown in Supplementary Figure S1.
Admixture
The SAC population is a 5-way admixed population with ancestral contributions from Bantu-speaking African populations, KhoeSan, Europeans, and South and East Asians (Chimusa et al., 2013; Salie et al., 2015). To avoid confounding during association testing the ancestral components are included as covariates (Daya et al., 2014a). Admixture was estimated for the autosome (chromosomes 1–22) and the X chromosome separately using the software ADMIXTURE (v1.3) (Alexander et al., 2009) and reference genotyping data for 5 ancestral populations. The reference populations used to infer ancestry were European (CEU) and South Asian (Gujarati Indians in Houston, Texas and Pathan of Punjab) extracted from the 1000 Genomes Phase 3 data (Sudmant et al., 2015), East Asian (Han Chinese in Beijing, China), African (Luhya in Webuye, Kenya, Bantu-speaking African, Yoruba from Nigeria) and San (Nama/Khomani) (Uren et al., 2016; Martin et al., 2017). Due to the limited number of individuals available for each reference population the SAC data had to be divided into 21 groups to equal the number of individuals per reference population. The number of individuals per reference population and admixed population has to be kept consistent in order to maximize the accuracy of the admixture results by not over-representing one particular population in the analysis. Therefore admixture inference was done separately for each of the 21 SAC groups, referred to as running groups. Each running group was analyzed five times at different random seed values. The results for each individual were averaged across the five runs in order to obtain the most accurate ancestry estimations [Shringarpure et al., 2016). Four ancestral components (African, San, European, and South Asian (Salie et al., 2015)] were included as covariates in the logistic regression association testing with the smallest component (East Asian) excluded in order to avoid complete separation of the data.
Association Analysis
SNP Based Association Analysis
Autosomal TB association testing was done with sex-stratified and combined datasets using the additive model in PLINK (version 1.71) (Purcell et al., 2007) in order to detect sex-based differences. TB association testing for the X chromosome were done separately in males and females using XWAS (version 2) and the results were combined using Stouffers method in order to obtain a combined association statistic (Chang et al., 2014; Gao et al., 2015). A sex-differentiated test was conducted for the X chromosome using the XWAS software to test for significant differences in genetic effects between males and females. SNP based association testing (sex-stratified or not) compares the frequency of alleles between cases and controls to determine if a specific allele co-occurs with a phenotype (TB) more often than would be expected by chance. The sex-differentiation test on the other hand compares the effect size (OR) of a variant between the sexes to determine if a variant has a different effect on risk between the sexes. The sex-differentiation test is explained in more detail by Chang et al. (2014). X chromosome inactivation states were also included in the association testing as covariates using a method developed by Wang et al. (2014). To include inactivation states in the association analysis the most likely state was determined for each SNP. A variant can either be inactivated, or it can be skewed toward the deleterious or normal allele or the variant can escape inactivation. To determine which of the four states is most probable the likelihood ratio for each one was calculated and the inactivation state that maximized the likelihood ratio was applied to the SNP in question. This was done for each variant as inactivation states vary along the X chromosome [for a detailed description see Wang et al. (2014)]. Ancestry, sex and age were included in the analyses as covariates where applicable. Information on other risk factors known to influence TB susceptibility such as smoking and alcohol consumption was not available for this study cohort and could not be included as covariates. Multiple testing correction was done using the SimpleM method (Gao et al., 2010), which adjusts the significance threshold based on the number of SNPs that explains 95% of the variance in the study cohort. This method is less conservative than Bonferroni correction and is a close approximation of permutation results in a fraction of the time. For the autosome the genome-wide significance threshold was set to 5.0e-8 (Panagiotou and Ioannidis, 2012).
Gene Based Association Analysis
Gene-based association testing groups SNPs together and thus decreases the multiple testing burden and increase power to detect an association. Gene-based association testing was done using the XWAS v2 scripts, which were implemented using the Python2 (version 2.7.10) and R programming environment [version 3.2.4, (R Development Core Team, 2013)] and R packages corpcor and mvtnorm. Reference files for the known canonical genes on the X chromosome for human genome build 37 were included in the XWAS v2 software package and used to group variants and p-values by gene (Chang et al., 2014; Gao et al., 2015). Bonferroni correction was used to adjust for multiple testing instead of SimpleM, as all genes, unlike SNPs, are independent of each other in the context of association testing and as such the multiple test correction cannot be less than the number of genes tested.
Interaction Analysis
Genome-wide SNP interaction analysis was done using CASSI3 (v2.51). A joint effects model was implemented for a rapid overview of interactions of all variants across the genome (autosome and X chromosome). Variants from significant interactions were reanalyzed using a logistic regression approach with covariate correction, which would not be feasible for a genome-wide interaction analysis as it would be too computationally intensive. As there is no general consensus on the significance threshold for genome wide interaction analysis Bonferroni correction was used in order to avoid potential inflation of false positive results.
Discussion
In this GWAS we investigated TB susceptibility in the admixed SAC population, with a specific focus on sex-bias and the X chromosome. A sex-stratified QC protocol was applied to the data in order to conserve inherent differences between the sexes and all statistical analysis were conducted in a sex-stratified and combined dataset in order to fully assess the impact of sex on TB susceptibility and the male sex-bias it presents with. We found no significant associations on the autosome or X chromosome for both the sex-stratified and combined SNP and gene-based association testing. A few significant interactions were identified, but the impact of these on TB susceptibility is unclear and will require further investigation to validate and functionally verify.
For the combined autosomal SNP based association testing the only potential variant of interest is rs17410035 located in the DROSHA gene (Table 2) which is potentially involved in miRNA biogenesis and could impact TB susceptibility if immune related regulatory miRNA is affected. For the X-linked association testing the top association in males was in an uninformative pseudogene, while the female and combined analysis revealed the same variant, rs768568 located in the TBL1X gene (Table 3). The TBL1X protein has been shown to be a co-activator of NF-kB mediated transcription of cytokine coding genes, but the mechanism of activation is unclear (Park et al., 2016). NF-kB is a vital component of the proinflammatory signaling pathway and is involved in multiple immune pathways including TLRs (Lawrence, 2009), which have previously been shown to influence TB susceptibility (Schurz et al., 2015). Based on this one could extrapolate that variants in the TBL1X gene could affect activation and proinflammatory signaling of NF-kB, which could have a direct effect on the immune system and thus TB susceptibility. The direction of effect for this variant was the same in males and females (Table 3), but was less significant in males probably due to loss of power when analyzing haploid genotypes. For the variants identified in the sex differentiated analysis it is unclear how they could influence TB susceptibility as the top hit is located in a pseudogene. However, the sex differentiated test did reveal just how big the difference in effects can be between the sexes for a specific variant (Table 4). If these variants with opposite effects are not analyzed in a sex-stratified way then the effects would cancel each other out and any information on sex specific effects would be lost. The X-linked gene-based association test revealed no significant associations despite having more power than the SNP based association testing. A possible reason for this could be that Bonferroni correction was used and as this is very conservative possible associations could have been missed. When looking at the most significant associations (Table 5) however, it is unclear how the identified genes could be implicated in TB susceptibility.
The joint effects interaction analysis revealed several significant interactions, but as association results have been previously shown to be severely influenced by admixture (Daya et al., 2014b) only the results for the logistic regression analysis will be discussed here. A few variants were identified in the logistic interaction analysis that could impact TB susceptibility (Table 6). URI1 (rs1118924) is activated by TNFα and is involved in the NF-kB/CSN2/Snail pathway, CIITA (rs6498130) impacts expression of MHC class II proteins and receptors and rs35996537 (RNF125) and rs2040739 (RNF126) are both E3 ubiquitin ligase proteins which affect a multitude of cellular functions, such as apoptosis (Song et al., 2005) and protein degradation (Shin et al., 2015). NF-kB, TNFα, MHC class II, E3 ligases, apoptosis and T-cells have all been implicated in TB susceptibility and could collectively contribute by influencing the immune response (Hirsch et al., 1999, 2005; Torres et al., 2006; Fallahi-Sichani et al., 2012; Bai et al., 2013; Parandhaman and Narayanan, 2014; Shin et al., 2015; Franco et al., 2017). As TB is a complex disease all potential influential factors need to be considered and as such the interaction analysis cannot be ignored. Shortcomings of the interaction analysis are that they are very computationally intensive and suffer from a massive multiple test correction burden. Future research should thus focus on ways to prioritize variants for interaction analysis to decrease computation time as well as have sufficient sample size to minimize multiple test correction burden.
A previous GWAS in the SAC population found a significant association with TB susceptibility in the WT1 gene (rs2057178, OR = 0.62, p-value = 2.71e-6) (Chimusa et al., 2014). This association did not reach genome-wide significance in our study (OR = 0.75, p-value = 0.049). At the time of the GWAS by Chimusa et al., 2014) there were few African and KhoeSan (only 6 KhoeSan) individuals in the reference data used for imputation and the accuracy of imputation in this population was not known. As the identified variant (rs2057178) was imputed into the data it should have been validated in the SAC population using an appropriate genotyping approach. Secondly although the variant reached a significance threshold for the number of variants tested it did not reach genome wide significance threshold of 5.0e-8 (Panagiotou and Ioannidis, 2012). Finally, the GWAS performed by Chimusa et al. (2014) only contained 91 control individuals compared to 642 cases, which could affect the power of the study. Chimusa et al. (2014) were unable to replicate previous associations identified in the X-linked TLR8 gene (Davila et al., 2008). The two TLR8 variants in our data, rs3764880 (OR = 1.73, p-value = 3.1e-4) and rs3761624 (OR = 1.70, p-value = 3.94e-4) also did not show significant associations. While the haploid genotypes in males contributes to this, a second influential factor could be admixture. Chimusa et al. (2014) did not perform X chromosome specific admixture analysis, which could affect association testing of X-linked genes. Furthermore, only six KhoeSan reference individuals were available, which could affect the accuracy of admixture inference and severely affect the results. For our study 307 KhoeSan individuals were available, improving the admixture inference and could explain why stronger effects (higher OR) were detected for the TLR8 variants when compared to Chimusa et al. (2014). It is also important to note that using global ancestry components as covariates does not correct for ancestry at any specific locus and as a result each locus in this population could have up to five different ancestries. This could greatly reduce power and contribute to the lack of replication between studies. In order to address this future studies could incorporate local ancestry inference into the analysis in order to determine the number of ancestries at a locus of interest. Other candidate genes identified in previous GWAS studies were also separately analyzed here, but associations did not replicate (Online Supplementary Data Sheet S2).
We did not find any significant associations with TB susceptibility, but highlight the need for sex-stratified analysis. Closer inspection of the data revealed that a large number of SNPs with opposite direction of effects for not only the X chromosome, but the autosome too. Sex specific effects has previously been reported for autosomal variants associated with pulmonary function in asthma (Berhane et al., 2000). In the SAC population these opposite effects have previously been observed for X-linked variants in the TLR8 gene (Daya et al., 2013) and the same is observed in this study. Sex-stratified analysis should therefore be included in association studies and incorporated in the study design. This can be done by keeping the male to female ratio balanced in the cases and controls. It would also be prudent to do the power calculation for the males and females separately. This will ensure sufficient power for sex-stratified analysis and could elucidate informative sex specific effects. This study was done in a 5-way admixed population. As was observed for the interaction analysis including admixture components significantly changes the association results. Furthermore it was observed (unpublished results) that the ancestral distribution between the X chromosome and autosome are different (Figure 1), which is an indication of sex-biased admixture (Goldberg and Rosenberg, 2015; Shringarpure et al., 2016) and highlights the importance of including X chromosome admixture components for X-linked and sex-bias analysis. It is important to note here that the ancestral components in the SAC present with a very wide range (Figure 1) and all this variability could affect the power of association studies. It is therefore desirable to increase the sample size when analyzing admixed individuals. Alternatively, a meta-analysis can be conducted, including data from all five ancestral populations, or local ancestry inference could be included in the analysis.