Skip to main content

Integration of multi-omics data to unveil the molecular landscape and role of piRNAs in early-onset colorectal cancer

Abstract

Background

The incidence of early-onset colorectal cancer (EOCRC) (< 50 years) has been steadily rising, with a parallel increase in metastatic and invasive cases. To elucidate the molecular mechanisms underlying this aggressive phenotype, we performed comprehensive multi-omics profiling to delineate the distinct features of EOCRC, with a focus on key drivers of metastatic and invasive potential.

Methods

We initially characterized the genome, epigenome, and transcriptome of tumors from 515 (69 EOCRC and 446 late-onset CRC [LOCRC]) cases in The Cancer Genome Atlas. Key candidate molecules were further validated using RNA-seq and scRNA-seq data. Multi-omics profiling revealed PIWIL1/piRNA as a hallmark of EOCRC, with further validation through in vitro functional assays, transcriptomic profiling, and Kaplan-Meier survival analysis.

Results

EOCRC demonstrated a mutational landscape similar to that of LOCRC, with comparable oncogenic driver mutations and somatic copy-number alterations. However, EOCRC exhibited a higher frequency of deletion in chromosomes 6, 15, and 19 regions, along with metabolic reprogramming favoring aerobic glycolysis and lipid metabolism. Integrative transcriptomic and DNA methylation analyses identified six EOCRC-specific molecules, including PIWIL1. Notably, PIWIL1 was mainly expressed in epithelial cells, with lower expression in EOCRC versus LOCRC. Its downstream piRNAs (FR019019, FR019089, and FR132045) were also downregulated in EOCRC. Functional experiments demonstrated that FR019089/FR019019 overexpression suppressed migration and invasion. Clinically, low FR019089 levels correlated with significantly shorter progression-free and overall survival in EOCRC patients. Additionally, downstream pathways of FR019089 and FR019019 overexpression were enriched in anti-cancer-related signaling pathways.

Conclusions

Our multi-omics approach yields novel insights into the molecular underpinnings of EOCRC and we characterize the role of PIWIL1-associated piRNAs in modulating EOCRC metastasis and invasion. FR019089 shows promise as a prognostic biomarker with potential clinical utility in the risk stratification and management of EOCRC patients.

Graphical Abstract

Peer Review reports

Background

Colorectal cancer (CRC) ranks as the third most diagnosed cancer and the second leading cause of cancer-related mortality worldwide, representing a significant global health burden [1]. Despite a global decrease in CRC incidence over the last 40 years, early-onset colorectal cancer (EOCRC, defined by a diagnosed age under 50 years) has been steadily increasing at an alarming rate of 2–4% annually [2,3,4,5]. The EOCRC tends to be localized in the distal colon and rectum [6]. Owing to the lack of routine screening in younger patients and the inherently more aggressive features of EOCRC tumors, a larger proportion (85%) of EOCRC patients are symptomatically at diagnosis, with a significantly higher prevalence of advanced-stage disease than the overall CRC population (50%) [6, 7].

Environmental factors have been implicated in EOCRC pathogenesis, yet the molecular changes driving its widespread rise remain elusive [2, 8, 9]. Notably, EOCRC exhibits distinct molecular alterations compared to late-onset CRC (LOCRC). A few studies have attempted to profile the somatic mutational landscape of EOCRC patients, revealing that somatic mutations in oncogenic driver genes differ quite a bit between EOCRC and LOCRC [7, 10]. A significant mutation frequency in PIK3R1, PDGFRA, FLT3, and KDR gene mutations has been identified [6]. Another study integrated deconvolution or mutational signatures of EOCRC and epidemiologic data to uncover mutagenic processes that contribute to tumorigenesis across the age continuum [11]. However, existing multi-omics studies on EOCRC remain largely descriptive, with limited functional validation of key molecular drivers implicated by omics data.

Beyond genomics and epigenomics, transcriptomic studies have identified mRNA-based signatures for EOCRC [12,13,14]. Nevertheless, the precise involvement of key non-coding RNAs (ncRNAs), such as PIWI-interacting RNAs (piRNAs), in EOCRC remains poorly understood. Recently, evidence gathered has found that small ncRNAs (including piRNAs) are closely related to the occurrence and progression of cancers, highlighting their potential as diagnostic or prognostic biomarkers [15]. Piwi Like RNA-Mediated Gene Silencing 1 (PIWIL1), associated with a class of small ncRNAs called piRNAs of 20–35 nt, plays a crucial role in somatic cell function and increased CRC progression [16]. PIWIL1 overexpression correlates with poor tumor differentiation, infiltration, lymph node invasion, metastasis, and reduced overall and disease-free survival in CRC [17]. Nevertheless, detailed exploration of PIWILs and piRNAs in CRC, especially in EOCRC, remains limited, leaving the cellular pathways and molecular mechanisms involving somatic piRNAs poorly understood.

The study identified distinct genomic and epigenomic alterations in EOCRC through integrated multi-omics analysis, demonstrating the power of multi-omics approaches in elucidating age-specific tumor heterogeneity. Building on these findings, we investigated the functional mechanisms by which key piRNAs mediate the enhanced migratory and invasive properties in EOCRC.

Methods

Study design

This study integrated multi-dimensional omics data, including whole-genome sequencing, bulk/single-cell RNA-seq, and DNA methylation profiling. Genomic data were utilized to characterize the mutational status and somatic copy number variations (CNVs) related to EOCRC and LOCRC. Based on transcriptome data, we elucidated the feature of tumor microenvironment and metabolic shifts related to age of onset CRC. To probe into the EOCRC-specific molecular events, we analyzed the differentially expressed genes (DEGs) at transcriptional level and differentially methylated CpGs (DMCs) at the epigenetic level; pivotal genes were then validated through integrating the DNA methylation and gene expression data. Single-cell RNA-sequencing data were further analyzed to specify the cell types with expression of the key candidate genes in CRC tumors. Based on multi-omics findings, we investigated the functional role of PIWIL1/piRNA through in vitro experiments.

Patients and clinical tissue samples

