Identification of Prognostic Biomarkers and Molecular Targets Among JAK Family in Breast Cancer.

BACKGROUND: Janus kinases (JAKs) are a family of non-receptor tyrosine kinases involved in multiple malignancies. However, clinical values of JAKs as prognostic markers and potential mechanism as molecular targets in breast invasive carcinoma (BC) are not completely clarified.
METHODOLOGY: TIMER, UALCAN and GEPIA were used to assess the expression and methylation levels of JAKs in BC. Kaplan-Meier Plotter, bc-GenExMiner, SurvExpress, TRGAted, MethSurv, and SurvivalMeth were used to assess the multilevel prognostic significance of JAKs in breast cancer patients. And cBioPortal, TIMER, STRING, GeneMANIA, NetworkAnalysis, LinkedOmics, DAVID 6.8, and Metascape were applied for multilayer networks and functional enrichment analyses. Correlations between immune cell infiltrates/their gene markers and JAKs were evaluated by TIMER.
RESULTS: We first explored the expression and methylation level of JAKs in breast cancer and found significantly reduced JAK1 and JAK2 expression at mRNA and protein levels, significantly higher JAK3 protein expression, and significantly increased TYK2 expression at mRNA level but decreased at protein level. In addition, hypermethylation of JAK3 and TYK2 and hypomethylation of JAK1 were found in tumor samples. In terms of prognostic values of JAKs in BC patients, low transcriptional levels of JAK1, JAK2, JAK3, and TYK2 indicated worse OS/DMFS/PPS/RFS/DFS, inferior DFS, worse RFS, and shorter OS/DMFS/RFS, respectively. The mRNA signature analysis showed that high-risk group had unfavorable OS/RFS/MFS. Low JAK2 protein level indicated unfavorable DSS/PFS in BC patients. Five CpGs of JAK1, four CpGs of JAK2, 20 CpGs of JAK3, and 13 CpGs of TYK2 were significantly associated with prognosis in BC patients. The DNA methylation signature analysis also suggested worse prognosis in the high-risk group. For potential biological roles of JAKs, interaction analyses, functional enrichment analyses for biological process, cellular component, molecular function, and KEGG pathway analyses of JAKs and their neighbor genes in BC were conducted. Kinase targets, gene-miRNA interactions, and transcription factor-gene interactions of JAKs were also identified. Furthermore, JAKs were found to be significantly related to immune infiltrates as well as the expression levels of multiple immune markers in BC.
CONCLUSION: JAKs showed multilevel prognostic value and important biological roles in BC. They might serve as promising prognostic markers and possible targets in breast cancer.

Introduction

Breast cancer is the second most common cancer worldwide. It accounts for more than one-fifth of all female invasive tumors, and 16% of all cancers in women.1,2 Every year an estimated 1,700,000 new cases are diagnosed, this aggressive disease is the leading cause of cancer-related deaths under the age of 45.3 And according to the expression of estrogen receptor (ER), progesterone receptor (PR) and human epidermal growth factor receptor 2 (HER2), breast cancer is broadly divided into luminal ER positive (luminal A and luminal B), HER2 enriched, and basal-like molecular subtypes.4,5 Despite recent great progress in diagnostic methods and treatment technologies,6–8 long-term survival rate of these patients is still not satisfactory, especially in the metastatic setting (median survival less than 24 months). It would be beneficial to conduct more accurate risk assessment and optimize treatment for patients with breast cancer. Therefore, it is essential to develop and identify valuable prognostic markers and therapeutic targets for breast cancer.

JAK (Janus kinase) is a family of intracellular, non-receptor tyrosine kinases with the molecular weight between 120 kD and 140 kD. Structurally, they can be divided into two parts. At the C-terminal, there are two closely connected JAK active sample domains, and the closest to the C-terminal is the catalytic active reaction region of JAK kinase. At the N-terminal, there are five subdomains (A, B, C, D, and E), known as JH domain, which may be involved in molecular coupling of JAK kinase and cytokine receptor. JAKs can mediate the activation reaction of signal protein after binding of cytokines and receptors.9,10 The members of JAK family are JAK1, JAK2, JAK3, and TYK2. In the amino acid sequences of human JAK family proteins, JAK1, JAK2, and JAK3 have high homology with TYK2. JAKs are involved in autoimmune diseases and malignancies, so the targeting of specific JAK members can be used to treat related diseases.11,12 The activated JAK protein phosphorylates the STAT protein to facilitate the formation of STAT dimer, eventually activating the JAK/STAT pathway.13 Inhibition of the JAK1/STAT3 pathway can effectively inhibit ovarian tumor progression.14 JAK2 also promotes invasiveness and drug resistance in colorectal cancer cells through the JAK2/STAT3 signaling pathway.15 The JAK3 mutation is of great significance as a response marker for targeted therapy of JAK kinase or anti-PD1 immunotherapy.16 TYK2 is involved in tumor immune-surveillance in human cancers,17 and it is necessary for the survival of anaplastic large cell lymphoma cells by activating MCL1 expression.18

Although the JAK family members have been reported to be involved in the progression of various types of cancers, the clinical value of the entire JAK family remains poorly investigated in breast cancer. This study investigated the clinical value of JAK family members as prognostic markers and therapeutic targets in breast cancer, based on multiple public resources and reliable integrative bioinformatics analysis. The flow diagram of this study is shown in .

