High expression of JAM2 indicates better prognosis and immunotherapy response in breast cancer

In our study, multiple databases were used to explore the potential role and underlying mechanism of junctional adhesion molecule B (JAM2) in breast cancer (BRCA). The data of JAM2 was downloaded from The Cancer Cell Line Encyclopedia (CCLE), the Genotype-Tissue Expression (GTEx), The Cancer Genome Atlas (TCGA) and Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) databases. Receiver operating characteristic (ROC) curve analysis was performed to analyze the area under the curve (AUC) of JAM2 expression correlated with normal breast tissue and breast cancer tissue. Gene set enrichment analysis (GSEA) was used to identify the potential biological mechanisms of the JAM2. The expression of JAM2 mRNA were downregulated In most tumors, including BRCA, which may be due to the hypermethylated status. The AUCs, which were 0.929 and 0.887 by the logistic regression and random forest algorithms, indicated that JAM2 mRNA expression have good diagnostic value in BRCA. Univariate and multivariate analyses indicated JAM2 as an independent prognostic factor for the overall survival of BRCA patients in both the TCGA cohort (HR = 0.62, P = 0.034) and METABRIC cohort (HR = 0.77, P = 0.001). GSEA showed that multiple tumor pathways were suppressed in the JAM2 high expression group. The expression of JAM2 was most positively related to the epithelial-mesenchymal transition (EMT) score (r = 0.38; P <0.01) by the reverse-phase protein array (RPPA) analysis. Patients with high JAM2 expression may be more sensitive to immunotherapy. 18 chemotherapy drugs that patients in the JAM2 low expression group were more sensitive to were identified. Our results demonstrated the diagnostic and prognostic value of JAM2. Analysis of the molecular mechanisms indicates the potential role of JAM2 as a tumor suppressor, and high JAM2 expression may predict a better immunotherapy response in BRCA.


Introduction
Breast cancer (BRCA), one of the most prevalent malignancies in women, is one of the most common cancers worldwide along with lung and colon carcinoma and primarily contributes to malignancy-related mortality in women [1,2]. In 2020, approximately 276,480 new cases of BRCA will be diagnosed among women in the US, and 42,170 of these women will consequently die from it [3]. Mortality from BRCA has continuously decreased since 1989 due to improvements in treatment and early detection by mammography [4]. However, not all patients have benefitted uniformly from these advances, as mostly indicated by the distinct divergence in the application of adjuvant hormonal therapy due to the versatile classification of BRCA. Currently, BRCA is dominantly divided into 4 distinct subtypes according to the status of estrogen receptor (ER), progesterone receptor (PR), human epidermal growth factor receptor 2 (HER2) and Ki-67, and the treatment regimens for different subtypes vary dramatically [5]. Although advances have been dedicated to therapeutic regimens and the discovery of new drugs, the outcomes of BRCA patients with metastasis remain obstinately poor and have hardly improved. In fact, distant metastasis is responsible for over 90% of BRCA-related mortality [6]; hence, identifying metastatic targets for BRCA is of great interest and is essential.
Junctional adhesion molecule B (JAM2), a member of the immunoglobulin superfamily (IgSF), is located on the membrane of eukaryotic cells and enriched at cell junctions in epithelial and endothelial cells [7]. Based on its biofunction, JAM2 has been reported to participate in multiple biological processes, including tight junction assembly, regulation of paracellular permeability, angiogenesis, tumor metastasis and cell proliferation [8][9][10]. Structurally distinguished from other JAM family members, JAM2 could specifically interact with JAM-C and integrin α4β1 and subsequently form multimer interactions with integrin counter-receptors [11]. Indeed, it was reported that diminishing the interactions of JAM-B and JAM2 significantly reduced monocyte numbers in the abluminal compartment by increasing reverse transmigration [12].
Moreover, Ebnet et al identified the interaction between PAR-3 and the JAM2 and JAM-C complex at its first PDZ domain [13]. Additionally, JAM2 was reported to negatively regulate pro-angiogenic pathways by inhibiting VEGF-induced ERK1/2 activation [14]. Importantly, JAM2 has been demonstrated to play an important role in the metastasis of multiple tumor cells. Hajjari et al showed that the JAM2 and JAM-C complex could activate the expression of the actin filament-associated protein (AFAP) gene and that the dysregulation of JAM2 and JAM-C potentially contributed to the progression of gastric adenocarcinoma tumors [15]. In addition, the interaction of JAM2 and JAM-C could activate c-Src, a proto-oncogene, and therefore regulate cell invasion and migration in glioma [16]. Interestingly, overexpression of JAM2 inhibited VEGFinduced angiogenesis and thus attenuated the progression of lung carcinoma in vivo [17].
However, the expression levels of JAM2 vary dramatically among different tumor types, and the role of JAM2 in breast cancer has rarely been evaluated.

