A FLOT1 host regulatory allele is associated with a recently expanded Mtb clade in patients with tuberculosis

Genome-to-genome study identifies human-Mtb genomic interactions

We identified an association between an intronic human variant, rs3130660, located in the 6p21 region on chromosome 6 and a phylogenetic marker (Position 271640) of a subclade of L2 Mtb (OR = 1.39, 95% CI: 1.33-1.44, P = 9.24 × 10⁻¹⁵, Fig. 2B). Individuals with each rs3130660 A allele were ∼40% more likely to be infected with Mtb strains carrying the interacting Mtb variant marker. We noted that the associated Mtb variant is in tight linkage with 78 of the 2,298 common Mtb variants (Pearson r² > 0.8, Table S2).

To better understand these 79 Mtb variants, we constructed a maximum likelihood phylogenetic tree using WGS of 1,556 collected isolates. In this analysis, we looked at 13,981 Mtb variants which included both the 2,298 common variants along with 11,683 variants that passed stringent quality control. We classified the strains into well-defined Mtb lineages, including L2 and L4.1-8 using previously defined lineage defining markers(16). We observed that previously defined lineages fell within distinct branches of the phylogenetic tree (Fig. 2C). Strikingly, we observed that all 79 Mtb SNPs were present in the same L2 subclade of the phylogenetic tree. We henceforth define the clade by the presence of the g2g Mtb variant (Position 271640) as the g2g-L2 clade.

Human allele associated with Mtb diversity colocalizes with FLOT1 expression in lung and other tissue types

The strongest association in the human genome is an intronic variant (rs3130660) of the FLOT1 gene (Fig. 3A). FLOT1 is a lipid raft-associated scaffolding protein that plays a role in membrane trafficking and phagosome maturation(17, 18). We investigated whether this SNP modulates the expression levels of nearby genes by performing an expression quantitative trait loci (eQTL) analysis. Using the Genotype Tissue Expression (GTEx release v8(19)) database in lung, we observed that this SNP was strongly associated with FLOT1 expression (effect size = 0.28, P = 2.22 × 10⁻¹⁶, Fig. 3B, Fig. S3). Each rs3130660 A allele corresponds to an 1.32-fold increase in FLOT1 expression in the lung. The most strongly associated variant with FLOT1 expression in GTEx was rs3132610 which is in tight linkage with rs3130660 (r² = 0.90 in 1KG European population(20), r² = 0.77 in this Peruvian cohort). To understand whether the rs3130660 allele corresponds to the signal reported in GTEx, we applied a colocalization analysis using the coloc software(21). We note that there are differences in genetic ancestry between our study and GTEx, which are predominantly individuals of European ancestry, and this might reduce the power to detect colocalizing signals across datasets(22). We also note that there are multiple independent alleles influencing FLOT1 expression(19), which might further reduce power to detect colocalization (Fig. 3B). Despite this, we observed strong evidence of colocalization in the lung (posterior probability = 0.98). Consistent with the colocalization analysis, when examining 14,976 SNPs overlapping with eQTL summary statistics of FLOT1, we observed strong correlation between association statistics of FLOT1 expression in lung and infection with g2g-L2 strain (Pearson r = 0.89, Fig. S4). We then looked at other cell and tissue types and assessed colocalization signals between FLOT1 and other 49 protein coding genes within the 700kb window of rs3130660 (Fig. 3C, Table S3). Using 108 RNA sequencing datasets included in the eQTL catalog release 4(23), we found evidence of colocalization between the identified interacting host SNP and FLOT1 eQTL association in 15 cell and tissue types (posterior probability > 0.8), including T cells (posterior probability = 0.99) and adipose (posterior probability 0.96, Fig. 3D, Fig. S5). These results provided strong evidence for colocalization between the host variant associated with the g2g-L2 strain and host variants associated with FLOT1 expression levels in multiple tissues and cell-types. This suggests that FLOT1 may be the causal human gene interacting with the Mtb genome.

• Download figure
• Open in new tab

FIG 3. Association in the human host FLOT1 locus.

