OR Oncology Reports 1021-335X 1791-2431 D.A. Spandidos 10.3892/or.2022.8408 OR-48-05-08408 Articles Intra-tumor heterogeneity of cancer stem cell-related genes and their potential regulatory microRNAs in metastasizing colorectal carcinoma UrhKristian 1 ZidarNina 1 TomažičAleš 2 3 BoštjančičEmanuela 1 Institute of Pathology, Faculty of Medicine, University of Ljubljana, 1000 Ljubljana, Slovenia Department of Abdominal Surgery, Division of Surgery, University Medical Centre Ljubljana, 1000 Ljubljana, Slovenia Faculty of Medicine, University of Ljubljana, 1000 Ljubljana, Slovenia Correspondence to: Dr Emanuela Boštjančič, Institute of Pathology, Faculty of Medicine, University of Ljubljana, Korytkova 2, 1000 Ljubljana, Slovenia, E-mail: emanuela.bostjancic@mf.uni-lj.si 11 2022 14 09 2022 48 5 193 06042022 13072022 Copyright © 2022, Spandidos Publications 2022

Intra-tumor heterogeneity (ITH) is related to cancer progression, therapy resistance and recurrences, and is one of the challenging fields in cancerogenesis research. Cancer stem cells (CSC) are thought to be crucially involved in the pathogenesis of several cancer types, including colorectal carcinoma (CRC), and associated with ITH. In the present study, the expression gradient of four genes related to CSC (L1TD1, SLITRK6, ST6GALNAC1 and TCEA3) and their potential regulatory microRNAs (miRNAs) were investigated in the central part and invasive front of the primary tumor, as well as in lymph node and liver metastases. In total, 63 formalin-fixed paraffin-embedded biopsy samples of primary tumor (central part, invasive tumor front), as well as lymph node and liver metastases from 19 patients with CRC, were analyzed. The expression of selected genes (L1TD1, SLITRK6, ST6GALNAC1 and TCEA3) and miRNAs (miR-199a-3p, miR-425-5p, miR-1225-3p, miR-1233-3p and miR-1303) was evaluated using reverse transcription-quantitative PCR. Significant differences in expression were identified for all investigated genes in lymph node metastasis, but not in the liver metastases. All investigated miRNAs were significantly differentially expressed in lymph node metastasis, and miR-199a-3p, miR-425-5p and miR-1233-3p in liver metastasis. Furthermore, a negative correlation between the expression of miR-199a-3p and expression of its potential target gene SLITRK6 was confirmed. The present results provide further evidence that expression of CSC-related genes and their potential regulatory miRNAs contribute to ITH in CRC, lymph node and liver metastasis. The SLITRK6 gene and its regulatory miRNA miR-199a-3p are promising for further validation in functional studies to deepen the present understanding of the regulation of CSC-related genes in CRC.

 intra-tumor heterogeneity cancer stem cell-related genes miRNAs metastases colorectal carcinoma Slovenian Research Agency P3-0054 J3-1754 This research was funded by the Slovenian Research Agency (research core funding grant no. P3-0054; project funding grant no. J3-1754; PhD research funding for KU). Introduction

One of the major issues in the field of oncology is tumor heterogeneity. Tumor heterogeneity is described as differences between tumors of the same type in different patients (inter-tumor heterogeneity), as well as between cancer cells in a single tumor of one patient [intra-tumor heterogeneity (ITH)]. These differences may be morphological, physiological and/or genetic and may result in differences in progression, metastasis and response to treatment (1,2).

In sporadic colorectal cancers (CRC), inter-tumor heterogeneity is widely present, resulting in cases without any identified well-known or with numerous different genetic aberrations, suggesting that there are still undiscovered mechanisms (3,4). One of the challenging fields in CRC progression is ITH, which is frequently a source of variability in tumors. ITH may be characterized by temporal (variations within a given tumor over time) and spatial (variations in distinctive regions of a tumor) differences in genetic mutations, epigenetic regulation and expression of coding and non-coding genes, which may be specific to each individual patient. ITH may be clonal and related to different types of genomic instability that may also explain certain differences between the primary tumor and its metastasis. ITH is also closely related to cancer progression, aggressiveness, therapy resistance and recurrences (4).

Tumor progression is thought to rely on a minority of cells in a given tumor recognized as cancer stem cell (CSC)-like cells, which are capable of self-renewal and differentiation. It is thought that CSCs represent the basis for tumor growth and metastatic spread (3), since they have a higher propensity towards invasion, suggesting that they are enriched in all stages of metastasis. CSCs may be found as a sub-population on the invasive tumor front as well as based on genes related to CSCs, supporting this hypothesis (5–7). Certain studies reported the potential involvement of CSCs in the progression of CRC (8,9). Genes associated with CSC features may be promising prognostic and therapeutic markers. It has been previously indicated that CSC-associated molecular profiles may predict tumour regeneration and disease relapse after conventional therapy in patients with CRC (9–14). Several potential anti-CSC targeted drugs have emerged in previous studies, with some of them making their way to the clinic (15). Furthermore, studying microRNAS (miRNAs) regulating selected genes is a promising therapeutic approach (16), with miRNAs having anti- or pro-metastatic effects, while they may also be considered as potential biomarkers for metastasis (17).

Clonal selection of CSCs is considered the main mechanism underlying differences between primary tumors and metastasis, suggesting that ITH is present in metastatic spread. Polyclonal seeding of CSCs and inter-metastatic exchange of cancer clones are also possible and may contribute to ITH (4). In addition, in most CRC cases, the sub-clonal origin of the local lymph node metastases is thought to be different from that of distant metastases (18), where liver metastases is most commonly due to anatomical factors related to the portal circulation (19), resulting in a complex picture of potential expression patterns for CSC-related genes and regulatory miRNAs.

