Study of Diffusion Weighted Imaging Derived Diffusion Parameters as Biomarkers for the Microenvironment in Gliomas

Objectives To explore the efficacy of diffusion weighted imaging (DWI)-derived metrics under different models as surrogate indicators for molecular biomarkers and tumor microenvironment in gliomas. Methods A retrospective study was performed for 41 patients with gliomas. The standard apparent diffusion coefficient (ADCst) and ADC under ultra-high b values (ADCuh) (b values: 2500 to 5000 s/mm2) were calculated based on monoexponential model. The fraction of fast diffusion (f), pseudo ADC (ADCfast) and true ADC (ADCslow) were calculated by bi-exponential model (b values: 0 to 2000 s/mm2). The apparent diffusional kurtosis (Kapp) was derived from the simplified diffusion kurtosis imaging (DKI) model (b values: 200 to 3000 s/mm2). Potential correlations between DWI parameters and immunohistological indices (i.e. Aquaporin (AQP)1, AQP4, AQP9 and Ki-67) were investigated and DWI parameters were compared between high- and low-grade gliomas, and between tumor center and peritumor. Receiver operator characteristic (ROC) curve and area under the curve (AUC) were calculated to determine the performance of independent or combined DWI parameters in grading gliomas. Results The ADCslow and ADCuh at tumor center showed a stronger correlation with Ki-67 than other DWI metrics. The ADCst, ADCslow and ADCuh at tumor center presented correlations with AQP1 and AQP4 while AQP9 did not correlate with any DWI metric. Kapp showed a correlation with Ki-67 while no significant correlation with AQPs. ADCst (p < 0.001) and ADCslow (p = 0.001) were significantly lower while the ADCuh (p = 0.006) and Kapp (p = 0.005) were significantly higher in the high-grade than in the low-grade gliomas. ADCst, f, ADCfast, ADCslow, ADCuh, Kapp at the tumor center had significant differences with those in peritumor when the gliomas grade became high (p < 0.05). Involving ADCuh and Kapp simultaneously into an independent ADCst model (AUC = 0.833) could further improve the grading performance (ADCst+ADCuh+Kapp: AUC = 0.923). Conclusion Different DWI metrics fitted within different b-value ranges (low to ultra-high b values) have different efficacies as a surrogate indicator for molecular expression or microstructural complexity in gliomas. Further studies are needed to better explain the biological meanings of these DWI parameters in gliomas.


