The impact of JAK/STAT inhibitor ruxolitinib on the genesis of lymphoproliferative diseases

Background/aim Ruxolitinib, a JAK/STAT signaling pathway inhibitor targeted drug, has been approved for the controlling of disease symptoms and splenomegaly in patients with myeloproliferative neoplastic diseases. Recently, it has been proposed that ruxolitinib-induced JAK/STAT pathway inhibition in myelofibrosis is associated with an elevated frequency of aggressive B-cell lymphomas. However, the biological basis and significance of this pharmacobiological adverse event is unknown. The aim of this bioinformatics study is to detect any possible confounding effects of ruxolitinib on the genesis of lymphoproliferative disorders. Materials and methods The gene expression data were retrieved from the E-MTAB-783 Cancer Genome Project database. Gene expression data for all available genes in 26 cell lines belonging to various types of lymphomas were chosen for use in this in silico analysis. Results We identified genes that were significant in developing resistance to ruxolitinib in lymphoma cell lines. Conclusion Based on the results of our present study, ruxolitinib may potentially lead to the pathological expression of the transcription factors important in lymphoma genesis, neoplastic commitment on the progenitor lymphoid cells, inhibition of repressor transcriptions protective for lymphoma development, inhibition of apoptosis, promotion of neoplastic proliferation, transcriptional activation, and proliferation of malignant neoplastic B cells.


Introduction
Ruxolitinib, a novel oral JAK 1 and 2 inhibitor, was newly certified as a revolutionary treatment for patients suffering from intermediate/high risk myelofibrosis (1). It was demonstrated that the exposure of cutaneous T-cell lymphoma cell lines triggered a dose-dependent inhibition of cell proliferation by a mechanism that impacted the control of DNA synthesis, with a significant decline in the basal levels of phospho-STAT3 (2).
Furthermore, the ruxolitinib-mediated decrease in circulating proinflammatory cytokine levels is connected to its inhibitory effects on JAK1-mediated signaling; this may make ruxolitinib a tempting therapeutic agent for MYD88 mutation-positive lymphoma patients, as this mutation leads to cell-autonomous activation of JAK 1 and JAK 2 (3). JAK/STAT signaling has been detected to be vital in the mammary gland, lymphocytes, adipocytes, neuronal cells, cardiomyocytes, hepatocytes, eye cells, and stem cells. Deregulation in the JAK/STAT pathway is suggested as a cause of numerous diseases including the emergence and the advancement of cancer cells (4).
The JAK/STAT pathway controls embryonic development and is also a part of the regulation of many processes such as stem cell maintenance, hematopoiesis, and the inflammatory response (5). STATs have critical as well as nonredundant functions in lymphocyte development. They impact cell fate decisions like differentiating naive T cells, regulate the intensity and period of inflammatory responses, and contributes to pathogenic mechanisms in chronic inflammatory diseases. STAT3 governs the differentiation of naive T cells into the regulatory and inflammatory T-cell lineages. It also controls cell growth, apoptosis, and the transcription of inflammatory genes, and it contributes to the development of chronic inflammatory diseases, as well as malignant and neurodegenerative diseases (5). STATs are also known to control counteracting cellular events by contributing to the induction of apoptosis, as well as in differentiation and stem cell maintenance. STAT3 signaling throughout mammary gland involution induces epithelial cell death. In opposition to this, persistent and enhanced STAT3 activation is an initiator of tumorigenesis (5). Moreover, JAK/STAT signaling leads to immunosuppression and controls inflammatory responses, obesity, stem cell maintenance, and the premetastatic niche development. These effectors and the direct and mediated mechanisms of JAK/STAT signaling in and on tumors cells improved our outlooks for JAK/STAT-based cancer therapeutics (5).
The aim of this bioinformatics study is to detect any possible confounding effects of ruxolitinib on the genesis of lymphoproliferative disorders. Recently, it has been proposed that ruxolitinib-induced JAK/STAT pathway inhibition in myelofibrosis is associated with an elevated frequency of aggressive B-cell lymphomas (6). Elucidation of this unique drug acting on multiple biological downstream paths is important since it is already accepted as a standard medication for the real-life management of patients with chronic myeloproliferative disorders, namely myelofibrosis and polycythemia vera (PV) (7,8). The management of such patients is based on randomized clinical trials such as the COMFORT-I, COMFORT-II, and RESPONSE trials and real-life administration data (9)(10)(11). However, chronic myeloproliferative disorders are relatively benign diseases in comparison to the acute myeloproliferation. Therefore, there is a general safety concern for the long-term administration of ruxolitinib, particularly for the emergence of secondary malignancies such as lymphomas (9,10). Hence, the findings of this study can underline the complication of lymphoma development with the probable relative long-term usage of ruxolitinib as a JAK/STAT signal transduction inhibitor (6).

