Identification of the role of mono-ADP-ribosylation in colorectal cancer by integrated transcriptome analysis

Our previous study clarified the carcinogenic properties of arginine-specific mono-ADP-ribosyltransferase 1 (ART1), which results in a critical post-translational modification that changes the structure and function of proteins and is widely involved in important processes. This study provides, for the first time, a comprehensive transcriptomic analysis of colorectal cancer cells with ART1 silencing by Illumina RNA-Seq and related verification experiments. Lentiviral infection was used to construct a CT-26 cell line with stable knockdown of the ART1 gene, and a whole transcriptome sequencing technique was performed to identify differentially expressed genes (DEGs). GO and KEGG classification/enrichment analyses and verification experiments were performed to determine the role of ART1 in the progression of colorectal cancer. A total of 5552 DEGs, GO functions and KEGG pathways with the highest enrichment, various SNPs, and diverse splicing patterns were identified. Importantly, knockdown of ART1 affected the splicing of certain key genes related to tumor cell growth and downregulated the expression of the key gene PTBP1 for alternative splicing. The overall attenuation of the endoplasmic reticulum unfolded protein response (UPR) signaling pathway caused by the inhibition of mono-ADP-ribosylation of GRP78 could disrupt UPR signaling, leading to the occurrence of apoptosis to impede tumorigenesis. ART1, which is clustered in organelles, may promote the development of colorectal cancer by participating in a variety of new mechanisms, including endoplasmic reticulum stress regulation and alternative splicing, and may be a good clinical drug target for targeted therapy of CRC.


Introduction
Colorectal cancer (CRC) is the second leading cause of cancer-related deaths worldwide. In both sexes combined, CRC ranks third (6.1%) in terms of incidence and second (9.2%) in terms of mortality according to the status report on the global burden of cancer issued by the International Agency for Research on Cancer (IARC) based on GLOBO-CAN [1]. The 5-year survival rate is approximately less than 10% when the tumor is diagnosed at advanced stages. Targeted therapy is currently a vital treatment that can prolong the survival of patients with advanced stages. Unfortunately, there are still patients for whom no alternative targeted drugs are available. Additionally, adverse reactions to targeted drugs limit their long-term therapeutic effects, and some new breakthroughs are expected. Therefore, exploration of more pivotal and specific key molecules as therapeutic targets for this malignancy and elucidation of the molecular mechanisms of CRC by expanding basic research are urgently needed.
ADP-ribosylation is a common modification of macromolecules, such as proteins and nucleic acids, and widespread modification is of increasing interest because of its association with key biological functions and cellular processes, including transcription, DNA damage repair, stress responses, and cell signaling [2]. ADP-ribosylation is the enzymatic transfer of a single or multiple ADP-riboses from NAD + to acceptor molecules, termed MARylation and PARylation, respectively [3,4]. Compared to that on PARylation, 111 Page 2 of 14 research on the functions of MARylation is lacking, especially the identification of mono-ADP-ribosylated proteins. However, over time, MARylation has gradually become important due to an increasing amount of evidence showing that MARylation serves as a reversible post-translational modification that can be read by readers and removed by MAR-specific erasers. Furthermore, mono-ADP-ribosyltransferase-catalyzed ribosylation of protein substrates usually lead to protein inactivation, providing a mechanism to inhibit protein functions in both physiological and pathological conditions. ART1, which is widely found anchored in the cell membrane and localized in the cytoplasm, is the only enzyme reported to catalyze arginine-specific mono-ADP-ribosylation in humans and mice. In addition, in the past 5 years, our research group has conducted extensive research on mono-ADP-ribosylation and shown that ART1 and PARP-1 have a synergistic effect on CDDP-induced apoptosis of CT26 cells. The mechanism may be related to ART1 regulation of the RhoA/ROCK pathway, thereby affecting the expression of NF-kB and PARP-1 [5], which indicates that poly-ADP-ribosyltransferase is regulated by mono-ADP-ribosyltransferase in CRC [6]. Mono-ADP-ribosylation may be the basis of poly-ADP-ribosylation, which has an important role. In addition, roles of ART1 in the regulation of cell proliferation, apoptosis, tumor angiogenesis, and DNA repair have been described in our past research; among them, ART1 is a vital cancer-promoting factor that is significantly increased in colorectal cancer [7,8]. While ART1 is a dangerous carcinogen in colorectal cancer, we hope to find more differentially expressed genes through high-throughput RNA sequencing and collect functional enrichment and signal pathway enrichment information to explore the molecular mechanism of ART1 in the development of colorectal cancer, which will contribute to targeted therapy for colorectal cancer.
CT-26 cells share molecular features with undifferentiated and refractory human colorectal carcinoma cells [9]. Therefore, as a widely used clinical and basic research mouse tumor cell line, CT-26 is an optimal model for research on biological mechanisms or preclinical evaluation of targeted therapy and immunotherapy. Here, we used lentiviral transfection to establish a CT-26 cell line with stable knockdown of ART1 to carry out high-throughput sequencing and related enrichment analysis.

