Colorectal Adenoma Subtypes Exhibit Signature Molecular Profiles: Unique Insights into the Microenvironment of Advanced Precancerous Lesions for Early Detection Applications

1. Introduction

Colorectal cancer (CRC) encompasses a heterogeneous group of malignancies arising from the epithelial cells lining the colon and rectum, representing a substantial burden on public health, with over 1.9 million new cases and 935,000 deaths in 2020 (1). CRC is the third most diagnosed cancer and the second leading cause of cancer-related mortality globally. Understanding the intricate molecular mechanisms driving the initiation, progression, and metastasis of CRC is crucial for early detection and the development of targeted therapeutic interventions to improve clinical outcomes (2).

The etiology of CRC is multifactorial, with both genetic and environmental factors playing pivotal roles. Whilst hereditary syndromes, such as familial adenomatous polyposis (FAP) and Lynch syndrome account for a subset of cases (3), most of CRC cases occur sporadically and are strongly influenced by lifestyle factors including diet, physical inactivity, and smoking (4). Alarmingly, recent years have shown a dramatic increase in patients with early-onset CRC (5).

The adenoma-carcinoma sequence, proposed by Vogelstein (6), delineates a stepwise progression of CRC from small dysplastic lesions via Advanced Precancerous Lesions (APLs) to the development of invasive carcinomas. Genetic alterations play pivotal roles in driving this progression, particularly frame shift single nucleotide mutations as in Kirsten rat sarcoma viral oncogene homolog (KRAS) (7) and tumor protein p53 (TP53) (8) genes or both copy number variations (CNVs) and single nucleotide mutations as in adenomatous polyposis coli (APC) (9), among others (10). Dysregulation of key signaling pathways, such as the Wnt/β-catenin, Ras-MAPK, and PI3K-Akt, further contributes to uncontrolled cellular proliferation characteristic of CRC (11). Notably, CRC and adenomas also exhibit significant levels of apoptosis, indicating a complex balance between cellular proliferation and cell death. This interplay, particularly in tumors with high-turnover phenotypes, may influence both prognosis and susceptibility to chemotherapeutic interventions, as suggested by emerging evidence (12) (13).

APLs can be broadly classified into several subtypes based on their histological features: tubular (TA), Villous/Tubulovillous (VA/TVA), and serrated adenomas (SA). Each subtype varies in its histopathological characteristics, genetic mutations, and risk of progression to colorectal cancer (14). TAs, the most common type, are characterized by a predominantly tubular glandular structure. Histologically, they exhibit elongated, tube-like glands lined by dysplastic epithelial cells. Tubular adenomas are generally smaller and have a lower risk of progressing to colorectal cancer than other subtypes; however, the risk increases with the size and degree of dysplasia of the adenoma (15).

Villous adenomas are less common but are associated with a higher risk of malignant transformation. These adenomas have a frond-like or villous architecture, with finger-like projections extending into the lumen of the colon. VAs tend to be larger and more often display high-grade dysplasia and display a greater propensity for transformation to invasive cancer (16).

Tubulovillous adenomas exhibit a combination of both tubular and villous features, comprising 25-75% villous structures. This hybrid nature places them at an intermediate risk of malignant transformation between that of tubular and villous adenomas respectively. The proportion of the villous component within these adenomas is correlated with risk of cancerous transformation (17).

SAs, encompassing traditional serrated adenomas (TSAs) and sessile serrated polyps/lesions (SSPs/SSLs), represent a distinct category characterized by a serrated or saw-toothed epithelial architecture. These adenomas are associated with the serrated pathway of colorectal carcinogenesis, which involves distinct genetic mutations such as BRAF or KRAS (18)(19). Serrated adenomas, particularly of the sessile type, are recognized for their potential to progress to colorectal cancer through the serrated neoplasia pathway, a process that is often marked by CpG island methylator phenotype (CIMP) and microsatellite instability (MSI) (20).

The CIMP phenotype is characterized by widespread hypermethylation of CpG islands in the promoter regions of genes, leading to their silencing. This epigenetic alteration is particularly associated with specific genetic mutations, such as those in BRAF, and is more prevalent in right-sided colon cancers (21). Microsatellite instability, on the other hand, arises from defects in the DNA mismatch repair (MMR) system, resulting in insertion and deletion mutations within the genome's microsatellite regions (22). Remarkably, MSI-high status is a hallmark of Lynch syndrome-associated cancers and a subset of sporadic CRCs (23). Both CIMP and MSI phenotypes contribute to the heterogeneity of tumor behavior, prognosis, and response to therapy in colorectal cancer (24).