Microarray gene expression and drug cytotoxicity data
The gene expression data used in this study were retrieved from the E-MTAB-783 Cancer Genome Project (CGP) database and RMA-normalized (12). The database contained data for 13,513 genes, amounting to 22,279 probe sets. The expression of these genes in 773 cell lines, belonging to a variety of cancer types, was available. Out of these, gene expression data for all available genes in 26 cell lines belonging to various types of lymphomas (B-cell lymphoma, Burkitt lymphoma, and Hodgkin lymphoma) were chosen for use in in silico analysis. The drug cytotoxicity data for ruxolitinib were obtained from the E-MTAB-783 CGP database for drug sensitivity (12).

Identification of the genes whose expressions correlated with ruxolitinib sensitivity
The previously retrieved and RMA-normalized gene expression data for 13,513 genes in 26 lymphoma cell lines (B-cell lymphoma, Burkitt lymphoma, and Hodgkin lymphoma) were correlated with drug cytotoxicity data for ruxolitinib. Using Pearson correlation analysis, genes that could be associated with ruxolitinib sensitivity were identified.
From this analysis, genes with the highest correlation were then used for further analysis. A list of genes containing all the genes that we were especially interested in investigating was compiled. The selection criteria of genes included clinicopathological correlation of lymphoma genesis and ruxolitinib administration. During the selection of genes, a clinician and molecular biologist worked together in order to determine clinically relevant and biologically important genes to underline the research endpoints. These genes are henceforth referred to as genes of interest in this study and are compiled in Table 1. This list had 200 genes; 155 belonged to the JAK/STAT pathway, 31 were cancer stem cell (CSC) markers, and 14 were immune response markers.
The highly correlated genes from the Pearson correlation analysis were then searched to see if any of them were genes of interest. Following this, the top 10 directly correlated genes and top 10 inversely correlated genes that were highly correlated were selected. Lymphoma cell lines were hierarchically clustered based on these genes, using Gene Cluster 3.0 with Euclidean distance as a similarity metric parameter for genes and arrays, and complete linkage.

Determination of differentially expressed genes within resistant and sensitive groups
We had gene expression data for 26 lymphoma cell lines. We decided to divide these cell lines into two groups based on IC 50 values for ruxolitinib to see if we could distinguish between cell lines that were sensitive to ruxolitinib and cell lines that were resistant to ruxolitinib.
A bar chart and a heat-map were generated in Excel from the IC 50 values to visualize the increase in IC 50 values across the 26 cell lines. A two-tailed, unpaired t-test was performed between cell lines with the lowest IC 50 values, the sensitive group (HDLM2, TUR, BC3), and cell lines with the highest IC 50 values, the resistant group (RAJI, DG75, RPMI6666, DB). P < 0.05 was determined to be an appropriate standard for selection of the two groups.
The 200 genes of interest, containing 155 JAK/STAT pathway genes, 31 CSC markers, and 14 immune response markers, were now analyzed within the sensitive group (HDLM2, TUR, BC3) and resistant group (RAJI, DG75, RPMI6666, DB). A two-tailed, unpaired t-test was used to verify which genes were differentially expressed in the two groups. A cut-off of 0.05 was used for the P-values.
Separate hierarchical clusters were first created for all 26 lymphoma cell lines and then for only the 7 cell lines in the two groups, based on the significant differentially expressed genes. The clustering was done using Gene Cluster 3.0 with Euclidean distance as a similarity metric parameter for genes and arrays, and complete linkage.