Materials and Methods

Expression and Methylation Levels of JAKs

The JAK expression and methylation levels were compared between tumor tissues of breast invasive carcinoma (BC) and normal tissues using the DiffExp module of TIMER ()19 and UALCAN () web resource.20 GEPIA ()21 was used for the multiple gene comparison and pathological stage plot of JAK family in BC.

Survival Outcomes and mRNA Levels of JAKs or Family Signature

Kaplan–Meier Plotter () 22 was used to obtain the associations between JAK mRNA levels and overall survival (OS)/distant metastasis-free survival (DMFS)/post-progression survival (PPS)/relapse-free survival (RFS) in breast cancer. The associations between JAK mRNA expression and prognosis in breast cancer cases experiencing different clinical conditions were also assessed using gene chips data of Kaplan–Meier Plotter. The prognostic value of JAK mRNA for OS and disease-free survival (DFS) was also analyzed using the RNA-seq data in bc-GenExMiner tool ().23 SurvExpress ()24 was applied for the prognostic analysis of JAK family signature in BC cases from Van De Vijver Nature 200225 cohort for the OS/recurrence-free survival (RFS)/metastasis-free survival (MFS) and TCGA cohort for OS.

Survival Outcomes and JAK Protein Levels

The survival analysis based on JAK protein levels was performed in TRGAted ().26 Only the JAK2 protein was available on this platform. We explored the relationships between JAK2 protein level and OS/disease-specific survival (DSS)/DFS/progression-free survival (PFS) in BC patients.

Survival Outcomes and DNA Methylation of JAK or Family Signature

A web tool named MethSurv ()27 was used to assess the prognostic value of different CpG methylation patterns of JAK1/JAK2/JAK3/TYK2 in BC patients. The effect of DNA methylation of JAK family signature on BC prognosis was analyzed using SurvivalMeth ().28

Genetic Alteration, Co-Expression, PPI Analyses, and Enrichment Network of JAKs

Genetic alteration of JAKs was available from cBioPortal for Cancer Genomics ().29,30 In addition, co-expression of JAKs was studied using the correlation module of TIMER ().19 Protein–protein interaction (PPI) network was obtained from STRING ()31 and GeneMANIA servers ().32 Visual analytics of enrichment network of the four members of JAK family was based on the NetworkAnalysis ().33

Functional Enrichment Analysis of JAKs and Their Most Frequently Altered Neighboring Genes

The top 50 genes with highest frequency associated with JAKs were available from cBioPortal for Cancer Genomics ().29,30 The integrated network data, including biological process (BP), cellular component (CC), and molecular function (MF) and KEGG pathway, were taken from DAVID 6.8 ().34,35 The functional enrichment analysis of JAKs and their neighboring genes was also conducted on Metascape website ().36

Kinase, miRNA, and Transcription Factor (TF) Targets of JAKs

The top 10 kinase targets of JAKs were also identified from LinkedOmics (),37 and the predicted miRNA or TFs and their connections with JAKs were available from NetworkAnalysis ().33

JAKs and Infiltrating Immune Cells

The “Gene” module of TIMER ()19 was used to visualize the correlation of JAK1, JAK2, JAK3, and TYK2 gene expression with immune infiltration levels in BC. The “SCNA” module provided the comparison of tumor infiltration levels among BC samples with different somatic copy number alterations (CNA) for JAK genes. The “Survival” module was used to explore the clinical relevance of tumor immune subsets in a multivariate Cox proportional hazard model.

Correlation Analysis Between JAK Expression and Immune Cell Markers

The correlation module of TIMER ()19 was used for the correlation analysis between JAK expression and immune cell markers together with the Spearman’s rho value and estimated statistical significance. The related gene markers of multiple immune cells were obtained from the related references38–40 and the CellMarker database ().41

Statistical Analysis

The survival curves related to OS, DMFS, PPS, RFS, DFS, DSS, PFS, and MFS of JAKs were drawn from Kaplan Meier plotter, bc-GenExMiner, SurvExpress, TRGAted, MethSurv, or SurvivalMeth, and the P-values with hazard ratios (HRs) were provided. In NetworkAnalysis, miRNA–gene interaction data were collected from TarBase, while transcription factor and gene ChIP-seq data were collected from ENCODE. Only peak intensity signal <500 and the predicted regulatory potential score <1 were used (using BETA Minus algorithm). In this study, TCGA_BC dataset was used in the “LinkInterpreter” module. Gene Set Enrichment Analysis (GSEA) tool was applied to explore the kinase target of JAKs in the BC, with a minimum number of genes (size) of 3 and a simulation of 500. Gene correlations were evaluated with Spearman correlation coefficients and P-values. A P-value lower than 0.05 was considered statistically significant.

Results

Expression of JAK mRNA and Protein, and Methylation Level in Breast Cancer