(A) Manhattan plot of the g2g association analysis in the FLOT1 locus. The y-axis shows the -log₁₀(p-values) obtained using a linear mixed model. Point colors represent the LD r² with the lead host variant (rs3130660) in this Peruvian cohort. (B) Manhattan plot of the FLOT1 eQTL. The x-axis shows the human genomic positions. The y-axis shows the -log₁₀(p-values) for the eQTL association between host variants and FLOT1 expression level in the lung (obtained in GTEx). Point colors represent the LD r² between either of the two eQTL SNPs of FLOT1 expression (rs3132610, circled in red and rs9461629 circled in blue) (C) Genomic locations of all genes in the FLOT1 locus in GRCh38. (D) A barplot of all posterior probabilities of shared causal variants between the g2g and FLOT1 eQTL signals in 13 cell and tissue types selected among 108 bulk RNA-sequencing datasets. Completed colocalization results in 108 datasets are provided in Table S3.

To test whether this result reflected social mixing patterns within the community or a true biological association, we repeated our analysis adjusting for year of TB diagnosis and Mtb population structure (as reflected by principal components), neither of which substantially affected our results (Fig. S6). Notably, the prevalence of the g2g-L2 clade remained consistent over time (Fig. S7) and over geographical space (Fig. S8), reinforcing our conclusion that the reported signals are independent from these covariates.

We considered the possibility that the observed association between a host allele and the Mtb g2g-L2 subclade was driven by host ancestry (Fig. S9) or Mtb lineage (Fig. S10). We examined host alleles for associations to the two common Mtb lineages in Peru (L2 and L4), and found the strongest association was between the same FLOT1 allele and the L2 lineage, but that this association was substantially weaker than with the g2g-L2 lineage (rs3130660, OR = 1.37, 95% CI: 1.29-1.44, P = 7.04 × 10⁻⁹, Fig. S11). A conditional analysis including g2g-L2 obviated this association (P = 0.97). In contrast, adding L2 as a covariate to our original model did not obviate the association wit g2g-L2 (P = 8.15 × 10⁻⁸). This suggested to us that the observed association between rs3130660 and L2 was driven by the g2g-L2 subclade rather than L2 more broadly. We also tested the effect of five human host ancestry principal components on Mtb genomic variants (1,267 tests) and the association of host ancestry principal components on two Mtb lineages. We observed no significant associations, suggesting that host genetic ancestry is not associated with the variations in the Mtb genome in Peru.

Mtb g2g-L2 subclade is recently expanded in Peru and transmits differently than other L2 strains

To interrogate the population dynamics of local Mtb L2 strains identified in the g2g analysis, we used 178 L2 isolates from this study and 77 previously collected L2 strains from the same population to reconstruct a maximum likelihood phylogenetic tree(24, 25). We showed that in this population there are three major clades: clade A consisting of the g2g-L2 isolates, and clade B and C (Fig. 4A). We merged the WGS of these L2 clade isolates with 1,000 global L2 isolates that were previously published in 112 bioprojects (Table S4). We found the L2 strains from Peru formed subclades that were largely separated from other global L2 isolates (Fig. S12). This pattern of restriction suggests that Mtb L2 strains were introduced to Peru and then diversified locally.

• Download figure
• Open in new tab

FIG 4. Genetic structure of the Mtb clade associated with the human FLOT1 locus.

(A) A phylogenetic tree of L2 constructed with 255 L2 Peruvian isolates. The g2g-L2 clade identified via the g2g analysis is highlighted in red (Clade-A), two other Peruvian subclades of L2 are highlighted in blue (Clade-B) and green (Clade-C) respectively. Estimated origin time (median value) for the ancestor strains and cluster rate when using 6-SNP distance of each marked L2 clade (Clade-A, B and C) are listed. (B) Histogram of pairwise minimum SNP distance to the closest neighbors within the marked clades. (C) Comparison of transmission cluster rate between the three marked clades when using 6-SNP (dark blue) and 12-SNP (light blue) distance as the threshold.