INTRODUCTION
Gliomas are the most common primary brain tumors in adults and according to World Health Organization (WHO) guidelines are classified into four grades, which reflect increasing malignancy and worse prognosis (1). Assays of histopathology and molecular pathology based on tumor samples obtained by resection or biopsy are the gold standard for determining the pathological grade and molecular subtype. Accurate assessment of glioma grade, as well as phenotype and genotype, is of potential importance for the optimization of personalized treatment. However, inherently high heterogeneity of gliomas means that biopsy or localized resection may not be representative of the tumor as a whole.
Improving the non-invasive characterization of glioma physiology or pathology might help to improve the imageguided biopsy and therapy. Considering different regimes of bvalue could control the degree of diffusion-weighting in the diffusion weighted imaging (DWI), different tissue properties reflected by the water diffusion could be encoded into DWI signals. Therefore, DWI is one of the potential tools to provide surrogate noninvasive imaging biomarkers for microenvironment in gliomas as the water diffusion coefficients could intermediately reflect the microstructure, perfusion or water exchange effects associated with transmembrane transport, such as facilitated diffusion (2,3).
As the cell proliferation and the water transportation are highly suspected to influence the water diffusion properties in gliomas at extracellular, intracellular or transcellular space, Ki-67 or aquaporin (AQP) subtypes (AQP1, AQP4, AQP9) were welcomed molecular targets quantified by different DWI diffusion metrics (4,5). Ki-67 is an immunohistochemical marker for the proliferation in gliomas which is known to correlate with tumor grading (4) and prognosis (6). AQPs provide a major pathway for the water transportation through cell membrane (7). And AQP subtypes, including AQP1, AQP4 and AQP9, are overexpressed in glioma cells and correlated with tumor grade and malignancy (8)(9)(10)(11)(12). AQP1, AQP4 and AQP9 were reported to be related to angiogenesis, invasion and peritumoral edema in gliomas. AQP1 predominantly locates in the perivascular space and it has been reported that increased AQP1 might induce vasogenic brain edema (13) and acceleration of cell migration and invasion (10,14). AQP4 are mainly expressed by astrocytes and its redistribution is thought to control water mobility at the blood-brain interface and progress along with blood-brain barrier disturbance and vascular proliferation (11,15,16). Increased expression of AQP9 has been observed in glioma tissue near vessels and might promote the invasion and motility of cells (17,18).
Several studies have demonstrated that some diffusion coefficients derived from intravoxel incoherent motion (IVIM), diffusion kurtosis imaging (DKI) or stretched-exponential models such as slow apparent diffusion coefficient (ADC slow ), axial kurtosis and heterogeneity index a could correlate with Ki-67 expression in gliomas (4,19,20), whilst the DKI-derived mean kurtosis, diffusion tensor imaging-derived mean diffusivity, the apparent ADC and ADC derived from ultrahigh b-value model showed different correlation tendency with AQP in gliomas (5,21). However, the levels of the correlations seemed to be diverse and few studies have compared the efficacy of different kinds of model-derived parameters for indicating Ki-67 and AQP in gliomas among the same dataset.
Considering more DWI models to select more rational or efficient diffusion surrogate biomarkers for indicating molecular expression or phenomenologically depicting the gliomas microenvironment is necessary at present. In the current study, the diffusion metrics derived from mono-exponential model, IVIM model, DKI model and ultra-high b-value model, and their efficacy for quantifying Ki-67, AQP1, AQP4, and AQP9 were analyzed and compared. Among these models, the monoexponential model could reflect apparent diffusion combining diffusion and perfusion; IVIM model could separate the slow diffusion in response to intra-or extracellular water molecules from the fast diffusion in response to tubular or vascular perfusion; DKI model could help provide general structure information and quantitative information about diffusion deviation from freely gaussian diffusion; and ultra-high b-value model could provide the diffusion information in the specific high b-value range which is potential to indicate transmembrane diffusion (2,3). By comprehensively considering these models covering gaussian or non-gaussian assumption and b-value compartment, we expected to contribute as many shreds of evidence as possible during the formation of DWI biomarkers for the gliomas microenvironment. In addition, as the AQPs distribution indicates the disease progression (9,22), the spatial distribution of the diffusion coefficients might noninvasively indicate physiological or pathological conditions related to the change of water exchange or transportation. Therefore, the same set of diffusion metrics from central and peritumoral regions of low-and high-grade gliomas were also analyzed and compared to help depict a more complete tumor microenvironment picture.

Patients Population
This retrospective study was approved by the local institutional review board, and written informed consent was obtained from each patient. A total of 45 patients with pathologically proven gliomas diagnosed between October 2014 and May 2016 were enrolled in the study. The inclusion criteria were as follows: (a) MRI examinations were performed on patients prior to treatments of tumors and (b) the pathological diagnoses and histological indices were acquired by surgical resections of gliomas. Four patients were excluded due to the presence of head movement artifacts in the DWI images. The final analysis was performed for a total of 41 patients with gliomas.

MRI Data Analysis
Two radiologists, blinded to the reports concerning tumor pathology, reviewed and analyzed all the MR images independently on a remote workstation. The radiologists independently drew three different regions of interests (ROIs) on the T2w echo-planar image for each tumor in the solid parts and regions within 1 cm peritumoral parts, respectively. Each MRI parameter in the solid parts or peritumoral region were determined by the averaged value of three ROIs, respectively. Areas of necrosis, hemorrhage and cerebrospinal fluid were excluded to ensure accurate measurements.

