1. Cancer Biology
  2. Chromosomes and Gene Expression

Single-PanIN-seq unveils that ARID1A deficiency promotes pancreatic tumorigenesis by attenuating KRAS-induced senescence

  1. Shou Liu
  2. Wenjian Cao
  3. Yichi Niu
  4. Jiayi Luo
  5. Yanhua Zhao
  6. Zhiying Hu
  7. Chenghang Zong  Is a corresponding author
  1. Department of Molecular and Human Genetics, Baylor College of Medicine, United States
  2. Genetics and Genomics Graduate Program, Baylor College of Medicine, United States
  3. Cancer and Cell, Biology Graduate Program, Baylor College of Medicine, United States
  4. Dan L Duncan Comprehensive Cancer Center, Baylor College of Medicine, United States
  5. McNair Medical Institute, Baylor College of Medicine, United States
https://doi.org/10.7554/eLife.64204
Share this article
Cite this article
  1. Shou Liu
  2. Wenjian Cao
  3. Yichi Niu
  4. Jiayi Luo
  5. Yanhua Zhao
  6. Zhiying Hu
  7. Chenghang Zong
(2021)
Single-PanIN-seq unveils that ARID1A deficiency promotes pancreatic tumorigenesis by attenuating KRAS-induced senescence
eLife 10:e64204.
https://doi.org/10.7554/eLife.64204
  2. Download BibTeX
  3. Download .RIS

Abstract

ARID1A is one of the most frequently mutated epigenetic regulators in a wide spectrum of cancers. Recent studies have shown that ARID1A deficiency induces global changes in the epigenetic landscape of enhancers and promoters. These broad and complex effects make it challenging to identify the driving mechanisms of ARID1A deficiency in promoting cancer progression. Here, we identified the anti-senescence effect of Arid1a deficiency in the progression of pancreatic intraepithelial neoplasia (PanIN) by profiling the transcriptome of individual PanINs in a mouse model. In a human cell line model, we found that ARID1A deficiency upregulates the expression of aldehyde dehydrogenase 1 family member A1 (ALDH1A1), which plays an essential role in attenuating the senescence induced by oncogenic KRAS through scavenging reactive oxygen species. As a subunit of the SWI/SNF chromatin remodeling complex, our ATAC sequencing data showed that ARID1A deficiency increases the accessibility of the enhancer region of ALDH1A1. This study provides the first evidence that ARID1A deficiency promotes pancreatic tumorigenesis by attenuating KRAS-induced senescence through the upregulation of ALDH1A1 expression.

Introduction

Pancreatic ductal adenocarcinoma (PDAC) was the third leading cause of cancer-related death in the United States in 2018 and is projected to become the second leading cause of cancer-related death by 2030 (Rahib et al., 2014). Besides the clear driver mutations of PDAC, which include CDKN2A, TP53, SMAD4, and KRAS, there are a large number of genes with low-frequency mutations (Waddell et al., 2015; Bailey et al., 2016; Cancer Genome Atlas Research Network, 2017). The statistical significance of a large number of low-frequency mutations suggests that they likely play functional roles in the promotion of tumor development. From the evolutionary perspective of tumor development, the probability of acquiring one of the mutations in low-frequency-mutation genes is significantly higher than the acquisition of one of the mutations in the very few high-frequency-mutation genes such as TP53, especially at the very early stage when the number of pre-lesion cells is limited. Therefore, it is reasonable to believe that a portion of the low-frequency mutations could be acquired at the early stage and that these mutations could play essential roles in tumorigenesis.

Among the genes with low-frequency mutations in PDAC, the top hit is ARID1A (a subunit of the SWI/SNF chromatin remodeling complex) with an 8% mutation rate (Figure 1—figure supplement 1). In addition to PDAC, ARID1A is also frequently mutated in other cancer types, with 45.2% mutation rate in ovarian cancer, 18.7% in gastric cancer, 18.6% in bladder cancer, 13.7% in hepatocellular cancer, 11.5% in melanoma, 9.4% in colorectal cancer, 8.2% in lung cancer, and 2.5% in breast cancer (Kadoch et al., 2013). Therefore, it is crucial to study the ubiquitous mechanisms for ARID1A-deficiency-facilitated tumorigenesis in various types of cancers.

ARID1A is involved in the regulation of many biological processes of cells, including differentiation, proliferation, and apoptosis (Wu and Roberts, 2013). In pancreatic cancer, recent studies have shown that ARID1A is necessary to maintain terminal differentiation of pancreatic acinar cells, and knockout of ARID1A results in the accelerated formation of acinar-to-ductal metaplasia (ADM) and then pancreatic intraepithelial neoplasia (PanIN) lesions (Zhang, 2018; Livshits et al., 2018; Wang et al., 2019). Although ARID1A depletion can prime acinar cells for early-stage PanIN lesion formation by facilitating shifts in cell identity (Livshits et al., 2018), the underlying mechanisms for the acceleration of ARID1A-deficiency-promoted PanIN progression remain elusive. Besides PanIN lesions, ARID1A deficiency could also facilitate intraductal papillary mucinous neoplasm (IPMN) formation through multiple pathways, including MYC-mediated protein synthesis (Wang et al., 2019) and SOX9/mTOR pathway (Kimura et al., 2018).

To dissect the mechanisms whereby Arid1a knockout drives PanIN progression, we applied the single-cell RNA-seq method to profile the transcriptome of individual early-stage PanIN lesions from Arid1a knockout and wildtype mice. Our results showed that Arid1a knockout could effectively reduce KRAS-induced senescence in PanIN lesions. It is important to point out that cellular senescence has been shown to be the major rate-limiting step in KRAS-driven PanIN progression (Morton et al., 2010; Serrano et al., 1997; Li et al., 2018). Therefore, with significant attenuation of senescence, Arid1a knockout can achieve significant acceleration of PanIN progression. Mechanistically, we found that aldehyde dehydrogenases play an essential role in attenuating senescence by scavenging the reactive oxygen species (ROS) induced by oncogenic KRAS.

Results

Individual PanIN lesion RNA-Seq unveils the potential player contributing to the attenuation of Kras-induced senescence in Arid1a knockout mice

To identify the mediators that contribute to Arid1a-deficiency-promoted PanIN progression, we followed the PanIN progression in conditional Arid1a knockout mice with mutant Kras (Arid1afl/fl;Lox-Stop-Lox-KrasG12D/+;Ptf1aCreERT/+, we abbreviate it as AKC) and mice without Arid1a knockout (Arid1a+/+;Lox-Stop-Lox-KrasG12D/+;Ptf1aCreERT/+, we abbreviate it as KC) (Figure 1—figure supplement 2A). Consistent with the findings from other groups (Zhang, 2018; Livshits et al., 2018; Wang et al., 2019), we observed that Arid1a knockout facilitated the progression of lesions from ADM to PanIN3 (Figure 1—figure supplement 2B,C). At as early as the 2-month time point, the percentages of ADM, PanIN-1, and PanIN-2 were 74, 26, and 0%, respectively, in KC mice versus 53, 47, and 0% in AKC mice. At the 6-month time point, the percentages of ADM, PanIN-1, PanIN-2, and PanIN-3 were 77, 23, 0, and 0%, respectively, in KC mice versus 16, 74, 9, and 0.5% in AKC mice (Figure 1—figure supplement 2C). We observed the decrease in the percentage of ADM and increase in the percentage of high-grade PanIN lesions, indicating the accelerated progression of lesions.

To profile the transcriptome of individual PanIN lesions, we combined laser capture microdissection (LCM) with a highly sensitive single-cell RNA-seq method (MATQ-seq) developed in our lab (Sheng et al., 2017; Figure 1A). In total, we dissected and profiled 20 lesions from two KC mice and 24 lesions from two AKC mice. We only dissected PanIN-1 and PanIN-2 lesions because their duct-like structures can be easily recognized for dissection from frozen sections (Figure 1—figure supplement 3A). We observed that the lesions dissected from KC and AKC mice were clearly separated in the multidimensional scaling (MDS) plot, and 861 differentially expressed genes (DEGs) (Supplementary file 2), including 532 genes upregulated and 329 genes downregulated in AKC, were identified, which suggests the successful transcriptome profiling of individual PanIN lesions by MATQ-seq (Figure 1B).

Figure 1 with 4 supplements see all
Download asset Open asset
Single pancreatic intraepithelial neoplasia (PanIN) lesion RNA-seq unveils the potential player contributing to the attenuation of Kras-induced senescence in Arid1a knockout mice.

(A) Experimental scheme of transcriptome profiling of single PanIN lesions. (B) Multidimensional scaling plot showed a clear separation of the transcriptome profiles of lesions from AKC mice (Arid1afl/fl;LSL-KrasG12D/+;Ptf1aCreERT/+) and KC mice (Arid1a+/+;LSL-KrasG12D/+;Ptf1aCreERT/+). 24 lesions from 2 AKC mice and 20 lesions from 2 KC mice were used for single-lesion RNA sequencing. (C) Hallmark gene sets that are significantly enriched between lesions from AKC and KC mice. (D, E) Enrichment plots of P53_PATHWAY and SENESCENCE_UP from gene set enrichment analysis.

With DEGs observed between AKC and KC lesions, we performed gene set enrichment analysis (GSEA) using Hallmark gene sets (Liberzon et al., 2015) to interrogate the pathways perturbed by Arid1a knockout. We observed that 25 gene sets were downregulated, and surprisingly, only two gene sets were upregulated in AKC lesions (Figure 1C, Figure 1—figure supplement 4, and Supplementary file 1, false discovery rate (FDR) < 0.1). It is worth noting that among the 27 gene sets 2 are specifically associated with Kras activation: KRAS_SIGNALING_UP (the gene set upregulated upon Kras activation) and KRAS_SIGNALING_DN (the gene set downregulated upon Kras activation). The gene set KRAS_SIGNALING_UP was downregulated while the gene set KRAS_SIGNALING_DN was upregulated (Figure 1C and Figure 1—figure supplement 4). This observation suggests that the activities of Kras signaling are partially impaired by Arid1a deficiency.

Furthermore, we observed that the Tp53 signaling pathway was suppressed in AKC lesions (Figure 1D). It has been well established that upregulation of the Tp53-related pathway is closely associated with apoptosis or senescence. ARID1A mutations have also been shown to be mutually exclusive with TP53 mutations in endometrial cancer (Wu et al., 2017). To determine whether Arid1a is involved in the regulation of apoptosis, senescence, or both, we further examined the activity of related pathways in Arid1a KO lesions. Interestingly, we found that the senescence-associated signaling pathway is significantly suppressed in lesions from AKC mice (Figure 1E). In contrast, the pathway activity associated with apoptosis was not significantly changed (Figure 1—figure supplement 3B). These observations led us to hypothesize that Arid1a deficiency could promote PanIN lesion progression through the attenuation of Kras-induced senescence.

Furthermore, senescent cells feature senescence-associated secretory phenotype (SASP), including high levels of inflammatory cytokines and immune modulators. With the attenuation of senescence promoted by Arid1a deficiency, we expected to observe reduced inflammatory response in the GSEA. Indeed, we observed that multiple signaling pathways associated with inflammation, including TNFα signaling, IL6–STAT3 signaling, IL2–STAT5 signaling, IFN-α signaling, and IFN-γ signaling, were significantly suppressed in Arid1a KO lesions (Figure 1C, Figure 1—figure supplement 4, and Supplementary file 1).

In vivo, ex vivo, and in vitro verification of the attenuation of Kras-induced senescence by Arid1a deficiency

To verify the effect of Arid1a deficiency on Kras-induced senescence, we performed senescence-associated beta-galactosidase (SA-β-Gal) staining on lesions from KC and AKC mice. SA-β-Gal-positive lesions were observed in five out of seven (71%) KC mice. In contrast, only one out of six (17%) AKC mice showed SA-β-Gal-positive lesions. Among the mice with SA-β-Gal-positive lesions, the percentage of SA-β-Gal-positive lesions in KC mice was about twice of that in AKC mice (Figure 2A,B). These data confirmed that Arid1a knockout indeed reduced Kras-induced senescence.

