CTCF Binding Site Mutations: Linking Topologically Associated Domains Dysregulation to Cutaneous Squamous Cell Carcinoma Progression

CTCF Binding Site Mutations: Linking Topologically Associated Domains Dysregulation to Cutaneous Squamous Cell Carcinoma Progression | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article CTCF Binding Site Mutations: Linking Topologically Associated Domains Dysregulation to Cutaneous Squamous Cell Carcinoma Progression Amarinder Thind, Bruce Ashford, Nicholas Shannon, Dario Strbenac, and 6 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-6844715/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Background Cutaneous squamous cell carcinoma (cSCC) is the most common lethal malignancy with metastatic potential. The high mutational burden in cSCC has made it difficult to understand the significance of variants in the noncoding and regulatory genome. This study presents the first investigation of mutations at CCCTC-binding factor binding sites (CTCFbs) of topologically associating domains (TADs) across defined stages of disease progression —primary tumours that have not metastasized, primary tumours that have metastasized, and lymph node metastasis. Results By integrating matched whole-genome sequencing and RNA sequencing data from the same tumors, we reveal that CTCFbs are mutation hotspots (~1,100 mutations/Mb) in cSCC, with mutation densities far exceeding genome-wide averages (170–250 mutations/Mb). This study is also the first to prove genome-wide association of CTCFbs mutations with gene expression changes in cSCC. We report the novel finding that CTCFbs that overlap with other regulatory elements such as promoters and untranslated regions exhibit even higher mutational densities than overall CTCFbs. A pattern of mutually exclusive TAD loop CTCFbs mutations was observed between non-metastasizing and metastasizing primary tumors. TAD loops with mutated CTCFbs, are associated with significant transcriptional changes in the genes within those TADs, implicating genes including HHIPL2 , LINC02870 , and GDA in cSCC progression. Conclusions Our findings highlight the functional relevance of noncoding mutations at CTCFbs in cSCC and suggest their potential influence as drivers of tumor progression and metastasis. This integrative genomic analysis with detailed examination of TADs provides a foundation for future studies into 3D genome dysregulation in skin cancers. CTCF binding sites Topologically Associated Domains Cutaneous squamous cell carcinoma Whole Genome Sequencing 3D genome architecture Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Background Cutaneous squamous cell carcinoma (cSCC) is the most common malignancy with a propensity for aggressive metastasis. Resectable advanced disease requires surgery and adjuvant radiation to achieve control. This malignancy is characterized by an extremely high tumor mutation burden (TMB) across the genome ( 1 , 2 ). Understanding the genomic events driving cSCC progression, especially the transition from localized to metastatic disease, is crucial for improving patient outcomes. Previous research has underscored the potential utility of using TMB to predict immunotherapy response and as a vital tool for understanding tumor evolution ( 3 ). However, the specific genomic regions and regulatory elements contributing to TMB variability in highly mutated tumours like cSCC remain poorly defined. Central to this investigation is CTCF, a protein that binds to the core genome sequence CCCTC, known as CTCF binding sites (CTCFbs) ( 4 ). CTCF plays diverse roles, including transcriptional regulation, insulation (blocking enhancer-promoter interactions), chromatin architecture regulation, and RNA splicing. For example, CTCF organizes chromatin into topologically associated domains (TADs), creating self-interacting genomic regions that facilitate interactions between regulatory elements and genes ( 5 ). TADs are fundamental to gene regulation, as they compartmentalize the genome, restricting interactions to within the TAD and preventing the spread of heterochromatin's silencing effects ( 6 , 7 ). We have previously shown that CTCFbs are prone to UV-induced mutations in 15 cSCC lymph node metastases. Interestingly, cSCC exhibits CTCFbs mutation at a rate seven times that seen in melanoma, also a UV-mediated cancer ( 2 ). Indeed, our previous work identified 1,979 genes within 422 affected TADs, including both tumor suppressor genes and oncogenes. These findings suggest that mutations in CTCFbs could significantly disrupt regulatory processes, leading to gene expression abnormalities associated with tumor progression. While mutations in CTCFbs have been observed in various cancers ( 8 , 9 ), the impact of CTCFbs mutations on gene expression during cSCC progression has not been investigated. In this study, we present a comprehensive analysis of mutation density and substitution patterns at CTCFbs within TADs in a clinically stage-defined cohort of cSCC samples. This includes high-risk non-metastasizing primary tumors (PNM), metastatic primary tumors (PM) and lymph node metastases (LNM). We compared this to the distribution of mutations among various genomic regions across the entire genome and examined the impact of CTCFbs variants on gene expression (by integrating whole genome and RNA sequencing) with the aim of identifying potential drivers of metastasis (Fig. 1 ). Results TMB and Substitution Types of the Overall Cohort The mean genomic TMB across all samples in this study was 199.46 mutations/Megabase (Mb), with a median of 145.56 mutations/Mb. The TMB values ranged from 7.59 to 999.72 SNVs/Mb, reflecting a broad and heterogeneous distribution (Fig. 2A). The standard deviation of 200.85 mutations/Mb further underscores the high intersample variability in mutation density. Among the samples, CSCC_0126-M1 presented the highest TMB (999.72 mutations/Mb), whereas CSCC_0157-P1 presented one of the lowest values (7.58 mutations/Mb). This variation suggests potential biological differences in mutation accumulation across tumor samples. Mutational density across Genomic Regions and their Background First, a comparison of different noncoding/coding genomic regions with each other was conducted to investigate whether certain known genomic regions present higher or lower mutational densities relative to each other (Fig. 2B). The average mutation density at the CTCFbs, which is based on a 20 bp motif, was greatest with a mutational rate of approximately 1100 mutations per Mb. In contrast, the average mutation rates for the other five regions ranged between 170 and 250 per Mb, with long noncoding RNAs (lncRNAs) regions exhibiting the second highest average mutation rate at ~ 260 per Mb (Fig. 2B). Second, we investigated the surrounding regions of each genomic region at various scaling factors to determine whether the TMB is elevated significantly for a specific region than in its surrounding region (Supplementary Fig. 1). The genomic regions tested included the 3' UTR, 5' UTR, promoter regions, lncRNA regions, CDS, and others. When these regions were compared with their surroundings, a sharp decrease in mutation density was observed at both ends of the CTCFbs (Fig. 2C; Supplementary Fig. 1). For the other regions, the decrease was more gradual, and the greater the distance from the region was, the greater the reduction in mutation rate, except for CTCFbs, which showed a distinct pattern (Supplementary Fig. 1). This observation may indicate unique regulatory or structural features at CTCFbs that influence mutation accumulation. A notably high mean mutation density near the CTCFbs centre (position 0) indicated differences in mutation accumulation, suggesting a specific mutation profile associated with these binding sites (Fig. 2C). Additionally, the substitution rates per Mb at the CTCFbs (± 20 bp) highlighted the distribution of different mutation types (Fig. 3). C > T substitutions are more prevalent at the core of CTCFbs (positions 3 and 4 on Anchor 5 and positions − 4 and − 5 for Anchor 3). CTCFbs and Regulatory Element Interactions CTCFbs that overlap with other genomic elements present a greater mutational density than do CTCFbs regions alone (Fig. 4). Among these, CTCFbs overlapping with lncRNAs display the highest TMB, surpassing those overlapping with coding sequences (CDSs), untranslated regions (UTRs), or promoter regions (Supplementary Table 1). This enrichment suggests that CTCF-bound coding and regulatory elements are particularly susceptible to tumor-associated mutations, potentially affecting gene expression, protein function, and chromatin architecture. Additionally, while the density of CTCF motif mutations is typically averaged over a 20 bp window, certain positions show markedly higher mutation rates. Specifically, as shown in Fig. 3, positions 3 and 4 on the CTCFbs motif of Anchor 5 and positions − 4 and − 5 on the motif of Anchor 3, show mutation densities ranging from approximately 5,000 to 7,000 mutations per Mb, indicating potential hotspots for functional disruption. CTCF loops across the entire cohort The lengths of the CTCF loops varied significantly, ranging from 48.9 kb to 1904 kb. A total of 667 (out of 903 TAD loops with CTCFbs in both anchor regions) were mutated in 72 samples. The number of mutated samples per loop ranges from 1 to 29, with frequently mutated loops could represent key regulatory elements or regions with critical roles in the pathogenesis of the disease. Substitution types This analysis revealed that C > T substitutions were dominant across most samples, with a few samples showing a notable proportion of T > G and T > C substitutions (Supplementary Fig. 2). Interestingly, two samples presented very few C > T substitutions, suggesting potential sample-specific mutational patterns or underlying biological factors influencing mutation processes. The predominance of C > T substitutions could reflect mechanisms such as deamination of methylated cytosines, a common mutational signature in skin cancers. Overall mutation density and substitution types across cSCC clinical groups Mutation density per sample and per group was compared among PNM, PM, and LNM groups using the Kruskal-Wallis test, which yielded a non-significant result, indicating no statistically significant differences. To further explore whether mutation patterns differ between clinical groups, we analysed the mean mutation density of substitution types across groups (Supplementary Fig. 3). The relative contribution of each substitution type varied slightly among the groups. In general PNM had a lower CTCFbs mutation density compared to PM and LNM. Differences in CTCFbs Mutations : Analysis of mutation patterns within the CTCFbs regions of TADs revealed no substantial differences across clinical groups. Supplementary Fig. 4 displays the distribution of mutation density per mb (adjusted for sample size and number of regions) across CTCFbs motifs. CTCFbs in Anchor5 (top panels) showed a enrichment of C > T transitions, particularly at positions 3 and 4, corresponding to the "G and G/T" sites of the motif (middle panel). In contrast, the CTCFbs in Anchor3 (bottom panels) exhibited a notable peak of C > T transitions at positions − 4 and − 3 across all the groups. Association between CTCF binding site mutations and the expression of harbouring and nearby genes As described in the methods, we grouped samples with and without mutations for each loop using 54 shared WGS/RNA-Seq samples. Only 284 loops had mutations in at least three samples, which was set as the minimum group size required for differential gene expression (DGE) analysis. Notably, 17 of these 284 loops contained significantly differentially expressed genes (DEGs) when samples groups with and without mutations (regardless of clinical subgroup) were compared (Table 1). Out of these 17 loops, 3 overlaps with 3’UTR regions. Table 1 Shows the 17 loops that harbour significant differentially expressed genes (shown) between samples with and without mutations regardless of cohort subgrouping. Differential expression was considered significant using a threshold of log₂FC ≥ 1 or ≤ − 1 and adjusted p-value (adj. P) < 0.05. Only DEGs meeting these criteria are shown. Loop No. of Significant DEGs DEGs loop_1057 6 ENSG00000227726, ENSG00000279459, ENSG00000226627, ENSG00000280089, ENSG00000171671, ENSG00000286708 loop_21 3 RPS27P20, ENSG00000281386, ENSG00000281097 loop_28 2 TNFRSF8, TNFRSF1B (ENSG00000028137) loop_310 2 MT1L, ENSG00000205364 loop_404 2 ENSG00000287920, ENSG00000290062 loop_364 2 BOC, ENSG00000288079 loop_213 2 ENSG00000285090, ENSG00000285964 loop_682 2 LINC02783, ENSG00000159339 loop_204 1 PIP loop_336 1 ENSG00000233633 loop_75 1 RGS6 loop_894 1 CFAP221 loop_156 1 KALRN loop_26 1 CASQ2 loop_1033 1 FAM78A loop_207 1 RAMP3 loop_697 1 CD163 Furthermore, to test whether the association was not by random chance, we performed a permutations test with 1000 iterations, as described in the Methods, the observed number of loops associated with significant DEGs was 17, which exceeded the maximum number observed in the null distribution (16; Supplementary Fig. 5). This yielded an empirical p -value < 0.001, supporting that the observed loop-DEG associations are unlikely to have occurred by random chance. Mutually Exclusive CTCFbs mutations between metastatic vs non-metastatic subgroups We observed that certain CTCFbs were uniquely mutated in primary tumors that metastasized (PM, n = 15) compared with those that had not metastasized (PNM, n = 16) as shown in Fig. 5A. Six loops are only mutated in the PNM group, whereas 4 are uniquely mutated to the PM group, suggesting that mutations at specific CTCFbs may be linked to the metastatic potential of cSCC. However, these mutations are not present in every sample within each clinical subgroup, highlighting the heterogeneity of mutation patterns across individual cases; for example, loop_45 is mutated in 5/16 (31%) PNM samples. Importantly, the genes associated with the mutated loops (Fig. 5B) are previously known to be linked to cancer progression. We further, investigated whether the loops with different CTCFbs mutations between subgroups groups (refer to the methods section) are associated with the differential expression of genes within the TAD loop being tested. For these specific groups, first we considered samples for which we had both WGS and RNA-Seq in our cohort PM (n = 14) vs PNM (n = 11). Using these two groups, we performed DGE analysis and only genes with a |log2FoldChange| > 1 and a p-value ≤ 0.01) were considered as significantly DEG. GDA and HHIPL2 are observed as significant in this comparison. Using Spearman’s rank correlation, we identified genes whose expression levels were significantly associated with the mutation status of the loops in which they reside. Several genes exhibited strong positive or negative correlations (Fig. 6; Supplementary Table 2), indicating differential expression linked to loop mutations. To complement this, AUC values were calculated, revealing high classification performance for some genes in distinguishing between mutated and non-mutated loops based on expression patterns (for details refer to Method section). This is done for each loop (Fig. 5A) by focusing on a particular stage. ENSG00000259093 , located within loop_70, showed a strong negative correlation (Spearman ρ = -0.51) and a high AUC of 0.92, indicating consistent downregulation in the presence of loop mutations. In contrast, LYRM9 in loop_95 exhibited a positive correlation (ρ = 0.58, AUC = 0.91), suggesting that upregulation was associated with mutated loops. Notably, a subset of gene–loop pairs demonstrated exceptionally high classification performance. For example, AKAP6 (loop_1083) achieved a perfect classification (AUC = 1.0), whereas LINC02068 (loop_284) reached an AUC of 0.94 with a strong positive correlation (ρ = 0.60). These findings suggest that gene expression within certain loops may serve as an effective proxy or biomarker for structural disruption caused by somatic mutations. Interestingly, some gene–loop pairs presented low correlations but moderate AUC values, which may indicate nonlinear or threshold-based expression responses, potentially mediated by epigenetic buffering or compensatory mechanisms. We also performed other group comparisons for dissimilarity in CTCFbs mutations; however, no significant biological conclusion was drawn from this analysis. Table 3 shows various comparisons and parameters used for these analyses and a list of significant DEGs found in each of those comparisons. Table 3 Other Group Comparisons for Dissimilarity in CTCFbs Mutations. This table summarizes the results of additional group comparisons performed to assess the dissimilarity in CTCFbs mutations across various tumor stages. While no significant biological conclusions were drawn from these analyses, we report the comparison parameters used and the differentially expressed genes (DEGs) identified in each group. The comparisons included different tumor stages (PNM vs PM, primary vs Met, and PM vs Met) along with specific anchor region considerations, loop filter criteria, and the list of significant DEGs for each comparison. Comparison Sample Groups (WGS/RNA-seq) Anchor Regions Loop Filtering Criteria Significant DEGs PNM vs PM 15/11 vs 16/14 5′ and 3′ anchors p-value < 0.3 & (sumPNM ≤ 1 or sumPM ≤ 1) GDA, HHIPL2 Primaries (PNM + PM) vs LNM 31/25 vs 41/29 5′ and 3′ anchors p-value < 0.3 & (sumPrimary ≤ 1 or sumLNM ≤ 1) LINC02884, LEFTY2, MARK2P13 PM vs LNM 16/14 vs 41/29 5′ and 3′ anchors p-value < 0.4 & (sumPM ≤ 1 or sumLNM ≤ 1) ENSG00000261184, LINC02884, ENSG00000259093 Discussion This study investigated the impact of mutations at the CTCFbs of TADs on gene expression that lies within TADs and their role in cSCC cancer progression. In line with the very high TMB reported in our previous studies in cSCC ( 1 , 2 , 10 ), the mean genome TMB across all samples was high compared with that of other cancers. We observed an exceptionally high mutation density at the core of the CTCF motif, with notable enrichment of C > T substitutions particularly. Mutations were significantly concentrated at the core nucleotides of the CTCF binding motif, suggesting that these regions are particularly vulnerable to DNA damage. This pattern aligns with known mutational mechanisms, including cytosine deamination at methylated CpG dinucleotides, which frequently results in C > T transitions—a common signature in cancer genomes ( 11 ). The high mutation rate at CTCFbs—especially those overlapping with regulatory elements like long non-coding RNAs, UTRs or coding sequences (Fig. 4)—suggests potential disruption of local chromatin structure or enhancer-promoter insulation. However, the relatively short lengths of the CTCFbs overlapping regulatory regions and the limited sample size in this study constrain our ability to draw strong statistical conclusions or biological interpretations. Despite these limitations, the observed trends reveal intriguing patterns that merit further investigation. We also examined whether the 17 CTCFbs-associated loops linked to significantly DEGs overlapped with other genomic regions. Interestingly, 3 of the 17 loops binding site motif overlap with 3’UTR regions (Supplementary-Table-1). These loops are loop_28, loop_894, loop_364 and contains cancer progression related genes such as BOC, TNFRSF8 (CD30), and TNFRSF1B ( 12 – 15 ). Theses overlap needs future more comprehensive analyses to validate and extend these preliminary observations. CTCF binding site mutations and gene dysregulation : Among the 17 DEG-associated loops, several harbor genes implicated in hallmark processes of cancer such as immune evasion, proliferation, and metastasis. MAP4K1 has been reported to act as either an oncogene or a tumor suppressor gene (TSG), depending on the signaling environment ( 16 ). Similarly, RPS6KA3 (also known as RSK2 ) demonstrates dual roles: while often overexpressed in cancers, it can suppress proliferation under certain conditions ( 17 , 18 ). ATP1B1 also acts as a tumor suppressor, with evidence supporting its role in inhibiting cancer cell proliferation and migration ( 19 , 20 ). In contrast, NOS2 displays mixed behavior, functioning either as a tumor suppressor or a tumor promoter depending on the biological context ( 21 , 22 ). TNFSF10 functions predominantly as a tumor suppressor, known for its strong pro-apoptotic activity in cancer cells ( 23 ). TNFRSF8 (CD30) is a member of the tumor necrosis factor receptor superfamily and serves as a tumor marker. It is found on the surface of specific cells, including certain immune cells and cancer cells( 13 ). CD30 is also the target of the FDA-approved therapeutic brentuximab vedotin (Adcetris)( 24 ). LINC02870 promotes triple negative breast cancer ( 25 ) and hepatocellular carcinoma progression ( 26 ). PIP gene expression decreases gradually with increasing stage and grade of breast cancer ( 27 ). These associations underscore the potential biological importance of mutations within the CTCFbs of TADs in modulating gene expression in cSCC. Differential Loop Mutations and Metastatic Potential When PM and PNM tumors were compared, specific CTCF loops (e.g., loop_45) were found to be recurrently mutated in one group but not in the other. Interestingly, while mutations in a given loop were observed in up to ~ 31% of samples within a group, suggesting a recurrent, yet heterogeneous pattern of loop disruption, the DEG analysis still revealed that HHIPL2 and GDA were significantly altered in the PM vs PNM comparisons (regardless of loop mutation or not). HHIPL2 has been previously associated with lung acinar adenocarcinoma ( 28 , 29 ), further supporting its role in cancer metastasis. In addition, GDA is known to play a direct role in skin carcinogenesis by interacting with several cytokines and growth factors ( 30 ). These findings suggest that CTCFbs loop mutations may serve as potential markers of metastatic potential, but that CTCFbs mutation alone may not be only driver of gene expression changes. However, mutations were not uniformly present across all samples in a group, which may reflect tumor heterogeneity or stochastic mutation patterns. In the comparison of primary (PM + PNM) vs LNM, LEFTY2 is found to be downregulated in LNM. LEFTY2 is a member of the TGF-β superfamily that functions as a tumor suppressor by inhibiting epithelial–mesenchymal transition (EMT) ( 31 ). Its downregulation or epigenetic silencing has been linked to increased stemness, invasion, and metastatic potential in several cancers, including endometrial cancer ( 32 ). Functional impact vs mutational noise : Although many loops harboured mutations, only a small subset were associated with significant changes in gene expression. This may reflect a combination of factors including the non-functional nature of some mutations, compensatory mechanisms, or limitations of bulk RNA-seq in capturing cell-type-specific regulatory changes ( 33 , 34 ). Moreover, given that only ~ 31% of samples in a group may harbor mutations in any given loop, the effect size on transcriptomic output may be diluted in population-level analyses. Statistical and Experimental Limitations : While permutation testing confirmed that the observed DEG-loop associations are unlikely to be due to random chance (empirical p = 0.00), several limitations remain. The sample size, particularly in subgroup comparisons (e.g., PM vs. PNM), reduces the power for detecting subtle effects. Bulk RNA-seq data also limits the ability to infer regulatory consequences in specific cell populations, especially for lineage-specific genes. Furthermore, while our data suggest that mutations at CTCFbs can dysregulate gene expression, direct mechanistic validation (e.g., through ChIP-Seq, CRISPR editing, or 3D chromatin conformation assays) is necessary to confirm causality. Biological and Clinical Implications This study supports the hypothesis that noncoding mutations, particularly those in architectural elements such as CTCFbs, can contribute to tumor progression through dysregulation of TAD loops and associated gene expression. Given the growing interest in targeting chromatin architecture and epigenetic dysregulation in cancer, these findings may inform new strategies for biomarker discovery or therapeutic targeting in cSCC. Conclusion In this study, we reveal that CTCFbs are hotspots of mutation accumulation across all clinical stages of cSCC, exhibiting a markedly higher TMB compared to other coding and non-coding elements. Notably, mutations at CTCFbs were not uniformly distributed and showed distinct substitution patterns, particularly enriched for C > T transitions, consistent with models of UV-associated carcinogenesis. Through integrative analysis of matched RNA-Seq and WGS data, we identified a subset of CTCFbs loops where mutations were significantly associated with altered expression of genes located within or near the corresponding TAD. Among these, genes including HHIPL2 , GDA , and LINC02884 —known to be involved in oncogenic pathways—were differentially expressed, linking CTCFbs mutations to functional consequences in cSCC pathogenesis. Permutation testing confirmed these associations were unlikely due to random chance. Furthermore, we observed that certain CTCF-bound loops exhibited differential mutation patterns between PNM and PM, suggesting a potential role for CTCFbs mutations in promoting metastatic behaviour. Overall, our findings underscore the central role of CTCFbs and their genomic context in modulating mutational landscapes and gene expression in cSCC. By uncovering how disruptions at these critical regulatory hubs may contribute to cancer development and progression, this study offers new insights to inform biomarker discovery and therapeutic targeting strategies for skin cancers. Methods Sample Collection: This study was conducted with approval from the Institutional Human Research Ethics Committee (HREC/15/RPAH/266). Patients with resectable metastatic cSCC and high-risk non metastasising primary tumors were prospectively identified by the treating surgeons prior to surgery. Clinicopathological data, including age, sex, extent of nodal metastases, histology, and immunosuppressive status, were collected (Table 4 , Supplementary Table 3). Fresh tumor tissue from nodal metastases (n = 41) was harvested during surgery and immediately snap-frozen. A total of 72/75 cSCC, WGS samples were included in the analysis (refer to QC section), comprising primary tumors (both metastatic (PM (n = 16)) and non-metastatic (PNM (n = 15)), as well as lymph node metastases (LNM (n = 41)). The metastatic cohort included matched PM and LNM samples for 13 patients. Further, RNA sequencing (RNA-Seq) data were collected for 54/72 samples, representing both primary and metastatic stages of the tumor. The PNM group had to meet the following criteria: absence of metastases at > 24 months follow-up after resection of the primary or negative sentinel lymph node biopsy at time of resection or histologically negative neck dissection. Whole genome and total RNA sequencing and QC: Tumor tissue sections were processed for DNA and RNA extraction (using Qiagen AllPrep DNA and RNA kits, Qiagen, Hilden, Germany) and for estimating tumor cellularity. Only samples with a tumor content of greater than 30% (range: 35–95%) proceeded to DNA quality control (QC). QC procedures included spectrophotometry (Nanodrop 2000, Thermo Fisher Scientific Inc.), and gel electrophoresis. WGS was performed by AGRF (Melbourne, Australia), Genome.One (Darlinghurst, Australia) on Illumina HiSeq X to a depth of ×30–45 for whole blood (germline DNA) and ×65–90 for tumor samples. on the Illumina NovaSeq 6000 platform (Illumina). The average sequencing coverage was 94.56× (range: 64–143) for tumor samples and 41.08× (range: 30–56) for blood samples. Of the 75 samples sequenced, 72 passed QC. The remaining 3 samples exhibited extreme GC bias. Initially, total RNA was QCed before sequencing was performed by spectrophotometer. Following sequencing, the raw RNA-Seq data underwent quality control using the bioinformatics tool FastQC (version 0.11.9; Andrews, 2010). Low-quality reads were removed using Trim Galore (version 0.4.5) ( 35 ). After quality filtering, a total of 54 RNA samples were retained for downstream analysis. WGS and RNA-Seq pre-processing Somatic Variant Analysis Using DRAGEN (v4.3.6) on ICA v2: Somatic variant analysis was performed using DRAGEN pipeline version 4.3.6 on the Illumina Connected Analytics (ICA) v2 environment using in house shell scripting. Tumor-normal paired analyses were carried out in three stages: alignment, variant calling, and integrative somatic analysis. Further information on the somatic calling is available at https://help.dragen.illumina.com/product-guides/dragen-v4.3/dragen-dna-pipeline/small-variant-calling/somatic-mode . (a) Alignment of Normal Samples : FASTQ files from normal samples were aligned using the DRAGEN Germline pipeline (v4.3.6) with the GRCh38 ( hg38-alt_masked.graph.cnv.hla.rna_v4.tar.gz ) reference genome. Small variant calling was enabled to generate a germline VCF, which was used as input for downstream somatic calling. (b) Alignment of Tumor Samples : Tumor sample FASTQ files were aligned using the DRAGEN Somatic pipeline (v4.3.6), using the same reference as for the normal samples. (c) Tumor-Normal Paired Analysis : Paired analysis was performed to identify somatic SNVs. Inputs included: Aligned tumor and normal BAM and BAI files, Germline SNV VCF from the normal alignment step, A systematic noise BED file (downloaded from: https://webdata.illumina.com/downloads/software/dragen/resource-files/sv-systematic-noise-baseline-collection-3.0.0.tar ). Post Somatic Calling Filtering is done as default and hard filtered . vcf files are used for further analysis. RNA-Seq data processing: RNA sequence reads were mapped with STAR version 2.7.10a ( 36 ) onto the GRCh38. On average, for each sample 70% of the reads aligned to the reference genome (average > 87 million reads/sample). Transcript abundance was measured in terms of read counts using the same annotation file used for the transcriptome assembly, leveraging the featureCounts ( 37 ), R Bioconductor package, default parameters. The count matrix was used as input for gene differential expression analysis. TAD loops with CTCF motifs and mutations at CTCFbs TAD loops containing CTCF motifs were identified in NHEK tissue using the same methodology as described in Mueller et al., 2019 ( 2 ). This involved utilizing chromosome conformation capture (Hi-C) TAD maps from NHEK ( 38 ) along with chromatin immunoprecipitation sequencing data from ENCODE (The ENCODE Project Consortium, 2012) ( 39 ). A 20-bp motif and a 20-bp position-weighted matrix ( 4 ) were applied to select binding sites with high CTCF-binding probability. TADs were filtered as shown in Fig. 1 . TADs were excluded if either anchor region contained more than one CTCF binding motif, or if the CTCF motifs were not in a convergent orientation, as this configuration is most strongly associated with CTCF binding ( 38 ). The genomic coordinates defining each TAD were determined as the 3′-end of the upstream motif and the 5′-end of the downstream motif. 903 TAAD loops with CTCFbs in both anchor regions were identified. Identification of Mutated CTCFbs: Binding sites genomic co-ordinates of all 903TAD-loops were intersected with our cohort’s WGS somatic mutations data using bedtools intersect (version V2.31.1;( 40 )) to obtain a mutational data matrix (samples vs loops) that contains information’s of mutated loops in each sample. Further, Genes lying with-TADs and 1000bp outside either side of the CTCFb-motif in the TADs were extracted using Biomart R package ( 41 ). Mutations in CTCF binding motif vs background: To assess the mutational enrichment in CTCF binding motifs, we defined the motif region as ± 10 bp around the center of each CTCF binding site. The background was defined as the flanking regions of ± 1kb relative to the binding site center. CTCFbs and their flanking 1kb regions were defined and merged into a GenomicRanges object. Mutation data was loaded and represented as a GenomicRanges object . Overlaps between mutations and the 2kb CTCF regions (± 1kb) were identified, and mutation positions were normalized to a ± 1kb scale centered on the CTCF binding site. For each position within the 2kb region, the mutation rate was calculated as follow: Let R(x) represent the mutation density at position x within the 2kb region: \(\:R\left(x\right)=\frac{\text{N}\left(x\right)}{\text{L}.\text{S}}\:\:\) × 10 6 Where: N(x) = Number of mutations observed at position x across all samples and loops × 2 (both anchors) L = Total number of loops × 2 S = Total number of samples 10 6 = Scaling factor for mutations per megabase (Mb) Mutation densities for each region type (CTCFbs and surrounding) were calculated for each sample by aggregating the total number of mutations across all motif regions (1,806), dividing by the total genomic length of these regions collectively, and scaling to mutations per megabase. Exploring Various Genomic Elements: We investigated several genomic elements, including 3' UTR, 5' UTR, promoter regions, long non-coding regions, and CDS regions, to assess (a) whether these elements tend to harbor more mutations compared to their surrounding regions, and (b) how mutational density varies due to differences in the region lengths within each element. To conduct these analyses, we developed an in-house script, which is available upon request via GitHub ( https://github.com/amarinderthind/CTCF_Cancer_study ). For the first hypothesis, we defined surrounding regions for each genomic element using different scaling factors (10, 50, 100, 200, 400, 600, and 800), which were applied based on the region length. The scaling factor was applied to calculate the upstream and downstream boundaries of each region, where the surrounding regions' length was determined by multiplying the element length by the scaling factor (Supplementary Fig. 1). For each scaling factor, the formula for defining the surrounding region was: Surrounding Region Length = Region Length × Scaling Factor For each region, we computed the number of mutations that overlap with the surrounding regions, and the mutation density was calculated by normalizing the counts to mutations per megabase (mutations/Mb). To account for varying sample sizes and genomic element lengths, we normalized the counts using the total number of CTCFbs motifs regions (903*2 = 1806) and total sample size (72). The scaling was applied as: $$\:\text{M}\text{u}\text{t}\text{a}\text{t}\text{i}\text{o}\text{n}\:\text{D}\text{e}\text{n}\text{s}\text{i}\text{t}\text{y}\:\left(\text{p}\text{e}\text{r}\:\text{M}\text{b}\right)=\frac{\text{}\text{N}\text{u}\text{m}\text{b}\text{e}\text{r}\:\text{o}\text{f}\:\text{M}\text{u}\text{t}\text{a}\text{t}\text{i}\text{o}\text{n}\text{s}}{\left(\text{R}\text{e}\text{g}\text{i}\text{o}\text{n}\:\text{L}\text{e}\text{n}\text{g}\text{t}\text{h}\right)\text{x}\:\text{n}\text{u}\text{m}\text{b}\text{e}\text{r}\:\text{o}\text{f}\:\text{s}\text{a}\text{m}\text{p}\text{l}\text{e}\text{s}}\:\text{x}\:\text{1000,000}\:$$ Additionally, mutations were annotated to identify whether they occurred within the target region (e.g., 3' UTR, promoter, etc.). We then compared mutation densities between the regions and their surrounding elements at each scaling factor. For visualization, mean mutation densities across all samples were calculated and compared for each scaling factor (10x, 50x, 100x, etc.), highlighting differences in mutation patterns across genomic elements and their surrounding regions. Mutational Substitution Types in CTCFbs motif and cSCC progression: To characterize mutational patterns across samples, we analyzed single nucleotide substitutions across whole genomes. Substitution types were classified into six categories (C > A, C > G, C > T, T > A, T > C, T > G). To compare substitution profiles across sample groups, mean mutation proportions were computed and assessed using Wilcoxon rank-sum tests. To assess differences in mutational profiles during cSCC progression, we compared motif-specific mutations across different disease stages. Association between CTCFbs mutations and the expression of harbouring and nearby genes As shown in Fig. 1 G, Differential gene expression (DGE) analyses were performed for each loop using data from the mutation's matrix and the RNA-Seq count matrix. For each loop, two comparison groups were created based on the presence or absence of observed mutations at CTCFbs to assess the association of these mutations with the expression of genes within the loop and nearby genes (within 1 kb). Loops with CTCFbs mutated in fewer than 3 samples were excluded from this DGE analysis and 54 shared samples of RNA-Seq/WGS were considered for this analysis. In total, 284 DGE analyses were conducted, one for each loop. A table containing the loops and their associated significantly differentially expressed genes (defined as log2FC >|1| and p-adjust < 0.05) was compiled. DGE is performed using DESeq2 R package (version 1.142.1; Bioconductor 3.18; R 4.3) ( 42 ), where design formula contains batch-correction as RNA-Seq samples where sequenced in years of time. Permutation Test To assess the statistical significance of the observed association between gene expression changes and mutations in CTCFbs, we performed a permutation test with 1,000 iterations. In each iteration, sample labels were randomly permuted while preserving the original group sizes, and differential gene expression (DGE) analysis was conducted using DESeq2 for genes associated with each of the 284 chromatin loops. Given that the number of genes per loop varied, we quantified, per iteration, the number of loops containing at least one significantly differentially expressed gene (adjusted p value < 0.05). This yielded a null distribution of loop counts expected under the assumption of no association. The observed number of significant loop-DEG pairs was 17, while the maximum number observed in the permuted (null) distributions was 16. None of the 1,000 permutations produced a value equal to or greater than the observed. Using the standard empirical p -value formula with a continuity correction: $$\:\text{p}=\frac{\text{r}+1}{\text{N}+1}\:\text{}=\frac{0+1}{1000+1}\text{}=\frac{1}{1001}\text{}\approx\:0.001$$ Where: \(\:\text{p}\) is the empirical p -value, representing the probability of observing a value as extreme as or more extreme than the actual value under the null hypothesis. r is the number of permutations in which the test statistic (e.g., number of significant loop-DEG pairs) was greater than or equal to the observed value ( 17 ). N is the total number of permutations performed. this result indicates that the likelihood of observing 17 or more significant associations under the null hypothesis is less than 0.1%, supporting the non-random nature of the observed loop-DEG associations. Associations of CTCF binding site mutations and progression of cSCC: To investigate the role of CTCF binding site mutations in the progression of cSCC, we studied different cSCC sub-cohorts, i.e. cSCC primary tumors (non-metastasizing, n = 16; and metastasizing tumor, n = 15) and LNM (n = 41) samples. To identify loops with significantly different mutational CTCFbs between groups, a Chi-square test was applied when the expected frequency in each cell was ≥ 5, while Fisher's exact test was used when the expected frequency in any cell was < 5. A confusion matrix was constructed using the mutation data matrix, where '0' indicates the absence of a mutation event and '≥1' indicates the presence of a mutation event. Various comparison of groups was performed as reported in Table 3 . After getting the DE loops, genes list was extracted to perform the RNA-Seq DGE analysis. To further assess whether gene expression was associated with the mutation status of the CTCFbs of TAD loop in which each gene resides, we calculated Spearman’s rank correlation coefficient (ρ) between gene expression (treated as a continuous variable) and loop mutation status (binary: mutated vs. non-mutated). Spearman’s ρ is used for testing general associations between a continuous and a binary variable. In addition, we computed the area under the receiver operating characteristic curve (AUC) to evaluate the classification performance of gene expression in distinguishing between mutated and non-mutated loops. These analyses were performed in R, using the pROC package for AUC computation and base functions for correlation analysis. Table 4 Demographic and clinicopathological features of patients with lymph node metastases (LNM), primary metastatic (PM), and primary non-metastatic (PNM) cutaneous squamous cell carcinoma (cSCC). n refers to the number of samples (total = 72) derived from 59 patients, including 13 matched sample pairs between LNM and PM. Variable LNM n = 41 PM n = 16 PNM n = 15 Mean age, years (range) 70.8 (30–92) 70.6 (51–92) 75.7 (57–91) Sex, n (%) Female 4 (9.25) 4 ( 25 ) 3 ( 20 ) Male 37 (90.25) 12(75) 12 (80) Site of primary tumour, n (%) cheek NA 3 (18.75) 5 (33.33) neck NA 1 (6.25) 0 lip NA 1 (6.25) 0 eyebrow NA 1 (6.25) 0 ear NA 1 (6.25) 3 ( 20 ) temple NA 2 (12.5) 1 (6.66) nose NA 1 (6.25) 1 (6.66) postauricular NA 2 (12.5) 1 (6.66) scalp NA 4 ( 25 ) 2 (13.33) face NA 0 1 (6.66) pre-auricular NA 0 1 (6.66) Site of metastasis, n (%) Neck 21 (51.2) 12 (75) NA Parotid and neck 4 (9.75) 1 (6.25) NA parotid 15 (36.57) 3 (18.75) NA perifacial 1 (2.43) 0 NA T-stage at surgery, n (%) 0 or unknown 1 NA 3 (18.75) 6 ( 40 ) 2 NA 7 (43.75) 4 (26.67) 3 NA 3 (18.75) 5 (33.33) 4 3 (18.75) 0 N stage at surgery, n (%) 0 0 0 1 2 (4.88) 1 (6.25) 2 8 (19.51) 1 (6.25) 3 18 (43.9) 4 ( 25 ) unknown 13 (31.71) 10 (62.5) Histopathological grading, n (%) 1 2 (4.88) 1 (6.25) 1 (6.66) 2 5 (12.19) 2 (12.5) 10 (66.67) 3 24 (58.54) 9 (56.25) 4 (26.67) unknown 10 (24.39) 4 ( 25 ) 0 Lympho-vascular infiltration (LVI), n (%) No 16 7 14 Yes 15 8 1 Unknown 10 1 0 Peri-neural invasion (PNI), n (%) no 20 12 12 yes 12 3 3 unknown 9 1 0 Abbreviations TAD: Topologically Associated Domains CSCC: Cutaneous Squamous Cell Carcinoma CTCFbs: CTCF Binding Site LNM: Lymph node metastases PM: Primary metastatic PNM: Primary non-metastatic WGS : Whole Genome Sequencing RNA-Seq : RNA Sequencing DEG : Differentially Expressed Gene DGE : Differential Gene Expression TMB : Tumor Mutational Burden Declarations Ethics approval and consent to participate: This study was conducted with approval from the Institutional Human Research Ethics Committee (UOW/ISLHD HREC14/397). Consent for publication: “Not applicable” Availability of data and materials: (a) Code Availability : All the in-house script used for the secondary analyses are available from GitHub at https://github.com/amarinderthind/CTCF_Cancer_study (b) Data availability : Data is available on suitable request to the corresponding author. Competing interests: The authors declare that they have no competing interests. Funding: This work was funded by the Illawarra Cancer Carers, Cancer Institute NSW translational program grant 2020/TPG2081, National Health and Medical Research Council Project Grant APP1181179, and Tour de Cure RSP-00244-19/20. Authors' contributions: Conceptualization: AT (lead), BA (supporting), MR (supporting); Data curation: AT (equal), DS (equal); Formal analysis: AT (leading); Investigation: AT (lead), BA (supporting), MR (supporting), NS (supporting); Methodology: AT (lead), DS (supporting); Visualization: AT (lead), AKP (supporting), Writing - original draft: AT; Writing - editing: AT, MR, BA; Writing – review: AT, BA, MR, SM, NS, DS; Funding acquisition: BA (equal), MR (equal), RG (equal), JC (equal). Acknowledgements: We wish to acknowledge the National Computational Infrastructure (NCI) of Australia for providing computational resources that contributed to these results. References Thind AS, Ashford B, Strbenac D, Mitchell J, Lee J, Mueller SA, et al. Whole genome analysis reveals the genomic complexity in metastatic cutaneous squamous cell carcinoma. Frontiers in oncology. 2022;12:919118. Mueller SA, Gauthier M-EA, Ashford B, Gupta R, Gayevskiy V, Ch’ng S, et al. Mutational patterns in metastatic cutaneous squamous cell carcinoma. Journal of investigative dermatology. 2019;139(7):1449-58. e1. Bulen BJ, Khazanov NA, Hovelson DH, Lamb LE, Matrana M, Burkard ME, et al. Validation of Immunotherapy Response Score as predictive of pan-solid tumor anti-PD-1/PD-L1 benefit. Cancer research communications. 2023;3(7):1335-49. Kim TH, Abdullaev ZK, Smith AD, Ching KA, Loukinov DI, Green RD, et al. Analysis of the vertebrate insulator protein CTCF-binding sites in the human genome. Cell. 2007;128(6):1231-45. Long HS, Greenaway S, Powell G, Mallon A-M, Lindgren CM, Simon MM. Making sense of the linear genome, gene function and TADs. Epigenetics & Chromatin. 2022;15(1):4. Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485(7398):376-80. Yang J, Corces VG. Chromatin insulators: a role in nuclear organization and gene expression. Advances in cancer research. 2011;110:43-76. Poulos RC, Thoms JA, Guan YF, Unnikrishnan A, Pimanda JE, Wong JW. Functional mutations form at CTCF-cohesin binding sites in melanoma due to uneven nucleotide excision repair across the motif. Cell reports. 2016;17(11):2865-72. Fang C, Wang Z, Han C, Safgren SL, Helmin KA, Adelman ER, et al. Cancer-specific CTCF binding facilitates oncogenic transcriptional dysregulation. Genome biology. 2020;21:1-30. Gupta R, Strbenac D, Satgunaseelan L, Cheung VK-Y, Narayanappa H, Ashford B, et al. Comparing genomic landscapes of oral and cutaneous squamous cell carcinoma of the head and neck: quest for novel diagnostic markers. Modern Pathology. 2023;36(8):100190. Damaschke NA, Gawdzik J, Avilla M, Yang B, Svaren J, Roopra A, et al. CTCF loss mediates unique DNA hypermethylation landscapes in human cancers. Clinical Epigenetics. 2020;12:1-13. Wang S, Wang Y, Hao L, Chen B, Zhang J, Li X, et al. BOC targets SMO to regulate the Hedgehog pathway and promote proliferation, migration, and invasion of glioma cells. Brain Research Bulletin. 2024;216:111037. Dumitru AV, Țăpoi DA, Halcu G, Munteanu O, Dumitrascu D-I, Ceaușu MC, et al. The polyvalent role of CD30 for cancer diagnosis and treatment. Cells. 2023;12(13):1783. Gao Y, Shi H, Zhao H, Yao M, He Y, Jiang M, et al. Single‐cell transcriptomics identify TNFRSF1B as a novel T‐cell exhaustion marker for ovarian cancer. Clinical and Translational Medicine. 2023;13(9):e1416. Van der Weyden C, Pileri S, Feldman A, Whisstock J, Prince H. Understanding CD30 biology and therapeutic targeting: a historical perspective providing insight into future directions. Blood cancer journal. 2017;7(9):e603-e. Ling Q, Li F, Zhang X, Mao S, Lin X, Pan J, et al. MAP4K1 functions as a tumor promotor and drug mediator for AML via modulation of DNA damage/repair system and MAPK pathway. EBioMedicine. 2021;69. Chan L-K, Ho DW-H, Kam CS, Chiu EY-T, Lo IL-O, Yau DT-W, et al. RSK2-inactivating mutations potentiate MAPK signaling and support cholesterol metabolism in hepatocellular carcinoma. Journal of Hepatology. 2021;74(2):360-71. Zheng K, Yao S, Yao W, Li Q, Wang Y, Zhang L, et al. Association between RSK2 and clinical indexes of primary breast cancer: a meta-analysis based on mRNA microarray data. Frontiers in genetics. 2021;12:770134. Shi J-l, Fu L, Ang Q, Wang G-j, Zhu J, Wang W-d. Overexpression of ATP1B1 predicts an adverse prognosis in cytogenetically normal acute myeloid leukemia. Oncotarget. 2015;7(3):2585. Acconcia F. Evaluation of the sensitivity of breast cancer cell lines to cardiac glycosides unveils atp1b3 as a possible biomarker for the personalized treatment of erα expressing breast cancers. International journal of molecular sciences. 2022;23(19):11102. Thomas DD, Wink DA. NOS2 as an emergent player in progression of cancer. Mary Ann Liebert, Inc. 140 Huguenot Street, 3rd Floor New Rochelle, NY 10801 USA; 2017. p. 963-5. Coutinho LL, Femino EL, Gonzalez AL, Moffat RL, Heinz WF, Cheng RY, et al. NOS2 and COX-2 Co-expression promotes cancer progression: a potential target for developing agents to prevent or treat highly aggressive breast cancer. International Journal of Molecular Sciences. 2024;25(11):6103. He W, Wang Q, Xu J, Xu X, Padilla MT, Ren G, et al. Attenuation of TNFSF10/TRAIL-induced apoptosis by an autophagic survival pathway involving TRAF2-and RIPK1/RIP1-mediated MAPK8/JNK activation. Autophagy. 2012;8(12):1811-21. Yi JH, Kim SJ, Kim WS. Brentuximab vedotin: clinical updates and practical guidance. Blood research. 2017;52(4):243-53. Wang X, Wang Q, Wang H, Cai G, An Y, Liu P, et al. Small protein ERSP encoded by LINC02870 promotes triple negative breast cancer progression via IRE1α/XBP1s activation. Cell Death & Differentiation. 2025:1-12. Guo M, Zhuang H, Huang J, Shao X, Bai N, Li M, et al. LINC02870 facilitates SNAIL translation to promote hepatocellular carcinoma progression. Molecular and Cellular Biochemistry. 2023;478(9):1899-914. Urbaniak A, Jablonska K, Podhorska-Okolow M, Ugorski M, Dziegiel P. Prolactin-induced protein (PIP)-characterization and role in breast cancer progression. American journal of cancer research. 2018;8(11):2150. Zou Y, Cao C, Wang Y, Zhou Y, Yao S, Zhang L, et al. Multi-omics consensus portfolio to refine the classification of lung adenocarcinoma with prognostic stratification, tumor microenvironment, and unique sensitivity to first-line therapies. Translational Lung Cancer Research. 2022;11(11):2243. Han WJ, He P. A novel tumor microenvironment-related gene signature with immune features for prognosis of lung squamous cell carcinoma. Journal of Cancer Research and Clinical Oncology. 2023;149(14):13137-54. Di Iorio P, Beggiato S, Ronci M, Nedel C, Tasca C, Zuccarini M. Unfolding new roles for guanine-based purines and their metabolizing enzymes in cancer and aging disorders. Frontiers in Pharmacology. 2021;12:653549. Mason JM, Xu H-P, Rao SK, Leask A, Barcia M, Shan J, et al. Lefty contributes to the remodeling of extracellular matrix by inhibition of connective tissue growth factor and collagen mRNA expression and increased proteolytic activity in a fibrosarcoma model. Journal of Biological Chemistry. 2002;277(1):407-15. Gao X, Cai Y, An R. miR-215 promotes epithelial to mesenchymal transition and proliferation by regulating LEFTY2 in endometrial cancer. International journal of molecular medicine. 2018;42(3):1229-36. Do C, Skok JA. Factors that determine cell type–specific CTCF binding in health and disease. Current Opinion in Genetics & Development. 2024;88:102244. Thind AS, Monga I, Thakur PK, Kumari P, Dindhoria K, Krzak M, et al. Demystifying emerging bulk RNA-Seq applications: the application and utility of bioinformatic methodology. Briefings in bioinformatics. 2021;22(6):bbab259. Krueger F. Trim Galore!: A wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files, with extra functionality for RRBS data. Babraham Institute. 2015. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15-21. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923-30. Rao SS, Huntley MH, Durand NC, Stamenova EK, Bochkov ID, Robinson JT, et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell. 2014;159(7):1665-80. Consortium EP. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489(7414):57. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841-2. Durinck S, Spellman PT, Birney E, Huber W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nature protocols. 2009;4(8):1184-91. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology. 2014;15:1-21. Additional Declarations The authors declare no competing interests. Supplementary Files supplementryfuguresv2.docx SupplementaryTable1overlappingregions.csv SupplementaryTable2correlatinoAUCPNMPM.csv SupplementaryTable3Clinical.xlsx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-6844715","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":468015891,"identity":"5d299087-700a-4be4-9480-4271db124b9d","order_by":0,"name":"Amarinder Thind","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Amarinder","middleName":"","lastName":"Thind","suffix":""},{"id":468015892,"identity":"67083329-e183-42e3-9f7e-88ebae351141","order_by":1,"name":"Bruce Ashford","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAzElEQVRIiWNgGAWjYBACAyidwA+hmUnQItlAshaDA8RqMZc+Y/jhQ01dnvGN5GMSDBXWiQ3sZwzwarHsyzGWnHHscLHZjbQ0CYYz6YkNPDn4tRic4TFj5mE7kLjtRo7ZDca2w4kNDMRo+fOvLnHzDJCWf0At/G+I0MLYxpy4QQKkpQGoRYKALZY9bMWSvX2HE2eceZb+I+FYunGbxLMCvFrMeZg3fvjxrS6xvz35sMGHGmvZfv7kDXi1oIIEIGYjQf0oGAWjYBSMAhwAAJ0XRbJAYevaAAAAAElFTkSuQmCC","orcid":"","institution":"","correspondingAuthor":true,"prefix":"","firstName":"Bruce","middleName":"","lastName":"Ashford","suffix":""},{"id":468015893,"identity":"95ea2e95-f781-420c-9faf-704ab255fde8","order_by":2,"name":"Nicholas Shannon","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Nicholas","middleName":"","lastName":"Shannon","suffix":""},{"id":468015894,"identity":"7ca54cc7-8b3d-4cea-8082-1f4d32c83fdb","order_by":3,"name":"Dario Strbenac","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Dario","middleName":"","lastName":"Strbenac","suffix":""},{"id":468015895,"identity":"19222540-52a8-4236-8cef-ef58bc61fc7e","order_by":4,"name":"Ann Katrin Piper","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Ann","middleName":"Katrin","lastName":"Piper","suffix":""},{"id":468015896,"identity":"f419cbe9-2db3-42d6-9b51-ed79a8b34dea","order_by":5,"name":"Jenny Mitchell","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Jenny","middleName":"","lastName":"Mitchell","suffix":""},{"id":468015897,"identity":"d623facd-687c-4813-97a1-a98b90dc5298","order_by":6,"name":"Jonathan Clark","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Jonathan","middleName":"","lastName":"Clark","suffix":""},{"id":468015898,"identity":"01f157d4-de5c-46c0-acdf-dcfa056311ba","order_by":7,"name":"Ruta Gupta","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Ruta","middleName":"","lastName":"Gupta","suffix":""},{"id":468015899,"identity":"c3d902b9-f70e-4cd3-a453-e59b95a21151","order_by":8,"name":"Simon Mueller","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Simon","middleName":"","lastName":"Mueller","suffix":""},{"id":468015900,"identity":"03ebaffc-1604-42cd-afc7-b4763a25887b","order_by":9,"name":"Marie Ranson","email":"","orcid":"","institution":"","correspondingAuthor":false,"prefix":"","firstName":"Marie","middleName":"","lastName":"Ranson","suffix":""}],"badges":[],"createdAt":"2025-06-07 21:49:50","currentVersionCode":1,"declarations":{"humanSubjects":true,"vertebrateSubjects":false,"conflictsOfInterestStatement":false,"humanSubjectEthicalGuidelines":true,"humanSubjectConsent":true,"humanSubjectClinicalTrial":false,"humanSubjectCaseReport":false,"vertebrateSubjectEthicalGuidelines":false},"doi":"10.21203/rs.3.rs-6844715/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-6844715/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":85276826,"identity":"487ae31d-8a6d-4b27-9eb1-81ffbc6028b4","added_by":"auto","created_at":"2025-06-24 07:28:49","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":609204,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIntegrated multiomics analysis of CTCFbs in cSCC.\u003c/strong\u003e \u003cstrong\u003e(A)\u0026nbsp;\u003c/strong\u003eCTCFbs were identified in normal human epidermal keratinocytes (NHEKs) via public ChIP-seq data and motif-based analysis. \u003cstrong\u003e(B)\u003c/strong\u003e Public Hi-C data from NHEK was used to extract TADs with convergent CTCF anchors. \u003cstrong\u003e(C)\u003c/strong\u003e Schematic representation of a chromatin loop regulated by CTCF. Loss of CTCF binding can lead to erosion of the chromatin loop and reduced promoter–enhancer contacts, resulting in altered gene expression. \u003cstrong\u003e(D)\u003c/strong\u003e The TAD loops containing CTCF motifs were filtered via the NHEK ChIP-seq data to retain loops with confirmed CTCF binding. \u003cstrong\u003e(E)\u003c/strong\u003e Whole-genome sequencing (WGS) data from 72 cSCC samples (tumor and matched blood) were analysed to identify somatic SNV mutations within the CTCFbs of the filtered loops. \u003cstrong\u003e(F)\u003c/strong\u003e Distribution of mutation burden across different genomic features, highlighting higher mutation rates in CTCFbs. \u003cstrong\u003e(G)\u003c/strong\u003e Integration of RNA-Seq data from 54 samples (matched with WGS) enabled differential gene expression (DGE) analysis via DESeq2 to assess the impact of CTCFbs mutations on nearby gene expression (±1 kb from the binding site). \u003cstrong\u003e(H)\u003c/strong\u003e A permutation test (1,000 iterations) was performed to evaluate whether the observed number of loops with significant differentially expressed genes (DEGs) could arise by chance. \u003cstrong\u003e(I)\u003c/strong\u003e Finally, the role of CTCFbs mutations in cSCC progression was explored by comparing the mutation profiles across primary nonmetastasizing (PNM, n=15), primary metastasizing (PM, n=16), and lymph node metastasis (LNM, n=41) samples. Statistical tests (Chi-square/Fisher’s exact) were used to identify loops uniquely mutated in specific groups and their associated DEGs.\u0026nbsp;Created in\u0026nbsp;https://BioRender.com\u003c/p\u003e","description":"","filename":"floatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/81841a330e04679ab4b3f25a.png"},{"id":85275914,"identity":"d712c053-0d8b-42d1-8a65-a1df83c125b9","added_by":"auto","created_at":"2025-06-24 07:20:51","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":165224,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDistribution of tumor mutation burden (TMB) across cohort. (A)\u003c/strong\u003e Boxplot representation of the TMB distribution across all samples. The boxplot highlights the interquartile range with the red triangle and blue square representing the mean and median TMB, respectively\u003cstrong\u003e. (B)\u003c/strong\u003e Mutation density per region (per Mb) across all samples, including CTCF. (\u003cstrong\u003eC)\u003c/strong\u003e The mean mutation density (per Mb) at each position calculated across 72 samples and 903 loops (1,806 regions, representing both sides of the loops) within ±1 kb from the CTCFbs centre, with the central position adjusted to 0.\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/1943b248449f7c8e3b7386f9.png"},{"id":85275900,"identity":"41df1bd5-3773-42d2-b928-b9979f385262","added_by":"auto","created_at":"2025-06-24 07:20:50","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":184516,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCTCFbs with motif and substitutions (overall cohort): (A)\u003c/strong\u003e Shows the substitution density per Mb at the center of the CTCFbs (±20 bp) of 5 sites. Mutations are categorized by substitution type (e.g., C\u0026gt;T, T\u0026gt;A) and normalized to account for the total number of samples (72) and the total number of loops (903). The x-axis represents the position relative to the center of the CTCF binding site. \u003cstrong\u003e(B)\u003c/strong\u003eSimilar to A this\u003cstrong\u003e \u003c/strong\u003eshows the substitution density per Mb at the center of the CTCFbs (±20 bp) of 3 sites across all samples.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/7d896daf4adb5cba49a1bd02.png"},{"id":85275893,"identity":"dd121557-a97b-4d41-a521-e502d3da4f04","added_by":"auto","created_at":"2025-06-24 07:20:50","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":84337,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eMutation overlaps and mutational density in CTCF-bound regions from 72 cSCC samples (all stages).\u003c/strong\u003e This bar chart depicts the number of mutations shared between CTCFbs and various genomic regions (UTRs, CDSs, lncRNAs, and promoters). The box plot shows the distribution of mutations per megabase, the median, quartiles, and potential outliers of the TMB values for each category. The numbers above each box indicate the total number of mutations observed in that category.\u003c/p\u003e","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/581cb366f35a057fe703eaf7.png"},{"id":85275904,"identity":"6646c9a8-a45f-40ff-bc50-4dc5c0eef183","added_by":"auto","created_at":"2025-06-24 07:20:50","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":331071,"visible":true,"origin":"","legend":"\u003cp\u003eComparative Analysis of CTCF Binding Site Mutations and Gene Expression between Primary Metastatic (PM) and Primary NonMetastatic (PNM) cSCC Tumors. \u003cstrong\u003e(A) Plot of loops with different CTCF binding site (CTCFbs) mutations\u003c/strong\u003e presents loops with notable differences in CTCFbs mutations between the \u003cstrong\u003ePM\u003c/strong\u003e (n=15) and \u003cstrong\u003ePNM\u003c/strong\u003e(n=16) groups. The count indicates the number of samples with mutations in each group for each loop. \u003cstrong\u003e(B) Table of genes associated with differentially mutated loops.\u003c/strong\u003e This panel displays the genes associated with the loops listed in panel A. Genes were identified within or near (1 kb flanking) the CTCFbs motifs. Only genes with assigned hgnc symbol are presented here. \u0026nbsp;\u003cstrong\u003e(C) Volcano plot of differential gene expression.\u003c/strong\u003e This panel shows a volcano plot illustrating the differential gene expression analysis for genes associated with the differentially mutated loops. The x-axis represents the log2 fold change in gene expression between the \u003cstrong\u003ePM\u003c/strong\u003e and \u003cstrong\u003ePNM\u003c/strong\u003e groups, and the y-axis represents the -log10 p-value. The red points indicate genes with significant differential expression (HHIPL2 and GDA). The blue dotted lines represent the thresholds for significance (log2-fold change \u0026gt; 1, p-value \u0026lt; 0.01).\u003c/p\u003e","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/6e4c2f67eb9cb6cf2a7eced1.png"},{"id":85275899,"identity":"acc1f326-c2a5-4021-80e5-a4ca4b212756","added_by":"auto","created_at":"2025-06-24 07:20:50","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":101716,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSpearman correlation vs. AUC.\u003c/strong\u003e This scatter plot illustrates the relationship between the Spearman correlation coefficient of individual genes expression with a binary outcome (mutated/non-mutated) and the AUC achieved when using each gene expression alone is used to predict that outcome. Each point represents a different gene, labelled with its identifier. The color of each point indicates the predictive power of the corresponding feature on the basis of its AUC value: green signifies strong predictive power (high AUC), orange indicates fair-good predictive power (moderate AUC), and red denotes poor predictive power (low AUC). The vertical dashed line at a Spearman correlation of 0 highlights features with no linear correlation, whereas the horizontal dashed line at an AUC of 0.7 represents a common threshold for acceptable predictive performance. This visualization helps assess the individual utility of features for classification tasks, revealing that strong correlations (both positive and negative) tend to be associated with higher AUC values, while weak correlations often correspond to lower AUC values.\u003c/p\u003e","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/1dc2ff80ffba36cd6253799b.png"},{"id":85278330,"identity":"4d750594-03c2-4812-acdf-d1ec3cf5fcc4","added_by":"auto","created_at":"2025-06-24 07:44:52","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":3419689,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/1bcaa4cd-822e-49f0-9f5d-1a300a5cd637.pdf"},{"id":85275890,"identity":"25d6d198-4f8f-4328-959f-1cf1097f692f","added_by":"auto","created_at":"2025-06-24 07:20:50","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":717125,"visible":true,"origin":"","legend":"","description":"","filename":"supplementryfuguresv2.docx","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/43fd48fac0b3715e711d4e61.docx"},{"id":85275913,"identity":"e6ba9099-bc11-41ff-9ac6-7e290f7bdc02","added_by":"auto","created_at":"2025-06-24 07:20:50","extension":"csv","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":10186,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTable1overlappingregions.csv","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/5269c81e6da65cbdd5308634.csv"},{"id":85275888,"identity":"895008cb-69b8-48e0-a6c3-4bcb6df26ae9","added_by":"auto","created_at":"2025-06-24 07:20:49","extension":"csv","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":3666,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTable2correlatinoAUCPNMPM.csv","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/faa6146d3c2e39658a82898d.csv"},{"id":85276828,"identity":"b1e8a4d9-c2de-46e6-8135-d096e4e05787","added_by":"auto","created_at":"2025-06-24 07:28:50","extension":"xlsx","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":16300,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTable3Clinical.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-6844715/v1/bc9e833520afebb257d6eca4.xlsx"}],"financialInterests":"The authors declare no competing interests.","formattedTitle":"\u003cp\u003e\u003cstrong\u003eCTCF Binding Site Mutations: Linking Topologically Associated Domains Dysregulation to Cutaneous Squamous Cell Carcinoma Progression\u003c/strong\u003e\u003c/p\u003e","fulltext":[{"header":"Background","content":"\u003cp\u003eCutaneous squamous cell carcinoma (cSCC) is the most common malignancy with a propensity for aggressive metastasis. Resectable advanced disease requires surgery and adjuvant radiation to achieve control. This malignancy is characterized by an extremely high tumor mutation burden (TMB) across the genome (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e). Understanding the genomic events driving cSCC progression, especially the transition from localized to metastatic disease, is crucial for improving patient outcomes. Previous research has underscored the potential utility of using TMB to predict immunotherapy response and as a vital tool for understanding tumor evolution (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e). However, the specific genomic regions and regulatory elements contributing to TMB variability in highly mutated tumours like cSCC remain poorly defined. Central to this investigation is CTCF, a protein that binds to the core genome sequence CCCTC, known as CTCF binding sites (CTCFbs) (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e). CTCF plays diverse roles, including transcriptional regulation, insulation (blocking enhancer-promoter interactions), chromatin architecture regulation, and RNA splicing. For example, CTCF organizes chromatin into topologically associated domains (TADs), creating self-interacting genomic regions that facilitate interactions between regulatory elements and genes (\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e). TADs are fundamental to gene regulation, as they compartmentalize the genome, restricting interactions to within the TAD and preventing the spread of heterochromatin's silencing effects (\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e, \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eWe have previously shown that CTCFbs are prone to UV-induced mutations in 15 cSCC lymph node metastases. Interestingly, cSCC exhibits CTCFbs mutation at a rate seven times that seen in melanoma, also a UV-mediated cancer (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e). Indeed, our previous work identified 1,979 genes within 422 affected TADs, including both tumor suppressor genes and oncogenes. These findings suggest that mutations in CTCFbs could significantly disrupt regulatory processes, leading to gene expression abnormalities associated with tumor progression. While mutations in CTCFbs have been observed in various cancers (\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e), the impact of CTCFbs mutations on gene expression during cSCC progression has not been investigated.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eIn this study, we present a comprehensive analysis of mutation density and substitution patterns at CTCFbs within TADs in a clinically stage-defined cohort of cSCC samples. This includes high-risk non-metastasizing primary tumors (PNM), metastatic primary tumors (PM) and lymph node metastases (LNM). We compared this to the distribution of mutations among various genomic regions across the entire genome and examined the impact of CTCFbs variants on gene expression (by integrating whole genome and RNA sequencing) with the aim of identifying potential drivers of metastasis (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e).\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\"\u003e\n \u003ch2\u003eTMB and Substitution Types of the Overall Cohort\u003c/h2\u003e\n \u003cp\u003eThe mean genomic TMB across all samples in this study was 199.46 mutations/Megabase (Mb), with a median of 145.56 mutations/Mb. The TMB values ranged from 7.59 to 999.72 SNVs/Mb, reflecting a broad and heterogeneous distribution (Fig. 2A). The standard deviation of 200.85 mutations/Mb further underscores the high intersample variability in mutation density. Among the samples, CSCC_0126-M1 presented the highest TMB (999.72 mutations/Mb), whereas CSCC_0157-P1 presented one of the lowest values (7.58 mutations/Mb). This variation suggests potential biological differences in mutation accumulation across tumor samples.\u003c/p\u003e\n\u003c/div\u003e\n\u003ch3\u003eMutational density across Genomic Regions and their Background\u003c/h3\u003e\n\u003cp\u003eFirst, a comparison of different noncoding/coding genomic regions with each other was conducted to investigate whether certain known genomic regions present higher or lower mutational densities relative to each other (Fig. 2B). The average mutation density at the CTCFbs, which is based on a 20 bp motif, was greatest with a mutational rate of approximately 1100 mutations per Mb. In contrast, the average mutation rates for the other five regions ranged between 170 and 250 per Mb, with long noncoding RNAs (lncRNAs) regions exhibiting the second highest average mutation rate at ~ 260 per Mb (Fig. 2B). Second, we investigated the surrounding regions of each genomic region at various scaling factors to determine whether the TMB is elevated significantly for a specific region than in its surrounding region (Supplementary Fig. 1). The genomic regions tested included the 3' UTR, 5' UTR, promoter regions, lncRNA regions, CDS, and others. When these regions were compared with their surroundings, a sharp decrease in mutation density was observed at both ends of the CTCFbs (Fig. 2C; Supplementary Fig.\u0026nbsp;1). For the other regions, the decrease was more gradual, and the greater the distance from the region was, the greater the reduction in mutation rate, except for CTCFbs, which showed a distinct pattern (Supplementary Fig.\u0026nbsp;1). This observation may indicate unique regulatory or structural features at CTCFbs that influence mutation accumulation.\u003c/p\u003e\n\u003cp\u003eA notably high mean mutation density near the CTCFbs centre (position 0) indicated differences in mutation accumulation, suggesting a specific mutation profile associated with these binding sites (Fig. 2C). Additionally, the substitution rates per Mb at the CTCFbs (± 20 bp) highlighted the distribution of different mutation types (Fig. 3). C \u0026gt; T substitutions are more prevalent at the core of CTCFbs (positions 3 and 4 on Anchor 5 and positions − 4 and − 5 for Anchor 3).\u003c/p\u003e\n\u003ch3\u003eCTCFbs and Regulatory Element Interactions\u003c/h3\u003e\n\u003cp\u003eCTCFbs that overlap with other genomic elements present a greater mutational density than do CTCFbs regions alone (Fig.\u0026nbsp;4). Among these, CTCFbs overlapping with lncRNAs display the highest TMB, surpassing those overlapping with coding sequences (CDSs), untranslated regions (UTRs), or promoter regions (Supplementary Table\u0026nbsp;1). This enrichment suggests that CTCF-bound coding and regulatory elements are particularly susceptible to tumor-associated mutations, potentially affecting gene expression, protein function, and chromatin architecture.\u003c/p\u003e\n\u003cdiv\u003eAdditionally, while the density of CTCF motif mutations is typically averaged over a 20 bp window, certain positions show markedly higher mutation rates. Specifically, as shown in Fig. 3, positions 3 and 4 on the CTCFbs motif of Anchor 5 and positions − 4 and − 5 on the motif of Anchor 3, show mutation densities ranging from approximately 5,000 to 7,000 mutations per Mb, indicating potential hotspots for functional disruption.\u003c/div\u003e\n\u003cp\u003e\u003cstrong\u003eCTCF loops across the entire cohort\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe lengths of the CTCF loops varied significantly, ranging from 48.9 kb to 1904 kb. A total of 667 (out of 903 TAD loops with CTCFbs in both anchor regions) were mutated in 72 samples. The number of mutated samples per loop ranges from 1 to 29, with frequently mutated loops could represent key regulatory elements or regions with critical roles in the pathogenesis of the disease.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eSubstitution types\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis analysis revealed that C \u0026gt; T substitutions were dominant across most samples, with a few samples showing a notable proportion of T \u0026gt; G and T \u0026gt; C substitutions (Supplementary Fig.\u0026nbsp;2). Interestingly, two samples presented very few C \u0026gt; T substitutions, suggesting potential sample-specific mutational patterns or underlying biological factors influencing mutation processes. The predominance of C \u0026gt; T substitutions could reflect mechanisms such as deamination of methylated cytosines, a common mutational signature in skin cancers.\u003c/p\u003e\n\u003ch3\u003eOverall mutation density and substitution types across cSCC clinical groups\u003c/h3\u003e\n\u003cp\u003eMutation density per sample and per group was compared among PNM, PM, and LNM groups using the Kruskal-Wallis test, which yielded a non-significant result, indicating no statistically significant differences. To further explore whether mutation patterns differ between clinical groups, we analysed the mean mutation density of substitution types across groups (Supplementary Fig.\u0026nbsp;3). The relative contribution of each substitution type varied slightly among the groups. In general PNM had a lower CTCFbs mutation density compared to PM and LNM.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDifferences in CTCFbs Mutations\u003c/strong\u003e: Analysis of mutation patterns within the CTCFbs regions of TADs revealed \u003cstrong\u003eno substantial differences\u003c/strong\u003e across clinical groups. Supplementary Fig. 4 displays the distribution of mutation density per mb (adjusted for sample size and number of regions) across CTCFbs motifs. CTCFbs in Anchor5 (top panels) showed a enrichment of C \u0026gt; T transitions, particularly at positions 3 and 4, corresponding to the \"G and G/T\" sites of the motif (middle panel). In contrast, the CTCFbs in Anchor3 (bottom panels) exhibited a notable peak of C \u0026gt; T transitions at positions − 4 and − 3 across all the groups.\u003c/p\u003e\n\u003ch3\u003eAssociation between CTCF binding site mutations and the expression of harbouring and nearby genes\u003c/h3\u003e\n\u003cdiv\u003e\n \u003cp\u003eAs described in the methods, we grouped samples with and without mutations for each loop using 54 shared WGS/RNA-Seq samples. Only 284 loops had mutations in at least three samples, which was set as the minimum group size required for differential gene expression (DGE) analysis. Notably, 17 of these 284 loops contained significantly differentially expressed genes (DEGs) when samples groups with and without mutations (regardless of clinical subgroup) were compared (Table 1). Out of these 17 loops, 3 overlaps with 3’UTR regions.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv\u003e\n \u003ctable id=\"Tab1\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv\u003eTable 1\u003c/div\u003e\n \u003cdiv\u003e\n \u003cp\u003eShows the 17 loops that harbour significant differentially expressed genes (shown) between samples with and without mutations regardless of cohort subgrouping. Differential expression was considered significant using a threshold of log₂FC ≥ 1 or ≤ − 1 and adjusted p-value (adj. P) \u0026lt; 0.05. Only DEGs meeting these criteria are shown.\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eLoop\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eNo. of Significant DEGs\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eDEGs\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_1057\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eENSG00000227726, ENSG00000279459, ENSG00000226627, ENSG00000280089, ENSG00000171671, ENSG00000286708\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_21\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e3\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eRPS27P20, ENSG00000281386, ENSG00000281097\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_28\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eTNFRSF8, TNFRSF1B (ENSG00000028137)\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_310\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eMT1L, ENSG00000205364\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_404\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eENSG00000287920, ENSG00000290062\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_364\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eBOC, ENSG00000288079\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_213\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eENSG00000285090, ENSG00000285964\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_682\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eLINC02783, ENSG00000159339\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_204\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003ePIP\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_336\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eENSG00000233633\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_75\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eRGS6\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_894\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eCFAP221\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_156\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eKALRN\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_26\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eCASQ2\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_1033\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eFAM78A\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_207\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eRAMP3\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eloop_697\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eCD163\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003c/table\u003e\n\u003c/div\u003e\n\u003cdiv\u003e\n \u003cp\u003eFurthermore, to test whether the association was not by random chance, we performed a permutations test with 1000 iterations, as described in the Methods, the observed number of loops associated with significant DEGs was 17, which exceeded the maximum number observed in the null distribution (16; Supplementary Fig. 5). This yielded an empirical \u003cem\u003ep\u003c/em\u003e-value \u0026lt; 0.001, supporting that the observed loop-DEG associations are unlikely to have occurred by random chance.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec8\"\u003e\n \u003ch2\u003eMutually Exclusive CTCFbs mutations between metastatic vs non-metastatic subgroups\u003c/h2\u003e\n \u003cp\u003eWe observed that certain CTCFbs were uniquely mutated in primary tumors that metastasized (PM, n = 15) compared with those that had not metastasized (PNM, n = 16) as shown in Fig. 5A.\u003c/p\u003e\n \u003cp\u003eSix loops are only mutated in the PNM group, whereas 4 are uniquely mutated to the PM group, suggesting that mutations at specific CTCFbs may be linked to the metastatic potential of cSCC. However, these mutations are not present in every sample within each clinical subgroup, highlighting the heterogeneity of mutation patterns across individual cases; for example, loop_45 is mutated in 5/16 (31%) PNM samples. Importantly, the genes associated with the mutated loops (Fig. 5B) are previously known to be linked to cancer progression.\u003c/p\u003e\n \u003cp\u003eWe further, investigated whether the loops with different CTCFbs mutations between subgroups groups (refer to the methods section) are associated with the differential expression of genes within the TAD loop being tested. For these specific groups, first we considered samples for which we had both WGS and RNA-Seq in our cohort PM (n = 14) vs PNM (n = 11). Using these two groups, we performed DGE analysis and only genes with a |log2FoldChange| \u0026gt; 1 and a p-value ≤ 0.01) were considered as significantly DEG. \u003cem\u003eGDA\u003c/em\u003e and \u003cem\u003eHHIPL2\u003c/em\u003e are observed as significant in this comparison.\u003c/p\u003e\n \u003cp\u003eUsing Spearman’s rank correlation, we identified genes whose expression levels were significantly associated with the mutation status of the loops in which they reside. Several genes exhibited strong positive or negative correlations (Fig. 6; Supplementary Table 2), indicating differential expression linked to loop mutations. To complement this, AUC values were calculated, revealing high classification performance for some genes in distinguishing between mutated and non-mutated loops based on expression patterns (for details refer to Method section). This is done for each loop (Fig. 5A) by focusing on a particular stage. \u003cem\u003eENSG00000259093\u003c/em\u003e, located within loop_70, showed a strong negative correlation (Spearman ρ = -0.51) and a high AUC of 0.92, indicating consistent downregulation in the presence of loop mutations. In contrast, \u003cem\u003eLYRM9\u003c/em\u003e in loop_95 exhibited a positive correlation (ρ = 0.58, AUC = 0.91), suggesting that upregulation was associated with mutated loops.\u003c/p\u003e\n \u003cp\u003eNotably, a subset of gene–loop pairs demonstrated exceptionally high classification performance. For example, \u003cem\u003eAKAP6\u003c/em\u003e (loop_1083) achieved a perfect classification (AUC = 1.0), whereas \u003cem\u003eLINC02068\u003c/em\u003e (loop_284) reached an AUC of 0.94 with a strong positive correlation (ρ = 0.60). These findings suggest that gene expression within certain loops may serve as an effective proxy or biomarker for structural disruption caused by somatic mutations.\u003c/p\u003e\n \u003cp\u003eInterestingly, some gene–loop pairs presented low correlations but moderate AUC values, which may indicate nonlinear or threshold-based expression responses, potentially mediated by epigenetic buffering or compensatory mechanisms.\u003c/p\u003e\n \u003cp\u003eWe also performed other group comparisons for dissimilarity in \u003cstrong\u003eCTCFbs\u003c/strong\u003e mutations; however, no significant biological conclusion was drawn from this analysis. Table 3 shows various comparisons and parameters used for these analyses and a list of significant DEGs found in each of those comparisons.\u003c/p\u003e\n \u003cdiv\u003e\n \u003ctable id=\"Tab2\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv\u003eTable 3\u003c/div\u003e\n \u003cdiv\u003e\n \u003cp\u003e\u003cstrong\u003eOther Group Comparisons for Dissimilarity in CTCFbs Mutations.\u003c/strong\u003e This table summarizes the results of additional group comparisons performed to assess the dissimilarity in CTCFbs mutations across various tumor stages. While no significant biological conclusions were drawn from these analyses, we report the comparison parameters used and the differentially expressed genes (DEGs) identified in each group. The comparisons included different tumor stages (PNM vs PM, primary vs Met, and PM vs Met) along with specific anchor region considerations, loop filter criteria, and the list of significant DEGs for each comparison.\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eComparison\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eSample Groups\u003c/p\u003e\n \u003cp\u003e(WGS/RNA-seq)\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eAnchor\u003c/p\u003e\n \u003cp\u003eRegions\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eLoop Filtering\u003c/p\u003e\n \u003cp\u003eCriteria\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eSignificant DEGs\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003ePNM vs PM\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e15/11 vs 16/14\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e5′ and 3′ anchors\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ep-value \u0026lt; 0.3 \u0026amp; (sumPNM ≤ 1 or sumPM ≤ 1)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eGDA, HHIPL2\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003ePrimaries (PNM + PM) vs LNM\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e31/25 vs 41/29\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e5′ and 3′ anchors\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ep-value \u0026lt; 0.3 \u0026amp; (sumPrimary ≤ 1 or sumLNM ≤ 1)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eLINC02884, LEFTY2, MARK2P13\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003ePM vs LNM\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e16/14 vs 41/29\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e5′ and 3′ anchors\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ep-value \u0026lt; 0.4 \u0026amp; (sumPM ≤ 1 or sumLNM ≤ 1)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eENSG00000261184, LINC02884, ENSG00000259093\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003c/table\u003e\n \u003c/div\u003e\n\u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eThis study investigated the impact of mutations at the CTCFbs of TADs on gene expression that lies within TADs and their role in cSCC cancer progression. In line with the very high TMB reported in our previous studies in cSCC (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e), the mean genome TMB across all samples was high compared with that of other cancers. We observed an exceptionally high mutation density at the core of the CTCF motif, with notable enrichment of C\u0026thinsp;\u0026gt;\u0026thinsp;T substitutions particularly. Mutations were significantly concentrated at the core nucleotides of the CTCF binding motif, suggesting that these regions are particularly vulnerable to DNA damage. This pattern aligns with known mutational mechanisms, including cytosine deamination at methylated CpG dinucleotides, which frequently results in C\u0026thinsp;\u0026gt;\u0026thinsp;T transitions\u0026mdash;a common signature in cancer genomes (\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e). The high mutation rate at CTCFbs\u0026mdash;especially those overlapping with regulatory elements like long non-coding RNAs, UTRs or coding sequences (Fig.\u0026nbsp;4)\u0026mdash;suggests potential disruption of local chromatin structure or enhancer-promoter insulation.\u003c/p\u003e \u003cp\u003eHowever, the relatively short lengths of the CTCFbs overlapping regulatory regions and the limited sample size in this study constrain our ability to draw strong statistical conclusions or biological interpretations. Despite these limitations, the observed trends reveal intriguing patterns that merit further investigation. We also examined whether the 17 CTCFbs-associated loops linked to significantly DEGs overlapped with other genomic regions. Interestingly, 3 of the 17 loops binding site motif overlap with 3\u0026rsquo;UTR regions (Supplementary-Table-1). These loops are loop_28, loop_894, loop_364 and contains cancer progression related genes such as \u003cem\u003eBOC, TNFRSF8 (CD30), and TNFRSF1B\u003c/em\u003e (\u003cspan additionalcitationids=\"CR13 CR14\" citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e). Theses overlap needs future more comprehensive analyses to validate and extend these preliminary observations.\u003c/p\u003e \u003cp\u003e \u003cb\u003eCTCF binding site mutations and gene dysregulation\u003c/b\u003e: Among the 17 DEG-associated loops, several harbor genes implicated in hallmark processes of cancer such as immune evasion, proliferation, and metastasis. \u003cem\u003eMAP4K1\u003c/em\u003e has been reported to act as either an oncogene or a tumor suppressor gene (TSG), depending on the signaling environment (\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e). Similarly, \u003cem\u003eRPS6KA3\u003c/em\u003e (also known as \u003cem\u003eRSK2\u003c/em\u003e) demonstrates dual roles: while often overexpressed in cancers, it can suppress proliferation under certain conditions (\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e, \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e). \u003cem\u003eATP1B1\u003c/em\u003e also acts as a tumor suppressor, with evidence supporting its role in inhibiting cancer cell proliferation and migration (\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e, \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e). In contrast, \u003cem\u003eNOS2\u003c/em\u003e displays mixed behavior, functioning either as a tumor suppressor or a tumor promoter depending on the biological context (\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e). \u003cem\u003eTNFSF10\u003c/em\u003e functions predominantly as a tumor suppressor, known for its strong pro-apoptotic activity in cancer cells (\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e). \u003cem\u003eTNFRSF8 (CD30)\u003c/em\u003e is a member of the tumor necrosis factor receptor superfamily and serves as a tumor marker. It is found on the surface of specific cells, including certain immune cells and cancer cells(\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e). \u003cem\u003eCD30\u003c/em\u003e is also the target of the FDA-approved therapeutic brentuximab vedotin (Adcetris)(\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e). \u003cem\u003eLINC02870\u003c/em\u003e promotes triple negative breast cancer (\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e) and hepatocellular carcinoma progression (\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e). \u003cem\u003ePIP\u003c/em\u003e gene expression decreases gradually with increasing stage and grade of breast cancer (\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e). These associations underscore the potential biological importance of mutations within the CTCFbs of TADs in modulating gene expression in cSCC.\u003c/p\u003e \u003cp\u003e \u003cstrong\u003eDifferential Loop Mutations and Metastatic Potential\u003c/strong\u003e \u003cp\u003eWhen PM and PNM tumors were compared, specific CTCF loops (e.g., loop_45) were found to be recurrently mutated in one group but not in the other. Interestingly, while mutations in a given loop were observed in up to ~\u0026thinsp;31% of samples within a group, suggesting a recurrent, yet heterogeneous pattern of loop disruption, the DEG analysis still revealed that \u003cem\u003eHHIPL2\u003c/em\u003e and \u003cem\u003eGDA\u003c/em\u003e were significantly altered in the PM vs PNM comparisons (regardless of loop mutation or not). \u003cem\u003eHHIPL2\u003c/em\u003e has been previously associated with lung acinar adenocarcinoma (\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e, \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e), further supporting its role in cancer metastasis. In addition, \u003cem\u003eGDA\u003c/em\u003e is known to play a direct role in skin carcinogenesis by interacting with several cytokines and growth factors (\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e). These findings suggest that CTCFbs loop mutations may serve as potential markers of metastatic potential, but that CTCFbs mutation alone may not be only driver of gene expression changes. However, mutations were not uniformly present across all samples in a group, which may reflect tumor heterogeneity or stochastic mutation patterns. In the comparison of primary (PM\u0026thinsp;+\u0026thinsp;PNM) vs LNM, \u003cem\u003eLEFTY2\u003c/em\u003e is found to be downregulated in LNM. \u003cem\u003eLEFTY2\u003c/em\u003e is a member of the \u003cem\u003eTGF-β\u003c/em\u003e superfamily that functions as a tumor suppressor by inhibiting epithelial\u0026ndash;mesenchymal transition (EMT) (\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e). Its downregulation or epigenetic silencing has been linked to increased stemness, invasion, and metastatic potential in several cancers, including endometrial cancer (\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e).\u003c/p\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003eFunctional impact vs mutational noise\u003c/b\u003e: Although many loops harboured mutations, only a small subset were associated with significant changes in gene expression. This may reflect a combination of factors including the non-functional nature of some mutations, compensatory mechanisms, or limitations of bulk RNA-seq in capturing cell-type-specific regulatory changes (\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e, \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e). Moreover, given that only\u0026thinsp;~\u0026thinsp;31% of samples in a group may harbor mutations in any given loop, the effect size on transcriptomic output may be diluted in population-level analyses.\u003c/p\u003e \u003cp\u003e \u003cb\u003eStatistical and Experimental Limitations\u003c/b\u003e: While permutation testing confirmed that the observed DEG-loop associations are unlikely to be due to random chance (empirical p\u0026thinsp;=\u0026thinsp;0.00), several limitations remain. The sample size, particularly in subgroup comparisons (e.g., PM vs. PNM), reduces the power for detecting subtle effects. Bulk RNA-seq data also limits the ability to infer regulatory consequences in specific cell populations, especially for lineage-specific genes. Furthermore, while our data suggest that mutations at CTCFbs can dysregulate gene expression, direct mechanistic validation (e.g., through ChIP-Seq, CRISPR editing, or 3D chromatin conformation assays) is necessary to confirm causality.\u003c/p\u003e \u003cp\u003e \u003cstrong\u003eBiological and Clinical Implications\u003c/strong\u003e \u003cp\u003eThis study supports the hypothesis that noncoding mutations, particularly those in architectural elements such as CTCFbs, can contribute to tumor progression through dysregulation of TAD loops and associated gene expression. Given the growing interest in targeting chromatin architecture and epigenetic dysregulation in cancer, these findings may inform new strategies for biomarker discovery or therapeutic targeting in cSCC.\u003c/p\u003e \u003c/p\u003e"},{"header":"Conclusion","content":"\u003cp\u003eIn this study, we reveal that CTCFbs are hotspots of mutation accumulation across all clinical stages of cSCC, exhibiting a markedly higher TMB compared to other coding and non-coding elements. Notably, mutations at CTCFbs were not uniformly distributed and showed distinct substitution patterns, particularly enriched for C \u0026gt; T transitions, consistent with models of UV-associated carcinogenesis.\u003c/p\u003e \u003cp\u003eThrough integrative analysis of matched RNA-Seq and WGS data, we identified a subset of CTCFbs loops where mutations were significantly associated with altered expression of genes located within or near the corresponding TAD. Among these, genes including \u003cem\u003eHHIPL2\u003c/em\u003e, \u003cem\u003eGDA\u003c/em\u003e, and \u003cem\u003eLINC02884\u003c/em\u003e—known to be involved in oncogenic pathways—were differentially expressed, linking CTCFbs mutations to functional consequences in cSCC pathogenesis. Permutation testing confirmed these associations were unlikely due to random chance.\u003c/p\u003e \u003cp\u003eFurthermore, we observed that certain CTCF-bound loops exhibited differential mutation patterns between PNM and PM, suggesting a potential role for CTCFbs mutations in promoting metastatic behaviour. Overall, our findings underscore the central role of CTCFbs and their genomic context in modulating mutational landscapes and gene expression in cSCC. By uncovering how disruptions at these critical regulatory hubs may contribute to cancer development and progression, this study offers new insights to inform biomarker discovery and therapeutic targeting strategies for skin cancers.\u003c/p\u003e "},{"header":"Methods","content":"\u003ch2\u003eSample Collection:\u003c/h2\u003e\u003cp\u003eThis study was conducted with approval from the Institutional Human Research Ethics Committee (HREC/15/RPAH/266). Patients with resectable metastatic cSCC and high-risk non metastasising primary tumors were prospectively identified by the treating surgeons prior to surgery. Clinicopathological data, including age, sex, extent of nodal metastases, histology, and immunosuppressive status, were collected (Table\u0026nbsp;\u003cspan refid=\"Tab3\" class=\"InternalRef\"\u003e4\u003c/span\u003e, Supplementary Table\u0026nbsp;3). Fresh tumor tissue from nodal metastases (n = 41) was harvested during surgery and immediately snap-frozen. A total of 72/75 cSCC, WGS samples were included in the analysis (refer to QC section), comprising primary tumors (both metastatic (PM (n = 16)) and non-metastatic (PNM (n = 15)), as well as lymph node metastases (LNM (n = 41)). The metastatic cohort included matched PM and LNM samples for 13 patients. Further, RNA sequencing (RNA-Seq) data were collected for 54/72 samples, representing both primary and metastatic stages of the tumor. The PNM group had to meet the following criteria: absence of metastases at \u0026gt; 24 months follow-up after resection of the primary or negative sentinel lymph node biopsy at time of resection or histologically negative neck dissection.\u003c/p\u003e\u003ch2\u003eWhole genome and total RNA sequencing and QC:\u003c/h2\u003e\u003cp\u003eTumor tissue sections were processed for DNA and RNA extraction (using Qiagen AllPrep DNA and RNA kits, Qiagen, Hilden, Germany) and for estimating tumor cellularity. Only samples with a tumor content of greater than 30% (range: 35–95%) proceeded to DNA quality control (QC). QC procedures included spectrophotometry (Nanodrop 2000, Thermo Fisher Scientific Inc.), and gel electrophoresis. WGS was performed by AGRF (Melbourne, Australia), Genome.One (Darlinghurst, Australia) on Illumina HiSeq X to a depth of ×30–45 for whole blood (germline DNA) and ×65–90 for tumor samples. on the Illumina NovaSeq 6000 platform (Illumina). The average sequencing coverage was 94.56× (range: 64–143) for tumor samples and 41.08× (range: 30–56) for blood samples. Of the 75 samples sequenced, 72 passed QC. The remaining 3 samples exhibited extreme GC bias. Initially, total RNA was QCed before sequencing was performed by spectrophotometer. Following sequencing, the raw RNA-Seq data underwent quality control using the bioinformatics tool FastQC (version 0.11.9; Andrews, 2010). Low-quality reads were removed using Trim Galore (version 0.4.5) (\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e). After quality filtering, a total of 54 RNA samples were retained for downstream analysis.\u003c/p\u003e\u003ch2\u003eWGS and RNA-Seq pre-processing\u003c/h2\u003e\u003ch2\u003eSomatic Variant Analysis Using DRAGEN (v4.3.6) on ICA v2:\u003c/h2\u003e\u003cp\u003eSomatic variant analysis was performed using DRAGEN pipeline version 4.3.6 on the Illumina Connected Analytics (ICA) v2 environment using in house shell scripting. Tumor-normal paired analyses were carried out in three stages: alignment, variant calling, and integrative somatic analysis. Further information on the somatic calling is available at \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://help.dragen.illumina.com/product-guides/dragen-v4.3/dragen-dna-pipeline/small-variant-calling/somatic-mode\u003c/span\u003e\u003cspan address=\"https://help.dragen.illumina.com/product-guides/dragen-v4.3/dragen-dna-pipeline/small-variant-calling/somatic-mode\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. \u003cb\u003e(a) Alignment of Normal Samples\u003c/b\u003e: FASTQ files from normal samples were aligned using the DRAGEN Germline pipeline (v4.3.6) with the GRCh38 \u003cb\u003e(\u003c/b\u003ehg38-alt_masked.graph.cnv.hla.rna_v4.tar.gz\u003cb\u003e)\u003c/b\u003e reference genome. Small variant calling was enabled to generate a germline VCF, which was used as input for downstream somatic calling. \u003cb\u003e(b) Alignment of Tumor Samples\u003c/b\u003e: Tumor sample FASTQ files were aligned using the DRAGEN Somatic pipeline (v4.3.6), using the same reference as for the normal samples. \u003cb\u003e(c) Tumor-Normal Paired Analysis\u003c/b\u003e: Paired analysis was performed to identify somatic SNVs. Inputs included: Aligned tumor and normal BAM and BAI files, Germline SNV VCF from the normal alignment step, A systematic noise BED file (downloaded from:\u003c/p\u003e\u003cp\u003e \u003cspan class=\"ExternalRef\"\u003e \u003cspan class=\"RefSource\"\u003ehttps://webdata.illumina.com/downloads/software/dragen/resource-files/sv-systematic-noise-baseline-collection-3.0.0.tar\u003c/span\u003e \u003cspan address=\"https://webdata.illumina.com/downloads/software/dragen/resource-files/sv-systematic-noise-baseline-collection-3.0.0.tar\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e \u003c/span\u003e). Post Somatic Calling Filtering is done as default and hard filtered .\u003cem\u003evcf\u003c/em\u003e files are used for further analysis.\u003c/p\u003e\u003ch2\u003eRNA-Seq data processing:\u003c/h2\u003e\u003cp\u003eRNA sequence reads were mapped with STAR version 2.7.10a (\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e) onto the GRCh38. On average, for each sample 70% of the reads aligned to the reference genome (average \u0026gt; 87\u0026nbsp;million reads/sample). Transcript abundance was measured in terms of read counts using the same annotation file used for the transcriptome assembly, leveraging the featureCounts (\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e), R Bioconductor package, default parameters. The count matrix was used as input for gene differential expression analysis.\u003c/p\u003e\u003cp\u003e \u003cstrong\u003eTAD loops with CTCF motifs and mutations at CTCFbs\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eTAD loops containing CTCF motifs were identified in NHEK tissue using the same methodology as described in Mueller et al., 2019 (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e). This involved utilizing chromosome conformation capture (Hi-C) TAD maps from NHEK (\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e) along with chromatin immunoprecipitation sequencing data from ENCODE (The ENCODE Project Consortium, 2012) (\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e). A 20-bp motif and a 20-bp position-weighted matrix (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e) were applied to select binding sites with high CTCF-binding probability. TADs were filtered as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e. TADs were excluded if either anchor region contained more than one CTCF binding motif, or if the CTCF motifs were not in a convergent orientation, as this configuration is most strongly associated with CTCF binding (\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e). The genomic coordinates defining each TAD were determined as the 3′-end of the upstream motif and the 5′-end of the downstream motif. 903 TAAD loops with CTCFbs in both anchor regions were identified.\u003c/p\u003e\u003ch2\u003eIdentification of Mutated CTCFbs:\u003c/h2\u003e\u003cp\u003eBinding sites genomic co-ordinates of all 903TAD-loops were intersected with our cohort’s WGS somatic mutations data using \u003cem\u003ebedtools intersect\u003c/em\u003e (version V2.31.1;(\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e)) to obtain a mutational data matrix (samples vs loops) that contains information’s of mutated loops in each sample. Further, Genes lying with-TADs and 1000bp outside either side of the CTCFb-motif in the TADs were extracted using \u003cem\u003eBiomart R package\u003c/em\u003e (\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e).\u003c/p\u003e\u003ch2\u003eMutations in CTCF binding motif vs background:\u003c/h2\u003e\u003cp\u003eTo assess the mutational enrichment in CTCF binding motifs, we defined the motif region as ± 10 bp around the center of each CTCF binding site. The background was defined as the flanking regions of ± 1kb relative to the binding site center. CTCFbs and their flanking 1kb regions were defined and merged into a \u003cem\u003eGenomicRanges\u003c/em\u003e object. Mutation data was loaded and represented as a \u003cem\u003eGenomicRanges object\u003c/em\u003e. Overlaps between mutations and the 2kb CTCF regions (± 1kb) were identified, and mutation positions were normalized to a ± 1kb scale centered on the CTCF binding site. For each position within the 2kb region, the mutation rate was calculated as follow:\u003c/p\u003e\u003cp\u003eLet R(x) represent the mutation density at position x within the 2kb region:\u003c/p\u003e\u003cp\u003e \u003cspan class=\"InlineEquation\"\u003e \u003cspan class=\"mathinline\"\u003e\\(\\:R\\left(x\\right)=\\frac{\\text{N}\\left(x\\right)}{\\text{L}.\\text{S}}\\:\\:\\)\u003c/span\u003e \u003c/span\u003e× 10\u003csup\u003e6\u003c/sup\u003e Where:\u003c/p\u003e\u003cul\u003e \u003cli\u003e \u003cp\u003eN(x) = Number of mutations observed at position x across all samples and loops × 2 (both anchors)\u003c/p\u003e \u003c/li\u003e \u003cli\u003e \u003cp\u003eL = Total number of loops × 2\u003c/p\u003e \u003c/li\u003e \u003cli\u003e \u003cp\u003eS = Total number of samples\u003c/p\u003e \u003c/li\u003e \u003cli\u003e \u003cp\u003e10\u003csup\u003e6\u003c/sup\u003e = Scaling factor for mutations per megabase (Mb)\u003c/p\u003e \u003c/li\u003e \u003c/ul\u003e\u003cp\u003eMutation densities for each region type (CTCFbs and surrounding) were calculated for each sample by aggregating the total number of mutations across all motif regions (1,806), dividing by the total genomic length of these regions collectively, and scaling to mutations per megabase.\u003c/p\u003e\u003ch2\u003eExploring Various Genomic Elements:\u003c/h2\u003e\u003cp\u003eWe investigated several genomic elements, including 3' UTR, 5' UTR, promoter regions, long non-coding regions, and CDS regions, to assess (a) whether these elements tend to harbor more mutations compared to their surrounding regions, and (b) how mutational density varies due to differences in the region lengths within each element.\u003c/p\u003e\u003cp\u003eTo conduct these analyses, we developed an in-house script, which is available upon request via GitHub (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/amarinderthind/CTCF_Cancer_study\u003c/span\u003e\u003cspan address=\"https://github.com/amarinderthind/CTCF_Cancer_study\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e ). For the first hypothesis, we defined surrounding regions for each genomic element using different scaling factors (10, 50, 100, 200, 400, 600, and 800), which were applied based on the region length. The scaling factor was applied to calculate the upstream and downstream boundaries of each region, where the surrounding regions' length was determined by multiplying the element length by the scaling factor (Supplementary Fig.\u0026nbsp;1). For each scaling factor, the formula for defining the surrounding region was:\u003c/p\u003e\u003cp\u003eSurrounding Region Length = Region Length × Scaling Factor\u003c/p\u003e\u003cp\u003eFor each region, we computed the number of mutations that overlap with the surrounding regions, and the mutation density was calculated by normalizing the counts to mutations per megabase (mutations/Mb). To account for varying sample sizes and genomic element lengths, we normalized the counts using the total number of CTCFbs motifs regions (903*2 = 1806) and total sample size (72). The scaling was applied as:\u003c/p\u003e\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:\\text{M}\\text{u}\\text{t}\\text{a}\\text{t}\\text{i}\\text{o}\\text{n}\\:\\text{D}\\text{e}\\text{n}\\text{s}\\text{i}\\text{t}\\text{y}\\:\\left(\\text{p}\\text{e}\\text{r}\\:\\text{M}\\text{b}\\right)=\\frac{\\text{}\\text{N}\\text{u}\\text{m}\\text{b}\\text{e}\\text{r}\\:\\text{o}\\text{f}\\:\\text{M}\\text{u}\\text{t}\\text{a}\\text{t}\\text{i}\\text{o}\\text{n}\\text{s}}{\\left(\\text{R}\\text{e}\\text{g}\\text{i}\\text{o}\\text{n}\\:\\text{L}\\text{e}\\text{n}\\text{g}\\text{t}\\text{h}\\right)\\text{x}\\:\\text{n}\\text{u}\\text{m}\\text{b}\\text{e}\\text{r}\\:\\text{o}\\text{f}\\:\\text{s}\\text{a}\\text{m}\\text{p}\\text{l}\\text{e}\\text{s}}\\:\\text{x}\\:\\text{1000,000}\\:$$\u003c/div\u003e\u003c/div\u003e\u003cp\u003eAdditionally, mutations were annotated to identify whether they occurred within the target region (e.g., 3' UTR, promoter, etc.). We then compared mutation densities between the regions and their surrounding elements at each scaling factor. For visualization, mean mutation densities across all samples were calculated and compared for each scaling factor (10x, 50x, 100x, etc.), highlighting differences in mutation patterns across genomic elements and their surrounding regions.\u003c/p\u003e\u003ch2\u003eMutational Substitution Types in CTCFbs motif and cSCC progression:\u003c/h2\u003e\u003cp\u003eTo characterize mutational patterns across samples, we analyzed single nucleotide substitutions across whole genomes. Substitution types were classified into six categories (C \u0026gt; A, C \u0026gt; G, C \u0026gt; T, T \u0026gt; A, T \u0026gt; C, T \u0026gt; G). To compare substitution profiles across sample groups, mean mutation proportions were computed and assessed using Wilcoxon rank-sum tests. To assess differences in mutational profiles during cSCC progression, we compared motif-specific mutations across different disease stages.\u003c/p\u003e\u003cp\u003e \u003cstrong\u003eAssociation between CTCFbs mutations and the expression of harbouring and nearby genes\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eAs shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eG, \u003cb\u003eDifferential gene expression (DGE)\u003c/b\u003e analyses were performed for each loop using data from the mutation's matrix and the RNA-Seq count matrix. For each loop, two comparison groups were created based on the presence or absence of observed mutations at CTCFbs to assess the association of these mutations with the expression of genes within the loop and nearby genes (within 1 kb). Loops with CTCFbs mutated in fewer than 3 samples were excluded from this DGE analysis and 54 shared samples of RNA-Seq/WGS were considered for this analysis. In total, 284 DGE analyses were conducted, one for each loop. A table containing the loops and their associated significantly differentially expressed genes (defined as log2FC \u0026gt;|1| and p-adjust \u0026lt; 0.05) was compiled. DGE is performed using DESeq2 R package (version 1.142.1; Bioconductor 3.18; R 4.3) (\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e), where design formula contains batch-correction as RNA-Seq samples where sequenced in years of time.\u003c/p\u003e\u003cb\u003ePermutation Test\u003c/b\u003e\u003cp\u003eTo assess the statistical significance of the observed association between gene expression changes and mutations in CTCFbs, we performed a permutation test with 1,000 iterations. In each iteration, sample labels were randomly permuted while preserving the original group sizes, and differential gene expression (DGE) analysis was conducted using DESeq2 for genes associated with each of the 284 chromatin loops. Given that the number of genes per loop varied, we quantified, per iteration, the number of loops containing at least one significantly differentially expressed gene (adjusted \u003cem\u003ep\u003c/em\u003e value \u0026lt; 0.05). This yielded a null distribution of loop counts expected under the assumption of no association. The observed number of significant loop-DEG pairs was 17, while the maximum number observed in the permuted (null) distributions was 16. None of the 1,000 permutations produced a value equal to or greater than the observed. Using the standard empirical \u003cem\u003ep\u003c/em\u003e-value formula with a continuity correction:\u003c/p\u003e\u003cdiv id=\"Equb\" class=\"Equation\"\u003e \u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e\n$$\\:\\text{p}=\\frac{\\text{r}+1}{\\text{N}+1}\\:\\text{}=\\frac{0+1}{1000+1}\\text{}=\\frac{1}{1001}\\text{}\\approx\\:0.001$$\u003c/div\u003e \u003c/div\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eWhere:\u003c/p\u003e\u003cul\u003e \u003cli\u003e \u003cp\u003e \u003cspan class=\"InlineEquation\"\u003e \u003cspan class=\"mathinline\"\u003e\\(\\:\\text{p}\\)\u003c/span\u003e \u003c/span\u003e is the empirical \u003cem\u003ep\u003c/em\u003e-value, representing the probability of observing a value as extreme as or more extreme than the actual value under the null hypothesis.\u003c/p\u003e \u003c/li\u003e \u003cli\u003e \u003cp\u003er is the number of permutations in which the test statistic (e.g., number of significant loop-DEG pairs) was greater than or equal to the observed value (\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e).\u003c/p\u003e \u003c/li\u003e \u003cli\u003e \u003cp\u003eN is the total number of permutations performed.\u003c/p\u003e \u003c/li\u003e \u003c/ul\u003e\u003cp\u003ethis result indicates that the likelihood of observing 17 or more significant associations under the null hypothesis is less than 0.1%, supporting the non-random nature of the observed loop-DEG associations.\u003c/p\u003e\u003ch2\u003eAssociations of CTCF binding site mutations and progression of cSCC:\u003c/h2\u003e\u003cp\u003eTo investigate the role of CTCF binding site mutations in the progression of cSCC, we studied different cSCC sub-cohorts, i.e. cSCC primary tumors (non-metastasizing, n = 16; and metastasizing tumor, n = 15) and LNM (n = 41) samples.\u003c/p\u003e\u003cp\u003eTo identify loops with significantly different mutational CTCFbs between groups, a Chi-square test was applied when the expected frequency in each cell was ≥ 5, while Fisher's exact test was used when the expected frequency in any cell was \u0026lt; 5. A confusion matrix was constructed using the mutation data matrix, where '0' indicates the absence of a mutation event and '≥1' indicates the presence of a mutation event. Various comparison of groups was performed as reported in Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e3\u003c/span\u003e. After getting the DE loops, genes list was extracted to perform the RNA-Seq DGE analysis.\u003c/p\u003e\u003cp\u003eTo further assess whether gene expression was associated with the mutation status of the CTCFbs of TAD loop in which each gene resides, we calculated Spearman’s rank correlation coefficient (ρ) between gene expression (treated as a continuous variable) and loop mutation status (binary: mutated vs. non-mutated). Spearman’s ρ is used for testing general associations between a continuous and a binary variable. In addition, we computed the area under the receiver operating characteristic curve (AUC) to evaluate the classification performance of gene expression in distinguishing between mutated and non-mutated loops. These analyses were performed in R, using the \u003cem\u003epROC package\u003c/em\u003e for AUC computation and base functions for correlation analysis.\u003c/p\u003e\u003cdiv class=\"gridtable\"\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e\u003ctable float=\"Yes\" id=\"Tab3\" border=\"1\"\u003e\u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 4\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003e\u003cem\u003eDemographic and clinicopathological features of patients with lymph node metastases (LNM), primary metastatic (PM), and primary non-metastatic (PNM) cutaneous squamous cell carcinoma (cSCC).\u003c/em\u003e n \u003cem\u003erefers to the number of samples (total = 72) derived from 59 patients, including 13 matched sample pairs between LNM and PM.\u003c/em\u003e\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e\u003ccolgroup cols=\"4\"\u003e\u003c/colgroup\u003e\u003cthead\u003e\u003ctr\u003e\u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003eVariable\u003c/p\u003e \u003c/th\u003e\u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003eLNM\u003c/p\u003e \u003cp\u003en = 41\u003c/p\u003e \u003c/th\u003e\u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003ePM\u003c/p\u003e \u003cp\u003en = 16\u003c/p\u003e \u003c/th\u003e\u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003ePNM\u003c/p\u003e \u003cp\u003en = 15\u003c/p\u003e \u003c/th\u003e\u003c/tr\u003e\u003c/thead\u003e\u003ctbody\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eMean age, years (range)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e70.8 (30–92)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e70.6 (51–92)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e75.7 (57–91)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eSex, n (%)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eFemale\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e4 (9.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e4 (\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e3 (\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eMale\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e37 (90.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e12(75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e12 (80)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eSite of primary tumour, n (%)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003echeek\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e3 (18.75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e5 (33.33)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eneck\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003elip\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eeyebrow\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eear\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e3 (\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003etemple\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e2 (12.5)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1 (6.66)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003enose\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1 (6.66)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003epostauricular\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e2 (12.5)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1 (6.66)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003escalp\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e4 (\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e2 (13.33)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eface\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1 (6.66)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003epre-auricular\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1 (6.66)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eSite of metastasis, n (%)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eNeck\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e21 (51.2)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e12 (75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eParotid and neck\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e4 (9.75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eparotid\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e15 (36.57)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e3 (18.75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eperifacial\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e1 (2.43)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eT-stage at surgery, n (%)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e0 or unknown\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e1\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e3 (18.75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e6 (\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e2\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e7 (43.75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e4 (26.67)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e3\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNA\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e3 (18.75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e5 (33.33)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e4\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e3 (18.75)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eN stage at surgery, n (%)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e0\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e1\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e2 (4.88)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e2\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e8 (19.51)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e3\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e18 (43.9)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e4 (\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eunknown\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e13 (31.71)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e10 (62.5)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eHistopathological grading, n (%)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e1\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e2 (4.88)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1 (6.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1 (6.66)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e2\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e5 (12.19)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e2 (12.5)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e10 (66.67)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003e3\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e24 (58.54)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e9 (56.25)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e4 (26.67)\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eunknown\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e10 (24.39)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e4 (\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eLympho-vascular infiltration (LVI), n (%)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eNo\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e16\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e7\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e14\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eYes\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e15\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e8\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eUnknown\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e10\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003ePeri-neural invasion (PNI), n (%)\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u0026nbsp;\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eno\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e20\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e12\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e12\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eyes\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e12\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e3\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e3\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eunknown\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e9\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e1\u003c/p\u003e \u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e\u003c/tr\u003e\u003c/tbody\u003e\u003c/table\u003e\u003c/div\u003e"},{"header":"Abbreviations","content":"\u003cp\u003e\u003cstrong\u003eTAD:\u0026nbsp;\u003c/strong\u003eTopologically Associated Domains\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCSCC:\u003c/strong\u003e Cutaneous Squamous Cell Carcinoma\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCTCFbs:\u003c/strong\u003e CTCF Binding Site\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eLNM:\u003c/strong\u003e Lymph node metastases\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePM:\u003c/strong\u003e Primary metastatic\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePNM:\u003c/strong\u003e Primary non-metastatic\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eWGS\u003c/strong\u003e: Whole Genome Sequencing\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eRNA-Seq\u003c/strong\u003e: RNA Sequencing\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDEG\u003c/strong\u003e: Differentially Expressed Gene\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDGE\u003c/strong\u003e: Differential Gene Expression\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eTMB\u003c/strong\u003e: Tumor Mutational Burden\u003c/p\u003e\n"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate: \u003c/strong\u003eThis study was conducted with approval from the Institutional Human Research Ethics Committee (UOW/ISLHD HREC14/397).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for publication: \u003c/strong\u003e\u0026ldquo;Not applicable\u0026rdquo; \u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAvailability of data and materials: \u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e(a) \u003cu\u003eCode Availability\u003c/u\u003e: \u003c/strong\u003eAll the in-house script used for the secondary analyses are available from GitHub at https://github.com/amarinderthind/CTCF_Cancer_study\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e (b)\u003c/strong\u003e \u003cstrong\u003e\u003cu\u003eData availability\u003c/u\u003e: \u003c/strong\u003eData is available on suitable request to the corresponding author. \u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests: \u003c/strong\u003eThe authors declare that they have no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding: \u003c/strong\u003eThis work was funded by the Illawarra Cancer Carers, Cancer Institute NSW translational program grant 2020/TPG2081, National Health and Medical Research Council Project Grant APP1181179, and Tour de Cure RSP-00244-19/20. \u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026apos; contributions:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eConceptualization: AT (lead), BA (supporting), MR (supporting); Data curation: AT (equal), DS (equal); Formal analysis: AT (leading); Investigation: AT (lead), BA (supporting), MR (supporting), NS (supporting); Methodology: AT (lead), DS (supporting); Visualization: AT (lead), AKP (supporting), Writing - original draft: AT; Writing - editing: AT, MR, BA; Writing \u0026ndash; review: AT, BA, MR, SM, NS, DS; Funding acquisition: BA (equal), MR (equal), RG (equal), JC (equal).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgements:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe wish to acknowledge the National Computational Infrastructure (NCI) of Australia for providing computational resources that contributed to these results.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eThind AS, Ashford B, Strbenac D, Mitchell J, Lee J, Mueller SA, et al. Whole genome analysis reveals the genomic complexity in metastatic cutaneous squamous cell carcinoma. Frontiers in oncology. 2022;12:919118.\u003c/li\u003e\n\u003cli\u003eMueller SA, Gauthier M-EA, Ashford B, Gupta R, Gayevskiy V, Ch\u0026rsquo;ng S, et al. Mutational patterns in metastatic cutaneous squamous cell carcinoma. Journal of investigative dermatology. 2019;139(7):1449-58. e1.\u003c/li\u003e\n\u003cli\u003eBulen BJ, Khazanov NA, Hovelson DH, Lamb LE, Matrana M, Burkard ME, et al. Validation of Immunotherapy Response Score as predictive of pan-solid tumor anti-PD-1/PD-L1 benefit. Cancer research communications. 2023;3(7):1335-49.\u003c/li\u003e\n\u003cli\u003eKim TH, Abdullaev ZK, Smith AD, Ching KA, Loukinov DI, Green RD, et al. Analysis of the vertebrate insulator protein CTCF-binding sites in the human genome. Cell. 2007;128(6):1231-45.\u003c/li\u003e\n\u003cli\u003eLong HS, Greenaway S, Powell G, Mallon A-M, Lindgren CM, Simon MM. Making sense of the linear genome, gene function and TADs. Epigenetics \u0026amp; Chromatin. 2022;15(1):4.\u003c/li\u003e\n\u003cli\u003eDixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485(7398):376-80.\u003c/li\u003e\n\u003cli\u003eYang J, Corces VG. Chromatin insulators: a role in nuclear organization and gene expression. Advances in cancer research. 2011;110:43-76.\u003c/li\u003e\n\u003cli\u003ePoulos RC, Thoms JA, Guan YF, Unnikrishnan A, Pimanda JE, Wong JW. Functional mutations form at CTCF-cohesin binding sites in melanoma due to uneven nucleotide excision repair across the motif. Cell reports. 2016;17(11):2865-72.\u003c/li\u003e\n\u003cli\u003eFang C, Wang Z, Han C, Safgren SL, Helmin KA, Adelman ER, et al. Cancer-specific CTCF binding facilitates oncogenic transcriptional dysregulation. Genome biology. 2020;21:1-30.\u003c/li\u003e\n\u003cli\u003eGupta R, Strbenac D, Satgunaseelan L, Cheung VK-Y, Narayanappa H, Ashford B, et al. Comparing genomic landscapes of oral and cutaneous squamous cell carcinoma of the head and neck: quest for novel diagnostic markers. Modern Pathology. 2023;36(8):100190.\u003c/li\u003e\n\u003cli\u003eDamaschke NA, Gawdzik J, Avilla M, Yang B, Svaren J, Roopra A, et al. CTCF loss mediates unique DNA hypermethylation landscapes in human cancers. Clinical Epigenetics. 2020;12:1-13.\u003c/li\u003e\n\u003cli\u003eWang S, Wang Y, Hao L, Chen B, Zhang J, Li X, et al. BOC targets SMO to regulate the Hedgehog pathway and promote proliferation, migration, and invasion of glioma cells. Brain Research Bulletin. 2024;216:111037.\u003c/li\u003e\n\u003cli\u003eDumitru AV, Țăpoi DA, Halcu G, Munteanu O, Dumitrascu D-I, Ceaușu MC, et al. The polyvalent role of CD30 for cancer diagnosis and treatment. Cells. 2023;12(13):1783.\u003c/li\u003e\n\u003cli\u003eGao Y, Shi H, Zhao H, Yao M, He Y, Jiang M, et al. Single‐cell transcriptomics identify TNFRSF1B as a novel T‐cell exhaustion marker for ovarian cancer. Clinical and Translational Medicine. 2023;13(9):e1416.\u003c/li\u003e\n\u003cli\u003eVan der Weyden C, Pileri S, Feldman A, Whisstock J, Prince H. Understanding CD30 biology and therapeutic targeting: a historical perspective providing insight into future directions. Blood cancer journal. 2017;7(9):e603-e.\u003c/li\u003e\n\u003cli\u003eLing Q, Li F, Zhang X, Mao S, Lin X, Pan J, et al. MAP4K1 functions as a tumor promotor and drug mediator for AML via modulation of DNA damage/repair system and MAPK pathway. EBioMedicine. 2021;69.\u003c/li\u003e\n\u003cli\u003eChan L-K, Ho DW-H, Kam CS, Chiu EY-T, Lo IL-O, Yau DT-W, et al. RSK2-inactivating mutations potentiate MAPK signaling and support cholesterol metabolism in hepatocellular carcinoma. Journal of Hepatology. 2021;74(2):360-71.\u003c/li\u003e\n\u003cli\u003eZheng K, Yao S, Yao W, Li Q, Wang Y, Zhang L, et al. Association between RSK2 and clinical indexes of primary breast cancer: a meta-analysis based on mRNA microarray data. Frontiers in genetics. 2021;12:770134.\u003c/li\u003e\n\u003cli\u003eShi J-l, Fu L, Ang Q, Wang G-j, Zhu J, Wang W-d. Overexpression of ATP1B1 predicts an adverse prognosis in cytogenetically normal acute myeloid leukemia. Oncotarget. 2015;7(3):2585.\u003c/li\u003e\n\u003cli\u003eAcconcia F. Evaluation of the sensitivity of breast cancer cell lines to cardiac glycosides unveils atp1b3 as a possible biomarker for the personalized treatment of er\u0026alpha; expressing breast cancers. International journal of molecular sciences. 2022;23(19):11102.\u003c/li\u003e\n\u003cli\u003eThomas DD, Wink DA. NOS2 as an emergent player in progression of cancer. Mary Ann Liebert, Inc. 140 Huguenot Street, 3rd Floor New Rochelle, NY 10801 USA; 2017. p. 963-5.\u003c/li\u003e\n\u003cli\u003eCoutinho LL, Femino EL, Gonzalez AL, Moffat RL, Heinz WF, Cheng RY, et al. NOS2 and COX-2 Co-expression promotes cancer progression: a potential target for developing agents to prevent or treat highly aggressive breast cancer. International Journal of Molecular Sciences. 2024;25(11):6103.\u003c/li\u003e\n\u003cli\u003eHe W, Wang Q, Xu J, Xu X, Padilla MT, Ren G, et al. Attenuation of TNFSF10/TRAIL-induced apoptosis by an autophagic survival pathway involving TRAF2-and RIPK1/RIP1-mediated MAPK8/JNK activation. Autophagy. 2012;8(12):1811-21.\u003c/li\u003e\n\u003cli\u003eYi JH, Kim SJ, Kim WS. Brentuximab vedotin: clinical updates and practical guidance. Blood research. 2017;52(4):243-53.\u003c/li\u003e\n\u003cli\u003eWang X, Wang Q, Wang H, Cai G, An Y, Liu P, et al. Small protein ERSP encoded by LINC02870 promotes triple negative breast cancer progression via IRE1\u0026alpha;/XBP1s activation. Cell Death \u0026amp; Differentiation. 2025:1-12.\u003c/li\u003e\n\u003cli\u003eGuo M, Zhuang H, Huang J, Shao X, Bai N, Li M, et al. LINC02870 facilitates SNAIL translation to promote hepatocellular carcinoma progression. Molecular and Cellular Biochemistry. 2023;478(9):1899-914.\u003c/li\u003e\n\u003cli\u003eUrbaniak A, Jablonska K, Podhorska-Okolow M, Ugorski M, Dziegiel P. Prolactin-induced protein (PIP)-characterization and role in breast cancer progression. American journal of cancer research. 2018;8(11):2150.\u003c/li\u003e\n\u003cli\u003eZou Y, Cao C, Wang Y, Zhou Y, Yao S, Zhang L, et al. Multi-omics consensus portfolio to refine the classification of lung adenocarcinoma with prognostic stratification, tumor microenvironment, and unique sensitivity to first-line therapies. Translational Lung Cancer Research. 2022;11(11):2243.\u003c/li\u003e\n\u003cli\u003eHan WJ, He P. A novel tumor microenvironment-related gene signature with immune features for prognosis of lung squamous cell carcinoma. Journal of Cancer Research and Clinical Oncology. 2023;149(14):13137-54.\u003c/li\u003e\n\u003cli\u003eDi Iorio P, Beggiato S, Ronci M, Nedel C, Tasca C, Zuccarini M. Unfolding new roles for guanine-based purines and their metabolizing enzymes in cancer and aging disorders. Frontiers in Pharmacology. 2021;12:653549.\u003c/li\u003e\n\u003cli\u003eMason JM, Xu H-P, Rao SK, Leask A, Barcia M, Shan J, et al. Lefty contributes to the remodeling of extracellular matrix by inhibition of connective tissue growth factor and collagen mRNA expression and increased proteolytic activity in a fibrosarcoma model. Journal of Biological Chemistry. 2002;277(1):407-15.\u003c/li\u003e\n\u003cli\u003eGao X, Cai Y, An R. miR-215 promotes epithelial to mesenchymal transition and proliferation by regulating LEFTY2 in endometrial cancer. International journal of molecular medicine. 2018;42(3):1229-36.\u003c/li\u003e\n\u003cli\u003eDo C, Skok JA. Factors that determine cell type\u0026ndash;specific CTCF binding in health and disease. Current Opinion in Genetics \u0026amp; Development. 2024;88:102244.\u003c/li\u003e\n\u003cli\u003eThind AS, Monga I, Thakur PK, Kumari P, Dindhoria K, Krzak M, et al. Demystifying emerging bulk RNA-Seq applications: the application and utility of bioinformatic methodology. Briefings in bioinformatics. 2021;22(6):bbab259.\u003c/li\u003e\n\u003cli\u003eKrueger F. Trim Galore!: A wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files, with extra functionality for RRBS data. Babraham Institute. 2015.\u003c/li\u003e\n\u003cli\u003eDobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15-21.\u003c/li\u003e\n\u003cli\u003eLiao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923-30.\u003c/li\u003e\n\u003cli\u003eRao SS, Huntley MH, Durand NC, Stamenova EK, Bochkov ID, Robinson JT, et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell. 2014;159(7):1665-80.\u003c/li\u003e\n\u003cli\u003eConsortium EP. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489(7414):57.\u003c/li\u003e\n\u003cli\u003eQuinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841-2.\u003c/li\u003e\n\u003cli\u003eDurinck S, Spellman PT, Birney E, Huber W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nature protocols. 2009;4(8):1184-91.\u003c/li\u003e\n\u003cli\u003eLove MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology. 2014;15:1-21.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"CTCF binding sites, Topologically Associated Domains, Cutaneous squamous cell carcinoma, Whole Genome Sequencing, 3D genome architecture","lastPublishedDoi":"10.21203/rs.3.rs-6844715/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-6844715/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003e\u003cstrong\u003eBackground\u003c/strong\u003e\u003cbr\u003e\nCutaneous squamous cell carcinoma (cSCC) is the most common lethal malignancy with metastatic potential. The high mutational burden in cSCC has made it difficult to understand the significance of variants in the noncoding and regulatory genome. This study presents the first investigation\u003cstrong\u003e \u003c/strong\u003eof mutations at CCCTC-binding factor binding sites (CTCFbs) of topologically associating domains (TADs) across defined stages of disease progression —primary tumours that have not metastasized, primary tumours that have metastasized, and lymph node metastasis.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eResults\u003c/strong\u003e\u003cbr\u003e\nBy integrating matched whole-genome sequencing and RNA sequencing data from the same tumors, we reveal that CTCFbs are mutation hotspots (~1,100 mutations/Mb) in cSCC, with mutation densities far exceeding genome-wide averages (170–250 mutations/Mb). This study is also the first to prove genome-wide association of CTCFbs mutations with gene expression changes in cSCC. We report the novel finding\u003cstrong\u003e \u003c/strong\u003ethat CTCFbs that overlap with other regulatory elements such as promoters and untranslated regions exhibit even higher mutational densities than overall CTCFbs. A pattern of mutually exclusive TAD loop CTCFbs mutations was observed between non-metastasizing and metastasizing primary tumors. TAD loops with mutated CTCFbs, are associated with significant transcriptional changes in the genes within those TADs, implicating genes including \u003cem\u003eHHIPL2\u003c/em\u003e, \u003cem\u003eLINC02870\u003c/em\u003e, and \u003cem\u003eGDA\u003c/em\u003e in cSCC progression.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConclusions\u003c/strong\u003e\u003cbr\u003e\nOur findings highlight the functional relevance of noncoding mutations at CTCFbs in cSCC and suggest their potential influence as drivers of tumor progression and metastasis. This integrative genomic analysis with detailed examination of TADs provides a foundation for future studies into 3D genome dysregulation in skin cancers.\u003c/p\u003e","manuscriptTitle":"CTCF Binding Site Mutations: Linking Topologically Associated Domains Dysregulation to Cutaneous Squamous Cell Carcinoma Progression","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-06-24 07:20:44","doi":"10.21203/rs.3.rs-6844715/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"e4b2b1e7-3928-4831-9c75-5656b83fe468","owner":[],"postedDate":"June 24th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2025-06-24T07:20:45+00:00","versionOfRecord":[],"versionCreatedAt":"2025-06-24 07:20:44","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-6844715","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-6844715","identity":"rs-6844715","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}