Multi-Omics Analysis of the Tumor Microenvironment in Liver Metastasis of Colorectal Cancer Identified FJX1 as a Novel Biomarker

Colorectal cancer incidence and mortality have increased in recent years, with more than half of patients who died of colorectal cancer developing liver metastases. Consequently, colorectal cancer liver metastasis is the focus of clinical treatment, as well as being the most difficult. The primary target genes related to colorectal cancer liver metastasis were via bioinformatics analysis. First, five prognosis-related genes, CTAG1A, CSTL1, FJX1, IER5L, and KLHL35, were identified through screening, and the prognosis of the CSTL1, FJX1, IER5L, and KLHL35 high expression group was considerably poorer than that of the low expression group. Furthermore, the clinical correlation analysis revealed that in distinct pathological stages T, N, and M, the mRNA expression levels of CSTL1, IER5L, and KLHL35 were higher than in normal tissues. Finally, a correlation study of the above genes and clinical manifestations revealed that FJX1 was strongly linked to colorectal cancer liver metastasis. FJX1 is thought to affect chromogenic modification enzymes, the Notch signaling system, cell senescence, and other signaling pathways, according to KEGG enrichment analysis. FJX1 may be a critical target in colorectal cancer metastasis, and thus has the potential as a new biomarker to predict and treat colorectal cancer liver metastases.


INTRODUCTION
Colorectal cancer is one of the most frequent types of cancer in the digestive tract (Huang et al., 2013). According to the American Cancer Society 2018 global cancer statistics, colorectal cancer has an incidence rate of 11%, which ranks third, and a mortality rate of 9%, which ranks second (Bray et al., 2018).
Advanced colorectal cancer tumor cells invade lymphatic vessels and blood vessels from the initial location or are implanted into other sites, causing the cancer to spread and become more serious, for example, metastasis has been observed in the peritoneum (Pretzsch et al., 2019), ovary (Zhou and Ding, 2021), brain (Stewart et al., 2018), liver (de la Pinta et al., 2021), lung (Hou et al., 2017), and bone (Alghandour et al., 2020). The liver is the most common distant metastatic organ of colorectal cancer, with liver metastases observed in 15-25% of patients at diagnosis (Engstrand et al., 2018). Patients with liver metastases who do not receive treatment have a median survival duration of only 6.9 months, and the 5-year survival rate for patients with unresectable liver metastases is low (House et al., 2010). Colorectal cancer is a complicated process involving several stages, phases, and genes (Li et al., 2016), however, our knowledge of the molecular pathways that contribute to colorectal cancer metastasis is limited. Given the significant morbidity and mortality associated with colorectal cancer liver metastases, it is critical to elucidate the disease pathophysiology and metastatic mechanisms to improve colorectal cancer treatment.
MiRNAs may be involved in colorectal cancer liver metastases. Exosome-transporting miRNAs released by colorectal cancer cells mediate tumor cell metastasis in the liver and miRNAs released by hepatocytes promote the formation of metastatic tumors (Balacescu et al., 2018). Due to the hepatocytes' regenerative abilities, a high DNA mutation rate, and cell proliferation rate, they quickly develop into hepatocyte lesions when inflammatory lesions are continuously stimulated. Also, tumor-associated macrophages (TAM) create a favorable environment for colorectal cancer spread. TAM increases angiogenesis around tissues, and accelerates tissue remodeling and tumor cell extravasation, producing ideal conditions for the creation of tumor tissue to help repair wounds (Cassetta and Pollard, 2020;Pan et al., 2020).
The Notch signaling pathway is involved in the activation of hepatic stellate cells, the development of liver fibrosis (Tao et al., 2020), the differentiation and development of TAM (Palaga et al., 2018), the activation of the WNT pathway, and other processes that promote the progression of liver cancer (Moeini et al., 2012). FJX1 affects tumor cell proliferation and differentiation as a downstream target gene of the Notch signaling system. FJX1 is a gene on chromosome 11 (11p13) (Chai et al., 2019) and the expression of FJX1 mRNA and protein are upregulated in colorectal tumor tissues compared to normal intestinal epithelial tissues . FJX1 may promote endothelial cell capillary formation in a hypoxia-inducible factor 1-dependent manner (Al-Greene et al., 2013;Liu et al., 2020), as well as act as a proto-oncogene that affects metastasis and recurrence in colon cancer patients. It was reported that the transcription level of FJX1 was inhibited by miR-1249, thereby reducing the cell proliferation, migration, and invasion of colorectal cancer (Dang and Zhu, 2020).
Bioinformatics analyses were applied to investigate target genes linked to colorectal and liver cancer metastasis. FJX1 was the most clinically important colorectal cancer prognostic gene and affects biological pathways such as chromogenic modifying enzymes and Notch signaling, as well as cellular senescence and other signaling pathways, according to a KEGG enrichment study (Cui et al., 2021). Therefore, FJX1 may be a new target for colorectal cancer liver metastasis.

