The patients diagnosed with PDAC were recruited from Changhai Hospital, Shanghai. All patients have signed a written informed consent form for collection and use of samples. A total of 41 pairs of tumors and their distal normal tissues were obtained from surgical specimens and frozen at -80°C immediately before nucleic acid extraction. All procedures have complied with the Code of Ethics of the World Medical Association (Declaration of Helsinki) and the guidelines of the institutional review committee of the Shanghai Institute of Nutrition and Health, Chinese Academy of Sciences.

Total RNA was isolated using the Qiagen RNeasy Kit (Qiagen, Germantown, MD) and treated with DNase I (New England Biolabs). Total DNA was extracted using the QIAamp DNA Mini Kit (Qiagen, Germantown, MD). The concentration was quantified by the NanoDrop One (Thermo Fisher Scientific, MA, USA), and the RNA quality was analyzed by the 2100 Agilent Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA). Libraries were generated using a standard library preparation kit (KAPA) according to the manufacturer’s instructions. Libraries were quantified on the Qubit4 (Thermo Fisher Scientific, MA, USA) with the High Sensitivity DNA Assay (Thermo Fisher Scientific, MA, USA). Samples were equimolar pooled and sequenced on the Illumina NovaSeq (Illumina, San Diego, CA, USA) in 150-bp paired-end mode.

