CYP1B1: A Novel Molecular Biomarker Predicts Molecular Subtype, Tumor Microenvironment, and Immune Response in 33 Cancers

Simple Summary Cytochrome P450 Family 1 Subfamily B Member 1 (CYP1B1) is a critical metabolic enzyme of melatonin. Although melatonin has been identified to exhibit tumor suppressing activity, the role and mechanism of the clinical and immunological characteristics of CYP1B1 in cancer remain unclear. We comprehensively explored the clinical and immunological characteristics of CYP1B1. We identified that the dysregulated expression of CYP1B1 was associated with clinical characteristics and a tumor immune microenvironment, which may provide a promising predictor and molecular target for clinical immune treatment. Abstract Background: Cytochrome P450 Family 1 Subfamily B Member 1 (CYP1B1) is a critical metabolic enzyme of melatonin. Although melatonin has been identified to exhibit tumor suppressing activity, the role and mechanism of the clinical and immunological characteristics of CYP1B1 in cancer remain unclear. Methods: In this study, RNA expression and clinical data were obtained from The Cancer Genome Atlas (TCGA) across 33 solid tumors. The expression, survival, immune subtype, molecular subtype, tumor mutation burden (TMB), microsatellite instability (MSI), biological pathways, and function in vitro and vivo were evaluated. The predictive value of CYP1B1 in immune cohorts was further explored. Results: We found the dysregulated expression of CYP1B1 was associated with the clinical stage and tumor grade. Immunological correlation analysis showed CYP1B1 was positively correlated with the infiltration of lymphocyte, immunomodulator, chemokine, receptor, and cancer-associated fibroblasts (CAFs) in most cancer. Meanwhile, CYP1B1 was involved in immune subtype and molecular subtype, and was connected with TMB, MSI, neoantigen, the activation of multiple melatonergic and immune-related pathways, and therapeutic resistance. Conclusions: Together, this study comprehensively revealed the role and mechanism of CYP1B1 and explored the significant association between CYP1B1 expression and immune activity. These findings provide a promising predictor and molecular target for clinical immune treatment.


Introduction
Melatonin, an endogenous hormone, is secreted by the pineal gland in response to biological rhythm [1,2]. The synthesis and metabolism of the hormone involves a series of biological pathways. In the process of melatonin synthesis, N-acetylserotonin