Figure 2 with 1 supplement see all
Download asset Open asset
In vivo, ex vivo, and in vitro verification of attenuation of Kras-induced senescence by Arid1a deficiency.

(A) Representative images of senescence-associated beta-galactosidase (SA-β-Gal) staining of frozen pancreatic sections from KC mice and AKC mice. (B) SA-β-Gal-positive lesions were counted at five random fields under the microscope in four KC mice and one AKC mouse, presented as percentages. (C) Representative images of SA-β-Gal staining of ex vivo culture from KC and AKC mice 1 month after administration of tamoxifen. (D) Quantification of the intensity of SA-β-Gal staining at five random fields under the microscope. (E) Colony formation assay of ARID1A knockout cells and wildtype human pancreatic Nestin-expressing (HPNE) cells with KRAS induction by doxycycline (6 µg/ml) for 15 days. (F) Quantification of colony number in panel (E). The colony formation assay was performed twice. Student’s t-test: **p<0.01; ***p<0.001. Scale bars: 200 µm.

To further verify the effects of Arid1a knockout on senescence, we performed an ex vivo culture experiment using acinar cells isolated from AKC and KC mice. SA-β-Gal staining was performed to examine the senescence of acinar cells. As shown in Figure 2C,D, due to flat morphology of senescent cells and the large cell-size variation, we cannot accurately quantify the number of SA-β-Gal-negative cells. Instead of using the percentage of senescence cells, here we quantified the senescence based on the intensity of SA-β-Gal staining of the positive cells. We observed that the intensity of SA-β-Gal staining in acinar cells from AKC mice was significantly less than that from KC mice. Importantly, this result suggests that the effect of Arid1a knockout on senescence is likely an intrinsic cellular response rather than a microenvironmental response such as accelerated clearance of senescent cells by immune cells (Pignolo et al., 2020; Baker et al., 2011).

Based on the results of the ex vivo experiment, we next tested whether we could observe similar effects of ARID1A deficiency on KRAS-induced senescence in in vitro cell lines. We established a human pancreatic Nestin-expressing (HPNE) cell line (an intermediary cell type formed during ADM) with inducible KRASG12D expression (iKRAS-HPNE cells). The induction of KRAS activity was verified by western blot of both KRAS and phosphorylated ERK (Figure 2—figure supplement 1A). Next, we knocked out ARID1A by CRISPR-Cas9 in iKRAS-HPNE cells. Two isogenic HPNE clones were used in the following experiments (Figure 2—figure supplement 1B), and the knockout efficiency of ARID1A was further confirmed by qRT-PCR (Figure 2—figure supplement 1C).

To examine the effect of ARID1A deficiency on KRAS-induced senescence, we first performed SA-β-Gal staining in ARID1A knockout (ARID1A-KO) and wildtype iKRAS-HPNE cells upon KRAS induction. We observed that ARID1A-KO cells showed a slightly lower percentage of SA-β-Gal-positive cells compared with wildtype iKRAS-HPNE cells (Figure 2—figure supplement 1D). The in vitro result of SA-β-Gal staining is less significant than what we observed both in vivo and ex vivo. One possible reason for this discrepancy is that in terms of evaluation of senescence SA-β-Gal staining may not be as effective in immortalized cell lines as in tissue samples (Dimri et al., 1995). Colony formation assay is an alternative method to measure senescence. Here, we performed a colony formation assay to evaluate the ability of cells to escape from KRAS-induced senescence. We observed that knockout of ARID1A drastically increased the number of cells escaping from KRAS-induced senescence (Figure 2E,F).

ARID1A knockout reduces KRAS activities and inflammatory response

