An Integrated Genome-wide CRISPRa Approach to Functionalize lncRNAs in Drug Resistance

Abstract

Resistance to chemotherapy plays a significant role in cancer mortality. To identify genetic units affecting sensitivity to cytarabine, the mainstay of treatment for acute myeloid leukemia (AML), we developed a comprehensive and integrated genome-wide platform based on a dual protein-coding and non-coding integrated CRISPRa screening (DICaS). Putative resistance genes were initially identified using pharmacogenetic data from 760 human pan-cancer cell lines. Subsequently, genome scale functional characterization of both coding and long non-coding RNA (lncRNA) genes by CRISPR activation was performed. For lncRNA functional assessment, we developed a CRISPR activation of lncRNA (CaLR) strategy, targeting 14,701 lncRNA genes. Computational and functional analysis identified novel cell-cycle, survival/apoptosis, and cancer signaling genes. Furthermore, transcriptional activation of the GAS6-AS2 lncRNA, identified in our analysis, leads to hyperactivation of the GAS6/TAM pathway, a resistance mechanism in multiple cancers including AML. Thus, DICaS represents a novel and powerful approach to identify integrated coding and non-coding pathways of therapeutic relevance.

Keywords: AML; AXL/GAS6; CRISPR; CRISPRa; TEM; cancer; cytarabine; drug-resistance; leukemia; lncRNA.

Figure 1. Identification of Protein-Coding and Noncoding Gene Biomarkers Correlated with Differential Ara-C Response

(A) Distribution of Ara-C drug sensitivities across 760 pan-cancer cell lines profiled by both CCLE and CTD² studies, quantified by their Z-scaled area under the dose response curve values after regressing out lineage-specific effects. See also Table S1. (B) Distribution of Z-scaled drug resistance-gene expression Pearson correlation values of all analyzed genes. Representative protein-coding and non-coding gene symbols enriched beyond a Z-score threshold of ± 1.16 are demarcated. See also Table S1. (C) Summary of gene set enrichment analysis (GSEA) of protein-coding genes ranked by drug resistance-gene expression correlation values using annotated KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways. See also Table S3. (D) Representative KEGG pathways from GSEA of protein-coding genes ranked by drug sensitivity-gene expression correlation values as shown in Figures 1B–1C. See also Table S3. (E) Pearson correlation distributions of gene pair expression levels in the cancer cell line panel across 997 sense-antisense cognate gene pairs and 5,000 random protein coding- lncRNA gene pairs. Wilcoxon rank-sum test: p < 2.2e-16. (F) Relationship of drug sensitivity-gene expression correlation values between protein coding-lncRNA gene pairs across 997 sense-antisense cognate gene pairs (left panel: Pearson’s R = 0.552, p < 2.2e-16) and 5,000 random gene pairs (right panel: Pearson’s R = 0.021, p = 0.1338).

Figure 2. CRISPRa Functional Screening of Coding Genes Modulating Ara-C Response

(A) Distribution of Ara-C IC50 values across a panel of AML cell lines. (B) Effect of BCL2 overexpression (Blue) or DCK knockdown on sensitivity to Ara-C in MOLM14 cells. Data are represented as mean ± SD, n = 3. (C) Schematic of CRISPRa pooled screening for the identification of genes whose activation modulate sensitivity to Ara-C in MOLM14 cells. (D) Volcano plot summarizing the global changes in sgRNA representation of protein-coding genes before and after 14 days of treatment with Ara-C. A subset of genes validated herein (red text) or previously annotated (black text) to modulate Ara-C sensitivity are labeled. A false discovery rate threshold of 0.339 was determined by receiver operating characteristic analysis (Figure S3F). Red - enrichment in the CRISPRa screening; blue depletion in the CRISPRa screening; open black circles - genes previously associated with differential Ara-C sensitivity and above the significance threshold; filled black points genes validated herein. See also Figure S3C-F, S3H, and Table S4. (E) Summary of gene set enrichment analysis (GSEA) of protein-coding genes ranked by CRISPRa screening using annotated KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways. See Table S3. (F) Disease-free survival association with expression levels of ZBP1, MUL1, and PI4K2A, genes enriched in both protein-coding CRISPRa screening and drug sensitivity-gene expression correlation analyses among patients treated with Ara-C therapy within the TCGA-LAML patient cohort. ZBP1: VST expression level cutoff = 6.13 (low, n = 42; high, n = 79), log-rank test: p-value = 0.0074. MUL1: VST expression level cutoff = 9.64 (low, n = 108; high, n = 13), log-rank test: p-value = 0.0033. PI4K2A: VST expression level cutoff = 7.23 (low, 36; high, n = 85), log-rank test: p-value = 0.038. (G) Ara-C efficacy measurements in MOLM14 cells expressing sgRNAs targeting ZBP1, MUL1, or PI4K2A based on normalized MTS reads following 48 hours of treatment. Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p<0.001 (H) Modulation of apoptotic response upon stable expression of sgRNAs targeting ZBP1, MUL1, or PI4K2A in MOLM14 cells. The percentage of apoptosis is determined by annexin V and propidium iodide (PI) staining of cells treated with 0.25 pM Ara-C for 72 hours. Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001 (I) Proliferation of unchallenged MOLM14 cells expressing sgRNAs targeting ZBP1, MUL1, or PI4K2A. Proliferation is quantified over four days (D1-D4). Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001

