Analysis of the Relationship between the Expression Levels of Neutrophil Gelatinase-Associated Lipocalin and Cytokine Genes in Bone Marrow

Background: Recently, various associations of NGAL with several hematological cancers have been reported. However, given that the regulation of NGAL gene expression by cytokines is tissue-specific, NGAL expression in relation to those of cytokine genes has not been analyzed in bone marrow (BM) tissue. The purpose of this study was to analyze the association between NGAL and 48 cytokine gene expression levels in mononuclear cells (MNCs) of BM at the time of diagnosis of hematological malignancy and to explore the expression pattern of NGAL and related cytokine genes in patients with hematological malignancies and controls. Methods: BM MNCs were isolated from 48 patients, who were classified as patients presenting myeloproliferative neoplasm, acute myeloid leukemia, myelodysplastic syndrome, and as controls. NGAL and cytokine genes were analyzed using NanoString. Data on hematological parameters were collected from medical records. Single and multiple regression analyses were performed to analyze relationships. Results: Normalized counts of 26 cytokine genes were related to NGAL normalized counts, while STAT3 and TLR4 normalized counts had the highest explanatory power. The following multiple regression model was developed: NGAL normalized counts=4316.825 + 9.056 × STAT3 normalized counts + 844.226 × IL5 normalized counts + 17.540 × TLR1 normalized counts - 28.206 × TLR2 normalized counts - 42.524 × IRAK4 normalized counts. In the multiple regression analysis, STAT3 and TLR4 normalized counts showed multicollinearity. NGAL, STAT3, IL5, and TLR4 normalized counts showed similar intergroup patterns. Conclusions: NGAL normalized counts was predicted by a multiple regression model, while they showed similar intergroup patterns to STAT3, IL5, and TLR4 normalized counts.