Data Collection
The expression, phenotype, and survival data were obtained from the UCSC Xena database (https://xenabrowser.net/, accessed on 22 May 2020). The lymphocyte, immunomodulator, chemokine, immune subtype, and molecular subtype were downloaded from the Tumor-Immune System Interactions (TISIDB) database (http://cis.hku.hk/TISIDB/index.php, accessed on 23 March 2019). The cancer-associated fibroblasts (CAFs) data were obtained from TIMER 2.0 (http://timer.cistrome.org/, accessed on 2 July 2020). The TCGA database was used to obtain the tumor mutation burden (TMB), microsatellite instability (MSI), and neoantigen. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013) and was approved by the Ethics Committee of The Sixth People's Hospital of Huizhou City (PJ2022MI-KJ038).

Differential Gene Expression Analysis
To identify the expression differences of CYP1B1 between tumor samples and normal tissues, the expression data of 33 cancers was downloaded from the UCSC Xena database.
The expression values were normalized by Transcripts Per Million (TPM) transformation. A distinction with a threshold of p < 0.05 was considered as having a significance.

Survival Analysis
The KM-plotter analysis of 33 cancer patients were examined using univariate COX regression analysis to determine the prognostic significance of CYP1B1. The forest plot was performed using the R software forest plot package. A log-rank test with a threshold of p < 0.05 was considered as having a significance.

Correlation Analysis between CYP1B1 and Tumor Immune System
To clarify the relation between CYP1B1 and the tumor immune system, the potential relationship among CYP1B1 and 28 lymphocytes, 69 immunomodulators (45 immunostimulators, 24 immunoinhibitors), 41 chemokines, and 18 receptors, immune subtypes and molecular subtypes were explored using the TISIDB database. The relationship of CYP1B1 and CAF was assessed by the TIMER 2.0 database. The correlation of CYP1B1 and TMB, MSI, and neoantigen were further evaluated. p < 0.05 was considered as having a significance.

Single-Sample Gene Set Enrichment Analysis (ssGSEA)
To identify the CYP1B1 activity in cancer, the single sample gene sets enrichment analysis (ssGSEA) between high CYP1B1 and low CYP1B1 expression was evaluated by using R software GSVA package. p < 0.05 was considered as having a significance.

Association between CYP1B1 and Drug Response
CCLE gene expression data were quantile normalized among all different cell lines for partial correlation, and then Z-score normalization was applied in each tissue to calculate the expression difference between High-Low (using the median as a cutoff) IC50 groups. The X-axis indicates the mean/median expression difference across tissues. Correlations of CYP1B1-drug associations after controlling for the tissue average expression were analyzed.

Gene Set Enrichment Analysis (GSEA)
Correlations with other genes and ordered genes based on findings were performed to discover the biological aspects of CYP1B1. The sorted gene list was used in GSEA analysis to see if highly linked genes clustered in genuinely functional pathways. The reference gene set was annotated gene set c5.go.v7.4.symbols.gmt and c2.cp.kegg.v7.4.symbols.gmt. As previously stated, FDRs of 0.05 and p-values of 0.01 were considered significant.

Cell Culture
The lung cancer cell line PC9, breast cancer cell lines MDA-MB-231 was obtained from the Sun Yat-sen University Cancer Center. Cells were maintained at 37 • , 5% CO 2 in 10% DMEM (Invitrogen, Carlsbad, CA, USA) or RPMI-1640 (Thermo Fisher Scientific, Waltham, MA, USA) supplemented with 10% fetal bovine serum.

Wound Healing Assay
Cells were cultured in 10% DMEM (Invitrogen, Carlsbad, CA, USA) or RPMI-1640 (Gibco, Australia), supplemented with 10% fetal bovine serum at 37 • , 5% CO 2 . The medium was removed, and the surface of the inoculated cells was scratched and marked with a 10 µL pipette tip. A sterile 200-L pipette tip was used to scrape the surface of the cell monolayer to create the wound. Under an inverted microscope, photographs were taken at time 0 and 24 h after cell scratching (Olympus Corporation). The area of each wound was calculated using Image J software (National Institutes of Health).

Transwell Assay
Transwell plates (8-m pores) were used for Transwell migration or invasion studies. Next, 5 × 10 4 (migration assay) or 1 × 10 5 (invasion assay) cells resuspended in serum-free medium were placed in the top chamber, either uncoated or covered with Matrigel (BD Biosciences). The culture medium in the lower compartment included 10% FBS. The cells were fixed and stained after being incubated for 12 or 24 h. Cells on the undersides of the filters were seen and counted at a magnification of 200 magnification.

Cell Proliferation Assay
The cancer cell line was seeded in 1000 cells per plate and cultured for two weeks at 37 • C in a 5% CO 2 incubator. The cells were washed twice, fixed for 15 min with 4% paraformaldehyde, then stained for 20 min at room temperature with 1% crystal violet solution. The number of visible colonies was determined.

In Vivo Metastasis Model
The Guangzhou Medical University Cancer Center's Institutional Animal Care and Use Committee authorized the animal operations (Guangzhou, China, G2022-050). First, 1 × 10 7 PC9 cells were injected into the footpads of mice for the tumor metastasis model. After 2 weeks, the mice were euthanized, and their footpad tumors and inguinal lymph nodes were removed after being intraperitoneally administered with melatonin (25 mg/kg) every other day for 7 consecutive days.

Statistical Analysis
All statistical analyses were performed using R software (version 4.0.3). The correlations between CYP1B1 and clinicopathological features were detected by Chi square test. Univariate analysis was used to estimate the prognostic value of CYP1B1. The correlation analysis was assessed by the Spearman rank test. p < 0.05 was considered statistically significant.

The Expression, Tumor Stage and Clinical Grade of CYP1B1 in Cancer
By investigating RNA-seq data from 33 cancers in the TCGA, we explored the differential expression of CYP1B1 between tumor samples and healthy samples. As shown in Figure 1A, CYP1B1 was significantly downregulated in most of the cancers.

Correlation between CYP1B1 and CAFs in the Tumor Microenvironment
CAFs, the most abundant stromal cells in the tumor microenvironment (TME), are associated with tumor cell growth, invasion, and metastasis, metabolic reprograming, immune escape, and therapeutic resistance [29][30][31]. Four algorithms, including EPIC, MCPCOUNTER, XCELL, and TIDE, were used to evaluate the correlation between the CAFs and the CYP1B1 expression level in 33 cancers. Cancers with consistent correlations of four algorithms were considered to be importantly associated with CAFs infiltration. As shown in Figure 3A Figure 3F-I). These results indicated that CAFs infiltration mediated by CYP1B1 were critical for cancer occurrence and progression in the TME.

