SEMA6B Overexpression Predicts Poor Prognosis and Correlates With the Tumor Immunosuppressive Microenvironment in Colorectal Cancer

Background: Semaphorin 6b (SEMA6B) is a member of the semaphorin axon-guidance family and has been demonstrated to both induce and inhibit tumor progression. However, the role of SEMA6B in colorectal cancer (CRC) has remained unclear. This study sought to explore the promising prognostic biomarker for CRC and to understand the expression pattern, clinical significance, immune effects, and biological functions of SEMA6B. Methods: SEMA6B expression in CRC was evaluated via multiple gene and protein expression databases and we identified its prognostic value through The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Correlations between SEMA6B expression and components of the tumor immune microenvironment were analyzed by packages implemented in R, Tumor Immune Estimation Resource (TIMER), Gene Expression Profiling Interactive Analysis (GEPIA), and Tumor-Immune System Interactions database (TISIDB). RNA interference was performed to silence the expression of SEMA6B to explore its biological roles in the colon cancer cell lines HCT116 and LoVo. Results: The messenger RNA (mRNA) level of SEMA6B and the protein expression were higher in CRC tissues than adjacent normal tissues from multiple CRC datasets. High SEMA6B expression was significantly associated with dismal survival. Multivariate Cox regression analysis demonstrated that SEMA6B was an independent prognostic factor for progression-free survival (PFS). The nomogram showed a favorable predictive ability in PFS. Functional enrichment analysis and the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm revealed that the gene cluster associated with the high SEMA6B group were prominently involved in immune responses and inflammatory activities. Notably, SEMA6B expression was positively correlated with infiltrating levels of CD4+ T cells, macrophages, myeloid-derived suppressor cells (MDSCs), regulatory T cells (Tregs), neutrophils, and dendritic cells. Moreover, SEMA6B expression displayed strong correlations with diverse marker sets of immunosuppressive cells in CRC. Integrative analysis revealed that immunosuppressive molecules and immune checkpoints were markedly upregulated in CRC samples with high SEMA6B expression. Furthermore, knockdown of SMEA6B in colon cancer cells significantly inhibited cell proliferation, migration, invasion and reduced the mRNA levels of immunosuppressive molecules. Conclusion: Our findings provide evidence that high SEMA6B expression correlated with adverse prognosis and the tumor immunosuppressive microenvironment in CRC patients. Therefore, SEMA6B may serve as a novel prognostic biomarker for CRC, which offers further insights into developing CRC-targeted immunotherapies.


Results:
The messenger RNA (mRNA) level of SEMA6B and the protein expression were higher in CRC tissues than adjacent normal tissues from multiple CRC datasets. High SEMA6B expression was significantly associated with dismal survival. Multivariate Cox regression analysis demonstrated that SEMA6B was an independent prognostic factor for progression-free survival (PFS). The nomogram showed a favorable predictive ability in PFS. Functional enrichment analysis and the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm revealed that the gene cluster associated with the high SEMA6B group were prominently involved in immune responses and inflammatory activities. Notably, SEMA6B expression was positively correlated with infiltrating levels of CD4 + T cells, macrophages, myeloidderived suppressor cells (MDSCs), regulatory T cells (Tregs), neutrophils, and dendritic cells. Moreover, SEMA6B expression displayed strong correlations with diverse marker sets of immunosuppressive cells in CRC. Integrative analysis revealed that