Gene-set enrichment analysis
The gene expression data for all 13,513 genes found in the E-MTAB-783 CGP database, and for the 7 cell lines in the sensitive group (HDLM2, TUR, BC3) and resistant group (RAJI, DG75, RPMI6666, DB), were put into GSEA 3.0 to determine which genes were significantly different within the two groups. Furthermore, we wanted to highlight the pathways that these genes were a part of in both these groups (12).

Network analysis
Out of the 200 genes of interest, the 25 genes that were differentially expressed were used as input for the GeneMANIA app (13) in Cytoscape (14) to generate a network based on coexpression and genetic interactions. That was done in order to find if the two sets have functionally similar genes related to each other and to determine the associated functions for different groups of genes in the network. The network served as a starting point for the pathway analysis approach in Reactome (15) to filter out regulatory pathways that are relevant for both these gene sets.

Results
We identified genes that were significant in developing resistance to ruxolitinib in lymphoma cell lines. We used a statistical approach that allowed us to correlate the IC 50 data for ruxolitinib and gene expression for desired genes and create a linear relationship between the two.

Genes correlated with ruxolitinib sensitivity
The gene expression data for all 13,513 genes in the E-MTAB-783 CGP database for lymphoma cell lines were analyzed (12). By doing Pearson correlation analysis with these genes, we were able to find the most highly correlated genes with ruxolitinib sensitivity. A cut-off of 0.4 was used for the R-values so that a large number of genes could be considered. With this cut-off, 800 genes were shown to be highly correlated with a direct relationship while 584 were shown to be highly correlated with an inverse relationship. Further analyses were done using these highly correlated genes. Of these highly correlated genes, we found a number of genes that were of interest to us. We had compiled a list of 200 genes for this study that we wanted to investigate. The list contained 155 JAK/STAT pathway genes, 31 CSC markers, and 14 immune response markers. Our analysis showed that 4 JAK/STAT pathway genes had a direct relationship with ruxolitinib resistance,  IL13RA1, IL13RA2, IL15, IL15RA, IL19, IL2, IL20RA, IL21, IL21R, IL22, IL22RA1,  IL23A, IL24, IL26, IL2RA, IL2RB, IL2RG, IL3, IL3RA, IL4, IL4R, IL5, IL5RA,  IL6, IL6ST, IL7, IL7R, IL9, IL9R, IRF9 *Genes that are a part of the JAK/STAT mediated signaling pathway. **Genes that are known to be cancer stem cell markers within lymphoma. ***Genes belonging to the immune response pathway that are identified as markers for lymphoma. while 8 JAK/STAT pathway genes showed an inverse relationship. Figure 1A shows one of the genes, CSF2RA, which has an inverse relationship with the development of ruxolitinib resistance. On the other hand, we also found a direct linear relationship between the expressions of 5 CSC markers and resistance to the drug. Figure 1B shows the correlation of one of these genes, CD45. Unfortunately, no immune response markers were found among the highly correlated genes. Following this, the ten highest directly correlated genes and the ten highest inversely correlated genes were hierarchically clustered. The results shown in Figure 2 show the presence of two distinct clusters within lymphoma cell lines.

Determination of sensitive and resistant groups within lymphoma cell lines
The 26 available lymphoma cell lines were divided into two groups, where one group was significantly more resistant to ruxolitinib than the other. Significance was determined by performing a two-tailed, unpaired t-test between the groups and accepting P-values of less than 0.05. IC 50 data for ruxolitinib were plotted for each cell line available. As a result, a comparative bar chart was created, as shown in Figure 3. This allowed for the clear determination of the level of sensitivity of each cell line to ruxolitinib. A heat-map was created to further highlight the range of sensitivity to ruxolitinib in lymphoma cell lines, also seen in Figure 3. Based on the bar chart and the heat-map, the cell lines were divided into two groups. The sensitive group (HDLM2, TUR, BC3) included the cell lines that corresponded to the lowest ruxolitinib IC 50 values, hence being the most sensitive to this drug. The resistant group (RAJI, DG75, RPMI6666, DB) included the cell lines that corresponded to the highest ruxolitinib IC 50 values, being the least sensitive and most resistant to action of this drug. The resistant group, compared to the sensitive group, was significantly more resistant to ruxolitinib, with a P-value of 0.001.

