Genome-wide gene expression tuning reveals diverse vulnerabilities of M. tuberculosis

Abstract

Antibacterial agents target the products of essential genes but rarely achieve complete target inhibition. Thus, the all-or-none definition of essentiality afforded by traditional genetic approaches fails to discern the most attractive bacterial targets: those whose incomplete inhibition results in major fitness costs. In contrast, gene "vulnerability" is a continuous, quantifiable trait that relates the magnitude of gene inhibition to the effect on bacterial fitness. We developed a CRISPR interference-based functional genomics method to systematically titrate gene expression in Mycobacterium tuberculosis (Mtb) and monitor fitness outcomes. We identified highly vulnerable genes in various processes, including novel targets unexplored for drug discovery. Equally important, we identified invulnerable essential genes, potentially explaining failed drug discovery efforts. Comparison of vulnerability between the reference and a hypervirulent Mtb isolate revealed incomplete conservation of vulnerability and that differential vulnerability can predict differential antibacterial susceptibility. Our results quantitatively redefine essential bacterial processes and identify high-value targets for drug development.

Keywords: Bayes Theorem; CRISPR-Cas Systems; Drug Development; Essential genes; Mass Spectrometry; Mycobacterium smegmatis; Mycobacterium tuberculosis; Vulnerability.

Graphical abstract

Figure 1

Genome-scale CRISPRi fitness profiling in Mtb (A) Experimental design to quantify Mtb gene vulnerability. (i) The Mtb CRISPRi library was built by cloning an sgRNA oligo array into an anhydrotetracycline (ATc)-inducible Sth1dCas9 vector. The library was designed to target all possible Mtb genes with sgRNAs of varying predicted knockdown efficiencies. (ii) Cultures were passaged for approximately 30 generations in the presence (CRISPRi on) or absence of ATc. At the indicated time points, genomic DNA was harvested and sgRNA targeting sequences amplified for next-generation sequencing. (iii) The relative fitness of individual strains was quantified by the sgRNA log2 fold change (L2FC) over time (+ATc/–ATc). Relative fitness values were then used to quantify three parameters that define target vulnerability: (1) maximum fitness cost, (2) sensitivity to partial knockdown, and (3) the phenotypic lag between the timing of CRISPRi induction and onset of a fitness defect. (B) Boxen plots (mean and quantiles) comparing time-dependent changes in L2FC values of sgRNAs targeting genes defined as Essential (n = 63,867) or Non-Essential (n = 29,609) by TnSeq and of control Non-Targeting sgRNAs (n = 1,658). sgRNAs targeting TnSeq Uncertain genes (n = 563) are not shown. (C) Hierarchical clustering of gene level depletion from the experiment described in (A). Each row represents a single targeted Mtb gene. (D) Bar chart showing the overlap between gene calls by TnSeq and CRISPRi. 42 genes in the Mtb genome cannot be called by either method. See also Figure S1 and Table S1.

Figure S1

