Identifying metastasis-initiating miRNA-target regulations of colorectal cancer from expressional changes in primary tumors

Colorectal cancer (CRC) is prevalent with high mortality, with liver metastasis contributing as a major factor that worsens the survival of patients. The roles of miRNAs in CRC have been elucidated, subsequent to recent studies that suggest the involvement of miRNAs in cancer biology. In this study, we compare the miRNA and gene expression profiles of primary tumors between two groups of patients (with and without liver metastasis) to identify the metastasis-initiating microRNA-target gene regulations. Analysis from 33 patients with metastasis and 14 patients without metastasis revealed that 17 miRNAs and their 198 predicted target genes are differentially expressed, where the target genes showed association with cancer progression and metastasis with statistical significance. In order to evaluate the clinical implications of the findings, we classified CRC patients of independent data into two groups based on the identified miRNA-target regulations, where one group was closer to primary tumors with metastasis than the other group. The comparison of survival showed statistically significant difference, thereby implying the roles of the identified miRNA-target regulations in cancer progression and metastasis. The identification of metastasis-initiating miRNA-target regulations in this study will lead to better understanding of the roles of miRNAs in CRC progression.

Scientific RepoRtS | (2020) 10:14919 | https://doi.org/10.1038/s41598-020-71868-0 www.nature.com/scientificreports/ with small RNA-sequencing to classify miRNAs that predict metastasis of colon cancers 6 . Pizzini et al. obtained miRNA and mRNA expression data with microarray profiling, and compared their profiles in primary and metastatic tumors, and normal colon tissues 7 . Röhr et al. analyzed high-throughput sequencing data of miRNAs and mRNA transcripts from normal colon and metastasis tumor samples 8 , and they evaluated the up/downregulations of miRNAs from individual patients to identify candidates of therapeutic targets. Although these studies identified miRNAs or their target regulations associated with metastasis, they have certain limitations. Lin et al. 5 and Beumer et al. 6 investigated the expression of miRNAs but did not consider gene expressions that can be potential targets of miRNAs. Pizzini et al. 7 analyzed the expressions of miRNAs and genes together, but their study has limited miRNA coverage as they used microarrays that do not fully represent the currently known miRNAs. Since Beumer et al. 6 used miRNA expression profiles from small RNA-seq, their study covered most of the known miRNAs, but they did not consider mRNA expression profiles from matched samples. Another study 8 compared miRNA and mRNA expression profiles from primary colon tumor tissues and liver metastatic tumor tissues, but was more focused on identifying differentiated miRNA-gene regulations between primary and metastatic tumors rather than identifying miRNA-gene regulations that initiate the liver metastasis of CRC.
In this study, we use miRNA expression data obtained from small RNA-seq, together with whole transcriptome sequencing (WTS)-based gene expression profiles from matched tumor tissues. By investigating miRNA and mRNA profiles from primary tumor tissues with and without accompanying liver metastasis, we aim to identify miRNAs and their potential target genes that show altered expressions and may contribute to the initiation of liver metastasis. Our study will provide more convincing knowledge about the roles of miRNA-gene regulations that initiate liver metastasis in CRC.

