Radiogenomics Profiling for Glioblastoma-related Immune Cells Reveals CD49d Expression Correlation with MRI parameters and Prognosis

Although there have been a plethora of radiogenomics studies related to glioblastoma (GBM), most of them only used genomic information from tumor cells. In this study, we used radiogenomics profiling to identify MRI-associated immune cell markers in GBM, which was also correlated with prognosis. Expression levels of immune cell markers were correlated with quantitative MRI parameters in a total of 60 GBM patients. Fourteen immune cell markers (i.e., CD11b, CD68, CSF1R, CD163, CD33, CD123, CD83, CD63, CD49d and CD117 for myeloid cells, and CD4, CD3e, CD25 and CD8 for lymphoid cells) were selected for RNA-level analysis using quantitative RT-PCR. For MRI analysis, quantitative MRI parameters from FLAIR, contrast-enhanced (CE) T1WI, dynamic susceptibility contrast perfusion MRI and diffusion-weighted images were used. In addition, PFS associated with interesting mRNA data was performed by Kaplan-Meier survival analysis. CD163, which marks tumor associated microglia/macrophages (TAMs), showed the highest expression level in GBM patients. CD68 (TAMs), CSF1R (TAMs), CD33 (myeloid-derived suppressor cell) and CD4 (helper T cell, regulatory T cell) levels were highly positively correlated with nCBV values, while CD3e (helper T cell, cytotoxic T cell) and CD49d showed a significantly negative correlation with apparent diffusion coefficient (ADC) values. Moreover, regardless of any other molecular characteristics, CD49d was revealed as one independent factor for PFS of GBM patients by Cox proportional-hazards regression analysis (P = 0.0002). CD49d expression level CD49d correlated with ADC can be considered as a candidate biomarker to predict progression of GBM patients.


TAM and T cell markers correlated with nCBV, ADC, volume and necrosis values. We performed
Pearson's correlation analysis to determine the association between immune cell markers and MRI parameters, and Fig. 2A, Supplementary Fig. 1 and Table 2 summarize the correlation. The expression levels of CD68, CSF1R, CD33 and CD4 showed significant positive correlations with nCBV values based on both region of interests (ROIs) from FLAIR and CE T1WI and that of CD11b had a significant positive correlation with nCBV values only from CE T1WI ( Fig. 2B; Case 1 and 2). We found significant negative correlations between the expression levels of CD49d and CD3e and ADC values from both FLAIR and CE T1WI, and those of CD33 and CD123, and CD25 were negatively correlated with ADC values from FLAIR and CE T1WI, respectively ( Fig. 2B; Case 3 and 4). Tumor volumes based on FLAIR or CE T1WI had significant negative correlations with the expression levels of CD123, CD49d and CD117, but no immune cell markers showed a significant correlation with tumor necrosis or necrosis ratio.

