RNA-seq evaluation. High-quality reads were aligned to mouse reference genome (NCBIM37) using HISAT2 software with higher accuracy and faster running speed. Between 96%-98% of reads were aligned to the reference genome, about 88%-91% of reads aligned to the reference genome in the unique method, all multiple mapped ratios in different three group <10%, the difference between total mapped ratio with unique mapped ratio<10% (Table.1), suggesting that the sequencing quality and mapping results all meet the requirements, also DNA or ribosome contamination is eliminated. Besides, we can also intuitively illustrate that there is no problem in the process of RNA extraction and cDNA library construction through statistical analysis of genome structure distribution and density distribution of sequences on chromosome. Statistics show that the percentage of sequencing sequences located in the exon region is the highest, while the location of introns and intergenic regions is relatively low (Fig 1A), thereby eliminating the contamination of immature mRNA in the three samples. Located the sequence to each chromosome on the genome and plot it, statistics reveal that there are no abnormal coverage positions in chromosomes (Fig 1B).
mRNA structural variation. SNP analysis results showed that there are 37896 and 42592 SNP sites in Non-transfection and GFP-vector respectively, the ART1 knockdown cells possessed 29418 SNP sites. The incidence of insertion or deletion (Indel) of small fragments at a certain position in the genome in three group were much lower than that of SNP (Fig 1C). Among this SNP sites, the number of transition mutations(A>G, C>T, G>A, T>C) was far larger than that of transversion (A>C, A>T, C>A, C>G, G>C, G>T, T>A, T>G). Also, in this detection, the majority of SNP sites occurred in CDS, intronic, ncRNA and 3’UTR regions, and less happened in splice, intergenic and 5’UTR regions (Fig 1D).
RNA-Seq and ASprofile software together can be used to predict different alternative splicing patterns and splicing sites. Among the ten possible splice forms, skipped exon (SKIP), alternative 5’first exon (transcription start site, TSS) and alternative 3’last exon (transcription terminal site, TTS) had the highest incidence (Fig 1F). Researches have demonstrated that many types of factors including polypyrimidine tract binding protein (PTB), CELF family (CUG-BP, Elav-like family), CALD1, can participate in the regulation of alternative splicing, and splicing variants of variable splicing factors which affect tumor cell growth and angiogenesis[10-12]. Interestingly, in this study, the prediction of possible alternative splicing events found that the knockdown of ART1 affected the occurrence of the splicing of SLC39a14, VEGF, Met, and the number of alternative splicing events in multiple alternative splicing factors (Mbnl1, Ceflf1, Cald1 and Ptbp3) has also changed significantly. (Supplement Table1&2.)
Identification of the DEGs. TPM is currently recognized as the most accurate transcript quantification formula. The number of genes uniquely expressed in each group and the number of genes shared between three groups are shown in a Venn diagram (TPM>0). (Fig 2A). In the subsequent DEGs and enrichment analysis, GFP-vector was selected as the control group. DEGseq(version 1.12.4), an R package to identify differentially expressed genes or subtypes from RNA-seq data in different samples. The DEGs were screened with criterions: Q-Value<0.05 and |Fold Change|>2. In the GFP-shART1 vs. GFP-vector, 1052 genes highly expressed in GFP-shART1 group, and 4500 down-regulated genes were identified (Fig 2B). In addition, the top 10 up-and down-regulated genes in GFP-shART1 group were listed in Tables.2. It should be pointed out that after knocking down ART1, PTBP1, which plays an important role in pre-mRNA splicing and in the regulation of alternative splicing events was significantly down-regulated, log2FoldChange=-20.604763, qValue<0.01. This also suggested that ART1 has a certain regulatory effect on alternative splicing.
GO term enrichment analysis of DEGs between GFP-vector and GFP-shART1 groups. Enrichment analysis can be performed on the significance of GO functions of differential genes, which can help us to identify the biological pathways most relevant to biological phenomena in differential gene sets. It is considered that the function is enriched when the corrected Q-value <0.05. The GO enriched analysis indicated that the most enriched GO terms in DEGs between GFP-vector and GFP-shART1 were “intracellular” (GO: 0005622) and “intracellular part” (GO: 0044424), which belongs to the cellular component, followed by “intracellular organelle” (GO: 0043229) and “organelle” (GO: 0043226) (Q-value<0.05) (Fig 2C); It means that the active position of DEGs caused by ART1 knockdown are mainly located in intracellular, especially in organelles. For molecular function, DEGs were significantly enriched in “binding” (GO: 0005488), followed by “heterocyclic compound binding” (GO: 0005488) (Q-value<0.05) (Fig 2D), besides, more than half of the top 30 enriched molecular functions listed in the bubble chart are related to various forms of bindings; For biological process, the highest enrichment of DEGs was “cellular metabolic process” (GO: 0044237), followed by “metabolic process” (GO:0008152) (Q-value<0.05) (Fig 2E). All of the above are in line with our previous research that ART1 is abundantly expressed in the endoplasmic reticulum, may increase glycolysis and energy metabolism in CT26 cells under high glucose conditions by regulating the Akt/c-Myc signaling pathway and the expression of glycolytic enzymes[13].
KEGG enrichment analysis. The KEGG pathway enrichment analysis of differentially expressed genes was performed using Kobas software. The DEGs were enriched significantly in 33 KEGG pathways (Q-value<0.05), and the top 30 signal pathways with the highest enrichment were shown in the bubble chart. Ubiquitin mediated proteolysis(ko04120), Protein processing in endoplasmic reticulum(ko04141), Ribosome biogenesis in eukaryotes(ko03008), RNA transport(ko03013), Endocytosis(ko04144) were ranked in the top five (Fig 2F).
Among them, what arouses our great interest is the mitogen activated protein kinases (MAPK) signaling pathway with 87 differential genes enriched (Q-value: 0.0002466). The extracellular-signal-regulated kinases (ERK1/2), the jun N-terminal kinase or stress activated protein kinases (JUK or SAPK), p38/MAPK and ERK5 are four major subfamilies of MAPK signaling[14]. In particular, researchers have recognized that the main regulator of colorectal cancer cell proliferation is the ERK MAPK signaling pathway. Moreover, the mutation of the KRAS proto-oncogene is an early event in the development of CRC. The high activity of Ras is accompanied by the increase of ERK activity, which causes colorectal cancer cells to escape from normal growth and differentiation control and has the ability to invade surrounding tissues and organs[15]. Kras also hemizygously mutated at p.G12D in CT26 cells, however, KEGG enrichment map showed that the knockdown of ART1 in CT26 cells reduced the expression of MEK1(log2FoldChange=-1.85), ERK (log2FoldChange=-2.17) and its downstream c-myc (log2FoldChange=-1.66). Also, knocking down ART1 significantly down-regulated mitogen-activated protein kinase kinase 4 and 7(MKK4 and MKK7) (log2FoldChange-2.36= and -3.34) which can activate the JNK double phosphorylation site to prevent the activation of the JUK signaling pathway, in addition, the transcription level of JNK also decreased (Supplement Figure1.). All of the above indicated that ART1silencing can prevent the proliferation of colorectal cancer cells to a certain extent. It needs to be added that the previous research which published in 2017 also suggested that ART1 serves a faciliatory role in the proliferation and migration of CT-26 cells, and this effect may associate with the factors downstream of FAK and RhoA, c-myc, c-fos and COX2[16, 17], which is also consistent with the results obtained by this RNA sequencing. However,the potential biological mechanism of the down-regulation of various key signaling molecules caused by ART1 inhibition needs further research.
ART1 is a mono-ADP-ribosyltransferase that is stably anchored to the cell membrane or abundantly expressed in the cytoplasm and specifically catalyzes arginine residues. In our previous research, ART1 was also found to be specifically expressed in the endoplasmic reticulum, which is combined with the result that more differential expressed genes obtained in the GO classification and enrichment analysis can be enriched into the cellular component of organelle, let us have a strong interest in protein processing in endoplasmic reticulum pathway (ko04141) obtained from KEGG enrichment analysis (Fig 3A). Sequencing results showed that ART1 silencing significantly downregulated the expression of its substrate protein glucose-regulated protein 78(GRP78/BIP) (log2FoldChange=-1.98), also inhibited the expression of endoplasmic reticulum transmembrane protein (PKP like ER protein kinase(PERK) (log2FoldChange=-1.26), activating transcription factor 6(ATF6) (log2FoldChange=-2.81) and inositol-enquiring enzyme 1(IRE1α) (log2FoldChange=-4.6) in the three important unfolded protein reaction(UPR) signal transduction pathways, suggesting that ART1 is likely to play a key role in the regulation of unmatured protein flow and the relief of endoplasmic reticulum load in endoplasmic reticulum stress environment. In addition, we also selected the top 10 functional categories with the highest enrichment degree and the differential genes related to this function, drew the network diagram of the significant enrichment function and gene interaction (Figure3B).
ART1 inhibition induced the apoptosis of CT-26 cells by overall repression of UPR reaction. Endoplasmic reticulum stress promotes the processing of misfolded or unfolded proteins accumulated in the cavity of the endoplasmic reticulum, in this way, the normal function of the cell is better maintained, which is conducive to cell survival, and the subsequent response triggered by this is the unfolded protein response. However, the unfolded protein response is characterized by two-way regulation in the process of processing cell growth and death. On the one hand, it can maintain cell survival. On the other hand, cell apoptosis will occur immediately when severe endoplasmic reticulum stress exceeds the load of the cell. More importantly, recent studies have found that if the unfolded protein response is fully inhibited, the cell will lose the ability to maintain the homeostasis of the endoplasmic reticulum and die[18, 19]. GRP78/Bip plays an important role in this adjustment process. Under normal physiological conditions, the N-terminals of the three receptor proteins, PERK (PEK like ER kinase), IRE1α (inositol requiring 1), and ATF6α (activating transcription factor 6), combine with GRP78/Bip to form a complex, thereby being in an inactive state and inhibiting signal transduction. When stress occurs, misfolded or unfolded proteins will accumulate to cause the complex to disaggregate, release PERK, IRE1α, ATF6α, and activate related pathways to promote the degradation of misfolded proteins and the correct folding of unfolded proteins. Emerging evidence indicated that UPR signaling components may play important roles in cancer since the upregulation of GRP78, PERK, IRE1α and ATF6 were detected in various clinical and tumor cell samples. The inhibition of the signal pathway mediated by a single component will cause the activation of the compensation mechanism of other components. Therefore, it is necessary to find the key points that can directly and comprehensively destroy the three related signal pathways in the UPR response of tumor cells, which can inhibit the survival of tumor cells.
It has been established that ART1 is widely distributed in endoplasmic reticulum, and its role in the endoplasmic reticulum is mainly to modify GRP78 by mono-ADP ribosylation to inactivate it[20, 4]. Recent studies demonstrated that ART1 was transiently up-regulated during ER stress induction with the treatment of dithiothreitol or thapsigargin[21]. Also, compared with normal tissues, the significantly higher expression of ART1 in colorectal cancer has been clarified in previous experiments There, combined with the results obtained by RNA sequencing, we speculated that the increase of ART1 in colorectal cancer may be also associated to the presence of more endoplasmic reticulum stress stimuli in tumor cells.
Our previous study confirmed that ART1 was abundantly expressed in endoplasmic reticulum. In this study, we detected the activation of the signaling pathways mediated by the three interactor proteins of GRP78. PERK, IRE1α and ATF6 all decreased significantly after the knockdown of ART1 by lentiviral transfection, in addition, XBP-1s and p-eif2α were downregulated by ART1 inhibition, suggesting that inhibiting the expression of ART1 can completely block the activation of the three signal pathways. Besides, the expression of GRP78 also decreases after the UPR signaling pathways were blocked, which consistent with the result obtained in RNA-seq analysis (Fig 4A). At the same time, real-time PCR were performed to detect the mRNA level PERK, IRE1α, ATF6 and GRP78, the degree of PERK, IRE1α, ATF6 and GRP78 all lower in GFP-shART1 group (Fig 4B). Interestingly, it was found that the expression of Cleaved-caspase3 increased and the expression of Bcl2 decreased (Fig 4C), indicating that knocking down ART1 may hinder cell growth by completely inhibiting the UPR signaling pathway in the endoplasmic reticulum, and at the same time, it affects the level of Bcl2 and makes the cells tend to apoptosis. Through the results of this experiment, it is found that tumor cell apoptosis not only occurs when the endoplasmic reticulum stress is overloaded, but also when the UPR response, an important regulatory signaling pathway for endoplasmic reticulum stress, is fully inhibited.