Methods
Sample preparation and expression profiling. Primary tumor tissues were obtained from 47 CRC patients and prepared for miRNA and mRNA expression profiling. Of the 47 patients, 14 showed no metastasis, while 33 patients had accompanying liver metastasis. The detailed demographic and clinical characteristics of the 47 patients are given in Table 1. The two groups of patients did not show statistically significant differences in ages, sex, tumor locations, and the degree of cancer cell differentiation (See Supplementary Table S1 online for detailed information). All samples were surgically acquired before any treatment. Both of small RNA-seq and whole transcriptome-sequencing were conducted for these 47 primary tumor samples on Illumina HiSeq 2000. Whole transcriptome-sequencing was done using the Illumina TruSeq RNA Sample Preparation Kit v2, and paired-end sequencing with a read length of 100 bp and targeted read depth of 50 million reads/sample was performed. On average, the Q30 sequence quality was 93.26% and 87.17% of data was mapped to genome. Small RNA-seq was conducted as single-end sequencing with a read-length of 50 bp and targeted read depth of 10 million reads/sample. On average, the Q30 sequence quality was 94.66% and 86.71% of data was mapped to genome. Sequences were aligned to the human reference genome (hg19) with the Burrows-Wheeler Aligner (https :// bio-bwa.sourc eforg e.net/). Expression values were normalized using RNA-Seq by Expectation-Maximization Identifying metastasis-initiating miRNA-target regulations. In order to identify miRNAs and genes that show differential expression between patients with and without liver metastasis, we compared the expressions of miRNAs and genes between primary tissues with liver metastasis and without metastasis. Differential expressions between the two groups were evaluated with raw transcript read counts using the edgeR 16 R package, one of the most popular software tools for identifying differentially expressed RNAs from sequencing data, and the resultant p values of statistical significance were adjusted with the Benjamini-Hochberg method 17 .
In this study, a miRNA or a gene is determined to be differentially expressed (DEmiRNA or DEG) if it shows greater than or equal to twofold change in the average expression level between the tumors with and without metastasis, and its adjusted p value is less than 0.05. A power analysis suggested more stringent criteria for differential expressions, especially for miRNAs, but we maintained these thresholds as more stringent criteria were too restrictive in our data configuration and considering the technical difficulties in practice of obtaining miRNA sequencing data of similar amount with RNA-seq data (see Supplementary Figs. S1-S3 online). Based on the DEmiRNAs and DEGs obtained, we identified differentially expressed miRNA-target regulations as metastasis-initiating, by utilizing predicted miRNA-target gene information. The predicted human miRNAtarget information from TargetScan v7.2 18 is used in this study, which encompasses 120,702 predicted human miRNA family-target regulations. For each DEmiRNA, its predicted target genes are identified among the DEGs, and these predicted regulations between the DEmiRNA and its target DEGs constitute miRNA-target regulations that contribute to the initiation of CRC liver metastasis. The overview of entire processes is illustrated in Fig. 1.
Evaluating the metastasis-specificity of identified miRNA-target regulations. We evaluated the metastatic association and statistical significance of the identified miRNA-target regulations. In order to evaluate the association with metastasis, we applied the functional enrichment test to the identified target genes to evaluate their role in biological functions related to metastasis. Using the MSigDB 19 , 22 hallmark gene sets (Supplementary Table S2 online) were selected as representative factors related to metastatic biological functions. Hypergeometric p values were evaluated for each gene set with the identified target genes from the identified metastasis-initiating regulations, and adjusted with the Benjamini-Hochberg method 17 . Adjusted p value www.nature.com/scientificreports/ less than 0.05 for a gene set is indicative of statistically significant enrichment of the function corresponding to the gene set in the target genes. The association of identified target genes to metastasis can be represented by the number of enriched metastasis-related functions from the 22 gene sets: higher number of associated metastasis-related functions in target genes is indicative of greater association with metastasis. We hypothesized that if the identified miRNA-target regulations represent true biology related to metastasis, the target genes should be significantly enriched with metastasis-related functions than by random chance, by showing higher metastasis-relatedness. To verify this hypothesis, we applied a random permutation approach to test if the evaluated metastasis-relatedness of the identified miRNA-target regulations is greater than that obtained by random chance. For each random permutation, the metastasis status of expression samples was randomly shuffled and re-assigned, and the entire process (from identifying metastasis-related miRNA-target regulations to evaluating the metastasis-relatedness) was conducted for the randomly permuted data. By repeating these random permutations for a fixed number of times (1,000 in this study) and examining the cases showing greater than or equal to the original metastasis-relatedness, we estimated the p value of the original metastasis-relatedness. More detailed descriptions of this permutation test are given in Supplementary Fig. S4 online with a schematic diagram. where µ is the average, and σ is the standard deviation of z-transformed expressions from the sample group. In these equations, Eq. (1) is applied when there is no overlap between the two representative ranges R woM and R M , and qs. (2) and (3) are applied when there is overlap between these ranges (see Supplementary Fig. S5 online for further descriptions with schematic illustrations). For each sample from the independent TCGA data, each miRNA or gene from identified regulations is categorized as having similar expression to that of a sample group from our data set, if its z-transformed expression level from the TCGA sample falls within the representative range from our data set for the same miRNA or gene. By repeating this process for all miRNAs and genes of a TCGA case, the number of miRNAs and genes that resemble the expression ranges of primary tumors with metastasis and without metastasis can be counted. A TCGA case is predicted to resemble the primary tumors with metastasis of our data set if it has significantly more miRNA-target genes that have expressions similar to primary tumors with metastasis than miRNAs-target genes that have expressions similar to primary tumors without metastasis (by comparing the number of genes using the chi-square test with a p value threshold = 0.05); such TCGA cases constitute a separate group G 1 (vice versa for TCGA cases resembling primary tumors without metastasis, duly categorized as group G 2 ). A TCGA case is excluded from the analysis if the case has significantly more miRNAs (by the chi-square test with a p value threshold = 0.05) that resemble the opposite sample group and contradicting the miRNA-target expression patterns (for example, a case categorized as G 1 by genes but G 2 by miRNAs). After categorizing the independent TCGA cases into two groups representing primary tumors with and without metastasis, the survival periods of patients were compared, to recognize potential relevance with the identified metastasis-initiating miRNA-target regulations. Log-rank test 21 is used for this survival comparison, and the multivariate Cox proportional hazards (Cox-PH) model is used to adjust hazard ratios for available clinical variables from TCGA.

Identification of metastasis initiating miRNA-target regulations.
Comparing the expression levels of miRNAs and genes of the two primary tumor sample groups, with and without metastasis, we identified 58 DEmiRNAs and 639 DEGs ( Supplementary Fig. S6 online and Fig. 2). Of these, 516 metastasis-initiating if µ woM < µ wM and µ woM + σ woM < µ wM − σ woM or µ woM > µ wM and µ woM − σ woM > µ wM + σ woM (1)  Fig. 4. We observed that miRNA-130b is down-regulated and ZEB1 is up-regulated in primary tumors with metastasis (Fig. 4a). It is known that the expression of ZEB1, a key gene associated with epithelial-mesenchymal transition, increases subsequent to deregulation of miR-130b 22 . The identified miR-132-HB-EGF regulation (Fig. 4b) has also been reported in prostate cancer cell lines 23 , where the HB-EGF, associated with cell migration and invasion, is inhibited by miR-132. We also identified the inhibition of FGF7 (a well-known fibroblast growth factor involved in tumor growth and invasion, and a potent epithelial cell-specific growth factor) by miR-15a (Fig. 4c); a previous study reported that the down-regulation of miR-15a progresses pancreatic cancer 24 . Considering the disturbed inhibition of FGF7 by miR-15a in primary tissues with metastasis, it is interesting that we also identified the miR16-FGFR regulation as a metastasis-initiating miRNA-target regulation. Our result shows down-regulation of miR16 and corresponding up-regulation of FGFR in primary tumors with metastasis (Fig. 4d), where FGFR is a receptor for FGF7 and well known to be related to proliferation and migration of endothelial cells and angiogenesis 25 . The regulation of miR-30b-SERPINE1 (Fig. 4e) is also reported to be related to gastric cancer tumorigenicity 26 , where restoring the regulation of SERPINE1 by miR-30b increases apoptosis of cancer cells in vivo. www.nature.com/scientificreports/ Some of the identified miRNA-target regulations from our results (Fig. 5) include genes previously identified as relevant to CRC, with experimental validations of their regulations provided in corresponding studies. FGFR1 is reported to be associated with cell proliferation, differentiation and angiogenesis in cancer 27 . Previous studies have shown that miR-424 inhibits FGFR1 [28][29][30] , and the regulation of miR-424-FGFR1 was predicted to be metastasis-initiating from our results with down-regulated miR-424 and up-regulated FGFR1 (Fig. 5a). The PDGFRB gene plays an important role in proliferation and migration of various tumors 31 , and studies have reported that PDGFRB represses miR-30b 32 , which is also seen in this study (Fig. 5b). AKT3 is an oncogene associated with the regulation and metastasis of tumor cell survival 33 . We identified that regulation of AKT3 by miR-15a, miR-15b and miR-16 is disturbed in primary tumors with metastasis (Fig. 5c), where the miRNAs are down-regulated and AKT3 is up-regulated; this disturbed regulation of AKT by miRNAs with expressional changes has also been reported previously 34,35 . Identified miRNA-target regulations have a statistically significant relationship with metastasis. In order to evaluate the metastasis-relatedness of the identified miRNA-target regulations, we performed the functional enrichment test, as described in the Methods section. Functional enrichment analysis revealed that eight of the 22 hallmark gene sets related to cancer progression and metastasis are enriched in the target genes of the identified miRNA-target regulations (Fig. 6a), and this amount of enriched gene sets reflects the degree of metastasis-relatedness. A random permutation approach was applied 1,000 times, to evaluate its statistical significance by analyzing data with randomly shuffled sample group labels (with metastasis and without metastasis). Among the 1,000 randomly permuted instances, there was no case with metastasis-relatedness higher than or equal to the metastasis-relatedness of the original data, corresponding to p value < 0.001 (Fig. 6b). This statistical significance indicates that the identified miRNA-target regulations are significantly associated with the status of liver metastasis. Different patients' survivals based on the patterns of the identified miRNA-target regulations. From the processes described in Methods, 426 primary tumor samples of TCGA CRC data were cat- www.nature.com/scientificreports/ egorized into two groups, based on the expression similarity with the identified metastasis-initiating miRNAtarget regulations; one group has greater similar expression patterns to primary tumors with metastasis than the other group. Survival comparison between the two groups using the log-rank test revealed significantly differing survival rates (p value = 0.009, Hazard ratio = 1.44 (CI 1.09-1.91)). We also investigated whether this difference in survival is independent of other clinical characteristics by using the multivariate Cox-PH model. The adjusted hazard ratio with gender, race and age was 1.35 (CI 1.06-1.90, p value = 0.012), thus representing this survival difference is independent of other clinical variables. The group of patients with expression patterns more similar to the metastasis-initiating miRNA-target regulations showed median survival of 761 days, whereas the median survival of the other group was 1,005 days (Fig. 7). This significant difference of prognosis from analyzing independent data by the identified miRNA-target regulations implies that the identified regulations play important roles in the progression of CRC.

Discussion
In the current study, we compared the expression levels of miRNAs and genes of primary tumors from CRC patients, with and without metastasis, to identify metastasis-initiating miRNA-target regulations. We identified 516 metastasis-initiating miRNA-target regulations, including 17 miRNAs and 198 target genes. Statistical Evaluated statistical significance of metastasis-relatedness of the identified miRNA-target regulations by random permutation test (1,000 random permutations of sample group labels). All 1,000 instances of identifying metastasis-initiating miRNA-target regulations from randomly permuted data sets show less metastasis-relatedness than that of the identified miRNA-target regulations from the original data set (equivalent to p value < 0.001).
Scientific RepoRtS | (2020) 10:14919 | https://doi.org/10.1038/s41598-020-71868-0 www.nature.com/scientificreports/ analysis revealed that the identified metastasis-initiating miRNA-target regulations are strongly related to cancer progression and metastasis, with statistically significant difference. We also evaluated the clinical significance of the identified miRNA-target regulations by separating primary tumor samples from the TCGA CRC data set, based on their resemblance to the identified regulations. From the survival comparison of the two categorized patient groups, the group that exhibited more resemblance to the identified metastasis-initiating miRNA-target regulations had shorter survival periods than the group showing lesser resemblance. This analysis reveals the important clinical implications of the identified miRNA-target regulations, especially considering that liver metastasis is a major factor that affects the survival of CRC patients. The identified metastasis-initiating miRNAtarget regulations include molecular mechanisms that were reported to be associated with cancer progression and metastasis in different cancer types by previous studies, as well as regulations with well-known genes for CRC. The patterns in expression changes of miRNAs and their regulated genes in our data analysis are also consistent with that of previous findings. Some other novel regulations identified from the current study are potential candidates for future experimental validations. While our results present promising metastasis-initiating miRNAtarget regulations, it has limitation of using only expression data for miRNAs and genes. Other genetic features that can contribute to the metastasis of colorectal cancer include DNA somatic mutations and copy number variations (CNVs). Hence, by evaluating somatic mutations and CNVs along with expressions for analysis, further coverages of regulations between miRNA and target genes can be considered in progressing tumors and initiating metastasis. Another limitation is that we used predicted miRNA-target relationship information in identifying target genes of miRNAs, including many miRNA-target relationships that have not yet been experimentally verified. Nevertheless, we believe that the identified miRNA-target regulations from this study can provide meaningful clues and guidance for researching and experimentally validating metastasis-initiating miRNA-target regulations for CRC.

Data availability
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.