Sampling from the border of the tumor, including the surrounding stroma and the sub-border in comparison to the central part, may provide distinctive information, since different areas of the same tumor may have different patterns of gene expression. Differences in expression patterns between the central part of primary CRC and its invasive tumor front, as well as between primary tumor and lymph node and liver metastasis, have been reported (4). A total of four CSC-related genes LINE1 type transposase domain containing 1 (L1TD1), SLIT and NTRK like family member 6 (SLITRK6), ST6 N-acetylgalactosaminide alpha-2,6-sialyltransferase 1 (ST6GALNAC1) and transcription elongation factor A3 (TCEA3) (20–27) and their potential regulatory miRNAs (miR-199a-3p, miR-425-5p, miR-1225-3p, miR-1233-3p and miR-1303) that were indicated to be involved in the cancerogenesis of CRC according to a previous study by our group (28) were now analyzed for spatial and temporal expression changes to investigate their involvement in the development of metastatic CRC.

 Materials and methods <sec> <title>Tissue samples

All tissue samples were fixed for 24 h in 10% buffered formalin prior to paraffin embedding. After this step, tissues were cut into 3–4 µm slices and stained with haematoxylin and eosin for routine histopathological examination. CRC specimens were histopathologically examined and classified according to the pathologic Tumor-Nodes-Metastasis system (29). For the purposes of the present study, representative paraffin blocks from the years 2006 and 2015–2019 were collected retrospectively from the archives of the Institute of Pathology, Faculty of Medicine, University of Ljubljana (Ljubljana, Slovenia).

In all cases, three representative tissue cores were punched using a 0.6-mm inner diameter needle. Tissue cores were taken from the central part of the tumor, invasive tumor front and from lymph node metastases, liver metastases or both. Sporadic CRC cases with corresponding lymph node/liver metastases were included in the study. Only samples that passed the subsequent RNA quality assessment were included in the study.

 RNA isolation from formalin-fixed paraffin-embedded (FFPE) tissue cores