We next performed a Bayesian coalescent analysis to estimate the introduction times of these three L2 subclades(26). These results suggest that the g2g-L2 (clade A) was introduced to Peru more recently than the other two L2 subclades B and C (61.9 versus 138.7 and 112.3 years ago, Fig. 4A). We then compared the minimum SNP distance between these isolates and their closest neighbors in our WGS dataset to estimate an upper bound of the divergence time since the last transmission event for each TB case. The minimum SNP distances were significantly shorter for the g2g-L2 strains (median 1 SNPs) than isolates in clades B and C (median 8 and 18 SNPs for clade B and C respectively, P = 6.85 × 10⁻¹³ using the Kolmogorov-Smirnov test, Fig. 4B). To further compare the transmissibility to these clades, we assessed the frequency of L2 strains within our cohort. We expect that locally diversified strains which are introduced earlier would have a larger population size(27, 28). Instead, we found that the most recently introduced g2g-L2 strains have a larger than expected population size. Among 255 participants infected with L2 in our cohort, we observed 150 of 255 (58.8%) of them were infected by g2g-L2 strains (clade A), which was much higher than the number of infections caused by clade B or C (18.0% and 6.7%). These observations suggest that the g2g-L2 strains may have undergone a recent, rapid expansion.

Our findings raise the possibilities that g2g-L2 strains are more transmissible than other L2 clades or are more adept at causing disease in the Peruvian population starting in the 1950s after their introduction. To better understand the differences in the transmission dynamics among the three local clades, we explored the distribution of potential transmission clusters where two TB cases were considered to belong to the same transmission cluster if they were separated by fewer than 6 SNPs. We found among 150 participants infected with the g2g-L2 clade, 139 (92.7%) had a SNP distance smaller than 6-SNPs compared to closest neighbors, suggesting they are potentially linked by direct transmission. This cluster rate is much higher than that of other L2 strains collected from the same region (Fig. S13) (92.7% versus 40%, 42/105). We reached the same conclusion when using a less stringent 12-SNP cut-off (Fig. 4C). These data suggest that the g2g-L2 clade is associated with a higher level of recent transmission in the Peruvian population.

Mtb g2g-L2 subclade displays differing redox metabolism from other L2 strains

To identify biological features that distinguish g2g-L2 strains, we phenotyped 23 g2g and 11 control L2 strains using a high throughput imaging-based phenotyping platform(29) (Fig. 5A, Fig. S14). We assessed seven functional phenotypes including cell morphology (length, width and area), DNA content with DAPI staining, the pattern and rate of peptidoglycan synthesis using NBD-amino-D alanine (NADA) incorporation, redox state using cyan fluorescent protein (CFP) autofluorescence, and intrabacterial lipid inclusions accumulation with Nile Red-staining. Among these phenotypes, we observed a significantly higher CFP signal in g2g compared to other L2 strains (Wilcoxon test, P =2.0 × 10⁻⁴, Fig. 5B-C), indicative of a shift in the redox state of F₄₂₀ towards a more oxidized state. F₄₂₀ is a redox active enzyme cofactor found in actinobacteria and methanogens that mediates low-reduction potential metabolic reactions(38). To determine if this functional difference in F₄₂₀ results in biologically meaningful consequences, we assessed susceptibility to pretomanid, an antibiotic pro-drug that requires the reduced form of F₄₂₀ to be activated(30). The g2g-L2 strains were more resistant to pretomanid, consistent with the measured shift towards oxidized F₄₂₀ (Wilcox test, P = 0.032, Fig. 5D).

• Download figure
• Open in new tab

FIG 5. Functional characterization of g2g-L2 Mtb strains.

(A) Schematic of the high throughput imaging-based phenotyping platform. 23 g2g-L2 and 11 other L2 Mtb strains were grown in 7H9 to a final OD600 of approximately 1.0. NADA (10uM) was added 16 hours prior to fixation with 2% paraformaldehyde. All samples were analyzed by microscopy to assess the population mean of cell length, cell width, cell area, DNA content using DAPI, intrabacterial lipid inclusions accumulation with Nile Red-staining, the pattern and rate of peptidoglycan synthesis using NADA incorporation, and redox state using CFP autofluorescence (oxidized form of F420 enzyme cofactor). (B) Violin plots for each phenotype described in (A), showing the distribution of each feature between the g2g and other L2 strains. Each point shows a representative sample for each individual strain; samples were assayed at minimum in duplicate (Fig. S14). The dotted lines inside each plot indicate the first and third quartiles; the dashed line indicates the median. Each point represents the average of two independent replicates per strain. A Wilcoxon test was performed to measure the difference between g2g and other L2 strains. The red dotted line indicates the Bonferroni corrected significance threshold after multiple testing (-log₁₀(0.05/7)). (C) Representative images of Nile Red, NADA, CFP, and DAPI signals in two representative g2g-L2 strains and two non-g2g-L2 strains are shown (strains marked as enlarged black squares in (B)). (D) Based on a down-selection from the microscopy data, 5 g2g and 5 other L2 Mtb strains (two replicates per strain) were incubated with pretomanid (0.078uM) for 7 days. Resazurin was added to each well, and cultures were incubated an additional 24 hours. Variability after pretomanid exposure is measured by percent reduction in resazurin as calculated from drug-free controls. We used a Wilcoxon test for comparing the means of g2g and other L2 strains (P=0.032).