We first explored the expression of JAKs and methylation level in tumor and normal samples (Figure 1A–C). As shown in Figure 1A, JAK1 and JAK2 mRNA level was significantly decreased, while TYK2 mRNA level was significantly increased in tumor tissues of BC compared with the corresponding normal tissues. Higher JAK3 mRNA level was found in tumor tissues, but statistical significance was not reached. Protein levels of JAK1, JAK2, and TYK2 were significantly lower in tumor samples, while JAK3 protein level was significantly higher in tumor samples (Figure 1B). Significant hypermethylation of JAK3 and TYK2 as well as hypomethylation of JAK1 were found in BC tumor samples (Figure 1C).

Figure 1

Expression of JAK mRNA (A) and protein (B), and methylation level (C) in the BC cases.

We also compared the relative expression levels of JAK mRNA in BC and found that among the JAK family members the relative expression of JAK1 was the highest, while JAK3 had the lowest expression (). There was a significant correlation between JAK2 level and pathologic stage of BC ().

Prognostic Value of Individual JAK mRNA Expression Level in Breast Cancer Patients

The associations between JAK mRNA expression and prognosis were analyzed using the data from gene chip in Kaplan–Meier Plotter (Figure 2A–P). Low expression level of JAK1 (P =1.4e-06) and TYK2 (P=0.0018) indicated significantly shorter OS in patients with breast cancer (Figure 2A and M), while low expression of JAK2 (P=0.096) and JAK3 (P=0.086) showed a modest association with a worse OS (Figure 2E and I).

Figure 2

Survival curves of OS (A, E, I, M), DMFS (B, F, J, N), PPS (C, G, K, O), and RFS (D, H, L, P) for the expression of JAKs in patients with breast cancer.

There were significant associations between DMFS and JAK1 (P=0.00014) and TYK2 (P=0.029) expression in breast cancer patients (Figure 2B and N). Only JAK1 (P=0.043) was significantly related to PPS in breast cancer patients (Figure 2C).

Low transcriptional levels of JAK1 (P=5e-05), JAK3 (P=0.0022), and TYK2 (P=1.9e-08) in breast cancer patients were significantly associated with worse RFS (Figure 2A, L and P), while a modest association was found between JAK2 and RFS (P=0.085) (Figure 2H).

In addition, we also assessed the prognostic value of JAKs using the data from RNA-seq in bc-GenExMiner. We found low mRNA level of JAK1 and JAK2, which indicated significantly shorter OS and worse DFS in patients with breast cancer (all P<0.001) ().

The Relationship Between JAK mRNA Expression and Clinical Variables in Breast Cancer

We further studied the association between JAK mRNA expression levels and clinical variables in breast cancer. We found that low expression of JAK1 was significantly related to reduced OS in patients with the following subtypes: ER-positive or -negative, PR-positive, TP53 wild-type, basal, luminal B, lymph node positive or negative, grade 2 or 3, and luminal androgen receptor subtype (). Low expression of JAK2 was also significantly related to shorter OS in patients with basal, luminal B, lymph node positive or grade 3 subtype, and inferior DMFS in patients with ER-negative, HER2-negative, basal, luminal B, or luminal androgen receptor subtype (). Low expression of JAK3 was also obviously related to worse OS in patients with basal, grade 3, and basal-like 1 subtype, and inferior DMFS in patients with basal subtype (). Low expression of JAK3 or TYK2 was related to worse OS in patients with basal and luminal B subtype (). The results of prognostic values of JAKs in clinicopathological subtypes regarding OS, DMFS, PPS, and RFS in breast cancer are detailed displayed in –.

Prognostic Value of JAK mRNA Signature in Breast Cancer

The JAK family signature was input for prognostic analysis in SurvExpress. All patients were assigned to the high/low-risk groups based on the score value with an optimal cutoff. The dataset from the cohort by Van De Vijver was applied. For OS (Figure 3A–C), distinct expression patterns of JAK1 and JAK2 were noticed between low- and high-risk groups. The high-risk group displayed a significantly shorter OS compared with the low-risk group (HR=1.91, 95% CI: 1.21–3.03, P=0.005469), which was also confirmed in the TCGA cohort (, HR=1.57, 95% CI: 1.11–2.21, P=0.01007). In addition, for RFS (Figure 3D–F), the expression levels of JAK1, JAK2, and TYK2 were higher in the low-risk group than in the high-risk group, and the high-risk group showed a significantly worse RFS than the low-risk group (HR=2.01, 95% CI: 1.36–2.99, P=0.0005106). For MFS (Figure 3G–I), significantly lower expression patterns of JAK1, JAK2, and TYK2 were noticed in the high-risk group, and the cases in the high-risk group had an inferior MFS compared with the low-risk group (HR=2.16, 95% CI: 1.45–3.21, P=0.0001424).

Figure 3

The prognostic values of JAK member signature in BC via SurvExpress platform. (A, D, G) The mRNA expressions of JAK1/JAK2/JAK3/TYK2 between the high- and low-risk groups; (B, E, H) The heat maps of mRNA expression of JAK1/JAK2/JAK3/TYK2; (C, F, I) Survival curves of JAK signature between the high- and low-risk groups.

Prognostic Value of JAK2 Protein Expression in Breast Cancer

The prognostic values of JAK protein were investigated in the reverse-phase protein arrays data of The Cancer Proteome Atlas (TCPA) Portal. Of note, only JAK2 was available (Figure 4A–D). In fact, the BC patients with low protein expression of JAK2 had unfavorable DSS (HR=0.305, P=0.0068) and PFS (HR=0.48, P=0.014) compared with those with high protein expression.