The quality of raw reads was assessed by FastQC (https://github.com/s-andrews/FastQC), and subsequently the low-quality data were filtered by fastp (https://github.com/OpenGene/fastp). The hg19/GRCh37 of the human genome was used as reference. RNA-seq reads were aligned onto the reference genome using STAR (https://github.com/alexdobin/STAR) (parameters: –outSAMstrandField intronMotif –outSAMattributes All –readFilesCommand zcat –outFilterType BySJout –outFilterMultimapNmax 1 –alignSJoverhangMin8 –alignSJDBoverhangMin 1 –outFilterMismatchNmax 999 –outFilterMismatchNoverLmax 0.04 –alignIntronMin 20 –alignIntronMax 1000000 –alignMatesGapMax 1000000). WGS reads were aligned by BWA (https://github.com/lh3/bwa). All alignments were duplicated by Picard MarkDuplicates (http://broadinstitute.github.io/picard/), sorted, and indexed using SAMtools (https://github.com/samtools/samtools). Computationally, RNA editing is identified as a single-nucleotide base change between DNA and RNA. The RNA and DNA variants in each sample were both detected by REDItools (19) (https://github.com/BioinfoUNIBA/REDItools2) under unstranded strategy. Then, the sites containing DNA variants were removed and the rest were considered as RNA editing sites. Meanwhile, we also used SPRINT (20) (https://github.com/jumphone/SPRINT) to validate these RNA editing sites detected by REDItools.

The number of A-to-I RNA editing sites is huge, but most of the sites exhibit low editing levels and low sample coverage. These editing events are often caused by non-specific overediting, whereas most of these editing sites are neutral or slightly harmful. This poses a huge challenge for detecting editing sites with high confidence. Therefore, we further obtain credible RNA editing sites by satisfying stringent requirements (base quality score >30, total reads for each site in each sample ≥10, 1 > editing level for each site in each sample ≥0.1, remove sites with DNA mutation and multiple RNA variants, remove SNP sites) ( Figure 1A ).

RNA editing sites were classified into three groups, shared, normal-specific, and tumor-specific group, which implicated that RNA editing can not only produce new cancer-specific changes but also reedit existing editing sites in normal tissues to a certain extent. For the shared group, all editing sites were occurred in both tumor and normal samples and filtered by sample coverage over 10. The editing site was considered significantly different with a Fisher’s exact p value less than 0.05 based on editing level. For two specific groups, the editing sites were occurred only in normal or tumor samples and filtered by sample coverage (≥4).

Read counts per gene were calculated for each sample using hisat2 (http://ccb.jhu.edu/software/hisat2) and StringTie (http://ccb.jhu.edu/software/stringtie/). Differential expression analysis was performed underlying the Bioconductor package Deseq2 version 1.32.0 (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html), and the different expression genes were identified by filtering under log₂|fold-change| ≥1 and padj <0.05.

We used Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) molecular pathway enrichment analysis for the genes with RNA editing sites between tumor and normal tissues through R package “clusterProfiler” (http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html). Statistical significance was defined as p.adjust < 0.05.

The consensus motif in the MEME suite of motif-based sequence analysis tools (version 5.3.3, https://meme-suite.org/meme/) was identified using two motif widths of 30 and 50.

We used the database of Search Tool for Retrieval of Interacting Genes (STRING) (ver. 11.0, http://www.string-db.org/), an online tool to analyze protein network characteristics and make protein network visualization, to evaluate the network of protein–protein interaction (PPI) and calculate K-means for clustering.

To predict the RNA secondary structure alteration derived from individual RNA editing events, we used the extended surrounding region (± 1000 nt around the editing site) to search for minimum free energy (MFE) and partition function using RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). Meanwhile, MXfold2 (https://github.com/mxfold/mxfold2) was also applied to predict the RNA secondary structure of RNA sequence surrounding the editing sites (± 250 nt).

We used rMATS (https://rnaseq-mats.sourceforge.net/index.html) to detect and analyze alternative splicing events, including skipped exon, alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), mutually exclusive exons (MXE), and retained intron (RI).

Cell Ranger (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger) and Seurat (https://satijalab.org/seurat/) were applied to cluster and identify the cell types based on known marker expression. The marker genes were PRSS1, CTRB1, CTRB2, and REG1B for acinar cells, MS4A1, CD79A, CD79B, and CD52 for B cells, AMBP, CFTR, and MMP7 for type1 ductal cells, KRT19, KRT7, TSPAN8, and SLPI for type 2 ductal cells, CHGB, CHGA, INS, and IAPP for endocrine cells, CDH5, PLVAP, VWF, and CLDN5 for endothelial cells, LUM, DCN, and COL1A1 for fibroblast cells, AIF1, CD64, CD14, and CD68 for macrophage cells, ACTA2, PDGFRB, and ADIRF for stellate cells, and CD3E, CD4, and CD8 T cells (21). The reads for each cell type were extracted and combined. RNA editing events for each type were identified by REDItools. Then, credible RNA editing sites were filtered by satisfying stringent requirements (total reads for each site in each sample ≥10, 1 > editing level for each site in each sample ≥0.1, remove SNP sites).

The survival analyses were calculated by GEPIA2 (http://gepia2.cancer-pku.cn/#index) based on TCGA-PAAD datasets and ONCOMINE (https://www.oncomine.org) integrated pancreatic cancer data from TCGA and GEO.

F_ST analysis, as a genetic signature of positive selection, was performed using VCFtools.

Statistical analysis was performed by GraphPad Prism 9 (GraphPad Software, San Diego, CA, USA) and R packages. A value of p < 0.05 was considered statistically significant. Data were expressed as mean ± SEM.

To better understand the relationship between the RNA editing and PDAC, we respectively delineated the RNA editing profiles of normal and tumor samples and found that PDAC showed a severe A-to-I RNA editing imbalance. We analyzed the composition of high confident editing events and found that these events were enriched in A-to-I alternation (95.11%), which showed significant differences between tumor and normal samples (REDItools, paired t-test, p = 0.012) ( Figure 1B ). The distribution of RNA editing events in chromosome was also analyzed, and a preference in chromosomes 17 and 19 was observed ( Figure 1C , Figure S1C ).

We then assessed whether RNA editing sites showed differential editing levels between tumor and normal samples. A total of 1,657,409 and 654,485 unique A-to-I RNA editing sites were identified by REDItools and SPRINT, respectively, in which 504,659 editing sites were overlapped. The average number of A-to-I editing sites in tumor (average = 81504) was slightly higher than that in normal tissues (average = 70812) (paired t-test, p = 0.051) ( Figure 1D , Figure S1D ). A percentage of 56.10% of the samples showed a hyper average editing level in tumors than those in normal tissues without difference overall (REDItools, paired t-test, p = 0.567) ( Figure 1E , Figure S1E ). We also observed several genes that have been reported to change the editing level in cancer and even drive early tumor invasion and metastasis, including FLNB, SLC22A3, and AZIN1, which had significantly different editing levels in PDAC.

We further performed gene annotation and found that the number of editing sites per gene in tumor tissue was significantly higher than that in normal tissue (paired t-test, p = 0.002) ( Figure 1F ). The distribution map showed that most genes have multiple editing sites ( Figure 1G ). The types of RNA region showed an enrichment in intronic (80.99%), intergenic (9.29%), 3′UTR (3.07%), and exonic (4.06%) ( Figure 1H ).

From a quantitative and qualitative perspective, the imbalance of RNA editing detected by REDItools was divided into three modes, shared group, tumor-specific group, and normal-specific group, based on the distribution characteristics of RNA editing events that occurred in tumor and normal samples. For A-to-I editing events, 13,042 and 7,297 events (belonging to 3,816 and 2,594 genes, respectively) occurred in the tumor-specific and normal-specific groups, respectively ( Figure 2A ), whereas 9,343 events (belonging to 2,812 genes) with significant differences occurred in the shared group ( Figure 2A ). These editing sites among three groups were recognized as different RNA editing sites (DREs), in which 95.68% were identified by SPRINT. The level of DREs mainly ranged from 10% to 30% ( Figure 2B ). In the shared group, there were more sites with high editing levels in tumor samples. Moreover, in the specific group, a low level of editing sites was predominant. Compared with the public REDIportal database (http://srv00.recas.ba.infn.it/atlas/), 58.91% of all identified events were known and 79.35%, 85.19%, and 96.36% were known in the tumor-specific group, normal-specific group, and shared group, respectively ( Figure 2C ). The relative abundance of DREs showed a preference for chromosomes 17 and 19 ( Figures 2D, E ). The distribution of DREs in each gene was consistent with the total result ( Figure 2F ). Meanwhile, the types of RNA region of DREs showed a consistent result with total editing sites ( Figures S1G, H ), and low enrichment of DREs was found in TSS ( FigureS1I ). Eight A-to-I DREs were identified as non-synonymous mutation, which might lead to the change in protein ( Table S1 ). Furthermore, the DREs in the shared group displayed a hyper mean editing level and a higher editing event count in tumor than normal tissues (paired t-test, editing level: p < 0.001, editing event count: p = 0.023) ( Figures 2H, I ). In addition, the profile of total 29,682 A-to-I DREs could distinguish tumor and normal tissues well ( Figure 2G ).

To investigate whether the genes with A-to-I DREs are involved in the physiological process of pancreas or tumorigenesis, we obtained 108 enriched pathways for total 5,319 genes with DREs by functional enrichment analysis. Pancreatic cancer, adherens junction, and MAPK signaling pathway were enriched in the tumor-specific group, whereas pancreatic secretion only existed in the normal-specific group ( FigureS2A ). In general, RNA editing events under different modes displayed different functional mechanisms.

To determine whether ADAR family members have effects on RNA editing imbalance in PDAC, we compared the expression level of ADAR1 and its relation with the editing level and the amount of editing sites.

The expression of ADAR1 in tumor was significantly higher than that in normal tissues, rather than ADAR2 and ADAR3 (paired t-test, ADAR1: p < 0.0001, ADAR2: p = 0.357; ADAR3: p = 0.303) ( Figure 3A ). Meanwhile, TCGA and GTEx data analysis showed that the expression profiles of the ADAR family were consistent with our results ( Figure S3A ). The A-to-I mean editing level was positively correlated with ADAR1 expression level, whereas it was lightly negative correlated with ADAR2 (spearman test, ADAR1: p < 0.0001, coef = 0.2174, ADAR2: p = 0.0285, coef = 0.059, ADAR3: p = 0.2689, coef = 0.01526) ( Figure 3B ). The global association detected between ADAR1 expression and average editing level was also observed at each editing site ( Figure 3C ). Meanwhile, only ADAR1 observed a significant association with prognosis (ADAR1, p = 0.042; ADAR2, p = 0.29; ADAR3, p = 0.066) ( Figure S3B ). We then explored the triplet codon characteristics of A-to-I editing sites and found sequence preference, especially upstream G underrepresentation and downstream G excess, which was consistent with the known ADAR1 motif ( Figure 3D ).

Then, we further classified samples with ADAR expression between paired tumor and normal tissues into ADAR1-upregulated and -downregulated groups to explore the function of ADAR1 on regulating RNA editing events. We found that the mean editing level was positively correlated with the amount of editing sites in the ADAR1 upregulated group. However, the opposite result was observed in the ADAR1-downregulated group ( Figure 3E ).

Above all, these results suggest that the differential expression of ADAR1 manipulate the A-to-I imbalance in tumor samples in PDAC.

In order to explore the associations between RNA editing and gene expression in PDAC, differential gene expression analysis based on RNA-Seq was performed and 2,022 different expression genes (DEGs) were identified, including 712 upregulated and 1,310 downregulated genes in tumor ( Figure 4A ). The DEGs were also found to distinguish tumor and normal tissues, but not better than the DREs did ( Figure 4B ). Pathway analysis showed that downregulated DEGs were enriched in pancreatic secretion, PPAR signaling pathway, protein digestion and absorption, and cAMP signaling pathway, whereas upregulated DEGs were enriched in ECM–receptor interaction and IL-17 signaling pathway ( Figures 4C, D ). Notably, in the pathogenesis of PDAC, the expression of tumor-related genes changed greatly.

In order to detect whether the RNA editing events affect gene expression, we focus on the DEGs with DREs. A total of 140 genes were filtered out ( Figure 5A ), which included 40 genes with 73 editing sites in the shared group, 78 genes with 370 A-to-I editing sites in the tumor-specific group, and 61 genes with 266 A-to-I editing sites in the normal-specific group. These differentially edited sites accompanied by different expressions were named eDREs, and these genes were named as edDEGs. The genes in the shared group were enriched in pancreatic secretion; focal adhesion; protein digestion and absorption; regulation of actin cytoskeleton; alanine, aspartate, and glutamate metabolism; and ECM–receptor interaction ( Figure 5B ). For genes in the tumor-specific group, they were enriched in ECM–receptor interaction, focal adhesion, central carbon metabolism in cancer, regulation of actin cytoskeleton, and others. For genes in the normal-specific group, KEGG pathways were enriched in the PPAR signaling pathway; alanine, aspartate, and glutamate metabolism; pancreatic secretion; and protein digestion and absorption. Furthermore, only the gene set in the tumor-specific group observed a significant association with prognosis (total 140 edDEGs: p = 0.1, shared group: p = 0.022, normal-specific group: p = 0.13, tumor-specific group, p = 0.001) ( Figure 5C ). These results suggested that RNA editing might take a role in the pathogenesis of PDAC.

Then, the expression levels of the 140 edDEGs were validated using TCGA-PAAD datasets (http://gepia2.cancer-pku.cn/#index). A set of 93 edDEGs were filtered out to showed similar expression levels with our results and showed a significant association with prognosis (p = 0.013) ( Figure 5D ). We then analyzed the association between the level of RNA editing and the level of gene expression and found that it could be either positive or negative. For example, the level of the RNA editing event in chromosome 4: 110846540 on EGF showed a significantly negative correlation with its expression level (p = 0.0122) ( Figure 5E ).

Among them, 41 edDEGs were found to be significantly associated with prognosis in TCGA-PAAD dataset, the combination of which also showed a significant association with prognosis (p < 0.001) ( Figure 5F ). Moreover, the PPI network of 41 edDEGs was divided into three clusters, in which ITGA2, ITGA3, ITGB6, and LAMC2, MUC4, and B3GNT3 were the nodes with most connectivity ( Figure 5G ). These genes were reported to be related with the pathogenesis of PDAC. ITGA2 can promote cell migration, invasion, metastasis, and chemoresistance in PDAC through ECM remodeling (22). LAMC2 can promote cancer progression and gemcitabine resistance through modulation of EMT and ATP-binding cassette transporters in PDAC (23). ITGA3 has been reported with an upregulated expression in PDAC (24). MUC4-promoted neural invasion is mediated by the axon guidance factor Netrin-1 in PDAC (25). In summary, RNA editing indeed affects the expression of important functional genes, suggesting an important role of RNA editing in the pathogenesis of PDAC.

It was reported that RNA editing could affect RNA stability by altering the RNA secondary structure and alternative splicing (26). Thus, we analyzed the RNA sequences and the corresponding secondary structure close to the editing site. We predicted the potential motif surrounding 177 eDREs derived from the 41 edDEGs by MEME and found the preference of GA-rich and CA-rich regions in sequences ( Figure 5H ). Furthermore, we also found that 42.67% of these eDREs might impact the RNA secondary structure based on RNAfold analysis based on minimum free energy (MFE) and partition function ( Table S2 ), which in turn might affect RNA stability and function. For example, the predicted RNA structure was changed from stem to loop at a GA-rich region due to the nearby RNA editing event in chromosome 11: 117986470 on TMPRSS4 ( Figure 5I ), which has been reported as an independent risk factor in PDAC to promote cell proliferation and inhibit apoptosis through activating ERK1/2 Signaling pathway (27).

A total of 3,647 significantly different alternative splicing events were identified between tumor and normal, of which 41.98% contained at least one RNA editing site ( Figure 5J , Table S3 ). These differentially edited sites accompanied by different alternative splicing events were named sDREs, and these genes were named edDSGs. Major sDREs participated in skipped exons. Among them, LAMC2, RAB27B, PHYHD1, MYH11, and CERS4 were also occurred in 140 edDEGs with co-changes in expression and editing. The alternative splicing events of LAMC2, RAB27B, MYH11, and CERS4 showed that once an RNA editing event occurred on the upstream intron, at least one downstream exon was shipped ( Figure 5L ). Exon 3 of CERS4 is a protein-coding region, and its loss because of A-to-I editing in intron2 would result in truncated protein, thereby affecting its protein function. Two kinds of sDREs in intron 2 of RAB27B lead to an alteration of the RNA secondary structure and also generated an abnormal transcript ( Figure 5L ). PHYHD1 displayed alternative 5′ splicing sites due to A-to-I editing in intron 10 leading to a truncated exon 10 ( Figure 5L ). We further observed the editing levels among various alternative splicing types and found that hyper editing events accounted for a large proportion of each component in the tumor ( Figure 5K ). These suggested that aberrant RNA editing in tumors may lead to aberrant variable splicing through an unstable RNA secondary structure.

These results suggest that RNA editing regulates gene expression and protein synthesis by affecting alternative splicing and the RNA secondary structure.

It was reported that positive selection is one of the mechanisms driving RNA editing (28, 29). In order to determine whether the RNA editing events suffer from positive selection in the tumorigenesis, we performed calculation of F_ST for each DREs. A percentage of 5.48% of DREs belonging to 1,005 genes showed the signal of positive selection ( Figure 6A ). These genes were enriched in regulation of actin cytoskeleton, focal adhesion, and adherens junction ( Figure 6B ), and 48 genes were immune-related genes. Furthermore, we screened out the positive selection signals of all immune genes ( Table S4 ). Chromosomes 19 and 17 had relatively more selection signals ( Figure 6C ), which was consistent with the enrichment of editing events in chromosomes 19 and 17. Meanwhile, an interesting result was found in which the positively selected DREs were mainly located in the normal- and tumor-specific groups. The hot genes with the most positive signal events in the normal-specific group were RBPJ, PFKFB3, CAMK1D, ACACB, and WWOX, whereas the hot genes in the tumor-specific group were SLC35F3, PCDH7, NF1, and RAB27B. Among them, RBPJ was an important transcriptional regulator in the Notch signaling pathway (30). PFKFB3 was required for cell-cycle progression and prevention of apoptosis (31). RAB27B expression is considered an independent prognostic marker for PDAC (32).

In summary, positive selection indeed has an effect on the RNA editing in PDAC and has a preference for the RNA editing sites in tumor immune-related genes.

We downloaded a single-cell RNA-seq data set (PRJCA001063) of pancreatic cancer containing 24 tumors and 11 normal pancreases and identified 10 major clusters including type 1 ductal, type 2 ductal, acinar, endocrine, endothelial, fibroblast, stellate, macrophage, T, and B cells, through principal component analysis based on gene expression. We then extracted the RNA editing events for each cluster. In general, the significant difference in editing level occurred in type 2 ductal, macrophage, and T cells, whereas the significant difference in editing sites occurred in type 1 ductal, type 2 ductal, macrophage, stellate, and T cells ( Figure 7A ). Among them, type 2 ductal cells greatly contributed to RNA editing events in pancreatic cancer tissues, whereas type 1 ductal cells had the largest number of RNA editing events in normal tissues. We then identified the DREs and their located genes and found that type 2 ductal, endothelial, fibroblast, stellate, macrophage, and T cells were dominated in tumor-specific editing sites, whereas normal-specific editing sites were dominated type 1 ductal cells ( Figure 7B ). The shared group showed a significantly higher editing level in type1 ductal cells and a significantly lower editing level in fibroblast, stellate, and macrophage in tumors, whereas in the specific group, a significantly higher editing level occurred in acinar, B, type 2 ductal, endocrine, macrophage, and T cells and a significantly lower editing level was displayed in stellate cells ( Figure 7C ). The Pearson analysis showed that it could be distinguished into normal and tumor based on the level of DREs ( Figure 7D ). Then, we compared the distribution of genes with RNA editing sites between RNA-seq of 41 paired PDAC samples and single-cell sequencing data of 35 PDAC ( Figure 7E ). There were 187 overlapped observed only in the tumor-specific group between RNA-seq and single-cell sequencing, of which 72.06% were detected in type 2 ductal cells.

In conclusion, the level and amount of RNA editing events were different among various cell types of PDAC, and type 2 ductal cells contributed the most to RNA editing in tumors.