Figure 3. CRISPRa Functional Screening of Noncoding Genes Modulating Ara-C Response

(A) Left panel: summary of the CaLR library design specifications, including lncRNA gene numbers, transcriptional start sites (TSS), and total sgRNA numbers. Right panel: relationships between coding genes and lncRNA genes for corresponding lncRNA classifications. See also Table S5. (B) Volcano plot summarizing the global changes in sgRNA representation of noncoding genes before and after 14 days of treatment with Ara-C. A subset of genes either validated herein to modulate Ara-C sensitivity (red text) or previously annotated in various cancer-related pathways (black text) are labeled. A false discovery rate threshold of 3.51e-5 was determined by analysis of nontargeting sgRNA negative controls at the transcript level (Figure S4H). Red points - enrichment in the CRISPRa screening; blue points - depletion in the CRISPRa screening; filled black points - genes validated herein. See also Figure S4E-I and Table S6. (C) Percentages of significantly enriched or depleted protein-coding or noncoding genes from CRISPRa screens detected in the TCGA-LAML patient samples. Chi-squared test: ***, p = 6.92e-3, (D) Gene expression level distributions of significantly enriched or depleted protein-coding or noncoding genes from CRISPRa screens detected in the TCGA-LAML patient samples. Wilcoxon rank-sum test: ***, p = 5.4e-7. (E) Guilt-by-association pathway annotation of enriched genes identified in the CaLR screen. KEGG pathway gene sets were used for this analysis.

Figure 4. Validation of CaLR Screening Results

(A) Fold change (FC) of MOLM14 cell viability treated with 0.25 pM Ara-C for 48 hours. Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01,***, p < 0.001. (B) Fold change (FC) of expression levels of targeted lncRNAs upon overexpression of enriched sgRNAs versus endogenous levels. Data are represented as mean ± SD, n = 3. (C) Ara-C efficacy measurements in MOLM14 cells expressing sgRNAs targeting indicating genes based on normalized MTS reads following 48 hours of treatment with the indicated concentrations of Ara-C. Data are represented as mean ± SD, n = 3, Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001 (D) Proliferation of unchallenged MOLM14 cells expressing sgRNAs targeting indicating genes. Proliferation is quantified over four days (D1-D4). Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. (E) Left panel: modulation of apoptotic response upon stable expression of sgRNAs targeting a panel of significantly enriched sgRNAs as determined through CaLR screening in MOLM14 cells. The percentage of apoptosis is determined by annexin V and propidium iodide (PI) staining of MOLM14 cells stably infected with individual sgRNAs and treated with 0.25 pM Ara-C for 72 hours. Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. Right panel: representative flow cytometry plots of annexin V/PI staining intensities corresponding to two sgRNAs promoting survival versus nontargeting control. (F) Immunofluroscence images (left panel) for DAPI and phospho-YH2A.X staining in MOLM14 cells stably infected with sgRNAs targeting the lncRNA genes shown, and treated with 25 pM Ara-C for 24 hours. Staining is quantified in the right panel. Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. (G) Disease-free survival association with expression levels of GAS6-AS2 and AC008073.2, genes enriched in both noncoding CRISPRa screening and drug resistance-gene expression correlation analyses among patients treated with Ara-C therapy within the TCGA-LAML patient cohort. GAS6-AS2: VST expression level cutoff = 3.38 (low, n = 92; high, n = 29), log-rank test: p-value = 0.035. AC008073.2: VST expression level cutoff = 4.39 (low, n = 93; high, n = 28), log-rank test: p-value = 0.0026.