Cell culture and lentivirus cell line establishment
The mouse colorectal cancer cell line CT26 isolated from Balb/c mice were provided by Professor Wei YQ of Sichuan University. The construction and packaging of lentiviruses were performed by GeneChem Corporation (Shanghai, China). Establishment of nontransfected, GFP vectortransfected, and GFP-shART1-transfected cell lines was described previously [6]. Briefly, CT-26 cells were infected with lentivirus for 48 h and selected in 2 µg/ml puromycin in culture medium. All groups of CT26 cells were grown in the recommended RPMI-1640 medium (HyClone, Logan, UT, USA) supplemented with 10% fetal bovine serum and antibiotic solution (100 µg/ml streptomycin and 100 U/ml penicillin, HyClone) at 37 °C and an atmosphere of 5% CO 2 . The cells in the three groups were lysed with RNAiso Plus (TaKaRa Bio, Inc.) and sent to Sangon Biotech Co., Ltd.

RNA isolation
Total RNA was extracted using a Total RNA Extractor (TRIzol) kit (B511311, Sangon, China) according to the manufacturer's instructions. DNA contamination was eliminated by DNase I treatment.

RNA-Seq read mapping
Clean reads were processed and aligned to the mouse reference genome (NCBIM37) using HISAT2 (http:// ccb. jhu. edu/ softw are/ hisat2/ manual. shtml, version 2.0) with default parameters. Hierarchical indexing and multiple alignment strategies were used to solve the short spliced-read alignment problem. Then, RSeQC (http:// rseqc. sourc eforge. net/; version 2.6.1) was used to statistically analyze the alignment results.

Analysis of gene structure and chromosome distribution
The homogeneity distribution and the genome structure were assessed by Qualimap (http:// www. quali map. org; version 2.2.1). BEDTools (http:// samto ols. sourc eforge. net; version 2.26.0) was used to statistically analyze the gene coverage ratio.

Genetic structural analysis
In the latest version of BCFtools (http:// samto ols. sourc eforge. net; version 1.5), SNP calling can be achieved by BCFtools instead of Samtools and BCFtools at the same time. SNPs (single-nucleotide polymorphisms) and indels (insertions or deletions) in each sample were extracted by SNP/InDel calling, and Snp Eff was used to determine the distribution of mutation sites on the genome structure. Additionally, quality value > 20 and coverage cutoffs > 8 were used as filtering criteria in this process. In this study, ASprofile software (download form http:// ccb. jhu. edu/ softw are/ ASpro file/; version 2.36) was used to extract and quantify and compare alternative splicing (AS) events from RNA sequence data. We can directly compare multiple transcripts of the same gene through asport software to identify different alternative splicing events.

Differential gene expression (DEG) analysis
We measured the relative abundances of transcripts through standardized RNA-Seq fragment counts with StringTie (http:// ccb. jhu. edu/ softw are/ strin gtie; version 1.3.3b). The unit of measurement is transcripts per million (TPM), which eliminates the influence of gene lengths and sequencing discrepancies to enable direct comparison of gene expression between samples. The calculation formula of TPM was as follows: The Venn Diagram package in R language (Venn Diagram R package; https:// CRAN.R-proje ct. org/ packa ge= VennD iagram; version 1.6.17) was used to visually show the similarity and overlap of the number of expressed genes between samples. DESeq2 (R package; https:// github. com/ mikel ove/ DESeq2; version 1.12.4) is an R package and was used to estimate variance-mean-dependence in count data from RNA sequencing assays and determine differentially expressed genes (DEGs) between two samples based on a model using the negative binomial distribution. Genes were considered significantly differentially expressed if Q-value < 0.001 and |fold change| > 2. Gene expression differences were visualized by scatter plots, MA plots, and volcano plots.

Functional analysis of differentially expressed genes (DEGs)
Functional enrichment analysis, including Gene Ontology (GO) and KEGG analyses, was performed to identify which DEGs, proteins, or other molecules were significantly enriched in GO terms or pathways. DAVID v 6.7 was used to map DEGs to the GO terms in the database followed by calculation of the number of genes in every term, and a hypergeometric test was performed to identify significantly enriched GO terms in the gene list from the background of the reference gene list. Cluster Profiler was used for functional enrichment analysis. GO terms annotation was visualized by a histogram. GO enrichment analysis results are displayed using a bubble chart. KOBAS software was used for enrichment analysis of KEGG pathways of differentially X i = total exon fragment∕reads L i = exon length KB expressed genes. In general, GO terms and KEGG pathways with Q-values < 0.05 were considered significantly altered.

Western blot analysis
Total protein was extracted with RIPA lysis buffer and cocktail (1:100). Equal amounts of protein from each group were loaded onto sodium dodecyl sulfate (SDS) polyacrylamide gels at 80 V for 0.5 h and then 100 V for 1 h and finally transferred to polyvinylidene fluoride (PVDF) blotting membranes. The membrane was blocked with 5% milk or bovine serum albumin (BSA) for 1 h at room temperature and then incubated with the primary antibody overnight at 4 °C. All antibody concentrations were strictly in accordance with the instructions. The membrane was removed the next day, washed three times, and then incubated with the secondary antibody for 1 h at room temperature. Protein bands were visualized by luminal chemiluminescence (Bio-Rad, Hercules, CA, USA). The gray value was detected and analyzed by ImageJ software.

Real-time polymerase chain reaction analysis
Briefly, quantitative real-time PCR of PERK, IRE1α, ATF6, and GRP78 was performed using SYBR Green Master Mix (TaKaRa Bio, Inc.) on a 7500 system (Applied Biosystems, Foster City, CA). The primers are shown as follows:

RNA-seq evaluation
High-quality reads were aligned to the mouse reference genome (NCBIM37) using HISAT2 software, which has high accuracy and a fast running speed. Approximately 96-98% of the reads were aligned to the reference genome, approximately 88-91% of the reads were aligned to the reference genome in the unique method, all multiple mapped ratios in three different groups were < 10%, and the difference between the total mapped ratio and the unique mapped ratio was < 10% ( Table 1), suggesting that the sequencing quality and mapping results all met the requirements, and DNA or ribosome contamination was eliminated. In addition, we intuitively illustrated that there was no problem in the RNA extraction and cDNA library construction through statistical analysis of the genome structure distribution and density distribution of sequences on the chromosomes. Statistics showed that the percentage of sequencing sequences located in the exon region is the highest, while that of introns and intergenic regions are relatively low (Fig. 1a), indicating that the contamination of immature mRNA in the three samples was eliminated. Locating the sequence to each chromosome on the genome and plotting it, we revealed that there were no abnormal coverage positions in the chromosomes (Fig. 1b).

mRNA structural variation
SNP analysis results showed that there were 37,896 and 42,592 SNP sites in nontransfected and GFP vector-transfected cells, respectively, and the ART1 knockdown cells had 29,418 SNP sites. The incidence of indels of small fragments at a certain position in the genome in the three groups was much lower than that of SNPs (Fig. 1c). Among these SNP sites, the number of transition mutations (A>G, C>T, G>A, T>C) was much larger than that of transversions (A>C, A>T, C>A, C>G, G>C, G>T, T>A, T>G). Additionally, in this study, the majority of SNP sites occurred in CDSs, intronic sequences, ncRNAs, and 3′UTRs and fewer occurred in splice sites, intergenic regions, and 5′UTRs ( 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. 1e). Studies have demonstrated that many types of factors, including polypyrimidine tract-binding protein (PTB), CELF family members (CUG-BP, Elav-like family), and CALD1(Caldesmon 1), can participate in the regulation of alternative splicing and splicing variants of variable splicing factors that affect tumor cell growth and angiogenesis [10][11][12]. Interestingly, in this study, the prediction of possible alternative splicing events showed that the knockdown of ART1 affected the occurrence of SLC39a14, VEGF, and Met splicing, and the number of alternative splicing events in multiple alternative splicing factors (Mbnl1, Ceflf1, Cald1, and Ptbp3) also changed significantly (Supplement Tables 1 and 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 DEG and enrichment analysis, the GFP vector was selected as the control group. DEGseq (version 1.12.4) is an R package to identify differentially expressed genes or subtypes from RNA-Seq data in different samples. The DEGs were screened with the following criteria: Q-value < 0.05 and |fold change|> 2. In the GFP-shART1 vs. GFP vector comparison, 1052 genes were highly expressed in the GFP-shART1 group, and 4500 downregulated genes were identified (Fig. 2b). In addition, the top 10 up-and down-regulated genes in the GFP-shART1 group are listed in Table 2. Notably, after knockdown of ART1 expression, PTBP1, which plays an important role in pre-mRNA splicing and in the regulation of alternative splicing events, was significantly downregulated (log2 fold change = − 20.604763, Q-value < 0.01). This finding also suggested that ART1 has a certain regulatory effect on alternative splicing.

GO term enrichment analysis of DEGs between the GFP vector and GFP-shART1 groups
Enrichment analysis can be performed on differentially expressed genes to determine the significance of GO functions, which can help us to identify the biological pathways most relevant to the biological phenomena in differentially expressed gene sets. The function is considered to be enriched when the corrected Q-value was < 0.05. The GO enrichment analysis indicated that the top enriched GO terms for the DEGs between the GFP vector and GFP-shART1 groups were "intracellular" (GO: 0005622) and "intracellular part" (GO: 0044424), which belong to the cellular component category, followed by "intracellular organelle" (GO: 0043229) and "organelle" (GO: 0043226) (Q-value < 0.05) (Fig. 2c), indicating that the DEGs induced by ART1 knockdown were mainly located intracellularly, especially in organelles. For molecular function, the DEGs were significantly enriched in "binding" (GO: 0005488), followed by "heterocyclic compound binding" (GO: 0005488) (Q-value < 0.05) (Fig. 2d), and more than half of the top 30 enriched molecular functions listed in the bubble chart were related to various forms of binding. For biological process, the category with 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 findings are in line with our previous research showing that ART1 is abundantly expressed in the endoplasmic reticulum and 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 are shown in the bubble chart. Ubiquitinmediated proteolysis (ko04120), protein processing in the endoplasmic reticulum (ko04141), ribosome biogenesis in eukaryotes (ko03008), RNA transport (ko03013), and endocytosis (ko04144) were the top five pathways (Fig. 2f). We were interested in the mitogen-activated protein kinase (MAPK) signaling pathway, with enrichment of 87 differentially expressed genes (Q-value: 0.0002466). Extracellular signal-regulated kinases (ERK1/2), 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, mutation of the KRAS proto-oncogene is an early event in the development of CRC. The high activity of Ras is accompanied by an increase in ERK activity, which causes colorectal cancer cells to escape normal growth and differentiation control and invade surrounding tissues and organs [15]. Kras was also hemizygously mutated at p.G12D in CT26 cells; however, the KEGG enrichment map showed that the knockdown of ART1 expression in CT26 cells reduced the expression of MEK1 (log2 fold change = − 1.85), ERK (log2 fold change = − 2.17), and its downstream factor c-myc (l log2 fold change = − 1.66). Additionally, knocking down ART1 expression significantly downregulated mitogen-activated protein kinase kinase 4 and 7 (MKK4 and MKK7) (log2 fold change = − 2.36 and − 3.34), which can activate the JNK double-phosphorylation site to prevent the activation of the JNK signaling pathway. In addition, the transcription level of JNK decreased (Supplement Fig. 1). All of the above results indicated that ART1 silencing can prevent the proliferation of colorectal cancer cells to a certain extent. Previous research published in 2017 also suggested that ART1 promotes proliferation and migration of CT-26 cells, and this effect may be associated with factors downstream of FAK and RhoA, such as c-myc, c-fos, and COX2 [16,17], which is also consistent with the results obtained by RNA sequencing. However, the potential biological mechanism of the downregulation 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; this finding combined with the result that more differentially expressed genes obtained in the GO classification and enrichment analysis can be enriched into the cellular component of organelle indicated the importance of protein processing in the 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) (log2 fold change = − 1.98) and also inhibited the expression of the endoplasmic reticulum transmembrane protein PKP-like ER protein kinase (PERK) (log2 fold change = − 1.26), activating transcription factor 6 (ATF6) (log2 fold change = − 2.81), and inositol-enquiring enzyme 1 (IRE1α) (log2 fold change = − 4.6) in the three important unfolded protein response (UPR) signal transduction pathways. These results indicated that ART1 is likely to play a key role in the regulation of immature protein flow and the relief of endoplasmic reticulum load under endoplasmic reticulum stress conditions. In addition, we selected the top 10 functional categories with the highest enrichment degree and the DEGs related to this function and drew a network diagram of the significantly enriched functions and gene interactions (Fig. 3b).

ART1 inhibition induced the apoptosis of CT-26 cells by overall repression of the UPR
Endoplasmic reticulum stress promotes the processing of misfolded or unfolded proteins that accumulate in the cavity of the endoplasmic reticulum; thus, the normal function of the cell is maintained, which is conducive to cell survival, and the unfolded protein response is subsequently triggered. However, the unfolded protein response is characterized by two-way regulation of cell growth and death. On the one hand, this process can maintain cell survival. On the other hand, cell apoptosis will occur immediately when severe endoplasmic reticulum stress exceeds the capacity 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 process. Under normal physiological conditions, the N-termini of the three receptor proteins, PERK, IRE1α, and ATF6α, combine with GRP78/Bip to form a complex, which is in an inactive state and inhibits signal transduction. When stress occurs, misfolded or unfolded proteins accumulate to cause the complex to disaggregate; release PERK, IRE1α, and ATF6α; and activate related pathways to promote the degradation of misfolded proteins and the correct folding of unfolded proteins. Emerging evidence indicates that UPR signaling components may play important roles in cancer since the upregulation of GRP78, PERK, IRE1α, and ATF6 expression was detected in various clinical and tumor cell samples. The inhibition of the signaling pathway mediated by a single component will activate a compensatory mechanism of other components. Therefore, it is necessary to identify the key points that can directly and comprehensively destroy the three related signaling pathways in the UPR of tumor cells, which can inhibit the survival of tumor cells.
ART1 is widely distributed in the endoplasmic reticulum, and its role in the endoplasmic reticulum is mainly to modify GRP78 by mono-ADP-ribosylation, inactivating it [4,20]. Recent studies have demonstrated that ART1 was transiently upregulated during ER stress induction with dithiothreitol or thapsigargin treatment [21]. Additionally, compared with that in normal tissues, the significantly higher expression of ART1 in colorectal cancer has been clarified in previous experiments. Combined with the results obtained by RNA sequencing, we speculated that the increase in ART1 in colorectal cancer may also be associated with the presence of more endoplasmic reticulum stress stimuli in tumor cells.
Our previous study confirmed that ART1 was abundantly expressed in the 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 signaling pathways. In addition, the expression of GRP78 decreased after the UPR signaling pathways were blocked, which was consistent with the  (Fig. 4a). Moreover, real-time PCR was performed to detect the mRNA levels of PERK, IRE1α, ATF6, and GRP78, and the levels of PERK, IRE1α, ATF6, and GRP78 were all lower in the GFP-shART1 group than the control group (Fig. 4b). Interestingly, it was found that the expression of cleaved caspase 3 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; moreover, it affects the level of Bcl2 and promotes cell apoptosis. Through the results of this experiment, we found that tumor cell apoptosis occurs not only when endoplasmic reticulum stress is overloaded but also when the UPR, an important regulatory signaling pathway for endoplasmic reticulum stress, is fully inhibited.

Discussion
Our previous research found that ART1 plays an important role in the proliferation, apoptosis, adhesion, movement, metastasis, and angiogenesis of mouse colon cancer cells. First, ART1 plays a facilitatory role in the proliferation and migration of CT26 cells and the potential mechanism associated with the downstream factors FAK, RhoA, c-Myc, c-fos, and COX-2 [10]. Additionally, under high-glucose conditions, ART1 may increase the generation of ATP and lactic acid by upregulating the AKT/mTOR/c-Myc pathway, thus stimulating the proliferation and inhibiting the apoptosis of CT-26 cells by increasing the expression of PKD1 and LDHA [22]. Second, ART1 could influence the apoptosis rate through regulation of the mRNA and protein expression of the critical factor Tubb3 [23]. In addition, ART1 has a positive correlation with VEGF and integrin  αVβ3 expression and is associated with the angiogenesis of colorectal carcinoma or upregulates HIF-1α via the PI3K/ AKT signaling pathway, thus promoting the expression of angiogenic factors [8]. Therefore, ART1 is highly overexpressed in colorectal carcinoma tissues and can accelerate many aspects of the development of this disease. Moreover, research on targeted proteins modified by ART1 has received increased attention. Integrin α7, human neutrophil peptide (HNP-1), platelet-derived growth factor-BB (PDGF-BB), fibroblast growth factor-2 (FGF-2), and GRP78/Bip are receptor proteins directly modified by ART1 and are related to tumor immunity, angiogenesis, cell proliferation, etc. [21,[24][25][26]. More importantly, the above research results are roughly consistent with the results of this transcriptome analysis, and RhoA (log2 fold change = − 1.04), c-Myc (log2 fold change = − 1.66), COX-2 (log2 fold change = − 3.81), and PKD1 (log2 fold change = − 4.53) all decreased significantly. Therefore, we hope to find more molecular mechanisms to confirm the carcinogenic properties of ART1 in colorectal cancer. Moreover, since ART1 mainly plays a role by modifying its substrate protein, we hope to find a signaling pathway related to the substrate protein directly modified by ART1 that may contribute to the targeted or combined treatment of colorectal cancer. This study provides the first comprehensive insight into the effect of ART1 silencing on the transcriptome of CT-26 colorectal cancer cells. We selected CT-26 cells that share molecular features with aggressive, undifferentiated, refractory human colorectal carcinoma cells, strictly adhered to good quality control standards of sequencing data, and then performed accurate reference sequence alignment. Using a whole transcriptome sequencing technique (RNA-Seq), we were able to identify the levels of differentially expressed genes with TPM values and subsequently performed GO and KEGG classification and enrichment analysis. Moreover, the numbers and forms of SNPs and diverse splicing patterns were subjected to diversified statistical analysis. The above findings provide a basis for analysis of ART1 in colorectal cancer from multiple perspectives and to explore the potential specific molecular mechanism.
TPM was used as a standardized value to standardize the number of read counts in transcripts and genes. However, compared with the classic FPKM (fragments per kilobase million) and RPKM (trans per million) algorithms, TPM has a different processing order; that is, the gene length is first considered, the sequencing depth is considered, and the calculated value can be directly used for comparison between samples. Generally, this data standardization algorithm is superior to FPKM and RPKM [27]. A total of 1052 highly expressed genes and 4500 downregulated genes in the GFP-shART1 group were identified. Among the top ten genes listed in this study that increased and decreased after ART1 knockdown, sorted by log2|fold change|, ART1 was selected, and ART1 silencing can significantly increase the expression of the protein arginine transferase prmt1. Prmt1, which is mainly distributed in the cytoplasm and nucleus, is associated with a variety of cancers and participates in various cellular processes, such as cell signal transduction, DNA damage repair, and transcriptional regulation. In particular, studies have suggested that PRMT1 is positively correlated with the proliferation and migration of colorectal cancer cells and can inhibit the apoptosis of some colorectal cancer cells. Therefore, whether the high expression of prmt1 caused by ART1 silencing shows the same effect in colorectal cancer is worthy of further investigation. Moreover, ART1 silencing can strongly increase the expression of GAPDH, which provides new ideas for the selection of internal controls in our daily research.
GO functional annotation and KEGG pathway enrichment analysis were used to screen the biological functions, cell localization, and signaling pathways affected by ART1. The differentially expressed genes obtained by knocking down ART1 expression were enriched in 33 signaling pathways when the Q-value was less than 0.05 as the standard, and KEGG enrichment analysis screened 73 pathways when the P-value was < 0.05 as the standard. The MAPK signaling pathway (ko04010) and protein processing in the endoplasmic reticulum (ko04141) are the two signaling pathways we will focus on. Fabrizio et al. found ADP-ribosylated GRP78/Bip on arginine residues largely distributed on cellular plasma membrane, in the lumen of ER, which is the result of the modification of the enzyme ART1 [21]. Furthermore, the information from UniProt (https:// www. UniPr ot. org) and our previous results indicated that there is a large amount of ART1 distributed in the endoplasmic reticulum. Researchers have also provided evidence that GRP78/BiP undergoes ADP-ribosylation and that ADP-ribosylation is a rapid post-translational mechanism for reversible inactivation of GRP78/BiP through interfering with allosteric coupling Fig. 3 Processing in endoplasmic reticulum pathway maps of DEGs and KOG enrichment network for KEGG enrichment analysis. a All the highlighted gene products in the figure belong to the genes annotated in this transcriptome. The rectangular nodes represent gene products (such as enzymes or some RNA regulators), and the circular nodes represent compounds (that is, substrates or products). The white background circle and the rectangular rectangle indicate other pathways associated with this pathway. The colors and positive and negative indicate the relative up-and-down relationship of genes in the two comparison samples. Red indicates upregulated genes. Reddish colors indicate higher gene expression, green indicates downregulated genes, and green colors indicate higher gene expression. Yellow indicates nondifferentially expressed genes. The color block size represents the proportion of the corresponding gene in the gene product. b The square represents the functional classification enriched by DEGs, the circle represents the names of related genes, and the connecting lines indicate the correlation between functions or genes (Color figure online) of two domains of GRP78/BiP and may play an important role in fast regulation of unfolded protein load [28,29]. Then, more direct evidence indicated that ADP-ribosylated GRP78/BiP provides a buffering system that balances the protein processing rates with those of protein synthesis. All of the above results suggested that ART1 is abundantly expressed in the endoplasmic reticulum and modifies GRP78 to inactivate it, which prompted us to investigate this signaling pathway in this transcriptome analysis. Therefore, we speculated that ADP-ribosylated GRP78/BiP will affect its downstream UPR and the regulation of endoplasmic reticulum homeostasis. Consistent with the sequencing results, our Western blot and RT-PCR results showed that ART1 silencing significantly downregulated the expression of its substrate protein GRP78/ BiP and inhibited the expression of the endoplasmic reticulum transmembrane proteins PERK, ATF6, and IRE1α. The downstream key signaling molecules XBP-1 and p-eif2α were also significantly inhibited, indicating that ART1 inhibition can completely block the activation of the unfolded protein response signaling pathway in the endoplasmic reticulum. More importantly, and consistent with the sequencing results, Bcl-2 (log2 fold change = − 1.95) decreased significantly, and cleaved caspase 3 increased significantly after intervention in ART1, suggesting the occurrence of apoptosis. Since members of all classes of the Bcl-2 family localize to the ER membrane and have been confirmed to influence ER homeostasis, we speculate that the unfolded protein response is fully inhibited, which increases the release of Ca 2+ in the ER, inhibits the expression of Bcl-2, and induces the occurrence of apoptosis. However, the detailed molecular mechanism is still unclear, and it remains to be defined whether increased ADP-ribosylated GRP78/BiP can affect the release of GRP78/BiP interactors and the downstream signaling pathways that participate in this process.
Through analysis, we not only found a different number of variable shear events in each group but also found that ART1 knockdown affected many types of factors, including PTBP, the CELF family, and CALD1, which can participate in the regulation of alternative splicing. PTBP1, which binds to the pyrimidine-rich region in the RNA sequence, interferes with the expression of the corresponding protein by affecting the positioning of the mRNA and variable splicing and ultimately triggers a change in the biological effect in the cell [30,31]. In the sequencing analysis results of this study, PTBP1 was affected by ART1 knockdown and showed a significant decrease in expression. Therefore, how the alternative splicing molecule PTBP1 regulated by ART1 affects the development of colorectal cancer is still the subject of future research.
In this study, CT26 cells with ART1 intervention were successfully sequenced. A total of 5552 differentially expressed genes were found and annotated. Enrichment analysis showed that ART1 may be involved in the development of colorectal cancer by participating in a variety of new mechanisms, including endoplasmic reticulum stress regulation and variable shear regulation. All of the above findings provide us with new ideas for understanding the carcinogenic properties of arginine mono-ADP-ribosylation in colorectal cancer.
Author contributions SXZ and JLD processed the cell culture, RNA extraction, and analysis of the sequencing result. SXZ wrote the manuscript. YPY and HJG were responsible for the process of validation experiments. YT, MX designed and conducted the experiments. QSL, ML also participated in analysis and interpretation of data of RNA-Seq. YLW was responsible for funding acquisition and manuscript modification.