Raw Data Source and Expression Difference Analysis
GEO provided data on in vivo models of organ-specific colon cancer metastasis, GSE64595. Total RNA from CD110 cells sorted from primary colonic tumors (CRC102-PT and CRC108-PT) and matching hepatic metastases (CRC102-LM and CRC108-LM) were analyzed and data from UCSC xena (https://xenabrowser.net/datapages/) were used to determine differential gene expression in unpaired samples of tumor and normal tissues (Barrett et al., 2013;Vivian et al., 2017). The expression differences were assessed using RNAseq data in level 3 HTSeq-fragments per kilobase per million (FPKM) format in the Colon adenocarcinoma (COAD) project of TCGA (https://portal.gdc.cancer.gov/). The RNAseq data was then transformed from FPKM to TPM (transcripts per million reads) format for log2 conversion using the R software (version 3.6.3) and the R package ggplot2 (version 3.3.3). The results of immunohistochemical staining were utilized to investigate the protein expression of various genes in the HPA database (Uhlen et al., 2015;Thul et al., 2017). The patient information in TCGA is summarized in Supplementary Table S1.

Clinical Significance
The differential genes were utilized to build the Least absolute shrinkage and selection operator (Lasso) model, which was used to search for prognostic molecules using the R packages glmnet (version 4.1-2) and survival (version 3.2-10). After tenfold cross-validation, five molecules with non-zero variables were found, and the risk factor map and ROC curve were created to prove the Lasso prognostic model's accuracy. R and R packages: survminer package (version 0.4.9) and survival package (version 3.2-10) were used for data analysis and cleaning based on TCGA expression and phenotypic data. Finally, numerous survival curves were produced to examine the significance of these prognostic genes in terms of patient survival. A nomogram was used to set scale scores to characterize the conditions of each variable in the multi-factor regression model so that the total score could be calculated to predict the probability of the event, and Calibration examined the fitting of the actual probability and the probability of the model prediction in various conditions based on multi-factor regression analysis. The Cox regression model was built with survival outcome and survival time as dependent variables to investigate the influence of various independent variables on survival (Liu et al., 2018).

Enrichment Analysis
The biological function was explored by single gene differential analysis of the TCGA database, and the patients Frontiers in Genetics | www.frontiersin.org July 2022 | Volume 13 | Article 960954 were divided into two groups according to their FJX1 expression, and then R (version 3.6.3) and R package: DESeq2 (version 1.26.0) were conducted for differential analysis. By simulating the effect of knockdown or overexpression, we explored the genes that FJX1 may affect, further inferring the functions or pathways that may be involved. GO/KEGG enrichment analysis and GSEA for 86 genes with differential significance with FJX1 with | log2(FC)|>2 and p. adj<0.05 were performed and the log2(FC) corresponding to the differential genes was used to calculate the z-score corresponding to each entry as follows:

Immune Microenvironment
The RNAseq data in TCGA were used to calculate the immune infiltration score using two methods. Single sample gene set enrichment analysis (ssGSEA) was conducted to calculate the degree of immune cell infiltration. The other method was estimate, the algorithm content has StromalScore, ImmuneScore, and ESTIMATEScore. In addition, two immune infiltration analysis methods, Ecotyper, and xCell, were used to further analyze the immune characteristics in the GEO dataset as a validation set. R package: GSVA (version 1.34.0) and ssGSEA (GSVA packet built-in algorithm) (Bindea et al., 2013;Hanzelmann et al., 2013).

Statistical Analysis
The results are performed with R software and R packages. Significant identification: ns, p ≥ 0.05; p, p < 0.05; pp, p < 0.01; ppp, p < 0.001.

Screening of Target Genes Related to Colorectal Metastatic Liver Cancer
A flow chart of the full-text logic is provided in Figure 1. Transcriptome data on colorectal cancer liver metastasis were obtained from the GEO database GSE64595. Group 1 consisted of six primary colorectal cancer samples and group 2 consisted of six colorectal cancer samples with liver metastasis. Variance Stabilizing Normalization (VSN) was applied to the data and principal component analysis (PCA) showed that there were significant differences among the 12 samples ( Figure 2A). The Uniform Manifold Approximation and Projection (UMAP) plot obtained from the dimensionality reduction analysis of the observational data also showed differences between the samples ( Figure 2B). In the normalized boxplots of the samples, the samples in the primary and metastatic cancer groups were well normalized ( Figure 2C). These results showed that the independent samples of these 12 primary and metastatic cancers were of good quality, which was statistically significant for grouping the cancer samples. In total, 89 differential genes satisfied the |log2(FC)|>1 & p. adj<0.05 thresholds ( Figure 2D), and the top 40 genes are shown in the heatmap ( Figure 2E). Prognostic Lasso analysis was used to construct and screen prognostic models, revealing that the expression of CTAG1A, CSTL1, FJX1, IER5L, and KLHL35 was higher in metastatic cancers than in primary cancers, indicating that the significantly higher expression of these five genes may drive cancer cell metastasis ( Figures 3A,B). With prognostic data from TCGA to map risk factors, CTAG1A was excluded from further analysis due to the lack of data. The risk factor diagram of the remaining four genes showed that the risk factors increased with transcription levels and poor prognostic outcomes in patients ( Figure 3C). The ROC curves further verified that the AUC of the four genes was greater than 0.5, indicating high diagnostic accuracy of the four genes ( Figure 3D).

Prognosis Analysis
To investigate the correlation between the expression of the four genes and prognosis, a survival curve was drawn based on the Frontiers in Genetics | www.frontiersin.org July 2022 | Volume 13 | Article 960954 clinical data of TCGA, showing that the prognosis of the group highly expressing CSTL1, FJX1, IER5L, and KLHL35 was significantly worse than that of the low expression group ( Figures 4A-D). The calibration plots to determine how well the model predicted the actual outcome showed the strong predictive power of these four genes to predict 1-, 3-, and 5year survival ( Figures 4E,F). A nomogram was then constructed to further evaluate patient prognosis ( Figure 4G).

Gene Transcription and Protein Expression
The protein expression in the paired and unpaired samples of TCGA colorectal cancer was also investigated. The Mann-Whitney U test (Wilcoxon rank-sum test) for unpaired samples showed the expression of CSTL1, FJX1, IER5L, and KLHL35 in tumor tissue was higher than that in normal tissue ( Figures 5A-H). Similarly, the paired samples t-test also showed that expression levels in tumor tissues were much higher than in normal tissues (paired and unpaired samples, P CSTL1 , P FJX1 , P IER5L , and P KLHL35 < 0.001). The immunohistochemical staining showed that the expression of FJX1 in colorectal cancer tissues was significantly different from that in normal colon tissues, while the protein level of KLHL35 was not significantly different ( Figures 5I-L). However, no immunohistochemical data were available for the remaining two genes. In addition, FJX1 was found to be localized in vesicles within the cells (Figures 5M-P).

Clinical Relevance
The clinical relevance of CSTL1, FJX1, IER5L, and KLHL35 in 521 colon cancer cases was assessed, showing the expression level of CSTL1 mRNA was higher than that in normal tissues, and the expression levels of PD, SD, PR, and CR were higher than those in normal tissues (Supplementary Figure S1). Similarly, the expression of FJX1 mRNA in pathological stage T1, T2, T3, and T4 of colorectal cancer was higher than that in normal tissue, and the expression of lymphoid-infiltrating tumor tissue was higher than that of non-infiltrating tissue, and the expression level of colon adenocarcinoma was higher than that of normal tissue. Expression levels were higher in PD, SD, PR, and CR stages for the primary treatment outcome than in normal tissues, and FJX1 expression levels were higher in DSS, OS, and PFI in the death group (Supplementary Figure  S2). The same situation appeared in the other two genes (Supplementary Figures S3, S4).
Promoter methylation levels of 313 colon adenocarcinoma samples and 37 normal tissue samples were analyzed from the TCGA database using the UALCAN platform ( Figures 6A-D), indicating that CSTL1 promoter methylation levels in primary malignancies were lower than normal tissue. The promoter methylation levels of FJX1, KLHL35, and IER5L in primary malignant tumors were higher than those in normal tissues. Subsequently, we compared the promoter methylation levels of four genes for TP53-mutated (n = 174) and non-TP53-mutated colon adenocarcinomas (n = 136) ( Figures 6E-H), showing that the promoter methylation level of CSTL1 was lower in tumor tissues with and without TP53 mutation than in normal tissues. The promoter methylation levels of FJX1, KLHL35, and IER5L were higher in TP53 and non-TP53-mutated tumor tissues than in normal tissues, and the methylation levels of FJX1 in non-TP53-mutated tumors were higher than those in TP53-mutated tumors. In addition, DNA methylation of the genes was validated in the EWAS database, the same as the UALCAN database ( Figures 6I-L). Finally, we performed a survival analysis of the promoter methylation levels of the different genes ( Figures 6M-P), showing that the survival rate of 148 samples with CSTL1 methylation levels greater than or equal to 0.863 was significantly higher than that of 143 samples with methylation levels less than 0.863. However, according to the hypomethylation of CSTL1 in the tumor, which was contrary to the survival curve, this gene was not considered a biomarker of colorectal cancer liver metastasis.

FJX1-Related Biological Pathways
By analyzing the correlation between CSTL1, FJX1, IER5L, and KLHL35 genes and clinical manifestations, FJX1 was the most clinically relevant prognostic gene, so follow-up analyses were performed to explore the pathways and biological functions affected by FJX1. First, single gene enrichment analysis was performed and FJX1 downstream target genes were obtained by mock knockout or overexpression of FJX1. Then, GO/ KEGG was performed on these downstream differential genes, and 122 GO pathway-enriched biological entries and 6 KGEE pathway-enriched biological entries were obtained ( Figures 7A-C). The biological processes involved mainly included nucleosome organization, chromatin assembly, and nucleosome assembly. The cellular composition mainly included the protein-DNA complex, DNA packaging complex, and the nucleus. The molecular function involves nucleosome DNA binding, taste receptor activity, and bitter taste receptor activity. The results of the KEGG enrichment analysis were mainly systemic lupus erythematosus, alcoholism, and taste transduction. Based on the enrichment analysis, the logFC corresponding to the extracted molecules was used to calculate the z-score corresponding to each entry. The z-scores for these entries were all negative, implying that these pathways may be inhibited ( Figures 7D-G). Meanwhile, enrichment analysis of the GSEA gene set showed that the expression of FJX1 significantly inhibited the functions of ten

Analysis of the Tumor Microenvironment
Differences in the immune-infiltrating microenvironment of primary and metastatic cancers based on data from TCGA-COAD revealed a positive correlation between FJX1 and the level of infiltration of most immune cells. As the transcription level of FJX1 increases, so does the level of immune cell infiltration (Supplementary Figure  S5). Cellular heterogeneity in the tumor microenvironment is an emerging area of research. xCell was used to perform cell-type enrichment analysis of the gene expression data from 64 immune and stromal cell types, revealing no significant differences in most immune or stromal cells in the six primary and six metastatic samples. Interestingly, the ImmuneScore was significantly increased in metastatic cancers and was almost absent in primary cancers. In addition, the B-cell infiltration score was also significantly increased as was the ImmuneScore, which was similar to the MicroenvironmentScore, suggesting that the coelevation of the microenvironment and ImmuneScore in metastatic cancer was caused by a marked elevation of B cells. B cells are the main immune cells of humoral immunity and the expression of FJX1 was proportional to the infiltration level of B cells. Furthermore, cDCs, iDCs, class-switched memory B cells, and memory B cells were all significantly increased in metastatic cancers compared to primary cancers, whereas the infiltration of MEPs and osteoblasts was significantly reduced. The ImmuneScore in TCGA was calculated using an estimation algorithm showing that the expression levels of the ImmuneScore and FJX1 were also positively correlated. Ten possible cancer ecotypes in each sample were calculated and five CE values were found to be consistent with our sample: CE1, CE2, CE6, CE9, and CE10. It can be seen that the cancer ecotypes of primary cancer are CE6 and CE10, while the cancer ecotypes of metastatic cancer are CE1, CE2, and CE9 (Supplementary Figures S6, S7).

DISCUSSION
Colorectal cancer is one of the most common cancers, with a high rate of morbidity and mortality (Brody, 2015;Castells, 2016;Nasseri and Langenfeld, 2017). Each year, more than 1.4 million new cases of colorectal cancer are diagnosed, resulting in around 700,000 deaths (Jin et al., 2020). Colorectal cancer can only be cured through surgical excision of the lesion but the current therapy benefit for advanced metastatic colorectal cancer is modest (Roncucci and Mariani, 2015). The liver is the most prevalent site of colorectal cancer metastatic spread and although simultaneous liver resection can help to delay the progression of colorectal Frontiers in Genetics | www.frontiersin.org July 2022 | Volume 13 | Article 960954 cancer, the prognosis for this group of patients is still dismal (Chakedis and Schmidt, 2018). The lack of early indicators to assess colorectal cancer liver metastases is of concern. We screened the differentially expressed genes with the Lasso prognostic model after analyzing the metastatic and primary data of colon cancer and found five genes with non-zero variables. Prognosis data for CTAG1A was missing, so it was excluded from the subsequent prognosis analysis. The survival curves of the remaining four genes were then verified, and all of them exhibited considerable predictive value. HR > 1 also suggested that they were risk factors for three prognosis groups of patients. The association between factors in the prediction model and the accuracy of the prognosis model was depicted by nomograms and calibration charts. The elevated expression of four genes in tumors was confirmed in both the unpaired and paired samples in TCGA. FJX1 protein levels in the tumor were higher than in normal tissue, while KLHL35 protein levels were undetectable in both tumor and normal tissues, and no immunohistochemistry data for the remaining two genes were available.
Regarding the methylation data of the four prognostic genes, except for the CSTL1 promoter area, the degree of methylation of the remaining three genes dramatically increased. Simultaneously, when CSTL1 was divided into subgroups based on the degree of methylation, the hypermethylation subgroup had a worse prognosis despite two major databases showing that the degree of methylation in the normal group was higher, so CSTL1 was also eliminated. Then, the univariate/multivariate Cox regression revealed that only the FJX1 subgroup was significant among the four prognostic genes, with the low expression group serving as the control, and the HR was 1.904, indicating that FJX1 expression was an important hazard factor for patients. Clinical analysis revealed the highest correlation between FJX1 and clinical factors, nevertheless, FJX1 was more significantly expressed in metastatic liver cancer than in the original colon cancer from GEO. The logFC in the volcanic diagram is greater than 1, indicating that FJX1 is most likely the determining factor in tumor spread. FJX1 was found to be localized in vesicles, indicating that it is likely to release cells along with the vesicles. In addition, FJX1 expression is directly proportional to the level of infiltration of most immune cells, and the lmmuneScore increased with high FJX1 expression, suggesting that an elevation in FJX1 could The z-scores for these entries were all negative, implying that these pathways may be inhibited.
Frontiers in Genetics | www.frontiersin.org trigger an inflammatory response that accelerates tumor growth. The infiltration scores and alterations of 64 types of immune cells were calculated using xCell analysis of GSE64595 expression profile data. The immunological score and the microenvironment score both increased dramatically according to the heat map. Th2 cells were inversely proportional to FJX1 expression in metastatic cancer, according to TCGA data. Surprisingly, the heat map also revealed that as FJX1 expression increased, the number of infiltrating Th2 cells decreased. Furthermore, early cancer carcinoma ecotypes include CE6 and CE10, whereas metastatic cancer ecotypes have CE1, CE2, and CE9, indicating that the tumor immune milieu becomes more complicated as FJX1 expression increases. CE1, CE2, and fibroblasts are microenvironmental cells that play a role in the establishment of the immune microenvironment in metastatic cancer. However, xCell shows that in metastatic cancer, there are more infiltrating fibroblasts. To summarize, these findings suggest that FJX1 stimulates the formation of tumor fibroblasts and is critical in the progression of metastatic cancer. FJX1 mRNA and protein are upregulated in human colorectal tumor epithelium compared to rectal adenomas, and high expression of FJX1 is linked to a poor patient prognosis (Al-Greene et al., 2013), in line with our observations. FJX1 is an angiogenesis regulator whose expression levels are elevated in a variety of cancers. FJX1 protein levels in the endometrium of women without endometriosis do not fluctuate during the menstrual cycle, however, during normal endometrial secretion, FJX1 levels in women with endometriosis were considerably greater than in those without endometriosis, implying that elevated FJX1 protein levels may play a role in endometriosis etiology (Chang et al., 2018).
The tumor microenvironment and colorectal cancer liver metastases are inextricably linked. Hepatic fibroblast function changes increase the establishment of metastatic foci of colorectal cancer cells in the liver, and high production of microRNA-10a in the tumor microenvironment can block the activation of hepatic fibroblasts, preventing colorectal cancer liver metastasis . Exosomes have a key role in the establishment of the pre-metastatic microenvironment and tumor cell-derived exosomes can boost the pre-metastatic microenvironment by stimulating the pre-metastatic microenvironment in mice models of pancreatic cancer liver metastasis and melanoma lung metastasis. Exosomes speed up tumor metastasis, and exosomes from intestinal cancer cells might stimulate the establishment of a pre-metastatic milieu in liver tissue, speeding up tumor spread (Lai et al., 2020). Also, hypoxic and acidic microenvironments are both linked to colorectal cancer liver metastasis (Larionova et al., 2020).

CONCLUSION
In conclusion, FJX1 may not only be a cause of colorectal cancer but also a risk factor for liver metastasis as a prognostic biomarker.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
JZ designed the study. JZ, YH, and HZ were involved in the database searches and statistical analyses. JZ and HZ were involved in the writing of the manuscript and its critical revision. JZ and HZ were responsible for the submission of the final version of the paper. All authors approved the final version. All authors agree to be accountable for all aspects of the work.

FUNDING
This work was funded and supported by the Project of Wannan Medical College (WK 2019F44) and the Health Research Project of Huai'an City (HAWJ201930).

ACKNOWLEDGMENTS
We are very grateful for data provided by databases such as TCGA, GEO.