We examined putative Mtb genes that might be driving this functional phenotype. To this end, we identified nonsynonymous changes in 49 variants among 79 common Mtb variants defining g2g-L2 (Table S2). We found many of these occurred in genes that overlap with redox-related pathways (Fig. S15). In a gene set enrichment analysis, g2g-L2 was enriched for nonsynonymous mutations associated with flavin adenine dinucleotide (FAD) binding compared to all non-g2g-L2 nonsynonymous mutations (Fisher’s exact test, P = 0.008). FAD is a redox coenzyme of several important reactions in metabolism. Genes that overlap with the FAD binding pathway (GO:0050660) include trxB2 (Thr2Asn), fadE3 (Tyr59His), fadE7 (Met364Ile), fadE34 (Ala371Val), and fadE35 (Asp417Asn). Three of these five mutations were predicted to attenuate protein function using SIFT(31) (SIFT score < 0.05, Table S2). Intriguingly we also observed a nonsynonymous mutation in ctpC (Asp633Asn) that has an attenuated function (SIFT score = 0.02). Previous study showed that cells lacking ctpC are both variable and sensitive to superoxide(32). Taken together, we found a cytological feature (oxidized form of F₄₂₀) that can differentiate g2g-L2 strains from other co-circulating L2 strains, which can manifest as reduced susceptibility to the anti-TB drug pretomanid and concordant of mutations that accumulated by g2g-L2.

Ethics Statement

We recruited 1,632 subjects from a large catchment area of Lima, Peru that included 20 urban districts and ∼3.3 million residents to donate a blood sample for use in our study. We obtained written informed consent from all the participants. The study protocol was approved by the Institutional Review Board of Harvard School of Public Health and by the Research Ethics Committee of the National Institute of Health of Peru.

Participant enrollment and follow-up

The study design and methods were previously described in detail(24, 44). Briefly, over the four-year period from 2009-2012, we identified patients ≥ 15 years of age who had received a diagnosis of pulmonary TB at any of 106 participating health centers. We confirmed the microbiological status of their disease with either a positive sputum smear or mycobacterial culture. We also recorded the index patients’ baseline smear status, HIV status, and drug-resistance profiles. Index cases with HIV infection or infected with multiple Mtb strains were excluded from the analyses.

Host genotyping, quality control and imputation

DNA samples were genotyped using a customized genotyping array (LIMAArray) based on whole-exome sequencing data from 116 active TB cases to optimize the capture of genome-wide genetic variation in Peruvian individuals(10). We excluded samples that were missing more than 5% of the genotype data, had an excess of heterozygous genotypes, and/or duplicated with identity-by-state > 0.9. We also excluded variants with a call rate less than 95%, with duplicated position markers, minor allele frequency < 1%, and marked deviation from Hardy-Weinberg equilibrium (excluded if P < 10⁻⁵). After these quality control metrics, the genotype data were pre-phased using SHAPEIT2(45). We then used IMPUTE2(46) to impute genotypes at untyped genetic variants using the 1000 Genomes Project Phase 3(20) as the reference panel.

We imputed eight classical HLA genes HLA-A, -B, -C, -DQA1, -DQB1, -DRB1, -DPA1, and -DPB1, amino acids and intergenic variants using HLA-TAPAS using a multi-ancestry reference panel which contains data from 21,546 unrelated individuals(47).

Mycobacterium tuberculosis variant calling and phylogeny construction