Figure 4

The prognostic value of JAK2 protein in BC. (A) The survival curve for OS. (B) The survival curve for DSS. (C) The survival curve for DFS. (D) The survival curve for PFS.

Prognostic Value of Single CpG of JAKs in Breast Cancer

The prognostic value of DNA methylation of JAK family in breast cancer was analyzed by MethSurv. The heat map of DNA methylation results is displayed in . Specifically, cg18227442 of JAK1, cg13629201 of JAK2, cg11085762 of JAK3, and cg11369662 of TYK2 showed the highest DNA methylation level. Overall, we found that five CpGs of JAK1, four CpGs of JAK2, 20 CpGs of JAK3, and 13 CpGs of TYK2 were significantly associated with prognosis in breast cancer patients (Table 1).

Table 1

The Prognostic Value of CpG in the JAK Family by MethSurv (P<0.05)

Gene-CpG	HR	LR Test P value
JAK1-Body-Open_Sea-cg25192081	0.515	0.0025
JAK1-5ʹUTR-Open_Sea-cg00420347	1.688	0.012
JAK1-5ʹUTR-OpenSea-cg15997411	0.62	0.018
JAK1-5ʹUTR-Island-cg09061759	1.593	0.038
JAK1-TSS1500-S_shore-cg07453451	1.592	0.029
JAK2-Body-Island-cg03390472	0.535	0.003
JAK2-Body-Island-cg13430118	0.587	0.0096
JAK2-Body-Island-cg13629201	1.831	0.0067
JAK2-1stExon-5ʹUTR-Island-cg15814898	0.59	0.024
JAK3-TSS1500-S_Shore-cg11249283	0.533	0.0035
JAK3-TSS200-S_Shore-cg02285920	0.605	0.017
JAK3-TSS200-S_Shore-cg03491584	0.618	0.015
JAK3-TSS200-S_Shore-cg17423207	0.556	0.0043
JAK3-TSS200-S_Shore-cg21988119	0.624	0.02
JAK3-Body-S_Shore-cg14934242	0.571	0.029
JAK3-Body-S_Shore-cg25623545	1.579	0.032
JAK3-5ʹUTR-1stExon-SShore-cg27035251	0.455	7.50E-05
JAK3-Body-N_Shore-cg19622422	0.542	0.0022
JAK3-5ʹUTR-Island-cg04916091	0.661	0.046
JAK3-5ʹUTR-Island-cg05315206	1.75	0.0051
JAK3-5ʹUTR-Island-cg05635754	0.499	0.00085
JAK3-5ʹUTR-Island-cg10357657	0.677	0.048
JAK3-Body-Island-cg05796838	0.572	0.01
JAK3-Body-Island-cg05950755	0.564	0.0084
JAK3-Body-Island-cg06655414	0.592	0.0084
JAK3-Body-Island-cg11085762	0.603	0.013
JAK3-Body-Island-cg27074995	0.487	0.0025
JAK3-Body-Island-cg27309110	0.598	0.0096
JAK3-Body-S_Shelf-cg12957407	0.581	0.027
TYK2-Body-Island-cg00590639	0.587	0.012
TYK2-Body-Island-cg11575644	0.527	0.005
TYK2-Body-Island-cg15662184	0.537	0.011
TYK2-Body-Island-cg18038269	0.553	0.004
TYK2-Body-Island-cg21480173	0.534	0.0031
TYK2-Body-Island-cg23306063	0.467	0.00014
TYK2-TSS1500-Island-cg06622468	0.55	0.0091
TYK2-TSS200-Island-cg16068727	1.765	0.0046
TYK2-5ʹUTR-1stExon-Island-cg18856252	0.617	0.022
TYK2-3ʹUTR-N_shelf-cg07515324	0.573	0.0058
TYK2-3ʹUTR-N_shelf-cg07515325	0.586	0.014
TYK2-3ʹUTR-N_shelf-cg15197202	0.586	0.014
TYK2-TSS1500-S_Shore-cg26265787	0.334	7.90E-05

Prognostic Value of the DNA Methylation of JAK Family Signature in Breast Cancer

The gene symbols of the JAK family members (JAK1, JAK2, JAK3, TYK2) were input for prognostic analysis in SurvivalMeth. Significant expression patterns of all JAKs were found between the low- and high-risk groups (Figure 5A). The heat map showed that DNA methylation level of JAK2 was the highest (Figure 5B), and the high-risk group displayed an unfavorable prognosis compared with the low-risk group (P=0.0060540538) (Figure 5C).

Figure 5

The prognostic value of the DNA methylation of JAK member signature in BC via SurvivalMeth server. (A) The methylation level of CpGs in the high- and low-risk group. (B) The heatmap of CpG methylation level. (C) The Kaplan–Meier curve of OS.

Genetic Alteration, Co-Expression, Physical Interaction Analyses, and Enrichment Network of JAKs

Comprehensive analyses of molecular characteristics of JAKs were further conducted. First, the genetic alterations of these members were analyzed. As displayed in Figure 6A, JAK1, JAK2, JAK3, and TYK2 were altered in 6%, 12%, 5%, and 9% of the 960 samples, respectively. Altered mRNA expression, amplification, and deep deletion of JAKs were also commonly found in these samples.