DWI Data Processing
DWI data were transferred to a workstation (Advantage Workstation 4.5; GE Healthcare) for processing. ADC st was calculated from b values of 0 and 1000 sec/mm 2 by using a monoexponential DWI model (23) as Equation (1): where S(b) represents the signal intensity in the presence of diffusion sensitization, and S(0) represents the signal intensity in the absence of diffusion sensitization. This model used the least square fit for linear fitting (24).
ADC slow was obtained from IVIM model with b values from 0 to 2000 sec/mm 2 (23) as Equation (2): Eq: (2) where f represents the fraction of fast diffusion component, ADC fast represents the pseudo-diffusion coefficient, and ADC slow represents the slow diffusion coefficient. The Levenberg-Marquardt fit was used for nonlinear fitting (24). ADC uh under ultra-high b values was calculated from b values of 2500, 3000, 4000 and 5000 sec/mm 2 by using the above monoexponential DWI model (25) as Equation (1).
Apparent diffusional kurtosis (K app ) was calculated from DKI model (26) according to Equation (3) with b-value ranging from 200 to 3000 sec/mm 2 : where D app (unit: ×10 -3 mm 2 /s) is the apparent diffusion coefficient fitted in low b-value range of 200-1000 sec/mm 2 , K app (unitless) is the apparent diffusional kurtosis which is fitted with b-value up to 3000 sec/mm 2 .

Immunohistochemistry Data Analysis
Specimens acquired from surgical resections were embedded in paraffin and immunohistochemistry stains were applied for quantification of AQP1, AQP4, AQP9 and Ki-67. Histological indices of AQP1, AQP4, AQP9 and Ki-67 were independently measured by two pathologists by using a HMIAS-2000 Medical Color Image Analysis System (Champion Image Engineering Co., Ltd., Wuhan, China). The pathologists independently placed three different ROIs in the solid parts of each tumor for each patient. The ROIs excluded the areas of necrosis and hemorrhage. The IODs were measured for AQP1, AQP4, AQP9 and Ki-67, which were averaged from three delineated ROIs.

Statistical Analysis
MedCalc software (version 19.0, MedCalc, Belgium) and R software (version 3.5.1) were used for statistical analyses.
Correlations between DWI-derived parameters and histological indices were computed by using the Pearson's correlation analysis. Steiger's Z-test was used for comparing each two correlation coefficients (27). A correlation coefficient (r) of 0.75-1.00 was set to indicate very good to excellent correlation; 0.50-0.74, moderate to good correlation; 0.25-0.49, fair correlation; and 0.24 or lower indicate little or no correlation (28,29). The Mann-Whitney U test was used for statistically comparing parameters between high-grade and low-grade gliomas. The significant difference of DWI-derived parameters between tumor center and peritumoral regions was analyzed by using the intragroup paired t-test or Wilcoxon matched-pairs signed rank test. A receiver operator characteristic (ROC) curve and area under the ROC curve (AUC) were calculated to evaluate the model performance in grading gliomas. The maximum Youden index was used to determine the threshold for calculating the sensitivity and specificity. AUC of the paired models were compared by Delong's test. The continuous net reclassification improvement (NRI), and integrated discrimination improvement (IDI) indices were analyzed to assess the added value of combined models (30). Interobserver agreement for each measurement was calculated by using an Intraclass Correlation Coefficient (ICC) with 95% confidence interval (CI). Two sides p values less than 0.05 were considered statistically significant. The R packages were mainly involved as follows: "icc" was used for ICC calculation by setting "twoway" and type of "agreement", "glmnet" package for logistic regression, "pROC" package for ROC analysis, "PredictABEL" was used for NRI and IDI evaluation, and "ggboxplot" for boxplot.

Differences in DWI-Derived Parameters and Histological Indices Between High-Grade and Low-Grade Gliomas
The significant differences in DWI-derived imaging parameters and histological indices between the high-grade and low-grade gliomas are listed in Table 1. Values of ADC st and ADC slow at the center of the tumor were significantly lower in high-grade gliomas than in low-grade gliomas, whereas the value of ADC uh at the center of the tumor was significantly increased in the high-grade gliomas than in the low-grade gliomas. Significant higher K app at tumor center were found in highgrade gliomas than in low-grade gliomas. Analysis of the histological indices showed that the expression of AQP1, The variables with abnormal distribution were depicted by mean ± SD. *P-value < 0.05 indicated statistical significance.
AQP4 and Ki67 were significantly higher in the high-grade than in the low-grade gliomas, whereas the expression of AQP9 showed no significant difference between the high-grade and low-grade gliomas. The MRI and immunohistochemistry staining of the high-grade and low-grade gliomas are illustrated in Figures 1, 2, respectively.

