Effect of HPSE and HPSE2 SNPs on the Risk of Developing Primary Paraskeletal Multiple Myeloma

Multiple myeloma (MM) is a plasma cell malignancy that is accompanied by hypercalcemia, renal failure, anemia, and lytic bone lesions. Heparanase (HPSE) plays an important role in supporting and promoting myeloma progression, maintenance of plasma cell stemness, and resistance to therapy. Previous studies identified functional single nucleotide polymorphisms (SNPs) located in the HPSE gene. In the present study, 5 functional HPSE SNPs and 11 novel HPSE2 SNPs were examined. A very significant association between two enhancer (rs4693608 and rs4693084), and two insulator (rs4364254 and rs4426765) HPSE SNPs and primary paraskeletal disease (PS) was observed. SNP rs657442, located in intron 9 of the HPSE2 gene, revealed a significant protective association with primary paraskeletal disease and lytic bone lesions. The present study demonstrates a promoting (HPSE gene) and protective (HPSE2 gene) role of gene regulatory elements in the development of paraskeletal disease and bone morbidity. The effect of signal discrepancy between myeloma cells and normal cells of the tumor microenvironment is proposed as a mechanism for the involvement of heparanase in primary PS. We suggest that an increase in heparanase-2 expression can lead to effective suppression of heparanase activity in multiple myeloma accompanied by extramedullary and osteolytic bone disease.


Introduction
Multiple myeloma (MM) is a malignant neoplasm of plasma cells. The disease is defined by the presence of 10% or more clonal plasma cells in the bone marrow (BM), accompanied by one or more myeloma-defining events (MDE) consisting of the CRAB criteria (hypercalcemia, renal failure, anemia, and lytic bone lesions), and three specific biomarkers: clonal bone marrow plasma cells ≥ 60%, serum free light chain (FLC) ratio ≥ 100, and more than one focal lesion in MRI [1,2]. A subclone of plasma cells, capable of growing outside the BM, may lead to extramedullary multiple myeloma (EMD) [3]. Extramedullary myeloma, presented at initial diagnosis (primary EMD) or recurrence (secondary EMD), includes two categories: i. Extramedullary disease (EMD), an aggressive form of MM with soft-tissue plasmacytomas resulting from hematogenous spread; ii. Paraskeletal disease (PS) characterized by the presence of soft-tissue plasmacytomas arising from direct growth of skeletal tumors with distraction of the bone cortex. Solitary plasmacytoma (SP) and plasma cell leukemia (PCL) are typically excluded from the definition of EMD [3,4]. Despite the implementation of new therapeutic agents, MM is still an incurable disease due to recurrent relapses in most patients [5]. of which are rs4693608 (enhancer) and rs4364254 (insulator) SNPs [30,31]. Combinations of rs4693608 and rs4364254 SNPs are significantly associated with high (HR), intermediate (MR) and low (LR) HPSE gene levels in normal leukocytes and mononuclear cells (MNCs) in both peripheral and umbilical cord blood [32,33]. The HR group includes all individuals with the rs4693608 AA genotype. The MR group consists of carriers AG-TT and AG-TC genotypes, and the LR group includes the genotypes GG-TT, GG-TC, GG-CC, and AG-CC. Moreover, a highly significant association was found between rs4693608 and the risk of acute and extensive chronic graft versus host disease (GVHD) after allogeneic stem cell transplantation (HSCT). Importantly, the discrepancy between the recipient and donor in these SNPs significantly affects the risk of acute GVHD [33,34]. Notably, the enhancer rs4693084 SNP and insulator rs28649799 SNP were shown to be associated with the yield of G-CSF-mediated CD34 + cell mobilization in normal donors [31,35]. In addition, rs4426765 is associated with heparanase gene expression in activated MNCs, as well as CD3 cell number and lymphocyte count in response to G-CSF-induced cell mobilization [31].
In the present study, the functional SNPs of the HPSE and HPSE2 genes were analyzed. For the first time, 11 SNPs of HPSE2 were examined. Our results show a very significant association between two enhancer (rs4693608 and rs4693084), and two insulator (rs4364254 and rs4426765) HPSE SNPs and primary paraskeletal disease (PS). SNPs rs4693608 and rs4364254 exhibited the most significant results. The frequencies of the LR genotype, previously shown to be associated with a low level of HPSE expression and a low risk of developing acute GVHD [32][33][34], were significantly higher in the group with total and paraskeletal disease compared with control patients. Significant correlations were found with rs4693608, rs4693084, and rs4426765 SNPs and multiple myeloma stage at diagnosis. SNP rs657442, located in intron 9 of the HPSE2 gene, revealed a significant protective association with primary PS disease and lytic bone lesions. The impact of rs658225 and rs754585 HPSE2 gene SNPs was noted in secondary extramedullary disease. The present study demonstrates the promoting (HPSE) and protective (HPSE2) role of regulatory gene elements in the development of paraskeletal disease and bone morbidity. We suggest that signal discrepancy between malignant myeloma cells and normal cells of the microenvironment provides a possible mechanism for the involvement of heparanase in primary EMD.