Beyond the genetic landscape, emerging evidence implicates the tumor microenvironment, immune system, and epigenetic modifications in CRC pathogenesis (25). Altogether, the heterogeneity observed in CRC tumors underscores the need for multiple approaches that consider the diverse molecular profiles exhibited by individual patients (26).

The primary objective of this study was to delineate a molecular signature that encompasses differentially methylated regions (DMRs), CNVs, and genomic mutations within subtypes of APLs to identify a molecular signature that distinctly separates APL subtypes. This paper reports the results of extensive molecular testing of a large cohort of well-characterized clinical samples from patients with diagnosed APLs.

2. Materials and Methods

2.1. Samples

In this multi-center retrospective study, a total of 166 fresh frozen (FF) tissue samples were obtained from different international biobank repositories, including TCBN (Caen, France), VCB (Melbourne, Victoria, Australia), Indivumed (Hamburg, Germany), IDIBAPS (Barcelona, Spain), LBIH Biobank (Liverpool, United Kingdome) and Dx Biosamples (San Diego, CA, United States).

96 and 46 paired tumor tissues and normal adjacent tissues from 48 APL and 23 stage I CRC were designedly selected for -omics detection. An independent set of APL samples were used for validation purposes.

Tissue samples with different histological background (sessile serrated lesions (SSLs), tubular adenomas (TAs) or villous/tubulovillous adenomas (VAs/TVAs)) were used to obtain DNA for whole genome bisulfite sequencing.

For our analysis, we define two samples as paired when extracted from the same patient though of differing tissue or sample type.

Demographic sample donor features are listed in Table 1:

2.2. Sample Preparation and Whole-Genome Bisulfite Sequencing (WGBS)

Genomic DNA (gDNA) from FF tissue was extracted using a DNeasy Blood & Tissue kit (Qiagen, Valencia, CA) according to a protocol from the manufacturer. Extracted gDNA was then fragmented into segments of about 400 bp with a Covaris S220 ultrasonicator. The extracted and shredded gDNA was bisulfite-converted with EZ DNA Methylation-Lightning kit (Zymo Research, Inc., Irvine, CA, USA). Sequencing libraries were prepared from the bisulfite converted DNA by using Accel-NGS Methyl-seq DNA library kit (Swift Biosciences [Ann Arbor, MI, USA]) and consequently sequenced to an average depth of 37.5x with NovaSeq6000 (Illumina, San Diego, CA, USA) equipment, using paired-end sequencing (2×150 bp). On average each sample had 335.58 M sequences.

2.3. Genomic Data and Methylation Processing

Quality checks and trimming were performed using FastQC (27) TrimGalore v0.4.5 (a wrapper tool around Cutadapt (28), which removed adapter sequences and poor-quality bases and reads).

The remaining high-quality reads were aligned to a bisulfite-converted human genome (Ensembl 91 assembly, hg38) using the Bismark Bisulfite Read Mapper (v0.20.0) (29), which makes use of Bowtie 1.2.1.1 alignment software (30). On average 89% of the sequences were aligned.

Methylation calls for every single C analyzed was performed by Bismark bismark_methylation_extractor script. For each CpG, Beta values (β) were calculated as:

β = CGmethylated / (CGmethylated + CGunmethylated),

where CGmethylated is the number of methylated cytosine and (CGmethylated + CGunmethylated) is the sum of methylated and unmethylated cytosine (total number of reads) at that position. The bisulfite conversion rate was estimated from Cytosine (C) in non-CpG context and was above 99%.

CpGs were first selected based on their coverage: CpGs not covered by at least 5 uniquely mapped reads in all samples were discarded for subsequent analyses such as those CpG with high beta difference variance, ending with 14848915 CpGs sites.

2.4. Analyses of Differentially Methylated Position/Region

To find differentially methylated regions, a pairwise analysis was performed using R in-house scripts, beta values were the input for all the analysis.

For each of the remaining CpGs, a paired Wilcox test was applied; p-values were corrected by Benjamin Hochberg (BH) algorithm and only those with a corrected p-value less than 0.05 and a methylation difference between the two groups of 0.2 were selected. Significant CpGs that fall in a SNP position and that have a population global minor allele frequency (MAF) more than 1% (dbSNP Build 155) were discarded.