Correlation Between DWI-Derived Parameters and Histological Biomarkers
The correlation analysis between DWI-derived imaging parameters and histological indices revealed that expression of AQP1, AQP4 and Ki67 were correlated with ADC st , f, ADC slow and ADC uh at the center of the tumor while had little correlation with those of peritumoral imaging parameters (see Table 2 and Figure 3). As shown in Table 2, the ADC st at the center of the tumor showed a small negative correlation with AQP1, AQP4 and Ki67 (p < 0.05). The f coefficient at the center of the tumor showed a small negative correlation with AQP1 and Ki67 (p < 0.05). The ADC slow at the center of the tumor showed a small negative correlation with AQP1, AQP4 and a moderate negative correlation with Ki-67 (p < 0.05). The ADC uh at the center of the tumor showed a small positive correlation with AQP1, AQP4 and Ki67. The DWI metric of Kapp at the center of the tumor showed weak positive and moderate negative correlation with Ki-67 expression (p < 0.05), while its correlation with AQPs was not obvious. There was no significant correlation between the expression of AQP9 and any of the imaging parameters (p > 0.05).
In addition, the correlation analysis results between each pair of DWI-metrics derived from tumor center were summarized in Table S2. The ADC st showed significant correlation (p < 0.05) with other DWI metrics (f, ADC slow , ADC uh and K app ), while ADC uh had no significant correlation (p > 0.05) with ADC slow and K app . The Steiger's Z test results by comparing different DWI metrics' correlations with the same histological biomarker were also summarized in Table S3. It indicated that there was no significant difference between ADC st and ADC slow for its respective correlation with an individual molecule. While significant differences existed between ADC uh and ADC st or ADC slow because of the inverse correlation tendency.

Differences in Imaging and Histological Biomarkers Between Tumor Center and Peritumoral Areas
The DWI-derived imaging parameters at the center of the tumor and peritumoral area in low-and high-grade gliomas were statistically analyzed as shown in Figure 4.
In summary, all the DWI-imaging parameters (including f, ADC st , ADC fast , ADC slow , ADC uh and K app ) at the tumor center had significant differences compared with those in peritumoral areas when the glioma grade became high. And the ADC uh and K app at tumor center were significantly higher than that of peritumor in both low-grade gliomas and high-grade gliomas.

Diagnostic Performance of Imaging and Histological Biomarkers for Differentiation Between Low-and High-Grade Gliomas
The ROC analysis for significant DWI imaging metrics measured from tumor center and histological biomarkers in discriminating high-grade gliomas from low-grade gliomas were conducted and compared as shown in Figure 5 and Table 3.  Among molecular biomarkers, Ki-67 presented a good diagnostic ability with the largest AUC and good sensitivity and specificity, which was followed by AQP1. AQP4 presented the best sensitivity whereas poor specificity was determined at a threshold of 0.14. For DWI imaging biomarkers, the ADC st had the largest AUC and relatively better sensitivity and specificity which was followed by ADC slow . The ADC uh and K app performed with a similar AUC of around 0.76. While ADC uh had a good sensitivity the same as the ADC st at a threshold of 0.24 and K app had a better specificity (80%) at a threshold of 0.512 compared with other independent imaging metrics. Meanwhile, we assessed the discrimination performance by combining two to three DWI metrics among ADC st , ADC slow , ADC uh and K app . From the aspect of AUC performance, the (ADC st + ADC uh : AUC = 0.895) and (ADC uh + K app : AUC = 0.882) represented better improvement among two-parameter-based combination models, and (ADC st + ADC uh +K app : AUC = 0.923) showed better improvement among three-parameter-based combination models. The paired AUC comparison results derived from the Delong's test were summarized in Table S4 and the continuous NRI and IDI were summarized in Table S5. Delong's test result showed that ADC uh at tumor center versus Ki-67 and AQP4