Correlation between CYP1B1 and TMB, MSI and Neoantigen
To understand the role of CYP1B1 in immunotherapy, the association between CYP1B1 and immunotherapy-related biomarkers (TMB, MSI and neoantigen) was further assessed. As shown in Figure 4, CYP1B1 expression was positively associated with the TMB in COAD (p = 0.0044; Figure 4A) and LIHC (p = 0.00011; Figure 4A), and a negative association was found in STAD (p = 0.038; Figure 4A). CYP1B1 expression positively correlated significantly with MSI in OV (p = 0.0019) and LIHC (p = 0.0038) in Figure 4B, while a negative association in KIRP (p = 0.03), SARC (p = 8.7 × 10 −5 ), STAD (p = 0.0074), SKCM (p = 0.05), CHOL (p = 0.0032) and HNSC (p = 0.029) was identified in Figure 4B. Similarly, CYP1B1 expression positively correlated significantly with neoantigen in LIHC (p = 0.0012) and READ (p = 0.046) in Figure 4C, while a negative association in STAD (p = 0.037) was revealed in Figure 4C. The results indicated that immunotherapy-related biomarkers mediated by CYP1B1 may play important roles in immune pathways. As the CYP1B1 was one of metabolic enzymes in melatonin, we further explored the melatonin-related pathways, including melatonin biosynthesis, circadian clock, entrainment of the circadian clock, entrainment of the circadian clock by photoperiod, and BMAL1_clock_NPAS2 circadian gene expression. As shown in Figure 4D and E, a high expression of CYP1B1 was associated with the inactivation of melatonin biosynthesis in LIHC (p < 0.0001) and STAD, associated with the activation of the circadian clock, entrainment of the circadian clock, entrainment of the circadian clock by photoperiod, and BMAL1_clock_NPAS2 circadian gene expression ( Figure 4D,E). These results indicated that melatonin metabolism mediated by CYP1B1 may be involved in the complexity and heterogeneity of TME in cancer.
To further uncover the role of CYP1B1 in TME, 10 classical immune-related pathways, including immunological synapse, innate immune response, immunoglobulin complex, immunoglobulin binding, type 2 response, humoral immune response, immune effector process, B cell mediated immunity, adaptive immune response, and T cell mediated immunity, were enriched by ssGSEA analysis. Compared with low expression of CYP1B1, a high expression of CYP1B1 was significantly associated with the activation of immune-related pathways in LIHC ( Figure 4F) and STAD ( Figure 4G).  . Two groups (High-expression and Low-expression) of pathway scores with boxplots using the Mann-Whitney U test. * p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001.

The Association of CYP1B1 Expression and Therapeutic Response
To identification the function of CYP1B1 in breast cancer, as shown in Figure 5A, GSEA analysis indicated that CYP1B1 was significantly associated with epithelial mesenchymal transition (NES = 2.625, FDR = 7.6 × 10 −10 ), angiogenesis (NES = 2.169, FDR = 9.2 × 10 −6 ), and regulation of the immune response (NES = 2.713, FDR = 3.1 × 10 −9 ) and circadian rhythm (NES = 0.822, FDR = 8.5 × 10 −1 ). Wound healing and transwell migration and invasion assays revealed that the overexpression CYP1B1 increased breast cancer migration and invasion ( Figure 5B and C). Colony assays identified that CYP1B1 also improved breast cancer proliferation ( Figure 5D). To explore the function of CYP1B1 in vivo, a tumor metastasis model was established and the immunohistochemical staining analysis of serial tissue sections showed that the expression of angiogenesis marker CD31, proliferation marker ki67, metastasis marker MMP9, and neutrophil infiltration marker LY6G mediated by CYP1B1 could be reversed by melatonin ( Figure 6).

