In a previous study, it was demonstrated that H3R117 is a mono-ADP-ribosylated site in Lovo cells, and Lovo cells with a point mutation resulting in H3R117A were successfully constructed [mutation of arginine at residue 117 of H3 to a non-ADP-ribosylated alanine (14)]. The untransfected Lovo cells (donated by Professor Wei-Xue Tang, Chongqing Medical University) and H3R117A Lovo cells were both lysed with RNAiso Plus (Takara Bio, Inc.) and sent to Sangon Biotech Co., Ltd. Both types of cells were cultured in DMEM (HyClone; GE Healthcare Life Sciences) with 10% FBS (HyClone; GE Healthcare Life Sciences) with 5% CO₂ at 37°C.

The raw image format obtained by Illumina Hiseq^™ (Illumina, Inc.) via Casava (http://support.illumina.com/sequencing/sequencing_software/casava.ilmn; version 1.8.2) base calling analysis was converted to raw sequenced reads, and termed raw data. The raw data quality value and other information were statistically analyzed, and the raw data quality of samples was visually evaluated using FastQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/; version 0.11.2). Trimmomatic (http://www.usadellab.org/cms/index.php?page=trimmomatic; version 0.36) was used to filter the raw data to obtain the clean data by removing n-bases and removing the adaptor sequence in reads: Forward, AGATCGGAAGAGCACACGTCTGAAC and reverse AGATCGGAAGAGCGTCGTGTAGGGA. Low quality bases were removed from reads 3′-5′ (q-value <20); low quality bases were removed from reads 5′-3′ (Q-value <20); the sliding window method (23) was used to remove the bases with a mass value <20 in the tails of reads (window size, 5 bp); the reads and their matched reads with a length of <35 nucleotides were removed; 10,000 sequences were randomly extracted from the clean data and compared with the NCBI NT database (http://ncbi.nlm.nih.gov/) using blastn. E-value ≤1×10⁻¹⁰ and comparison of results with a similarity >90% and coverage >80% were used to calculate species distribution and perform pollution detection.

HISAT2 (version 2.1.0) was used to align the effective sample data to the reference genome, and statistical comparisons were performed using RseQC (http://rseqc.sourceforge.net/; version 2.6.1). RSeQC was used to analyze the duplicate reads and the distribution of inserted fragments based on the alignment results. Qualimap (http://www.qualimap.org; version 2.2.1) was used to test the gene body coverage and analyze the structural distribution of the genome. BEDTools (http://code.google.com/p/bedtools; version 2.26.0) was used for statistical analysis of gene coverage and read distribution on chromosomes.

BCFtools (http://samtools.sourceforge.net; version 1.5) was used to perform single nucleotide polymorphism (SNP)/insertion-deletion (indel) calling, and the SNPs/indels in each sample were extracted. Subsequently, SnpEff (snpeff.sourceforge.net; version 2.36) was used to calculate the distribution of variation sites on the genome structure; a quality value >20 and coverage cutoffs >8 were used as the criteria for filtering. Alternative splicing analysis was used for fusion gene analysis using ASprofier (http://ccb.jhu.edu/software/ASprofile; version 2.36) and EricScript (http://ericscript.sourceforge.net; version 0.55).

In the RNA-seq analysis, the expression levels of genes were estimated by counting the reads located in the genomic regions or gene exon regions. In addition to the true expression levels of the gene, the read count was positively correlated with the length and sequencing depth of the gene. To make the estimated gene expression levels of different genes and different experiments comparable, the concept of transcripts per million (TPM) was introduced. The calculation formula of TPM was as follows:

StringTie (http://ccb.jhu.edu/software/stringtie; version 1.3.3b) and known gene models (24) were used to assess gene expression.

Biological replicates are necessary for any biological experiment. The correlation of gene expression levels between samples was used as an important index to test the reliability of the experiments and the rationality of sample selection. The closer the correlation coefficient was to 1, the higher the similarity of expression patterns between samples. A total of three correlation indexes were calculated: Pearson, Kendall and Spearman (R software; version 3.1.2).

DESeq2 (R package; https://github.com/mikelove/DESeq2; version 1.12.4) was used for gene expression difference analysis, and the results of the expression difference analysis were visualized. The DEGs were mapped to the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) protein interaction network database (http://string-db.org/), which was used to construct the protein interaction network. Based on the results of the differential analysis, Venn diagrams (VennDiagram R package; https://CRAN.R-project.org/package=VennDiagram; version 1.6.17) were drawn.

TopGO (R package https://bioconductor.org/packages/topGO/; version 2.24.0) was used for Gene Ontology (GO) enrichment analysis, and the significance of GO-directed acyclic graph was drawn. ClusterProfiler (R package https://guangchuangyu.github.io/software/clusterProfiler; version 3.0.5) was used for KEGG pathway and eukaryotic orthologous group (KOG) classification enrichment analysis. The association analysis network diagram was drawn based on the results of the gene function enrichment analysis.

The cells were lysed in cell lysis buffer for western blotting (Beyotime Institute of Biotechnology) supplemented with PMSF (100:1; Beyotime Institute of Biotechnology) to extract the total proteins. The concentrations of the proteins were detected using a bicinchoninic acid protein assay kit (Beyotime Institute of Biotechnology). Proteins (20 µg per lane) were loaded on an SDS gel (10%), resolved using SDS-PAGE and subsequently transferred to PVDF membranes (GE Healthcare Life Sciences), Membranes were blocked in 5% non-fat milk for 2 h at room temperature and incubated with primary antibodies against β-actin (ProteinTech Group, Inc; cat. no. 66009-1-Ig; 1:20,000) or β-catenin (Wanleibio Co., Ltd.; cat. no. WL0962a; 1:500) at 4°C overnight, and subsequently incubated with a HRP-conjugated secondary antibody [ProteinTech Group, Inc.; cat. nos. SA00001-1 (mouse) or SA00001-2 (rabbit); 1:5,000] for 2 h at room temperature. Signals were visualized using enhanced chemiluminescence reagent (EMD Millipore), and the protein level analysis was performed using ImageLab software (Bio-Rad Laboratories, Inc.; version 6.0.0).

Total RNA was extracted from cells using RNAiso Plus reagent (Takara Bio, Inc.). RNA was reverse transcribed to cDNA using an Evo M-MLV RT kit with gDNA Clean for qPCR (Hunan Accurate Bio-Medical Co., Ltd.), according to the manufacturer's protocol (37°C for 15 min, 85°C for 5 sec, then 4°C). The primers were synthesized by Sangon Biotech Co., Ltd. The sequences of the primers were as follows (5′-3′): Galactosylxylosyl protein 3-β-galactosyltransferase (B3GALT6) forward, CGACGCCTACGAAAACCTCA and reverse, GTCGTACACGTAGGACAGGC; N-acetyl galactosamine 4-sulfate 6-O-sulfotransferase (CHST15) forward, GCCACTCAATGCCATCCAGA and reverse, ATGGCAGGCTCGAGAACCAC; chondroitin 4-sulfotransferase 2 (CHST12) forward, CGACGAGTTTCTGGACAAG and reverse, GGACCAGTCGTAGCCTCT; and β-actin forward, AGCGAGCATCCCCCAAAGTT and reverse GGGCACGAAGGCTCATCATT. qPCR was performed using SYBR^® Premix Ex Taq^™ II, according to the manufacturer's protocol (Takara Bio, Inc.) as follows: 95°C for 30 sec, 95°C for 5 sec, 60°C for 30 sec for a total of 40 cycles, followed by 95°C for 10 sec, and melting curve analysis between 65 and 95°C (increments of 0.5°C) for 5 sec. β-actin was used as the reference to normalize the results and the data were analyzed using the 2^−ΔΔCq method (25).

Cells (1×10⁵/ml) were fixed in 4% paraformaldehyde for 20 min at room temperature and permeabilized using 0.1% Triton X-100 for 20 min. Samples were blocked in goat serum (Wuhan Boster Biological Technology, Ltd.; cat. no. AR009) for 30 min and subsequently incubated with the following primary antibodies: Anti-β-catenin (Wanleibio Co., Ltd.; cat. no. WL0962a; 1:200) and anti-mono-ADP-ribose binding reagent (EMD Millipore; cat. no. MABE1076; 1:200) at 4°C overnight. The following day, cells were washed in PBS and incubated with Cy3-conjugated secondary antibody (ProteinTech Group, Inc.; cat. no. SA00009-2; 1:200) in the dark for 1 h at room temperature. DAPI (Wuhan Boster Biological Technology, Ltd.) was used for nuclear staining (5 min at room temperature). The images were captured using ZOETM Fluorescent Cell Imager (Bio-Rad Laboratories, Inc.) and analyzed using ImageJ software (version 1.48; National Institutes of Health) The experiment was repeated three times.

Cells were seeded in 6-well plates. When the cell density reached 90%, the monolayers were scratched using a 200 µl sterile pipette tip, and the cells were cultured in serum-free medium for 24 h at 37°C. The wounded monolayers were imaged using a ZOETM Fluorescent Cell Imager (Bio-Rad Laboratories, Inc.) at the start of the experiment and after 24 h. The experiment was repeated three times and the images were analyzed using ImageJ version 1.48 (National Institutes of Health).

Transwell assays were performed using Transwell chambers (Merck KGaA) coated with Matrigel (Corning, Inc.). Serum-free cell suspension with 5×10⁵ cells were added to the upper chamber of the well at 200 µl/chamber, whereas supplemented medium with 10% serum was added to the lower chamber. After 24 h of incubation at 37°C, the cells were fixed with 4% paraformaldehyde for 30 min and stained with 0.1% crystal violet for 20 min at room temperature, followed by washing three times. Images were captured (magnification, ×200) under a light microscope (Nikon Corporation), the experiment was repeated three times, and the number of cells was counted using ImageJ version 1.48 (National Institutes of Health).

DEGs were determined using a negative binomial distribution, and enrichment analysis was performed using a hypergeometric distribution. DESeq2 was used to analyze differences in gene expression difference. SnpEff was used to calculate the distribution of variation sites in the genome structure, and the filters used were mass value >20 and coverage >8. Variable shear events were classified using ASprofile based on the predicted gene model for each sample. TopGO was used for GO enrichment analysis and ClusterProfiler was used for KEGG pathway and KOG enrichment analysis. All experiments were performed three times and the data are presented as the mean ± standard deviation. Statistical analysis was performed using GraphPad Prism version 5 (GraphPad Software, Inc.). An unpaired two-tailed Students' t-test was used to compare groups. P<0.05 was considered to indicate a statistically significant difference.

A previous study identified mono-ADP-ribosylation of H3R117 in Lovo cells using liquid chromatography with tandem mass spectrometry (LC-MS/MS) (14). In the present study, immunofluorescence experiments were performed to confirm the mono-ADP-ribosylation in these cells. Mono-ADP-ribosylation was notably present in Lovo cells, particularly in the nucleus (Fig. S1). These results demonstrated that nuclear mono-ADP-ribosylation is present in CRC.

FastQC was used to visually evaluate the quality of the raw data (Fig. 1) and Trimmomatic was used to filter the raw data to obtain the clean data. Data information prior to and following quality control showed Q20 >90 and Q30 >90% (Tables I and II). These results showed that the experimental data was of high quality, and the stability of the instrumental analysis systems was good and met the requirements for subsequent analysis.

Certain important quality metrics can only be measured when reads are matched with the reference genome and their positions are aligned to the annotation information, including sequencing saturation, distribution of reads with different genome characteristics and coverage uniformity of transcripts.

The reference genome was used as the reference sequence, HISAT2 was used to map the quality-controlled sequencing reads to the reference genome, and statistical comparisons were performed using RSeQC; ~97% of clean reads were mapped to the reference genome (Table III). Duplicate reads were defined as sequences with exactly the same base arrangement. Duplicate reads are primarily the result of PCR amplification during library construction. As shown in Fig. 1A and B, the overall trend of the curve showed a linear smooth extension and diffusion. As the horizontal axis increased, the corresponding vertical axis decreased, suggesting that the results were normal. Ideally, for RNA-seq, there should be independent sampling between sequencing reads, and the distribution of reads on all expressed transcripts should be homogenized. The analysis of gene coverage from 5′-3′ is an indicator used to evaluate the randomness of transcriptome sequencing (26). By mapping the results, the coverage of each transcriptional site was calculated, and the coverage rate of each transcriptional site was calculated. By calculating the coverage ratio of each gene in each sample, the proportion of all the samples sequenced completely detected and the proportion of genes not detected was calculated, and could be used to determine whether there were genes with specific expression patterns among the samples (Fig. 1C).

RNA-seq can be used detect various structural variations in mRNA. In the present study, SNPs and alternative splicing (AS) were primarily analyzed. SNP analysis showed that the untransfected cells possessed 60,362 SNP sites and 3,979 indel sites. The cells expressing the mutant H3 possessed 54,918 SNP sites and 3,876 indel sites (Fig. 2A). Of the SNPs, >50% were present in introns (Fig. 2D and E). SNPs included transition and transversion, and the majority of the SNPs in both groups of cells were transitions. A total of four transition types were detected: A>G, C>T, G>A and T>C. Among these, A>G and T>C were the most common. Transversion was relatively rarely observed, but eight types of transversion were detected: A>C, A>T, C>A, C>G, G>C, G>T, T>A and T>G (Fig. 2B and C). Between the two cell types assessed, the same gene exhibited different forms of AS, which greatly increases the capacity for gene encoding, as well as the complexity of gene expression and the diversity of protein function. Abnormal AS may produce potentially harmful or beneficial effects in the organism under certain environmental conditions. In the present study, ASprofile was used to classify AS according to the predicted gene models of each sample. The classification was as follows: SKIP, exon skipping; MSKIP, cassette exons; IR, retention of single; MIR, multiple introns; AE, alternative exon ends; TSS, alternative transcription start site; and TTS, alternative transcription termination site. The three types of AS most commonly observed in both cell groups were TSS, SKIP and TTS (Fig. 2F and G).

A total of three correlation indices were calculated in the scatterplot of the replicated correlation test: Pearson=0.8294, Kendall=0.8258 and Spearman=0.8911. These correlation indices were all close to 1, and the scatterplot showed that the majority of the points were centered around the diagonal line. These results indicated that the similarity of expression patterns between samples was high (Fig. 3).

In total, 58,174 differentially expressed genes were identified; among these, 2,324 genes were significantly different (q-value <0.05; fold change >2). Compared with the untransfected cells, 1,391 upregulated genes and 933 downregulated genes were identified. DEGs were mapped using the STRING protein interaction network database to construct protein interaction networks (Fig. 4A). Based on the results of the differential analysis, Venn diagrams were drawn. The Venn diagram showed the number of expressed genes that were common and unique between the two groups (TPM >0). A total of 21,669 genes were present in both groups; 2,400 genes were unique to the N group, and 3,848 genes were unique to the mutant cells (Fig. 4B).

GO is an internationally standardized functional classification system for genes, providing a set of regularly updated descriptions to comprehensively describe the properties of genes and gene products in organisms. There are three ontologies in GO: Molecular function, cellular component and biological process (Fig. 5A). According to the results of the GO analysis, it was identified that for biological process, DEGs were significantly enriched in ‘metabolic process’ and ‘negative regulation of biological process’ (q-value <0.05); for cell component, DEGs were significantly enriched in ‘organelle’, ‘organelle part’ and ‘membrane-enclosed lumen’ (q-value <0.05); for molecular function, DEGs were significantly enriched in ‘binding’ and ‘catalytic activity’ (q-value <0.05; Fig. 5B).

Clusters of Orthologous Groups of proteins (COG) and KOG are the NCBI annotation system based on direct homologous relationships of genes. Among these, COG targets prokaryotes and KOG targets eukaryotes. COG/KOG combines evolutionary relationships, homologous genes from different species are divided into different Ortholog clusters. Currently, there are 4,873 categories for COG and 4,852 categories for KOG. The genes from the same ortholog have the same functions. In the present study, KOG classification indicated that the DEGs were classed into 26 functional categories. The only categories that were statistically significant were ‘chromatin structure’ and ‘dynamics’ (q-value <0.05). A total of 43 genes involved in the functional category were identified. Compared with the untransfected group, 14 genes were downregulated and 29 genes were upregulated. The top four genes with a log2FoldChange were: Histone cluster 2 H3 pseudogene 2 (HIST2H3PS2) (log2FoldChange=22.65), histone cluster 2 H3 family member D (HIST2H3D) (log2FoldChange=5.35), CENPB DNA-binding domain containing 1 (CENPBD1) (log2FoldChange=2.84), histone cluster 1 H2A family member I (HIST1H2AI) (log2FoldChange=2.78; Fig. 5C).

KEGG is a relatively complete database of biological systems, integrating genomic, chemical and system functional information. The pathway database is the most commonly used sub-database (27). In RNA-seq, KEGG enrichment was used to analyze the DEGs, which were shown as scatter plots (Fig. 6A). KEGG enrichment was measured using the Rich factor, q-value and the number of genes enriched in this pathway. A total of 306 KEGG pathways were identified, and the top three most significant pathways are presented in Fig. 6B-D. These were ‘glycosaminoglycan biosynthesis-chondroitin sulfate/dermatan sulfate’ (ko00532), ‘ubiquitin mediated proteolysis’ (ko04120) and ‘bladder cancer’ (ko05219). Furthermore, the levels of expression of carbohydrate sulfotransferase 14 (CHST14), chondroitin 4-sulfotransferase 2 (CHST12), uronyl 2-sulfotransferase (UST), B3GALT6, chondroitin polymerizing factor 2 (CHPF2) were upregulated, whereas CHST15 and chondroitin polymerizing factor (CHPF) were downregulated in the glycosaminoglycan biosynthesis-chondroitin sulfate/dermatan sulfate pathway. RT-qPCR was used to confirm the changes in the top three genes based on the log2FoldChange values. These were CHST12 (log2FoldChange=3.61), B3GALT6 (log2FoldChange=2.34) and CHST15 (log2FoldChange=−2.41) (P=0.0075) (Fig. S2). The RT-qPCR results were consistent with the sequencing results. Similarly, the levels of expression of certain genes changed significantly in the ubiquitin-mediated proteolysis pathway and bladder cancer pathway.

To examine the function of mono-ADP-ribosylation in invasion and metastasis in CRC, Transwell invasion and wound healing assays were performed. The results demonstrated that compared with the untransfected cells, the cell invasion and metastatic capacity (P=0.0105) of the mutated cells were reduced (Fig. S3A and B). Previous studies have shown that the glycosaminoglycan biosynthesis-chondroitin sulfate/dermatan sulfate pathway is associated with invasion and metastasis. Thus, how this pathway was involved in tumor metastasis was assessed. Previous studies have reported that chondroitin sulfate/dermatan sulfate proteoglycans interact with Wnt, resulting in increased accumulation of β-catenin, which regulates tumor metastasis and other biological activities (28). In the present study, western blotting and immunofluorescence were performed to evaluate the expression of β-catenin, and the results showed that the expression of β-catenin was decreased significantly in the mutation group (P=0.0138 and P=0.0023; Fig. S3C and D, respectively). These results suggested that mono-ADP-ribosylation of H3R117 in Lovo cells regulates the biosynthesis of chondroitin sulfate/dermatan sulfate, and that this is involved in tumor metastasis by regulating Wnt/β-catenin signaling.