Regions were formed summing CpGs that have a maximum distance of 50 bases. Only regions with 3 or more CpGs and more than 60% of significant CpGs respect to total number of CpGs of the region were selected for further analysis.

2.5. Pathway Enrichment Analysis

To analyze the biological meaning of our DMRs and add an extra layer of validation, we performed a pathway enrichment with KEGG (31) pathways. As a first step, an extensive mapping of the DMRs was performed to find putative cis-regulatory regions related to specific genes using the rGREAT function (32). Afterwards, with the associated genes information, we proceeded to perform a pathway enrichment using ClusterProfiler (33), with a pvalue of 0.05, and limiting the search to Homo sapiens. Resulting pathways were visualized using R package enrichplot (34).

2.6. Copy Number Alteration Detection

To evaluate large-scale copy number alterations (CNAs) the tool ichorCNA (35) was used. Read counts for each APL, CRC and respective normal samples were calculated with the HMMcopy Suite (36): each genome was divided into 10-kb non-overlapping bins and aligned reads were counted based on overlap within each bin; then, read counts were normalized to correct for GC-content and mappability biases. Centromeres were filtered based on chromosome gap coordinates obtained from UCSC (University of California at Santa Cruz) for hg38.

Normal samples were used to create a reference data set. This data was further used to correct systematic biases arising from library construction, sequencing platform, and DNA-specific artifacts.

Copy number analysis was performed on autosomal chromosomes only, using ichorCNA default parameters.

2.7. Somatic Mutation Detection