Figure 6

The genetic alteration, co-expression, and PPI network analyses of JAKs. (A) Genetic alterations in JAKs in BC using cBioPortal. (B) Correlation heat map of JAKs in BC. (C) PPI network of JAKs using STRING. (D) Physical interaction network of JAKs using GeneMANIA.

Then, the potential co-expressions of JAKs were explored. There was a strong correlation between JAK1 and JAK2, and low correlations were found among other JAK members (Figure 6B). Moreover, protein/protein interactions (PPI) of JAKs were explored using STRING. As expected, four nodes and six edges were obtained in the PPI network (Figure 6C). The physical interaction analyses of JAKs were also performed with GeneMANIA (Figure 6D). Moreover, the enrichment network of JAKs was done with NetworkAnalysis, and the integrated interactive visualization with the top 10 highly enriched items of JAKs (P<0.05) is displayed in . The results revealed that the functions of these four members were mainly related to immune response, immune system, binding and kinase activity in GO items (–) and Th1 and Th2 cell differentiation, Th17 cell differentiation, and JAK-STAT signaling path in KEGG analysis ().

Moreover, the top 50 genes with the highest frequency associated with JAKs were available from cBioPortal. The data suggested that TP53, TTN, PIK3CA, SYNE2, CDH1, COL4A2, HMGB3, ACACA, HMGXB4, LY75, PRKDC, CUBN, MYO9A, ATM, ADGRG4, SDK1, FRMPD4, MUC12, GATA3, NEB, MUC16, RNF213, NHS, CTCF, DNAH11, OBSCN, MDN1, MAP3K1, ABCC1, ANKRD36, APLP1, BCAS2, BEND5, CACNA1S, CATSPERB, CCDC138, CCDC34, CCT4, CRX, CWC27, CYP4V2, DENND2D, DGCR8, DNAJB7, FADS1, FAM49A, FGFR2, GAS2, GOSR1, and GRIN2B were primarily connected with the modulation and function of JAKs in breast cancer. The integrated network using DAVID 6.8 with these 54 genes is shown in .

In the BP category, the most highly enriched GO items, including MAPK cascade, innate immune response, phosphatidylinositol-3-phosphate biosynthetic process, and response to antibiotic, might be associated with the tumorigenesis and progression of BC (). In the CC category, the most significantly highly enriched items included extrinsic component of cytoplasmic side of plasma membrane, nucleolus, and cytoskeleton (). In the MF category, the JAKs and their neighboring genes were mainly enriched in binding and kinase activities (). As expected in KEGG pathway analysis (), hsa04151 (PI3K-Akt signaling pathway), hsa04550 (signaling pathways regulating pluripotency of stem cells), hsa04630 (JAK/STAT signaling pathway), and hsa05200 (pathways in cancer) were significantly related to the tumorigenesis and development of breast cancer.

In addition, we also performed the functional enrichment analysis in Metascape website. and show that the functions of JAKs and their neighboring genes were mainly enriched in ST interleukin 4 pathway, PID IFNG pathway (IFN-gamma pathway), and regulation of innate immune response. To better understand the correlation between JAKs and BC, the PPI network and mCODE components were also analyzed. We found that biological function was mainly associated with ST interleukin 4 pathway, interleukin receptor SHC signaling, interleukin-2 family signaling, and interleukin-20 family signaling ().

Kinase Targets, Gene–miRNA Interaction, and Transcription Factor (TF)-Gene Interaction of JAKs

In addition, the top 10 kinase targets of JAKs were also identified from LinkedOmics (Table 2). We also further explored potential miRNAs and TFs of JAKs as well as their interactions using NetworkAnalysis (Figure 7). The predicted miRNAs or TFs are shown in Figure 7A and B. Namely, there were 75 miRNA targets predicted in JAK1, 25 miRNAs in JAK2, 10 miRNAs in JAK3, and 23 miRNAs in TYK2. There were 15 TF targets predicted in JAK1, 22 TFs in JAK2, 6 TFs in JAK3, and 27 TFs in TYK2. miRNAs or TFs were also found to synchronously correlate with all JAK members or at least with three of them.

Table 2

The Top 10 Kinase Targets of JAKs in BC (P<0.05)

JAK1	Gene Set	Description	Size	Leading Edge Number	P value
	Kinase_MAPK7	mitogen-activated protein kinase 7	30	11	0
	Kinase_MAPK3	mitogen-activated protein kinase 3	171	62	0
	Kinase_MAPK1	mitogen-activated protein kinase 1	197	67	0
	Kinase_PRKACG	protein kinase cAMP-activated catalytic subunit gamma	84	24	0
	Kinase_PRKACB	protein kinase cAMP-activated catalytic subunit beta	85	27	0.002519
	Kinase_PRKX	protein kinase, X-linked	86	24	0.002538
	Kinase_ABL1	ABL proto-oncogene 1, non-receptor tyrosine kinase	83	34	0.002577
	Kinase_CDK5	cyclin dependent kinase 5	66	21	0.002639
	Kinase_SYK	spleen associated tyrosine kinase	35	23	0.002793
	Kinase_NUAK1	NUAK family kinase 1	4	2	0.003165