Mtb DNA was extracted from cryopreserved cultures of sputum samples. Each sample was thawed and subcultured on LJ and a big loop of colonies were lysed with lysozyme and proteinase K to obtain DNA using CTAB/Chloroform extraction and ethanol precipitation. Samples were sequenced on an Illumina HiSeq 2500 or 4000 sequencer yielding paired-end reads of length 125bp which were aligned to the reference assembly, H37Rv NC_000962.3 (GenBank accession CP003248). We called variants using Pilon (v1.22)(48). Genome coverage was assessed using Samtools (v1.10). We excluded isolates with evidence of mixed infection using the barcode method(16). We assigned a variant call as missing if the valid depth of coverage at a specific variant was less than 12 reads, if the mean read mapping quality at the site did not reach 10, or if none of the alternative alleles account for at least 85% of the valid coverage.

The phylogenetic tree was constructed based on the WGS Mtb alignment. Variants occurring in genes with repetitive elements such as transposases, proline-glutamate (PE) or proline-proline-glutamate (PPE)(49) were excluded to avoid any inaccuracies in the alignment. After applying these filters, 13,981 Mtb variants were used to conduct a genetic distance matrix. We built a phylogenetic tree using distance-based neighbor-joining coded in the ape package(50).

To characterize the Peruvian L2 Mtb strains in the context of global strains, we randomly selected 1,000 whole-genome sequenced L2 strains from published studies in 112 bioprojects to represent the major phylogenetic structures of L2. We used these global strains together with 255 L2 Peruvian strains to reconstruct the global L2 phylogenetic tree.

We estimated divergence times via BEAST v1.8.0(26), using an uncorrelated lognormal relaxed clock that allows for tree branches to evolve at different rates. The XML input file was modified to specify the number of invariant sites in the Mtb genomes. For the Mtb substitution rate, we used a normal distribution with a mean of 4.6 × 10⁻⁸ substitution per genome per site per year (3.0 × 10⁻⁸ to 6.2 × 10⁻⁸, 95% highest polar density interval), which was calibrated by ancient DNA samples(51, 52). An uncorrelated log-normal distribution was used for the substitution rate and a constant population size for the tree priors. We ran three chains of 5 × 10⁻⁷ generations and sampled every 10,000 generations to assure independent convergence of the chains; we discarded the first 10% as a burn-in. Convergence was assessed using Tracer (v1.7.0)(53), ensuring all relevant parameters reached an effective sample size > 100.

Transmission cluster analysis

To define transmission clusters of the collected L2 Mtb isolates, we applied two SNP thresholds (6 and 12) to separate a patient isolate from that of at least one other patient in the cluster. The 6-SNP threshold was defined as the range of SNP distances between paired isolates of the same strain obtained at different times from relapsed TB patients(54). The 12-SNP threshold was a previously defined upper threshold of genomic relatedness noted within human hosts and between epidemiologically related human hosts(55).

Genome-to-genome association analysis

To systematically test for interaction between human and Mtb genomic variants at the genome-wide level, we performed mixed logistic regression modeling using GEMMA(56) version 0.94.1 assuming an additive genetic model for each host and bacterial variant. For each bacterial variant j (j = 1,…, 1,267) and host variant i (i = 1,…, 676,110), we test where y_j is an N-vector of binary labels indicating the absence or presence of a bacterial variant in the host samples i (n = 1,556); x_i is an N-vector of host genotypes; β_ji is the additive effect of the host variant i on the bacterial variant j; c_j is the j^th column of an N×C matrix of covariates (age and sex) including a column of 1s; and ω_j is a c-vector of the corresponding coefficients including the intercept. We used genetic relatedness matrix (GRM) as a random effect to correct for cryptic relatedness and population stratification between collected host samples. The GRM was obtained from an LD-pruned (r² < 0.2) genotype with minor allele frequency >1%.

Permutation Strategy

To determine an appropriate empirical genome-wide significance threshold for the g2g analysis accounting for population structure both in the Mtb and human genome, we conducted 200 simulations for each g2g analysis using a permutation procedure. To preserve the phenotype (defined by the Mtb variants) structure, for each interaction, we randomly permuted the presence/absence status for each of the 1,267 Mtb variants within the 1,556 individuals included in the study. In total, we ran 1,267 × 200 = 253,400 association tests to obtain an empirical genome-wide significance threshold.

