RNA-seq data collection of primary intestinal epithelial organoids undergoing ER stress or nutrient starvation
In order to gain a global understanding of how IECs are regulated by stress, RNA-seq analysis was first conducted on mRNA extracted from small intestinal organoid cultures from six mice either left unstimulated or treated with thapsigargin, a molecule that induces ER stress (Fig. 1a). In a second set of experiments, intestinal organoids from six other mice were used to analyze by RNA-seq the impact of nutrient starvation on IEC gene expression programs. For both stresses, we chose to analyze an early response (4 h stimulation), to concentrate on immediate cellular responses to these stresses and limit the impact of potential regulatory feedback loops that could alter gene expression programs in these organoids in a paracrine or autocrine-dependent manner. Moreover, analyzing early responses to these stresses limited the potential impact of cell death, which was undetectable under our experimental conditions at this early time of stimulation (data not shown).
Gene expression analysis was performed using the RPKM (Read Per Kilobase of transcript per Million) normalization method. A total of 21,553 genes were analyzed per condition, of which approx. 16,000 had detectable (>0) expression in organoids in the experimental conditions analyzed (Fig. 1b and Additional file 1: Table S1). In agreement with the role of the intestinal epithelium in nutrient absorption and host defense, the genes most expressed in our organoid cultures were associated with metabolism (including Aldoa, Aldob, Mt1, Mt2, Apoa1, Ldha) and innate immunity (Gm15284, Defa24, Lyz1) (Additional file 2: Table S2). Global analysis showed that, for each experimental condition, the triplicate cultures of organoids displayed overall very robust consistency in gene expression, with average correlation coefficients (standard deviation between triplicates divided by the mean for each gene analyzed) ranging between 0.338 to 0.455 (Fig. 1b). The coefficient of determination (R2) for both datasets (CTR_1 versus Thapsi_1 and CTR_2 versus Starv_2) was close to 1 (0.9956 and 0.9982, respectively), thus indicating that little experimental noise was generated in our assays and that most genes were not extensively altered at the transcript level by the stimulations (Fig. 1b-d). Nevertheless, 651 and 245 genes were found to be upregulated (>2.5x) and down-regulated (<0.4x) in thapsigargin-stimulated organoids, respectively, while 240 and 165 genes were found to be upregulated (>2.5x) and down-regulated (<0.4x) in nutrient starved organoids, respectively (Fig. 1e and Additional file 2: Table S2). On average, genes expressed in intestinal organoids were upregulated 1.35x and 1.10x by thapsigargin and nutrient starvation, respectively (Fig. 1f).
Transcript level analysis of murine intestinal epithelial organoids undergoing ER stress
We first analyzed the transcriptional response to thapsigargin in primary intestinal organoids. As expected, the transcriptional landscape shaped by thapsigargin treatment displayed a very strong (p = 1.6 × 10-12) ER stress-associated signature, and Gene Ontology (GO) analysis revealed that GO #0034976 (“Response to endoplasmic reticulum stress”) was the most significantly associated GO group among genes upregulated >2.5 fold by thapsigargin. ER stress related genes such as Atf3, Chac1, Thbs1, Ddit3, Derl3, Herpud1 and Hspa5 were upregulated more than 2.5 fold (Fig. 2a). Furthermore, the ER stress and UPR-associated transcription factor Xbp1, which is alternatively spliced upon ER stress, clearly underwent splicing to generate the exon 4 Δ26bp ER stress-specific isoform (Fig. 2b), and its expression was upregulated 2.5 fold. Thus, thapsigargin stimulation resulted in a strong induction of ER stress-associated responses in intestinal organoids.
Analysis of the 651 genes induced by thapsigargin over 2.5 fold revealed that, in addition to the core ER stress reprogramming at the transcript level, genes associated with cellular signaling (Group 1), gene expression regulation (Group 4), metabolism (Group 5) and the cytoskeleton (Group 3) were upregulated, as well as genes encoding channels and transporters (Group 2), cell surface molecules, receptors and extracellular matrix proteins (Group 6) and secreted factors (Group 7) (Fig. 2c). Therefore, ER stress induction by thapsigargin triggers a global transcriptional adaptation in intestinal organoids that affects multiple cellular processes and pathways.
Interestingly, a number of inflammation-associated genes were upregulated and were found mainly in Groups 6-7. These include the chemokines Cxcl10, Cxcl5, Cxcl1 and Ccl20, the cytokines Il23a and Tnf, the acute-phase response factor Saa3 and other genes encoding secreted mediators associated with inflammation, such as Fgf21 and Lcn2, as well as genes associated with the cellular adaptation to an inflammatory milieu (Hmox1) or innate immunity (Ifitm1, Nos2, Mx2, Defb40) (Fig. 2d). Thapsigargin stimulation also upregulated the expression of Fut1, a gene responsible for the fucosylation of extracellular matrix glycoproteins, an event that has recently been shown to play a major protective role in the response to bacterial pathogens in the small intestine [15]. Thus, in the intestinal epithelium, the cellular response to the ER stress inducer thapsigargin bears striking similarities to the inflammatory response against microbial infection.
Thapsigargin stimulation also upregulated genes associated with chemotaxis, migration and locomotion as well as genes encoding growth factors or growth factors receptors (Fig. 2d), which were found in Groups 6-7. Such growth factors included Hgf, Egf, Ctgf, Gdf15, Manf, Vegfa, Vegfc and Vgf. This might be of physiological importance given the known contribution of ER stress in the modulation of cancer microenvironment and tumor growth [16].
Global transcriptional reprogramming commonly requires upregulation of a core group of transcription factors (TFs) as part of the first induction wave of so-called “immediate early genes”. In agreement, genes associated with Gene Expression Regulation (Group 4) represented the largest group of genes up-regulated by thapsigargin (N = 145) and among those genes, a remarkable number (N = 46) of TFs were identified (Fig. 2d). Moreover, TFs were also found among the genes that are down-regulated by thapsigargin, including Mycn and Myb as well as the Notch target Hes5. Interestingly, Notch-dependent genes and in particular Hes family members play key roles in intestinal cell fate decision and proliferation, suggesting a potential role of ER stress signaling in IEC lineage commitment through Hes5 down-regulation (see also below). Together, these immediate-early TFs likely play a critical role in shaping the global cellular transcriptional adaptation to ER stress.
Thapsigargin stimulation also down-regulated the expression of a number of genes in intestinal organoids (N = 245 for genes down-regulated >2.5 fold = regulated <0.4 fold) (Fig. 2e). Of particular significance, we noticed that genes involved in metabolism (Group 5) were over-represented when compared to upregulated genes, suggesting that thapsigargin stimulation results in potent down-regulation of major cellular metabolic pathways (Fig. 2f). Moreover, among the genes associated with cellular signaling (Group 1), we noticed that a number encode for proteins associated with cell cycle regulation (Fig. 2f), suggesting a direct impact of thapsigargin on cellular proliferation (see also below, Fig. 4). Together, this analysis identified the core transcriptional program regulated by ER stress in primary intestinal organoids.
Transcript level analysis of murine intestinal epithelial organoids undergoing nutrient starvation
Primary murine organoids require a complex mixture of growth factors and nutrients to proliferate and to undergo differentiation of their stem cells into the different absorptive or secretory lineages of a functional intestinal epithelium. We chose to perform a short-term nutrient starvation by replacing the normal organoid culture medium with a Krebs-Ringer bicarbonate (KRB) buffer that lacks growth factors and AAs, which should result in acute inhibition of mTOR signaling and induction of the GCN2-dependent arm of the ISR. Under these conditions, although starvation resulted in the regulation of substantially fewer genes than thapsigargin (see above Fig. 1e), we noticed that the relative proportions of upregulated genes from Groups 1 to 7 were comparable between the thapsigargin and nutrient starved organoids (Additional file 3: Figure S1A). Similar to thapsigargin-treated organoids, a number of TFs, genes associated with inflammation, chemotaxis and growth factor signaling were identified in nutrient starved organoids (Additional file 3: Figure S1B). GO #0006935 (chemotaxis) was the most significantly associated GO term for genes upregulated >2.5 fold by nutrient starvation. With regards to genes down-regulated by nutrient starvation, 41.8 % (69/165) were genes whose function was uncharacterized (Additional file 3: Figure S1C), which is much greater than the percentage of uncharacterized thapsigargin-repressed genes (22.9 %) (Additional file 3: Figure S1D), suggesting that a significant portion of the cellular response to nutrient starvation in intestinal organoids involves poorly characterized pathways and processes. Nevertheless, among the genes with an attributed function, we noticed that much fewer were associated with cell cycle regulation as compared to genes down-regulated by thapsigargin (Additional file 3: Figure S1E), suggesting that ER stress and nutrient starvation affect cell cycle regulation in a distinct manner (see also Fig. 4 below). Finally, we observed that a number of genes associated with cellular innate immunity (Card11, Ifitm6, Mx2, Nlrp10, Nlrp1b, Nos3 and Noxa1) were downregulated by nutrient starvation, which contrasts with results obtained in thapsigargin-stimulated cells.
Together, we provide a comprehensive analysis of the modulation of the transcriptional landscape by nutrient starvation and reveal global similarity/dissimilarity patterns as compared to the regulation occurring in cells undergoing ER stress.
Identification of a common transcriptional signature for ER stress and nutrient starvation
We aimed to take our analysis one step further and to identify the common transcriptional signature associated with these two ISR-triggering cellular stresses. Focusing first on upregulated genes, we found that a core group of 90 genes were upregulated by both stresses (Fig. 3a). Considering the pool sizes (N = 651 and N = 240 for thapsigargin stimulation and nutrient starvation, respectively), this intersection was found highly significant (p = 1.44×10-101), thus showing that a strong common transcriptional signature exists for genes upregulated by these two stresses. It must be noted that these common genes represented a smaller fraction of the thapsigargin upregulated pool (90/651 = 13.8 %) than the nutrient starvation upregulated pool (90/240 = 37.5 %), thus showing that in our experimental conditions, thapsigargin induced a wider spectrum of stress-specific genes than nutrient starvation. Consistent with this, a large majority of the genes induced by thapsigargin were not significantly modulated by nutrient starvation (Fig. 3b).
When upregulated genes were separated in functional categories, it appeared that most of the genes upregulated by both stresses were found in Groups 4 (“Gene expression”), 6 (“Cell surface molecules, receptors and extracellular matrix proteins”), 7 (“Secreted”) and 1 (“Cellular signaling”) (Fig. 3c-d). Among the 90 common genes, those that were most upregulated by both conditions included TFs (found in Group 4) and factors associated with inflammation and/or chemotaxis (in Groups 6/7) (Fig. 3c and e).
The analysis of down-regulated genes revealed that few genes (N = 18) were found to be down-regulated by both stresses (Additional file 4: Figure S2A-B), and these genes were either of uncertain function (7/18) or had unrelated functions (Additional file 4: Figure S2C-D).
Together, these results demonstrate the existence of a common transcriptional signature associated with two distinct ISR-inducing stresses. While a core group of 11 TFs likely promotes this transcriptional reprogramming, genes associated with inflammation and chemotaxis were unexpectedly also part of this common transcriptional signature.
Effect of cellular stress on intestinal stem cells and proliferation
In primary intestinal organoids, cell proliferation is dependent on the activity of ISCs and ISC-derived transit-amplifying cells that express the marker Lgr5 and occupy an important fraction of the intestinal crypt [17]. The above results suggested that ER stress and nutrient starvation affected the expression of cell cycle and/or proliferation genes differentially in intestinal organoids (see Fig. 2 and Additional file 3: Figure S1). Therefore, to identify the overall impact of thapsigargin versus nutrient starvation on ISCs, we analyzed how these stresses affected the transcript levels of the 151 ISC-enriched genes identified recently by sorting Lgr5+ intestinal epithelial cells [18]. Interestingly, the majority of ISC genes displayed reduced expression following thapsigargin treatment (fold induction <1), while the opposite result was observed in nutrient starved cells (Fig. 4a). Overall, expression of ISC genes was not upregulated (1.06 fold) by thapsigargin while this treatment upregulated gene expression 1.34 fold genome-wide (Figs. 1f and 4b). In contrast, ISC genes were significantly more upregulated (1.41 fold) than all genes (1.1 fold) by nutrient starvation (Figs. 1f and 4b). Despite these differences, we observed that a group of seven genes were similarly modulated by both stresses: Tnfrsf19, Wwtr1, Vav3, Esrrg and Notch1 were down-regulated by thapsigargin and nutrient starvation, while Car12 and Rasa3 were upregulated. These seven genes might thus represent the core group of genes involved in the cellular adaptation of ISCs to stress.
Finally, we aimed to directly test the differential effect of thapsigargin and nutrient starvation on cellular proliferation in intestinal organoids by measuring EdU incorporation in organoids undergoing stress for 4 h. In agreement with the ISC gene expression data, we observed that thapsigargin stimulation resulted in potent reduction of EdU incorporation in organoids, indicative of a block in proliferation, while nutrient starvation had the opposite effect (Fig. 4d-e). Thus, while ER stress affects ISC gene expression and proliferation in a way that is expected for a cellular stress and is in line with previous results [6], short-term nutrient starvation appeared to trigger an apparently paradoxical pro-proliferative response.
Analysis of the AS landscape in enteroids undergoing ER stress and nutrient deprivation
One key aim of this study was to gain insight into the role of AS on gene regulation during metabolic stress – an idea that has been understudied in mammalian systems. In order to investigate global changes in AS occurring upon nutrient starvation and ER stress, AS events were analyzed from the RNA-seq dataset using criteria that has previously been described [13]. Applying a stringent cutoff (p ≥ 0.90, DBS = 0.1, RN = 5, See
Methods
) to identify high confidence AS events, we identified a group of 85 AS events that were significantly upregulated by thapsigargin treatment and 42 AS events induced by nutrient starvation (Fig. 5a, Additional file 5: Figure S3A and Additional file 6: Table S3). When classifying the AS events by type of splicing event, the largest group during ER stress corresponded to exon skipping (S/I) events, followed by an approximately 25 % proportion being alternative 5’ or 3’ donor/acceptor changes, with complex events (C3, C2, IR-C) representing the smallest percentage. During nutrient starvation however, we observed a more even distribution among S/I events and events involving full intron retention (IR-S), as IR-S and complex intron retention (IR-C) events comprised more than half of the observed AS events. To rule out that the increase in AS events could be attributed to simply an upregulation in transcription, we analyzed the expression levels of the genes found to undergo AS upon ER stress or nutrient starvation. Only a handful of genes were found to be upregulated at the transcript level (Casp4, Clk4, Cdkn2aip, Gpcpd1, Alg12, Taf1a, and 1110021L09Rik for ER stress and Glipr1 and Maff for nutrient starvation), while the rest of the genes underwent no significant change in expression (Fig. 5b, Additional file 5: Figure S3B).
Interestingly, closer examination of the thapsigargin-regulated AS events displaying the highest confidence score in our analysis (p ≥ 0.95, DBS = 0.15, RN = 10) revealed that the vast majority of those events were predicted to produce a reading frame shifting caused by exon inclusion, skipping or alternative 5’/3’ splice site usage (Fig. 5c). Upon nutrient starvation conditions however, we observed that the proportion of AS events that resulted in frame shifted reading frames (29 out of 42 events = 69 %) was roughly that which would be expected by chance (2/3 = 66.6 %) (Additional file 5: Figure S3C). The predicted stress-induced AS events that met our cutoff criteria were visualized using Integrated Genome Viewer (IGV) (Fig. 5d and Additional file 5: Figure S3D).
Gene ontology analysis of the most enriched gene family among the genes undergoing AS during both ER stress and nutrient deprivation revealed a significant enrichment for genes involved in mRNA splicing (p = 1.4x10-8), including the U snRNA maturation-associated factor Lsm7, the SR protein Srsf7 and the Smn paralog Smndc1 (Fig. 5e and Additional file 5: Figure S3E), providing further support for our results that identified a key role for metabolic stress stimuli in the dynamic regulation of the splicing machinery [19].
Increased expression of PTC-containing AS isoforms of RNA processing/splicing genes upon metabolic stress
As mentioned above, GO analysis of the genes undergoing both thapsigargin- and nutrient deprivation-upregulated AS revealed a significant enrichment in genes associated with mRNA splicing regulation. Using the cut-off criteria outlined earlier, we observed a significant overlap (p = 9.6 × 10-8) of the AS events in splicing factors during nutrient starvation and ER stress and found five genes (Hnrnpd, Ogt, Srsf7, Ivns1abp and Hnrpdl) to be spliced identically during either stress (Fig. 6A). These events were analyzed by compared percent spliced-in (PSI) values, as well as Sashimi plots generated via IGV, and were all validated by RT-PCR analysis RNA splice variant specific primer pairs (Fig. 6B and primer list in Table 1). Interestingly, all of these events involved the partial or full retention of intronic sequences that resulted in the introduction of an in-frame premature stop codon, producing a transcript to be degraded by nonsense mediated decay (NMD) according to the 55-bp rule [20]. Furthermore, these NMD-inducing events are evolutionarily conserved, as identified by other groups [9].
Splicing of detained introns functions to regulate gene expression during cellular stress
A recent study identified a pool of mRNAs containing retained introns (referred to as detained introns or DIs in this study) that were spliced post-transcriptionally under stress conditions [21]. Interestingly, splicing of these DIs often resulted in a frameshifting event that generated potential NMD splice variants. Strikingly, several DI-containing genes reported in human cells by Boutz et al. were orthologs of genes found in our list (Eny2, Ogt, Srsf7, Rbm39, Zfp326, Tra2a, Cdc16, Ivns1abp, Tex10, Hnrpdl, Slc35b1, Ufd1l and Clk4), suggesting that ER stress might represent a physiological means to regulate splicing of DIs. Notably, all five of the AS events mentioned above occurring upon both ER stress and nutrient deprivation represented DIs identified by Boutz et al. to flank premature termination codon (PTC)-containing cassette sequences. Furthermore, extended manual analysis using IGV software revealed that this trend was consistent in RNA-binding proteins (Rbm39), U1 snRNP related (Snrnp70), U2 snRNP related (Sf3b1), core SR proteins (Srsf3, Srsf6, Srsf5), and other SR proteins (Srrm1, Tra2a), among others (Additional file 7: Figure S4). In agreement with the hypothesis that thapsigargin-regulated AS events corresponded to DIs, several AS events identified were predicted to generate potential NMD variants.
Boutz et al. identified that the SR protein kinase CLK4 and its homologue CLK1 possess PTC-containing DIs (introns 3 and 4) that are spliced out to generate a functional translated protein (Fig. 7a). This provides a means of autoregulation and serves as the upstream regulator of splicing of the downstream DIs, as treatment with the CLK inhibitor CB19 affected the splicing of >300 DIs [21]. Consistent with the previous observations that the “fully spliced” (DI-lacking) CLK1/4 transcripts encode a stable mRNA transcript, we observed an increase in expression of both CLK4 and CLK1 upon thapsigargin treatment (Fig. 7b). Furthermore, we observed an increase in CLK4 protein upon both thapsigargin and nutrient starvation (Fig. 7c). Interestingly, we observed significant splicing of introns 3 and 4 in both CLK1 and CLK4 upon thapsigargin treatment, analyzed by both PSI values and IGV Sashimi plots (Fig. 7d-e). Notably, there was no significant increase in CLK1/4 mRNA expression, nor did we observe any change in PSI levels of intron 3 or 4 upon nutrient starvation, suggesting the presence of another upstream mechanism to regulate DI splicing upon starvation.
Nutrient starvation and ER stress induce a subset of distinct AS events
In addition to a common signature of AS induced in enteroids during nutrient deprivation and ER stress, we noticed that most AS events (80/85 and 37/42 for ER stress and nutrient starvation, respectively) were stress-specific. Notably, the vast majority of these AS events that were unique to ER stress or nutrient starvation occurred in genes involved in diverse biosynthetic pathways. For example, AS events induced by stress occurred in various families of enzymes, such as cysteine-type peptidases (ie. Usp15, Usp40), mannosyltransferases (ie. Alg9, Alg12), isomerases (ie. Rpe, Trub2) and serine/threonine kinases (ie. Bco2, Pla2g2a). Interestingly, 3 of the identified AS targets in thapsigargin-stimulated cells are factors previously associated with ER stress response, including Casp4, Slc35b1 and Ufd1l [22–24] (Fig. 8a, Additional file 8: Figure S5A). Furthermore, among the genes observed to undergo AS upon nutrient deprivation, several were involved enzymatic processes. Examples of these genes include the ferric-chelate reductase Frrs1 and the NAD(P) transhydrogenase Nnt, which both undergo increased exon skipping upon starvation, and the 5’-nucleotidase Nt5c3, a gene for which we observed increased retention of a portion of intronic sequence (Fig. 8b, Additional file 8: Figure S5B). These results suggest that AS of genes involved in the response to stress contributes to a regulatory program associated with the cellular response to ER stress and nutrient deprivation.
In silico prediction of the use of alternative translation start sites to generate novel isoforms
We noticed that many AS events induced by thapsigargin affected exons towards the 5’ end of the mRNAs, and database search (Aceview) identified the existence of potential alternative ATG (aATG) as potential translation initiation sites for 12 genes (Casp4, Zfp326, Acot8, Ivns1abp, Tra2a, Rpl27a, G630016D24Rik, Slc35b1, Tex10, Rbm39, U2af1, BC024814) among the 32 genes with the highest confidence AS (p ≥ 0.95, DBS = 0.15, RN = 10, data not shown), which could support the expression of shorter variants with truncated N-terminal ends (Additional file 9: Figure S6). All of these genes have been reported by Wilson et al. to utilize a novel aATG in both mouse and humans to maintain the proper frame of reading. Interestingly, many of the isoforms that would be generated upon usage of an alternative aATG have prominent changes in either protein domains (Fig. 8c) or protein localization (Fig. 8d). For example, upon skipping of the Casp4 98 bp exon 3, a novel variant using a previously annotated aATG just after the splicing event would generate an isoform lacking the entire caspase activation and recruitment domain (CARD domain). Splicing of the solute carrier Slc35b1 in a manner that generates a 25 bp alternative 3’ event induced by thapsigargin and the 5’-nucleotidase Nt5c3, which undergoes an intron retention event of 55 bp upon nutrient starvation, resulted in the loss of 1 or more transmembrane domains and likely affects the subcellular localization of these proteins. Lastly, we predict that the thapsigargin-induced Alt3’ splicing of the ER stress gene Ufd1l and Alt5’ splicing of C1 family peptidase family gene Tinag would also utilize previously annotated aATGs to generate protein isoforms that are re-localized from the mitochondria to the cytosol, and the extracellular space to the nucleus, respectively [25].