Figure 5. GAS6-AS2 Promotes Drug Resistance In Vitro and In Vivo

(A) Integration of drug resistance-gene expression correlative analysis and forward genetic screenings identifies seven sense-antisense gene pairs which pass all significance thresholds, a higher number than expected by chance alone (Chi-squared test: p = 9.85e-7). (B) Fold change (FC) of MOLM14 cell viability treated with 0.25 pM Ara-C for 48 hours. Cells expressing individual sgRNAs targeting GAS6-AS2. Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. (C) Pearson correlation between cell viability versus GAS6-AS2 expression level for each of the 8 sgRNAs targeting GAS6-AS2. (D) Ara-C efficacy measurements in MOLM14 cells expressing sgRNAs #1 and #3 targeting GAS6-AS2 based on normalized MTS reads following 48 hours of treatment. Data are represented as mean ± SD, n = 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. (E) Left panel: representative flow cytometry data of MOLM14 cells expressing either control or GAS6-AS2-targeting sgRNAs, treated with 25 pM Ara-C for 24 hours and labeled with viability (propidium iodide (PI)) and apoptotic (annexin V) markers. Right panel percentage of apoptosis determined from quantification of staining results. Data are represented as mean ± SD, n > 3, Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. (F) Competition assay between populations of MOLM14 control-Blue and MOLM14 GAS6- AS2-Red following 25 pM Ara-C treatment. Left panels: representative flow cytometry plots. Right panel: ratios between red and blue cells over time. Data are represented as mean ± SD, n > 3. Welch two sample t-test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. (G) Schematic of an orthotopic xenograft competition assay between control (blue) and GAS6-AS2 (Red) MOLM14 cells with Ara-C treatment. (H) Ratios of control (blue) versus GAS6-AS2 (Red) MOLM14 cells from bone marrow of mice treated and analyzed at day 17 as outlined in Figure 5G. (I) Representative flow cytometry results of cells harvested from mouse bone marrow 17 days following transplantation and treatment with vehicle or Ara-C for 5 days.

Figure 6. GAS6-AS2 Activates GAS6/TAM Signaling

(A) Pearson correlation between GAS6-AS2 and GAS6 expression levels following GAS6- AS2 activation. Data are represented as mean of triplicate measurements. (B) Pearson correlation between GAS6-AS2 and GAS6 expression levels across the 760 cancer cell lines analyzed (Figure 1A-B). (C) Pearson correlation between GAS6-AS2 and GAS6 expression levels in AML patient samples. (D) Western blot analysis of differential GAS6/TAM signaling activation in response to individual control or GAS6-AS2 sgRNA overexpression. (E) Pearson correlation between GAS6-AS2 and AXL expression levels in AML patient samples. (F) Pearson correlation between GAS6-AS2 and AXL expression levels across the 760 cancer cell lines analyzed (Figure 1A-B). (G) Expression levels of GAS6-AS2, GAS6, and AXL in MOLM14 and K562 cell lines. (H) Ara-C efficacy measurements in MOLM14 and K562 cell lines, based on normalized MTS reads following 48 hours of treatment with the indicated concentrations of Ara-C. Data are represented as mean ± SD, n = 3.

Figure 7. GAS6-AS2 Demonstrates Trans-Regulation of AXL

(A) Fold change (FC) of GAS6-AS2, GAS6, and AXL in response to GAS6-AS2 knockdown via ASO in K562 cells. Data are represented as mean ± SD, n = 3. Welch two sample t- test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. (B) Modulation of Ara-C response upon GAS6-AS2 knockdown via ASO in K562 cells. Data are represented as mean ± SD, n = 3, Welch two sample t-test: *, p < 0.05. **, p < 0.01,***, p < 0.001. (C) Methylation of CpG islands in the HEK293T AXL promoter following modulation of GAS6-AS2 expression. n = 12, Chi-square test: *, p < 0.05. **, p < 0.01, ***, p < 0.001. (D) Gene ontology analysis of coding genes clustered with GAS6-AS2 as determined by k- means clustering (cluster #6 in Figure S7D). (E) Drug sensitivity-gene expression Pearson correlation values of DNA methyltransferases. Genes enriched beyond a Z-score threshold of ± 1.16 are colored in red. See also Figure 1B. (F) Distribution of FPKM-normalized transcript abundances associated with DNMT1 versus IgG. (G) Model summarizing the mechanism by which GAS6-AS2 regulates GAS6/TAM signaling.