A total of 16 healthy individuals and 15 pancreatic cancer patients were recruited between April 2015 to August 2019 with age range of 20 to 70 years for PC patients and 20 to 55 years for normal individuals. In both PC and normal individuals the female to male ratio was about 3:2. Surgical tissue and plasma samples of patients with confirmed PC (pancreatic ductal adenocarcinoma) and not undergoing any chemotherapy were obtained from the Institute of Postgraduate Medical Education & Research and the Chittaranjan National Cancer Institute (both Kolkata, India). The study was approved by the Institutional Ethics Committee (INST/IEC/2015/218 and IPGME&R/IEC/2022/318 for Institute of Postgraduate Medical Education & Research-Research Oversight Committee and CNCI-IEC-SG2-2023-69 for Chittaranjan National Cancer Institute-Institutional Ethics Committee). Written informed consent was procured from all participants prior to the study. A total of 5 ml peripheral venous blood was collected in vacutainer tubes (BD Biosciences) before routine surgery and plasma samples were processed as previously described (17). Normal plasma samples were collected from healthy individuals with no history of pancreatic disease and were processed in the same way. Tumour and adjacent normal pancreatic tissue (>5 cm from tumour margin) were collected from patients with PC. The samples were stored at -80˚C until use. Simultaneously, resected specimens were processed for histopathological assessment to confirm malignant or benign nature (Table SI).