Isolation of total RNA from three tissue cores per region was performed using a MagMax FFPE DNA/RNA Ultra kit (Applied Biosystems; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol with a modification. Protease digestion was performed overnight at 56°C with shaking for 15 sec at 300 rpm every 4 min on an Eppendorf ThermoMixer®C (Eppendorf SE). A NanoDrop®−1000 Spectrophotometer (Thermo Fisher Scientific, Inc.) was used to determine the concentration and assess the quality of the isolates at the wavelengths of 260, 280 and 230 nm.

 RNA quality assessment

Reverse transcription (RT) of RNU6B, a housekeeping small nuclear RNA gene, was used as the quality control, followed by amplification using quantitative real-time PCR (qPCR) and TaqMan methodology (Thermo Fisher Scientific, Inc.). All of the samples included in the study had passed this quality control step. Positive and negative amplification of RNU6B was in positive correlation with positive and negative amplification of GAPDH (100 bp) initially used as quality control in previous studies (data not shown) (30,31). In addition, TaqMan primers and probes that amplify and detect PCR products <100 bp in length were chosen, as indicated in Table I.

 Efficiency testing

A pre-designed mixture of probes and primers specific for miRNAs or target genes (mRNAs) was used. Three pools of RNA samples were created, obtained from primary colorectal tumors, lymph node metastases and liver metastases prior to qPCR. The obtained cDNA of miRNAs and pre-amplified cDNA of mRNAs was diluted in four steps, ranging from 5-point dilution to 625-point dilution, and the probes were tested for qPCR efficiency. The qPCR efficiency reactions were performed on a RotorGene Q (Qiagen GmbH) in triplicate. The efficiency was calculated as follows: E=10(−1/slope)-1 (32).

 RT of miRNAs

Looped primers for specific RT of miRNAs and a MicroRNA TaqMan RT kit (Applied Biosystems; Thermo Fisher Scientific, Inc.) were used as per the manufacturer's protocol for RT. RNU6B and miR-1247b were used as reference genes (RGs). miRNAs were tested relative to the geometric mean of the expression of RNU6B and miR-1247b (Table I). The RT reaction mix (10 µl) was prepared with 10 ng of total RNA sample, 1.0 µl of MultiScribe Reverse Transcriptase (50 U/µl), 1.0 µl of RT Buffer (10X), 0.1 µl of dNTP (100 mM), 0.19 µl RNAase inhibitor (20 U/µl) and 2.0 µl of RT primer (5X). The reaction conditions were: 16°C for 30 min, 42°C for 30 min and 85°C for 5 min.

 qPCR of miRNAs

qPCR for miRNAs was performed in a PCR mixture (10 µl) containing 5.0 µl TaqMan 2X FastStart Essential DNA Probe Master (Roche Diagnostics), 0.5 µl TaqMan assay 20X and 4.5 µl RT products diluted 100-fold. qPCR was performed on a RotorGene Q (Qiagen GmbH) in duplicate, as follows: Initial denaturation at 95°C for 10 min and 40 cycles of 15 sec at 95°C (denaturation) and 60 sec at 60°C (primer annealing and elongation). The signal was collected at the endpoint of every cycle.

 RT of mRNAs

The mRNAs summarized in Table I were analyzed relative to the geometric mean of the RGs IPO8 and B2M. They were reverse transcribed using a OneTaq RT-PCR Kit (New England BioLabs, Inc.) using random primers according to the manufacturer's instructions. RT reactions were performed with 3.0 µl (60 ng) of total RNA and 1.0 µl of Random Primer Mix incubated at 70°C for 5 min. The 10 µl RT mixture included 5.0 µl of M-MuLV Reaction Mix, 1.0 µl of M-MuLV reverse transcriptase and 4.0 µl of reaction mix after random priming. The reaction conditions were as follows: 25°C for 5 min, 42°C for 60 min and 80°C for 4 min.

 Pre-amplification and qPCR of mRNAs

Following RT, pre-amplification was performed using a TaqMan PreAmp Master Mix (Applied Biosystems; Thermo Fisher Scientific, Inc.) in a 10-µl reaction according to the manufacturer's protocol. The obtained PreAmp reaction was diluted 5-fold in all cases, except when investigating lymph-node metastases, where it was diluted 25-fold. For the qPCR, 4.5 µl of the diluted sample was used in a 10-µl reaction volume with 5.0 µl of 2X FastStart Essential DNA Probe Master Mix (Roche Diagnostics) and 0.5 µl of TaqMan 20X probe. The thermal cycling conditions were as follows: 50°C for 2 min, initial denaturation at 95°C for 10 min and 40 cycles of denaturation at 95°C for 15 sec and annealing at 60°C for 1 min. All qPCR analyses were performed on a Rotor Gene Q (Qiagen GmbH) in duplicate. The signal was collected at the endpoint of each cycle.

 Statistical analysis

The results were presented as relative gene expression. All quantification cycle values (Cqs) were corrected for PCR efficiencies and the expression of the gene of interest was calculated relative to a geometric mean of RGs, named ΔCq using the ∆∆Cq method (32). In CRC samples, expression differences in mRNAs and miRNAs were compared between the central part and invasive front, lymph node metastases or liver metastases, using ΔCq and the Wilcoxon Signed-Rank test. For all of the investigated correlations/associations, Spearman rank-order correlation was used. An additional Bonferroni correction was performed for the investigated comparisons after the Wilcoxon Signed Rank test and comparisons that failed to pass the adjusted α-value were mentioned accordingly. Statistical analysis of data was performed using SPSS version 24 (IBM Corporation). Differences were considered significant at P≤0.05.

 Results <sec> <title>Patients and tissue samples

A total of 19 patients with CRC were included in the study, namely seven patients with lymph node but not liver metastasis (N+ M0 group; mean age, 76.0±13.5 years; age range, 54–91 years; males/females, 6:1), three patients without lymph node but with liver metastases (N0 M+ group; mean age, 72.0±6.1 years; age range, 68–79 years; males/females, 2:1) and nine patients with both lymph node and liver metastases (N+ M+ group; mean age, 69.3±16.5 years; age range, 31–88 years; male/female, 5:4). Demographic data and clinicopathological features for each patient are also available in previously published work (33). In total, 63 tissue samples from 19 patients with CRC, with lymph node metastases and/or liver metastases, were analyzed. The invasive tumor front, lymph node and liver metastases were compared to the corresponding central part of the primary tumor.

 Expression of mRNAs

When comparing the invasive tumor front to the central part of the primary tumor, ST6GALNAC1 and TCEA3 exhibited a statistically significant difference in expression prior to Bonferroni correction; downregulation in this case. Regarding comparisons between the central part of the primary tumor and lymph node metastases, all genes investigated were downregulated with statistical significance. No statistically significant differences were observed when comparing liver metastases to the central part of the primary tumor for the investigated genes (Fig. 1).

 Expression of miRNAs

When comparing the invasive tumor front to the central part of the primary tumor, miR-199a-3p exhibited a statistically significant upregulation only prior to Bonferroni correction. Regarding comparisons between the central part of the primary tumor and lymph node metastases, all miRNAs investigated were upregulated, with statistical significance. In liver metastases, significantly upregulated expression was observed when comparing to the central part of the primary tumor for miR-199a-3p, miR-425-5p and miR-1233-3p, as indicated in Fig. 2.

A heatmap of the expression changes of all investigated genes and miRNAs in each patient is additionally available in Table SI. Fold changes and P-values for all investigated comparisons between groups on gene/miRNA expression are available in Table SII for gene expression and Table SIII for miRNA expression.

 Correlations between the investigated mRNAs and miRNAs

Spearman coefficients of the correlation revealed negative correlations of a weak or moderate nature for all significant correlations between the investigated genes and miRNAs. miR-1303 did not correlate significantly with any of the investigated genes. TCEA3 did not correlate significantly with any of the investigated miRNAs. L1TD1 exhibited a significant, negative moderate correlation with miR-199a-3p, miR-425-5p, miR-1225-3p and miR-1233-3p. SLITRK6 had a weak negative significant correlation with its proposed regulatory miRNA miR-199a-3p. In addition, ST6GALNAC1 had weak negative significant correlations with miR-199a-3p, miR-425-5p, miR-1225-3p and miR-1233-3p. ∆Cq comparisons between the potential regulatory miRNAs miR-199a-3p and miR-425-5p and target gene SLITRK6 are presented in Fig. 3. Additional expression comparisons for L1TD1, TCEA3 and their potential regulatory miRNAs are available in Fig. S1. Results including P-values are summarized in Table II.

As presented in Table III, significant positive Spearman correlation coefficients were obtained between all of the investigated miRNAs. The correlations were either weak (e.g. miR-425-5p and miR-1303, rs=0.2-0.39), moderate (e.g. miR-199a-3p and miR-1303, rs=0.40-0.59), strong (e.g. miR-199a-3p and miR-1225-3p, rs=0.6-0.79) or very strong (e.g. miR-199a-3p and miR-425-5p, rs>0.8).

 Discussion

In the present study, the expression of four genes related to CSC and CSC-like properties were validated in CRC tissue samples, obtained from the central part of primary tumors, invasive tumor front, lymph node and liver metastases. The investigated genes (L1TD1, SLITRK6, ST6GALNAC1, TCEA3) were previously identified as differentially expressed between normal mucosa, adenomas and CRC using bioinformatics analysis of publicly available microarray data (34), and validated to be involved in CRC carcinogenesis together with their potential regulatory miRNAs (28).

Regarding the expression of the investigated mRNAs, two different patterns were observed. The first pattern was observed in lymph node metastases compared to the central part of the primary tumor, where all investigated genes were downregulated. The second one was observed in liver metastases, where none of the investigated genes was differentially expressed when compared to the central part of the primary tumor. By contrast, their regulatory miRNAs exhibited variable expression, except in the case of lymph node metastases when compared to the central part of the tumor, where all of the investigated miRNAs exhibited the opposite expression trend to their target mRNAs.

In lymph node metastases, when compared to the central part of the primary tumor, all investigated genes were downregulated. Bioinformatics analysis indicated that a higher expression of L1TD1 in CRC was associated with longer disease-free survival (35), confirming the negative trend of expression of L1TD1 in relation to invasiveness observed in the present study. Small inhibitory RNA-mediated silencing of ST6GALNAC1 in gastric cancer cells was previously reported to lead to reduced growth, migration and invasion of cancer cells (36), whereas its overexpression enhanced their metastatic ability (37). High levels of ST6GALNAC1 were observed in ovarian CSCs and silencing of ST6GALNAC1 was indicated to reduce cell proliferation, migration, invasion, self-renewal ability and tumorigenicity (38). For TCEA3, a previous study by our group reported significantly different expression between CRC without and with lymph node metastases, which suggested a role in lymph node metastasis development (28). TCEA3 has also been associated with gastric cancer, in which high expression has been associated with better prognosis, lower proliferation of carcinoma cells and induction of apoptosis (39). The results thus suggest that higher expression of these genes is necessary for cancer progression to lymph nodes. To the best of our knowledge, the present study was the first to report that the SLITRK6 gene was downregulated in lymph node metastases of CRC, indicating its possible role in metastasizing CRC, which may be further investigated in future studies.

In liver metastases, none of the investigated genes exhibited differential expression when compared to the central part. This observation is not surprising, since in a previous study, gene expression profiling using microarrays clearly distinguished normal colon mucosa and normal liver from primary CRC and liver metastases, respectively; the authors observed moderate variations in expression of most differentially expressed genes, or their dysregulation limited to one individual (40). The observation, that primary CRC and liver metastases have a similar expression pattern was further supported by high-throughput transcriptome sequencing (41). This finding suggests that expression changes consistently occur during CRC development, but only a small number of them may be associated with metastatic progression to the liver (40). In liver metastases, only two different subpopulations of CSCs were identified based on the expression of OCT4, one in the peritumoral stroma and the other in tumor nests (42). The same group reported that the only two populations of CSCs in primary tumor may also be stratified by OCT4 expression (43). The present expression analysis supports this observation of a small number of CSC subclones in primary tumor and liver metastases.

Based on 213 archival biopsy samples investigating genetic changes between primary tumor, lymph node and distant metastases from 17 patients, it has been recently indicated that in 65% of cases of lymph node and liver metastases, they arise from independent sub-clones, whereas in 35% of cases, they share a common origin (18). The same group further confirmed that lymph node and distant metastases develop through different evolutionary mechanisms, with a higher inter-lesion heterogeneity of lymph node metastases (44). Furthermore, the KRAS mutation status in CRC confirmed a much lower level of concordance when comparing primary tumor and matched lymph node metastases; however, a high concordance rate was observed between primary and matched distant (e.g., liver) metastases (3). The observations of the present study further support the hypothesis that not only on the chromosomal, genetic and epigenetic levels, but also on the expression level, there is a higher ITH between primary tumor and lymph node metastases than between primary tumor and liver metastases, further confirming different developmental pathways of lymph node and liver metastases.

In lymph node as well as in liver metastases, all investigated miRNAs were up-regulated when compared to the center of the primary tumor, except miR-1303, which was downregulated in liver metastases. miR-199a-3p was indicated to target stemness and mitogenic-related pathways to suppress the expansion and tumorigenic capabilities of prostate CSCs in vitro (45). Its expression was observed to significantly differ between metastatic and low metastatic groups of patients with uveal melanoma (46). It was not able to promote tumor cell proliferation in melanoma; however, it may regulate metastatic invasion of melanoma, angiogenesis and endothelial cell recruitment (47). An analysis of clinical and pathologic data revealed that a higher miR-199a-3p expression contributed to more advanced lymphatic invasion and lymph node metastases, as well as liver metastases, in CRC (48). miR-425-5p facilitates epithelial-mesenchymal transition and extracellular matrix degradation and promotes hepatocellular carcinoma cell metastasis (49). Of note, it has been indicated that miR-1225-5p suppresses gastric cancer invasion and metastases (50) and that it inhibits apoptosis of pancreatic cancer cells (51). Furthermore, overexpression of miR-1233-3p promoted the migration and invasion of human breast cancer cells (52). Downregulation of miR-1303 inhibited the proliferation, migration and invasion of gastric cancer cells (53) and it also suppressed the proliferation, migration and invasion of prostate cancer cells (54). In addition, miR-1303 was one of the five miRNAs, where the expression signature was an independent predictor of poor metastasis-free survival of patients with breast cancer (55). It was also upregulated in a bioinformatics analysis of microarray data between primary CRC tumors and liver metastases (56). In non-small cell lung cancer, high expression of miR-1303 was associated with TNM stage, lymph node metastasis and shorter survival time, and its overexpression in H1299 and A549 cells promoted cell proliferation, migration and invasion (57). These data suggest that all five genes and all investigated miRNAs are involved in the formation of CRC metastases in lymph nodes, but only three out of five miRNAs are involved in the formation of liver metastases of CRC.

An advantage of the present study is the use of the punching technique, which enabled us to obtain tissue from the locations of interest, determined by microscopic analysis by a pathologist. A limitation of the present study is the relatively small sample size, as cases that fit the inclusion criteria are not very common. Another limitation is that only samples that successfully passed the initial quality control and samples with stable expression of the RGs were selected for further analysis, thus limiting the number of included samples. As another limitation, patients were sub-grouped as patients with only lymph node metastases (N+ M0), patients with only liver metastases (N0 M+) and patients with both lymph node and liver metastases (N+ M+), resulting in a small number of patients with only liver metastases. However, such cases are rare and it is difficult to collect sufficient sample numbers; therefore, these results should be interpreted with caution. Furthermore, there was an unequal male/female ratio; however, there are no public data available that investigated genes and miRNAs are differentially regulated in different genders in CRC. Another limitation of the present study is that the results were not validated further through a functional study using a CRC cell line and experimental animals. Finally, the study lacks direct validation of a regulatory role of the investigated miRNAs on the analyzed genes; however, inverse co-expression analysis may provide indirect support for the predicted regulation.

In conclusion, the present study analyzed the expression patterns of the L1TD1, SLITRK6, ST6GALNAC1 and TCEA3 genes and their potential regulatory miRNAs in tissues obtained using a punching technique to enable the validation of ITH of these genes in CRC. The central part of the primary tumor, the invasive tumor front, as well as lymph node and liver metastases were compared. Of note, all genes and all miRNAs were differentially expressed in lymph node metastases. However, none of the investigated genes were differentially expressed in liver metastases, whereas the majority of miRNAs were. The present results thus indicate a role of all of the investigated genes in the development of lymph node metastases, but for none of them in the development of liver metastases. As CSCs are involved in treatment resistance and disease recurrence, analysis of CSC markers has prognostic and therapeutic potential (9–14). Their potential regulatory miRNAs are easily delivered in vivo due to their small size. Synthetic miRNAs may therefore be administered systematically and may thus serve as therapeutic agents in the future (16,17). Future perspectives include validation of obtained results through a functional study both in vivo and in vitro.

 Supplementary Material Supporting Data Supporting Data Supporting Data Acknowledgements

Not applicable.

 Availability of data and materials

The ∆∆Cq data that support the findings of the present study are available in Table SI. The full ∆Cq data are available from the corresponding author upon reasonable request. Demographic data and clinicopathological features for each patient are available upon reasonable request from the authors or in previously published work by our group (33).

 Authors' contributions

Conceptualization of article: KU, EB and NZ. Methodology: KU and EB. Data acquisition: KU and AT. KU and EB checked and approved the authenticity of the raw data. Formal analysis: KU. Original draft preparation, KU and EB. Review and editing: EB, NZ and AT. Visualization: KU. Supervision: EB, AT and NZ. All authors have read and agreed to the published version of the manuscript.

 Ethics approval and consent to participate

The present study was performed according to the tenets of the Declaration of Helsinki. The study is retrospective, observational, performed on tissue samples that were obtained during routine diagnostic/therapeutic procedures, consisting of either excision or resection. Sufficient tissue was available for routine analysis and research; furthermore, tissue is still available for any additional routine analysis in the future. Prior to excision, or resection, informed consent was obtained for the routine surgical procedure. The National Medical Ethics Committee (Ljubljana, Slovenia) approved the study and any further requirement for informed consent for participation in scientific studies was waived (approval no. 0120-54/2020/7).

 Patient consent for publication

Further need for consent for participation in scientific studies was waived by the National Medical Ethics Committee (approval number 0120-54/2020/7) due to the nature of the study.

 Competing interests

The authors declare that they have no competing interests.

 References MarusykAPolyakKTumor heterogeneity: Causes and consequencesBiochim Biophys Acta18051051172010 HanahanDWeinbergRAHallmarks of cancer: The next generationCell144646674201110.1016/j.cell.2011.02.013 BlankARobertsDEIIDawsonHZlobecILugliATumor heterogeneity in primary colorectal cancer and corresponding metastasesDoes the apple fall far from the tree? Front Med (Lausanne)52342018 StantaGBoninSOverview on clinical relevance of intra-tumor heterogeneityFront Med (Lausanne)585201810.3389/fmed.2018.0008529682505 PragerBCXieQBaoSRichJNCancer stem cells: The architects of the tumor ecosystemCell Stem Cell244153201910.1016/j.stem.2018.12.009 GrilletFBayetEVilleronceOZappiaLLagerqvistELLunkeSCharafe-JauffretEPhamKMolckCRollandNCirculating tumour cells from patients with colorectal cancer have cancer stem cell hallmarks in ex vivo cultureGut6618021810201710.1136/gutjnl-2016-311447 ToloudiMApostolouPChatziioannouMPapasotiriouICorrelation between cancer stem cells and circulating tumor cells and their valueCase Rep Oncol44454201110.1159/000324403 SaikiYIshimaruSMimoriKTakatsunoYNagaharaMIshiiHYamadaKMoriMComprehensive analysis of the clinical significance of inducing pluripotent stemness-related gene expression in colorectal cancer cellsAnn Surg Oncol1626382644200910.1245/s10434-009-0567-5 DyllaSJBevigliaLParkIKChartierCRavalJNganLPickellKAguilarJLazeticSSmith-BerdanSColorectal cancer stem cells are enriched in xenogeneic tumors following chemotherapyPLoS One3e2428200810.1371/journal.pone.000242818560594 GiampieriRScartozziMLoretelliCPivaFMandolesiALezocheGDel PreteMBittoniAFaloppiLBianconiMCancer stem cell gene profile as predictor of relapse in high risk stage II and stage III, radically resected colon cancer patientsPLoS One8e72843201310.1371/journal.pone.007284324023782 Merlos-SuárezABarrigaFMJungPIglesiasMCéspedesMVRossellDSevillanoMHernando-MomblonaXda Silva-DizVMuñozPThe intestinal stem cell signature identifies colorectal cancer stem cells and predicts disease relapseCell Stem Cell8511524201110.1016/j.stem.2011.02.020 ColakSZimberlinCDFesslerEHogdalLPrasetyantiPRGrandelaCMLetaiAMedemaJPDecreased mitochondrial priming determines chemoresistance of colon cancer stem cellsCell Death Differ2111701177201410.1038/cdd.2014.3724682005 LombardoYScopellitiACammareriPTodaroMIovinoFRicci-VitianiLGulottaGDieliFde MariaRStassiGBone morphogenetic protein 4 induces differentiation of colorectal cancer stem cells and increases their response to chemotherapy in miceGastroenterology140297309201110.1053/j.gastro.2010.10.00520951698 LottiFJarrarAMPaiRKHitomiMLathiaJMaceAGanttGAJrSukhdeoKDeVecchioJVasanjiAChemotherapy activates cancer-associated fibroblasts to maintain colorectal cancer-initiating cells by IL-17AJ Exp Med21028512872201310.1084/jem.2013119524323355 ZeunerATodaroMStassiGDe MariaRColorectal cancer stem cells: From the crypt to the clinicCell Stem Cell15692705201410.1016/j.stem.2014.11.012 WangVWuWMicroRNA-based therapeutics for cancerBiodrugs231523200910.2165/00063030-200923010-0000219344188 WangDLiuJHuoTTianYZhaoLThe role of microRNAs in colorectal liver metastasis: Important participants and potential clinical significancesTumour Biol391010428317709640201710.1177/1010428317709640 NaxerovaKReiterJGBrachtelELennerzJKvan de WeteringMRowanACaiTCleversHSwantonCNowakMAOrigins of lymphatic and distant metastases in human colorectal cancerScience3575560201710.1126/science.aai851528684519 Gonzalez-VillarrealCAQuiroz-ReyesAGIslasJFGarza-TreviñoENColorectal cancer stem cells in the progression to liver metastasisFront Oncol101511202010.3389/fonc.2020.01511 WongRCIbrahimAFongHThompsonNLockLFDonovanPJL1TD1 is a marker for undifferentiated human embryonic stem cellsPLoS One6e19355201110.1371/journal.pone.001935521559406 NärväERahkonenNEmaniMRLundRPursiheimoJPNästiJAutioRRasoolODenessioukKLähdesmäkiHRNA-binding protein L1TD1 interacts with LIN28 via RNA and is required for human embryonic stem cell self-renewal and cancer cell proliferationStem Cells30452460201210.1002/stem.1013 EmaniMRNärväEStubbAChakrobortyDViitalaMRokkaARahkonenNMoulderRDenessioukKTrokovicRThe L1TD1 protein interactome reveals the importance of post-transcriptional regulation in human pluripotencyStem Cell Reports4519528201510.1016/j.stemcr.2015.01.01425702638 SandbergCJVik-MoEOBehnanJHelsethELangmoenIATranscriptional profiling of adult neural stem-like cells from the human brainPLoS One9e114739201410.1371/journal.pone.011473925514637 FarahaniEPatraHKJangamreddyJRRashediIKawalecMRao ParitiRKBatakisPWiechecECell adhesion molecules and their relation to (cancer) cell stemnessCarcinogenesis35747759201410.1093/carcin/bgu04524531939 AmbrizXde LanerollePAmbrosioJRThe mechanobiology of the actin cytoskeleton in stem cells during differentiation and interaction with biomaterialsStem Cells Int20182891957201810.1155/2018/2891957 OgawaTHirohashiYMuraiANishidateTOkitaKWangLIkeharaYSatoyoshiTUsuiAKuboTST6GALNAC1 plays important roles in enhancing cancer stem phenotypes of colorectal cancer via the Akt pathwayOncotarget8112550112564201710.18632/oncotarget.22545 ParkKSChaYKimCHAhnHJKimDKoSKimKHChangMYKoJHNohYSTranscription elongation factor Tcea3 regulates the pluripotent differentiation potential of mouse embryonic stem cells via the lefty1-nodal-smad2 pathwayStem Cells31282292201310.1002/stem.128423169579 UrhKŽlajpahMZidarNBoštjančičEIdentification and validation of new cancer stem cell-related genes and their regulatory microRNAs in colorectal cancerogenesisBiomedicines9179202110.3390/biomedicines902017933670246 BrierleyJDGospodarowiczMKWittekindCTNM Classification of Malignant Tumours8th editionWiley-BlackwellOxford, UK2017 RankovićBZidarNŽlajpahMBoštjančičEEpithelial-mesenchymal transition-related microRNAs and their target genes in colorectal cancerogenesisJ Clin Med81603201910.3390/jcm8101603 ŽlajpahMHauptmanNBoštjančičEZidarNDifferential expression of extracellular matrix-related genes DCN, EPHA4, FN1, SPARC, SPON2 and SPP1 in colorectal carcinogenesisOncol Rep42153915482019 LathamGJNormalization of microRNA quantitative RT-PCR data in reduced scale experimental designsMethods Mol Biol6671931201010.1007/978-1-60761-811-9_2 PavličAUrhKŠtajerKBoštjančičEZidarNEpithelial-mesenchymal transition in colorectal carcinoma: Comparison between primary tumor, lymph node and liver metastasesFront Oncol11662806202110.3389/fonc.2021.662806 HauptmanNBoštjančičEŽlajpahMRankovićBZidarNBioinformatics analysis reveals most prominent gene candidates to distinguish colorectal adenoma from adenocarcinomaBiomed Res Int20189416515201810.1155/2018/9416515 ChakrobortyDEmaniMRKlénRBöckelmanCHagströmJHaglundCRistimäkiALahesmaaREloLLL1TD1-a prognostic marker for colon cancerBMC Cancer19727201910.1186/s12885-019-5952-231337362 TamuraFSatoYHirakawaMYoshidaMOnoMOsugaTOkagawaYUemuraNAriharaYMuraseKRNAi-mediated gene silencing of ST6GalNAc I suppresses the metastatic potential in gastric cancer cellsGastric Cancer198597201610.1007/s10120-014-0454-z25532910 OzakiHMatsuzakiHAndoHKajiHNakanishiHIkeharaYNarimatsuHEnhancement of metastatic ability by ectopic expression of ST6GalNAcI on a gastric cancer cell line in a mouse modelClin Exp Metastasis29229238201210.1007/s10585-011-9445-122228572 WangWYCaoYXZhouXWeiBZhanLSunSYStimulative role of ST6GALNAC1 in proliferation, migration and invasion of ovarian cancer stem cells via the Akt signaling pathwayCancer Cell Int1986201910.1186/s12935-019-0780-7 LiJJinYPanSChenYWangKLinCJinSWuJTCEA3 attenuates gastric cancer growth by apoptosis inductionMed Sci Monit2132413246201510.12659/MSM.895860 KoehlerABatailleFSchmidCRuemmelePWaldeckABlaszykHHartmannAHofstaedterFDietmaierWGene expression profiling of colorectal cancer and metastases divides tumours according to their clinicopathological stageJ Pathol2046574200410.1002/path.1606 LeeJRKwonCHChoiYParkHJKimHSJoHJOhNPark doYTranscriptome analysis of paired primary colorectal carcinoma and liver metastases reveals fusion transcripts and similar gene expression profiles in primary carcinoma and liver metastasesBMC Cancer16539201610.1186/s12885-016-2596-327461012 HumphriesHNWickremesekeraSKMarshRWBraschHDMehrotraSTanSTItinteangTCharacterization of cancer stem cells in colon adenocarcinoma metastasis to the liverFront Surg476201810.3389/fsurg.2017.0007629404335 MunroMJWickremesekeraSKPengLMarshRWItinteangTTanSTCancer stem cell subpopulations in primary colon adenocarcinomaPLoS One14e0221963201910.1371/journal.pone.022196331491003 ReiterJGHungWTLeeIHNagpalSGiuntaPDegnerSLiuGWassenaarECEJeckWRTaylorMSLymph node metastases develop through a wider evolutionary bottleneck than distant metastasesNat Genet52692700202010.1038/s41588-020-0633-2 LiuRLiuCZhangDLiuBChenXRycajKJeterCCalhoun-DavisTLiYYangTmiR-199a-3p targets stemness-related and mitogenic signaling pathways to suppress the expansion and tumorigenic capabilities of prostate cancer stem cellsOncotarget75662856642201610.18632/oncotarget.10652 WorleyLALongMDOnkenMDHarbourJWMicro-RNAs associated with metastasis in uveal melanoma identified by multiplexed microarray profilingMelanoma Res18184190200810.1097/CMR.0b013e3282feeac618477892 PenchevaNTranHBussCHuhDDrobnjakMBusamKTavazoieSFConvergent multi-miRNA targeting of ApoE drives LRP1/LRP8-dependent melanoma metastasis and angiogenesisCell15110681082201210.1016/j.cell.2012.10.028 WanDHeSXieBXuGGuWShenCHuYWangXZhiQWangLAberrant expression of miR-199a-3p and its clinical significance in colorectal cancersMed Oncol30378201310.1007/s12032-012-0378-6 FangFSongTZhangTCuiYZhangGXiongQMiR-425-5p promotes invasion and metastasis of hepatocellular carcinoma cells through SCAI-mediated dysregulation of multiple signaling pathwaysOncotarget83174531757201710.18632/oncotarget.15958 ZhengHZhangFLinXHuangCZhangYLiYLinJChenWLinXMicroRNA-1225-5p inhibits proliferation and metastasis of gastric carcinoma through repressing insulin receptor substrate-1 and activation of β-catenin signalingOncotarget746474663201610.18632/oncotarget.6615 ZhongRLiSFangKYangLWangLmicroRNA-1225 inhibit apoptosis of pancreatic cancer cells via targeting JAK1Cell Cycle189901000201910.1080/15384101.2019.160812730990343 LuMWuYZengBSunJLiYLuoJWangLYiZLiHRenGCircEHMT1 inhibits metastatic potential of breast cancer cells by modulating miR-1233-3p/KLF4/MMP2 axisBiochem Biophys Res Commun526306313202010.1016/j.bbrc.2020.03.08432209259 ZhangSJFengJFWangLGuoWDuYWMingLZhaoGQmiR-1303 targets claudin-18 gene to modulate proliferation and invasion of gastric cancer cellsDig Dis Sci5917541763201410.1007/s10620-014-3107-5 LiuBZhouWJiangHXiangZWangLmiR-1303 promotes the proliferation, migration and invasion of prostate cancer cells through regulating the Wnt/β-catenin pathway by targeting DKK3Exp Ther Med1847474757201931772644 LereboursFCizeron-ClairacGSusiniAVacherSMouret-FourmeEBelichardCBrainEAlberiniJLSpyratosFLidereauRBiecheImiRNA expression profiling of inflammatory breast cancer identifies a 5-miRNA signature predictive of breast tumor aggressivenessInt J Cancer13316141623201310.1002/ijc.2817123526361 TorresSGarcia-PalmeroIBartoloméRAFernandez-AceñeroMJMolinaECalviñoESeguraMFCasalJICombined miRNA profiling and proteomics demonstrates that different miRNAs target a common set of proteins to promote colorectal cancer metastasisJ Pathol2423951201710.1002/path.4874 ChenJJiangTYuBLiTZhaoPYuanLQiJUpregulation of microRNA-1303 is a potential prognostic marker of non-small cell lung cancerCancer Biomark28439446202010.3233/CBM-20146132508320

∆Cq data obtained by normalization to the geometric mean of B2M and IPO8 reference genes for the investigated mRNAs between central parts of the primary tumor and the invasive tumor front, corresponding lymph node metastases and/or liver metastases. The statistical significances pertain to comparisons of lymph node metastases or liver metastases with corresponding central parts of the primary tumor. The whiskers of the boxplots represent the values outside of the upper and lower quartiles, up to the maximum and minimum values, excluding outliers, of ∆Cq data in patients for each gene (n=19). **P≤0.01; ***P≤0.001; ◦, outlier; Cq, quantification cycle.

∆Cq data obtained by normalization to the geometric mean of RNU6B and miR-1274b reference genes for the investigated miRNAs between central parts of the primary tumor and the invasive tumor front, corresponding lymph node metastases and/or liver metastases. The statistical significances pertain to comparisons of lymph node metastases or liver metastases with corresponding central parts of the primary tumor. The whiskers of the boxplots represent the values outside of the upper and lower quartiles, up to the maximum and minimum values, excluding outliers, of ∆Cq data in patients for each miRNA (n=19). **P≤0.01; ***P≤0.001; ◦, outlier; Cq, quantification cycle; miRNA, microRNA.

Correlation between the expression of (A) SLITRK6 and miR-199a-3p and (B) SLITRK6 and miR-425-5p for each of the investigated patients (n=19). Cq, quantification cycle; miRNA, microRNA.

Probes used for miRNA and mRNA quantification using reverse-transcription quantitative PCR.

Gene/miRNA name Probe ID/miRNA ID Sequence (probe sequence or mature miRNA sequence in 5′-3′ direction)
B2M Hs99999907_m1 GTTAAGTGGGATCGAGACATGTAAG
IPO8 Hs00183533_m1 GGGGAATTGATCAGTGCATTCCACT
L1TD1 Hs00219459_m1 TTTTTCGCCAGGCACCAAGGCACAG
SLITRK6 Hs00536106_s1 TTTCCATGGACTGGAAAACCTGGAA
ST6GALNAC1 Hs01027885_m1 AGGAGGCCTTCAGACGACTTGCCCT
TCEA3 Hs00957468_m1 GAAATCGAAGATCATATCTACCAAG
hsa-mir-199a-3p 002304 ACAGUAGUCUGCACAUUGGUUA
hsa-mir-425-5p 001516 AAUGACACGAUCACUCCCGUUGA
hsa-miR-1225-3p 002766 UGAGCCCCUGUGCCGCCCCCAG
hsa-mir-1233-3p 002768 UGAGCCCUGUCCUCCCGCAG
hsa-mir-1274b 002884 UCCCUGUUCGGGCGCCA
hsa-mir-1303 002792 UUUAGAGACGGGGUCUUGCUCU
RNU6B 001093 CGCAAGGATGACACGCAAATTCGTGAAGCGTTCCATATTTTT

RNU6B, B2M, IPO8 and hsa-mir-1274b were used as reference genes. miRNA/miR, microRNA.

Significant Spearman correlation coefficients and corresponding P-values for comparisons between the investigated genes and miRNAs.

Gene miR-199a-3p miR-425-5p miR-1225-3p miR-1233-3p miR-1303
L1TD1 −0.474 (P=0.001) −0.538 (P<0.001) −0.443 (P=0.004) −0.496 (P<0.001) /
SLITRK6 −0.259 (P=0.048) / / / /
ST6GALNAC1 −0.358 (P=0.006) −0.266 (P=0.050) −0.385 (P=0.006) −0.397 (P=0.002) /
TCEA3 / / / / /

miRNA/miR, microRNA; /, not significant.

Significant Spearman correlation coefficients and corresponding P-values for investigated comparisons between the miRNAs.

miRNA miR-199a-3p miR-425-5p miR-1225-3p miR-1233-3p miR-1303
miR-199a-3p 1 0.823 (P<0.001) 0.653 (P<0.001) 0.633 (P<0.001) 0.460 (P<0.001)
miR-425-5p 1 0.587 (P<0.001) 0.592 (P<0.001) 0.314 (P=0.025)
miR-1225-3p 1 0.886 (P<0.001) 0.699 (P<0.001)
miR-1233-3p 1 0.714 (P<0.001)
miR-1303 1

miRNA/miR, microRNA; /, not significant.