Discussion
In this study, we applied radiogenomics profiling to the association between quantitative MRI features and expression levels of immune cell markers in GBM patients, which was also correlated with PFS. We found that myeloid lineage immune cell markers (e.g., CD163, CD63, CD83 and CD49d) were more dominant than lymphoid lineage markers in GBMs. In terms of CBV values based on FLAIR or CE T1WI, the expression levels of CD68, CSF1R, CD33, CD4 and CD11b showed positive correlations. In addition, the expression levels of CD49d, CD3e, CD33, CD123 and CD25 were found to have negative correlations with ADC values based on FLAIR or CE T1WI. The expression levels of CD123, CD49d and CD117 had negative correlations with tumor volumes. PFS was also correlated with myeloid cell markers (e.g., CD11b, CD123, CD33, CD163, CD63 and CD49d) and lymphoid cell marker genes (e.g., CD25 and CD8). Among them, CD49d was the most important prognostic marker regardless of genetic status of GBM.
TAMs are known to have a role in tumor angiogenesis; therefore, we assumed that TAMs can affect the increase in nCBV values, which can be simultaneously correlated with TAM markers (e.g., CD11b, CD68 and CSF1R). Recent studies have elaborated an important role of CSF1R-colony stimulating factor-1 receptor in the brain TAM; for example, inhibition of CSF1R either reduces 28 or depolarizes TAMs 29 . In addition, CD33, a MDSC marker, has high correlation with nCBV values. The suppressive activity of MDSC is one of the most prevalent mechanisms of immune evasion in patients 30 , which occurs as a consequence of the aberrant myelopoiesis that arises in cancer. MDSCs are mobilized during tumorigenesis and infiltrate developing tumors, where they promote tumor angiogenesis 31 and disrupt major immunosurveillance mechanisms. Among the lymphoid lineage markers, CD4 is significantly positively correlated with nCBV values. As we described above, nCBV seem to be related to the function of TAMs, which induce angiogenesis in tumors. Regulatory T cell and type 2 helper T cells (T H 2) marked by CD4 express Interleukin-10 32 , Interleukin-4 and Interleukin-13 33 , which are cytokines for M2-type macrophage polarization 11,34,35 . In other words, TAMs can be induced primarily by pro-tumorigenic roles through suppressive immunosurveillance and anti-inflammatory cytokines, which can cause increase in nCBV resulting from angiogenesis, one of the key roles of TAMs in GBM 10,23 .
CD49d, CD33, CD123, CD3e, and CD25, which mark bone marrow-derived cells (MDSCs) dendritic cells (DC), T helper cells, cytotoxic T cells and regulatory T cells, respectively, are negatively correlated with ADC values. Dendritic cells are known as antigen presenting cells derived from the bone marrow, which induces innate and adaptive immune responses 36 . When tumors occur, immune cells, including MDSCs, which are the precursor of DCs, macrophages and granulocytes 9,12,13 , DCs and T cells, are recruited from bone marrow and converge on the tumor side by the immune system 13 . We believe that tumor cellularity, measured by ADC value, can be affected by immune cell infiltration as well as tumor tissue itself. In addition, myeloid cells seem to play a more important role than lymphoid cells in tumor cellularity. In addition, CD123, CD117, and CD49d, which mark DCs, mast cells and bone marrow-derived cells 8 , respectively, are negatively correlated with tumor volume. Before they are recruited by tumors, DCs and mast cells induce innate and adaptive immune responses to regress tumors and prevent relapse 9 . CD49d was revealed as a marker of hematopoietic bone marrow-derived cells 37 . DCs and mast cells activate T cells to remove tumors 8,12 , which can be explained by the finding that these immune cells recruit in the early stage of gliomagenesis.
In terms of PFS, as the previous study has been reported that high CD49d expression is associated with poor survival in chronic lymphocytic leukemia 38 , it is remarkable that CD49d was discovered an independent biomarker regardless of any other molecular characteristics, even though TAM markers are prominent among the immune cell population and correlate with several imaging features. CD49d is also known as integrin α4, forming a heterodimer with integrin β1 (VLA4), which is the specific marker of hematopoietic cells, namely, bone marrow-derived cells. Bowman RL et al. demonstrated that microglia specifically repress integrin α4 (CD49d), enabling its utility as a discriminatory marker between microglia and bone marrow-derived macrophages in primary and metastatic disease in both mouse and human 37 . Several studies reported that the brain-resident microglia and the infiltrating monocytes/macrophages of blood are the major glioma-associated inflammatory cells that constitute the tumor microenvironment 39,40 . Particularly, a recent report 41 and a clinical study 42 revealed that monocytes/macrophages, but not microglia and lymphocytes, are the most predominant TAMs in GBM. Monocytes, cells of the myeloid lineage, are released during inflammation and differentiate into macrophages to maintain immune homeostasis 43 . A previous study suggested that circulating monocytes are cytotoxic to tumor cells 44 ; however, when monocytes reach the tumor mass, the tumor molecular milieu induces differentiation to new cell types per tumor requirement 45 . A recent report indicated that reduction in the tumor-promoting effects of monocytes/macrophages in GBM can be considered as an adjuvant treatment for glioma 46 . However, the fate of the GBM-adhered monocytes/macrophages and their effect on GBM growth are still obscure. Adhesion molecules are known to mediate cell-cell interactions, particularly between immune cells and target tissues. Vascular cell adhesion molecule-1 (VCAM-1) is an inducible adhesion molecule that facilitates tight attachment to the monocyte/macrophage-associated integrin α4β1 (VLA4) 47 . VCAM-1 expressed on the surface of tumors interacts with VLA4 on monocytes/macrophages, which promotes tumor invasion, angiogenesis and metastasis 10 . In short, higher CD49d expression level in GBM is thought to induce TAMs to adhere to GBM, increasing cellularity and resulting in a decrease in ADC on MRI, attenuated tumor aggressiveness and poorer progression of patients. The previous studies have been reported that MRI parameters (e.g., nCBV and ADC) are used for evaluation of GBM patients' overall survival (OS) and PFS and they are related to poor progression of GBM. For instance, high nCBV and low ADC values are results of poor PFS 48 . Therefore, our findings related to immune cell markers and MR imaging parameters can follow antecedent researches.
MRI protocol. All patients underwent conventional, DWI and DSC perfusion MRI using a 3 T scanner (Verio; Siemens Healthcare Sector) with a 32-channel head coil. The conventional MRI included T1-weighted imaging (T1WI), such as transverse spin-echo imaging, before and after contrast enhancement or multi-planar reconstructed transverse, coronal imaging with a sagittal three-dimensional magnetization prepared rapid acquisition gradient echo (3D-MPRAGE) sequence before and after contrast enhancement, and transverse T2-weighted imaging (T2WI) with turbo spin-echo sequences and FLAIR images. Contrast-enhanced (CE) DWI was performed with a single-shot, spin-echo, echo-planar imaging (EPI) sequence in the axial plane before the injection of contrast material with b-values of 0 and 1000 sec/mm 2 , a TR of 6300 ms, a TE of 92 ms, an FA of 180°, a matrix of 240 × 240, an FOV of 240 × 240, a section thickness of 3 mm and an NEX of 3. DWI was acquired in three orthogonal directions and combined into a trace image. ADC maps were calculated on a voxel-by-voxel basis with the software that was incorporated into the MRI unit using these data.
The transverse DSC perfusion MRI was obtained with single-shot, gradient-echo, echo-planar sequences during the intravenous administration of gadobutrol at a concentration of 0.1 mmol/kg of body weight at a rate of 4 ml/sec using a power injector (Spectris; Medrad). A 30-ml bolus injection of saline followed at the same injection rate. For each section, 60 images were acquired at intervals equal to the TR. The parameters were as follows: TR, 1500 ms; TE, 30 ms; FA, 90°; matrix, 128 × 128; section thickness, 5 mm; intersection gap, 1 mm; FOV, 240 × 240 mm; sections, 15-20; voxel size, 1.875 × 1.875 × 5 mm 3 ; pixel bandwidth, 1563 Hz; and total acquisition time, 1 minute 30 seconds.
Image post-processing and data analysis. The conventional MR images, ADC maps, and DSC PWI were digitally transferred from the picture archiving and communication system workstation to a personal computer for further analysis. The relative CBV (rCBV) was obtained with a dedicated software package (nordicICE; Nordic Imaging Lab, Bergen, Norway) that applied an established tracer kinetic model to the first-pass data 49,50 . First, realignment was performed to minimize patient motion during the dynamic scans. A gamma-variate function, which approximates the first-pass response as it would appear in the absence of recirculation, was used to fit the 1/T2* curves to reduce the effects of recirculation. To reduce the contrast agent leakage effects, the dynamic curves were mathematically corrected by using leakage correction package available on the dedicated software 51 . After the elimination of recirculation and leakage of the contrast agent, rCBV was computed with numeric integration of the curve. To minimize variances in rCBV in an individual patient, the pixel-based rCBV maps were normalized by dividing every rCBV value in a specific section by the rCBV value in the unaffected white matter, and finally normalized rCBV (nCBV) maps were generated 52 .
Using a dedicated software package (nordicICE), co-registrations between the structural images (e.g., FLAIR images and CE T1WI) and the nCBV and ADC maps were performed based on geometric information stored in the respective data sets. The differences in the slice thickness between images were corrected automatically by re-slicing and co-registration based on the underlying structural images. The nCBV and ADC maps were displayed as color overlays on the both FLAIR images and CE T1WI.
One neuroradiologist (S.H.C. with 16 years of brain MR imaging experience) who was blinded to the clinical data drew polygonal ROIs that contained the entire enhancing lesions in each section of the co-registered images. Areas of necrosis, hemorrhage, or non-tumor macro-vessels that were evident on the CE T1WI were excluded from the ROIs. Then, the ROIs of T2 high SI lesions, regardless of contrast enhancement, were also defined on each transverse FLAIR image, avoiding the cystic and necrotic regions and the macrovessels. Because the ROI placement was conducted on the nCBV and ADC map co-registered with structural images, the margin of the  lesions could be defined with confidence. The entire volume of contrast-enhancing lesions, T2 high SI lesions, and necrosis, which was defined as a hypointense area without contrast enhancement on CE T1WI within the mass on the FLAIR images, was calculated. The data acquired from each section were summed to derive the voxel-by-voxel ADCs and nCBVs for the entire tumor extent based on both CE T1WI and FLAIR images by using nordic ICE. The ADC and nCBV histograms were plotted with ADC and nCBV on the respective x-axis with a bin size of 3 × 10 −5 mm 2 /sec and 0.1, respectively, whereas the y-axis was expressed as a percentage of the total lesion volume by dividing the frequency in each bin by the total number of analyzed voxels. For further quantitative analysis, the cumulative number of observations in all bins up to the specified bin was mapped on the y-axis as a percentage in the cumulative histograms. The 5th percentile point for ADC (5% ADC) and 95th percentile point for nCBV (95% nCBV) were derived (the Xth percentile point is the point at which X% of the voxel values that form the histogram are found to the left of the histogram) 53,54 . Immunohistochemistry staining. Antibodies for IHC analysis included mouse anti-human Integrin alpha4 (CD49d) antibody (sc-365209, Santa Cruz Biotechnology, Inc.) and rabbit anti-human Glial fibrillary acidic protein (GFAP) antibody (ab7260, abcam). Paraffin sections (4 μm) were dewaxed and rehydrated, Antigen retrieval was performed in a microwave by placing the sections in epitope retrieval solution (0.01 M citrate buffer, pH 6.0) for 20 minutes, then incubated in 3% hydrogen peroxide for 20 min at room temperature. The paraffin sections were then blocked with 3% BSA for 30 min, stained with antibodies for 1 hour at room temperature, washed with wash buffer (S3006, Dako) and stained with secondary antibody (A-11008, Alexa Fluor 488 by Invitrogen TM ; A-11032, Alexa Fluor 594 by Invitrogen TM ) for 30 min at room temperature. DAPI was used to stain the nucleus.
Statistical Analysis. All statistical analyses were performed using two commercial software programs (MedCalc version 17.2, MedCalc Software). A P value < 0.05 was considered statistically significant. Kolmogorov-Smirnov's test was used to determine whether the variables followed normal distribution. Non-parametric data are presented as median and interquartile range (IQR, range from the 25th to the 75th percentile), and parametric data are shown as the mean ± standard deviation. Pearson's correlation analysis for parametric data was performed for the correlation between the expression level of immune cell markers and quantitative imaging parameters.
The progression-free survival (PFS) was assessed using the Kaplan-Meier method according to the expression level of immune cell markers, which were compared using log-rank tests. GBM progression was defined according to RANO criteria 55 . We only recorded the first progression. PFS was calculated from the date of surgery to that of GBM progression, death, the last confirmation of no evidence of disease, or the most recent follow-up examination. Patients without an event were censored at the date of the most recent follow-up, regardless of whether they were scheduled for future follow-up or they had been lost to follow-up. Eight patients who expired from progression-unrelated conditions (e.g., infarction and infection) were excluded from PFS analysis. To determine thresholds in each immune cell marker expression level for PFS, receiver operating curve analysis was used. Multivariate analysis was performed using the Cox proportional hazards model, which was adjusted for the prognostic factors including the expression level of immune cell markers and IDH1 mutation status.