Differentially expressed genes within resistant and sensitive groups
The E-MTAB-783 CGP database for gene expression in cancer cell lines was used with a focus on the lymphoma cell lines that were part of the resistant and sensitive groups (12). Analysis of the comparative gene expression between these two groups revealed that 904 genes were differentially expressed, with P-values of less than 0.05. The results further showed that out of these 904 genes, 15 belonged to the list of 200 genes of interest that we had compiled for this study. These 15 genes are displayed in Table 2 along with their classifications; 9 of these genes are involved in the JAK/STAT pathway while 6 are CSC markers. These genes and their clinical significance are shown in Table 3. This allowed us to form a relationship between the expression of genes of interest and ruxolitinib sensitivity. To further demonstrate this relationship, differentially expressed genes for the JAK/STAT pathway and CSC markers were hierarchically clustered, as shown in Figure 4 and Figure 5, respectively. These figures can be seen having two distinct clusters.

Gene-set enrichment analysis
Gene-set enrichment analysis was then performed between the sensitive and resistant groups of the lymphoma cell lines. The analysis was performed with all the genes available in the E-MTAB-783 CGP database (12). The results highlighted the genes that showed significant differences between the two groups and provided us with insight about which pathways these significant genes were involved in.   The cell cytotoxicity data for ruxolitinib available in the E-MTAB-783 CGP database were compared for each lymphoma cell line. The bar chart above shows the vast range of IC 50 for lymphoma cell lines, which allowed us to create two distinct groups. The selected groups were significant with a P-value of 0.001 (20). Table 4 and Table 5 respectively show the pathways that were found to be statistically noteworthy, with P-values of less than 0.05, within the sensitive and the resistant groups.

Network analysis
Our previous analysis had revealed that 15 of 200 genes of interest were differentially expressed in the sensitive and resistant groups. Nine of these belonged to the JAK/ STAT pathway while the other 6 were CSC markers. These 15 genes were input into the GeneMANIA app in Cytoscape to generate a network based on coexpression and genetic interactions. Cytoscape tools were used to differentiate the input genes from both sets from the ones that GeneMANIA predicts as likely to share the same function based on their interactions, as seen in Figure 6. From the network, it is distinctly clear that the two sets of genes are connected with one another. Furthermore, the results of GeneMANIA analysis provided a list of associated functions for different groups of genes in the network, as shown in Table 6. The integration of different    genes together with their function allows for a complete view of potential regulatory mechanisms occurring in a biological process. Reactome analysis was then used to group genes based on their common pathways, demonstrated in Table 7.