To overcome the bisulfite conversion problem, which can cause cytosine to appear as thymine during sequencing and complicate the detection of true mutations, a pipeline based on EpiDiverse/snp (https://github.com/EpiDiverse/snp) (37) was created. In summary, bam files were first preprocessed using samtools (sort, calmd, and index), and then a double-masking procedure, which manipulates specific nucleotides and base quality (BQ) scores on alignments, was applied.

For each sample pair, pileup files were created with samtools. The somatic function of VarScan2 (38)(39) was used to call somatic single nucleotide variants (SNVs) and small insertions/deletions (INDELs) from the comparison between AA or CRC and NATs. Only positions that were present in both files and met the minimum coverage in both files were compared. Only mutations identified by the tool as somatic were annotated by SnpEff (40). To summarize, analyze, and visualize the mutations, the VCF output from VarScan2 and SnpEff was converted into Mutation Annotation Format (MAF) files and analyzed using the R Bioconductor package Maftools (41).

2.8. CpG Island Methylator Phenotype (CIMP)

To predict the CIMP status of each sample, we utilized CpGs derived from promoter regions that exhibited a high standard deviation (SD >0.2) of methylation levels in tumor tissues and a low methylation level (mean β <0.05) in normal tissues. Consensus clustering analysis based on the K-means algorithm was performed using the R package ConsensusClusterPlus (42). Given the limited number of available samples for DNA methylation profiling in this study, we classified samples with hypermethylation patterns as CIMP-H, while all other samples were categorized as CIMP-L.

2.9. Microsatellite Instability (MSI) Status Prediction

Due to the converted data, standard short tandem repeat DNA sequence search methods are unable to detect microsatellite instability (MSI). Therefore, in this study, the MSI status was determined by examining potential frameshift mutations in genes involved in the mismatch repair (MMR) pathway and employing a prediction method called microarray DNA methylation tumor tissue profiling based MSI estimation (MSIMEP) (43).

2.10. Principal Component Analysis (PCA)

Principal component analysis (PCA) reduces the dimensionality of multivariate data, to two or three that can be visualized graphically with minimal loss of information. PCA was performed on epigenomics data to illustrate the differences between APL histologies using the prcomp function in R package stats. The R package factoextra was used to visualize the PCA output.

3. Results

We obtained 166 fresh frozen samples, fully characterized clinically and histologically using standardized criteria (described in the Methods section and Supplementary Table 1). To ensure robustness and generalizability, samples were divided into a discovery dataset of 48 APL samples paired with normal adjacent tissue (NAT), a validation dataset of 26 non-paired APL samples and a third dataset comprised of 23 paired CRC samples (Table 1). The discovery dataset was used to explore and identify the APL signature, performing differential methylation analysis and CNV/mutation analysis, while the validation dataset was used exclusively for assessing our signature performance in the agnostic identification of APLs using the signature method. Finally, the CRC dataset was used to compare and understand the differences between the precancerous state versus the invasive cancer state. Research pipeline is resumed in Figure 1.

3.1. DMR Analysis Shows Unique Methylation Patterns Between APL’s Subtypes and CRC

We conducted a thorough discovery process for DMR between APL subtypes, aiming to identify specific changes in methylation patterns across the genome to discern differences between subtypes. The analysis began with a pairwise profiling of different APL subtypes and their corresponding Non-Adjacent Normal Tissue (NAT) as controls, followed by comparison among the three histology types for both APL and NAT. The findings from these analyses are summarized in Table 2.

Starting from 122,348 DMRs, divided into 6,263 hypermethylated and 116,050 hypomethylated, we focused on the hypermethylated unique regions from each APL to extract relevant biological information.

As illustrated in Figure 2a, marked as red circles, we pinpointed the unique hypermethylated regions arising downstream for the identification of representative signatures for each stage in the pre-malignant evolution as can be seen in the Venn diagram (Figure 2b). We found 2,003 hypermethylated regions across APL subtypes, which were absent in the normal adjacent tissues. It should be noted that DMRs are found not to be correlated with age (average Spearman correlation of 0.06) or sample origin (average Spearman correlation of 0.09).

After analyzing the genomic distribution of the DMRs, we found a global presence of methylated regions in all chromosomes. In Figure 2b, a Circos representation of the human genome, external rings represent the chromosomes from 1 to 22, the next two rings are the density in a dot or density plot of DMRs in each part of the chromosomes sequence. Of particular note, there is significant enrichment of DMRs on Chromosome 7.

As expected, CRC samples were found to display unique methylation patterns different from those found in APL subtypes and NATs (Figure S1).

In totality, these results suggest that, during the evolution of malignant transformation, there exist specific methylation patterns that characteristically define each stage of APLs and, ultimately, CRC.

3.2. DMR Signature Enables Sample Stratification According to APL Subtype

Next, to reveal features intrinsic to each APL subtype and to reduce the number of variables that do not participate in the phenotype, we performed a Principal Component Analysis (PCA). DMRs were effectively able to separate APL samples into three groups of high-similarity components characterized by histological subtype in both discovery and validation sets (Figure 3a). These results suggest that DMRs can differentiate between APL subtypes.

To confirm that the signature differentiates between pre-cancerous lesions and of statistical significance, we performed a Kruskal-Wallis test for feature contributions of the first dimension between the three subtypes (Figure 3b). The analysis demonstrates significant differences between each histological APL (p=4.4e-05).

These results demonstrate that our signature can separate control samples from APLs and CRC and furthermore, consistently stratifies samples by APL subtype.

3.3. CNAs and Mutation Analysis Suggests Synergistic Mechanisms Complementary to Epigenetic Control in APLs

As CRC typically develops through a multistep process involving the accumulation of genetic and epigenetic alterations (44), we chose to support the DMR analysis with an analysis of Copy Number Alterations (CNA). The CNA signals were mostly present in TA or VA/TVA samples (Figure 3a). The most frequent signals could be seen in chromosomes 7, 12, 19 and 20, while the majority of CRC samples exhibited CNA in chromosomes 7, 8, and 20, suggesting that different processes exist between the emergence of APLs and progression to CRC. As in DMRs, chromosome 7 is also highly enriched with CNVs (Figure 4a, barplot).

Since APC and KRAS are the most frequently mutated driver genes in the adenoma-to-carcinoma sequence (6), they could be expected to demonstrate alterations more frequently, and MUC6 and MUC3A mutations were more frequent that KRAS mutations, principally in TA and VA/TVA (Figure 4b). However, SSL exhibits alterations in neither KRAS nor APC.

The integration of CNA and mutation analysis with differential methylation analysis provides a comprehensive view of the molecular landscape of APLs. Our analysis suggests that there is a clear distinction between driver events and passenger events in the different histological subtypes of APL: TA and VA/TVA present the strongest signal in copy number variation. Conversely, SSL presents the weakest signal in CNA, as well as no hits in KRAS or APC, suggesting a different driver event which may maintain the phenotype with epigenetic alterations. Early-stage colorectal cancer (CRC) exhibits a restructured mutation profile compared to advanced premalignant lesions (APLs), with notable alterations in MUC12 and FSIP2, observed in 43% and 39% of cases, respectively. Additionally, significant copy number aberrations (CNAs) were detected, particularly involving chromosomes 7 and 20. These findings suggest that both CNAs and mutations in specific genomic regions contribute to the differences observed among the different histological sub-types and CRC

3.4. Pathway Enrichment Analysis Reveal Redundancy in the Circuitry of Transcriptional Regulation Cellular Fate and Proliferation Signaling in APL Progression

To investigate the biology underlying the differential methylation signal identified, we performed a Gene Ontology (GOs) and pathway enrichment analysis based on DMRs obtained from the APLs that revealed distinctive subtype dysregulation in multiple biological pathways (Figure 5).

In the SSL vs. TA comparison (Figure 5a), we found several morphogenic and cell fate processes such as cell-cell adhesion, cell fate commitment and specification, and dorsal/ventral pattern formation. Moreover, when comparing SSL vs. VA/TVA (Figure 5b), we found slightly different processes but still related to cell fate, such as pattern specification process, regionalization cell fate commitment, gland development, plasma membrane adhesion molecules, anterior/posterior pattern specification and cell fate specification.

Finally, a TA vs. VA/TVA comparison (Figure 5c), showed enrichment in processes such as cell fate commitment, regulation of epithelial cell proliferation, miRNA transcription and regulation of stem cell proliferation, among other terms related to tissue reorganization.

These results suggested that genes affected by the DMRs in each subtype showed unique traits but also a redundancy at pathway level that support cellular proliferation, morphogenesis, and tissue reorganization. These molecular events suggest that cellular reorganization and the potential to re-organize tissue structures are critical to the development of APLs, which is consistent with the known clinical evolution of APLs and CRC.

Oncological pathway enrichment analysis using the somatic mutations found in APLs and early-stage CRC (Figure 5d) elucidate interconnected signaling pathways that help to alter the normal colonic structure and, in the case of CRC, contribute to an invasive phenotype. Crosstalk between dysregulated pathways, such as NOTCH, WNT, Hippo, and RTK-KRAS underscore the complexity of molecular events driving lesion progression. Interestingly, PI3K, MYC, TGF-Beta, cell cycle, and TP53 pathways are more affected in CRC than in the APLs. It is important to note that there are substantial differences among subtypes: SSL exhibits a lower a lower proportion of their pathway-associated genes altered; TA and VA/TVA show a higher enrichment in TGF-beta signaling but similar to normal tissue in other pathways; and VA/TVA is the only subtype that shows enrichment in the cell cycle pathway.

Collectively, these results provide a comprehensive molecular portrait of advanced precancerous lesions (APLs), highlighting the convergence of transcriptional misregulation and oncogenic signaling in their progression from normal tissue to APL and then to CRC.

3.5. DMR Phenotype in APLs Is Not CIMP and MSI Related

To address the question of whether the phenotype related to DMRs is related to CIMP status, we performed an analysis of specific methylation in promoter regions, classifying those regions with hypermethylation patterns as CIMP-H, while the samples without clear hypermethylation or hypomethylated were classified as CIMP-L (Supp. Figure 2a). Contrary to expectations, the DMR phenotype demonstrated independence from classical CpG island methylator phenotype. The absence of a strong correlation suggests the existence of alternative mechanisms driving differential methylation in advanced precancerous lesions.

To assess microsatellite instability status (MSI), we used the MSIMEP method (Supp. Figure 2b). Our results revealed no significant association with the identified DMR phenotype. The epigenetic alterations observed in APLs appear to be distinct from the genomic instability associated with MSI.

These results suggest that the phenotype of APLs is driven by molecular events unrelated to MSI or CIMP events.

4. Discussion

We present evidence of signatures related to aberrant methylated regions, chromosomal aberrations, and somatic mutations in advanced precancerous lesions.

The identification and analysis of methylated pattern differences between APL subtypes provides valuable insights into the epigenetic alterations that precede CRC. This study highlights significant DMRs across various subtypes of APLs, emphasizing their potential role in the progression from normal colonic mucosa to malignancy. Our analysis suggests that epigenetic changes are subtype-specific and may contribute to the heterogeneity observed in CRC.

The hypermethylation of genes involved in cell cycle regulation, DNA repair, and apoptosis observed in our study suggests that these pathways are crucially disrupted in the early stages of colorectal pre-neoplasia. For example, the methylation of the CDKN2A gene, a well-documented event in colorectal adenomas (45), was confirmed in our SSL vs. VA/TVA results, reinforcing its role in early tumorigenesis. Other genes as IRX4 are known to be silenced epigenetically and down expressed on pancreatic cancer (44) or AMER2 (also called FAM123A), known to be a negative regulator of WNT signaling and hypermethylated in precancerous stages in gastric carcinomas (46).

Our study identified distinct DMR signatures that enable the stratification of samples according to their APL subtype. Subsequent unsupervised clustering of these DMRs demonstrated a clear separation of samples into their respective subtypes, indicating the robustness of these methylation markers. Our findings were validated using an independent cohort. These DMR signatures not only facilitate the precise stratification of APL samples but also hold potential as biomarkers for early detection and personalized treatment strategies in clinical settings.

The analysis of CNVs revealed notable differences among the premalignant lesion subtypes. CNVs are known to contribute to genomic instability, a hallmark of cancer (47), and their presence in premalignant lesions underscores their role in early carcinogenesis (48).

SSLs present with lower gains compared to VAs/TVAs, focused on chromosomes 7, 12, 13 and 20, and the lowest losses compared with other APLs.

Here, it is important to highlight that chromosome 7 is enriched in CNVs for TA and VA/TVA as well as for DMRs (Figure 2b) as a possible transformation pathway has been reported before in CRC and adenomas (49). Interestingly, all the Mucin genes included in our mutation analysis are encoded on chromosome 7. This family of glycoproteins has been implicated in chronic inflammatory states and the promotion of oncogenesis in numerous malignancies, especially CRC (50).

These findings can inform the development of subtype-specific diagnostic and therapeutic strategies. Further research is needed to understand the functional impact of these CNVs and their potential as biomarkers for early detection.

Somatic mutations are pivotal for the initiation and progression of CRC and are known to be linked to epigenetic modifications (51). Our study identified several recurrent mutations across different premalignant lesion subtypes, providing insight into the molecular events that may drive their progression to invasive malignancy.

KRAS, TP53, and APC mutations were predominantly found in VA/TVA, while SSL showed mutations on APC and TP53. This aligns with known molecular pathways involved in the adenoma-to-carcinoma sequence in VAs/TVAs but represents a finding previously unreported in the serrated pathway.

SSLs exhibited a higher mutational burden compared to VA/TVA samples, suggesting a more heterogeneous clonal evolution. This could imply that SSLs have a higher potential for malignant transformation due to their diverse mutational landscape. Of note, our results showed that the SSL phenotype is not associated with CIMP or MSI. The independence of this epigenetic signature suggests alternative pathways driving methylation changes and emphasizes the importance of exploring distinct molecular mechanisms in the context of precancerous lesion progression.

The identification of these somatic mutations provides crucial insights into the molecular heterogeneity of premalignant lesions.

The analysis of biological pathways has uncovered notable redundancy in the transcriptional regulation circuitry that governs cellular fate and proliferation signaling in the APLs. This redundancy, manifested through multiple overlapping pathways, ensures robust control of critical cellular processes despite genetic and epigenetic alterations. Such redundancy may facilitate malignant transformation by providing alternative routes for sustaining proliferative and survival signals. These findings underscore the complexity of transcriptional networks in CRC progression and highlight potential challenges and opportunities for therapeutic interventions targeting these redundant pathways to effectively disrupt the progression from premalignant lesions to invasive carcinoma. Further, understanding these differences could aid in developing targeted prevention strategies and improve the precision of early detection methods.

5. Conclusions

Our integrated analysis of differentially methylated regions (DMRs), copy number variations (CNVs), and somatic mutations across various subtypes of advanced premalignant lesions of the colon provides a comprehensive understanding of the molecular alterations involved in early colorectal carcinogenesis. These results emphasize the potential of these molecular markers as valuable biomarkers for early detection and as targets for preventive strategies.

We are aware that our study has several limitations that should be addressed in future research. First, the validation of these biomarkers requires larger, independent cohorts, ideally analyzed in a double-blinded manner to reduce potential bias. Second, functional analyses are needed to investigate the biological interplay between methylation changes and chromosomal aberrations, which may help clarify their roles in the progression of colorectal cancer.