Multiple omics datasets were obtained from two databases, namely the TCGA [18] (N = 515) and the Gene Expression Omnibus (GEO) datasets under accession code GSE39582 [19, 20] (N = 585), GSE17536 [21, 22] (N = 177), and GSE17537 [21, 23] (N = 55). The genomics, transcriptomics, epigenomics, and clinical information of 515 European individuals with CRC (a combination of COAD and READ) from the TCGA were downloaded. To validate the expression of candidate biomarkers for EOCRC, expression profiles from three GEO datasets were extracted. We downloaded scRNA-seq data (GSE132465) [24, 25], comprising 23 primary tumor samples from EOCRC patients and 10 matched normal mucosa samples. The datasets mentioned above were downloaded from several websites (https://portal.gdc.cancer.gov/; http://www.cbioportal.org/; https://www.ncbi.nlm.nih.gov/geo/).

The acquisition of all tissue samples utilized in this research was conducted following the approval granted by the Biomedical Ethics Committee of Zhejiang University (ZGL202306-3), with each patient providing written informed consent. A total of 18 EOCRC patients who were firstly diagnosed with CRC were recruited from the First Affiliated Hospital of Zhejiang University School of Medicine between 1 st Mar, 2023 and 19 th Feb, 2024. Eighteen pairs of colorectum tumor and control biopsy specimens were obtained during laparoscopic surgery. The control tissue was collected at least 2 cm from the edge of the CRC tissue. Moreover, all CRC and control tissues were examined and evaluated by at least three pathologists. The research was conducted following recognized ethical guidelines.

Genomic data processing

The mutation annotation format file with aggregated mutation of CRC cases was obtained from the TCGA portal. Candidate somatic single-nucleotide variants were identified using Mutect algorithms. To identify significantly mutated genes, we used the “maftools” R package to analyze the MAF files of the TCGA mutation data. The “lollipopPlot” function was applied to display mutation types in the two different subgroups, and results were visualized using Lollipop charts. GISTIC software was widely applied to detect both broad and focal (potentially overlapping) recurring events. We used GISTIC 2.0 [26] to identify genes exhibiting significant amplification or deletion. The parameter threshold was defined as amplification or deletion length > 0.1 and P < 0.05.

Identification of differential genes and differentially methylated CpG sites

We defined protein-coding genes differentially expressed between EOCRC and LOCRC tissue samples using the “edgeR” package on log2-transformed gene expression levels (quantified by transcript per million). Additionally, to account for differences in the distribution of baseline characteristics between EOCRC and LOCRC, the inverse probability of treatment weighting (IPTW) using propensity scores was performed to balance baseline covariates associated with each group [27]. Gender and stage at diagnosis were included in the IPTW process as covariates. By this approach, all patients in each group contributed to the final analysis with a specific calculated weight. Weighted gene expression data was analyzed to support the stability of DEGs results without balancing baseline covariates. Genes with an absolute log2-fold change (FC) ≥ 1 and false discovery rate (FDR) < 0.05 in comparison were considered DEGs. GO analysis was conducted using the enrichGO function implemented in the R package “clusterProfiler” (v. 3.8.1) [28].

DNA methylation array data (Illumina Infinium HumanMethylation450 BeadChip, HM450) were obtained from the GDC Data Portal, which interrogated 485,764 cytosine positions of the human genome, out of which 482,421 positions (99.3%) were CpG dinucleotides [29]. We applied the following criteria for quality control: (i) probes with mean detection P value ≥ 0.01 in > 5% of samples were removed from all samples; (ii) probes on the X or Y chromosome were removed; (iii) probes overlapping with single-nucleotide polymorphisms were removed; (iv) probes mapped to multiple sites in the human genome were removed. Finally, 308,242 probes were kept for further analysis. The CpG probe annotation file was downloaded from the ENCODE Project database (http://genome.ucsc.edu/ENCODE/downloads.html). Each CpG probe was annotated with the corresponding gene [30]. The package “limma” (implemented in “ChAMP” packages [31]) was used to identify DMC sites between EOCRC and LOCRC patients. A linear model was conducted to calculate the P value for differential methylation, and CpG probes meeting the criteria of β value ≥ 0.1 and FDR < 0.05 were defined as significant differential CpGs.

Correlation analysis between DNA methylation and gene expression

The correlations between the methylation levels of DMCs and the expression levels of their corresponding genes were examined using Spearman’s rank correlation. To obtain the DNA methylation driver genes, we considered FDR < 0.05 and the absolute Spearman rank correlation coefficient |r|> 0.30 as a statistically significant correlation. We further considered the consistency of direction between the methylation level of DMCs and the corresponding gene expression level for a more rigorous screening. DMCs located in promoters were selected in this work, and we applied the following criteria: (i) Wilcoxon rank-sum test was used to detect DNA methylation differential genes with P < 0.05. (ii) Negative correlations (Spearman r < − 0.30, FDR < 0.05) between DMCs and their corresponding genes.

piRNA transcriptome data collection and processing

To identify the expression levels of candidate piRNAs between EOCRC and LOCRC tissues, we obtained the piRNA expression matrix of 341 CRC patients (44 EOCRC patients and 297 LOCRC patients) from the TCGA dataset, which was characterized by Martinez et al. [32]. piRNAs with an absolute log2FC ≥ 1 and FDR < 0.05 were considered EOCRC-specific piRNAs using the “edgeR” package.

Pathway enrichment and immune signature analysis

GSEA was performed using the “GSVA” package [33]. Enrichment results with a criterion of FDR < 0.05 were considered significant. Gene sets were derived from the Broad Institute’s Molecular Signatures Database (MsigDB). In addition, the fatty acid-associated genes and glycolysis-associated genes were obtained from the hallmark, reactome, and kegg gene sets listed in MsigDB. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the R package “clusterProfiler.” The CIBERSORT platform (https://cibersort.stanford.edu/) was used to estimate immune cell proportions based on the LM22 reference profile. Computing gene set scores for representative immune signatures (immune checkpoint, immunosuppression modulators for TC-6, immune response to tumor cells, etc.) was performed using the “ssGSEA” algorithm.

Single-cell transcriptome analysis

We followed the pipeline to process the dataset [24]. Briefly, the sequencing data was aligned with the human reference genome (GRCh38) and processed using the CellRanger 2.1.0 pipeline (10 × Genomics). The original gene expression matrix was filtered using the “Seurat” package [34] based on the cell counts with > 1000 unique molecular identifier (UMI) counts, > 200 genes and < 6000 genes, and > 20% mitochondrial gene expression in UMI, which filtered gene expression matrix with 67,296 cells. After defining the global cell type, cells with several genes exceeding the outliers were removed to eliminate potential doublets, which generated a final gene expression matrix with 63,689 cells available in this study. The GSE132465 dataset contained annotated cell types from each sample. To eliminate cells of ambiguous identity, we determined diverse cell types in the GSE132465 dataset, and the “SingleR” package was used for initial clustering and cell type identification. T cells, B cells, epithelial cells, stromal cells, myeloid cells, and endothelial cells were especially recognized by t-distributed stochastic neighbor embedding (t-SNE) in the “Seurat” R package. The raw counts were then normalized and transformed to the natural log scale for subsequent analysis.

Cell culture and treatment

Human CRC cell lines (SW620, SW480, DLD-1, HT29, HCT116, HCT8, and RKO) and human normal intestinal epithelial cell line (NCM460) were purchased from the American Type Culture Collection (Manassas, VA, USA). SW620, SW480, DLD-1, HT29, HCT116, HCT8, and NCM460 were grown in complete RPMI- 1640 medium (VivaCell Biosciences, Shanghai, China). RKO was cultured in complete DMEM medium (VivaCell Biosciences). All media were supplemented with 10% FBS (Gibco, Gaithersburg, MD) and penicillin/streptomycin (NCM Biotech, Suzhou, China). Cells were grown in a humidified atmosphere at 37℃ with 5% CO2. FR019089 mimic, FR019089 inhibitor, FR019019 mimic, and FR019019 inhibitor were purchased from GenePharm (Shanghai, China). Transfection of mimic and inhibitor was conducted using GenMute siRNA Transfection Reagent (SignaGen Laboratories, Jinan, China).

RNA isolation, qRT-PCR, RNA-seq, and analysis

Total RNA from cells and tissues was isolated using the Trizol reagent (Invitrogen, Carlsbad, CA, USA). For qRT-PCR, total RNA (500 ng) was used to synthesize cDNA via the miRNA 1 st Strand cDNA Synthesis Kit (Vazyme Biotech, Nanjing, China, #MR101), and the miRNA Universal SYBR qPCR Master Mix (Vazyme Biotech, #MQ101) was used to conduct the amplification reactions according to protocols. The primers for qRT-PCR are listed in Additional file 1: Table S1.

For RNA-seq, total RNA was sequenced on an Illumina Novaseq™ 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the vendor’s recommended protocol. Raw Illumina sequencer output was converted to FASTQ format and was aligned to the Homo sapiens GRCh38 using the HISAT2 (https://ccb.jhu.edu/software/hisat2). The StringTie was used to estimate the expression levels of all transcripts by calculating FPKM. The EdgeR package was applied to analyze DEGs (absolute log2FC ≥ 1 and FDR < 0.05). GSEA was performed using “GSVA” packages [33]. Enrichment results with a criterion of qvalue < 0.05 were considered significant. Volcano plots were created using the “pheatmap” package. RNA sequencing data is available upon request.

CCK8 assay and colony formation assay

For the CCK8 assay, 2000 cells were plated in a 96-well microplate. After 0, 24, 48, and 72 h of proliferation, 10 μL of CCK8 (Biosharp, Hefei, China) was added to each well for 1 h, and cell viability was detected at 450 nm. For the short-term colony formation assay, 5000 cells were seeded in 12-well plates. After 1 week of proliferation, cells were fixed in 4% methanol and stained with crystal violet to count cell number. Three independent experiments were performed, and ImageJ software was used for quantification.

Transwell migration and invasion assay

Transwell and matrigel chamber plates (Corning Costar, New York, USA) were utilized to measure cell motility and invasiveness. 150,000 cells for RKO and 70,000 cells for HCT8 were loaded into a transwell cultured with serum-free media on the upside of the membrane, and a complete medium with 10% FBS was added at the bottom of the insert. Cells were led to adhere and transmigrate through the membrane at 37 °C, 5% CO2. After 28 h for RKO and 17 h for HCT8, transmigrated cells on the lower surface of the membrane were fixed in 4% methanol and stained with crystal violet. Images of cells on the lower face of the filter were captured in fields at × 10 magnification. Independent experiments were performed in triplicate. ImageJ software was used for quantification.

Statistical analysis

The Wilcoxon rank-sum test, t-test, and Kruskal-Wallis test were used for quantitative data comparisons, and Fisher’s exact test and \(\chi^2\) test were used to assess the difference in qualitative data. To investigate the prognostic potential of candidate piRNAs, Kaplan–Meier analysis and log-rank test were used to estimate and compare OS and PFS rates among different subgroups stratified by piRNA expression (high versus low level). A log-rank test P < 0.05 was used to define differences in survival time. For wet experiments, statistical analyses were carried out using GraphPad Prism software (version 8.0). For two- or three-group comparisons, Student’s t test or one-way ANOVA followed by Dunnett’s test was used and the figure legends specified the statistical tests used for each data. Graphs represent mean ± standard deviation (SD) unless otherwise stated. All P-values were two-sided, and P < 0.05 was considered statistically significant.

Results

Overview of clinical features and genomic alterations of EOCRC

In the TCGA cohort, a group of CRC patients provided tumor tissues for multi-omics profiling (Fig. 1A), and a total of 69 (13.40%) EOCRC and 446 (86.60%) LOCRC patients with summarized demographic characteristics in Additional file 1: Table S2 were included in this study. Patients with EOCRC were more frequently initially diagnosed with stage III (40.58% vs. 27.80%) or stage IV (20.29% vs. 13.68%). The prevalence of low microsatellite instability (MSI-low) tumors was higher among patients with EOCRC (23.19% vs. 15.93%), whereas the difference in MSI-high status between EOCRC and LOCRC was not observed (13.04% vs. 13.90%). Next, mutational analysis revealed that APC, TP53, and KRAS were the most frequently mutated genes (Fig. 1B, C), consistent with prior studies [7, 10]. Although APC mutations were less common in EOCRC than in LOCRC (67% vs. 76%), their mutation types and distribution differed markedly between subgroups (Additional file 2: Figure S1). CNV analysis showed amplifications in chromosomes 8, 13, and 17 in both groups, while deletions in chromosomes 6, 15, and 19 were more prevalent in EOCRC (Fig. 1D, E). Collectively, these results highlight distinct genomic alterations and molecular characteristics in EOCRC.

Fig. 1
figure 1

Population features and genomic alterations of EOCRC and LOCRC. A A total of 69 EOCRC and 446 LOCRC patients are enrolled in the TCGA cohort. The study diagram in the TCGA cohort is comprised of a comparison in somatic mutation and CNVs, profiling of transcriptome and DNA methylation, selection of diagnostic or prognostic biomarkers, and clinical variation. B,C Somatic mutations of the top mutated genes among different subtypes in EOCRC (B) and LOCRC (C). Genes are ordered by the total mutation frequencies. D,E The distribution of CNVs features across chromosomes for EOCRC (D) and LOCRC (E). TCGA, The Cancer Genome Atlas; EOCRC, early-onset colorectal cancer; LOCRC, late-onset CRC; CNV, copy number variation

Shift of immune microenvironment and cellular metabolism in EOCRC

The immune system is widely recognized as a critical factor in CRC onset and progression [35], and the tumor immune microenvironment has emerged as an independent biomarker of immunotherapy efficiency. To elucidate the feature of the tumor microenvironment of EOCRC, we utilized the single sample gene set enrichment analysis (ssGSEA) to quantify the abundance of 22 immune cell types. The tumor region of EOCRC was significantly enriched with memory B cells, whereas decreasingly enriched in NK cells and M1 macrophages (Fig. 2A). The significant pathways included pathways related to cancer-associated pathways and immune response (Fig. 2B). Interestingly, further GSEA analysis revealed upregulated cancer-associated pathways, including chemokine signaling pathway, JAK-STAT signaling pathway, IL-17 signaling pathway, and bacterial-related pathways (FDR < 0.05; Fig. 2C, D and Additional file 1: Table S3). In contrast, the downregulated pathways were mainly enriched in metabolic pathways, including ascorbate and aldarate metabolism, glycine, serine, and threonine metabolism, oxidative phosphorylation, cholesterol metabolism, and glycosphingolipid biosynthesis pathways (FDR < 0.05; Fig. 2E, F and Additional file 1: Table S3). Altogether, these results indicate the existence of an immune-inflammatory microenvironment in the tumor region, which is consistent with previous studies [36]. To assess metabolic differences, we analyzed fatty acid oxidation and glycolysis, finding significantly elevated expression of associated genes in EOCRC (Fig. 2G, H), indicating a shift toward aerobic glycolysis and lipid metabolism that may drive aggressive tumor biology.

Fig. 2
figure 2

Enrichment pathways of EOCRC. A The cell proportions in tumor microenvironment of EOCRC and LOCRC groups. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001, significant statistical differences between the two subgroups were assessed using the Wilcoxon test. B Volcano plot indicates the comparison at transcriptomic levels between EOCRC and LOCRC through ssGSEA enrichment scores. Upregulated pathways in the EOCRC group are shown in red and downregulated pathways are shown in blue. C–F GSEA results for the upregulated (C,D) and downregulated (E,F) pathways in the EOCRC group and LOCRC group, respectively (FDR < 0.05). G,H Quantification of selected genes of the fatty acids pathway (G) and glycolytic pathway (H) in the TCGA cohort (log2 intensity levels). Each boxplot shows the median and interquartile range (25 th–75 th percentiles). *P < 0.05, **P < 0.01, ***P < 0.001, significance was determined using a Wilcoxon rank-sum test. EOCRC, early-onset colorectal cancer; LOCRC, late-onset CRC; ns, not significant; GSEA, gene-set enrichment analysis; NES, normalized enrichment score

Differentially expressed PIWIL1 associated with tumorigenesis of EOCRC

We then probed into the EOCRC-specific molecular events at both transcriptional and epigenetic levels. Transcriptomic analysis identified 168 DEGs in LOCRC versus EOCRC (40 upregulated and 128 downregulated; Fig. 3A). Pathway enrichment revealed the top 27 significant gene ontology (GO) terms (P < 0.05, Fig. 3B), including cell function and the nervous system pathways. More DEGs were analyzed using weighted gene expression data (81 upregulated and 180 downregulated; Additional file 1: Table S4), which contained all the molecules analyzed above. At the epigenetic level, DNA methylation profiling detected 241 significant DMCs (230 hypermethylated and 11 hypomethylated) in LOCRC compared to EOCRC (Fig. 3C). Genes harboring these DMCs were enriched in 12 GO pathways (P < 0.05, Fig. 3D), involving carboxylic acid, monocarboxylic acid, and organic acid transport.

Fig. 3
figure 3

Gene expression and epigenetics difference of EOCRC and LOCRC tumors. A Volcano plot illustrates the differentially expressed genes in LOCRC. Upregulated genes in LOCRC and EOCRC patients are shown in red and blue, respectively. B Gene Ontology enrichment analysis for the differential expression genes. C Volcano plot indicates differentially methylated CpGs between EOCRC and LOCRC. Up-methylated CpGs in LOCRC and EOCRC patients are shown in red and blue, respectively. D Gene Ontology enrichment analysis for the CpG-mapped genes. EOCRC, early-onset colorectal cancer; LOCRC, late-onset CRC

The role of DNA methylation in gene regulation is multifaceted, but methylation in promoters is often considered the hallmark of gene repression [37]. Through integrated analysis of transcriptome and DNA methylation data, we consistently identified six key modules: PIWIL1, pyroglutamylated RFamide peptide receptor (QRFPR), carbohydrate sulfotransferase 8 (CHST8), GATA binding protein 4 (GATA4), membrane-spanning 4-domains A1 (MS4A1), and Wnt family member 3 A (WNT3A). Significant negative correlations were observed between promoter DMCs and gene expression levels (PIWIL1, r = − 0.54, P = 4.14e − 27; CHST8, r = − 0.20, P = 1.76e − 04; QRFPR, r = − 0.39, P = 5.83e − 14; GATA4, r = − 0.44, P = 1.07e − 17; MS4A1, r = − 0.13, P = 0.014; WNT3A, r = − 0.10, P = 0.065; Additional file 2: Figure S2 and Additional file 1: Table S5). These results highlight that in a variety of molecular processes, the promoter DNA hypermethylation, and hypomethylation regulate gene expression, which may be underneath the initiation and progression of EOCRC.

PIWIL1 belongs to the Piwi clade of the Argonaute family, mainly expressed in germ cells but absent in somatic tissues [38]. Overexpression of PIWIL1 is common to several tumor types (e.g., hepatocellular carcinoma (HCC) [39]), and its aberrant expression has been associated with tumorigenesis and poor prognosis [40]. Therefore, we first validated its expression pattern in CRC using TCGA and three external datasets (GSE39582, GSE17536, and GSE17537). As shown in Fig. 4A–D, the transcription level of PIWIL1 was significantly lower in EOCRC compared with LOCRC (P < 0.05). Further comparison of PIWIL1 expression in EOCRC, LOCRC, and healthy tissues using the TCGA data and clinical specimens (a total of 18 pairs of colorectum tumor and control specimens; baseline information of the CRC patients is shown in Additional file 1: Table S6) found that PIWIL1 was upregulated in both EOCRC and LOCRC tissues (Fig. 4E–G). Single-cell analysis of the tumor microenvironment revealed predominant PIWIL1 expression in epithelial cells, with significantly higher levels in LOCRC (Fig. 4H–M), particularly in CMS1 and CMS3 subtypes (Additional file 2: Figure S3). These findings suggest PIWIL1 as a key regulator of molecular heterogeneity between EOCRC and LOCRC.

Fig. 4
figure 4

Validation of differentially expressed PIWIL1 in different datasets. A–D The expression of PIWIL1 in EOCRC and LOCRC patients from the TCGA (A), GSE39582 (B), GSE17536 (C), and GSE17537 (D) datasets. E–G Differential expression analysis of PIWIL1 between normal and EOCRC (E) or LOCRC (F) tissues from the TCGA cohort and clinical specimens (G). Significance was determined using a Wilcoxon rank-sum test. H,I UMAP plots of unsupervised clustering result of tumor subclusters in EOCRC (H) and LOCRC (I) from the GSE132465 dataset. The tSNE plots show single cells of CRC tissues across six major cell types. J,K The expression pattern of PIWIL1 in single-cell level in EOCRC (J) and LOCRC (K) tissues. Cells expressing PIWIL1 are colored red in the tSNE plots using the normalized value. L–M Bubble plots of the average and percent expression of PIWIL1 in different cell subtypes of EOCRC (L) and LOCRC (M) tissues. EOCRC, early-onset colorectal cancer; LOCRC, late-onset CRC; tSNE, t-distributed stochastic neighbor embedding

FR019089 and FR019019 as anti-oncogenic piRNAs were involved in the progression and metastasis of EOCRC

Since PIWI proteins are central to piRNA biogenesis, we then examined the expression patterns of piRNAs in EOCRC and LOCRC samples. Three differentially expressed piRNAs (FR019019, P = 8.20e − 07; FR019089, P = 5.32e − 05; FR132045, P = 3.02e − 05; Fig. 5A–C) were identified as lower-expressed in EOCRC compared with LOCRC. These piRNAs can be retrieved in piRBase, piRNAdb, and RNAcentral (IDs and RNA sequences in Additional file 1: Table S7). To investigate their uncharacterized roles in CRC, we firstly examined their basal levels in SW620, SW480, DLD-1, HT29, HCT116, HCT8, RKO, and NCM460 cells. FR019089 and FR019019 showed relatively consistent expression patterns in the cell lines mentioned above (Additional file 2: Figure S4 A-C). Furthermore, results of qRT-PCR confirmed the under-expression of FR019089 and FR019019 in CRC tissues compared to the control specimens (Fig. 5D,E), indicating FR019089 and FR019019 as EOCRC-related anti-oncogenic piRNAs that may be involved in CRC progression and metastasis.

Fig. 5
figure 5

FR019089 and FR019019 are EOCRC-specific piRNAs. A–C Differentially expressed FR019089 (A), FR019019 (B), and FR132045 (C) between EOCRC and LOCRC patients from the TCGA datasets. D,E qRT-PCR detection of FR019089 (D) and FR019019 (E) in a subset of 18 cancerous and control tissues from patients of colorectal cancer. U6 was used as an internal control. Statistical analyses were performed by two-tailed Student’s t test. EOCRC, early-onset colorectal cancer; LOCRC, late-onset CRC

FR019089 and FR019019 inhibit migration and invasion in CRC cells

To further investigate the molecular function of FR019089 and FR019019, RKO and HCT8 cells were transfected with the mimics and inhibitors for FR019089 and FR019019. The knockdown and overexpression efficiencies were verified by qRT-PCR (Additional file 2: Figure S4 D-G). The CCK8 assay showed that FR019089 and FR019019 did not significantly affect cell multiplication, except for a modest increase in HCT8 proliferation upon FR019019 overexpression (Additional file 2: Figure S5). However, clone formation assays did not corroborate this proliferative effect, suggesting FR019019 may not sustain long-term growth advantage. We then interrogated whether they may regulate cell migration and invasion in CRC cells. As illustrated in Fig. 6A,B, the amount of migrated and invasive cells was remarkably cut down after mimic-FR019089 or -FR019019 treatment compared with the mimic-NC group. In contrast, the number of migrated and invasive cells in the inhibitor-FR019089 or -FR019019 group obviously rose (Fig. 6C,D). Collectively, our data show that FR019089 and FR019019 suppress CRC cell motility without altering proliferation, suggesting their role in inhibiting metastatic potential.

Fig. 6
figure 6

FR019089 and FR019019 attenuate migration and invasion in colorectal cancer cells and the prognostic significance of FR019089 in CRC patients. A,B RKO (A) and HCT8 (B) cells were transfected with either mimic-NC, -FR019089, or -FR019019, transwell assay was tested for the capability of migration and invasion. C,D RKO (C) and HCT8 (D) cells were transfected with either inhibitor-NC, -FR019089, or -FR019019, and transwell assay tested for the capability of migration and invasion. The experiment was repeated thrice independently with similar results. Bar = 100 μm. The quantitative results are presented as the mean ± SD. Statistical analyses were performed by one-way ANOVA with Dunnett’s test. **P < 0.001, ***P < 0.0001 vs. the mimic-NC or inhibitor-NC group. E–H The overall survival and progression-free survival analysis was performed by the Kaplan–Meier test and log-rank method in EOCRC (E,F) and LOCRC (G,H) patients, respectively. EOCRC, early-onset colorectal cancer; LOCRC, late-onset CRC

Enhanced cell invasion and migration promote cancer metastasis, leading to poor prognosis. We therefore assessed the association of FR019089 and FR019019 with patient survival. Kaplan–Meier analysis revealed that low FR019089 expression was significantly associated with shorter progression-free survival (PFS) and overall survival (OS) in EOCRC (Fig. 6E,F), but not in LOCRC patients (Fig. 6G,H; Additional file 2: Figure S6). Taken together, these results identify FR019089 as a promising prognostic biomarker specific to EOCRC.

Downstream regulated pathways of FR019089 and FR019019

Growing evidence has shown that piRNAs function by forming specific piRNA silencing complexes that bind a diverse spectrum of downstream target genes [41, 42]. We therefore performed RNA-seq to explore putative downstream targets of FR019089 and FR019019. RKO cells were transfected with mimic-FR019089, -FR019019, or negative control. Differential expression analysis revealed 136 upregulated and 13 downregulated genes under FR019089 overexpression, versus 36 upregulated and 7 downregulated genes for FR019019 (Fig. 7A,B). GSEA analysis indicated that mimic-FR019089 transfection promoted the epithelial cell apoptotic process, the positive regulation of programmed cell death, and negative regulation of stem cell proliferation (Fig. 7C and Additional file 1: Table S8), while FR019019 overexpression induced the epithelial cell apoptotic process and the secondary metabolic process, and negative regulation of canonical Wnt signaling pathway (Fig. 7 and Additional file 1: Table S9). Several downstream regulatory molecules (TNF, RIPK3, GATA3, and PPARG) of FR019089 and FR019019 were selected for verification by qRT-PCR. TNF and RIPK3 are involved in triggering programmed cell death, whereas GATA3 and PPARG are strongly associated with suppressing cell proliferation [43,44,45]. The results of qRT-PCR confirmed that high levels of FR019089 and FR019019 promoted the expression of TNF, RIPK3, GATA3, and PPARG (Additional file 2: Figure S7). These data suggest that FR019089 and FR019019 share partially overlapping functions, with enriched pathways linked to reduced cancer aggressiveness [46,47,48].

Fig. 7
figure 7

Identification of downstream pathways of FR019089 and FR019019. RKO cells were treated with mimic-FR019089, -FR019019, or -negative control (n = 4 for each group). A,B Volcano plots illustrate the differentially expressed genes between mimic-negative control with -FR019089 (A) or -FR019019 (B). Upregulated genes are shown in red and downregulated genes are shown in blue. C,D Results of gene-set enrichment analysis for significant pathways related to aggressiveness of colorectal cells in mimic-FR019089 (C) and mimic-FR019019 (D) group (qvalue < 0.05)

Discussion

Our findings utilized multi-omics integration to unveil genomic alterations, transcriptomic and epigenomic characteristics, as well as tumor microenvironment and metabolic pathways specific to EOCRC. Notably, integrative analysis of transcriptomic and DNA methylation profiles identified six candidate genes potentially implicated in EOCRC pathogenesis. We observed specific downregulation of PIWIL1, FR019089, and FR019019 in EOCRC, which correlated with increased migration and invasion, as well as poorer survival. Furthermore, in vitro experiments demonstrated the ability of FR019089 and FR019019 to inhibit cell migration and invasion, and subsequent RNA-seq confirmed the downstream anti-oncogenic signaling pathways mediated by high expression of FR019089 and FR019019. These findings provide novel insights into the PIWI/piRNA regulatory axis underlying the aggressive clinical phenotype of EOCRC, and we highlight that FR019089 has potential application in personalized treatment strategies.

In contrast to the declining trend of CRC-related mortality in adults 50 years or older, the overall incidence and mortality of EOCRC have been steadily rising [2, 12]. Emerging evidence supports that EOCRC represents a distinct clinical entity from LOCRC in terms of genetic predisposition, pathological characteristics, and tumor biology. Familial, hereditary, and sporadic EOCRC are three subgroups of EOCRC. Family history and hereditary conditions account for about 30% of EOCRC [49]. Familial adenomatous polyposis and Lynch syndrome are typical ones. For sporadic EOCRC, both environmental factors and genetic susceptibility are responsible for its occurrence, but the related germline genetic variation is unknown. Currently, Wang et al. [41] identified 49 independent genetic loci significantly associated with EOCRC risk, and Alkhateeb et al. [50] found recurring missense mutations in several genes (AXIN2, ALK, CDKN1A, MAP3K1, ROS1, EPCAM, KDM5A, and AURKA) in over 50% of Canadian samples in addition to previously known colorectal cancer-associated genetic variants. Moreover, Laskar et al. [51] sought out two novel risk loci (1p34.1 and 4p15.33) for EOCRC and found novel evidence of probable causal associations for modifiable lifestyle with increased EOCRC risk. EOCRC patients exhibit several common clinic-pathological features, including left colon localization of the primary tumor, more advanced stage at diagnosis, poorer tumor cell differentiation, higher prevalence of signet ring cell histology, and MSI-high due to germline mutations in the DNA mismatch repair [52,53,54,55]. Although a lower percentage of patients were considered to have EOCRC in our study (13.40%), we observed more aggressive features of EOCRC, including a higher frequency of stage III and IV at diagnosis and more prevalence of microsatellite instability tumors. These features might be the consequence of a lack of systematic screening in younger patients and more aggressive intrinsic tumor biology [55].

Migration and invasion are two facets of tumor malignancy. Apart from genetic selection, several traits can explain the enhanced aggressiveness of EOCRC. The metabolic shift of tumors to aerobic glycolysis (the Warburg effect) is a well-established hallmark of cancer, and metastatic cells have been shown to increase the uptake, synthesis, and utilization of lipids as a fuel source [56]. Our findings corroborated this paradigm, revealing significant upregulation of both fatty acid oxidation and glycolytic pathways in EOCRC. Moreover, epigenetic regulation, encompassing DNA methylation, histone modifications, and non-coding RNAs, provides a versatile mechanism for dynamic adaptation to microenvironmental inputs [56]. In this study, six EOCRC-specific molecules (PIWIL1, QRFPR, GATA4, CHST8, MS4A1, and WNT3A) exhibited DNA methylation-mediated regulation, consistent with prior reports [57,58,59,60]. These molecules likely participate in tumor microenvironment modulation and metabolic reprogramming [61,62,63,64]. For instance, PIWIL1 overexpression has been reported to accelerate the growth of HCC tumors through increasing oxygen consumption and energy production via fatty acid metabolism; additionally, PIWIL1-overexpressing HCC cells can attach myeloid-derived suppressor cells (MDSCs) into the tumor microenvironment, which in turn initiate the expression of immunosuppressive cytokine IL10 [39]. However, the precise mechanisms by which PIWIL1 orchestrates metabolic rewiring and immune microenvironment remodeling in CRC remain to be fully elucidated.

PIWIL1 is a member of the PIWI protein family [38]. We observed that PIWIL1 is expressed predominantly in epithelial cells and especially in CMS1 and CMS3. There are marked differences in the intrinsic biological underpinnings between four consensus molecular subtypes: CMS1, hypermutated, microsatellite unstable, strong immune activation; CMS2, epithelial, chromosomally unstable, marked WNT, and MYC signaling activation; CMS3, epithelial, evident metabolic dysregulation; and CMS4, prominent TGF-β activation, stromal invasion, and angiogenesis [65]. Analysis of the CRC tumor microenvironment and cellular metabolism revealed the upregulation of chemokine signaling, antigen processing and presentation, and IL-17 signaling pathways in both EOCRC and LOCRC; meanwhile, a remarkable shift in energy metabolism happened in CRC tissues. Importantly, these changes correlated with elevated PIWIL1 expression in CMS1/CMS3 subtypes, identifying it as a potential therapeutic target.

PIWI proteins and their mediated epigenetic modifications play an essential role in colorectal carcinogenesis and progression [66]. Our study found an increased level of PIWIL1 in CRC tissues relative to normal colonic mucosa. Although PIWIL1 has been reported to drive cancer metastasis in piRNA-independent mechanisms [67], it also functionally cooperates with piRNAs to shape the clinicopathological characteristics of CRC [16]. Emerging evidence suggests the clinical significance of piRNAs in CRC and that genetic variations in piRNAs may modulate CRC susceptibility [14, 68]. However, despite these advances, there is no research on piRNAs specifically related to EOCRC and their association with clinical and pathological phenotypes within this context. Here, novel EOCRC-specific piRNAs (FR019089 and FR019019) were identified, both of which displayed anti-migration and anti-invasion effects, and FR019089 exhibited the potential as a therapeutic molecule. Collectively, our findings contribute to a deeper understanding of the mechanisms of piRNAs in EOCRC tumorigenesis and further underscore the translational value of piRNAs.

However, this study still has several limitations. Firstly, although we endeavored to maximize the inclusion of EOCRC cases to enhance analytical rigor, these results could represent part of the omics characteristic of EOCRC patients due to the older diagnosed age of CRC patients in TCGA and GEO datasets. Secondly, although not reaching statistical significance, we observed potentially important molecular alterations in EOCRC, including lower mutation frequency and copy number variations (deletions/amplifications). Thirdly, the study lacked large-sample and cross-ethnicity external validation of EOCRC in independent datasets, as well as molecular mechanism investigations into PIWIL1/FR019089 and PIWIL1/FR019019, and research on DNA modification. Finally, no significant difference in EOCRC vs. LOCRC for target piRNAs was observed due to the small sample size of CRC patients. Therefore, more in-depth profiling focusing on EOCRC patients and further exploration of mechanisms elaborating EOCRC pathophysiology are of great necessity in future research.

Conclusions

Our findings characterize the molecular landscape of EOCRC and implicate FR019089 and FR019019 of piRNA as potential modulators influencing the metastatic and invasive characteristics of EOCRC. FR019089 exhibits potential as a predictive biomarker for EOCRC prognosis. To the best of our knowledge, these data provide the first evidence of the role of piRNAs in EOCRC, potentially supporting their use as prognostic biomarkers in EOCRC patients. These findings could be instrumental in advancing precision cancer therapy.

Data availability

The data that support the findings of this study are publicly available as described in the Method section or available from the corresponding author upon reasonable request. Analysis codes used to generate the results were deposited to the GitHub repository (https://github.com/yunyun-1996/EOCRC.git).

Abbreviations

CRC:

Colorectal cancer

EOCRC:

Early-onset colorectal cancer

LOCRC:

Late-onset colorectal cancer

ncRNA:

Non-coding RNA

piRNA:

PIWI-interacting RNA

PIWIL1:

Piwi Like RNA-Mediated Gene Silencing 1

MSI:

Microsatellite instability

ssGSEA:

Single sample gene set enrichment analysis

DEGs:

Differentially expressed genes

DMCs:

Differentially methylated CpGs

QRFPR:

Pyroglutamylated RFamide peptide receptor

CHST8:

Carbohydrate sulfotransferase 8

GATA4:

GATA binding protein 4

MS4A1:

Membrane spanning 4-domains A1

WNT3A:

Wnt family member 3A

NES:

Normalized enrichment score

PFS:

Progression-free survival

OS:

Overall survival

GEO:

Gene Expression Omnibus

FC:

Fold change

FDR:

False discovery rate

KEGG:

Kyoto Encyclopedia of Genes and Genomes

t-SNE:

t-Distributed stochastic neighbor embedding

UMI:

Unique molecular identifier

SD:

Standard deviation

GO:

Gene ontology

IPTW:

Inverse probability of treatment weighting

References

  1. Siegel RL, Miller KD, Wagle NS, Jemal A. Cancer statistics, 2023. Ca: Cancer J Clin. 2023;73(1):17–48. https://doiorg.publicaciones.saludcastillayleon.es/10.3322/caac.21763.

    Article  PubMed  Google Scholar 

  2. Akimoto N, Ugai T, Zhong R, Hamada T, Fujiyoshi K, Giannakis M, Wu K, Cao Y, Ng K, Ogino S. Rising incidence of early-onset colorectal cancer - a call to action. Nat Rev Clin Oncol. 2021;18(4):230–43. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41571-020-00445-1.

    Article  PubMed  Google Scholar 

  3. Giannakis M, Ng K. A common cancer at an uncommon age. Science (American Association for the Advancement of Science). 2023;379(6637):1088–90. https://doiorg.publicaciones.saludcastillayleon.es/10.1126/science.ade7114.

    Article  CAS  Google Scholar 

  4. Zhao J, Xu L, Sun J, Song M, Wang L, Yuan S, Zhu Y, Wan Z, Larsson S, Tsilidis K, et al. Global trends in incidence, death, burden and risk factors of early-onset cancer from 1990 to 2019. Bmj Oncology. 2023;2(1):e49. https://doiorg.publicaciones.saludcastillayleon.es/10.1136/bmjonc-2023-000049.

    Article  Google Scholar 

  5. Mauri G, Patelli G, Crisafulli G, Siena S, Bardelli A. Tumor, “age” in early-onset colorectal cancer. Cell. 2025;188(3):589–93. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.cell.2024.12.003.

    Article  CAS  PubMed  Google Scholar 

  6. Atikukke G, Alkhateeb A, Porter L, Fifield B, Cavallo-Medved D, Facca J, Elfiki T, Elkeilani A, Rueda L, Misra S. P-370 Comprehensive targeted genomic profiling and comparative genomic analysis to identify molecular mechanisms driving cancer progression in young-onset sporadic colorectal cancer. Ann Oncol. 2020;31:S209–10. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.annonc.2020.04.452.

    Article  Google Scholar 

  7. Willauer AN, Liu Y, Pereira AAL, Lam M, Morris JS, Raghav KPS, Morris VK, Menter D, Broaddus R, Meric Bernstam F, et al. Clinical and molecular characterization of early-onset colorectal cancer. Cancer-Am Cancer Soc. 2019;125(12):2002–10. https://doiorg.publicaciones.saludcastillayleon.es/10.1002/cncr.31994.

    Article  CAS  Google Scholar 

  8. Liu PH, Wu K, Ng K, Zauber AG, Nguyen LH, Song M, He X, Fuchs CS, Ogino S, Willett WC, et al. Association of obesity with risk of early-onset colorectal cancer among women. Jama Oncol. 2019;5(1):37–44. https://doiorg.publicaciones.saludcastillayleon.es/10.1001/jamaoncol.2018.4280.

    Article  PubMed  Google Scholar 

  9. Nguyen LH, Liu PH, Zheng X, Keum N, Zong X, Li X, Wu K, Fuchs CS, Ogino S, Ng K, et al. Sedentary behaviors, TV viewing time, and risk of young-onset colorectal cancer. Jnci Cancer Spect. 2018;2(4):pky73. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/jncics/pky073.

    Article  Google Scholar 

  10. Lieu CH, Golemis EA, Serebriiskii IG, Newberg J, Hemmerich A, Connelly C, Messersmith WA, Eng C, Eckhardt SG, Frampton G, et al. Comprehensive genomic landscapes in early and later onset colorectal cancer. Clin Cancer Res. 2019;25(19):5852–8. https://doiorg.publicaciones.saludcastillayleon.es/10.1158/1078-0432.CCR-19-0899.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Koh G, Degasperi A, Zou X, Momen S, Nik-Zainal S. Mutational signatures: emerging concepts, caveats and clinical applications. Nat Rev Cancer. 2021;21(10):619–37. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41568-021-00377-7.

    Article  CAS  PubMed  Google Scholar 

  12. Nakamura K, Hernandez G, Sharma GG, Wada Y, Banwait JK, Gonzalez N, Perea J, Balaguer F, Takamaru H, Saito Y, et al. A liquid biopsy signature for the detection of patients with early-onset colorectal cancer. Gastroenterology. 2022;163(5):1242–51. https://doiorg.publicaciones.saludcastillayleon.es/10.1053/j.gastro.2022.06.089.

    Article  CAS  PubMed  Google Scholar 

  13. Furuhashi S, Bustos MA, Mizuno S, Ryu S, Naeini Y, Bilchik AJ, Hoon D. Spatial profiling of cancer-associated fibroblasts of sporadic early onset colon cancer microenvironment. Npj Precis Oncol. 2023;7(1):118. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41698-023-00474-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Weng W, Liu N, Toiyama Y, Kusunoki M, Nagasaka T, Fujiwara T, Wei Q, Qin H, Lin H, Ma Y, et al. Novel evidence for a PIWI-interacting RNA (piRNA) as an oncogenic mediator of disease progression, and a potential prognostic biomarker in colorectal cancer. Mol Cancer. 2018;17(1):16. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12943-018-0767-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Chen S, Ben S, Xin J, Li S, Zheng R, Wang H, Fan L, Du M, Zhang Z, Wang M. The biogenesis and biological function of PIWI-interacting RNA in cancer. J Hematol Oncol. 2021;14(1):1–93. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13045-021-01104-3.

    Article  CAS  Google Scholar 

  16. Sellitto A, Geles K, D’Agostino Y, Conte M, Alexandrova E, Rocco D, Nassa G, Giurato G, Tarallo R, Weisz A, et al. Molecular and functional characterization of the somatic PIWIL1/piRNA pathway in colorectal cancer cells. Cells (Basel, Switzerland). 2019;8(11):1390. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/cells8111390.

    Article  CAS  Google Scholar 

  17. Sun R, Gao CL, Li DH, Li BJ, Ding YH. Expression status of PIWIL1 as a prognostic marker of colorectal cancer. Dis Markers. 2017;2017:1204937. https://doiorg.publicaciones.saludcastillayleon.es/10.1155/2017/1204937.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Network CGA. Comprehensive molecular characterization of human colon and rectal cancer. Nature. 2012;487(7407):330–7. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nature11252.

    Article  CAS  Google Scholar 

  19. Marisa L, de Reynies A, Duval A, Selves J, Gaub MP, Vescovo L, Etienne-Grimaldi MC, Schiappa R, Guenot D, Ayadi M, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. Plos Med. 2013;10(5):e1001453. https://doiorg.publicaciones.saludcastillayleon.es/10.1371/journal.pmed.1001453.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Marisa L, de Reynies A, Duval A, Selves J, Gaub MP, Vescovo L, Etienne-Grimaldi MC, Schiappa R, Guenot D and Ayadi M, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. GEO. 2013. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39582.

  21. Smith JJ, Deane NG, Wu F, Merchant NB, Zhang B, Jiang A, Lu P, Johnson JC, Schmidt C, Bailey CE, et al. Experimentally derived metastasis gene expression profile predicts recurrence and death in patients with colon cancer. Gastroenterology. 2010;138(3):958–68. https://doiorg.publicaciones.saludcastillayleon.es/10.1053/j.gastro.2009.11.005.

    Article  CAS  PubMed  Google Scholar 

  22. Smith JJ, Deane NG, Wu F, Merchant NB, Zhang B, Jiang A, Lu P, Johnson JC, Schmidt C, Bailey CE, et al. Experimentally derived metastasis gene expression profile predicts recurrence and death in patients with colon cancer. GEO. 2010. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17536.

  23. Smith JJ, Deane NG, Wu F, Merchant NB, Zhang B, Jiang A, Lu P, Johnson JC, Schmidt C and Bailey CE, et al. Experimentally derived metastasis gene expression profile predicts recurrence and death in patients with colon cancer. GEO. 2010. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17537.

  24. Lee HO, Hong Y, Etlioglu HE, Cho YB, Pomella V, Van den Bosch B, Vanhecke J, Verbandt S, Hong H, Min JW, et al. Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer. Nat Genet. 2020;52(6):594–603. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41588-020-0636-z.

    Article  CAS  PubMed  Google Scholar 

  25. Lee HO, Hong Y, Etlioglu HE, Cho YB, Pomella V, Van den Bosch B, Vanhecke J, Verbandt S, Hong H and Min JW, et al. Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer. GEO. 2020. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132465.

  26. Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12(4):R41. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/gb-2011-12-4-r41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Kool R, Dragomir A, Kulkarni GS, Marcq G, Breau RH, Kim M, Busca I, Abdi H, Dawidek M, Uy M, et al. Benefit of neoadjuvant cisplatin-based chemotherapy for invasive bladder cancer patients treated with radiation-based therapy in a real-world setting: an inverse probability treatment weighted analysis. Eur Urol Oncol. 2024;7(6):1350–7. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.euo.2024.01.014.

    Article  PubMed  Google Scholar 

  28. Yu G, Wang L, Han Y, He Q. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics (Larchmont, NY). 2012;16(5):284–7. https://doiorg.publicaciones.saludcastillayleon.es/10.1089/omi.2011.0118.

    Article  CAS  Google Scholar 

  29. Sandoval J, Heyn H, Moran S, Serra-Musach J, Pujana MA, Bibikova M, Esteller M. Validation of a DNA methylation microarray for 450,000 CpG sites in the human genome. Epigenetics-Us. 2011;6(6):692–702. https://doiorg.publicaciones.saludcastillayleon.es/10.4161/epi.6.6.16196.

    Article  CAS  Google Scholar 

  30. Bibikova M, Barnes B, Tsan C, Ho V, Klotzle B, Le JM, Delano D, Zhang L, Schroth GP, Gunderson KL, et al. High density DNA methylation array with single CpG site resolution. Genomics. 2011;98(4):288–95. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.ygeno.2011.07.007.

    Article  CAS  PubMed  Google Scholar 

  31. Tian Y, Morris TJ, Webster AP, Yang Z, Beck S, Feber A, Teschendorff AE. ChAMP: updated methylation analysis pipeline for Illumina BeadChips. Bioinformatics. 2017;33(24):3982–4. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/btx513.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Martinez VD, Vucic EA, Thu KL, Hubaux R, Enfield KS, Pikor LA, Becker-Santos DD, Brown CJ, Lam S, Lam WL. Unique somatic and malignant expression patterns implicate PIWI-interacting RNAs in cancer-type specific biology. Sci Rep-Uk. 2015;5:10423. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/srep10423.

    Article  Google Scholar 

  33. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/1471-2105-14-7.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36(5):411–20. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nbt.4096.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Fletcher R, Wang YJ, Schoen RE, Finn OJ, Yu J, Zhang L. Colorectal cancer prevention: Immune modulation taking the stage. Bba-Rev Cancer. 2018;1869(2):138–48. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.bbcan.2017.12.002.

    Article  CAS  Google Scholar 

  36. Du M, Gu D, Xin J, Peters U, Song M, Cai G, Li S, Ben S, Meng Y, Chu H, et al. Integrated multi-omics approach to distinct molecular characterization and classification of early-onset colorectal cancer. Cell Rep Med. 2023;4(3):100974. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.xcrm.2023.100974.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Baylin SB. DNA methylation and gene silencing in cancer. Nat Clin Pract Oncol. 2005;2(S1):S4–11. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/ncponc0354.

    Article  CAS  PubMed  Google Scholar 

  38. Peng JC, Valouev A, Liu N, Lin H. Piwi maintains germline stem cells and oogenesis in Drosophila through negative regulation of Polycomb group proteins. Nat Genet. 2016;48(3):283–91. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/ng.3486.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Wang N, Tan HY, Lu Y, Chan YT, Wang D, Guo W, Xu Y, Zhang C, Chen F, Tang G, et al. PIWIL1 governs the crosstalk of cancer cell metabolism and immunosuppressive microenvironment in hepatocellular carcinoma. Signal Transduct Tar. 2021;6(1):86. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41392-021-00485-8.

    Article  CAS  Google Scholar 

  40. Dong P, Xiong Y, Konno Y, Ihira K, Xu D, Kobayashi N, Yue J, Watari H. Critical roles of PIWIL1 in human tumors: expression, functions, mechanisms, and potential clinical implications. Front Cell Dev Biol. 2021;9:656993. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fcell.2021.656993.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Watanabe T, Lin H. Posttranscriptional regulation of gene expression by Piwi proteins and piRNAs. Mol Cell. 2014;56(1):18–27. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.molcel.2014.09.012.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Zhong F, Zhou N, Wu K, Guo Y, Tan W, Zhang H, Zhang X, Geng G, Pan T, Luo H, et al. A SnoRNA-derived piRNA interacts with human interleukin-4 pre-mRNA and induces its decay in nuclear exosomes. Nucleic Acids Res. 2015;43(21):10474–91. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/nar/gkv954.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Annibaldi A, Meier P. Checkpoints in TNF-induced cell death: implications in inflammation and cancer. Trends Mol Med. 2018;24(1):49–65. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.molmed.2017.11.002.

    Article  CAS  PubMed  Google Scholar 

  44. Lan T, Li H, Zhang D, Xu L, Liu H, Hao X, Yan X, Liao H, Chen X, Xie K, et al. KIAA1429 contributes to liver cancer progression through N6-methyladenosine-dependent post-transcriptional modification of GATA3. Mol Cancer. 2019;18(1):186. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12943-019-1106-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Savic D, Ramaker RC, Roberts BS, Dean EC, Burwell TC, Meadows SK, Cooper SJ, Garabedian MJ, Gertz J, Myers RM. Distinct gene regulatory programs define the inhibitory effects of liver X receptors and PPARG on cancer cell proliferation. Genome Med. 2016;8(1):74. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13073-016-0328-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Lin Z, Wan AH, Sun L, Liang H, Niu Y, Deng Y, Yan S, Wang QP, Bu X, Zhang X, et al. N6-methyladenosine demethylase FTO enhances chemo-resistance in colorectal cancer through SIVA1-mediated apoptosis. Mol Ther. 2023;31(2):517–34. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.ymthe.2022.10.012.

    Article  CAS  PubMed  Google Scholar 

  47. Prasad CP, Chaurasiya SK, Guilmain W, Andersson T. WNT5A signaling impairs breast cancer cell migration and invasion via mechanisms independent of the epithelial-mesenchymal transition. J Exp Clin Canc Res. 2016;35(1):144. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13046-016-0421-0.

    Article  CAS  Google Scholar 

  48. Zhao H, Ming T, Tang S, Ren S, Yang H, Liu M, Tao Q, Xu H. Wnt signaling in colorectal cancer: pathogenic role and therapeutic target. Mol Cancer. 2022;21(1):144. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12943-022-01616-7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Hofseth LJ, Hebert JR, Chanda A, Chen H, Love BL, Pena MM, Murphy EA, Sajish M, Sheth A, Buckhaults PJ, et al. Early-onset colorectal cancer: initial clues and current views. Nat Rev Gastro Hepat. 2020;17(6):352–64. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41575-019-0253-4.

    Article  Google Scholar 

  50. Alkhateeb A, Atikukke E, El Keilani A, Elfiki TA, Orlando J, Cavallo-Medved D, Rueda L, Fetz A, Ravid-Einy A, Porter LA, et al. Comparative genomic analysis to identify a signature of sporadic colorectal cancer development in young adults. J Clin Oncol. 2024;42(3_suppl):208. https://doiorg.publicaciones.saludcastillayleon.es/10.1200/JCO.2024.42.3_suppl.208.

    Article  Google Scholar 

  51. Laskar RS, Qu C, Huyghe JR, Harrison T, Hayes RB, Cao Y, Campbell PT, Steinfelder R, Talukdar FR, Brenner H, et al. Genome-wide association study and Mendelian randomization analyses provide insights into the causes of early-onset colorectal cancer. Ann Oncol. 2024. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.annonc.2024.02.008.

    Article  PubMed  Google Scholar 

  52. Di Leo M, Zuppardo RA, Puzzono M, Ditonno I, Mannucci A, Antoci G, Russo RA, Patricelli MG, Elmore U, Tamburini AM, et al. Risk factors and clinical characteristics of early-onset colorectal cancer vs. late-onset colorectal cancer: a case-case study. Eur J Gastroen Hepat. 2021;33(9):1153–60. https://doiorg.publicaciones.saludcastillayleon.es/10.1097/MEG.0000000000002000.

    Article  Google Scholar 

  53. Archambault AN, Su Y, Jeon J, Thomas M, Lin Y, Conti DV, Win AK, Sakoda LC, Lansdorp-Vogelaar I, Peterse EFP, et al. Cumulative burden of colorectal cancer-associated genetic variants is more strongly associated with early-onset vs late-onset cancer. Gastroenterology. 2020;158(5):1274–86. https://doiorg.publicaciones.saludcastillayleon.es/10.1053/j.gastro.2019.12.012.

    Article  CAS  PubMed  Google Scholar 

  54. Siegel RL, Miller KD, Goding Sauer A, Fedewa SA, Butterly LF, Anderson JC, Cercek A, Smith RA, Jemal A. Colorectal cancer statistics, 2020. Ca: Cancer J Clin. 2020;70(3):145–64. https://doiorg.publicaciones.saludcastillayleon.es/10.3322/caac.21601.

    Article  PubMed  Google Scholar 

  55. Martinez-Perez D, Vinal D, Pena-Lopez J, Jimenez-Bou D, Ruiz-Gutierrez I, Martinez-Recio S, Alameda-Guijarro M, Rueda-Lara A, Martin-Montalvo G, Ghanem I, et al. Clinico-pathological features, outcomes and impacts of COVID-19 pandemic on patients with early-onset colorectal cancer: a single-institution experience. Cancers. 2023;15(17):4242. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/cancers15174242.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Gerstberger S, Jiang Q, Ganesh K. Metastasis. Cell. 2023;186(8):1564–79. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.cell.2023.03.003.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Chen Z, Yang HJ, Lin Q, Zhu MJ, Yu YY, He XY, Wan XP. Estrogen-ERalpha signaling and DNA hypomethylation co-regulate expression of stem cell protein PIWIL1 in ERalpha-positive endometrial cancer cells. Cell Commun Signal. 2020;18(1):84. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12964-020-00563-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Faenza M, Benincasa G, Docimo L, Nicoletti GF, Napoli C. Clinical epigenetics and restoring of metabolic health in severely obese patients undergoing batriatric and metabolic surgery. Updates Surg. 2022;74(2):431–8. https://doiorg.publicaciones.saludcastillayleon.es/10.1007/s13304-021-01162-9.

    Article  PubMed  Google Scholar 

  59. Santoro ML, Santos CM, Ota VK, Gadelha A, Stilhano RS, Diana MC, Silva PN, Spíndola LMN, Melaragno MI, Bressan RA, et al. Expression profile of neurotransmitter receptor and regulatory genes in the prefrontal cortex of spontaneously hypertensive rats: relevance to neuropsychiatric disorders. Psychiat Res. 2014;219(3):674–9. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.psychres.2014.05.034.

    Article  CAS  Google Scholar 

  60. Hellebrekers DM, Lentjes MH, van den Bosch SM, Melotte V, Wouters KA, Daenen KL, Smits KM, Akiyama Y, Yuasa Y, Sanduleanu S, et al. GATA4 and GATA5 are potential tumor suppressors and biomarkers in colorectal cancer. Clin Cancer Res. 2009;15(12):3990–7. https://doiorg.publicaciones.saludcastillayleon.es/10.1158/1078-0432.CCR-09-0055.

    Article  CAS  PubMed  Google Scholar 

  61. Pavlasova G, Borsky M, Seda V, Cerna K, Osickova J, Doubek M, Mayer J, Calogero R, Trbusek M, Pospisilova S, et al. Ibrutinib inhibits CD20 upregulation on CLL B cells mediated by the CXCR4/SDF-1 axis. Blood. 2016;128(12):1609–13. https://doiorg.publicaciones.saludcastillayleon.es/10.1182/blood-2016-04-709519.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Mudd TJ, Lu C, Klement JD, Liu K. MS4A1 expression and function in T cells in the colorectal cancer tumor microenvironment. Cell Immunol. 2021;360:104260. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.cellimm.2020.104260.

    Article  CAS  PubMed  Google Scholar 

  63. Chartrel N, Picot M, El MM, Arabo A, Berrahmoune H, Alexandre D, Maucotel J, Anouar Y, Prevost G. The neuropeptide 26RFa (QRFP) and its role in the regulation of energy homeostasis: a mini-review. Front Neurosci-Switz. 2016;10:549. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fnins.2016.00549.

    Article  Google Scholar 

  64. Ye D, Liu H, Zhao G, Chen A, Jiang Y, Hu Y, Liu D, Xie N, Liang W, Chen X, et al. LncGMDS-AS1 promotes the tumorigenesis of colorectal cancer through HuR-STAT3/Wnt axis. Cell Death Dis. 2023;14(2):165. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41419-023-05700-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Guinney J, Dienstmann R, Wang X, de Reyniès A, Schlicker A, Soneson C, Marisa L, Roepman P, Nyamundanda G, Angelino P, et al. The consensus molecular subtypes of colorectal cancer. Nat Med. 2015;21(11):1350–6. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nm.3967.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Lin Y, Zheng J, Lin D. PIWI-interacting RNAs in human cancer. Semin Cancer Biol. 2021;75:15–28. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.semcancer.2020.08.012.

    Article  CAS  PubMed  Google Scholar 

  67. Li F, Yuan P, Rao M, Jin CH, Tang W, Rong YF, Hu YP, Zhang F, Wei T, Yin Q, et al. piRNA-independent function of PIWIL1 as a co-activator for anaphase promoting complex/cyclosome to drive pancreatic cancer metastasis. Nat Cell Biol. 2020;22(4):425–38. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41556-020-0486-z.

    Article  CAS  PubMed  Google Scholar 

  68. Chu H, Xia L, Qiu X, Gu D, Zhu L, Jin J, Hui G, Hua Q, Du M, Tong N, et al. Genetic variants in noncoding PIWI-interacting RNA and colorectal cancer risk. Cancer-Am Cancer Soc. 2015;121(12):2044–52. https://doiorg.publicaciones.saludcastillayleon.es/10.1002/cncr.29314.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We would like to thank Nan Zhou from the Core Facilities, Zhejiang University School of Medicine for their technical support.

Funding

This research was supported by the National Natural Science Foundation of China (82204019), the Zhejiang Provincial Clinical Research Center for CANCER (2022E50008, 2024ZY01056), the Noncommunicable Chronic Diseases-National Science and Technology Major Project (2024ZD0520100), and the Key R&D Program of Zhejiang (2023C03049).

Author information

Authors and Affiliations

Authors

Contributions

Conceptualization: XL, KD, HZ; Methodology: SZ, LY, JZ, QX, JS, LW, YZ, YL; Investigation: SZ, LY, JZ; Visualization: SZ, LY; Supervision: XL, KD, HZ, ET, MD; Writing—original draft: SZ, LY, JZ; Writing—review & editing: XL, KD, HZ. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Honghe Zhang, Kefeng Ding or Xue Li.

Ethics declarations

Ethics approval and consent to participate

The acquisition of all tissue samples utilized in this research was conducted following the approval granted by the Biomedical Ethics Committee of Zhejiang University (ZGL202306-3), with each patient providing written informed consent.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

12916_2025_4074_MOESM1_ESM.xlsx

Additional file 1: Table S1. Sequences of qPCR primers. Table S2. Clinical and tumor characteristics of patients with early-onset and late-onset colorectal cancer. Table S3. GSEA analysis of the different expression genes. Table S4. Differentially expressed genes identified in the TCGA database after using inverse probability of treatment weighting method to balance baseline covariates associated with EOCRC and LOCRC. Table S5. Consistent differential methylation sites and transcription genes identified in the TCGA database. Table S6. Baseline information of CRC patients. Table S7. RNA sequences of differentially expressed piRNAs and their IDs in piRBase, piRNAdb and RNAcentral. Table S8. GSEA analysis of RNA-seq results (mimic-FR019089 vs mimic-NC). Table S9. GSEA analysis of RNA-seq results (mimic-FR019019 vs mimic-NC).

12916_2025_4074_MOESM2_ESM.docx

Additional file 2: Figure S1. Mutation in APC gene in CRC subgroups. Figure S2. The Spearman correlation of DNA methylation with gene expression. Figure S3. Expression pattern of PIWIL1 in different cell types and CMSs of age-of-onset colorectal cancer. Figure S4. Basal expression of FR019089, FR019019, and FR132045 in colon cell lines and efficiency verification of mimic and inhibitor. Figure S5. Effects of FR019089 and FR019019 on cellular proliferation. Figure S6. The prognostic significance of FR019019 in CRC patients. Figure S7. FR019089 and FR019019 promote the expression of TNF, RIPK3, GATA3, and PPARG.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhou, S., Yu, L., Zhao, J. et al. Integration of multi-omics data to unveil the molecular landscape and role of piRNAs in early-onset colorectal cancer. BMC Med 23, 250 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12916-025-04074-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12916-025-04074-2

Keywords