Study Population
Two hundred and seventy eight patients were included in the study. Nineteen patients were diagnosed with monoclonal gammopathy (MGUS) and 28 with smoldering multiple myeloma (SMM). Six patients with MGUS and 4 with SMM progressed to active multiple myeloma (MM). Clinical and disease characteristics, including disease stage by the International Staging System (ISS), bone involvement, renal failure, extramedullary disease, beta-2-microglobulin, creatinine, albumin, calcium, and hemoglobin levels were collected from medical records. The definition of renal failure is based on an estimated glomerular filtration rate (eGFR) or creatinine clearance (CrCl) of 29-15 mL/min/1.73 m 2 . Patient characteristics are presented in Table 1. The study was approved by the Sheba Medical Center Ethics Committee and the Israeli Ministry of Health (#4247) and all subjects gave their written informed consent.

SNPs Analysis
Genomic DNA was extracted from peripheral blood applying the Wizard ® Genomic DNA Purification kit (Promega, Madison, WI, USA). Functional SNPs of the HPSE and HPSE2 genes were genotyped using Real-Time SNP Assay (Bio Search Technologies, Novato, CA, USA) and TaqMan SNP Genotyping Assay (Thermo Fisher Scientific, Pleasanton, CA, USA). PCR reactions were performed using an ABI PRISM 7700 sequence detector (Applied Biosystems, Warrington, UK) according to the manufacturer's instructions (Quanta, Gaithersburg, MD, USA).

DNA Constructs
A fragment of the HPSE gene enhancer (440 bp), located in intron 2 and including rs4693083, rs4693608, rs4693084 and rs4693609 SNPs, was cloned via PCR using PCR II-TOPO vector (Invitrogen by Life Technologies, Carlsbad, CA), as previously described in detail [30]. The multiple myeloma cell line RPMI-8226, isolated from a patient with plasmacytoma, was seeded into six-well culture dishes according to Lonza instructions (AMAXA biosystems, Lonza, Germany). Cells were transfected with 2 µg of each DNA construct using the Ingenio Electroporation Kit (Mirus Bio, Madison, WI, USA) and Nucleofector TM (AMAXA biosystems, Lonza, Germany). Twenty-four hours after electroporation, the cells proceeded to RNA analysis.

Purification of Total RNA, Generation of cDNA, and Real-Time Quantitative RT-PCR
Total RNA was obtained from RPMI-8226 cell line using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Complementary DNA (cDNA) was synthesized using 1µg total RNA applying the qScript cDNA Synthesis kit (Quanta, Gaithersburg, MD, USA). PerfeCTa SYBR Green Fast Mix ROX (Quanta, Gaithersburg, MD, USA) Real-time PCR reactions to amplify the HPSE and HPSE2 genes and the RPL-19 housekeeping gene were performed using an ABI PRISM 7700 sequence detector (Applied Biosystems, Warrington, UK) in a total volume of 20 µL. The primer sequences for HPSE and RPL-19 genes were previously published [36]. HPSE2 expression was determined using primers 5'CAGGGCATTGATGTCGTGATAC3' and 5'GCCAGTAGTCTG-GTAATGGGTTA3'. The experiments were performed in triplicates and a standard curve was created. A dissociation curve was generated after each experiment to confirm that one product was amplified. Comparative CT method (∆∆CT) was used to quantify HPSE and HPSE2 gene expression, and relative quantification (RQ) was calculated as 2 −∆∆CT .

Statistical Analysis
All calculations were performed using the NCSS software (NCSS, Kaysville, Utah, USA). A p-value ≤ 0.05 was considered statistically significant. Genotype and allele frequencies of the SNPs were calculated by direct counting. Possible differences between samples in the frequency of each SNP genotype or allele were assessed using the χ 2 test. The Yates correction was used when analyzing the interaction of two or more SNPs. Due to the low frequency of allele A of SNP rs657422, association analysis was carried out using Fisher's exact test.

Analysis of HPSE and HPSE2 Gene SNPs among Multiple Myeloma Patients
The retrospective study included 278 multiple myeloma patients (174 males and 104 females) and 211 controls. Two hundred and thirty one patients were diagnosed with multiple myeloma, 28 with SMM and 19 with MGUS. The median age was 60 years (range 27-90). All statistical analyzes, except for association with the patient's diagnosis (MM, SMM and MGUS), were performed in 259 patients who developed active MM. Patient characteristics are presented in Table 1. Six HPSE gene SNPs were used in the case vs control analyses. Two functional SNPs rs4693608 and rs4693084 are located in the HPSE gene enhancer, three functional SNPs (rs4426765, rs28649799 and rs4364254) are mapped in the HPSE gene insulator, and rs4693602 is situated in the 3'UTR ( Figure 1). SNPs of the HPSE2 gene were analyzed for the first time. Eleven SNPs that could potentially be functional polymorphisms or be associated with functional SNPs were selected ( Figure 1). A correlation analysis was performed between the HPSE and HPSE2 SNPs and the patient's diagnosis (MM, SMM, and MGUS). The frequency of the TT genotype (rs4364254) was slightly higher in patients with MM compared to the control group of healthy individuals, but statistical analysis showed only a trend toward correlation (Supplementary Table S1). No associations were observed for rs4364254 allele frequencies (Supplementary Table S1). In addition, no correlations were found between the HPSE2 gene SNPs and the patient's diagnosis.
control analyses. Two functional SNPs rs4693608 and rs4693084 are located in the HPSE gene enhancer, three functional SNPs (rs4426765, rs28649799 and rs4364254) are mapped in the HPSE gene insulator, and rs4693602 is situated in the 3'UTR ( Figure 1). SNPs of the HPSE2 gene were analyzed for the first time. Eleven SNPs that could potentially be functional polymorphisms or be associated with functional SNPs were selected ( Figure 1). A correlation analysis was performed between the HPSE and HPSE2 SNPs and the patient's diagnosis (MM, SMM, and MGUS). The frequency of the TT genotype (rs4364254) was slightly higher in patients with MM compared to the control group of healthy individuals, but statistical analysis showed only a trend toward correlation (Supplementary Table S1). No associations were observed for rs4364254 allele frequencies (Supplementary Table S1). In addition, no correlations were found between the HPSE2 gene SNPs and the patient's diagnosis.

Association between HPSE Gene SNPs and Stage of Multiple Myeloma Patients at Diagnosis
The influence of HPSE and HPSE2 SNPs on the stage of multiple myeloma at diagnosis was analyzed in 151 patients. The stage was determined according to the ISS [37]. Stage I was diagnosed in 25% of patients, stage II in 43.7%, and stage III in 30.5%. The association between HPSE gene SNPs and MM staging is shown in Table 2. Significant correlations were found with enhancer rs4693608 and rs4693084 SNPs and insulator rs4426765 SNP. The most significant results were obtained in comparing stage I and stage II/ III. The frequencies of the AA genotype were higher in patients diagnosed with stage II and III (15.4% in stage I vs. 40.9% and 32.6% in stages II and III, respectively). The frequencies of the GG genotype were higher in patients diagnosed with stage I vs II and III (25.6% at stage I vs. 12% and 15.2% at stages II and III, respectively). The p-value for genotype frequency was p = 0.022 and the p-value for allele frequency was p = 0.0082 (Table  2). No correlations with HPSE2 gene SNPs were found.

Association between HPSE Gene SNPs and Stage of Multiple Myeloma Patients at Diagnosis
The influence of HPSE and HPSE2 SNPs on the stage of multiple myeloma at diagnosis was analyzed in 151 patients. The stage was determined according to the ISS [37]. Stage I was diagnosed in 25% of patients, stage II in 43.7%, and stage III in 30.5%. The association between HPSE gene SNPs and MM staging is shown in Table 2. Significant correlations were found with enhancer rs4693608 and rs4693084 SNPs and insulator rs4426765 SNP. The most significant results were obtained in comparing stage I and stage II/ III. The frequencies of the AA genotype were higher in patients diagnosed with stage II and III (15.4% in stage I vs. 40.9% and 32.6% in stages II and III, respectively). The frequencies of the GG genotype were higher in patients diagnosed with stage I vs II and III (25.6% at stage I vs. 12% and 15.2% at stages II and III, respectively). The p-value for genotype frequency was p = 0.022 and the p-value for allele frequency was p = 0.0082 (Table 2). No correlations with HPSE2 gene SNPs were found.

Association between HPSE Gene SNPs and Primary Bone-Related Disease in Multiple Myeloma
Fifty patients with primary extramedullary disease were identified. Forty-one of 50 were diagnosed with paraskeletal disease, one with EMD from hematogenous spread and two had both bone-depended and bone-independent EMD. One patient was diagnosed with solitary plasmacytoma, and five developed plasma cell leukemia. Table 3. Calculations were performed for all individuals diagnosed with any type of EMD and separately for patients who developed paraskeletal disease. Four functional SNPs revealed highly significant associations ( Table 3). The frequency of the GG genotype and the G allele of rs4693608 SNP was higher in patients with EMD compared to the control group (32% versus 13% for the GG genotype, 18% versus 36.4% for the AA genotype, p = 0.003, and 57% versus 38.3%, p = 0.001 for the G allele, respectively, Table 3). Similar results were obtained for paraskeletal disease (p = 0.018 for genotype comparison and p = 0.006 for allele comparison) ( Table 3). All three insulator SNPs were associated with EMD, rs4364254 being the most prominent SNP ( Table 3). The frequency of the CC genotype and the C allele of rs4364254 SNP was higher in patients with primary total EMD and patients with paraskeletal disease (20% versus 8.6%, p = 0.003 for the CC genotype; 43% versus 24.4%, p = 0.00032 for allele C in total EMD; 18.6% versus 8.6%, p = 0.019 for genotype CC; and 40.7% versus 24.4%, p = 0.003 for allele C in paraskeletal disease, respectively, Table 3). In our previous studies, we investigated the interaction between enhancer and insulator HPSE gene SNPs [31,32,34]. We allocated all possible HPSE genotype combinations into three groups (HR, MR, and LR) correlating with high, intermediate and low heparanase mRNA expression levels and with high, intermediate and low risk of acute GVHD, respectively [31,32,34]. Given the results presented in Table 3, we can conclude that rs4693608 and rs4364254 exhibited the most significant association with EMD. Therefore, correlations were analyzed in HR (genotypes AA-TT and AA-TC), MR (genotypes AG-TT and AG-TC), and LR (genotypes GG-TT, GG-TC, GG-CC and AG-CC) groups. Frequency of the LR genotype was significantly higher in patients with total and paraskeletal disease compared with patients without EMD disease (44% for total EMD and 41.9% for paraskeletal disease versus 16.1% in the other patients, χ 2 = 16.036, p = 0.00033 for total EMD and χ 2 = 14.32, p = 0.002 for paraskeletal disease after Yates' correction) ( Table 3).

Association between HPSE2 Gene SNPs and Primary Extramedullary Disease (EMD) in Multiple Myeloma
Eleven HPSE2 SNPs were analyzed. Only one SNP, rs657442, located in intron 9, demonstrated a significant correlation with primary EMD (Table 4). Of note, all patients (100%) diagnosed with EMD were carriers of the CC genotype, while in the control group, only 86.4% were carriers of this genotype (p = 0.0065 for total EMD and p = 0.011 for paraskeletal disease). The A allele frequency was 7.1% in the control group and null in the EMD group (p = 0.038 for total EMD and p = 0.0066 for paraskeletal disease, Table 4).  The interactions between HPSE genotypes and the significant HPSE2 SNP rs657442 are presented in Table 5. A protective effect of the A allele was observed. Five patients with the LR genotype, identified in association with the risk of EMD, were carriers of the A allele for SNP rs657442 and had neither primary nor secondary EMD (Table 5). Primary EMD was not observed in 8 carriers of the A allele and MR genotype. Three of them developed secondary EMD, one with hematogenous spread and two with paraskeletal disease. The frequency of the A allele is low, and there are many variations with null frequencies (Table 5). Therefore, the impact of HPSE-HPSE2 gene interaction was assessed by allocation of all possible variants into three groups with high-, intermediate-and low-risk. In group A, we hypothesized that the protective effect of the A allele could downgrade the risk of developing EMD from high to intermediate, and from intermediate to low. Therefore, group A1 included the HR-CC, HR-CA, HR-AA, and MR-CA genotypes. The LR-CA variant was added to group A2. Since we could not assess the magnitude of the A allele-mediated protective effect against EMD in group B, we assumed that the protective effect was high in all subjects with the A allele. As a result, the LR-CA and MR-CA genotypes were included in group B1. Chi-square analysis was performed for both groups, and after Yates' correction, the summarized effect was found to be highly significant compared to the p-values for each SNP alone (p = 0.000015 for group A and p = 0.0000098 for group B) ( Table 5).

Association between HPSE2 Gene SNPs and Secondary Extramedullary Disease (EMD) in Multiple Myeloma
Thirty-eight patients developed secondary EMD. Seven of them had soft-tissue plasmacytomas with hematogenous spread, 27 had bone-related EMD, and 4 had both types of EMD. No correlation with HPSE gene SNPs was observed. Analysis of the association between HPSE2 gene SNPs and secondary EMD is shown in Table 6. Some associations were identified with rs658225 (intron 3), rs621644 (intron 9), and rs754585 (intron 11) ( Figure 1 and Table 6). Notably, the G allele of rs658225 SNP exhibited a protective effect. The G allele frequency in the control group (18%) was higher compared to patients with secondary EMD (8.1%) (p = 0.03, Table 6). In addition, the TT genotype of rs754585 SNP disclosed a protective effect (33.3% in the control group vs. 13.5% in the secondary EMD group, p = 0.027). However, the frequency of the CC genotype was the same in both groups (Table 6). An opposite effect was observed in heterozygous individuals (64.9% in the secondary EMD group vs. 41.7% in the control group).
The interaction between these SNPs was assessed (Table 7). Since the number of patients with secondary EMD was relatively small, the analysis was carried out not only in the group of secondary EMD but also for all patients who revealed EMD. The same trend was observed both in total and secondary EMD compared with the control group. The protective effect of the G allele was manifested in the possessors of the AG genotype for rs658225 (Table 7).   The influence of SNP rs754585 was demonstrated in carriers of the AA genotype. The risk of secondary EMD in patients with the AA-TT genotype was low compared with the control group (10.5% vs. 22.5%), while in individuals with the AA-CT genotype the risk of secondary EMD was higher (52.6% in the secondary EMD group vs. 23.3% in the control group). No differences were found among AA-CC carriers. Due to the small number of individuals with the AG-TT genotype, it was impossible to confirm an addictive apologetic effect of rs658225 and rs754585 HPSE2 SNPs (Table 7). Additional studies are needed in a larger cohort of patients with EMD.

Association between HPSE2 Gene SNPs and Bone Morbidity in Multiple Myeloma Patients
One hundred and fifty-eight multiple myeloma patients (74.9%) exhibited one or more osteolytic bone lesions, while 53 patients (25.1%) had no lytic bone disease. Analysis of the functional HPSE gene SNPs did not reveal a correlation with bone involvement. Eleven SNPs of the HPSE2 gene were genotyped among patients with multiple myeloma. A highly significant correlation with bone morbidity was observed with rs657442 SNP (intron 9), a SNP with a relatively rare frequency of the A allele. Only one individual was genotyped as homozygous for the A allele. Therefore, the following comparison is related to the frequencies of the CA genotype and the A allele. Only 5.7% of patients with the CA genotype had bone involvement compared to 20.8% of patients who had no detectable bone disease (p = 0.001 for genotype frequency, and p = 0.00016 for allele frequency, Table 8). Collectively, HPSE2 gene SNP rs657442 revealed a protective effect.

Modification of HPSE and HPSE2 Gene Expression after Transient Transfection of RPMI-8226 MM Cells with a Fragment of the Enhancer
In a previously published study, we found an active HPSE enhancer in intron 2 of the HPSE gene [30]. To determine the functional effects of enhancer SNPs, we applied a luciferase reporter gene with a minimal promoter to measure enhancer activity. Briefly, RPMI-8226 cells were transiently transfected with DNA constructs in the sense and antisense directions or with an empty vector [30]. As an additional control, RPMI-8226 cells were subjected to an electric pulse. The HPSE gene enhancer fragment (440 bp) included rs4693083, rs4693608, rs4693084, and rs4693609 SNPs. We have previously reported [30] that the 440 bp fragment of the HPSE gene exhibits enhancer activity in both the sense and antisense directions. In the present study, we evaluated the effect of the enhancer fragment transfection on the expression of HPSE and HPSE2 and the results are presented in Figure 2. HPSE expression was increased in all the transfected samples, while the level of the HPSE gene remained unchanged in the control samples.  [30]. DNA constructs were prepared in sense and antisense direction. The HPSE gene enhancer fragment (440 bp) included rs4693083, rs4693608, rs4693084, and rs4693609 SNPs. Vr-A construct contains allele A of rs4693608 and allele G of rs4693084. Vr-B fragment includes allele G of rs4693608 and allele T of rs4693084. Samples of RPMI-8226 cells before and after transfection were analyzed by quantitative RT-PCR. RQ (relative quantification) values are shown on the y-axis. HPSE gene expression levels were elevated in enhancer fragment-transfected samples compared to controls, while HPSE2 gene expression was decreased to undetectable levels in the same samples. Samples treated with AMAXA electrical pulse alone or transfected with empty plasmid were used as controls. No differences were observed between untreated and treated samples. Data shown are of representative experiments performed in triplicates. Differences in RQ among various constructs were evaluated by t-test. p value of ≤ 0.05 was regarded as statistically significant. * p value of ≤ 0.05. ND-not detected.
The expression level of the HPSE2 gene is low in most malignant cell lines [19,20,21,22,23,24,25,26,27,28,29]. Likewise, RMPI-8226 cells expressed a very low, but detectable level of HPSE2 and there was no effect to electrical pulse and transfection with an empty plasmid. Transfection of the HPSE enhancer fragment resulted in an undetectable level of HPSE2 expression (Figure 2), suggesting that the HPSE gene enhancer not only regulates HPSE gene expression but also has an opposite effect on the expression of the HPSE2 gene. The undetectable level of the HPSE2 gene makes it impossible to evaluate the role of specific HPSE enhancer functional SNPs in co-regulation of both genes.

Discussion
Association studies open up great prospects for elucidating the genetic basis of diseases. Moreover, functional SNP-based studies are effective because they focus on SNPs with the highest prior probability of being associated [38]. Functional SNPs located  [30]. DNA constructs were prepared in sense and antisense direction. The HPSE gene enhancer fragment (440 bp) included rs4693083, rs4693608, rs4693084, and rs4693609 SNPs. Vr-A construct contains allele A of rs4693608 and allele G of rs4693084. Vr-B fragment includes allele G of rs4693608 and allele T of rs4693084. Samples of RPMI-8226 cells before and after transfection were analyzed by quantitative RT-PCR. RQ (relative quantification) values are shown on the y-axis. HPSE gene expression levels were elevated in enhancer fragment-transfected samples compared to controls, while HPSE2 gene expression was decreased to undetectable levels in the same samples. Samples treated with AMAXA electrical pulse alone or transfected with empty plasmid were used as controls. No differences were observed between untreated and treated samples. Data shown are of representative experiments performed in triplicates. Differences in RQ among various constructs were evaluated by t-test. p value of ≤0.05 was regarded as statistically significant. * p value of ≤0.05. ND-not detected.
The expression level of the HPSE2 gene is low in most malignant cell lines [19][20][21][22][23][24][25][26][27][28][29]. Likewise, RMPI-8226 cells expressed a very low, but detectable level of HPSE2 and there was no effect to electrical pulse and transfection with an empty plasmid. Transfection of the HPSE enhancer fragment resulted in an undetectable level of HPSE2 expression (Figure 2), suggesting that the HPSE gene enhancer not only regulates HPSE gene expression but also has an opposite effect on the expression of the HPSE2 gene. The undetectable level of the HPSE2 gene makes it impossible to evaluate the role of specific HPSE enhancer functional SNPs in co-regulation of both genes.

Discussion
Association studies open up great prospects for elucidating the genetic basis of diseases. Moreover, functional SNP-based studies are effective because they focus on SNPs with the highest prior probability of being associated [38]. Functional SNPs located in regulatory elements, such as enhancers and insulators, help not only to indicate the involvement of a particular gene in the development of common diseases but also to explore the possible mechanisms of these enhancers or insulators in pathological processes [39,40]. Moreover, these regions can be used as targets for subsequent treatments [41].
In the present study, five previously identified and characterized functional HPSE gene SNPs were analyzed among patients with multiple myeloma. A highly significant association was found between the low-frequency alleles of rs4693608, rs4693084, rs4426765 and rs4364254 SNPs and an increased risk of primary paraskeletal disease. A previously published approach to studying the interaction between HPSE enhancer and insulator SNPs was applied [31,32,34]. The LR genotype, which correlates with low expression of the HPSE gene, was found to be associated with primary EMD.
Heparanase-2 appears to function as a natural inhibitor of the heparanase enzyme and a tumor suppressor [19,24,42]. While heparanase enhances the pathogenesis of cancer, heparanase-2 plays a protective role. It appears that the heparanase/Hpa2 ratio impacts tissue hemostasis [43,44] and the balance between cancer progression and suppression [19,24,42].
The present study aimed to identify and characterize possible functional SNPs located in the HPSE2 gene. Importantly, the rs657442 SNP, genotyped at intron 9, showed a protective effect on primary paraskeletal disease and bone-related morbidity. According to the EN-CODE Registry of candidate cis-Regulatory Elements (cCREs), rs657442 SNP is located in a region with a distal enhancer-like signature (https://genome.ucsc.edu/cgi-bin/hgTracks?d b=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&v irtMode=0&nonVirtPosition=&position=chr10%3A98567418%2D98567918&hgsid=150243 2943_OAIfBAS927s66nGeaamrwOHsFwZi, accessed on 20 May 2020). Additional three HPSE2 SNPs (rs658225, rs621644 and rs754585) revealed some correlation with secondary EMD. The rs658225 polymorphism was mapped to intron 3 and, according to the ENCODE Registry, is localized to the CTCF-related insulator (https://genome.ucsc.edu/cgi-bin/hgTr acks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=def ault&virtMode=0&nonVirtPosition=&position=chr18%3A43420931%2D43421131&hgsid= 1502451125_iiE7HfeK7puYOrj2ozO3noTl2oPi, accessed on 20 May 2020). Intron 3 of the HPSE2 gene is a large intron that contains two pseudogenes and a gene encoding for MIR6507. It can be assumed that SNPs located in intron 3 may play an important role in the regulation of HPSE2 and other processes associated with MIR6507, and the RPL7P36 and ARL5AP2 pseudogenes. The rs624644 SNP has also been mapped to a region with a distal enhancer-like signature and could potentially be functional. SNP rs754585 was found to be in complete linkage disequilibrium with many polymorphic SNPs that are located in the 3'UTR region that can influence the expression and function of the HPSE2 gene.
Heparanase is one of the main regulators of aggressive tumor behavior [7]. It drives angiogenesis, tumor growth, and metastasis [9,10,14,[45][46][47][48]. Soluble heparanase is detected in the bone marrow of patients with myeloma associated with a high microvessel density [9]. Our results indicate that patients with genotypes AA-AA and AG-AC (rs4693608-rs4426765) had an increased risk of being diagnosed with stage II and III of the disease. Our previous studies revealed that the A allele is associated with higher heparanase levels in healthy individuals [32,33] and a higher risk of developing GVHD [34]. Andersen et al. identified a significant correlation between the A allele of rs4693608 and the risk of vertebral fractures at diagnosis [49]. Taken together, it is conceivable that the HPSE gene enhancer plays an important role in progression of the disease before diagnosis, and that the risk of multiple myeloma being diagnosed at later stages of the disease is higher in individuals with the AA genotype.
The interaction between myeloma cells and the bone marrow microenvironment promotes signaling cascades and activates chemotaxis and adhesion of myeloma cells to the bone marrow [50]. MM progression and extramedullary dissemination require continuous MM cell egress and spread into new extramedullary sites, regulated by dynamic changes in the expression of chemokine receptors and adhesion molecules. Thus, increased CXCR4/CXCL12 chemokine signaling in MM cells promotes epithelial to mesenchymal transition (EMT)-like phenotype with enhanced invasive properties [50][51][52][53][54]. Consequently, the CXCR4-mediated EMT-like program involves E-cadherin downregulation, resulting in decreased adhesion and enhanced egression of MM cells [54]. Furthermore, down-regulation of adhesion molecules (VLA-4, p-selectin, CD56, CD44), and membraneembedded tetraspanins as well as up-regulation of heparanase are part of the mechanism regulating extramedullary spread in MM [3,18,53].
Heparanase is a multitasking protein whose activation promotes, among other effects, the release of various cytokines, chemokines and growth factors in the tumor microenvironment [9,55,56]. Moreover, nuclear heparanase in activated human T-lymphocytes modifies the activity of inducible immune response genes by association with euchromatin and regulation of histone H3 methylation [57]. ChIP-on-chip data revealed that heparanase bound to regulatory regions of myeloma-related genes (i.e., CXCL12) [57], encodes for SDF-1, CXCR4, CXCL2, CXCR3, CCL27, CX3CL1, and LCN2 [58][59][60]. The HPSE2 gene is also on this list [57], indicating that the interaction between HPSE and HPSE2 may occur at both the gene and protein levels [23]. This was confirmed our results. Transfection of HPSE gene enhancer to myeloma cells resulted in increased expression of the HPSE gene and subsequent downregulation of the HPSE 2 gene.
Unexpectedly, the present results indicate that MM patients with the LR genotype disclosed an increased risk of developing primary paraskeletal disease. Studies performed over the past two decades have shown that heparanase plays a major role in modulating the bone marrow microenvironment to support the progression of multiple myeloma [6]. High heparanase activity in myeloma cells correlates with myeloma growth, bone marrow angiogenesis, and osteolytic bone disease [6,9,11,12]. Our previous results showed that transfection of the HPSE enhancer and insulator fragments into RPMI-8226 and CAG human myeloma cells led to the highest level of luciferase activity compared to other similarly transfected cell types [30,31]. Thus, one would expect that the risk of developing a more aggressive disease in patients with HR genotype would be higher in comparison to carriers of the LR genotype. Indeed, a significant correlation was found between the rs4693608 AA genotype, which constitutes the bulk of the HR genotype, and a more progressive disease (stages II and III) at the time of diagnosis.
To explain how the LR genotype increases the risk of paraskeletal disease, we propose the following mechanism, which is summarized in Figure 3. In our previous published studies [33,34], we found that the effect of a discrepancy between recipient and donor in HPSE gene SNPs leads to an increased risk of developing acute GVHD. The same effect was now observed in primary paraskeletal disease. To simplify, we adopt the concept of the "fight-or-flight" response to stress, first described by Walter Cannon [61]. Animals respond to stress with a set of predictable reactions known as "fight-or-flight" [62]. The function of the autonomic nervous system (ANS) in mediating the "fight-or-flight" response is well studied [63]. Briefly, the ANS is a major modulator of immune response and immune effector cells by activation of pro-and anti-inflammatory functions [64]. In the hematopoietic system, the ANS regulates stem cell niche homeostasis, and fine-tunes the inflammatory response [64]. Based on this approach, we named the HR genotype "fight" and the LR genotype "flight". Our previous studies revealed that the HR genotype exhibited high HPSE mRNA expression, low plasma heparanase level, significant HPSE up-regulation, and most importantly, a higher index threshold after LPS treatment [32,33]. In contrast, the LR genotype disclosed low HPSE gene expression level, high plasma heparanase level, minor HPSE up-regulation, and no effect or downregulation of HPSE in response to LPS treatment [32,33]. A lower index threshold was observed ( Figure 3A). Downstream heparanase signaling affects a variety of molecular mechanisms including transcription of immunologically related genes and miRNAs, and release of various cytokines, chemokines, and growth factors from the extracellular matrix [9][10][11]57]. Moreover, heparanase is an important modulator of exosome biogenesis via the syndecan-syntenin-ALIX pathway, affecting the level of exosomal cargo and the regulation of neurotransmission [65]. ALIX pathway, affecting the level of exosomal cargo and the regulation of neurotransmission [65]. We suggest that the regulatory elements of the HPSE gene are involved in cell responsiveness to activating signals. After HSCT, engrafted donor cells (LR genotype), are exposed to the patient cells (HR genotype) and pro-inflammatory environment. When the stimulation threshold is exceeded, such donor cells with the LR "flight" genotype react with recipient cells harboring the HR "fight" genotype, resulting in aberrant activation and increased risk of developing acute GVHD increased ( Figure 3B) [33,34].
We explain the identified risk of developing primary paraskeletal disease in patients with the LR genotype using a previously published discrepancy model [33,34]. Although aggressive myeloma cells exhibited the LR genotype, these cells have acquired a "fight" phenotype with high heparanase activity. In other words, while myeloma cells express HR-like phenotype, the normal microenvironment continues to express the LR genotype, leading to differences in signaling responses. When the signaling stimulation threshold is exceeded, aberrant activation of normal cells in the microenvironment results in bone cortex rupture and paraskeletal disease progression ( Figure 3C). Previous studies have shown that heparanase expression by myeloma cells significantly increases local and systemic osteolysis [11]. In addition, heparanase can also alter the fate of osteoblast progenitor cells, directing cellular differentiation toward adipocytes rather than osteoblasts [66]. We suggest that the regulatory elements of the HPSE gene are involved in cell responsiveness to activating signals. After HSCT, engrafted donor cells (LR genotype), are exposed to the patient cells (HR genotype) and pro-inflammatory environment. When the stimulation threshold is exceeded, such donor cells with the LR "flight" genotype react with recipient cells harboring the HR "fight" genotype, resulting in aberrant activation and increased risk of developing acute GVHD increased ( Figure 3B) [33,34].
We explain the identified risk of developing primary paraskeletal disease in patients with the LR genotype using a previously published discrepancy model [33,34]. Although aggressive myeloma cells exhibited the LR genotype, these cells have acquired a "fight" phenotype with high heparanase activity. In other words, while myeloma cells express HR-like phenotype, the normal microenvironment continues to express the LR genotype, leading to differences in signaling responses. When the signaling stimulation threshold is exceeded, aberrant activation of normal cells in the microenvironment results in bone cortex rupture and paraskeletal disease progression ( Figure 3C). Previous studies have shown that heparanase expression by myeloma cells significantly increases local and systemic osteolysis [11]. In addition, heparanase can also alter the fate of osteoblast progenitor cells, directing cellular differentiation toward adipocytes rather than osteoblasts [66]. The effect of discrepancy in signaling between myeloma cells and the normal microenvironment, and the aberrant activity of cells in the bone marrow niche contribute to osteolysis, bone cortex distraction and various manifestation of bone-associated extramedullary disease.