DISCUSSION
In the current study, DWI parameters metrics derived from low to ultra-high b values based on mono-exponential model, IVIM model, DKI model and ultra-high b-value model were analyzed from three aspects: (1) their association with expression of histological molecular biomarkers including Ki-67, AQP1, AQP4, and AQP9; (2) if they could be taken as noninvasive surrogated indicators for depicting microenvironment as tumor progression; (3) their preoperative diagnosing ability to differentiate low-and high-grade gliomas.  Our histological results were consistent with previous studies that the AQP1, AQP4 and Ki-67 expression were significantly increased in the high-grade gliomas compared with the lowgrade gliomas (31,32). Our dataset also showed that ADC st , ADC slow decreased and ADC uh , K app increased significantly in high-grade gliomas compared with low-grade ones.
The alteration of DWI metrics might be explained by molecular origins. The data presented in our study demonstrated that different DWI models fitted within different b-value regimes could influence the derived DWI metrics' correlation level or tendency with histological biomarkers expression. Most of the DWI metrics derived from tumor center had correlations with Ki-67 which indicated tumor cell proliferation-related diffusion restriction effect (33). The ADC slow and ADC uh showed stronger correlations with Ki-67 than other DWI metrics in the current study. It was consistent with previous studies that ADC slow eliminated the effect of microcirculation and could better reflect the water restriction from tumor density and extracellular volume or deposition changes during tumor proliferation (4,34). It has been also reported that ADC value obtained from high-b value (b = 3000 sec/mm 2 ) had a stronger correlation with Ki-67 index, which was in agreement with the current result for ADC uh (35). ADC uh might give more chance to indirectly reflect the slower diffusion component involved in Ki-67-related proliferation.
The ADC st , ADC slow and ADC uh presented correlations with AQP1 and AQP4 in the current study, while AQP9 had no correlation with any DWI metrics. The ADC st and ADC uh presented stronger correlations with AQP1 expression. For AQP4, the ADC st and ADC slow showed stronger correlations.
With the knowledge of slower water transportation speed (36) resulted from the AQPs, previous studies indicated that ADC derived in the high-b value range might be one potential marker for AQPs expressions in gliomas (5) or other diseases (25,37). However, our results showed that the ADC uh was not completely independent from other DWI metrics such as ADC st to independently indicate AQPs expression. And there were already a few studies that reported a negative correlation between ADC st and AQP4 during brain injury or ischemia, which were hypothesized to be associated with decreases in the extracellular space caused by cell swelling (28,38). The water diffusion condition was originated from the combined contribution of the water exchange effect by AQP overexpression and the proliferation effect by Ki-67. Upregulations of AQP1 and AQP4 could enhance migration and invasion of glioma cells. In the migration or invasion process, the AQPs might be involved in or play key roles in coordinated cell volume changes (39), reduced tumor cell adhesion with surrounding cells (40), degradation of extracellular matrix (41) and rapid transmembrane water fluxes during lamellipodia formations of glioma cells (14). In addition, AQP1 upregulation is associated with angiogenesis (32) and tumorassociated edema formation (42) and AQP4 redistribution might be functional in the reabsorption of excess cerebral fluid during vasogenic edema (38,43,44). These processes might influence the water diffusion in the tissue and transmembrane or microcirculation. In addition, mean kurtosis, as one potential DWI metric measuring the degree of diffusion hindrance or restriction, has been reported to correlate with AQP4 expression in gliomas (21). In the current dataset, the K app showed a larger value in high-grade gliomas which was in accordance with previous studies (45,46). However, it only showed a significant weak positive correlation with Ki-67 which indirectly reflects proliferation-related diffusion barriers (47). No significant correlation with AQP1 or AQP4 was found for K app in the current research. As parameters derived from models involving high-b value, ADC uh and K app had different correlation levels with AQPs which might be resulted from the different b-value range for fitting. But the current dataset could not provide sufficient explanations for such deviation. More pathological evidences and a larger sample size should be supplemented to reveal the actual mechanism for AQP-related diffusion effect.
The statistical differences between the DWI metrics derived from tumor center and peritumoral area were analyzed to explore if DWI metrics could be useful to noninvasively describe the glioma microenvironment. In the current study, as gliomas grade increased, the DWI-derived parameters in the tumor center deviated more significantly from those in the peritumoral area. The ADC st at the center of the tumor decreased significantly than that of the peritumoral area in high-grade gliomas which might be a good indicator for the enhanced cellularity in the solid tumor center (48). However, the ADC st might ignore the perfusion effect resulted from tumor vasculature. Therefore, bi-exponential IVIM model was used to distinguish the perfusion and diffusion effects. As expected, the ADC fast increased, ADC slow decreased and K app increased respectively at tumor center compared with those at the peritumoral area in the high-grade gliomas whereas these differences were not significant in the low-grade gliomas. Firstly, the decrease of "true" diffusion coefficient ADC slow and K app in the tumor center might represent the more complicated or heterogenous tumor microstructure, which was in accordance with proliferation-related Ki-67 expression (49). Secondly, as ADC fast (b < 200 sec/mm 2 ) was able to reflect effective perfusion index (50), the increase of ADC fast at tumor center might be well correlated to perfusion-related angiogenesis of gliomas, which was consistent with the previous study (24). It has been reported that the growth of malignant glioma (WHO grade III and IV) is dependent on new neovascularization which may result in increased permeability, blood flow and transport properties (51,52). The relative cerebral blood volume increased more frequently in the high-grade gliomas than in the low-grade gliomas (53). Moreover, the AQP1 upregulation might also contribute to the angiogenesis in high grade gliomas (41). The volume fraction of fast pool f is mainly affected by the fraction of capillaries and microcirculations. But the f parameter at tumor center and peritumor showed an inverse tendency compared with perfusion-related ADC fast , in which the f became smaller at tumor center as glioma grade increase. Such controversy between f and ADC fast existed in a previous report (24). The possible reason for this is that the IVIM model doesn't involve the influence of echo time to f. The T2 value of brain tissues including white or grey matter under 3.0T MRI is greatly deviated from blood's T2 value. Such deviation might be the reason for the uncertainty of f value. Although the current study showed correlations for AQP expression with ADC uh , ADC st or ADC slow , we still cannot infer that the AQP1 or AQP4 distribution in the tumor center might be denser than peritumor as the immunohistochemical results were only obtained for the solid tumor part. In addition, it has been reported that AQP1 up-regulation is associated with angiogenesis in gliomas and is predominantly located perivascularly and in areas of tumor infiltration whereas distant from the tumor center (32). AQP4 was also reported to have higher expression in both tumor and peritumor than in normal tissues in gliomas, but the degree of peritumoral edema only positively correlates with the expression level of AQP4 in peritumor (9,22). Therefore, further research is warranted to clarify the spatial distribution of AQPs expressions or tissue arrangement by using pathological, molecular imaging or functional imaging methods and contribute more explanations for DWI metrics.
Finally, the diagnostic performance of the DWI parameters (ADC st , ADC slow, ADC uh and K app ) in differentiating high-grade and low-grade gliomas. We demonstrated that the ADC st had largest AUC 0.833 which was followed by ADC slow. While ADC uh had a good sensitivity the same as the ADC st . Although there have been reported that conventional ADC st performed well in grading gliomas (47,54), the added value of other DWI metrics derived from IVIM, ultra-high-b value or DKI model is still worthy of study (5,19,55). The current research revealed that increased NRI and IDI indices could be obtained by introducing ADC uh and K app into independent ADC st in grading gliomas.
Comprehensive analysis of the relationship between DWIderived parameters and indices of histopathology has demonstrated the potential of establishing imaging biomarkers reflecting the tumor microenvironment in gliomas. However, further research is required to address several limitations of the present study. Firstly and most importantly, whereas imaging ROIs were placed both at the center of the tumor and in the peritumor, histology assays were not performed separately for these two areas and tumor heterogeneity may have resulted in selection bias. In future studies techniques should be applied so that the imaging and histology refer to the same sampling locations. In addition, by considering the limited sampling area for pathological methods, molecular imaging and functional imaging such as perfusion imaging might be potential methods to help explain the biological significance of DWI metrics. Secondly, the number of patients was relatively small and statistical analysis could not be reliably applied to study different glioma subtypes. As well as studying a larger cohort, additional molecular biomarkers, including proteins and genes, should be included in future studies. Finally, treatment planning must include consideration of regions of glioma infiltration and study of the microenvironment in and beyond the peritumor region is important.
In conclusion, different DWI metrics fitted within different bvalue ranges (from low to ultra-high b values) could act with different efficacy as surrogate indicator for molecular expression or microstructural complexity in gliomas. Further studies which associate pathology or physiology with imaging performance are needed to better explain the biological meanings for these DWI parameters in gliomas.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board of Henan Provincial People's Hospital. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
YB, TL, LC, HG, WW, GZ, SL, HL, and MW designed the overall study, and performed the experiments. YB and MW analyzed imaging data. LW and LK performed histological studies. YB, TL, and MW interpreted the results, and wrote the manuscript. NR revised the language of the manuscript. All authors contributed to the article and approved the submitted version.