JAK2	Kinase_LYN	LYN proto-oncogene, Src family tyrosine kinase	50	29	0
	Kinase_LCK	LCK proto-oncogene, Src family tyrosine kinase	43	26	0
	Kinase_SYK	spleen associated tyrosine kinase	35	21	0
	Kinase_FGR	FGR proto-oncogene, Src family tyrosine kinase	12	8	0
	Kinase_FYN	FYN proto-oncogene, Src family tyrosine kinase	66	28	0
	Kinase_EGFR	epidermal growth factor receptor	46	20	0
	Kinase_HCK	HCK proto-oncogene, Src family tyrosine kinase	23	13	0
	Kinase_MYLK	myosin light chain kinase	10	5	0
	Kinase_MYLK3	myosin light chain kinase 3	10	5	0
	Kinase_MYLK4	myosin light chain kinase family member 4	10	5	0
JAK3	Kinase_LCK	LCK proto-oncogene, Src family tyrosine kinase	43	27	0
	Kinase_LYN	LYN proto-oncogene, Src family tyrosine kinase	50	25	0
	Kinase_HCK	HCK proto-oncogene, Src family tyrosine kinase	23	12	0
	Kinase_ZAP70	zeta chain of T-cell receptor associated protein kinase 70	12	4	0
	Kinase_FGR	FGR proto-oncogene, Src family tyrosine kinase	12	6	0
	Kinase_FYN	FYN proto-oncogene, Src family tyrosine kinase	66	21	0
	Kinase_SYK	spleen associated tyrosine kinase	35	16	0
	Kinase_FER	FER tyrosine kinase	8	4	0
	Kinase_PKM	pyruvate kinase, muscle	4	2	0
	Kinase_CSK	c-src tyrosine kinase	4	2	0
TYK2	Kinase_IKBKB	inhibitor of nuclear factor kappa B kinase subunit beta	28	13	0.005263
	Kinase_PRKD3	protein kinase D3	8	4	0.009259
	Kinase_CHUK	conserved helix-loop-helix ubiquitous kinase	27	9	0.018717
	Kinase_NEK9	NIMA related kinase 9	4	2	0.019512
	Kinase_NUAK1	NUAK family kinase 1	4	2	0.0199
	Kinase_MKNK1	MAP kinase interacting serine/threonine kinase 1	7	2	0.020725
	Kinase_BCR	BCR, RhoGEF and GTPase activating protein	12	5	0.02649
	Kinase_FGFR1	fibroblast growth factor receptor 1	7	4	0.03681
	Kinase_MET	MET proto-oncogene, receptor tyrosine kinase	8	4	0.038674
	Kinase_RPS6KB1	ribosomal protein S6 kinase B1	23	5	0.042553

Figure 7

The predicted miRNAs and TFs of JAKs and their interaction networks. (A) The predicted networks of miRNAs. (B) JAKs and the predicted networks of TFs and JAKs. Red circle, JAK family; square, predicted miRNA or TFs; line, predicted interactions.

Infiltrating Immune Cell and JAK Expression in Breast Cancer

Given the close relationship between immunological features and prognosis in cancer, we further investigated the correlation between immune infiltrates and JAK expression in BC (Figure 8). Overall, JAK1, JAK2, and JAK3 expression was highly positively associated with all six immune infiltrates (B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells) (all P<0.05). A significant positive correlation was found between TYK2 expression and four kinds of immune cell infiltration, including B cells, CD4+ T cells, neutrophils, and dendritic cells (P< 0.05). In different breast cancer subtypes, including basal-like, HER2+, and luminal, the correlations between JAK expression and the six types of infiltrating immune cells were not all the same (Figure 8).

Figure 8

The correlations between JAKs and immune cell infiltration in BC. (A) JAK1 and immune cell infiltration. (B) JAK2 and immune cell infiltration. (C) JAK3 and immune cell infiltration. (D) TYK2 and immune cell infiltration.

We also compared the tumor infiltration levels among BC with different somatic copy number alterations for JAKs (). Moreover, the multivariate Cox proportional hazard model was applied for JAKs and six tumor-infiltrating immune cells in BC and the subtypes. Overall, macrophages (P=0.059) and JAK2 (P=0.075) were marginally associated with the clinical outcome of all BC patients (). In addition, macrophages (P=0.010), B cells (P=0.028), and dendritic cells (P=0.042) were significantly related to the clinical outcome in basal, HER2, and luminal subtype, respectively (–).

Correlation Analysis Between JAK Expression and Various Markers of Immune Cells

To investigate the effects of JAK expression on tumor-infiltrating immune cells, we analyzed the correlations between JAK1/JAK2/JAK3/TYK2 expression and various markers of immune cells via public databases. The innate and adaptive immune cells investigated included monocytes, tumor-associated macrophages (TAMs), M1 and M2 macrophages, neutrophils, natural killer cells, dendritic cells, CD8+ T cells, B cells, and T cells (general), as well as the subsets of T cells (Th1, Th2, Tfh, Th17, Tregs, and exhausted T cells). The correlation was also adjusted for tumor purity due to its influences on the immune infiltration analysis.