We measured the distributions of the minimum p-values of the variants (P_min) for each interaction of the permutation. We defined an empirical genome-wide significance threshold, - log₁₀(P_sig), as the 95th percentile (1-α) of -log₁₀(P_min).

Colocalization analysis using eQTLs

We integrated our GWAS results with cis-eQTL data using a Bayesian method (coloc v3.2-1)(57). We evaluated whether the GWAS and eQTL associations best fit a model in which the associations are due to a single shared variant (summarized by the posterior probability) as opposed to different regulating variants. We used gene expression datasets from the eQTL catalog(23) release 4. We tested pairwise colocalization between 108 bulk RNA-sequencing expression datasets in 69 distinct cell and tissue types and GWAS from the most significant genome-to-genome pair (position 271640 as phenotype). We used GWAS and all variant-gene cis-eQTL associations tested in each tissue (including non-significant associations) of genes within a 800kb window around the top association host SNP (rs3130660). A posterior probability of colocalization >0.8 was considered to be a strong evidence for a causal gene.

Bacterial strains and culture conditions

All 23 g2g and 11 non-g2g Mtb L2 strains were drug sensitive and HIV negative. They were grown shaking at 37°C. Cultures were grown in 7H9 media (Middlebrook 7H9 salts with 0.2% glycerol, 10% OADC [oleic acid, albumin, dextrose, catalase]). Cell density (OD600) was determined by a spectrophotometer.

Microscopy imaging and analysis

Mtb cultures were grown to OD600 of approximately 1.0, then fixed with 2% paraformaldehyde (PFA) for 1 hour. The fluorescent D-amino acid NADA (3-[(7-Nitro-2,1,3-benzoxadiazol-4-yl)amino]-D alanine hydrochloride) was added at a final concentration of 1mM 16 hours prior to fixation. All samples were seeded onto molded 1.8% agarose in phosphate buffered saline (PBS) with DAPI (2.5μg/mL) and Nile red (0.1μg/mL). Samples were incubated at 37°C for one hour prior to imaging.

Samples were imaged with a Plan Apo 100× 1.45 NA objective using a Nikon Ti-E inverted, widefield microscope equipped with a Nikon Perfect Focus system with a Piezo Z drive motor, Andor Zyla sCMOS camera, and NIS Elements (v4.5). Semi-automated imaging was carried out using a customized Nikon JOBS script to locate imaging fields of interest, 24 images were taken for each strain. Cell segmentation and quantification was performed using our previously published pipeline, MOMIA(29).

Minimum inhibitory concentration (MIC) determination

Pretomanid (PA-824) MIC value for each strain was determined following the microplate-based Alamar Blue assay method as previously described(58). 96-well plates containing two-fold serial dilutions of the antibiotic Pretomanid (PA-824) in 7H9 media were prepared, creating a concentration gradient from 40 to 0.039 µM. A drug-free control well was included. All Mtb strains were grown to mid-log phase and diluted to OD600 = 0.003 into each well of the prepared plates. Plates were incubated at 37°C without shaking in sealed plates. After 7 days, a freshly prepared mixture of 0.15mg/ml resazurin (Sigma Aldrich) diluted in PBS and 20%

Tween-80 was added to each well. Plates were incubated for 24 hours, after which absorbance was read at 570 nm, with a reference wavelength of 600 nm. Percent reduction of resazurin was normalized for each strain by subtracting the OD570 of the minimum value wells (40µM drug condition) and scaling by the reduction measurement of the no-drug control. Data are representative of two independent replicates per strain.

Pathway analysis of interacting Mtb genes

We used Variant Effect Predictor(59) (v.105) to first annotate the Mtb genome. Among 79 highly correlated (Pearson r²>0.8) common Mtb variants that are associated with the human genome, 49 were annotated as nonsynonymous (missense) variants. We defined Mtb genes that have one of these 49 missense variants as g2g-L2 genes. We used the DAVID(60) Functional Annotation Clustering function to group g2g-L2 genes into a similar functional network. We used PANTHER(61) to test for functional overlaps with known pathways. We calculated a Fisher’s exact p-value by constructing a 2×2 contingency table with g2g-L2 genes (n = 48) and all other Mtb genes with a common missense mutation in the g2g analysis (n = 898) overlapping the flavin adenine dinucleotide binding pathway included in the gene ontology category (GO:0050660).