Discussion
Melatonin has been considered a promising anti-cancer drug in the treatment of cancer [32][33][34]. However, the molecular and clinical characteristics of the melatonergic metabolic enzyme CYP1B1 in cancer still remain unknown. In this study, we comprehensively investigated the clinical and immunological pattern of CYP1B1 determined from RNA-seq data across TCGA pan-cancer. The results indicated that the dysregulated expression of CYP1B1 was correlated with clinical and immunological characteristics in cancer and could be a promising predictor and molecular target for clinical immune treatment.
CYP1B1, one member of the cytochrome P450 family, has been identified in tumorigenesis [35,36]. The dysregulated expression of CYP1B1 was explored between tumor and normal tissues. As a result of the abnormal expression in cancer, CYP1B1 has been defined as a candidate tumor antigen [37]. Interestingly, CYP1B1 has also been exploited as a molecular target for immunotherapy, owing to the restricted expression profile in normal tissues [38]. Differential CYP1B1 expression was also observed in different clinical stages and tumor grades, and a high expression of CYP1B1 was also observed in renal cell carcinoma, which was related to advanced grades and late stages [39]. In vitro and vivo function experiments identified that CYP1B1 could increase tumor progression and metastasis. The results indicated that CYP1B1 was correlated with clinical characteristics in cancer and could be a prognostic predictor in cancer.
In the present study, the infiltration of lymphocyte, immune regulators, CAFs, immune subtype, and molecular subtype was mediated by CYP1B1. By using genetic study of the melatonergic microenvironment across 14 solid tumors, Lv et al. also identified that the melatonin catabolic enzymes, including CYP1B1, were associated with TMB and prognosis [40]. In vitro experiments found that the CYP1B1 expression was significantly correlated with B7-H3 expression in colorectal cancer. Nevertheless, in vivo experiments revealed that HLA-A*0201 could be processed and presented by CYP1B1-specific cytotoxic T lymphocytes (CTLs) [41,42]. The neutrophil infiltration mediated by CYP1B1 in the TME could be reversed by melatonin in vivo. Based on these findings, we speculated that CYP1B1 may regulate the immune cell infiltration in the tumor microenvironment, which could act as a molecular marker for clinical therapy.
TMB, MSI, and neoantigen have been identified as biomarkers for predicting the response of immune checkpoint inhibitors in cancer [43,44]. Our results found that CYP1B1 expression was associated with TMB, MSI, and neoantigen in pan-cancer, especially in LIHC and STAD. ssGESA analysis CYP1B1 was associated with melatonergic and immunerelated pathways and therapeutic response. Recent studies showed that CYP1B1 was involved in the drug resistance of tumor cells, such as paclitaxel and docetaxel [45][46][47]. The clinical trial of CYP1B1-directed vaccination identified that patients with solid and hematologic tumors can benefit from anti-CYP1B1 immunity [48]. Thus, targeting CYP1B1 may be a useful way for the development of anticancer treatment.
However, several issues need to be further explored. First, the clinical and immunological characteristics of CYP1B1 in cancer were based on the public database; the roles of CYP1B1 still needs to be further verified by multi-center data. Second, the function of CYP1B1 has been clarified in vitro and vivo, but the function of CYP1B1 knock-down or its inhibition need to be further explored and the potential mechanism of CYP1B1 in TME remains unclear. Third, although the CYP1B1 has promising predictive power for ICI in 11 clinical cohorts, its predictive power still needs to be validated on a larger number of immune cohorts. Finally, the immunological CYP1B1 in pan-cancer needs additional investigation to clarify its function and processes.

Conclusions
In conclusion, we comprehensively analyzed the clinical and immunological characteristics of CYP1B1 across 33 solid tumors. Our results identified that the dysregulated expression of CYP1B1 was associated with the clinical stage, tumor grade, immune cell infiltration, TMB, MSI, neoantigen, activation of multiple melatonergic and immune-related pathways, and therapeutic resistance. Targeting CYP1B1 might be a promising predictor and molecular target for clinical treatment.  Figure S2. The negative correlation between CYP1B1 expression and tumor grade in LIHC. Figure S3. The correlations between CYP1B1 expression and immunomodulator. (A) The spearman correlations between CYP1B1 expression and immunostimulator in cancer. (B) The spearman correlations between CYP1B1 expression and immunoinhibitor in cancer. The heatmap represents rho value. Red color means positive correlation, blue color means negative correlation. Four associations were showed by dot plots. Figure S4. The correlations between CYP1B1 expression and chemokine and receptor. (A) The spearman correlations between CYP1B1 expression and chemokine in cancer. (B) The spearman correlations between CYP1B1 expression and receptor in cancer. The heatmap represents rho value. Red color means positive correlation, blue color means negative correlation. Four associations were showed by dot plots. Figure S5. The IC50 between high and low CYP1B1 expression in normal tissue.