Data collection
We downloaded 21

Gene set enrichment analysis (GSEA)
The BRCA samples were separated into high and low JAM2 expression groups based on the optimal cutoff value of JAM2 by the surv_cutpoint function of the survminer R package (https://CRAN.R-project.org/package=survminer). The clusterProfiler R package was used to perform GSEA [20]. The annotated gene set files "c2.cp.kegg.v7.1.entrez.gmt" and "h.all.v7.1.entrez.gmt" were used as references. The significance threshold was set as p value < 0.05.

Estimation of immunotherapy and chemotherapy response
The estimated scores, immune scores and stromal scores of BRCA were analyzed by the estimate R package [21]. The Immune Cell Abundance Identifier (ImmuCellAI) (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/), which uses a gene set signaturebased method to precisely estimate the infiltration score of 24 immune cell types, was used to predict the immunotherapy response (anti-PD1 or anti-CTLA4 therapy) of BRCA in the TCGA and METABRIC cohorts [22]. The chemotherapy response was estimated by the half maximal inhibitory concentration (IC50) of each BRCA patient from the TCGA cohort through the pRRophetic R package [23].

Statistical analysis
All statistical analyses were performed using R version 4.0.0 (2020-04-24). For continuous variables, the Kruskal-Wallis H test was performed. Univariate and multivariate Cox regression analyses were used to identify the predominant prognostic factors of overall survival (OS) (P < 0.05). Kaplan-Meier survival curves were compared using the log-rank test. The ggplot2 R package (https://CRAN.Rproject.org/package=ggplot2) was used to plot the correlation bubble diagram. P < 0.05 (two-sided) was considered statistically significant.

The mRNA expression of JAM2 from pan-cancer analyses of different databases
Three different databases were used to explore the mRNA expression of JAM2 in different tumors and normal tissues. First, the data of each tumor cell line were downloaded from the CCLE database, and the expression of JAM2 was observed in 21 different tissues ( Fig 1A). In addition, the GTEx database was used to analyze the expression levels of JAM2 in 31 normal tissues, as shown in Fig 1B. Interestingly, compared to that in other tumor cells, the expression of JAM2 in BRCA cell lines was relatively low, but compared to that in other normal tissues, the expression of JAM2 in BRCA cells was relatively high. Considering the small number of normal samples in TCGA, we integrated the normal tissues from the GTEx database into the TCGA database to analyze the differential expression of JAM2 in 27 tumor and normal tissues ( Fig 1C). Except for six tumors (marked in blue), the expression of JAM2 was low in all other tumors (marked in yellow) compared to normal tissues.  Table S1). The results indicated that the expression of JAM2 was downregulated in BRCA tissues compared to normal tissues in the TCGA cohort (Fig 2A). Moreover, the expression of JAM2 was altered in subtypes of breast cancer; among them, it was highly expressed in the claudin-low subtype and expressed at low levels in the basal subtype ( Fig 2B). The same result was confirmed in the METABRIC database ( Fig 2C,   2D). In addition, to identify the diagnostic and prognostic value of JAM2 in BRCA, receiver operating characteristic (ROC) curve analysis was performed to analyze the area under the curve (AUC) of JAM2 expression correlated with normal breast tissue and breast cancer tissue by two algorithms. As shown in Fig 2E, the AUCs were 0.929 and 0.887 by the logistic regression and random forest algorithms, respectively, in the TCGA cohort, indicating that JAM2 has good clinical diagnostic value. Furthermore, the patients were separated into either the JAM2 high expression group or the JAM2 low expression group based on the optimal cutoff value, which was derived from the surv_cutpoint function of the survminer R package. The results indicated that patients with low JAM2 mRNA expression had a lower OS rate (Fig 2F; P = 0.016). Next, the METABRIC dataset was used to validate the results, and the same results were produced (Fig 2G, 2H).
Univariate and multivariate analyses were also performed to identify the prognostic value of JAM2. Multivariate Cox regression analysis was conducted by utilizing the significant prognostic factors identified in the univariate analysis (P < 0.05). The expression of JAM2 was determined to be an independent prognostic factor (hazard ratio (HR) = 0.62, P = 0.034) in the TCGA cohort ( Fig 3A). Furthermore, tumor stage (HR = 3.06, P < 0.01) and age (HR = 2.29, P < 0.01) were also independent prognostic predictors in the TCGA cohort ( Fig 3A). The ER, PR and HER2 statuses were also associated with BRCA prognosis; however, there were too many patients with incomplete in the TCGA cohort. Therefore, the METABRIC dataset, which included BRCA patients and comprehensive clinical information, was used to validate the predictive prognostic value of JAM2. In the METABRIC cohort, after taking into account factors such as ER, PR and HER2 statuses, high JAM2 mRNA expression (HR = 1.88, P = 0.004) was also an independent prognostic factor by multivariate analysis ( Fig 3B).

Hypermethylation of JAM2 in BRCA
To identify the downregulated mechanism of JAM2 in BRCA, multiple methods were used to analyze the methylation status. Three methods, the beta value of the probe with the largest standard deviation (Fig 4A), the mean beta value of all probes within the gene (Fig 4B) and the beta value of the probe with the smallest correlation coefficient with mRNA expression (strongest negative correlation) (Fig 4C), were used to represent the methylation levels of JAM2, and all methods demonstrated that JAM2 was significantly more highly methylated in BRCA tissues than in normal tissues (P < 0.05, Fig 4A, 4B and 4C). Pearson correlation analysis of 4 DNA methyltransferases showed that the expression of JAM2 was negatively related to the expression of DNMT1 (r = -0.12, P <0.01), DNMT3B (r = -0.14, P <0.01) and EHMT2 (r = -0.19, P <0.01) (Fig 4D-G). In addition, the high beta value of the cg01975706 probe, which was the most negatively correlated with the expression of JAM2 (Fig 4H), indicated a worse survival rate (Fig 4I). The above information suggests that the low expression of JAM2 in BRCA may be due to hypermethylation.
GSEA was conducted to identify the potential biological mechanisms of the JAM2 high expression group and the JAM2 low expression group in BRCA. The MSigDB results of c2.cp.kegg.v7.1.entrez.gmt indicated that the gene signatures of "pathways in cancer", "leukocyte transendothelial migration", "TGF beta signaling pathway", "jak/stat signaling pathway", and "basal cell carcinoma" were significantly suppressed in the JAM2 high expression group (Fig 5E) ( Table S2). The MSigDB results of h.all.v7.1.entrez.gmt indicated that the gene signatures of "inflammatory response", "TNFA signaling via NFKB", "epithelial mesenchymal transition", "KRAS signaling down", and "estrogen response early" were significantly suppressed in the JAM2 high expression group (Fig 5F) (Table S3). Overall, these findings strongly indicate the potential role of JAM2 as a tumor suppressor in BRCA.

The use of reverse-phase protein array (RPPA) to explore the potential function of JAM2 in BRCA
The RPPA pathway scores, which correlate with tumor lineage and indicate the protein expression signatures of pathway activity, were obtained by previously published literature in TCGA [24]. The analysis of the correlation between the expression of JAM2 and the pathway scores of RPPA indicated that JAM2 was most positively related to the epithelial-mesenchymal transition (EMT) score (r = 0.38; P <0.01) and most negatively associated with the proliferation score (r = -0.56; P <0.01) ( Fig 6A). In addition, all pathway scores were significantly related to the expression of JAM2 except for the hormone a score, apoptosis score and PI3K/Akt score (Fig 6A) ( Table S4). Our analysis also indicated that the pathway scores for EMT, Ras/MAPK and RTK were significantly higher in the high JAM2 expression group. However, the pathway scores for tumor purity, proliferation, cell cycle, DNA damage response and TSC/mTOR were significantly lower in the high JAM2 expression group (Fig 6B-I).
Furthermore, correlation analysis was also performed between the expression of JAM2 and 226 tagged proteins in RPPA. The results indicated that 59 proteins were positively correlated with JAM2 (P < 0.01) and 72 proteins were negatively correlated with JAM2 (P < 0.01) (Fig 7A) ( Table S5). The three proteins most positively and negatively associated with JAM2 are shown in Fig 7 B-G.

Relationships between JAM2 expression and neoantigens, tumor mutational burden (TMB) and microsatellite instability (MSI)
Most cancers are associated with mutations and genomic instability, such as chromosomal instability (CIN) or MSI [25]. Neoantigens are nascent antigens encoded by mutant genes in tumor cells, mainly produced by point mutations, deletion mutations, gene fusions, etc., that are not the same as proteins expressed by normal cells [26].
Neoantigens can be utilized to identify tumor-specific T cells or to generate vaccines [27]. The TMB was estimated by analyzing somatic mutations, including single-nucleotide variants (SNVs) and insertions/deletions (INDELs) per megabase of the panel sequences examined, and can be a potential predictive biomarker for immunotherapy [28]. Therefore, the MSI scores, number of neoantigens and TMB in each TCGA sample were calculated, and significant negative correlations were found between JAM2 expression and MSI (r = -0.08; P = 0.013), neoantigens (r = -0.18; P < 0.001) and TMB (r = -0.09; P = 0.007) (Fig 8A-C).

Immunotherapy and chemotherapy responses of high-and low-JAM2 patients in BRCA
Increasing reports indicate that the tumor immune microenvironment plays an important role in tumor development. Therefore, the estimate scores, immune scores and stromal scores were analyzed by the estimate R package. The correlation analysis indicated that the expression of JAM2 was positively associated with these genes (all P <0.05) (Fig 8D-F). Tumor-infiltrating lymphocytes are independent predictors of lymph node status and survival in cancer patients. Twenty-eight immune infiltration cell scores and 6 immune infiltration cell scores were downloaded from The Cancer Immunome Atlas (https://tcia.at/) [29] and Tumor IMmune Estimation Resource (TIMER) [30], respectively. The results indicated that except for the CD4 T cell score, all immune infiltration cell scores were positively correlated with JAM2 expression (Fig 8G-H). Common immune checkpoint genes were also collected to analyze the relationship between JAM2 expression and immune checkpoint gene expression. The correlation bubble diagram showed that most of the immune checkpoint genes were significantly positively correlated with JAM2 ( Fig 8I).
Immunotherapy using immune checkpoint blockade targeting CTLA-4 and PD-1 has emerged as a promising strategy for the treatment of various malignancies [31].
Therefore, ImmuCellAI [22], which is a gene set signature-based method, was used to estimate the clinical response to immune checkpoint blockade in the TCGA and METABRIC cohorts. The results indicated that JAM2 expression was higher in those predicted to be sensitive to immunotherapy than in those who were not sensitive ( Fig   9A, 9C), and the immune infiltration scores were also different between the high and low JAM2 expression groups in the TCGA and METABRIC cohorts (Fig 9B, 9D).
The IC50 value of each BRCA patient in the TCGA cohort was predicted by the pRRophetic R package [23]. BRCA patients with a low expression of JAM2 have a worse prognosis. The 18 chemotherapy drugs that patients in the JAM2 low expression group were more sensitive to are shown in Fig 10 (Table S6). These results may provide clinical treatment strategies for future studies.

Discussion
JAM2 belongs to a subfamily of immunoglobulins that regulate tight junction assembly and leukocyte recruitment during inflammation, angiogenesis, cell migration and proliferation [32]. More research has shown that JAM2 is differentially expressed in different types of cancers. Kok-Sin et al found that the expression of JAM2 was downregulated due to its hypermethylated status in colorectal cancer [33]. Similarly, the methylation of JAM2 was increased, and its expression was decreased in esophageal squamous cell carcinoma [34]. However, in some tumors, JAM2 seems to play a Hajjari et al found that JAM2 was significantly upregulated in gastric tumor samples compared with adjacent normal tissues [15].
In our study, multiple databases were used to explore the potential role and underlying mechanism of JAM2 in BRCA. The expression of JAM2 was relatively low in BRCA cell lines compared to other cancer cell lines from the CCLE database (Fig 1A), and the expression of JAM2 was relatively high in normal breast tissue compared to other normal tissues in the GTEx database ( Fig 1B). Interestingly, except for cholangiocarcinoma (CHOL), glioblastoma multiforme (GBM), acute myeloid leukemia (LAML), low-grade glioma (LGG), liver hepatocellular carcinoma (LIHC) and pancreatic adenocarcinoma (PAAD), downregulation was observed in all other tumors compared with that in normal tissues through analyses of the GTEx and TCGA databases ( Fig 1C).
DNA methylation is an imperative epigenetic mechanism that regulates several biological processes. Aberrant DNA methylation in cancer is typically characterized by the promoter-specific hypomethylation of genes, particularly cancer suppressor genes, leading to gene silencing [37,38]. Therefore, the methylation status of JAM2 in BRCA was analyzed by multiple methods to identify the downregulated mechanism of JAM2 in BRCA. Finally, similar to that in esophageal squamous cell carcinoma and colorectal cancer, we found that the expression of JAM2 was downregulated, which may be due to the hypermethylated status of the JAM2 gene in BRCA (Fig 4A-4I).
The role of JAM2 in tumorigenesis and development has also been investigated in previous research. Huang JY et al suggested that JAM2 and JAM-C may potentially be involved in gastric adenocarcinoma tumor progression [39]. Qi LF also observed that downregulated JAM2 could inhibit glioma migration in vitro by blocking the NOTCH pathway [40]. An interesting finding suggested that overexpression of JAM-B is responsible for suppressing VEGF-induced angiogenesis in a Down syndrome trisomy 21 mouse model, thereby restraining the tumor effects in a mouse model of lung cancer [17]. JAM2 is involved in BRCA progression, but the underlying mechanism remains unclear [41]. To identify the potential molecular mechanisms of JAM2 in BRCA, 16 oncogenes and 17 TSGs associated with OS and JAM2 were identified.
GSEA showed that the gene signatures of "epithelial mesenchymal transition", "pathways in cancer", "jak/stat signaling pathway" and "basal cell carcinoma" were significantly suppressed in the JAM2 high expression group. These findings strongly indicate the potential role of JAM2 as a tumor suppressor in BRCA.
The analysis of proteomic data is one of the commonly used tools in the field of tumor pathophysiology, and the RPPA technique is widely used in the TCGA project to study the tumor proteome [24]. The three proteins most positively and negatively associated with JAM2 were caveolin 1, collagen VI, and PEA15 and MSH6, ASNS, and RBM15.
Furthermore, Ludwig et al found that JAM2 plays a role in T lymphocyte rolling and firm adhesion by interacting with α4β1 [42]. Therefore, the estimate scores, immune scores and stromal scores were analyzed by the estimate R package, and the expression of JAM2 was positively associated with the abovementioned scores (all P <0.05). In addition, except for the CD4 T cell score, all immune infiltration cell scores were positively correlated with JAM2 expression in different databases. We estimated the clinical response to immune checkpoint blockade in the TCGA and METABRIC cohorts and found that JAM2 expression was higher in those predicted to be sensitive to immunotherapy than in those who were not sensitive. These results suggest that JAM2 plays a broad role in BRCA inflammatory infiltration and that high expression of JAM2 may predict a better immunotherapy response.
Finally, the pRRophetic R package was used to predict the IC50 of each BRCA patient in the TCGA cohort. Eighteen chemotherapy drugs that patients in the JAM2 low expression group were more sensitive to were found to provide clinical treatment strategies in the following study.

Conclusions
Our results demonstrated the diagnostic and prognostic value of JAM2. Analysis of the molecular mechanisms indicates the potential role of JAM2 as a tumor suppressor, and high JAM2 expression may predict a better immunotherapy response in BRCA.

Data Availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest
The authors declare that they have no competing interests.

Funding Statement
The present study was supported by grants from the National Natural Science Foundation of China (nos. 81772979 and 81472658).