We found that the JAK1 expression significantly correlated with the expression of gene markers in TAMs and Treg in BC patients after adjusting for tumor purity (Figure 9). The correlations between the expression level of JAK2 and the expression of gene marker sets of TAMs, natural killer cells, dendritic cells, CD8+ T cells, T cells (general), Th1, Tfh, Th17, Treg, and T cell exhaustion were still significant after adjustment for purity (). JAK3 expression was significantly associated with TAMs, M2, natural killer cells, dendritic cells, CD8+ T cells, T cells (general), and Th1 and T cell exhaustion in BC after adjustment for purity (). Regarding the correlation between TYK2 and immune cells in BC cases, significant associations were found for dendritic cells, CD8+ T cells, T cells (general), Th2, Treg, and T cell exhaustion (). The results strongly indicated that JAKs might influence prognosis of BC through regulating immune response.

Figure 9

Correlation between the expression of JAK1 and marker genes of infiltrating immune cells in BC using the TIMER database. (A) Correlation between the expression of JAK1 and immune molecular genes. NS >0.05; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. (B) The scatter plots of correlation between JAK1 expression and the gene markers of TAMs and Treg in BC.

Discussion

Janus kinases are a family of intracellular, nonreceptor tyrosine kinases, which comprise four members (JAK1, JAK2, JAK3 and TYK2). JAKs were initially identified using polymerase chain reaction and verified by the human genome project.10,42–45 JAKs possess two near-identical phosphate-transferring domains and mainly exertive functions through transducing cytokine-mediated signals via the JAK-STAT pathway.46 JAKs have been reported to be implicated in the pathogenesis of inflammatory and immune disorders as well as malignant tumors.47–49 And the clinical significance of JAKs in breast cancer was also reported. For example, JAK1 expression inversely correlates with tumor size, lymph node status, and TNM of breast cancer patients.50 And somatic mutations of JAKs (including JAK1, JAK2 and JAK3) have been commonly found in breast cancer and have potential implications for clinical management.51–53 TYK2 plays important roles in growth and metastasis of breast cancer.54,55 JAK inhibitors, such as CP-690,550 (JAK3 target), Lestaurtinib (multikinase target), and AZ-01/AZ-60 (JAK2 target) are currently under development for the treatment of autoimmune diseases and hematological malignancy.46,56,57 JAKs are closely associated with cytokine receptors, such as IFNγ, IL-2, IL-4 and IL-10,58,59 and have important immunologic and inflammation functions.60 In addition, JAK expression levels were also associated with the therapeutic effect in cases treated with anti-PD1 immunotherapy. The anti-PD1 therapy could improve the survival time of mice with NK/T-cell neoplasms and the level of JAK2 was greatly increased in mice treated with anti-PD1.61 Another study found that in patients with non-small cell lung cancer (NSCLC), JAK3 mutations upregulated PD-L1 expression, which explains their benefit from anti-PD1 treatment.16 However, the prognostic significance and biological functions of JAKs in breast cancer have not been comprehensively analyzed.

First, we examined the expression levels of JAKs in the breast cancer. Compared to normal tissues, tumor samples showed significantly decreased JAK1 and JAK2 expression at both mRNA and protein levels, significantly higher JAK3 expression at protein level, and significantly increased TYK2 expression at mRNA level but decreased at protein level. Besides, JAK3 and TYK2 showed hypermethylation, while JAK1 showed hypo-methylation in BC tumor samples. Moreover, in BC tissues, the relative transcriptional level was the highest for JAK1 and the lowest for JAK3, and there was a significant correlation between JAK2 level and BC pathologic stage. These data suggested that JAKs might play multiple important roles in the breast carcinoma tumorigenesis and progression.

Then, we provided systematic prognostic landscape of the expression level of JAK mRNA, protein, and DNA methylation status in breast cancer via multiple databases. First, we found that, through the gene chips data from K-M plotter, low transcriptional level of JAK1 was significantly related to worse OS/DMFS/PPS/RFS, and a modest association was found between worse OS/RFS and decreased JAK2 mRNA level. JAK3 mRNA level displayed a modest association with shorter OS, and a significant relationship with a worse RFS. Low TYK2 level indicated worse OS/DMFS/RFS. Then, we further assessed the association between JAK mRNA expression levels and OS/DMFS/PPS/RFS in breast cancer cases experiencing different clinical conditions. The results using the RNA-seq data of bc-GenExMiner tools showed that breast cancer patients with low mRNA level of JAK1 or JAK2 had significantly shorter OS and worse DFS. Second, we found low JAK2 protein expression group had unfavorable DSS and PFS. Third, we found that five CpGs of JAK1, four CpGs of JAK2, 20 CpGs of JAK3, and 13 CpGs of TYK2 were significantly related to prognosis in breast cancer patients. The mRNA signature analysis showed unfavorable OS/RFS/MFS in the high-risk group of BC patients, and the DNA methylation signature analysis indicated that the high-risk group had worse survival than the low-risk group. These results revealed strong evidence for powerful predictive abilities of JAKs as promising prognostic factors in BC patients.

As breast cancer emerges and progresses by a multistep and multifactorial process involving progressively accumulated genetic, epigenetic, and microenvironmental alterations,62–65 we further performed comprehensive analyses of the molecular characteristics, interaction networks, and the immune microenvironment.