Genome-scale CRISPRi fitness profiling in Mtb H37Rv, related to Figures 1 and 2, STAR Methods and Table S1 (A) Histogram depicting the number of sgRNAs per gene in the Mtb CRISPRi library (RLC12; Addgene #163954). (B) Next generation sequencing quality-control metrics for the Mtb CRISPRi library. The “Plasmid” column depicts metrics for the RLC12 plasmid library following cloning and isolation from E. coli. The “H37Rv Mtb” column depicts library metrics following transformation and expansion in Mtb H37Rv. Skew ratio represents the ratio between top and bottom 10% of sgRNA counts. (C-F) Correlation heatmap of the triplicate screens depicted in Figure 1A. Panel (C) depicts the correlation between non-targeting sgRNAs in the –ATc cultures; panel (D) depicts the correlation between non-targeting sgRNAs in the +ATc cultures; panel (E) depicts the correlation between TnSeq essential gene targeting sgRNAs in the –ATc cultures; panel (F) depicts the correlation between TnSeq essential gene targeting sgRNAs in the +ATc cultures. G, generation. (G) Boxen plots comparing time-dependent changes in sgRNA L2FC values (mean ± quantiles) comparing –ATc to Input (i.e., generation 0). sgRNAs are grouped according to whether they target genes defined as Essential by TnSeq (n = 63,867 sgRNAs) or Non-Targeting sgRNAs (n = 1,658). ns, not significant. (H) Density plot to detect potential new “bad-seed” sequences. The plot shows the L2FC (+ATc/–ATc at generation 24.3) of all sgRNAs targeting non-essential genes (dashed line), and sgRNAs targeting non-essential genes that contain the indicated sgRNA seed sequences (defined as the five PAM-proximal nucleotides of the sgRNA targeting sequence) displaying the strongest depletion from the library. See STAR Methods “Estimate of the “bad-seed” effect” for more detail. (I) Violin plot showing the behavior of sgRNAs containing the strongest “bad-seed” sequences identified for SpydCas9 (Cui et al., 2018). Only sgRNAs targeting a CRISPRi non-essential gene were analyzed. sgRNAs with a PAM-proximal ‘ACCCA’ sequence (n = 24) show some evidence for target-independent depletion (i.e., "bad-seed" behavior). Dot and error bars represent mean and SD. ^∗p = 0.021; ns, not significant. (J) Heatmap showing the behavior of mismatched sgRNAs in the competitive fitness experiment depicted in Figure 1A. ΔL2FC represents the difference in depletion between essential gene-targeting sgRNAs with perfectly matching targeting sequences and the corresponding mismatched sgRNAs. Mismatched sgRNAs contain mismatches between the sgRNA targeting sequence and the gene target at the indicated position (x axis; 22 is the sgRNA nucleotide furthest from the PAM). Mismatched sgRNAs were not designed but were the result of errors during library synthesis or cloning. (K) Frequency of ATc-resistant colonies that occur after transformation of four unique sgRNAs targeting the essential genes gyrB (ms0005), dnaE1 (ms3178), mmpL3 (ms0250), and pptT (ms2648) in Msmeg. Dots represent transformations performed in biological duplicate; error bars indicate median ± 95% CI. CFU, colony forming unit; NT, non-targeting. (L) Table summarizing the mutations observed in the CRISPRi plasmid in independent ATc-resistant colonies. All but two isolates show unique deletions, duplications, or an inversion (all generically marked as Δ to indicate lack of CRISPRi functionality) within the sgRNA, Cas9, or both. WT, wild-type; TetR, Tet repressor protein; oriE, E. coli origin of replication. (M) Line plot showing all sgRNAs targeting dnaA (rv0001) in the Mtb H37Rv CRISPRi fitness experiment. “Flatliner” sgRNAs of presumed CRISPRi-resistant subpopulations are indicated in green. See STAR Methods for details. (N and O) Distribution of sgRNA depletion slopes (β_e) for sgRNAs targeting essential genes (n = 63,867 sgRNAs) stratified by targeted PAM sequence (N) or sgRNA targeting sequence length (O). Black dots and lines show the median and 25%–75% percentiles. Dot and error bars represent mean and SD. NT, non-targeting.

Figure 2

Identification of features that dictate sgRNA strength (A) Two-line model fits for three different sgRNAs targeting mmpL3. (B) Bar plot showing the regression coefficients (mean ± SEM) for each sgRNA feature identified by the linear model, colored by feature type. All features were represented by more than 500 sgRNAs except for the 20%–30% GC (n = 18) and 90%–100% GC (n = 458) bins. (C) Comparison of measured versus linear model predicted CRISPRi activity (mean ± SEM) of 29 sgRNAs against a Renilla luciferase target in Msmeg; sgRNAs are color coded from blue (strength = 0) to red (strength = 1). The green dot indicates a control non-targeting sgRNA. RLU, relative light unit. (D) Line plot showing the behavior of sgRNAs targeting the essential gene mmpL3 and non-essential gene clgR. sgRNAs are color coded by predicted strengths as in (C). Circles represent our sequencing limit of detection. Triangles represent the point of observation of rare CRISPRi-resistant subpopulations, beyond which sgRNA L2FC values are not plotted (see STAR Methods for details). See also Figure S1 and Data S1.

Figure 3

A quantitative framework to predict gene vulnerability to transcriptional silencing (A) Description of the logistic curve parameters used to model gene level vulnerability. The x axis depicts the linear model predicted strengths of gene-targeting sgRNAs. The y axis depicts the fitness cost of individual sgRNAs estimated by the two-line model from Figure 2A. See details in STAR Methods. (B) Logistic curve fit to all sgRNAs (dots) targeting mmpL3. The black line represents the mean logistic curve and range (gray) from 5,000 parameter samples. Mean parameter estimates and their 95% highest density interval (HDI) are indicated. (C) Logistic regression fits for four example genes of differing vulnerability along with their corresponding VI. Lines represent fits generated by the sampling procedure with the dark line representing the mean fit. (D) Circos plot showing all targeted Mtb H37Rv genes (dots) with their VI. Genes in the upper quartile of vulnerability are depicted as red dots (filled red, confident VI; unfilled red, low-confidence VI). Genes encoding the targets of first line TB therapy (rpoB, inhA, and embAB) are highlighted by blue dots. The outer ring represents the gene-level L2FC value at 28.8 generations. The inner purple lines represent decreasing VI values, with the most vulnerable genes located closest to the center of the circle. The phthiocerol dimycocerosates (PDIM)/phenolic glycolipid (PGL) locus (gray) contains no vulnerable genes. See also Figure S2 and Data S2.

Figure S2

Individual vulnerability model parameters are gene specific, and vulnerability is not correlated with gene expression levels, related to Figure 3 and Data S1 (A) Histogram showing the per-gene Spearman correlation between the rate of depletion $\left({\beta }_e\right)$ estimated from the Bayesian vulnerability model and the predicted strength for targeting sgRNAs. All CRISPRi essential genes with confident vulnerability calls in H37Rv (n = 552) are included in this analysis. (B) Histogram of the vulnerability indices estimated from 5,000 parameter samples for mmpL3 (rv0206c), embB (rv3795), rv2477c, and moeB1 (rv3206c). The vulnerability index 95% credible regions are depicted by dashed lines. (C) Histogram showing the potential influence of the CRISPRi polar effect on vulnerability. The difference in vulnerability index between any downstream gene and its respective upstream gene in the operon is depicted (VI downstream gene – VI upstream gene; n = 657 comparisons). Dashed line depicts the mean difference in VI (mean, 1.658). (D) Histogram of vulnerability indices for genes predicted to be essential by CRISPRi and with confident vulnerability calls, highlighting genes predicted to have an essential domain according to TnSeq (DeJesus et al., 2017). (E) Violin plot depicting the vulnerability index for different groups of genes: all CRISPRi essential genes with confident vulnerability calls (All Ess; n = 552), genes predicted to have an essential domain (Domain Ess; n = 26), genes without an essential domain (Not Domain Ess; n = 526), and genes in the top (n = 138) and bottom (n = 138) quartiles of vulnerability index. Dot and error bars represent mean and SD. Significance (p-value) is calculated using a two-sided t-test. (F-H) Scatterplot of gene vulnerability ratios and/or individual gene parameter estimates. Only confident vulnerability index estimates are shown (see main text for details). (F) depicts the relationship between $\gamma$ and $M$; (G) depicts the relationship between ${\beta }_{max}$ and $M$; (H) depicts the relationship between ${\beta }_{max}$ and $\gamma$. (I and J) Scatterplot showing the relationship between gene mRNA levels as quantified by RNaseq (I) or protein levels as quantified by mass spectrometry (J) (Schubert et al., 2015) and gene vulnerability. Only confident vulnerability index estimates are shown (see STAR Methods for details).

Figure S3

Genome-scale CRISPRi in Msmeg, related to Figure 4, Table S1, and STAR Methods (A) Histogram depicting the number of sgRNAs per gene in the Msmeg CRISPRi library (RLC11; Addgene #163955). The library targets 6,642 of the 6,679 annotated Msmeg genes. (B) Next generation sequencing quality-control metrics for the Msmeg CRISPRi library. The “Plasmid” column depicts metrics for the RLC11 plasmid library following cloning and isolation from E. coli. The “Msmeg” column depicts library metrics following transformation and expansion in Msmeg. Skew ratio represents the ratio between top and bottom 10% of sgRNA counts. (C and D) Correlation heatmap of the triplicate screens performed in Msmeg depicting TnSeq essential gene (Dragset et al., 2019) targeting sgRNAs in the –ATc (C) and +ATc (D) cultures. (E) Boxen plots comparing time-dependent changes in sgRNA L2FC values targeting genes defined as Essential (n = 27,702 sgRNAs) and Non-Essential (n = 120,429) by TnSeq (Dragset et al., 2019). Mean L2FC (solid line) and quantiles beyond the 25^th and 75^th percentiles are shown (boxes). Also depicted are control Non-Targeting sgRNAs (n = 7,421). (F) Hierarchical clustering of gene level depletion from the Msmeg CRISPRi fitness screen. Each row represents a single targeted Msmeg gene. (G) Bar chart showing the overlap between gene calls by TnSeq (Dragset et al., 2019) and CRISPRi. 73% of TnSeq essential calls (291 of 401) are shared with CRISPRi. (H) Histogram showing the per-gene Spearman correlation between the rate of depletion $\left({\beta }_e\right)$ estimated from the Bayesian vulnerability model and the predicted strength for targeting sgRNAs. All CRISPRi essential genes with confident vulnerability calls in Msmeg are included in this analysis. (I) Growth kinetics of the hypomorphic sgRNAs (mean ± SD) shown in Figure 4C. The linear model predicted sgRNA strengths are listed in parentheses next to each gene name. All strains were grown for 15 generations in the presence or absence of ATc and then used to seed cultures for the time-course experiment shown here. Growth for 15 generations ± ATc ensures all strains have reached steady-state growth in response to CRISPRi target gene knockdown. NT, non-targeting. (J) Quantification of target gene mRNA levels by qRT-PCR (biological triplicates; mean ± SEM) of the hypomorphic strains depicted in Figure 4C. (K) Effect of titrating the ATc concentration (range 0-500 ng/mL) on growth (mean ± SD) of the indicated strains from Figure 4C. These strains encode either a non-targeting (NT) sgRNA or a strong sgRNA (predicted strength range, 0.94 – 1.00) against the indicated target. Strains are color coded by vulnerability as in Figure 4D.

Figure 4

Vulnerability predictions correlate with the magnitude of target knockdown needed to reduce bacterial fitness (A) Scatterplot of the linear model coefficients (as in Figure 2B) for Mtb H37Rv (x axis) and Msmeg (y axis). (B) Mean logistic regression fits for the indicated Msmeg genes of varying vulnerability. (C) Phenotypic consequences of hypomorphic (hypo) and strong knockdown of the genes depicted in (B). Predicted sgRNA strengths (P.S.) are listed next to each sgRNA and are color coded according to the scale in Figure 2D. The percent increase in strain doubling time (Dt) of each hypo sgRNA compared with a non-targeting control (95% confidence interval [CI]) was quantified at steady-state growth (Figure S3I). nm, not measured. (D) Quantification of target gene protein levels (mean ± SD) by label-free mass spectrometry (+ATc) of the 6 hypo strains depicted in (C). qRT-PCR quantification of target gene mRNA levels for the same strains is depicted in Figure S3J. See also Figure S3, Table S1, and Data S2.

Figure 5

Pathway analysis identifies differentially vulnerable processes in mycobacteria (A) Heatmap of fitness cost (scaled ${\beta }_e$) as a function of increasing sgRNA strength. Each row represents a single Mtb gene for which a high-confidence VI is available. (B) Table depicting evolutionary conservation between Mtb and eight other bacterial species. For the most vulnerable (VUL; n = 138) and invulnerable (INV; n = 138) H37Rv Mtb genes, the frequency with which a homolog was identified (“genes with homolog”) and the average amino acid similarity (“average similarity of homologs”; % ± SEM) are reported. For the four bacterial species for which genome-wide essentiality calls are available, conservation of essentiality (%) is also listed. M. smeg, M. smegmatis; M. abs, M. abscessus; C. glut, C. glutamicum; B. sub, B. subtilis. ^∗∗∗∗p < 0.0001. ns, not significant. (C) Bubble plot of the enriched (p < 0.05) PATRIC subclasses for the top quartile VUL and bottom quartile INV Mtb and Msmeg (Msm) genes. Conserved subclass enrichment is depicted in bold type. The star represents subclasses where some or all of the corresponding Msmeg homologs are non-essential (Figure S4C), which, for the purposes of this analysis, were considered INV. (D) Logistic regression curves of the indicated Mtb gene groups. Each colored line represents a single gene. The solid black line represents the locally estimated scatterplot smoothing (LOESS) fit of the individual mean logistic regressions. (E) Detailed view of the different vulnerabilities of Mtb genes involved in DNA replication. Genes are color coded by their VI. Darker shades of purple indicate higher vulnerability. The density scale represents the fraction of CRISPRi essential genes with confident VI calls. Figure adapted from (Yao and O’Donnell, 2010). ^∗, low-confidence call. See also Figures S4 and S5 and Data S2.

Figure S4

Evolutionary conservation of vulnerability, related to Figure 5 and Data S2 (A) Violin plot depicting the gene level dN/dS ratios (ω) estimated by GenomegaMap (Wilson and CRyPTIC Consortium, 2020) for five groups of genes: all analyzed Mtb genes (All; n = 3,979), TnSeq non-essential genes (Non-Ess; n = 3,271), TnSeq essential genes (Ess; n = 624), and genes in the top (Vulnerable, n = 138) and bottom (Invulnerable; n = 138) quartiles of vulnerability index. Dot and error bars represent mean and SD. Significance (p-value) is calculated with a two-sided t-test. (B) Sankey plot showing vulnerability conservation between Msmeg and Mtb. 92% of (78 of 85) vulnerable Msmeg genes (VUL; upper quartile) have an Mtb homolog that also ranks in the upper or middle quartile of vulnerability. 82% (70 of 85) of the invulnerable Msmeg genes (INV; lower quartile) have an Mtb homolog that ranks in the lower or middle quartile of vulnerability. (C) Logistic regression curves of the indicated Mtb H37Rv and Msmeg gene groups (PATRIC subclasses indicated above the species name) starred in Figure 5C. Each colored line represents the mean logistic regression curve for a single gene. The solid black line represents the LOESS fit of all logistic regressions. Note that several invulnerable Mtb genes are non-essential in Msmeg. (D) Logistic regression fits and summary LOESS fit for the indicated gene groups (PATRIC subclass) that synthesize the three main mycobacterial cell envelope components. The dashed line is a reference to the solid line of the mycolic acid biosynthesis genes. (E) Vulnerability estimates for the Mtb coenzyme A biosynthetic pathway. Genes are color coded as in Figure 5E. Shown are the mean logistic regression curves for each gene. The solid black line represents the LOESS fit for the indicated gene group. (F) Logistic regression fits for the drug targets rpoB (rv0667) and def (rv0429c). Lines represent fits generated by the sampling procedure with the dark line representing the mean fit. (G) Expression-fitness relationships for genes involved in the cytoplasmic steps of peptidoglycan synthesis between E.coli, B. subtilis (adapted from Hawkins et al., 2020), and Mtb H37Rv, Mtb HN878 and Msmeg. Each colored line represents the expression-fitness relationship for a single gene in the indicated group; the solid black line represents the LOESS fit for the indicated gene group. Only genes with confident vulnerability calls are shown for Mtb and Msmeg. (H) Comparison of the expression-fitness relationships of nadD and nadE between E.coli, B. subtilis (adapted from Hawkins et al., 2020), and Mtb H37Rv, Mtb HN878 and Msmeg.

Figure S5

tRNA synthetases are choke points in Mtb translation, related to Figure 5 and Data S2 Vulnerability estimates for Mtb H37Rv amino acid metabolic genes and tRNA synthetases. Only genes that are CRISPRi essential and have a vulnerability call with high confidence are shown. Genes are color coded as in Figure 5E. The density scale in the figure legend represents the fraction of CRISPRi essential genes with certain vulnerability calls.

Figure S6

Genome-scale CRISPRi in Mtb HN878, related to Figure 6 and STAR Methods (A) Read depth plot of Mtb HN878 whole genome sequencing mapped to the H37Rv genome (GenBank: NC_018143). The 248 SNPs affecting 664 sgRNAs of our CRISPRi library are indicated in red. Significant decreases and increases in read depth mark a genomic deletion and duplication, respectively, in our HN878 clone and are highlighted in gray. (B) Next generation sequencing quality-control metrics for the Mtb HN878 CRISPRi library. The “Plasmid” column depicts metrics for the RLC12 plasmid library following cloning and isolation from E. coli. The “Mtb HN878” column depicts library metrics following transformation and expansion in Mtb HN878. Skew ratio represents the ratio between top and bottom 10% of sgRNA counts. (C) Correlation heatmap of the triplicate screens performed in Mtb HN878 depicting TnSeq essential gene (DeJesus et al., 2017) targeting sgRNAs in the –ATc cultures. (D) Correlation heatmap of the triplicate screens performed in Mtb HN878 depicting TnSeq essential gene (DeJesus et al., 2017) targeting sgRNAs in the +ATc cultures. (E) Scatterplot of the linear model coefficients (as in Figure 2B) for Mtb H37Rv (x axis) and Mtb HN878 (y axis). (F) Histogram showing the per-gene Spearman correlation between the rate of depletion $\left({\beta }_e\right)$ estimated from the Bayesian vulnerability model and the predicted strength for targeting sgRNAs. All CRISPRi essential genes with confident vulnerability calls in HN878 are included in this analysis. (G) Liquid growth assay (mean ± SD) using the sgRNAs targeting two differentially essential genes used in Figure 6B. NT, non-targeting. (H) Quantification of target gene mRNA levels by qRT-PCR (n = 6 technical replicates; mean ± SEM) following CRISPRi silencing of rv2017 and rv2228c in H37Rv and HN878. Gene expression levels were normalized to the non-targeting control for each strain. (I) Quantification of target gene mRNA levels by qRT-PCR (technical triplicates of biological duplicates) of cydABCD and qcrCAB in H37Rv and HN878. For each gene, HN878 expression levels were compared to H37Rv (control).

Figure 6

Differential VI predicts strain-specific susceptibility to antibacterial agents (A) Bar chart showing the overlap between CRISPRi gene essentiality calls in H37Rv and HN878. (B) CRISPRi knockdown of two genes predicted to be essential in H37Rv and non-essential in HN878. NT, non-targeting. (C) Correlation between VI in H37Rv and HN878 for all genes (black) and CRISPRi essential genes for which high-confidence VI calls are available (blue). (D) Histogram showing the normalized differential VI between HN878 and H37Rv for genes with a high-confidence call in both strains. Quartiles are delineated by a dotted line. (E) Logistic regression fits for accD6 in H37Rv (black) and HN878 (turquoise). Lines represent fits generated by the sampling procedure with the dark line representing the mean fit. (F) Phenotypic consequences of accD6 knockdown. The optical density 600 (OD₆₀₀) L2FC (+ATc/–ATc; mean ± SD) was calculated for three accD6 sgRNAs (1–3) and a non-targeting control sgRNA in H37Rv and HN878. Strains were pre-treated with ATc for 3 days prior to starting the depicted time course. (G) Bubble plot of the enriched (p < 0.05) PATRIC subclasses for genes more VUL in HN878 versus H37Rv. The star represents a subclass where some or all of the corresponding H37Rv homologs are non-essential, which, for the purposes of this analysis, were considered INV. (H–L) Effect of rifampicin (H), ethambutol (I), isoniazid (J), Q203 (K), and ND-10885 (L) on growth (mean ± SD) of H37Rv and HN878. (M) Gene-level L2FC measurements for cydABCD and inhA from the H37Rv and HN878 CRISPRi screens at ~29 generations. (N and O) Effect of novobiocin (N) and SPR719 (O) on growth (mean ± SD) of H37Rv and HN878. See also Figure S6 and Data S2.