Introduction
Neutrophil gelatinase-associated lipocalin (NGAL), a member of the lipocalin superfamily, is a glycoprotein originally isolated from the secondary granules of human neutrophils [1]. While NGAL is useful as a biomarker especially in kidney injury, it is also expressed in other tissues in response to various pathological conditions, such as ischemia, tissue injury, and cancer [1][2][3]. Being a pleiotropic mediator of various inflammatory processes [4,5], NGAL acted at different levels as a regulator of multiple responses [6]. Recently, association of NGAL with several hematologic cancers have been reported [1,[5][6][7][8].
A few studies measuring NGAL protein levels in the bone marrow (BM) of patients with hematological Ivyspring International Publisher cancer reported that NGAL levels were significantly lower in patients with acute myeloid leukemia (AML) and myelodysplastic syndromes (MDS) than in controls [8][9][10]. Given that neutrophilic precursors are the major source of NGAL synthesis and these exist primarily in the BM but not the peripheral blood, such studies have enabled an accurate comparative analysis of NGAL levels according to the hematological malignant disease entity [8].
Nevertheless, NGAL expression has not been analyzed in BM aspirate cells of patients with hematological malignancies. In particular, although NGAL expression is regulated by various types of inflammatory cytokines, it should be noted that the regulation of NGAL expression is tissue-specific [1]. Accordingly, analysis of the relationship between NGAL and inflammatory cytokines in BM cells is important for elucidating the mechanism of the regulation of NGAL in the BM.
The purpose of this study was to identify the expression profile of NGAL and 48 cytokine genes in mononuclear cells (MNCs) of BM aspirate at the time of the diagnosis of a hematological malignancy using the nCounter system [11] and analyze the relationship between NGAL expression and those of 48 cytokine genes. Additionally, we explored the expression patterns of NGAL and related cytokine genes in patients with hematological malignancies and controls.

Sample collection and preparation
The study was approved by the Institutional Review Board of Korea University Ansan Hospital and performed in accordance with the Declaration of Helsinki. Patients were enrolled between May 2018 and July 2019. Informed consent was obtained from all participants (n=48). Aliquots of leftover BM aspirate samples were collected from patients who underwent BM examination for the diagnosis of hematological malignancies. Patients were classified into myeloproliferative neoplasm (MPN), AML, and MDS groups based on the WHO diagnostic criteria, after the BM smear and pathological review. The control group (n=9) consisted of patients with lymphoma without BM involvement (n=7) or normocellular marrow without hematological malignancy (n=2) [8,9]. None of the patients manifested symptoms of active infections, inflammatory diseases, or kidney failure [7].

MNC isolation
Isolation of MNCs from BM aspirates was performed using Lymphoprep in SepMate tubes (STEMCELL Technologies, Vancouver, Canada) [12]. First, 4.5 mL of Lymphoprep (density gradient medium; density: 1.077 g/mL; STEMCELL Technologies) was added to 15 mL of SepMate tubes and pipetted through the central hole of the insert. The BM aspirate sample was diluted with an equal volume of phosphate-buffered saline (PBS) and mixed gently, then pipetted down the side of the (vertically held) tube. The tubes were then centrifuged (1200 × g, 10 min) and the top layer containing the enriched MNCs was poured into a new tube (15-mL conical polypropylene tube; Becton Dickinson), and PBS was added, resulting in a total volume of 14.5 mL. After centrifugation (300 × g, 8 min), the supernatant was carefully removed, and 2 mL of PBS was added to the infranatant (including BM MNCs). The solutions were mixed well, and the mixture was transferred to 1.5-mL Eppendorf tubes and centrifuged (400 × g, 5 min). The supernatant was then removed and the infranatant and isolated BM MNCs were stored at -80 °C until mRNA extraction.

mRNA extraction from BM MNCs
Frozen cell pellets were thawed slowly in a pre-chilled bead bath. Then, 1 mL Trizol (#79306-200mL; Qiagen, Valencia, CA, USA) was added to the thawed cell pellets. The mixture was briefly vortexed and incubated for 5 min at 20-24 °C, and 200 μL chloroform (#34854-1L; Sigma Aldrich, St. Louis, MO, USA) was added. The mixture was shaken for 15 s and incubated for 5 min at 20-24 °C. After centrifugation (13,000 × g, 15 min, 4 °C), the aqueous phase was transferred to new 1.5-mL tubes and an equal volume of isopropanol was added. The tubes were inverted manually four times, incubated at 20-24 °C for 10 min, and centrifuged (13,000 × g, 10 min, 4 °C). The supernatant was discarded. The pellets were washed twice with 1 mL of 70% ethanol and centrifuged (7,500 × g, 5 min, 4 °C). RNA pellets were air-dried. RNA was resuspended in 10-20 μL RNase-free water, transferred to new 1.5-mL tubes, and RNA quality was determined using a 2100 Bioanalyzer capillary electrophoresis system (Agilent Technologies, Santa Clara, CA).
A reaction mixture containing RNA (5 μL; 100-300 ng), Master Mix (8 μL; reporter CodeSet and hybridization buffer), and capture probeSet (2 μL) was prepared, and the solution was mixed and spun down. It was incubated at 65 °C in a thermocycler (Bio-Rad Laboratories Inc., Hercules, CA, USA) for 16 h (maximum hybridization, 48 h). Samples were then transferred to the preparation station (NanoString Technologies, Inc.) with a prepared nCounter Master Kit and cartridge. After binding the sample to the cartridge, 12 lanes per run of the nCounter prep station were run for approximately 2.5-3 h. Cartridges were transferred to the Digital Analyzer (NanoString Technologies Inc.) for analysis and scanned on a digital analyzer at 555 fields of view.
As quality checks for raw data, we confirmed the following: (i) Imaging quality control (QC): Field of View>75%; (ii) Binding Density QC: 0.1<Binding Density<2.25; (iii) Positive Control Linearity QC: R2>0.95; (iv) Limit of Detection QC: 0.5 fM positive control probe>2 standard deviation + mean of the negative controls. Data were normalized using the control probes from Codeset. Normalized expression values of NGAL and 48 cytokine genes were used for the statistical analysis.

Clinical data collection
Baseline demographic data and data about hematological parameters, including hemoglobin levels, white blood cell (WBC) counts, neutrophil and platelet counts, and C-reactive protein (CRP) levels, were collected from medical records. The estimated glomerular filtration rate (eGFR) was calculated using the Chronic Kidney Disease Epidemiology Collaboration equation.

Statistical analyses
Quantitative data are presented as median [quartile1 (Q1), Q3] values. To ascertain that the data were normally distributed, the Shapiro-Wilk test was used. For the statistical analysis of patient demographic features and laboratory parameters, either a parametric or non-parametric test was performed to obtain reliable statistical results. In general, statistical power is higher in parametric tests. However, when collected samples deviated from a normal distribution, non-parametric tests lead to stronger relative power than their parametric counterparts. Consequently, better protection against type II errors relative to parametric tests was ensured. The Kruskal-Wallis H test, a non-parametric counterpart of one-way ANOVA, was conducted on the sample data in which normality was untenable.
For pairwise comparisons when a significant difference was found among groups, Scheffé's post-hoc test or Dunn's pairwise test was applied. In the parametric post-hoc analysis, we applied Scheffé's method because the test is robust for different sample sizes between groups.
A simple regression analysis was performed to analyze the relationship between NGAL expression levels and other cytokine gene levels. Multiple regression analysis was performed to simultaneously analyze the relationship between NGAL gene counts and all statistically significant cytokine gene normalized counts.
To control the inflated Type I error caused by performing multiple comparisons and multiple simple regressions, a Benjamini-Hochberg false discovery rate (FDR) correction was applied [13][14][15][16][17][18]. The correction method adjusts the p values in a simpler but more efficient way than the traditional Bonferroni correction. The adjusted alphas were determined to be 0.018 for pairwise comparisons and 0.026 for multiple simple regressions. The procedure to determine the adjusted alpha has been detailed [13].
In the multiple regression analysis, multicollinearity was also analyzed to identify closely related independent variables. When independent variables presented variance influence factor (VIF) values>10, the independent variables were considered to have multicollinearity. The predictive accuracy of the multiple regression models was assessed using Akaike's information criterion (AIC) and adjusted R 2 values. Statistical significance was set at p<0.05. Statistical analyses were performed using SPSS version v25.0 (IBM SPSS Statistics, Armonk, NY, USA).

Patient characteristics
BM aspirates were collected from 48 patients with a median age of 62 years (range 22-89 years). All 48 patients were at the initial diagnosis stage. Patient demographics, hematological parameters, and normalized counts of NGAL are presented in Table 1. Six variables exhibited significant differences among the four groups. For numerical data, the significant results of the pairwise comparison tests between the groups are shown in Supplementary Table S1.  Table 2 presents the normalized counts of each of the 48 cytokine genes.

Multiple regression analysis of the normalized counts of NGAL and those of other cytokine genes
Using normalized counts of NGAL as the dependent variable and 26 significant predictors (normalized counts of 26 cytokine genes) identified in the simple regression analysis as independent variables, multiple regression analysis was performed. Two multiple regression models were developed ( The adjusted R 2 and AIC values are listed in Table 3. The AIC value was lower for Model 1 than for Model 2, and the adjusted R 2 value for Model 1 was higher than that for Model 2. This suggested that model 1 had higher predictive accuracy for NGAL normalized counts than model 2 (Table 3) [19,20]. However, multiple regression analysis showed multicollinearity between TLR4 (VIF=17.277) and STAT3 (VIF=22.683), two independent variables in model 1 (Table 3). Thus, model 2, in which TLR4 was removed from the predictor variables, was interpreted as more appropriate for predicting NGAL normalized counts than model 1. Because TLR4 and STAT3 normalized counts exhibited multicollinearity, a scatter plot between TLR4 and STAT3 normalized counts was created (Figure 1).
The relationship between the two variables in BM was further analyzed by simple regression analysis, which showed the following relationship:  The median (Q1, Q3) of IL5 normalized counts in the control group was 5.57 (1.52, 11.02). MPN (n=20) had the highest IL5 normalized counts [7.25 (3.80, 11.96)]. The MPN group showed statistically higher IL5 normalized counts than the AML group ( Figure  2C). The AML group exhibited the lowest IL5 normalized counts [1.69 (1.04, 4.57)].

Discussion
We analyzed the association of cytokines' expressions, using BM MNCs at the diagnostic time of myeloid malignancies, which represent clonal diseases of hematopoietic stem cells (HSCs) [21]. BM MNCs include not only hematopoietic progenitor cells at different stages of maturation (drived by HSCs) but also mesenchymal stromal cells (MSCs) comprising the HSC niche [22]. HSC niches are local tissue microenvironments that maintain and regulate stem cells, being often located near trabecular bone and created partly by MSCs and endothelial cells [23]. Recent studies showed that myeloid malignancies are caused by different genetic and epigenetic changes of HSCs and functional changes in HSC niche cells such as MSCs, while the interplay by various kinds of cytokines between HSCs and the niche play an important role [24]. Especially the interplay between both of them were reported to change during disease evolution [24]. In this context, our study displayed the association of cytokines in MNCs (including HSCs and MSCs) involved in the interplay at the diagnostic time of myeloid malignancy, especially presenting the association of NGAL with the relatively well known cytokine genes such as STAT3, TLR4, IL5, TLR1 and TLR2.
Demographic and clinical data showed that AML and MDS groups had significantly lower Hb levels than the control (Table 1; Supplementary Table  S1). Hb, neutrophil counts, platelet counts, and NGAL normalized counts were significantly lower in the AML and MDS groups than in the MPN group (Table  1; Supplementary Table S1). NGAL normalized counts showed a pattern similar to BM NGAL protein levels in previous studies [8,10]. CRP levels were significantly higher in the AML group than those in the control group (Table 1;  Supplementary Table 1). Although CRP levels (protein levels measured in peripheral blood, Table 1) were higher in the AML group, there was no significant difference between the AML and control groups with respect to CRP expression in BM MNCs (Supplementary Table 2). This is because plasma CRP is mainly produced in hepatocytes, but CRP gene expression in our study was measured in BM mononuclear cells [25]. In line with that, in our study, CRP gene showed very low normalized counts among all 49 genes (Table 2), and did not show significant differences between the groups (Supplementary Table 2).
Multiple regression model 2 (Table 3) for predicting NGAL normalized counts included five independent factors, among which, STAT3 normalized counts had the highest explanatory power (R 2 =0.799, Table 2). In sequence, IL5 (R 2 =0.342, Table  2), TLR1 (R 2 =0.303, Table 2), TLR2 (R 2 =0.193, Table 2), and IRAK4 (R 2 =0.170, Table 2) normalized counts exhibited explanatory power. In a previous study using human macrophages, STAT3 was an important transcription factor for NGAL expression [26]. Additionally, NGAL has been suggested to induce the growth and proliferation of breast cancer cells, lung cancer cells, and hepatoblatoma [26]. Similarly, our study based on human BM MNCs showed that NGAL normalized counts had a statistically high association with STAT3 normalized counts. Additionally, NGAL and STAT3 normalized counts showed similar patterns based on disease entity, especially in the MPN and MDS groups (Figure 2 A and B). Further studies including more human BM MNC specimens in several hematological malignancy entities could identify the association between NGAL and STAT3 more clearly, helping to elucidate their developmental mechanisms.
In multiple regression model 2, the TLR4 normalized count was removed as an independent factor because it exhibited multicollinearity with the STAT3 normalized count. However, the TLR4 normalized count was an independent factor with higher explanatory power (R 2 =0.833, Table 2) for the NGAL normalized count than the STAT3 normalized count (R 2 =0.799, Table 2) in simple regression analysis. Additionally, the pattern of NGAL normalized counts in the disease and control groups was more similar to that of TLR4 normalized counts than to that of STAT3 normalized counts (Figure 2 A,  B, and D). Although the relationship between TLR4 and NGAL has been reported in previous studies, it was found in organs or tissues other than BM, such as the urinary system and A549 cells [1,27]. In the urinary system, TLR4 expression is crucial for the secretion of urinary NGAL but does not directly affect NGAL expression [27]. In A549 cells (adenocarcinomic human alveolar basal epithelial cells), the transfection of both TLR4 and its cofactor MD2 led to the upregulation of NGAL expression [1]. Studies using human macrophages and monocytes showed an association between TLR4 and STAT3, but not a direct relationship of TLR4 to NGAL [28,29]. Previous studies based on human monocytes and macrophages suggested that TLR4 does not affect NGAL expression as directly as STAT3 [26,28,29]. Accordingly, TLR4 normalized counts, but not STAT3 normalized counts, were finally removed as an independent factor in our multiple regression model. In multiple regression model 2, IL5 normalized counts had the second highest explanatory power (R 2 =0.342, Table 2), and if the pattern of the STAT3 normalized counts was combined with the pattern of IL5 normalized counts, the combined pattern became more similar to that of the NGAL normalized counts (Figure 2 A, B, and C).
There was no sign of severe infections in AML group of our study. Nevertheless, AML group showed high CRP and low NGAL expression (Table  1), which seems contradictory. As for high CRP levels at the diagnostic time of AML, AML patients with high WHO performance score were reported to have higher CRP levels than those with low score (patients with performance score 3/4 had 6.81 mg/dL higher CRP compared to patients with performance score 0), while in that previous study patients who had a bacteraemic episode within 30 days of AML diagnosis were excluded [30].
However, WHO performance scores of AML patients in our study could not be identified accurately from medical records. Only given that four, two and one AML patients had over 70 years of age, PML-RARA rearrangement and brain hemorrhage, respectively, it could be speculated that AML group in our study would have had relatively high WHO scores. This possibly could lead to high CRP in AML group of our study.
As for low NGAL, a previous study using BM supernatant showed that AML group had low NGAL and high CRP levels [8]. That study suggested that the reason to low NGAL levels in AML group would be due to neutrophils or neutrophilic precursors as synthesis and storage sites of NGAL being suppressed by increased leukemic cells [8].
Nevertheless, as the limitation of our study, since MNCs were isolated from the BM via density, neutrophils would not have been necessarily preserved in the buffy coat and this could have biased the NGAL expression level in our study (more lymphocytes and less myeloid cells). To prevent such bias, future studies measuring the expression of cytokines and analyzing their association at single cell level are needed, which will also help elucidate the association of cytokines in myeloid malignancy more accurately and specifically.
In conclusion, to the best of our knowledge, the present study analyzed for the first time, the association between the normalized counts of NGAL and those of the cytokine genes in human BM MNCs. First, simple regression analysis identified 26 cytokine gene normalized counts that were statistically significantly related to NGAL normalized counts, although STAT3 and TLR4 normalized counts had the highest explanatory power. Stepwise multiple regression analysis was used to develop a multiple regression model as follows: NGAL normalized counts=4316.825 + 9.056 × STAT3 normalized counts + 844.226 × IL5 normalized counts + 17.540 × TLR1 normalized counts -28.206 × TLR2 normalized counts -42.524 × IRAK4 normalized counts. Multiple regression analysis showed that STAT3 and TLR4 normalized counts exhibited multicollinearity and a statistically high association. STAT3, IL5, and TLR4 normalized counts showed similar patterns to NGAL normalized counts in hematological malignancy and control groups. Future studies are warranted to identify and confirm these relationships more accurately and specifically in hematological malignancy diseases, which would help elucidate their developmental mechanisms.