INTRODUCTION
Colorectal cancer (CRC) is a major cause of cancer mortality around the word. Approximately 20% of patients occurred metastases at diagnosis (Garborg et al., 2013;Siegel et al., 2014). Although many advances in systemic therapy and liverdirected treatments made thus far, 5 years survival rate is only 12-14% in patients with metastatic CRC (Siegel et al., 2017;Bray et al., 2018). Over the past few years, full recognition of the complex interactions between cancer cells and the immune system has led to a rapid development in immunotherapeutic approaches. Immunotherapeutic strategies include immune checkpoint inhibitors, cancer vaccines, adoptive cell transfer, oncolytic viral therapy, and carcinoembryonic antigen (CEA) T-cell bispecific antibodies, which focus on selectively enhancing the host immune system to fight cancer (Kakimi et al., 2017;Szeto and Finley 2019;Wrobel and Ahmed 2019).
Immune checkpoint inhibitors currently represent the main domain of immunotherapy and have achieved clinical benefits to patients with advanced cancer including renal cell carcinoma, non-small cell lung cancer and malignant melanoma (Reck et al., 2016;Atkins et al., 2017;Gettinger et al., 2018;Ribas and Wolchok 2018). Recent success in using antibodies against various immune checkpoints such as cytotoxic T-lymphocyteassociated antigen 4 (CTLA-4), programmed death 1 (PD-1), and programmed death ligand 1 (PD-L1) for cancer immunotherapy has brought this approach being implemented as a new treatment modality for CRC, especially in terms of targeting the microsatellite instability-high (MSI-H) phenotype (Le et al., 2015;Overman and McDermott 2017;Overman et al., 2018;Le et al., 2020). However, clinical immunotherapeutic trials have revealed that anti-CTLA-4 monoclonal antibodies (mAbs) yield unsatisfactory clinical efficacies in unselective CRC patients (Chung et al., 2010), and anti-PD-L1 mAbs and anti-PD-1 have shown little or no response rates in metastatic CRC (mCRC) (Le, Uram and Wang 2015;Overman and McDermott 2017). Although there is clear clinical evidence for a therapeutic role of immune checkpoint inhibitors in deficient mismatch repair (dMMR) or MSI-H mCRC, the majority of mCRC patients with proficient MMR (pMMR) or microsatellite stable (MSS) phenotypes do not benefit from this type of immunotherapy (Koi and Carethers 2017;Ciardiello and Vitiello 2019;Ganesh et al., 2019;Judge et al., 2020;Morse et al., 2020). Furthermore, previously described molecular features, such as immunoscore, PD-1, PD-L1, MSI, mutational load, and consensus molecular subtypes have not been identified in predicting responses to immune checkpoint inhibitors based on immunotherapy. (Emambux and Tachon 2018;Sveen et al., 2020). Therefore, it is required to discover novel biomarkers with latent prognostic value and screen immune-based therapeutic targets for CRC patients.
Semaphorin family members were initially characterized as axon-guidance factors with functions in axonal navigation, but have subsequently also been linked to the pathology of various diseases, such as cancer, immune disease and neurodegenerative disease (Műzes and Sipos 2014;Neufeld and Mumblat 2016;Franzolin and Tamagnone 2019). Accumulated studies have shown that some semaphorins-including semaphorin 3E (SEMA3E), SEMA4D, SEMA5A, SEMA6D, and SEMA7A-play vital roles in tumorigenesis and tumor development by promoting angiogenesis and tumor-cell migration, as well as the epithelial-mesenchymal transition (EMT); in contrast, SEMA3A, SEMA3B, and SEMA3F exhibit tumor-inhibitory effects (Neufeld et al., 2012;Neufeld and Mumblat 2016;Gurrapu and Tamagnone 2019). Mechanisms that account for the diversity of semaphorin signaling responses in different cellular contexts can profoundly affect these different biological activities. SEMA6B, a member of the semaphorin axonguidance family, has recently been investigated in terms of human SEMA6B gene expression and its roles in cancer. In breast cancer tissues, the SEMA6B promoter undergoes abnormal methylation, and downregulation of SEMA6B messenger RNA (mRNA) has been found in tumor samples (D'Apice and Costa 2013; Kuznetsova et al., 2007). In CRC patients, miR-30b could mediate axon guidance and is significantly negatively correlated with SEMA6B (Coebergh van den Braak et al., 2018). Among different human cell lines, high levels of SEMA6B mRNA have been observed in MCF-7 breast adenocarcinoma cells, and these levels have been found to be downregulated by 9-cis-retinoic acid, an anti-proliferative and differentiation-promoting agent (Murad and Collet 2006). Functionally, SEMA6B has been found to exert complex roles in the development and progression of tumors such as breast cancer (D'Apice and Costa 2013; Murad and Collet 2006), glioblastoma (Kigel and Rabinowicz 2011), gastric cancer (Ge and Li 2013), and testicular cancer (Ji and Wang 2020). However, the prognostic value of SEMA6B in CRC and the relationship between SEMA6B and immune responses remain elusive.
At present study, we assessed SEMA6B expression and clarified its potential prognostic value in CRC patients using The Cancer Genome Atlas (TCGA), Human Protein Atlas (HPA), and Gene Expression Omnibus (GEO) databases. Moreover, we explored the underlying biological functions and relevant pathways of SEMA6B and investigated correlations of SEMA6B with a variety of tumor-infiltrating immune cells (TIICS) as well as tumor-immunity status via comprehensive bioinformatic analyses. Taken together, our present findings may help to uncover prominent immunoregulatory roles of SEMA6B in the CRC microenvironment, and provide a promising biomarker and target for CRC diagnosis and immunotherapy.