Total RNA enriched with small RNAs was isolated from the plasma samples using miRNAeasy Serum/Plasma advanced kit (Qiagen GmbH). Briefly, 200 µl plasma sample was centrifuged at high speed of 16,000 x g for 5 min at 4˚C to remove any cellular debris or particulate matter that may interfere with downstream RNA isolation. Next, QIAzol lysis reagent was followed by vortexing and phase separation after adding chloroform as per the manufacturer's instructions. The aqueous phase was then separated, and ethanol precipitation of RNA was performed followed by passage through the column provided with the kit. Column-bound RNA was washed and eluted using the buffer solution, according to the manufacturer's instructions. Quantification was performed using a multi-channel spectrophotometer (ND 8000; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. Additionally, total RNA was isolated from tissue using QIAzol (Qiagen GmbH) and PureLink RNA mini kit (Ambion; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. The quality of isolated total RNA was determined using Agilent RNA 6000 Nano chips in a 2100 Bioanalyzer (Agilent Technologies, Inc.) and NanoDrop spectrophotometer (Thermo Fisher Scientific, Inc.). Quantification was performed using Qubit and the Quant-iT RNA assay kit broad range (Thermo Fisher Scientific, Inc.). Total RNA samples with RNA integrity number >7 were selected for small RNA library preparation and Illumina sequencing.

The quality of isolated small RNA was checked using small RNA chips in a 2100 Bioanalyzer (Agilent Technologies, Inc.) and quantitation was performed using a Qubit Fluorometer. Small RNA sequencing library preparation was performed using Illumina^® TruSeq^® Small RNA Library Prep kit (Illumina, Inc.; cat. no. RS-200-0012) according to the manufacturer's instructions. A total of 10 ng isolated small RNA was used for library preparation. The first step was to ligate adapters to 3' and 5' ends of the RNA molecule. Subsequently reverse transcription and amplification were performed to generate a cDNA library, using the reagents provided in kit following manufacturer's instructions. Gel purification step that selects bands 145-160 bp long was performed to prepare the final small RNA sequencing library for clustering and sequencing. The quality of small RNA sequencing libraries was checked using high sensitivity D1000 screen tape in a 2200 TapeStation (Agilent Technologies, Inc.) and final library quantification was performed using a Qubit Fluorometer (Thermo Fisher Scientific, Inc.). Single end 1X 50 bp sequencing of pooled libraries was performed in a Novaseq 6000 (Illumina, Inc.).

Initial quality control and visualization of small RNA sequencing data were performed using FastQC (version 0.12.0) (18) and MultiQC (v1.24) (19). Adapter trimming (Illumina TruSeq small RNA adapters) was performed using the TrimGalore tool (v0.6.10) (20). Sequence reads of a length of 24-35 nucleotides were retained in the analysis. Poor quality reads (Phred score <20) were filtered out.

Alignment of reads to the reference genome and quantification of piRNA expression. The filtered sequencing reads were aligned to the human genome reference (hg19) using Bowtie2 aligner (Version 2.5.1) (21). Quality-checking of the sequencing alignment data was performed using SAMtools (v1.21) (22), Sambamba (v0.5.0) (23) and Qualimap (v2.3) (24). Aligned sequencing reads were overlapped with other small RNA sequencing information from the DASHR (v2.0) database (in BED file format) (25) to filter out other small ncRNA. Raw piRNA expression counts were quantified using the Featurecounts tool (v1.6.0.3) (26), with a piRNAdb annotation file (version 1.7.5; reference genome, hg19) (27) in GTF format. For normalization of the count data, ‘estimateSizeFactors’ was used to calculate size factors for each sample, using the Median ratios method and normalized counts were obtained using ‘counts’ of DESeq2.

Analysis of differential piRNA expression. The R package DEseq2 (version 1.12.3) (28) was used to identify DE piRNAs. Wald test was used for assessing statistical significance with adjusted P-value <0.1 as the threshold and -log2FC >0.58 or <-0.58 for up- and downregulated piRNA, respectively. Visualization of DE piRNAs was performed through heatmap and volcano plots, using R packages pheatmap (10.32614/CRAN.package.pheatmap) and dplyr (10.32614/CRAN.package.dplyr). A detailed schematic of the data analysis pipeline and quality filtering is shown in Fig. S1.

Differentially expressed piRNAs and protein coding genes were used to predict the target genes for piRNAs using miRanda target prediction algorithm (29). Briefly, the FASTA sequences of DE piRNAs and gene coding transcripts were retrieved from piRNAdb (version 1_7_5; hg19 reference; pirnadb.org/) and Ensembl databases(Homo_sapiens.GRCh37.cdna.all.fa;hg19reference; ensembl.org/index.html), respectively, and miRanda was run using alignment score ≥170 and a binding energy ≤-20 kcal/mol. Next, the piRNA-target gene pairs were filtered based on degree of sequence complementarity in the primary (2-11 nucleotides) and secondary seed site (12-21 nucleotides) with no mismatch and wobble base pairing within the primary seed site and only one mismatch and no wobble base pairing within the secondary seed site, as reported previously with minor modifications (30). Additionally, targets for each piRNA were manually verified using piRNADb (27).

To understand the functional aspects of targets of DE piRNAs, gene set enrichment analysis (GSEA) was used. Gene Ontology (GO) (geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (genome.jp/kegg/) were used to perform the enrichment analysis. KEGG database integrates genomic information to explore metabolic pathways, genetic information processing and cellular processes, whereas GO analysis involves the computational examination of gene sets to identify and categorize biological processes (BPs), cellular components (CCs) and molecular functions (MFs). Enrichr was used for GO and KEGG analysis (31,32). P<0.05 was considered to indicate significant enrichment. The Cancer Genome Atlas-Pancreatic Adenocarcinoma (TCGA-PAAD) dataset was used through Gene Expression Profiling Interactive Analysis; gepia.cancer-pku.cn/).

Most piRNAs in the genome originate from 25-35-bp-long discrete loci termed ‘piRNA clusters’, which serve a key role in the silencing of TEs (33). proTRAC command-line tool (34) was used with default parameters (sliding window size, 5,000 bp; sliding window increment, 1,000 bp; minimum fraction of hits with 1T(U) or 10A, 0.75; minimum size of piRNA cluster, 1,000 bp) to identify piRNA clusters in the tissue and plasma samples. Overlap of identified clusters with known repetitive elements in the genome was determined with Repeatmasker database (reference genome, hg19) annotation (35) and BEDtools package (36).

Wald test was used for assessing statistical significance with adjusted P-value <0.1 as the threshold and -log2FC >0.58 or <-0.58 in R package DEseq2 was used to identify DE-piRNAs. Pearson's correlation test. In GO enrichment analysis, hypergeometric distribution mathematic model was used to obtain the P-value of the Pathways. In KEGG pathways, multiple Benjamini and Hochberg testing was performed.

Small RNA sequencing was performed to obtain a median of ~52.45 million reads (range, 28.32-65.68 million) in tumor tissues and adjacent normal samples. A median of 15.79 million reads (range, 3.29-92.69 million) was obtained in the case of the plasma samples (PC and CP cases and respective normal samples). Detailed metrics and quality control information of small RNA sequencing in tissue and plasma samples of patients with pancreatitis and PC are provided in Table SII.

In pancreatic tumor tissues compared with normal tissue, 30 piRNAs were significantly upregulated and six piRNAs were significantly downregulated (Table I; Fig. 1A and C).

There were 16 normal samples and 15 samples from patients with PC. One piRNA was up- and 10 were downregulated (Table SIII; Fig. 1B and D).

Matched plasma and tissue specimens were evaluated for piRNA expression profiles in four PDAC plasma samples. Normalized expression levels of 20 piRNAs were positively correlated in plasma and tissues of patients with PDAC (Pearson's correlation coefficient >0.3; Table SIV). This suggests a common trend of piRNA deregulation between biospecimens and increases the chances of tumour tissue piRNA alteration being reflected in plasma. A total of three piRNAs, hsa-piR-23246, hsa-piR-32858 and hsa-piR-9137, were highly correlated between tumour tissue and plasma of patients with PDAC and may have an important role in disease development.

piRNAs have been found to operate in a manner similar to miRNA, as evident from studies, regulating gene expression through complementary base pairing and exhibiting an inverse correlation with target mRNA expression (9,10,37). By using Miranda algorithm and piRNAdb database to predict the targets of DE piRNAs, 413 mRNA targets for six downregulated DE piRNAs and 1,984 targets for 30 upregulated DE piRNAs were identified (Table SV).

To determine the role of deregulated piRNAs in patients with PDAC, pathway analysis using predicted genes of DE piRNAs was performed. The top 10 GO classifications of BP, CC and MF were used (Fig. 2A and B). Using targets of upregulated piRNAs, top GO processes included ‘monoatomic cation transport’ and ‘maturation of SSU-rRNA’ in the BP subgroup, ‘pre-ribosome, small subunit precursor’ and ‘90S pre-ribosome’ in the CC and ‘mRNA 5' UTR binding’ and ‘histone demethylase activity’ in the MF subgroup (Table SVI, Table SVII, Table SVIII, Table SIX, Table SX and Table SXI). Similarly, using downregulated piRNA and targets, top GO processes were identified as cAMP-mediated signaling regulation, regulation of transcription and regulation of arginine-histamine methylation in BP group. Rough endoplasmic reticulum membrane and RISC complex were the top CC subgroup, while histone demethylase activity and promoter-specific chromatin binding were the top GO processes in the MF subgroup. Additionally, 154 KEGG pathways were identified among dysregulated piRNAs and their targets (Tables SXII and SXIII). Pathway enrichment analysis showed that several key pathways such as pathways of glycolysis and gluconeogenesis, pathways of bile secretion, pathways in cancer, glucagon signaling pathway, ‘insulin signaling pathway’ and MAPK signaling pathway were enriched (Fig. 3A and B) (Tables SXII and SXIII). These results indicated alteration of multiple malignancy-specific pathways in PC.

proTRAC was used to predict genomic location enriched with piRNA clusters in PC and normal samples. A total of 262 clusters were identified in plasma samples (Table SXIV) containing 40 piRNAs. Among these, 18 were exclusive to PC patients, 15 were common to both patient and normal samples and seven were unique to normal samples. A total of 34 clusters was identified in the patient genome while analyzing the tissue samples. Out of these, 17 clusters showed high normalized count reads. From these high-density clusters, 10 clusters demonstrated the presence of 25 enriched piRNAs (Table SXV). Similarly, within normal tissue samples, 24 clusters were identified alongside 12 highly concentrated clusters. One cluster was the origin of six piRNAs, of which hsa-piR-32859, hsa-piR-22269, hsa-piR-20792, hsa-piR-32002 and hsa-piR-15181 were DE.

Although there was an asymmetric distribution throughout the chromosomes, clusters of piRNAs were more prevalent in chromosomes 9, 10, 18, 20 and Y (Table SXV; Fig. 4A). There were no clusters detected on chromosomes 1, 4, 5, 6, 7, 8, 12, 14, 19, 22 and X. Furthermore, the nucleotide preference among the piRNAs was also investigated. Specifically, piRNAs encoded from the plus strand exhibited the strongest preference for 1U compared with those from the minus strand (9,10). However, a higher bias for 10A was observed in the minus compared with the plus strand (Fig. 4B). This indicated the importance of orientation in acquiring or determining the type of biogenesis pathway of piRNA production. Based on this genomic feature, orientation is key in determining the type of pathway by which piRNA is produced. Fig. S2 shows a representative image of piRNA cluster visualization.

TEs or jumping genes occupy 45% of the entire human genome. TEs such as long Interspersed Nuclear Elements and SINE (Short Interspersed Nuclear Elements) are small repetitive sequences and easily enter or jump into any position of genome (38). This type of insertion leads to mutations and contributes to cancer development (38). To understand how piRNAs silence TEs, their origin was analyzed from the genome cluster and targeted coordinates in those regions were identified (Table SXV; Fig. 5) (39-41).

There was a higher density of piRNA origin found in tumour samples compared with normal samples. Additionally, density of piRNA clusters present in tumour samples was higher in SINE repeats. Long Terminal Repeats regions that harbour piRNAs were identified only in tumour samples. In normal samples, piRNAs were mapped to the DNA transposons like Tc1/mariner, DNA/TcMar-Tigger, DNA/hAT-Charlie and DNA/hAT-Tip100.

Pre-operative plasma samples from eight patients with CP were used and circulating piRNA expression pattern was compared with that of healthy control individuals (n=16). A total of four upregulated piRNAs and 15 downregulated piRNAs were found in the plasma of patients with CP (Table SXVI; Fig. S3). piRNAdb was used to identify potential targets of those piRNA (Table SXVII). Most of the target genes contribute to chronic inflammation in multiple organs (Table II). The results of the present study suggested similar changes of piRNAs in pancreatic tissues of patients with CP that cause deregulation of target genes contributing to chronic inflammation. Additionally, the involvement of these target genes in PC was investigated using TCGA-PAAD. Almost all the genes were upregulated in pancreatic tumour tissue, indicating that the inflammatory condition in CP may also be present in pancreatic tumour tissue and promote carcinogenesis.