Available from cBioPortal for Cancer Genomics, alterations of JAK1, JAK2, JAK3, and TYK2 were found in 6%, 12%, 5%, and 9% of the BC cases, respectively. Altered mRNA expression, amplification, and deep deletion were also commonly observed. In addition, a strong correlation was found between JAK1 and JAK2, and low correlations were found among other JAK members. These data confirmed that accumulated genetic alterations are involved in the occurrence and progression of breast tumors.

The related molecular biology functions of JAKs were also analyzed. The PPI and the enrichment network of JAKs showed that JAK family was mainly related to immune response, differentiation of immune cells (Th1, Th2, and Th17 cells) and JAK-STAT signaling path. Th1, Th2 and Th17 cells, the subsets of CD4+ T-cell, execute diverse important immune functions.66 It is known that CD4⁺T cells are essential to achieve an effective immune response to pathogens, and the differentiation of naive CD4⁺T cells mainly depends on the cytokine microenvironment, the extracellular matrix and transcription factors.66–71 In breast cancer, Th1 produces IFN-γ and IL-12, exerting anti-tumor activity. Th2 produces IL-4 and IL-10 to promote tumor activity.72 The Th1/Th2 balance is closely related to tumor immunity in breast tumor.72,73 Th17 cells, as the most important factor involved in the inflammation, play significant roles in multiple diseases by production of IL17 cytokine, including breast cancer.74,75 And the JAK-STAT signaling pathway is involved in various fundamental processes such as cell division and apoptosis, stem cell niche maintenance, inflammation and immunity, as well as tumor formation and cancer drug resistance.76–80 The functional enrichment analysis of JAKs and their most frequently altered neighboring genes was also conducted. The mainly highly enriched items included MAPK cascade, innate immune response, binding and kinase activities in BP/CC/MF categories, and hsa04151 (PI3K-Akt signaling pathway), hsa04630 (JAK-STAT signaling pathway), and hsa05200 (pathways in cancer) in KEGG pathway analysis. They play critical roles in diverse cellular functions, including cell growth, proliferation, and survival, and have also been implicated in therapeutic resistance, angiogenesis, tumor migration, invasion, and metastasis.81–85 Moreover, we also sought to characterize the kinase targets, potential miRNAs, and TFs of JAKs as well as their interactions. The top 10 kinase targets of JAKs were identified, and there were 75 miRNA targets predicted in JAK1, 25 miRNAs in JAK2, 10 miRNAs in JAK3, and 23 miRNAs in TYK2. Fifteen TF targets were predicted in JAK1, 22 TFs in JAK2, 6 TFs in JAK3, and 27 TFs in TYK2. The interactions of JAK-miRNA and JAK-TF and the potential kinase targets might be involved in the formation and development of breast tumors.

Moreover, the correlation between JAKs and immune cell infiltration in BC was also assessed. JAK1, JAK2 and JAK3 expression were found to be highly positively associated with six immune infiltrates (B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils and dendritic cells), and TYK2 expression was highly positively related to infiltrations of B cells, CD4+ T cells, Neutrophils, and Dendritic cells. In different breast cancer subtypes, including basal-like, HER2+ and luminal, the correlations between JAKs expression and immune cells infiltrates were also analyzed. We also explored the correlation between somatic CNA of JAKs and abundance of immune infiltrates. In addition, Macrophage, Dendritic cells, and B_cell infiltrates could be an independent factor for clinical survival in basal, luminal, and Her2 cases, respectively. We also performed correlation analysis between JAKs expression and various markers of immune cells. After adjusting for tumor purity, we found that JAK1 expression was significantly associated with the expression of gene markers of TAMs and Treg in BC patients. JAK2 expression correlated with the expression of gene marker sets of TAMs, natural killer cells, dendritic cells, CD8+ T cells, T cells (general), Th1, Tfh, Th17, Treg, and T cell exhaustion. JAK3 expression was significantly related to TAMs, M2, natural killer cells, dendritic cells, CD8+ T cells, T cells (general), Th1, and T cell exhaustion. Significant correlation was found between TYK2 and some of the immune cells (dendritic cells, CD8+ T cells, T cells [general], Th2, Treg, and T cell exhaustion) in BC. These results above indicated JAKs might have an impact on survival of BC patients through regulating immune response and tumor-infiltrating lymphocytes (TILs).

In previous studies, TILs were demonstrated to be closely related to clinical outcomes of BC patients, a systematic evaluation of these cell infiltrations may be able to guide the prognosis and appropriate treatment of breast cancer.86–88 In this work, our findings further showed the expression of JAKs was associated to the infiltration of immune/inflammatory cells. JAKs could be used as prognostic factors and immunotherapy targets in breast cancer.

Our results provided new insights into the identification of JAKs as promising drug targets and powerful prognostic markers in breast cancer. However, some limitations to our study should be acknowledged. First, there might be potential biases caused by the confounders from bioinformatics databases. Second, the findings in our work are based on the data from public databases, and the exact mechanisms need to be validated in future in vitro or in vivo experiments.

In conclusion, the observations of this work suggest a potential role of JAKs as prognostic factors and possible targets in breast cancer.