Data Resources
RNA-sequencing data for the TCGA-colon adenocarcinoma (TCGA-COAD) and TCGA-rectal adenocarcinoma (TCGA-READ) cohorts, including 638 CRC samples and 51 normal tissue samples, were downloaded from public databases (https://portal.gdc.cancer.gov/). Corresponding clinicopathological characteristics for each patient-including age, gender, race, tumor location, disease type, tumor stage, tumor-node-metastasis (TNM) classification, venous invasion, lymphatic invasion, pretreatment CEA level, and survival information-were also retrieved from the TCGA data portal.
Only patients with both survival information and expression data were included in the present study. Another mRNA expression profile for 308 normal tissue samples was obtained in transcriptsper-million (TPM) format from the Genotype-Tissue Expression (GTEx) project (https://www.gtexportal.org/home/datasets), which is another large-scale repository cataloging gene expression from healthy individuals. Then, Ensembl gene IDs were mapped to human gene SYMBOL in terms of GENCODE V22 annotations for human datasets through the use of R/Bioconductor packages. As the raw mRNA sequence datasets from TCGA were normalized in terms of fragments per kilobase million (FPKM) via log 2 (FPKM+1), these datasets were scaled to a total depth of 10 6 fragments per sample and were interpreted as TPM in order to more easily compare the proportion of reads that was aligned to a given gene in each sample. Subsequently, any gene with a mean expression of ≤0.3 across all samples was deleted from the final mRNA expression matrices for subsequent analysis. Six independent datasets from the GEO database were used for external validation in the present study, including GSE41258, GSE44076, GSE37182, GSE20842, GSE83889, and GSE39582, together with survival information. A normalized expression matrix from GEO database was applied directly for the analyses. The protein expression levels of SEMA6B in clinical specimens from CRC patients and normal tissues were examined using immunohistochemical data from the HPA database (http://www.proteinatlas.org/). Since the data used in the present study were provided by TCGA and GEO, informed consent or ethical approval was not required. Furthermore, the present study fully adhered to all TCGA publication guidelines.

Survival Analysis
Kaplan-Meier (KM) curves were plotted to compare overall survival (OS), disease-free survival (DFS), progression-free survival (PFS), and relapse-free survival (RFS). These curves were generated with an optimum cut-off value for SEMA6B mRNA expression using the survfit function from the R package 'survminer', and a log-rank test was conducted to compare differences between survival status. Univariate and multivariate analyses of Cox proportional-hazards regression models were performed to obtain hazard ratios (HRs) with 95% confidence intervals (CIs) and statistical significance; the results were illustrated using a forest plot via GraphPad Prism 8.0. PFS-related nomogram models were established based on the multivariate Cox regression results. Calibration curves were drawn, and the concordance index (C-index) was computed to assess the prediction power of the nomogram.
Additionally, the prognostic values of SEMA6B expression in breast, esophageal, stomach, liver, lung, and ovarian cancers were assessed by the best cut-off values via Kaplan-Meier plotter (www.kmplot.com). HRs with 95% CIs and log-rank p values were also computed on the Kaplan-Meier plotter web page.

DNA Mutation and Methylation Analyses
To investigate the regulation of expression associated with the expression profile of SEMA6B, DNA mutation and methylation analyses were explored via online databases. Specifically, somatic mutation information was identified by the cBioPortal platform (www.cbioportal.org), which is a comprehensive web resource for exploring, visualizing, and analyzing multidimensional cancergenomic data. Methylation changes in SEMA6B in CRC and adjacent normal tissues were compared using UALCAN (http:// ualcan.path.uab.edu/index.html) (Chandrashekar et al., 2017) and Wanderer (http://maplab.imppc.org/wanderer/) (Díez-Villanueva et al., 2015) databases, which are web tools that can be employed to analyze DNA methylation profiles and gene expression from TCGA.

ROC Analysis
Receiver operating characteristic (ROC) analysis was used to evaluate the diagnostic accuracy for both OS and PFS; areas under the curve (AUCs) as well as p values were calculated via SPSS 25.0 software.

Identification of Differentially Expressed Genes
In accordance with the optimum cutoff value in KM survival analysis for OS, patients were classified into two groups (low and high SEMA6B expression) across TCGA datasets. Linear models were used to screen differentially expressed genes (DEGs) between these two groups by using the R package, "limma". The threshold for identifying DEGs was set as the false discovery rate (FDR)-adjusted p value <0.01 and absolute value of log2 (fold change) ≥ 1. With the R package, "gghplot2", a volcano plot was generated to visualize fold changes and t-test criteria.

GO and KEGG Pathway Enrichment Analyses of DEGs
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using 1,789 overexpressed DEGs were performed using R software with the aid of the "clusterProfiler" package. Biological processes (BP) and molecular functions (MF) were included in the GO enrichment analysis. Only terms with an FDR-adjusted p value <0.01 were considered to be statistically enriched. The top-15 enriched terms ordered by q value, from small to large, are shown in the corresponding plot.
Gene Set Enrichment Analysis of DEGs GSEA (version 4.1.0) was used to evaluate correlations between SEMA6B expression (high vs low) using the TCGA dataset. The annotated gene set was c2. cp.kegg.v6.2. symbols.gmt. Standard settings with 1,000 runs of gene permutations were employed for each analysis to determine the enriched pathways. Normalized enrichment scores (NES) and FDR-adjusted p values were obtained to indicate significantly enriched gene sets and pathways.

Gene Set Variation Analysis and Functional Annotation
To investigate the difference on biological pathways and processes according to the expression patterns of SEMA6B, GSVA was employed with the "GSVA" R package. GSVA is a non-parametric unsupervised method to explore the variation of pathway activity over samples (Hänzelmann et al., 2013). The annotated gene set was also "c2. cp.kegg.v6.2. symbols" downloaded from the MSigDB database (https://www.gseamsigdb.org/gsea/index.jsp). An adjusted p value ≤0.05 was considered statistically significance. A heatmap was drawn to display the enriched score value of each sample using the "pheatmap" R package.

Generation of Immune and Stromal Scores
The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm (Yoshihara and Shahmoradgoli 2013) was applied to calculate immune scores, stromal scores, ESTIMATE scores, and tumor purity for each patient from the TCGA dataset via R software loaded with the "estimate" package with default parameters.

TIMER Database Analysis
The Tumor Immune Estimation Resource (TIMER) is a webbased platform for systematic analysis of immune infiltrates across diverse cancer types from the TCGA (http://timer. comp-genomics.org/) (Li et al., 2017), which adopts a deconvolution of previously published computational approaches to infer TIICs from gene expression profiles. In the present study, correlations between SEMA6B expression and TIICs (B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells) were investigated separately in TCGA-COAD and TCGA-READ. Meanwhile, correlations between SEMA6B expression and gene markers of TIICs were also analyzed via the "Gene_Corr" module.

Tumor-Immune System Interactions Database Analysis
To further determine the relationship between SEMA6B mRNA expression and the abundance of TIICs, the Tumor-Immune System Interactions database (TISIDB) was used to determine the correlation between SEMA6B mRNA expression and tumorinfiltrating lymphocytes (TILs). TISIDB is a web portal for assessing tumor and immune system interactions, which integrates multiple heterogeneous data types (http://cis.hku.hk/ TISIDB/index.php) (Ru et al., 2019).

Gene Correlation Analysis
The Gene Expression Profiling Interactive Analysis (GEPIA2) online database (http://gepia2.cancer-pku.cn/#index) (Tang et al., 2019) was used to determine correlations between SEMA6B mRNA expression and gene markers of TIICs in CRC and adjacent normal tissues. Pearson correlation coefficients between SEMA6B expression and immunoinhibitor gene sets were visualized using R software with the "corrplot" package.

Heatmaps and Hierarchical Clustering Analysis
The "complexHeatmap" package from Bioconductor was used to plot heatmaps in terms of the expression of immunosuppressive gene sets as well as immune and stromal scores in different subgroup samples with R software. Specifically, Z-score normalization was applied in the expression dataset matrix, and then Euclidean distance was used to determine hierarchically clustered.

Analysis of Cell Viability
The mRNA expression level of SEMA6B in different colon cell lines was examined using the Cancer Cell Line Encyclopedia (CCLE) website (https://sites.broadinstitute.org/ccle). According to the survey, we selected two cell lines with high SEMA6B levels for subsequent research, i.e., HCT116 and LoVo cells.
A cell-counting kit 8 (CCK8) was used to determine cell proliferation. HCT116 and LoVo cells were plated and cultured in a 96-well plate with 1,500 cells per well, and then were interfered with by SEMA6B-and NC-siRNA for 0, 24, 48, or 72 h. After the interference, the supernatant was removed, and 100 μl DMEM or DMEM F12K was added in the presence of 10 μl of CCK8, and the cells were incubated at 37°C for 4 h. Then, the absorbance at 450 nm was assayed. potency of cells in the migration and invasion stages. For the invasion assay, 2 × 10 5 cells in serum-free medium were added to the upper chamber. Then, 500 μl of 10% FBS DMEM or DMEM F12K was added to the lower chamber, and the number of cells that had migrated after 48 h was quantified by counting five random fields under a microscope (IX70; Olympus Corp., Tokyo, Japan). Similar methods were adopted for the migration assay, except that 1 × 10 5 cells were added to the upper chamber without the Matrigel coating. Five random fields were counted for each chamber.

Statistical Analysis
Data analysis and visualization were performed using R software (version 4.0.0) with appropriate packages, as well as with SPSS25.0 (IBM Corp., Armonk, New York, United States) and GraphPad Prism 8.0 (GraphPad Software Inc., San Diego, CA, United States). Pearson's correlations were used to analyze pairwise relationships between mRNA levels of different genes. For continuous variables, differences among groups were analyzed using one-way analysis of variance (ANOVA) for comparisons of more than two groups, and via t-tests for comparisons between only two groups. Chi-squared tests (χ 2 ) were used to evaluate associations between SEMA6B expression and clinicopathological parameters. A two-sided p < 0.05 was considered statistically significant.

SEMA6B Expression is Upregulated in CRC Tissues
Using RNA-sequencing data from the TCGA-CRC dataset and GTEx project, we found that SEMA6B mRNA expression in CRC tissues (n 638) was significantly higher than that in normal colorectal tissues (n 359; p < 0.001). Meanwhile, upregulated SEMA6B mRNA levels in CRC tissues were also validated in GSE datasets, including GSE41258 (p < 0.001), GSE44076 (p < 0.001), GSE37182 (p < 0.001), GSE20842 (p 0.003), and GSE83889 (p 0.015). ( Figure 1A). Moreover, immunohistochemical staining obtained from the HPA database demonstrated that protein expression levels of SEMA6B were consistent with their transcriptional levels in comparison with those in normal tissues ( Figure 1B). We further detected the function of methylation in regulating the expression of SEMA6B and found that the DNA methylation levels of SEMA6B were dramatically downregulated in CRC tissues compared with those in normal colorectal samples (p < 0.001) using the UALCAN web tool (Supplementary Figure S1A). As shown in Supplementary Figure S1B, data from the Wanderer web tool were similar (p < 0.05; normal n 44, tumor n 400), with most of the SEMA6B probes in the 450 methylation array exhibiting significant differences between CRC and normal specimens. Otherwise, few mutational or somatic copy-number alterations of SEMA6B were observed in CRC tissues (Supplementary Figure S2).

High SEMA6B Expression Predicts Poor Prognosis of CRC Patients
We next assessed associations of SEMA6B expression with clinicopathological features of CRC patients using the TCGA-CRC dataset. As shown in Figure 2, the expression levels of SEMA6B were significantly associated with venous invasion, T stage, MSI, KRAS mutation, and CMS. However, SEMA6B expression in CRC was not significantly correlated with other clinicopathological characteristics, including age, gender, tumor site, lymphatic invasion, M stage, N stage, pathological stage, and patient status.
To evaluate the prognostic value of SEMA6B, CRC patients were divided into high and low SEMA6B expression groups according to the optimal cut-off value of SEMA6B levels. Kaplan-Meier survival curve analysis showed that CRC patients with high SEMA6B expression had shorter OS (p 0.040), 5 year survival (p 0.042), DFS (p 0.002), and PFS (p < 0.001) compared to those with low SEMA6B expression ( Figure 3A). To validate these findings, we analyzed the prognostic significance of SEMA6B using another GEO cohort (GSE39582). In line with results in the   DFS for GSE17538 dataset. Forest plots visualizing the hazard ratios (HRs) and 95% confidence intervals (CI) for each variable. Differences with p < 0.05 were considered significant.
Frontiers in Molecular Biosciences | www.frontiersin.org December 2021 | Volume 8 | Article 687319 7 TCGA dataset, the high SEMA6B expression group exhibited poor survival (OS: p 0.022; RFS: p 0.050) ( Figure 3B). Subsequently, the Kaplan-Meier plotter database was used to analyze the prognostic potential of SEMA6B in different cancers. Interestingly, high SEMA6B expression levels were associated with poor prognosis of OS in lung adenocarcinoma and lung squamous cell carcinoma (Supplementary Figure S3), while SEMA6B expression had less influence on the prognoses within other types of cancers (Supplementary Figure S3).
Furthermore, to explore the clinical prognostic significance of SEMA6B in CRC, Cox regression analysis was performed to determine PFS and OS. Univariate Cox regression analysis showed that SEMA6B, venous invasion, N stage, pretreatment CEA, T stage, and M stage were significantly associated with PFS in CRC patients sourced from the TCGA dataset (p < 0.05; Figure 4A). Meanwhile, multivariate Cox regression analysis showed that SEMA6B, T stage, and M stage were independent prognostic factors for PFS among CRC patients (p < 0.05; Figure 4A). For OS, univariate Cox regression analysis showed that SEMA6B, age, venous invasion, pretreatment CEA, T stage, N stage, and M stage have statistical significance; however, SEMA6B expression was not an independent prognostic factor in the multivariate analysis for OS (p > 0.05; Supplementary Figure S4). Moreover, in the GSE17538 validation CRC dataset, SEMA6B was an independent prognostic factor for DFS in multivariate Cox regression analysis (p < 0.05; Figure 4B).

Prognostic Nomogram Models Based on SEMA6B for CRC
According to the results of multivariate Cox regression analysis, the independent prognostic signature of SEMA6B, M stage, and N stage was enrolled to establish a nomogram model to predict the 3-and 5 year PFS probabilities of each patient for clinical practice ( Figure 5A). The C-index of this model was 0.72 (95% CI, 0.70-0.75). Moreover, the calibration curve revealed that the 3-and 5 year PFS rates predicted by the nomogram were highly consistent with the actual observation outcomes ( Figure 5B). These results demonstrated that the nomogram models show a favorable prognostic ability for predicting PFS. In addition, ROC analysis demonstrated that the AUC values for OS or PFS of the prediction model-including pathological M stage, N stage, T stage, and SEMA6B expression-were significantly improved from 0.639 to 0.759 for OS and from 0.641 to 0.719 for PFS ( Figures 5C,D), indicating the additive predictive value of SEMA6B compared with that of known prognostic factors.

SEMA6B-Related Biological Processes
According to the aforementioned results, SEMA6B may play an essential role in the biological functions of CRC. To clarify the The calibration curve of the nomogram when predicting 3-and 5 years PFS probability. ROC curves for predicting OS (C) and PFS (D) in CRC patients via prognostic models from the TCGA dataset. The x-axis indicates the false-positive rate, which is presented as '1-Specificity'. The y-axis indicates the true-positive rate, which is designated as "Sensitivity".
Frontiers in Molecular Biosciences | www.frontiersin.org December 2021 | Volume 8 | Article 687319 biological roles of SEMA6B expression in CRC, DEGs were detected between high and low SEMA6B expression groups. Volcano plot analysis identified 1,794 DEGs ( Figure 6A). Among them, 1,789 genes were upregulated and 5 genes were downregulated ( Figure 6A). Then, the biological functions of these DEGs were explored by KEGG signaling pathway and GO annotation analysis. The top-15 pathways from the KEGG results showed that the most significantly altered pathways in the SEMA6B high expression group were those involving cytokine-cytokine receptor interactions, chemokine signaling pathways, rheumatoid arthritis, viral protein interactions with cytokines and cytokine receptors, and complement and coagulation cascades ( Figure 6B). Furthermore, GSEA was conducted to determine the potential mechanism underlying SEMA6B in CRC. According to NES, we selected the most highly enriched signaling pathways. As shown in Figure 6C, the genes in the SEMA6B high expression group were mainly enriched in inflammatory activities and immune-related processes such as leukocyte transendothelial migration, chemokine signaling pathways, cytokine receptor interactions, and complement cascades. Furthermore, tumor-aggressiveness feature-related gene sets including those involving the Janus kinase-signal transducer and activator of transcription (JAK-STAT) signaling pathway, cancer pathways, and vascular endothelial growth factor (VEGF) signaling pathway were also significantly enriched in the high SEMA6B expression group. Meanwhile, the results of GO analysis revealed that many biological functions of these DEGs were primarily associated with inflammatory responses and adaptive immune responses from categories of biological processes and molecular functions ( Figure 7A), respectively. To determine the biological behaviors based on distinct SEMA6B profiles in CRC patients, GSVA was performed. As expected, we found that the SEMA6B high-risk group was markedly enriched in tumor cell proliferation, immune response, EMT-related pathways, such as the JAK-STAT signaling pathway, antigen processing and presentation, natural killer cell mediated cytotoxicity, the chemokine signaling pathway, cytokinecytokine receptor interaction, ECM receptor interaction, and focal adhesion ( Figure 7B).  Supplementary Table S1. As shown in Figure 8A, SEMA6B expression was positively correlated with stromal   Given that TILs are a prognostic indicator for CRC (Reissfelder et al., 2015;Berntsson et al., 2017), associations between SEMA6B gene expression and TILs infiltration levels in CRC were analyzed via the TIMER database. The results revealed that SEMA6B was negatively correlated with tumor purity in COAD (r −0.39, p 2.86E-16) and READ (r −0.39, p 1.94E-06), whereas it was strongly positively correlated with infiltrating levels of CD4 + T cells (r 0.504, p 2.65E-27), macrophages (r 0.465, p 4.53E-23), neutrophils (r 0.542, p 5.37E-32), and dendritic cells (r 0.545, p 1.90E-32) in COAD ( Figure 8B). Similar trends in terms of correlations between SEMA6B and TILs infiltrating levels were also observed in READ ( Figure 8B). Additionally, we also investigated the relationship between the expression levels of SEMA6B and immune infiltration in TISIDB database so as to further verify the role of SEMA6B in TME. Spearman's correlation analysis illustrated that SEMA6B was strongly related to immune infiltration in COAD (Supplementary Figure S5A) and READ (Supplementary Figure S5B To further determine correlations between SEMA6B and TILs, we analyzed relationships between SEMA6B and marker genes of different immune cells in COAD and READ via the TIMER and GEPIA databases. After correlation adjustments by purity, the findings demonstrated that expression levels of most marker sets of CD4 + T cells, Tregs, exhausted T cells, monocytes, tumor-associated macrophages (TAMs), M2 macrophages, and dendritic cells presented strong correlations with SEMA6B expression in COAD and READ (Table 1). In the GEPIA database, the analyses revealed that SEMA6B was strongly related to TILs infiltration in CRC tissues compared to that in normal colorectal samples, especially in terms of important signatures of pro-cancer immune cells, such as chemokine (C-C motif) ligand (CCL)-2, CD68, interleukin 10 (IL10) of TAMs, CD163, V-set and immunoglobulin domaincontaining (VSIG4), membrane-spanning 4-domain subfamily A, member 4A (MS4A4A) of M2 macrophages, forkhead box protein P3 (FOXP3), transforming growth factor beta (TGFβ), C-C chemokine receptor 8 (CCR8), signal transducer and activator of transcription 5B (STAT5B) of Tregs, PD-1, CTLA4, and T-cell immunoglobulin mucin 3 (TIM-3) of T cell exhaustion ( Table 2).

SEMA6B Overexpression is Indicative of the Tumor Immunosuppressive Microenvironment
To further determine the influence of SEMA6B on the tumor microenvironment, we investigated correlations between SEMA6B and expression of genes negatively regulating the cancer-immunity cycle, which consists of a cycle of processes involving eradication of cancer by the immune system (Chen and  Mellman 2013). Gene signatures were downloaded from the Tracking Tumor Immunophenotype website (http://biocc. hrbmu.edu.cn/TIP/index.jsp) (Xu et al., 2018). Unsupervised hierarchical clustering analyses revealed that genes involved in the negative regulation of the cancer-immunity cycle were mostly upregulated in the high SEMA6B expression group ( Figure 9A). A correlation coefficient heatmap showed that SEMA6B had a significant positive correlation with immunosuppressive molecules, including colony stimulating factor 1 receptor (CSF1R), PD-L2, FOXP3, CD86, and TGFβ1 ( Figure 9B). As expected, the expression levels of immune checkpoints-such as PD-L1, PD-1, CTLA-4, T cell immunoreceptor with Ig and ITIM domains (TIGIT), TIM-3, and lymphocyte activation gene 3 (LAG-3)-in the high SEMA6B expression group were significantly higher than those in the low SEMA6B expression group ( Figure 9C). These results suggest that SEMA6B plays a vital role in immune escape via promoting the tumor immunosuppressive microenvironment.  The mRNA level of SEMA6B in different colon cancer cell lines was investigated through CCLE datasets. Bar plots show that HCT116, SW480, LoVo, SW620, GP2D, and LS513 are the top six colon cancer cell lines with the highest SEMA6B expression levels ( Figure 10A). Therefore, HCT116 and LoVo cell lines were chosen for further study. QRT-PCR analysis revealed that SEMA6B silencing can not only significantly reduce the expression of SEMA6B, but also decrease the mRNA levels of immunosuppressive molecules, EDNRB, TGFB2, IL1B, IL6ST, BTLA, PD-L1, LGALS1, ICAM1, HAVCR2 in the selected cells lines compared to the corresponding controls ( Figure 10B). Knockdown of SEMA6B significantly reduced the proliferation in both cell lines according to the results of the CCK-8 assay ( Figure 10C). Furthermore, the cell numbers that migrated through the membrane were significantly reduced in the SEMA6Bsilenced groups of both cell lines according to the results of the Transwell migration assay ( Figure 10D). Transwell invasion assay also showed that SEMA6B silencing significantly decreased the degree of invasiveness in both selected colon cancer cell lines ( Figure 10E). These data demonstrate that SEMA6B knockdown reduces the growth and progression of colon cancer cells as well as suppresses the formation of an immunosuppressive microenvironment.

DICUSSION
Accumulating evidence has indicated that the intra-tumoral immune contexture (i.e., type, functional orientation, density, and location of immune cells) of solid tumors may be a dominant determinant of clinical outcomes (Becht et al., 2016;Fridman et al., 2017). Despite CRC being one of the major cancer types for which new immune-based cancer treatments are currently in development, the prognosis of patients with advanced CRC remains poor (Emambux and Tachon 2018;Ciardiello and Vitiello 2019). Therefore, it has been necessary to further identify immune-related biomarkers and better elucidate the underlying molecular mechanisms of CRC to improve patient prognosis and guide the development of immunomodulation for CRC treatments. In the present study, we focused on SEMA6B, a member of the semaphorin axon-guidance family, which exhibits immune functions related to the control of cellular movements and cell-cell communication (Tawarayama et al., 2010;Koivisto et al., 2020). Given that our knowledge of the role of SEMA6B in cancers is limited, in the present study, we aimed at investigating its prognostic value in CRC and revealed its associated biological functions by performing a comprehensive analysis of population databases. Here, we found that the mRNA expression levels of SEMA6B were significantly increased in CRC compared with those in normal colorectal tissues using the TCGA and GEO databases. A similar trend in the protein expression levels of SEMA6B was also observed in the HPA database. In addition, our data demonstrated that gene hypomethylation was influential in upregulation of SEMA6B expression in CRC. These findings imply that enhanced SEMA6B expression had considerable effects on CRC progression in a manner that may be due to epigenetic alterations. Of note, SEMA6B levels were related with the disease type, venous invasion, T stage, MSI, KRAS mutation, and CMS of CRC patients. Moreover, high SEMA6B expression was found to predict worse survival in all cohorts, and was further shown to be an independent prognostic factor of PFS in CRC patients. The nomogram based on SEMA6B, M stage, and N stage could facilitate accurate prediction of the 3-and 5 year PFS probabilities in CRC patients. To our knowledge, this is the first study to report a consistent association between increased SEMA6B levels and poor prognosis in CRC. These results indicate a malignant biological influence of high SEMA6B levels in CRC.
Studies based on the roles of SEMA6B have also been reported for other malignant human tumors. In breast cancer, high levels of SEMA6B in human MCF-7 cells exhibit an in vitro invasive capacity, and show potential as a key regulator of tumor progression (Murad and Collet, 2006). Recently, SEMA6B has been shown to exhibit higher expression levels in testicular cancer tissues than in normal tissues and is considered to be a predictor of poor prognosis in patients with testicular germ-cell tumors (Ji and Wang, 2020). On the contrary, SEMA6B has been reported to be downregulated in human glioblastoma cells upon prolonged treatment with the anti-tumor action of all-trans retinoic acid (ATRA) (Correa et al., 2001). Increasing evidence has indicated that SEMA6B is related to macrophages and correlates with a favorable prognosis in glioma patients (Sun et al., 2019). Thus, these conflicting findings suggest that SEMA6B may play differential roles in different kinds of human cancers, and that discrepancies between SEMA6B expression and prognostic values may derive from underlying mechanisms pertinent to specific biological properties in various cancers. Exploring and revealing the mechanisms of SEMA6B in CRC may facilitate the discovery of novel therapeutic approaches for CRC patients.
Biological pathway analysis and functional enrichment analysis in our present study illustrated that immune-related processes, inflammatory activities, and cancer signaling pathways were significantly enriched in the high SEMA6B expression group. Previous studies have indicated that immune infiltrating levels in tumor sites influence prognosis and the response rate of chemotherapy and immunotherapy in CRC patients (Waniczek et al., 2017). The ESTIMATE algorithm was first reported by Yoshihara et al. to assess immune-cell infiltration and the presence of stromal cells based on gene expression data (Yoshihara and Shahmoradgoli, 2013). In the present study, we revealed that high SEMA6B expression was positively associated with higher stromal scores, immune scores, and ESTIMATE scores but was negatively association with tumor purity. Another noteworthy finding of the present study was that SEMA6B expression was correlated with diverse immune infiltration levels in CRC. Our results showed that there were moderate-to-strong positive relationships between SEMA6B expression levels and infiltration levels of macrophages, MDSCs, NK cells, Tregs, and Th1 cells, as well as significantly positive correlations between infiltrating levels of CD4 + T cells, neutrophils, and dendric cells and SEMA6B expression in CRC. Likewise, the relationships between gene markers of different immune cells and SEMA6B expression are suggestive of a pivotal role of SEMA6B in regulating the tumor immune microenvironment. In addition to identifying markers of CD8 + T-cell activation, we also found that NK cells were more active in tumors with high SEMA6B expression, which indicates that these tumors may have an antitumor immune microenvironment (Qu et al., 2018). However, the above-mentioned immune cells do not serve as the key effectors that destroy tumor cells in CRC patients. This phenomenon may be due to the following possibilities. First, immunosuppressive cells such as Tregs, M2 macrophages, and MDSCs (Qu et al., 2018;Xiang and Gilkes 2019) are known to dominate the immune microenvironment in tumors with high SEMA6B expression. Second, SEMA6B has the potential to induce CD8 + T-cell exhaustion, as reflected by positive correlations with T-cell exhaustion markers. TIM-3, a crucial surface protein on exhausted T cells (Anderson 2012), was highly correlated with SEMA6B expression in CRC in the present study. Third, SEMA6B expression represents weak correlations with gene markers of M1 macrophages, whereas M2 macrophage indicators exhibited strong correlations in the present study. Previous studies have underscored the dualistic roles of macrophages in tumors: the M1 phenotype exhibits proinflammatory and tumoricidal activities, whereas the M2 phenotype exhibits anti-inflammatory activities and exerts functions in immunosuppression and malignant progression of tumors (Biswas and Mantovani 2010;Mantovani and Sica 2010). Our present results suggest a potential regulatory role of SEMA6B in polarization of TAMs with the M2 phenotype for evading immune surveillance. The cancer-immunity cycle is a series of functional stepwise events required to obtain an efficient control of cancer growth by the immune system (Chen and Mellman 2013;Horton et al., 2019). This process consists of the following seven steps (Chen and Mellman 2013): 1) releasing of cancer cell antigens; 2) cancer-antigen presentation; 3) priming and activation; 4) trafficking of immune cells to tumors; 5) infiltration of immune cells into tumors; 6) recognition of cancer cells by T cells; and 7) killing of cancer cells. Immune checkpoint molecules-such as PD-1, PD-L1, and CTLA4-can inhibit the development of an active immune response by acting primarily at the level of T-cell development and proliferation (step 3) (Chen and Mellman 2013;Chen and Mellman 2017). Negative gene sets in the cancer-immunity cycle-including PD-1, PD-L1, LAG-3, TIM-3, and TIGIT-can have an inhibitory function that primarily acts to modulate active immune responses in the tumor bed (step 7) (Chen and Mellman 2013;Chen and Mellman 2017). Therefore, tracking tumor immunophenotypes may be essential for further elucidating the mechanisms of cancer immunity and progressing the development of biomarkers of responses to immunotherapy. Our present study demonstrated that in CRC patients with high SEMA6B expression, the genes involved in the prevention of T-cell priming and the induction of tolerance were increased, and immune checkpoints-such as PD1, PD-L1, CTLA-4, LAG3, TIM-3, and TIGIT-were also upregulated in this group. Hence, accumulating evidence suggests that SEMA6B may contribute to tumor development by both attenuation of tumor-specific cytotoxic T-cell responses and induction of an immunosuppressive state (Kigel and Rabinowicz 2011;Ge and Li 2013;Cui and Bian 2021).
We conducted SEMA6B knockdown experiments to further confirm the biological role of SEMA6B in two colon cancer cell lines, i.e., HCT116 and LoVo. The results demonstrate that SEMA6B silencing enables a reduction in proliferation, migration and invasion in vitro; meanwhile, the mRNA expression levels of immunosuppressive molecules were also diminished in silenced colon cancer cells by qRT-PCR. Gui et al. reported that siRNA-mediated knockdown of SEMA6B weakened gallbladder cancer cell invasion and migration (Cui and Bian 2021). Ge et al. also discovered that SEMA6B may promote gastric cancer invasion and metastasis and represents a reliable biomarker for the clinical diagnosis and therapy of gastric cancer (Ge and Li 2013). Overexpression of SEMA6B in BHK-21 cells could promote cell proliferation, and the inhibition of SEMA6B signaling suppresses tumor formation for the sake of abrogation of bFGF and VEGF signaling (Kigel and Rabinowicz 2011). Overall, our findings reveal that SEMA6B may play a key role in regulating CRC progression and helps shape the immunosuppressive microenvironment. However, further studies are required to validate the biological functions of SEMA6B in CRC.
There were some limitations of our present study. First, only transcriptomic expression of SEMA6B expression with clinical data was analyzed to predict prognosis in CRC from public databases; thus, our results should be validated in larger sample size. Second, this was a retrospective study with selection biases inherent in the cohorts; thus a prospective study is also needed. Third, the molecular mechanisms of SEMA6B in CCR remain unclear, despite a series of functional annotations and enrichment analysis being investigated. Thus, further study is required to recognize the potential biological mechanisms of SEMA6B using additional experimental approaches.
In conclusion, our present study explored the clinical value and biological processes of SEMA6B, using CRC samples from the TCGA and GEO databases on a large scale. Our data revealed Frontiers in Molecular Biosciences | www.frontiersin.org December 2021 | Volume 8 | Article 687319 that high SEMA6B expression was significantly correlated with cancer progression, poor survival, and immune infiltration in patients with CRC. Moreover, high SEMA6B expression was correlated with increased immunosuppressive cellular infiltration and the expression of immune checkpoints. Moreover, in vitro cell studies validate that overexpress SEMA6B may promote proliferation and metastasis in two colon cancer cell lines, and help to foster an immunosuppressive microenvironment. Our present findings therefore offer novel insights into how SEMA6B affects prognosis and the immune microenvironment in CRC and suggests that SEMA6B may serve as a potential biomarker for developing immunotherapeutic strategies for assessing the efficacy and responsiveness of CRC treatments.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.