Mycobacterium abscessus pathogenesis identified by phenogenomic analyses

Abstract

The medical and scientific response to emerging and established pathogens is often severely hampered by ignorance of the genetic determinants of virulence, drug resistance and clinical outcomes that could be used to identify therapeutic drug targets and forecast patient trajectories. Taking the newly emergent multidrug-resistant bacteria Mycobacterium abscessus as an example, we show that combining high-dimensional phenotyping with whole-genome sequencing in a phenogenomic analysis can rapidly reveal actionable systems-level insights into bacterial pathobiology. Through phenotyping of 331 clinical isolates, we discovered three distinct clusters of isolates, each with different virulence traits and associated with a different clinical outcome. We combined genome-wide association studies with proteome-wide computational structural modelling to define likely causal variants, and employed direct coupling analysis to identify co-evolving, and therefore potentially epistatic, gene networks. We then used in vivo CRISPR-based silencing to validate our findings and discover clinically relevant M. abscessus virulence factors including a secretion system, thus illustrating how phenogenomics can reveal critical pathways within emerging pathogenic bacteria.

Fig. 1. Multidimensional phenotyping of M. abscessus.

a, Phenotypic variability of clinical M. abscessus isolates was assessed across multiple dimensions (described in Methods) including: planktonic growth (assessed by serial OD measurement) in a range of different carbon sources; MIC of a range of clinically relevant antibiotics assessed on day 3 (MIC early) and day 11 (MIC late) to quantify intrinsic and inducible drug resistance; macrophage infection (4 h post infection), intracellular replication (2 days post infection) and death (2 days post infection) quantified using high-content imaging of differentiated THP-1 cells incubated with tdTomato-expressing clinical isolates; survival and immune response of Drosophila melanogaster infected with clinical isolates; and clinical outcomes (lung function decline and clearance of M. abscessus from sputum samples) of infected patients. Ami, amikacin; Cla, clarithromycin; Clo, clofazimine; FEV1, forced expiratory volume; Fox, cefoxitin; Lin, linezolid; OD-AUC, area under the OD curve; qPCR, quantitative polymerase chain reaction. b, Pearson correlation coefficients within and across phenotypic groups shown as a matrix, with two-sided non-significant (unadjusted P > 0.05) associations shown in white. Source data

Fig. 2. Phenotypic groups.

a,b, Clustering of clinical isolates, using k-means (pink, blue, green) or t-SNE, based on experimentally observed phenotypes only, reveals similar groups that differ in their clinical outcomes. c, Distribution of specific phenotypes across the three phenotypic groups (bacterial growth: P = 3.6 × 10⁻⁷; clarithromycin resistance: P = 2.5 × 10⁻¹¹⁵; intracellular bacterial replication P = 1.8 × 10⁻¹⁰; macrophage death: P = 5.1 × 10⁻¹⁰; Drosophila survival: P = 8.5 × 10⁻¹¹). AU, arbitrary units. d, Maximum likelihood phylogenetic tree of M. abscessus subspecies and corresponding phenotypic groups. e, Nearest phylogenetic neighbours most commonly belong to the same phenotypic group. P values were calculated using a two-tailed chi-squared test or one-way analysis of variance, as appropriate. Source data

Fig. 3. Integrating computational structural modelling into GWAS.

a, Genome-wide associations were performed for all phenotypes with the top variants extracted (up to five per association) and ordered using hierarchical clustering (red, linear model; blue, mixed model; P values were calculated using a two-sided Wald test). Pairwise R² measurements of the identified genetic variants (grey scale) show extensive genome-wide linkage (LD). b, To identify causal variants and overcome genome-wide linkage, the functional impacts of genetic variants were classified as having high effects (large deletions, frameshifts, start/stop alterations; red), moderate effects (inframe insertions/deletions; blue and green) and low effects (synonymous and intergenic variants; grey). The impacts of missense mutations were estimated using proteome-wide computational structural modelling with variants considered as having high (red), moderate (blue) or low (green) functional effects based on terciles of the change in protein stability, estimated using mCSM. c, Manhattan plot of the mixed model GWAS analysis of 264,122 genetic variants for intracellular M. abscessus replication (two-sided Wald test) with the threshold for multiple hypothesis testing (black line). Several loci in mbtD, including four missense mutations, were identified as potential mechanisms relevant for intracellular M. abscessus survival. (Inset) Three-dimensional structural model of MbtD with the high-effect missense mutation Ile256Thr shown in red. d, MbtD knockout mutants complemented with wild-type or identified MbtD variants had similar growth rates in broth culture but replicated differently within THP-1 cells. Experiments were performed in triplicate on at least three separate occasions and are presented as mean ± s.e.m. Conditions were compared with a two-sided unpaired t-test. Source data