Discussion
The JAK-STAT pathway may play an important role in the pathogenesis of lymphoma (2,16). In humans, the activator of the JAK-STAT family comprises four JAKs (JAK 1, JAK 2, JAK 3, and tyrosine kinase 2 (TYK 2)) and seven STATs (STAT 1, STAT 2, STAT 3, STAT 4, STAT 5a, STAT 5b, STAT 6) proteins (5). JAK-STAT proteins play a key role in regulating lymphoid homeostasis and immunity. This includes the maintenance of the balance between the development of T helper 1 (Th1) and T helper 2 (Th2), T-cell response, regulatory T (Treg) cells, and the function of memory CD8+ cells in myeloid and lymphoid development (5). JAK-STAT signaling is perversely activated in lymphoma by multiple mechanisms, together with inappropriate autocrine and paracrine cytokine stimulation. Genetic and epigenetic changes of negative regulators of JAK/STAT signaling, such as loss of function, SOCS1 mutations, and deletions of the protein tyrosine phosphatase 2 (PTPN2), may also cause  deregulated JAK/STAT initiation in Hodgkin lymphoma, primary mediastinal large B-cell lymphoma, diffuse large B-cell lymphoma, follicular lymphoma, and Peripheral T-cell lymphoma (13). As mentioned, one of the foremost pathways taking part in the signal transductions of a wide array of these cytokines is the JAK/STAT cascade. The JAK/ STAT pathway is crucial for many significant biological processes, including broad immune and hematopoietic cell functions (1). Acquired mutations may influence the JAK/STAT pathway by activating members of the JAK and STAT families directly, inactivating proteins whose typical function is to deactivate the JAKs, or establishing autocrine signaling loops that drive JAK-mediated multiplication. Some of the mutations that are found to be linked to lymphoma development directly target genes that encode elements of the JAK/STAT pathway. JAK 1 and STAT 3 mutations show this relationship between lymphoma development and the JAK/STAT pathway (2). However, even in the presence of STAT or JAK mutations, the whole functional cytokine receptors, JAK and STAT, were required to maintain activation and malignant cell proliferation (17). STAT 3 and STAT 5b mutations are present in aggressive lymphomas emerging from natural killer cells as well as the mutated STATs that are related to increased tyrosine phosphorylation. They supply a growth benefit to natural killer cells, which can be partially inhibited by a JAK 1 and JAK 2 inhibitor (4). Preliminary studies considering the efficacy of JAK inhibitors illustrate that therapeutic agents targeting the JAK/STAT signaling pathway can be used to treat patients with lymphoma (3). On the other hand, ruxolitinibinduced JAK/STAT pathway inhibition in myelofibrosis is associated with an elevated frequency of aggressive B-cell lymphomas (6). Based on the results of our present in silico study, there must be some concern for the development of lymphoproliferative neoplastic disease in a given patient under ruxolitinib since the drug affects numerous genes that have a clear impact on the pathobiology of lymphomas. Table 3 summarizes the list of genes present in lymphoid cancer stem cells that were affected by ruxolitinib usage and their clinical correlations. Thus, ruxolitinib may potentially lead to the pathological expression of the transcription factors important in lymphomagenesis, neoplastic commitment on the progenitor lymphoid cells, inhibition of repressor transcriptions protective for lymphoma development, inhibition of apoptosis, promotion of neoplastic proliferation, transcriptional activation, and proliferation of malignant neoplastic B cells (Table 3). Hence, we fully agree with Porpaczy et al. that the detection of a preexisting B-cell clone may identify individuals at risk for lymphoma development and any ruxolitinib candidate should undergo a bone marrow biopsy procedure for the evaluation of clonal lymphoid baseline proliferation (6). This suggestion is particularly helpful for early stage myelofibrosis or PV since those patients potentially have decades for survival (7,18,19). Ruxolitinib is administered for PV just for the control of symptoms, but if lymphoma due to ruxolitinib complicates the clinical picture, the survival could be decreased from about two decades to several months in a given PV patient. Ruxolitinib is a JAK/STAT signaling pathway inhibitor targeted drug with predictable pharmacobiological actions. The main function of the JAK/STAT signaling pathway is cellular proliferation in health and disease. The drug has been approved for the treatment of patients with high-or intermediate-risk myelofibrosis (MF) with symptomatic splenomegaly. The development of lymphoma due to ruxolitinib might be acceptable for a high-risk advanced stage of MF. However, based on the results of our current study, the administration of ruxolitinib may not be rational for symptom control only in early-stage MF, low-risk MF, and PV patients for the potential pathobiological risk of lymphoma development.
The clinical relevance of the novel results of our study adds support to the concerns of Porpaczy et al. (6). If pharmacobiological aspects of ruxolitinib are related to the danger of the development of lymphoma in patients with chronic myeloproliferative disorders, then decades of usage of the drug could be harmful at least for the patient subpopulation with already present preneoplastic bone marrow lymphoid follicles. Clinicians dealing with the management of chronic myeloproliferative disorders should be aware of this fact and should clinicopathologically check lymphoid neoplastic disorders before and during the administration of ruxolitinib. Of course, myeloid and lymphoid neoplasia are different diseases, but long-term administration of JAK/STAT inhibitors may be located at a crossroads to complicate lymphoma in a patient with chronic myeloproliferative disease.
The development of any drug from bench side to the clinic is very difficult, painful, and expensive. Sometimes unexpected pharmacobiological adverse effects could complicate the treatment strategy with novel drugs. Therefore, a well-developed scientific strategy is absolutely necessary for the design of clinical studies at all critical points. Ruxolitinib is the first clinically useful targeted therapy in Ph*-negative MPNs. In the interest of a clinically advanced approach to the study of this drug, the progression of lymphoma development associated with ruxolitinib in addition to the disease risk profile and scoring of MPNs should be included. Table 7. The resulting pathways from Reactome analysis and the number of genes for each of them ordered by P-value.