With the consistent in vivo, ex vivo, and in vitro observation of the anti-senescence effects of ARID1A deficiency, we next investigated the molecular mechanisms by which ARID1A knockout promotes the escape from KRAS-induced senescence. We first performed RNA-seq on ARID1A-KO (clone #2) and wildtype HPNE cells with or without KRAS induction. We observed a clear separation between ARID1A-KO and wildtype HPNE cells under both conditions in the MDS plot (Figure 3A).

Figure 3 with 4 supplements see all
Download asset Open asset
ARID1A knockout upregulates aldehyde dehydrogenase (ALDH) expression.

(A) Multidimensional scaling plot demonstrated clear separation between the transcriptome profiles of ARID1A-KO human pancreatic Nestin-expressing (HPNE) cells and wildtype cells with or without KRAS induction. RNA sequencing was performed with three biological repeats. (B) Volcano plot of differentially expressed genes between ARID1A knockout cells and wildtype cells with KRAS induction. (C) Venn diagram showing the upregulated genes (upper) and downregulated genes (bottom) that are shared between cells with (gray) or without (blue) KRAS induction. (D) ALDH1A1 mRNA levels quantified by sequencing data are significantly different between ARID1A-KO cells and wildtype cells with (left) or without (right) Kras induction. CPM: count per million reads. (E) Western blot for ALDH1A1 expression in ARID1A-KO cells and wildtype cells with KRAS induction. (F) mRNA level of Aldh3a1 in KC and AKC lesions based on pancreatic intraepithelial neoplasia (PanIN)-seq data. APM: amplicon per million reads. (G) IHC staining against ALDH3A1 in KC and AKC lesions. Scale bars: 200 µm. (H) Comparison of ALDH3A1 levels between KC and AKC lesions based on the intensity of staining in (G). H-score was calculated by counting the number of lesions with different levels of staining intensity at four random fields under the microscope. Student’s t-test: ***p<0.001; ****p<0.0001.

Next, we performed GSEA using the Hallmark gene sets. Interestingly, we found the gene set KRAS_SIGNALING_UP was downregulated in ARID1A-KO cells compared with wildtype cells (Figure 3—figure supplement 1A,B), which is consistent with our data from PanIN lesions. To verify this observation, we examined the activity of ERK, a classical downstream effector of KRAS signaling, in ARID1A-KO and wildtype cells. As shown in Figure 3—figure supplement 2, the phosphorylation of ERK upon KRAS induction in ARID1A-KO cells was significantly reduced in comparison with the wildtype cells. The consistent observation between in vivo and in vitro indicates that ARID1A deficiency partially impairs the activities of KRAS signaling.

Next, we performed an interaction test on the gene expression data to identify the genes that have different responses to KRAS activation depending on ARID1A status (Supplementary file 4). GSEA results showed that six signaling pathways were remarkably repressed in ARID1A-KO HPNE cells, and four of them are involved in inflammatory response (Figure 3—figure supplement 3A,B), which is consistent with what we observed in PanIN lesions from AKC mice (Figure 1C).

Furthermore, we examined the expression of a classic marker of cellular senescence: CDKN1A. We observed that the activation of CDKN1A expression upon KRAS induction was also significantly reduced in the ARID1A-KO cell line (Figure 3—figure supplement 3C,D), which indicates that our HPNE cell line model successfully recapitulated the senescence phenotypes observed in the mouse model.

ARID1A knockout significantly upregulates the expression of aldehyde dehydrogenase (ALDH) family members

To identify the underlying players that promote the attenuation of cellular senescence, we analyzed the DEGs between wildtype and the ARID1A-KO cell line with or without KRAS induction (Figure 3B, Figure 3—figure supplement 4A, and Supplementary file 3). For the first clone with KRAS induction, we identified 125 upregulated genes and 165 downregulated genes between the wildtype and the ARID1A-KO line (Figure 3B). To exclude the genes whose expression changes may be associated with mutant KRAS signaling, we compared the list of DEGs under two conditions: with or without KRAS induction. As shown in Figure 3C, for the upregulated genes, 57 out of 125 genes (46%) are shared between the two conditions and the expected number of random overlap is 1.79. For the downregulated genes, 54 out of 165 genes (33%) are shared and the expected number of random overlap is 3.38. These results indicate that these genes are mainly dependent on ARID1A deficiency.

Among the DEGs between ARID1A-KO and wildtype HPNE cells, ALDH1A1 exhibits the significant change in differential gene expression for both conditions: with or without KRAS induction (Figure 3B and D, Figure 3—figure supplement 4A, and Supplementary file 3). To rule out any potential clonal bias, we also performed RNA-seq on a second clone (clone #11). We observed that ALDH1A1 was also significantly upregulated in the second clone under both conditions (Figure 3—figure supplement 4B–D and Supplementary file 3). The upregulation of ALDH1A1 in ARID1A-KO cells was further verified by both qRT-PCR (Figure 3—figure supplement 4E) and western blot (Figure 3E). Considering that ALDH1A1 has been shown to participate in the clearance of ROS (Raha et al., 2014) and ROS are vital mediators of KRAS-induced senescence (Storz, 2017), we hypothesize that ALDH1A1 is the gene that mediates the effect of ARID1A deficiency on KRAS-induced senescence.

Next, we examined our PanIN-seq data to evaluate the expression of Aldh1a1 and other members of the ALDH family. Interestingly, we observed that Aldh3a1 is significantly upregulated in lesions from AKC mice (Figure 3F and Supplementary file 2), which was further confirmed by immunohistochemistry (IHC) staining (Figure 3G,H). This result suggests that in different species different types of ALDH family proteins can be used to mediate the attenuation of Kras-induced senescence in Arid1a-deficient cells.

ARID1A KO facilitates escape from KRAS-induced senescence via ALDH1A1

Given the important role of ALDH in ROS clearance, a high level of ALDH could also be important for the development of KRAS-driven PDAC. Here, we analyzed the expression of ALDH family members in normal pancreas and PDAC samples (Bailey et al., 2016; GTEx Consortium, 2013). In normal pancreas tissues, we mainly observed the expression of ALDH1A1 (Figure 4—figure supplement 1A), with different cell types exhibiting different expression levels of ALDH1A1 (Figure 4—figure supplement 1B). Since the tumor cells are mainly epithelial cells, we only compared PDAC data to pancreatic ductal cells to avoid the confounding factors caused by the cell type difference. As shown in Figure 4—figure supplement 1B, there are four subclusters of normal ductal cells. The average expression level of ALDH1A1 in normal pancreatic ductal cells (clusters 1–3) is less than 50. We excluded cluster 4 since the ALDH1A1-positive cells are indicative of the ductal stem cell population (Rovira et al., 2010). In contrast to the expression levels in normal ductal cells, we observed that in 63% of PDAC samples, the expression levels of ALDH1A1 are higher than 50 TPM, and in 10% of samples, the expression levels are higher than 200 TPM (Figure 4 and Figure 4—figure supplement 1C). Furthermore, we examined the mutation levels in ALDH1A1. We observed that only 0.2% of the patients (1 out of 576 patient samples from two cohorts [Bailey et al., 2016; Cancer Genome Atlas Research Network, 2017]) acquired mutations in ALDH1A1 (Figure 4B). This observation further supports our hypothesis that ALDH1A1 plays an important role in KRAS-driven PDAC development.

Figure 4 with 2 supplements see all
Download asset Open asset
ARID1A knockout facilitates escape from KRAS-induced senescence by upregulating ALDH1A1 expression.

(A) Heatmap of the expression levels of aldehyde dehydrogenase (ALDH) family members in pancreatic ductal adenocarcinoma (PDAC) patients (Bailey et al., 2016). (B) Mutation rates of ALDH1A1, ALDH1A3, and ALDH3A1 in PDAC patients. The mutation data from these two studies (Bailey et al., 2016) and TCGA (Cancer Genome Atlas Research Network, 2017) were used for analysis. (C) Colony formation assay of ARID1A knockout cells and wildtype cells with KRAS induction, treated with and without ALDH inhibitor DEAB. (D) Quantification of colony number in panel (C). The colony formation assay was performed twice. (E) Colony formation assay of ARID1A knockout cells expressing shRNA targeting ALDH1A1 and scramble shRNA control. (F) Quantification of colony number in panel (E). The colony formation assay was performed twice. (G) Measurement of the reactive oxygen species (ROS) level using an H2DCFDA-based ROS detection assay kit. Percentage of positive cells measured by flow cytometry. (H) Measurement of the ROS levels in the ARID1A-KO cells with ALDH1A1 knockdown or scramble shRNA. (I) Working model for ARID1A-deficiency-promoted pancreatic intraepithelial neoplasia (PanIN) lesion progression via inhibition of ROS production. Student’s t-test: *p<0.05; **p<0.01; ***p<0.001.

Next, to validate the essential role of ALDH1A1 in promoting the escape of cells from KRAS-induced senescence, we performed a colony formation assay in HPNE cells with and without N,N-diethylaminobenzaldehyde (DEAB, a pan-inhibitor of ALDH) treatment. We observed that inhibition of ALDH1A1 activity significantly decreased the number of colonies formed in ARID1A knockout cells; in contrast, no significant changes were observed in the wildtype cells (Figure 4C,D). To rule out the unknown effects of DEAB on HPNE cells, we also performed a colony formation assay on ARID1A-KO HPNE cells with and without ALDH1A1 knockdown. The knockdown efficiency was verified by qRT-PCR (Figure 4—figure supplement 2). We also observed that the colony number in ARID1A-KO cells with ALDH1A1 knockdown was significantly less than that without ALDH1A1 knockdown (Figure 4E,F), which is consistent with the results of the ALDH inhibitor experiment.

Furthermore, we examined the levels of ROS production in ARID1A-KO cells and wildtype cells. We observed that the fraction of ROS-positive cells in ARID1A-KO iKRAS-HPNE cells was significantly less than in wildtype cells, regardless of KRAS induction (Figure 4G). To verify the role of ALDH1A1 in reducing ROS, we measured the ROS level in ARID1A-KO cells with ALDH1A1 knockdown. We observed that the ROS level in ALDH1A1 knockdown cells was significantly higher than that observed in the cells with scramble shRNA (Figure 4H), indicating that ALDH1A1 is responsible for the clearance of ROS in HPNE cells with ARID1A knockout. From the results above, we conclude that ARID1A knockout can effectively reduce cellular ROS levels by upregulating the expression of ALDH1A1, which then leads to significant attenuation of KRAS-induced senescence and acceleration of PanIN progression (Figure 4I).

ARID1A KO alters gene expression mainly by modulating the chromatin accessibility of distal regulatory elements

As one of the DNA binding subunits of the SWI/SNF complex, ARID1A can regulate gene expression by modulating the chromatin accessibility for transcription factor (TF) binding, recruitment of coactivators/corepressors, and facilitation of chromatin looping required to approximate promoters with distal enhancers (Wu and Roberts, 2013). To investigate how ARID1A deficiency enhances the expression of ALDH1A1, we performed ATAC-seq on ARID1A-KO and wildtype HPNE cells. The distribution of the fragment sizes and the distance from transcription start sites for both wildtype and ARID1A-KO are shown in Figure 5—figure supplement 1A–F (Supplementary file 5). As shown in Figure 5A, the Spearman correlation coefficients between knockout and wildtype cells are significantly lower than those between replicate samples in both groups (Figure 5B and Figure 5—figure supplement 1G), suggesting that the changes in DNA accessibility were robustly captured.

Figure 5 with 6 supplements see all
Download asset Open asset
ARID1A knockout activates transcription of the ALDH1A1 gene by increasing the accessibility of its enhancer region.

(A) Spearman correlation coefficients of the read counts in peaks between ARID1A-KO human pancreatic Nestin-expressing (HPNE) cells and wildtype cells. ATAC sequencing was performed with three biological repeats. (B) The average fold change of chromatin accessibility in the promoters and enhancers of differentially expressed genes identified in RNA-seq. Mann–Whitney–Wilcoxon test: ****p<0.0001. (C, D) The scatter plot of the read counts in each peak between ARID1A-KO cells and wildtype cells for promoter (C) and distal regions (D). The peaks with significantly increased read density in ARID1A-KO cells compared with wildtype cells are colored in red. The peaks with significantly decreased read density are colored in blue. (E) Heatmap of the gained or lost distal regions between ARID1A-KO cells and wildtype cells. (F) The top five transcription factor (TF) binding motifs significantly enriched in the enhancer ATAC peaks with increased accessibility in ARID1A-KO cells. (G) The top five TF binding motifs significantly enriched in the enhancer ATAC peaks with decreased accessibility in ARID1A-KO cells.

Next, we examined the accessibility of the promoter and enhancer regions of the genes that are differentially expressed between ARID1A-KO cells and wildtype cells. We separated the DEGs into two groups based on the fold change. The genes with positive fold-change values are the genes upregulated in ARID1A-KO cells (upregulated genes group in Figure 5B), and the genes with negative fold-change values are the genes downregulated in ARID1A-KO cells (downregulated genes group in Figure 5B). We plotted a scatter plot of read counts in peaks between wildtype and ARID1A-KO for promoters (Figure 5C) and enhancers (Figure 5D) and observed that the number of peaks affected by ARID1A deficiency in the distal regulatory regions is significantly larger than the number of peaks in promoter regions.

The heatmap of reads for the differential peaks is shown in Figure 5E, and the average read density profiles are shown in Figure 5—figure supplement 2A,B. For the differential peaks of promoters, the heatmaps of reads and the average read density profiles are shown in Figure 5—figure supplement 2C,D. We also performed analysis for the distribution of the differential peaks (Figure 5—figure supplement 3A,B) and a functional enrichment analysis (Figure 5—figure supplement 3C,D) using the GREAT algorithm (McLean et al., 2010). We observed consistent enrichment in enhancer regions. Significant interactions between the SWI/SNF complex and distal regulatory regions have also been observed in colon cancer (Mathur et al., 2017).

Next, we examined the association between DEGs and the peaks with differential accessibility. We observed that the number of DEGs is associated with peak changes with statistical significance for both promoters and enhancers (Figure 5—figure supplement 4A–D). We also noticed that the number of DEGs associated with peak changes in enhancer elements is significantly larger than the number of genes associated with peak changes in the promoter regions. This result further supports that ARID1A knockout alters gene expression mainly by modulating the chromatin accessibility of the enhancer elements.

Considering the observation that AR1D1A deficiency impairs the activities of KRAS signaling pathways based on the GSEA of transcriptome data (KRAS_SIGNALING_UP in Figure 3—figure supplement 1A), next we examined the chromatin accessibility of the genes involved in KRAS signaling. We observed that in comparison to the wildtype cells, chromatin accessibility decreased in ARID1A-KO cells (Figure 5—figure supplement 5). This observation suggests that AR1D1A deficiency impairs the activities of the KRAS signaling pathways by partially impairing the chromatin accessibility of the genes downstream of the KRAS pathways.

We next performed motif enrichment analysis for the differential ATAC peaks. We separated the ATAC peaks with significant changes between ARID1A-KO and wildtype cells into four groups: distal peaks with increased accessibility in ARID1A-KO cells, distal peaks with decreased accessibility in ARID1A-KO cells, promoter peaks with increased accessibility in ARID1A-KO cells, and promoter peaks with decreased accessibility in ARID1A-KO cells. We then performed motif enrichment analysis by using the AME algorithm (McLeay and Bailey, 2010) on each group of ATAC peaks. Interestingly, we found that the binding motifs of SRY, FOX family, CDX1, SOX5, etc., were enriched in the ATAC peaks with significantly increased accessibility in ARID1A-KO cells (Figure 5F, Figure 5—figure supplement 6A,B and Supplementary file 7). The binding motifs of the FOS-JUN family, NFE2, NF2L1, etc., were significantly enriched in the ATAC peaks with significantly decreased accessibility (Figure 5G, Figure 5—figure supplement 6B, and Supplementary file 7), which is consistent with the results of a recent study showing ARID1A as a co-factor of AP-1 (Sen et al., 2019).

ARID1A KO activates transcription of the ALDH1A1 gene by increasing the accessibility of its enhancer region

To understand how the gene expression of ALDH1A1 is altered by ARID1A knockout, we examined the accessibility of the promoter and distal regulatory elements for the ALDH1A1 gene in both ARID1A-KO and wildtype cells. Interestingly, we found that there was a significant increase in accessibility in 9 out of 11 peaks at the distal regions when we compared ARID1A-KO cells with wildtype cells (Supplementary file 6). To further verify these functional regions, we compared the landscape of two well-known enhancer markers (H3K27ac and H3K4me1) in seven highly expressed ALDH1A1 cell lines (ALDH1A1high cell lines) (Mei et al., 2017) with our ATAC-seq peaks (Supplementary file 8). We observed that five enhancer loci clearly overlapped with the active enhancer markers (Figure 6A and Figure 6—figure supplement 1A,B).

Figure 6 with 2 supplements see all
Download asset Open asset
ARID1A knockout activates transcription of the ALDH1A1 gene by increasing the accessibility of its enhancer region.

(A) The ATAC-seq tracks and H3K4me1/H3K27ac ChIP-seq tracks of the distal regions of the ALDH1A1 gene. The ChIP-seq tracks from different cell lines are labeled in different colors and overlaid. The figure with separated tracks is shown in Figure 6—figure supplement 1. The ALDH1A1high cell lines include A549, LOUCY, A673, 22RV1, VCAP, K562, and HepG2. The ALDH1A1low cell lines include HCT116, MCF7, Panc1, PC9, and DOHH2. The ChIP-seq data were obtained from the ENCODE database (ENCODE Project Consortium, 2012). (B) Read counts in the five enhancer peaks of the ALDH1A1 gene in ARID1A-KO human pancreatic Nestin-expressing (HPNE) cells and wildtype cells. p-value: *p<0.05; ***p<0.001; ****p<0.0001. (C) Transcription factors (TFs) that regulate ALDH1A1 expression by binding to its enhancer loci. Left panel: the list of TFs that can potentially bind to the enhancer regions of ALDH1A1 based on the ChIP-seq data from seven ALDH1A1high cell lines and five ALDH1A1Low cell lines. The TFs whose binding events are preferentially detected in the datasets from ALDH1A1high cell lines (Fisher’s test, p<0.05) are marked in red. The TFs that are observed in more than 50% of the datasets from ALDH1A1high cell lines but do not have enough datasets for statistical tests are marked in blue. The TFs whose binding events are observed in more than 25% but less than 50% of the datasets from ALDH1A1high cell lines are marked in black. The TFs that are not expressed in HPNE cell lines were removed. Middle panel: expression of ALDH1A1 in EP300 knockdown cells and cells with scramble shRNA. Right panel: expression of ALDH1A1 in NR3C1 knockdown cells and cells with scramble shRNA.

We also plotted out the landscape of H3K27ac and H3K4me1 in five lowly expressed ALDH1A1 cell lines (ALDH1A1low cell lines) (Mei et al., 2017), and we did not observe any significant overlap between our ATAC-seq peaks and H3K27ac/H3K4me1-enriched regions (Figure 6A). For the five enhancer loci, we observed a consistent increase in accessibility in ARID1A knockout cells in comparison with wildtype cells (Figure 6B). Overall, these results confirm that ARID1A deficiency upregulates ALDH1A1 expression by increasing the accessibility of the associated enhancer elements.

To identify proteins that could potentially bind to these five enhancer loci, we examined the TF ChIP-seq datasets from seven ALDH1A1high cell lines and five ALDH1A1low cell lines (Mei et al., 2017). We counted the binding events of each TF in the five candidate enhancer peaks for both ALDH1A1high cell lines and ALDH1A1low cell lines (Supplementary file 9). We list the TFs that could bind to the enhancer regions in Figure 6C. The TFs whose binding events were preferentially detected in the datasets from ALDH1A1high cell lines (Fisher’s test, p<0.05) are indicated in red and include EP300 and NR3C1. For the TFs that do not have enough datasets for a statistical test, we indicated the TFs whose binding events are observed in more than 50% of ALDH1A1high cell line datasets in blue, and the TFs whose binding events are observed in more than 25% but less than 50% of ALDH1A1high cell line datasets in black (Figure 6C, left panel). The TFs that are not expressed in HPNE cell lines were removed.

To verify the regulation of ALDH1A1 expression by the TFs identified above, we knocked down EP300 and NR3C1, respectively, in ARID1A-KO HPNE cell line (Figure 6—figure supplement 2A,B). In Figure 6C (middle and right panels), we observed that both knockdowns could significantly impair the transcription of the ALDH1A1 gene. Collectively, these data also indicate that ARID1A deficiency promotes active transcription of ALDH1A1 by altering genome accessibility and therefore allowing the binding of EP300 and/or NR3C1 to the corresponding enhancer loci.

Discussion

As a bona fide oncogenic driver of pancreatic cancer, the mechanism by which mutant KRAS counteracts the oncogenic stress that it induces is critically important for the progression of precursor lesions (Storz, 2017). It has been shown that KRAS can upregulate the expression of multiple oxidoreductases via NRF2 (nuclear factor, erythroid derived 2, like 2) to counteract oncogenic stress (DeNicola et al., 2011). Meanwhile, KRAS can also enhance NADPH production by reprogramming the metabolism of glutamine to assist oxidoreductases in scavenging ROS (Son et al., 2013). However, the observation that a high percentage of senescent PanIN lesions occur in the pancreas of KC mice suggests that the mechanisms described above likely have limited effects in reducing Kras-induced ROS in neoplastic lesions. In the early stage of PDAC development, it is plausible that new mutations can be acquired to help cells to reduce ROS and escape the oncogenic Kras-induced senescence.

In this study, we successfully adapted a high-sensitivity total-RNA-based single-cell RNA-seq method, MATQ-seq, to profile the transcriptome of single lesions. We showed that transcriptome profiling of individual lesions is not only technically feasible but also it can provide important insights about the potential mechanisms of tumor progression. The transcriptome profiling of early lesions directly led us to unveil the effects of ARID1A knockout in attenuating KRAS-induced senescence and identify the important roles of ALDH1A1 in mitigating the ROS stress induced by oncogenic KRAS.

It is worth noting that ARID1A has also been linked to the regulation of ROS through other pathways. Sun et al. found that ARID1A overexpression causes an increase in ROS by activating transcription of cytochrome P450 enzymes (CYP450) at the initiation stage of liver cancer (Sun et al., 2017). Interestingly, Ogiwara et al. found that ARID1A deficiency results in elevated ROS by inhibiting the transcription of SLC7A11 (a transporter gene required for the import of cystine and the production of glutathione) in ovarian cancer (Ogiwara et al., 2019). Therefore, the effects of ARID1A deficiency on ROS are likely tissue-specific.

Besides reducing ROS levels through the upregulation of ALDH family proteins, we also observed that ARID1A deficiency directly suppresses the activities of KRAS signaling. Our in vitro and in vivo data showed that the activation of Hallmark gene sets KRAS_SIGNALING_UP and KRAS_SIGNALING_DN is remarkably impaired upon ARID1A knockout. Consistently, the activation of several inflammation-related signaling pathways upon Kras induction is inhibited in PanIN lesions and HPNE cells with ARID1A deficiency compared with wildtype counterparts, including TNFα signaling, IL6–STAT3 signaling, IL2–STAT5 signaling, IFN-α signaling, and IFN-γ signaling. Both reduction of KRAS activities and inflammatory cytokine response are also consistent with our observation of the increased attenuation of Kras-induced senescence upon ARID1A knockout.

Overall, experimental evidence shows that mutant KRAS signaling needs to be tuned to efficiently drive tumorigenesis: too much signaling could lead to cell growth arrest, while too little signaling could limit proliferation (Li et al., 2018). To find the ‘sweet spot’ of KRAS activation, KRAS could harness its own activities by changing the mutation pattern (including mutation position and type of substitution) and seeking cooperation with other genes (e.g., mutagenesis of tumor suppressors) through the course of passive selection during tumor development (Li et al., 2018). Overall, our study demonstrates the important roles of epigenetic mutations in harnessing KRAS activities and promoting tumorigenesis. Future studies of such effects in other frequently mutated epigenetic genes are greatly desired.

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Genetic reagent (Mus musculus)Ptf1aCreERTJackson LabNo: 019378
Genetic reagent (M. musculus)Lox-Stop-Lox-KrasG12DJackson LabNo: 008179
Genetic reagent (M. musculus)Arid1afloxPMID:18448678Dr. Zhong Wang lab (University of Michigan)
Cell line (Homo sapiens)hTERT-HPNEDr. Jennifer Baily lab (UT Health)RRID:CVCL_C466
Recombinant DNA reagentpInducer20-KRASG12DDr. Haoqiang Ying lab (MD Anderson Cancer Center)Inducible expression of KRASG12D
Recombinant DNA reagentpL-CRISPR.EFS.tRFPAddgeneRRID:Addgene_57819Pol III-based sgRNA expression backbone
Recombinant DNA reagentpGIPZOpen BiosystemsPol III-based shRNA expression backbone
Sequence-based reagentsgARID1AThis papersgRNACAGCGGTACCCGATGACCAT
Sequence-based reagentshALDH1A1This papershRNAGGAGTGTTTACCAAAGACATT
Sequence-based reagentshEP300-1This papershRNACGGCAAACAGTTGTGCACA
Sequence-based reagentshEP300-2This papershRNAAGCTACTGAAGATAGATTA
Sequence-based reagentshNR3C1-1This papershRNACCAACGGTGGCAATGTGAA
Sequence-based reagentshNR3C1-2This papershRNAAGCTGTAAAGTTTTCTTCA
Commercial assay or kitROS detection assay kitBioVisionCat# K936-250
Commercial assay or kitSenescence β-Galactosidase Staining KitCell Signaling TechnologyCat# 9860
Chemical compound, drugALDEFLUOR diethylaminobenzaldehyde (DEAB) reagent, 1.5 mM in 95% ethanolStemcell Technologies IncCat# 01705
Antibody(Rabbit polyclonal) anti-ALDH1A1AbcamAbcam Cat# ab23375, RRID:AB_2224009WB(1:400)
Antibody(Rabbit polyclonal) anti-Aldh3a1AbcamAbcam Cat# ab76976, RRID:AB_1523110IHC(1:100)
Antibody(Rabbit polyclonal) anti-RasAbcamAbcam Cat# ab180772, RRID:AB_2884935WB(1:500)
Antibody(Mouse monoclonal) anti-β-ActinSigma-AldrichSigma-Aldrich Cat# A1978, RRID:AB_476692WB(1:4000)
Antibody(Rabbit monoclonal) anti- Phospho-p44/42 MAPK (Erk1/2)Cell Signaling TechnologyCell Signaling Technology Cat# 4370, RRID:AB_2315112WB(1:1000)

Mice

All animal experiments in this study were performed in accordance with a protocol approved by IACUC of Baylor College of Medicine. In this study, the following mice strains were generated: Arid1afl/fl; Lox-Stop-Lox-KrasG12D/+;Ptf1aCreERT/+, Arid1afl/+;Lox-Stop-Lox-KrasG12D/+;Ptf1aCreERT/+, and Arid1a+/+;Lox-Stop-Lox-KrasG12D/+;Ptf1aCreERT/+. The mouse with the target allele type of Lox-Stop-Lox-KrasG12D/+ was ordered from Jackson Laboratory (Stock No: 008179 | LSL-K-ras G12D). The promoter of Ptf1a is used to drive the expression of inducible Cre in adult acinar cells (Stock No: 019378 | Ptf1aCreERT/+, denoted as C). Removal of the floxed exon 8 of the Arid1a (denoted as Afl) leads to loss of Arid1a transcript (Gao, 2008). Oncogenic Kras mutant allele is silenced by a floxed STOP cassette (Lox-Stop-Lox-KrasG12D/+, denoted as K). The adult mice (6–8 weeks) of the above genotypes were administered with tamoxifen (75 mg/kg/day for five consecutive days) to induce efficient ablation of Arid1a and activation of oncogenic Kras in pancreatic acinar cells.

Cell lines

Request a detailed protocol

We received the HPNE cell line (RRID:CVCL_C466) from Dr. Jennifer Bailey’s lab (UT Health Center), and the line is cultured in media recommended by ATCC. The identity of the cell line has been authenticated by STR profiling, and the mycoplasma contamination in cell lines has been routinely tested. HPNE with inducible KRASG12D/+ (iKRAS-HPNE) was generated by transduction of plasmid pInducer20-KRASG12D into parental HPNE cells. After G418 selection for 15 days, the survived cells were then used for single-cell expansion in 96-well plate. For ARID1A knockout HPNE cells, isogenic clone of iKRAS-HPNE was infected with lentivirus packaged with pL-CRISPR.EFS.tRFP-sgARID1A. After infection for 5 days, the RFP-positive cells were sorted and single-cell expansion was performed. ARID1A knockout was confirmed by Sanger sequencing of the isogenic clones. The guide RNA sequence used for ARID1A knockout is CAGCGGTACCCGATGACCAT. For ALDH1A1 knockdown cells, two ARID1A knockout HPNE clones were infected by lentivirus packaged with plasmid pGIPz-ALDH1A1. After infection for 5 days, the GFP-positive cells were sorted. Knockdown efficiency was confirmed by RT-PCR. The shRNA sequence targeting ALDH1A1 is GGAGTGTTTACCAAAGACATT. For knockdown of C/EBPβ, EP300, or NR3C1, HPNE cells with ARID1A knockout were transfected with plasmid pGIPz-EP300 or pGIPz-NR3C1. After infection for 5 days, the GFP-positive cells were sorted. Knockdown efficiency was confirmed by RT-PCR. The shRNA sequences targeting NR3C1 are CCAACGGTGGCAATGTGAA and AGCTGTAAAGTTTTCTTCA; and targeting EP300 are CGGCAAACAGTTGTGCACA and AGCTACTGAAGATAGATTA.

Primary acinar isolation and culture

Request a detailed protocol

Pancreata isolated from KC mice and AKC mice 1 month after 5-day tamoxifen administration were rinsed twice in cold 1× HBSS buffer, then minced into small pieces and digested with digestion buffer (HBSS buffer with 10 mM HEPES and 0.5 mg/ml collagenase and 0.25 mg/ml trypsin inhibitor) for 20–30 min at 37°C. During the incubation. the tissue was pipetted every 5 min. After washing twice with washing buffer (HBSS buffer with 5% FBS), the digested tissue was resuspended with media (Waymouth media with 2.5% FBS and 0.25 mg/ml trypsin inhibitor and 100 U/ml Penicillin-Streptomycin), filtrated with 100 µm strainer and seeded into 10 cm dishes overnight at 37°C to remove fibroblasts and ductal cells. The unattached acinar cells were then transferred into collagen-coated plates for growth. To activate Kras expression, 25 ng/ml EGF was added into the media for 5 days. Cells were then used for SA-β-Gal staining, RT-PCR, and western blot.

Western blot

Request a detailed protocol

Western blot was performed using the standard protocol. Antibodies used in this study include ALDH1A1 (Abcam Cat# ab23375, RRID:AB_2224009), ALDH3A1 (Abcam Cat# ab76976, RRID:AB_1523110), KRAS (Abcam Cat# ab180772, RRID:AB_2884935), phosph-Erk1/2 (Cell Signaling Technology Cat# 4370, RRID:AB_2315112), and β-actin (Sigma-Aldrich Cat# A1978, RRID:AB_476692).

Colony formation assay

Request a detailed protocol

HPNE cells (3 × 104/well) were seeded into 6-well plates. The cells were treated with doxycycline (6 µg/ml for 15 days) with and without DEAB (1.5 µM for 30 days, Stemcell Technologies Inc 01705). The media was changed every 2 days. When the colonies were large enough, the cells were fixed and stained with crystal violet.

ROS measurement

Request a detailed protocol

HPNE cells were treated with doxycycline (6 µg/ml for 5 days). ROS level was measured by Flow Cytometry using ROS Detection Assay Kit (BioVision, Cat# K936-250).

SA-β-Gal staining

Request a detailed protocol

SA-β-Gal staining was performed on slides of freshly frozen tissues or cells using Senescence β-Galactosidase Staining Kit (Cell Signaling Technology, Cat# 9860). Total and SA-β-Gal-positive lesions or cells were counted at random fields under the microscope, and positive rates were calculated. For quantification of SA-β-Gal staining of primary acinar cells, due to the difficulty of recognizing the nuclei, average optical density (OD) was used to quantify the intensity of SA-β-Gal staining. 8-bit images were adjusted for white balance and color-deconvoluted using Feulgen light green vector in ImageJ. The average gray values of the green channel were measured. OD was calculated using the following formula: OD = log10 (255/gray value).

PanIN quantification

Request a detailed protocol

Tissues were fixed in 4% paraformaldehyde overnight, processed, and embedded in paraffin. Paraffin-embedded sections were subjected to hematoxylin and eosin staining (H&E staining). ADM, mPanIN-1A, mPanIN-1B, mPanIN-2, and mPanIN-3 were quantified using ImageJ for morphometric analysis based on scanned H&E slides at 20× magnification. Grades of lesions were determined based on the characteristic morphology criterion (Gopinathan et al., 2015) with consulting of pathology core. PanIN-1A: flat epithelium composed of columnar cells with basally oriented nuclei. PanIN-1B: identical to PanIN-1A lesions but exhibit papillary or basally pseudostratified architecture. PanIN-2: show mild nuclear abnormalities, including loss of polarity, nuclear enlargement, nuclear crowding, and nuclear pleomorphism. PanIN-3: show a predominantly papillary or micropapillary architecture with abnormal cribriforming, budding, and luminal necrosis; more severe cytological atypia, such as loss of nuclear polarity, dystrophic goblet cells, nuclear irregularities, and macro nucleoli. In case the lesion harbors more than one PanIN grade, the lesion was graded based on the component with the highest grade. Numbers of lesions of different grades were counted for at least five fields of view. The area of tissue was measured for each field of view. Lymph nodes of the pancreatic area were excluded. Numbers of lesions and tissue areas were summed up to calculate lesion number per area.

IHC quantification

Request a detailed protocol

For quantification of IHC results against ALDH3A1, H-score method was used. In brief, staining intensity (not stained: 0; weakly stained: +1; moderately stained: +2; or strongly stained: + 3) was determined for each lesion of interest in the field. The H-score was calculated by the following formula: 3 × percentage of strongly stained cells + 2 × percentage of moderately stained cells + 1 × weakly stained cells, giving a range of 0–300.

Bulk RNA-seq

Request a detailed protocol

HPNE cells were treated with doxycycline (6 µg/ml) for 5 days. RNA samples were prepared using the standard protocol for Trizol. mRNA was enriched using NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB, E7490), and the library was prepared using the NEBNext Ultra II RNA Library Prep Kit for Illumina (NEB, E7770). All libraries were sequenced on Illumina Nextseq500 platform. Reads were aligned to hg19 assembly of the human genome by STAR aligner (Dobin et al., 2013), and transcripts counting was performed by HTseq-count (Anders et al., 2015). Differential gene expression analysis was performed by using edgeR (Robinson et al., 2010) with a cutoff of FDR at 0.05. To identify the genes with differential response to oncogenic KRAS in KO and WT cells, we also performed the interaction analysis in edgeR.

Analysis of ALDH1A1 expression in normal pancreas and PDAC

Request a detailed protocol

The expression profiles of ALDH genes in normal pancreas were obtained from GTEx database. The expression level of ALDH1A1 in different cell types in normal pancreas was obtained from Human Protein Atlas database. The PDAC RNA-seq data were from ICGC-PACA-AU cohort. The raw count data were downloaded from https://dcc.icgc.org/https://dcc.icgc.org/https://dcc.icgc.org/https://dcc.icgc.org/.

ATAC-seq experiment

Request a detailed protocol

ATAC-seq was performed following the protocol of Howard Chang’s lab (https://www.nature.com/articles/nmeth.4396) with slight modifications. In brief, 5 × 104 cells were lysed with ATAC-Resuspension Buffer (RSB) containing 0.1% NP40 and 0.1% Tween-20. After incubation on ice for 3 min, the cell lysates were washed by RSB with 0.1% Tween-20. The cell lysates were then incubated with transposition mixture at 37°C for 30 min. After amplification, the transposed fragments were purified with magnetic beads. Finally, 4 ng fragments were used for the generation of the library. All libraries were sequenced on Illumina Nextseq500 platform.

ATAC-seq data analysis

Request a detailed protocol

Reads were then mapped to the hg19 assembly by Bowtie2 (Langmead and Salzberg, 2012) after removing the adaptor sequence. The quality control of ATAC-seq data was performed by using the ATACseqQC R package (Ou et al., 2018). Next, the mapped reads from three technical replicates of each genotype were combined for the peak calling by MACS2 (Zhang et al., 2008). Peaks from wildtype samples and ARID1A-KO samples were combined to get a union peak set. All the peaks were then annotated by HOMER (Heinz et al., 2010). HTseq-count (Anders et al., 2015) was used for read counting. edgeR package (Robinson et al., 2010) was then used for normalization between different samples and for peak differential analysis. The read density profiles of the differential peaks were plotted by deepTools (Ramírez et al., 2016). The motif enrichment analysis was performed by using AME.

ENCODE data analysis

Request a detailed protocol

The expression levels of ALDH1A1 in the cell lines in ENCODE were obtained from the Cancer Cell Line Encyclopedia (CCLE) database (file name: file-CCLE_RNAseq_genes_rpkm_20180929.gct.gz; Ghandi et al., 2019). Seven cell lines with RPKM >10 were categorized as ALDH1A1high cell lines. Five cell lines with RPKM <0.5 were selected as ALDH1A1low cell lines. The A549 H3K27ac and H3K4me1 ChIP-seq data were reanalyzed by using hg19 assembly as described above. The peak files of H3K27ac/H3K4me1 from the other cell lines were directly obtained from the ENCODE database. The ALDH1A1-related differential ATAC peaks, which are overlapped with the H3K27ac/H3K4me1 peaks in at least four out of seven ALDH1A1high cell lines, were characterized as functional enhancer regions.

PanIN-seq

Request a detailed protocol

We used the MMI CellCut platform to perform LCM. 40–50% laser power was used with the cutting speed of 18 µm/s to dissect microscopic lesions. 1.6 µl of MATQ-seq lysis buffer (1 μl of 0.2% Triton X100 [Sigma-Aldrich], 0.4 μl of primer mix, 0.12 μl dNTP, 0.05 μl 0.1 M DTT [Life Technologies], and 2 U RNaseOUT [Life Technologies]) was added onto the isolation cap where the dissected tissue was attached (MMI, Prod. No. 50206). We used a pipette tip to scrape the laser-dissected tissue into the lysis buffer and then pipetted the lysis buffer into the tube. Sample tubes were then placed on a thermocycler and incubated at 72°C for 3.5 min, followed by 1 min incubation on ice. 2.4 µl of MATQ-seq first strand synthesis buffer (0.8 μl 5× First Strand Buffer [Life Technologies], 0.2 μl 0.1 M DTT, 4 U RNaseOUT [Life Technologies], 30 U Superscript III [Life Technologies], and 1.15 μl RNase-free water) was then added. The reverse transcription program was same as MATQ-seq. After reverse transcription, the residual primers were then digested by using T4 polymerase at 37° for 40 min and 75° for 20 min. RNA was then digested by using RNase-H and RNase-If at 37° for 15 min and 72° for 15 min. Following that, dC-tailing and second strand synthesis were performed as described in MATQ-seq. The library prep for PanIN-seq samples were same as MATQ-seq. All libraries were sequenced on Illumina Nextseq500 platform.

PanIN-seq data analysis

Request a detailed protocol

The raw sequencing data trimming and barcode retrieval were performed as previously described (Sheng et al., 2017). The reads were mapped to the genome MM10 using STAR with the following parameters: --outFilterMismatchNoverLmax 0.05 --outFilterMatchNmin 16 --outFilterScoreMinOverLread 0 --outFilterMatchNminOverLread 0. We used Gencode annotation release mm10 (GRCm38.p4) for transcript annotation. Unique barcode counting and gene expression-level quantification were performed as previously described with a few modifications: the mapping position of the reads was included as part of the identity of the corresponding barcodes; only reads mapped to the exon region were used for gene expression-level quantification. The genes with APM >2 in at least five samples were retained for differential gene expression analysis, which was performed by using edgeR with a cutoff of FDR at 0.05.

Functional enrichment analysis

Request a detailed protocol

Functional enrichment analysis was performed by using GSEA with the MSigDB Hallmark gene sets, and the senescence-related gene sets from the MSigDB curated gene sets. For PanIN-seq, the permutation of phenotype label was used to calculate the p-value.

Data availability

Sequencing data have been deposited in GEO under GSE160444.

The following data sets were generated
    1. Liu S
    2. Cao WJ
    3. Niu YC
    4. Luo JY
    5. Zhao JH
    6. Zong C
    7. Hu Z
    (2020) NCBI Gene Expression Omnibus
    ID GSE160444. ARID1A deficiency promotes pancreatic tumorigenesis by reducing KRAS-induced senescence.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE160444
The following previously published data sets were used
    1. The International Cancer Genome Consortium
    (2016) European Genome-phenome Archive (EGA)
    ID EGAS00001000154. Genomic analyses identify molecular subtypes of pancreatic cancer.
    https://www.ebi.ac.uk/ega/studies/EGAS00001000154
    1. The Cancer Genome Atlas Research Network
    (2017) GDC Data Portal
    ID TCGA-PAAD. Integrated Genomic Characterization of Pancreatic Ductal Adenocarcinoma.
    https://portal.gdc.cancer.gov/projects/TCGA-PAAD

References

Decision letter

  1. Maureen E Murphy
    Senior and Reviewing Editor; The Wistar Institute, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The authors have identified the importance of ALDH1A1 in preventing oncogene-induced senescence in precursor lesions of pancreatic cancer. Using a combination of sequencing, computational analysis and experimental work they find that ALDH1A is negatively regulated by ARID1A and that loss of ARID1A results in its upregulation and a decrease in intracellular ROS. This highlights a mechanism whereby PanINs prevent senescence and drive a pro-growth phenotype.

Decision letter after peer review:

Thank you for submitting your article "Single-PanIN-seq Unveils that ARID1A Deficiency Promotes Pancreatic Tumorigenesis by Attenuating KRAS Induced Senescence" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by Maureen Murphy as the Senior and Reviewing Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission. Please be advised that there was considerable discussion about your work, and that considerable revisions are requested by the reviewers. You may wish to create a plan of attack for the revisions and send this to the Senior Editor, who can share it with reviewers.

As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is "in revision at eLife". Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

Summary:

The authors have identified the importance of ALDH1A1 in preventing oncogene-induced senescence in precursor lesions of pancreatic cancer. By a combination of sequencing, computational analysis and experimental work they find that ALDH1A is negatively regulated by ARID1A and that loss of ARID1A results in its upregulation and a decrease in intracellular ROS. This highlights another new mechanism that is exploited by PanINs to prevent senescence and drive a pro-growth phenotype.

Essential revisions:

1. The conclusion that "a highly complex trans-differentiation occurs following the knockout of Arid1a" is not supported by the evidence provided in Figure S3. It is unclear what panels B and C are referring to (no labels) and there is no scale given. The authors also use the RNA-sequencing data to conclude that "the activities of Kras signaling are partially suppressed by ARID1A deficiency." This is a major point and should be supported by experimental evidence in their cell models.

2. Although Figure S5B shows that genetic manipulation of ARID1A was successful, there is no evidence by western blot or PCR that ARID1A has truly been knocked out in their HPNE cells. This model was used for a large portion of the data and must be characterized in a convincing way.

3. The evidence for senescence in vitro is not strong. Several times throughout, a colony formation assay is used as a measure of senescence, rather than SA-b-gal, molecular markers or secreted factors. It is unclear why the authors measured SA-b-gal staining intensity rather than %-positive cells. In particular, the evidence that ARD1A regulates Kras-induced senescence in HPNE cells is lacking.

4. It is unclear if the ROS phenotype is mediated by ALDH1A1. Therefore, the conclusion that "ARID1A knockout can effectively reduce cellular ROS level by upregulating the expression of ALDH1A1" is not supported.

5. There is no functional experimental evidence provided for how ARID1A regulates the expression of ALDH1A1. This is a major point of novelty in the paper and should therefore be supported by experimental evidence.

6. In the patients with pancreatic ductal adenocarcinoma, typically one sees point mutations, InDel or structure variations in the cancer genes; can the authors list the mutations in ARID1A in their samples to see if there are any recurrent mutations (i.e, visualized by lollipop plot, using ICGC-PACA-AU or ICGC-PACA-CA whole genome sequencing dataset, https://dcc.icgc.org/)? Can the authors perform an eQTL analysis (for example using FastQTL) to see if any mutations in ARID1A and +/- 1Mb cis-mutations are associated with the expression of ALDH1A1? The mutations and gene expression datasets could be readily downloaded from ICGC-PACA-AU or ICGC-PACA-CA (ICGC data portal).

7. The authors describe differences in proportions of distinct populations of ADM, PanIN-1, PanIN-2, PanIN-3 populations. At no point in their introduction/results do they describe the PanIN-1/2/3 populations. This would be a good addition for clarity. In addition there are no details on how these were quantified in the methods. In the introduction the authors state that previously ARID1A KO has been found to increased ADM – however Figure S1 shows the opposite that ARID1A KO leads to decreased ADM – while this reflects a trajectory of PanIN progression – not including any background on this is confusing and needs clarification.

8. How many genes in total were differentially expressed following MATQ-seq? Were lowly expressed genes removed? Was there any requirement that a gene was found to be expressed in at least N samples etc? These kind of filters are standard for scRNA-seq analysis and would seem to be appropriate for this type of data.

9. Using multiple pairwise comparisons is a restrictive way of analyzing this data – the authors should consider re-analyzing this data using an interaction test – trying to identify genes that respond differently to KRAS induction dependent on whether ARID1A is knocked out or not. Or consider an approach based on clustering of expression across differentially expressed genes to identify patterns of response across all four conditions. One prediction of the model the authors are proposing should be that SASP genes are upregulated in the KRAS vs WT but blunted in KRAS-ARID1A vs WT – is this observed?

10. Figure 3F: Figure 3F – please add a title containing the gene name – having the y-axis as ABM is confusing at first glance. Any supplementary table/data to show that ALDH1A1 was not in the list of DEG in AKC mice? Or to show that ALDH3A1 was? Is there a difference in regulation between species (i.e. both ALDH1A1 and ALDH3A1 are expressed, but only ALDH1A1 responds) or that only ALDH1A1 is expressed?

11. There is no comprehensive analysis linking the results from the in vitro MATQ-seq in mouse and the in vitro work in human. There is no discussion of whether orthologous genes are seen to up/downregulated in both. Linking together the results from the AKC vs KC lesions with KRAS-ARID1A-KO vs KRAS-Wildtype HPNE cells would allow the results from both models to be compared – which would strengthen the results presented.

12. Figure 4A: What is the expression of ALDH family members in the normal pancreas? This is really needed to make a comparison and statement as presented by the authors. In addition, Figure 4A is really showing high expression of ALDH1A1 and ALDH2 and not all "ALDH family proteins".

13. Figure 4H: The authors have demonstrated no evidence that in PanINs KRAS is a target of ARID1A? Was this observed in their analysis of HPFE cells? The origin of this link is not apparent.

14. ATAC-seq analysis: there is no analysis of the gained or lost regions at all. It is not stated in the text as to the numbers of gained/lost nor is there any analysis of the genes they are near (i.e. using GREAT or something similar). Are the sites that have significant changes in accessibility enriched for specific TFBSs etc.… i.e. other papers have suggested that ARID1A is a co-factor for AP-1 etc.…. (https://clinicalepigeneticsjournal.biomedcentral.com/articles/10.1186/s13148-019-0690-5).

15. Figure 5: Figure 5A: the text uses the phrase that Spearmans correlation coefficients are significantly lower.….. this is implying a hypothesis test – please rephrase. Figure 5B/C: needs to be clarified on how it shows "changes in accessibility of regulatory regions were significantly correlated with the alterations in gene expression levels in ARID1A-KO cells". How many of the 311 gained promoter ATAC regions overlap with genes which are upregulated (and vice versa)?

16. Supp. Table S4: There are 12 peaks reported in this table of which 10 are significantly "gained" this is different to the numbers reported in the manuscript. Please clarify/correct.

17. Figure S10: What is significance of identifying TFs? What are the significance of these TFs? Have they appeared in previous publications? Are any of these changing expression in ARID1A-KO RNA-seq? This is an underdeveloped analysis of what looks like good and interesting data.

18. Not all of the results from the differential expression analyses are available as Supplemental tables – this should be fixed.

19. In several parts of the manuscript – significance testing was carried out on only two data points (there can be no reliable estimate of variance using only two data points) using parametric methods – this is likely giving a false impression of significance – please either increase N or use a more appropriate test.

https://doi.org/10.7554/eLife.64204.sa1

Author response

Essential revisions:

1. The conclusion that "a highly complex trans-differentiation occurs following the knockout of Arid1a" is not supported by the evidence provided in Figure S3. It is unclear what panels B and C are referring to (no labels) and there is no scale given.

After more careful thinking, we agree with the reviewer’s critique. It is difficult to prove the transdifferentiation processes in vivo. In the revision, we deleted the related figure and descriptions.

The authors also use the RNA-sequencing data to conclude that "the activities of Kras signaling are partially suppressed by ARID1A deficiency." This is a major point and should be supported by experimental evidence in their cell models.

To provide the experimental support as suggested, we performed the western blot of ERK on ARID1A-KO and wild type HPNE cells and found that the activation of ERK signaling by KRAS was impaired in ARID1A-KO HPNE cells. We have added the new result to Figure 3—figure supplement 2.

2. Although Figure S5B shows that genetic manipulation of ARID1A was successful, there is no evidence by western blot or PCR that ARID1A has truly been knocked out in their HPNE cells. This model was used for a large portion of the data and must be characterized in a convincing way.

Following the reviewer’s suggestion, we have provided the qRT-PCR result (Figure 2—figure supplement 1C) to demonstrate that the expression of ARID1A is indeed significantly reduced in our HPNE cells with ARID1A knockout.

3. The evidence for senescence in vitro is not strong. Several times throughout, a colony formation assay is used as a measure of senescence, rather than SA-b-gal, molecular markers or secreted factors. It is unclear why the authors measured SA-b-gal staining intensity rather than %-positive cells.

First, we would like to point out that we have performed SA-b-gal staining on in vitro cells. However, we observed that the difference between ARID1A-KO and WT cells was not as significant as in vivo (Figure 2—figure supplement 1D). We have included this staining result in our first version. We attributed this result to the potential inefficiency in using SA-b-gal staining to identify the senescent HPNE.

The reason that we chose to measure SA-b-gal staining intensity in ex vivo experiment (Figure 2C-D) is as following: we found out that it was difficult to count the SA-b-gal negative cells due to the flat morphology of senescent cells and the large cell-size variations. Without accurate counting of negative cells, we cannot achieve the accurate calculation of the percentage of positive cells. We have added additional explanation in the revised manuscript to clarify the reason. To further clarify, for in vivo measurement, we still use %-positive analysis.

In particular, the evidence that ARD1A regulates Kras-induced senescence in HPNE cells is lacking.

To address the criticism for lack of evidence, we added the following experimental data in the revised manuscript: we examined the expression levels of CDKN1A in HPNE cell line model since CDKN1A is classical marker of cellular senescence. In the RNA-seq data, we observed that the expression of CDKN1A upon KRAS induction is also significantly suppressed in ARID1A-KO cells in comparison to wild type cells (Figure 3—figure supplement 3C). This result was further verified by qRT-PCR, as shown in Figure 3—figure supplement 3D. This new result provides additional experimental evidence supporting that ARID1A regulates KRAS-induced senescence in HPNE cells.

4. It is unclear if the ROS phenotype is mediated by ALDH1A1. Therefore, the conclusion that "ARID1A knockout can effectively reduce cellular ROS level by upregulating the expression of ALDH1A1" is not supported.

To provide the direct support for the role of ALDH1A1 in reducing cellular ROS, we added the ALDH1A1 knockdown experiment in the revision. The knockdown efficiency is shown by qRTPCR in Figure 4—figure supplement 2. We measured Kras-induced ROS level in ARID1A knockout HPNE cells with or without ALDH1A1 knockdown. As a result, ROS level in ALDH1A1 knockdown cells was significantly higher than in the cells with scramble shRNA upon Kras induction, indicating that ALDH1A1 is responsible for the clearance of ROS in HPNE cells with ARID1A knockout. This new result has been added to Figure 4H.

5. There is no functional experimental evidence provided for how ARID1A regulates the expression of ALDH1A1. This is a major point of novelty in the paper and should therefore be supported by experimental evidence.

In our first version of manuscript, we showed that ARID1A knockout could significantly increase the chromatin accessibility at the enhancer regions of ALDH1A1 (Figure 6A) in our ATAC-seq analysis. By analyzing the published ChIP-seq datasets of different TFs, we have identified the candidate factors that are involved in the regulation of ALDH1A1 expression. In particular, EP300 and NR3C1 could bind to two enhancer loci of ALDH1A1 (Figure 6C, left panel).

Here, to verify the predicted TFs based on ATAC-seq data as well as address this criticism for lack of experimental evidence on the role of ARID1A in regulating the expression of ALDH1A1, we have knocked down EP300 and NR3C1, respectively, in ARID1A-KO HPNE cells (Figure 6figure supplement 1). As shown in Figure 6C (middle and right panels), knockdown of EP300 or NR3C1 significantly impairs the transcription of ALDH1A1 gene. As a control, we also showed that the expression of EP300 and NR3C1 are not changed between ARID1A-KO cells and wild type cells (Supplementary File 3). Therefore, these new experimental data support the following regulatory mechanism: ARID1A deficiency promotes the accessibility of the enhancer regions of ALDH1A1, which leads to the binding of EP300 and NR3C1 to the ALDH1A1 enhancers and the activation of ALDH1A1 expression.

6. In the patients with pancreatic ductal adenocarcinoma, typically one sees point mutations, InDel or structure variations in the cancer genes; can the authors list the mutations in ARID1A in their samples to see if there are any recurrent mutations (i.e, visualized by lollipop plot, using ICGC-PACA-AU or ICGC-PACA-CA whole genome sequencing dataset, https://dcc.icgc.org/)?

To address this question, we have analyzed the mutation data from two cohorts (ICGC-PACA-AU and TCGA) and plotted the ‘Lollipop’ plots by using the MutationMapper web tool at http://www.cbioportal.org/. First, we didn’t find the recurrent mutations in the coding regions of ARID1A in these PDAC patients (Figure 1—figure supplement 1). It is worth noting that most of these mutations generate truncated proteins, which further support the roles of loss-of-function mutations in ARID1A.

Can the authors perform an eQTL analysis (for example using FastQTL) to see if any mutations in ARID1A and +/- 1Mb cis-mutations are associated with the expression of ALDH1A1? The mutations and gene expression datasets could be readily downloaded from ICGC-PACA-AU or ICGC-PACA-CA (ICGC data portal).

To test the eQTL based analysis, we have summarized the somatic mutations identified in ICGCPACA-AU and ICGC-PACA-CA datasets. In total, we identified 1510 somatic SNVs within +/- 1Mb regions of ARID1A. However, only 16 of them are recurrent mutations detected in multiple donors (15 mutations detected in two donors and 1 mutation detected in three donors). In the eQTL analysis, the SNVs in query need to be detected in multiple donors to reduce person-to-person variation and warrant statistical power. Here, considering the low number of recurrent mutations based on the two PDAC somatic SNV datasets, we conclude that we do not have statistical power to identify the mutations that are associated with the expression of ALDH1A1. Furthermore, we also examined whether there is significant enrichment of variants around the potential enhancer elements of ARID1A gene based on ATAC-seq data. However, we did not find any significant clustering of variants, which prevents us from performing meaningful eQTL analysis.

7. The authors describe differences in proportions of distinct populations of ADM, PanIN-1, PanIN-2, PanIN-3 populations. At no point in their introduction/results do they describe the PanIN-1/2/3 populations. This would be a good addition for clarity. In addition there are no details on how these were quantified in the methods.

As suggested, we added the details on lesion quantification in the method section in Supplementary Materials (Section PanIN quantification, Page 4).

In the introduction the authors state that previously ARID1A KO has been found to increased ADM – however Figure S1 shows the opposite that ARID1A KO leads to decreased ADM – while this reflects a trajectory of PanIN progression – not including any background on this is confusing and needs clarification.

We agree with the reviewer’s assessment. We have added additional description in the Main text (last paragraph, Page 4). Here, we would like to emphasize that in the previous study by Livshits et al. (Reference 8), the authors observed increase of ADM in adult mice two weeks after Arid1a was knocked out. In our study, we examined the effects of Arid1a knockout on PanIN progression at much later stage, from 2 months to 6 months after Arid1a knockout. Indeed, as the reviewer pointed out, this reflects a trajectory of PanIN progression.

8. How many genes in total were differentially expressed following MATQ-seq? Were lowly expressed genes removed? Was there any requirement that a gene was found to be expressed in at least N samples etc? These kind of filters are standard for scRNA-seq analysis and would seem to be appropriate for this type of data.

To address this critique, we have added the number of the differentially expressed genes into our manuscript (Supplementary File S2). The reviewer is also right. As a standard procedure, we do have a hard cutoff to remove lowly expressed genes when we performed DEG analysis on MATQseq. We only retained the genes with APM > 2 in at least five samples for DEG analysis. We apologize for missing some of the analysis information in our first version and we have added the details in the revised manuscript.

9. Using multiple pairwise comparisons is a restrictive way of analyzing this data – the authors should consider re-analyzing this data using an interaction test – trying to identify genes that respond differently to KRAS induction dependent on whether ARID1A is knocked out or not. Or consider an approach based on clustering of expression across differentially expressed genes to identify patterns of response across all four conditions. One prediction of the model the authors are proposing should be that SASP genes are upregulated in the KRAS vs WT but blunted in KRAS-ARID1A vs WT – is this observed?

We thank the reviewer for this suggestion. We have performed an interaction test by using edgeR to identify the genes that respond differentially to KRAS induction depending on the mutation status of ARID1A and we have added the results in the revised manuscript. To summarize, with the interaction test, we identified 156 genes whose responses to KRAS activation are significantly different between ARID1A-KO and wild type HPNE cells (FDR < 0.05) (Supplementary File S4). Next, we ranked the genes based on their differences in the response to KRAS induction (measured by the fold change in the output of edgeR DEG test) between ARID1A-KO and WT cells. Gene set enrichment analysis (GSEA) was then performed on the ranked gene list (Figure 3—figure supplement 3A-B). As a result, we observed that the activation of 4 inflammation-related pathways upon Kras induction (including interferon-α response, interferon-γ response, inflammatory response and TNF-α signaling via NF-κB pathways) are significantly suppressed in ARID1A-KO HPNE cells compared to WT cells, which suggests that these pathways are blunted in ARID1A-KO cells upon KRAS induction.

10. Figure 3F: Figure 3F – please add a title containing the gene name – having the y-axis as ABM is confusing at first glance.

To address this comment, we have added the gene name to Figure 3F. We have also modified the y-axis label as “gene expression levels (APM)” and the definition of the unit APM is given in the caption of the figure.

Any supplementary table/data to show that ALDH1A1 was not in the list of DEG in AKC mice? Or to show that ALDH3A1 was?

Yes, we included the full list of DEG table of our PanIN-seq data in Supplementary File S2 to show that Aldh3a1, but not Aldh1a1, is significantly upregulated in AKC mice.

Is there a difference in regulation between species (i.e. both ALDH1A1 and ALDH3A1 are expressed, but only ALDH1A1 responds) or that only ALDH1A1 is expressed?

Yes, there are species-based differences. ALDH3A1 is not expressed in HPNE cell line, but it is expressed in mouse model, vice versa for ALDH1A1 in HPNE cell line. Therefore, we do not have a scenario that both ALDH1A1 and ALDH3A1 are expressed, but only ALDH1A1 responds. But the species difference still indicates the regulation could be different. In our regulation analysis with ATAC-seq data and experimental verification, we mainly focus on the regulation of ALDH1A1 for human samples.

11. There is no comprehensive analysis linking the results from the in vitro MATQ-seq in mouse and the in vitro work in human. There is no discussion of whether orthologous genes are seen to up/downregulated in both. Linking together the results from the AKC vs KC lesions with KRAS-ARID1A-KO vs KRAS-Wildtype HPNE cells would allow the results from both models to be compared – which would strengthen the results presented.

To address this critique, we have compared the results from the in vivo mouse model to the in vitro human cell line model. We only find a limited number of genes that are significantly downregulated or upregulated in both (10 genes, the expected number of shared genes is 5.3). We speculate that the reasons for the limited number of orthologous genes sharing similar changes include cell type difference (the cell that we sequenced by PanIN-seq exist in the intermediate state of cell identity changes, while the HPNE has the fully established identity as epithelia cells) and the species difference.

Despite the limited number of orthologous genes sharing similar changes, we found that the perturbed pathways are quite similar between these two models. Among the 11 pathways significantly downregulated (FDR < 0.05) in ARID1A-KO cells under KRAS induction, 8 of them are also significantly downregulated (FDR < 0.05) in AKC lesions compared to KC lesions (the expected number is 3.74) (the shared pathways are labeled with asterisk in Figure 3—figure supplement 1). Here we didn’t compare the upregulated pathways in both models since the number of pathways that are upregulated in AKC lesions is small (only two pathways with FDR <0.1).

Overall, this new analysis suggests that while there are gene expression differences between these two model systems (mouse versus human, in vivo versus in vitro, and the neoplasia versus epithelia like cells), however, the general effects of ARID1A deficiency on the lesion cells are consistent across species as well as between in vivo and in vitro systems. We have clarified these points in the revised manuscript and added the corresponding analysis in Figure 3—figure supplement 1.

12. Figure 4A: What is the expression of ALDH family members in the normal pancreas? This is really needed to make a comparison and statement as presented by the authors.

To address this question, we have investigated the expression of ALDH family members in the normal pancreas and found that ALDH1A1 is the main subtype expressed in the normal pancreas (Figure 4—figure supplement 1A). With the published single cell RNA-seq data, we also observed that the expression of ALDH1A1 varies across different cell types, and the endocrine cells have the highest expression level, followed by exocrine glandular cells (Figure 4—figure supplement 1B).

As suggested by the reviewers, we also compared the expression of ALDH1A1 in PDAC to normal pancreas. Since the tumor cells are mainly epithelial cells, we only compared PDAC data to pancreatic ductal cells to avoid the confounding factors caused by the cell type difference. As shown in Figure 4—figure supplement 1B, there are four subclusters of ductal cells. The average expression level of ALDH1A1 in normal pancreatic ductal cells (cluster 1,2 and 3) is less than 50. we exclude the cluster 4 since ALDH1A1 positive cells in ductal cells correspond to the stem cell population. In PDAC patient cohort, we observed that 63% of samples have the expression level of ALDH1A1 larger than 50 TPM and 10% of samples have ALDH1A1 expression higher than 200 TPM (Figure 4—figure supplement 1C). These observations strongly suggest that ALDH1A1 plays important roles in a significant portion of PDAC tumors.

In addition, Figure 4A is really showing high expression of ALDH1A1 and ALDH2 and not all "ALDH family proteins".

We have modified the main text based on the reviewers’ suggestions to make the description more accurate.

13. Figure 4H: The authors have demonstrated no evidence that in PanINs KRAS is a target of ARID1A? Was this observed in their analysis of HPFE cells? The origin of this link is not apparent.

Based on our PanIN-seq data, we concluded that the activities of KRAS signaling pathway are clearly suppressed in AKC lesions compared to KC lesions. In term of studying the regulation, due to the limited number of cells in PanIN lesion, we were not able to perform ATAC-seq for the lesion cells, which limits our conclusion on whether KRAS is a direct target of ARID1A. On the other hand, since we also observed the similar effect on KRAS pathways upon ARID1A knockout in the cell line model, we can use our cell line data to investigate the potential regulatory relationship between ARID1A and KRAS signaling pathway.

To do so, we analyzed the changes of chromatin accessibility the genes in KRAS signaling pathway (the genes in Hallmark_KRAS_Signaling_UP gene set from MSigDB) by using the ATAC-seq of our HPNE cell data. As a result, we found that the fold-changes of chromatin accessibility of the genes involved in KRAS signaling pathway are significantly lower than that of randomly selected genes (Figure 5—figure supplement 5). This result suggests that ARID1A engages in the regulation of the genes involved in KRAS signaling pathway by modifying the chromatin accessibility of their regulatory regions. We have modified the description in our revised manuscript to explain this potential regulation mechanism (Page 13, first paragraph, Main text).

14. ATAC-seq analysis: there is no analysis of the gained or lost regions at all. It is not stated in the text as to the numbers of gained/lost nor is there any analysis of the genes they are near (i.e. using GREAT or something similar).

Since we focused on the regulatory function of ARID1A on the expression of ALDH1A1 in our manuscript, so we did not include the detailed analysis of ATAC-seq data in our first version. To address this critique, we also added the result of the global analysis of our ATAC-seq data using GREAT. Consistent with Figure 5C-D, the peaks with significantly changed accessibility between ARID1A-KO and WT cells are enriched at the distal regions (these new results are added to Figure 5—figure supplement 3A-B). The number of peaks with differential accessibility are all presented in the new figures. We also added the result of the GO enrichment analysis in GREAT in Figure 5—figure supplement 3C-D.

Are the sites that have significant changes in accessibility enriched for specific TFBSs etc.… i.e. other papers have suggested that ARID1A is a co-factor for AP-1 etc.…. (https://clinicalepigeneticsjournal.biomedcentral.com/articles/10.1186/s13148-019-0690-5).

To address this comment, we performed the motif enrichment analysis by using AME function in MEME. Briefly, we separate the peaks with significant changes into four types: distal peaks with increased accessibility, distal peaks with decreased accessibility, promoter peaks with increased accessibility and promoter peaks with decreased accessibility. We found the peaks with increased accessibility (including both promoter peaks and distal peaks) are enriched with FOX family TFs, SRY, CDX1, SOX5, TEAD4 etc. (Figure 5F, Figure 5—figure supplement 5A and Supplementary File S7). On the other hand, the peaks with decreased accessibility are enriched with FOS/JUN TF family, NFE2, NF2L1 etc., (Figure 5G, Figure 5—figure supplement 5B and Supplementary File S7). These results are consistent with previous reports as mentioned by the reviewers. We also added the suggested reference in the revised manuscript.

15. Figure 5: Figure 5A: the text uses the phrase that Spearmans correlation coefficients are significantly lower.….. this is implying a hypothesis test – please rephrase.

As suggested, we have rephrased our description. Furthermore, we provide the statistical test for Spearman’s correlation coefficients in Figure 5—figure supplement 1G.

Figure 5B/C: needs to be clarified on how it shows "changes in accessibility of regulatory regions were significantly correlated with the alterations in gene expression levels in ARID1A-KO cells".

To address this comment, we have added more detailed descriptions about the procedures in the revised manuscript (Page 12, the second paragraph). In details, we separated the DEGs identified in RNA-seq between ARID1A-KO and WT cells into two groups based on the fold-change (Figure 5B). The genes with positive log-fold change values are the genes upregulated in ARID1A-KO cells (Up-regulated genes group in Figure 5B) and the genes with negative log-foldchange values are the ones downregulated in ARID1A-KO cells (down-regulated genes group in Figure 5B). Next, we calculated the fold-change of chromatin accessibility for each of these DEGs and plotted the distribution of fold-changes of chromatin accessibility for these two groups of genes. As a result, we observed that the average fold-change of chromatin accessibility of upregulated genes is significantly higher than that of down-regulated genes (Figure 5B), which supports our statement.

How many of the 311 gained promoter ATAC regions overlap with genes which are upregulated (and vice versa)?

We performed the suggested analysis to examine the association between differentially expressed genes and the peaks with differential accessibility for both promoters and enhancers. As shown in Figure 5—figure supplement 4A-D, we observed that DEGs are significantly associated with the promoter peaks and the distal peaks with significantly differential accessibility. Furthermore, we noticed that the number of differentially expressed genes associated with the enhancer elements is larger than the number of genes associated with the promoters. This result is consistent with the observation that ARID1A knockout alters gene expression mainly by modulating the chromatin accessibility of the enhancers elements. A new paragraph describing this result (Page 12, last paragraph, Main text) is added in the revised manuscript.

16. Supp. Table S4: There are 12 peaks reported in this table of which 10 are significantly "gained" this is different to the numbers reported in the manuscript. Please clarify/correct.

Our description is probably not very clear. In the new Supplementary File S6, we listed all the ATAC-seq peaks that are linked to ALDH1A1 including both distal peaks and promoter peaks. In the Main text, we only described the peaks in the distal regions, therefore there are 11 peaks. We have added one sentence to clarify this point.

17. Figure S10: What is significance of identifying TFs? What are the significance of these TFs? Have they appeared in previous publications? Are any of these changing expression in ARID1A-KO RNA-seq? This is an underdeveloped analysis of what looks like good and interesting data.

Based on the literature search, we did not find previous reports showing that NR3C1 and EP300 are involved in the regulation of ALDH1A1 in cellular senescence. Here, we unveil how ARID1A mediate the gene expression of ALDH1A1 by these two genes for the first time, which shows the significance of identifying these TFs involved in ARID1A regulation.

There are no expression changes for both of the two genes between wildtype and ARID1A KO. To deepen our study, we also knocked down EP300 and NR3C1, respectively, in ARID1A-KO HPNE cells (Figure 6—figure supplement 1). As shown in Figure 6C (middle and right panels), knockdown of EP300 or NR3C1 indeed significantly impairs the transcription of ALDH1A1 gene.

18. Not all of the results from the differential expression analyses are available as Supplemental tables – this should be fixed.

In our first version of manuscript, we only included the genes that reach the statistical significance (FDR < 0.05). To address this critique, we have now provided the full table of DEG results in the original manuscript in Supplementary File S2.

19. In several parts of the manuscript – significance testing was carried out on only two data points (there can be no reliable estimate of variance using only two data points) using parametric methods – this is likely giving a false impression of significance – please either increase N or use a more appropriate test.

We apologize for our description that may cause the confusion about our experimental data and the statistics. For the figures that the reviewer mentioned with two data points, including Figure 2F, Figure 4D and Figure 4G, the two data points corresponds to two clones (clone #2 and clone #11), but for each clone, we did perform 3 biological replicates. We have clarified the description on the corresponding statistical tests.

https://doi.org/10.7554/eLife.64204.sa2

Article and author information

Author details

  1. Shou Liu

    Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, United States
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, validation, visualization, Writing – original draft, Writing – review and editing
    Contributed equally with
    Wenjian Cao and Yichi Niu
    Competing interests
    No competing interests declared

  2. Wenjian Cao

    1. Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, United States
    2. Genetics and Genomics Graduate Program, Baylor College of Medicine, Houston, United States
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing
    Contributed equally with
    Shou Liu and Yichi Niu
    Competing interests
    No competing interests declared

  3. Yichi Niu

    1. Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, United States
    2. Genetics and Genomics Graduate Program, Baylor College of Medicine, Houston, United States
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing
    Contributed equally with
    Shou Liu and Wenjian Cao
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4376-7792

  4. Jiayi Luo

    1. Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, United States
    2. Cancer and Cell, Biology Graduate Program, Baylor College of Medicine, Houston, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9873-0671

  5. Yanhua Zhao

    Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, United States
    Contribution
    Project administration
    Competing interests
    No competing interests declared

  6. Zhiying Hu

    Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, United States
    Contribution
    Project administration
    Competing interests
    No competing interests declared

  7. Chenghang Zong

    1. Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, United States
    2. Genetics and Genomics Graduate Program, Baylor College of Medicine, Houston, United States
    3. Cancer and Cell, Biology Graduate Program, Baylor College of Medicine, Houston, United States
    4. Dan L Duncan Comprehensive Cancer Center, Baylor College of Medicine, Houston, United States
    5. McNair Medical Institute, Baylor College of Medicine, Houston, United States
    Contribution
    Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – original draft, Writing – review and editing
    For correspondence
    czong@bcm.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8337-8038

Funding

NIH Office of the Director (1DP2EB020399)

  • Chenghang Zong

Robert and Janice McNair Foundation (McNair Scholarship)

  • Chenghang Zong

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are grateful to the McNair family for their support. CZ was also supported by the NIH Director’s New Innovator Award (1DP2EB020399). We thank Dr. Zhong Wang (University of Michigan) for providing the Arid1a floxed mice, Dr. Jennifer Bailey (UT Health Center) for providing the HPNE cell line, and Dr. Haoqiang Ying (MD Anderson Cancer Center) for providing inducible KRAS plasmids. We thank Kuanwei Sheng, Ejune Chen, Jinxiang Yuan, and other Zong lab members for their help in this project. We thank Dr. Christophe Herman for his proofreading and helpful discussion.

Ethics

This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols (#AN-6434) of Baylor College Medicine. Every effort was made to minimize suffering.

Senior and Reviewing Editor

  1. Maureen E Murphy, The Wistar Institute, United States

Version history

  1. Received: October 20, 2020
  2. Accepted: May 12, 2021
  3. Accepted Manuscript published: May 13, 2021 (version 1)
  4. Version of Record published: June 14, 2021 (version 2)

Copyright

© 2021, Liu et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,796
    Page views
  • 269
    Downloads
  • 6
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Shou Liu
  2. Wenjian Cao
  3. Yichi Niu
  4. Jiayi Luo
  5. Yanhua Zhao
  6. Zhiying Hu
  7. Chenghang Zong
(2021)
Single-PanIN-seq unveils that ARID1A deficiency promotes pancreatic tumorigenesis by attenuating KRAS-induced senescence
eLife 10:e64204.
https://doi.org/10.7554/eLife.64204

Categories and tags

Research organisms

Further reading

    1. Cancer Biology

    ONC201/TIC10 enhances durability of mTOR inhibitor everolimus in metastatic ER+ breast cancer

    Elena Farmaki, Aritro Nath ... Andrea H Bild
    Research Article

    The mTOR inhibitor, everolimus, is an important clinical management component of metastatic ER+ breast cancer (BC). However, most patients develop resistance and progress on therapy, highlighting the need to discover strategies that increase mTOR inhibitor effectiveness. We developed ER+ BC cell lines, sensitive or resistant to everolimus, and discovered that combination treatment of ONC201/TIC10 with everolimus inhibited cell growth in 2D/3D in vitro studies. We confirmed increased therapeutic response in primary patient cells progressing on everolimus, supporting clinical relevance. We show that ONC201/TIC10 mechanism in metastatic ER+ BC cells involves oxidative phosphorylation inhibition and stress response activation. Transcriptomic analysis in everolimus resistant breast patient tumors and mitochondrial functional assays in resistant cell lines demonstrated increased mitochondrial respiration dependency, contributing to ONC201/TIC10 sensitivity. We propose that ONC201/TIC10 and modulation of mitochondrial function may provide an effective add-on therapy strategy for patients with metastatic ER+ BCs resistant to mTOR inhibitors.

    1. Cancer Biology
    2. Computational and Systems Biology

    High-content microscopy reveals a morphological signature of bortezomib resistance

    Megan E Kelley, Adi Y Berman ... Gregory P Way
    Research Article

    Drug resistance is a challenge in anticancer therapy. In many cases, cancers can be resistant to the drug prior to exposure, i.e., possess intrinsic drug resistance. However, we lack target-independent methods to anticipate resistance in cancer cell lines or characterize intrinsic drug resistance without a priori knowledge of its cause. We hypothesized that cell morphology could provide an unbiased readout of drug resistance. To test this hypothesis, we used HCT116 cells, a mismatch repair-deficient cancer cell line, to isolate clones that were resistant or sensitive to bortezomib, a well-characterized proteasome inhibitor and anticancer drug to which many cancer cells possess intrinsic resistance. We then expanded these clones and measured high-dimensional single-cell morphology profiles using Cell Painting, a high-content microscopy assay. Our imaging- and computation-based profiling pipeline identified morphological features that differed between resistant and sensitive cells. We used these features to generate a morphological signature of bortezomib resistance. We then employed this morphological signature to analyze a set of HCT116 clones (five resistant and five sensitive) that had not been included in the signature training dataset, and correctly predicted sensitivity to bortezomib in seven cases, in the absence of drug treatment. This signature predicted bortezomib resistance better than resistance to other drugs targeting the ubiquitin-proteasome system. Our results establish a proof-of-concept framework for the unbiased analysis of drug resistance using high-content microscopy of cancer cells, in the absence of drug treatment.

    1. Biochemistry and Chemical Biology
    2. Cancer Biology

    FOXP2 confers oncogenic effects in prostate cancer

    Xiaoquan Zhu, Chao Chen ... Yanyang Zhao
    Research Article Updated

    Identification oncogenes is fundamental to revealing the molecular basis of cancer. Here, we found that FOXP2 is overexpressed in human prostate cancer cells and prostate tumors, but its expression is absent in normal prostate epithelial cells and low in benign prostatic hyperplasia. FOXP2 is a FOX transcription factor family member and tightly associated with vocal development. To date, little is known regarding the link of FOXP2 to prostate cancer. We observed that high FOXP2 expression and frequent amplification are significantly associated with high Gleason score. Ectopic expression of FOXP2 induces malignant transformation of mouse NIH3T3 fibroblasts and human prostate epithelial cell RWPE-1. Conversely, FOXP2 knockdown suppresses the proliferation of prostate cancer cells. Transgenic overexpression of FOXP2 in the mouse prostate causes prostatic intraepithelial neoplasia. Overexpression of FOXP2 aberrantly activates oncogenic MET signaling and inhibition of MET signaling effectively reverts the FOXP2-induced oncogenic phenotype. CUT&Tag assay identified FOXP2-binding sites located in MET and its associated gene HGF. Additionally, the novel recurrent FOXP2-CPED1 fusion identified in prostate tumors results in high expression of truncated FOXP2, which exhibit a similar capacity for malignant transformation. Together, our data indicate that FOXP2 is involved in tumorigenicity of prostate.