Fig. 4. CC-DCA for assessing genome-wide epistasis.

CC-DCA was used to identify co-evolving variants among ~10¹² potential variant combinations of 2,366 clinical M. abscessus isolates. a, Circos plot of the M. abscessus chromosome showing the 100,000 top variant-to-variant couplings with a distance of >100 bp (black lines), coupling density (green; range 0–56,307 couplings per 5 kb) and SNP density (red; range 14–1,961 SNPs per 5 kb). b, Significant variant–variant couplings identified through DCA were pooled to gene–gene couplings. Whereas variant–variant couplings indicate the total number of co-evolutionary signals within a single gene, gene–gene couplings reflect the number of putative gene interactions. NR, non-ribosomal. c, Networks of co-evolving (and therefore probably epistatic) genes based on the 1,000 strongest DCA-derived gene–gene couplings ranked by coupling number, colour coded by functional class. The strength and number of couplings are shown by edge colour and thickness respectively. d, Example of a highly coupled gene network (highlighted by a circle in c) involving components of the mycobactin biosynthesis pathway e, CRISPR-induced transcriptional repression of several genes within this cluster demonstrates impaired mycobacterial survival within macrophages. Experiments were done in triplicate (with three guides per gene) on at least three separate occasions, are presented as mean ± s.e.m., and compared with the empty vector control using a two-sided unpaired t-test, **P < 0.001, ***P < 0.0001 (MbtD knockdown d2: P = 0.0002; MbtE kd d1: P = 0.0006, d2: P = 4.3 × 10⁻⁵; MbtF kd d1: P = 0.0001, d2: P = 1.3 × 10⁻⁵; MAB_4071 kd d1: P = 0.008; MAB_4072 kd d1: P = 0.0002). Source data

Fig. 5. Integrating GWAS and DCA to reveal the genetic networks of in vivo virulence in M. abscessus.

a, Representative image of Drosophila melanogaster infected with M. abscessus (magenta) resembles mycobacterial infection in other organisms (independently repeated over five times), with infection of phagocytes (green) and formation of granuloma-like structures (inset). b, Genome-wide association (using a linear model and applying Wald test statistics) reveals a putative secretion system protein and a peptide synthetase to be highly associated with Drosophila survival. The black horizontal line marks the multiple hypothesis testing threshold based on the number of independent variants. HP, hypothetical protein. c, Both variants align to clinical isolates with long survival, including a dominant circulating clone, within the subspecies M. a. abscessus. d, Deletion in MAB_0471 was associated with persistent respiratory infection in CF patients (two-sided unpaired t-test). e, CRISPR–dCas9 knockdown of MAB_0471 and MAB_3317 (unlike the essential gene yidC) did not affect growth in liquid culture (left) but in vivo silencing did lead to prolonged survival of infected Drosophila, as shown by Kaplan–Meier survival analysis (log-rank test, P = 7.6 × 10⁻¹⁷) generated from data from at least 18 infected flies per bacterial strain. f, Epistatic gene network, derived from DCA outputs, revealed direct coupling of MAB_0471 with other putative secretion system proteins including MAB_0472 and a distant connection to the peptide synthetase MAB_3317. g, In vivo silencing of MAB_0472 replicated virulence attenuation (log-rank test, P = 3.6 × 10⁻¹²). Source data