Structural network alterations in patients with nasopharyngeal carcinoma after radiotherapy: A 1-year longitudinal study

This longitudinal study explored the changed patterns of structural brain network after radiotherapy (RT) in patients with nasopharyngeal carcinoma (NPC). Diffusion tensor imaging (DTI) data were gathered from 35 patients with NPC at four time points: before RT (baseline), 0∼3 (acute), 6 (early delayed), and 12 months (late-delayed) after RT. The graph theory was used to characterize the dynamic topological properties after RT and the significant changes were detected over time at the global, regional and modular levels. Significantly altered regional metrics (nodal efficiency and degree centrality) were distributed in the prefrontal, temporal, parietal, frontal, and subcortical regions. The module, that exhibited a significantly altered within-module connectivity, had a high overlap with the default mode network (DMN). In addition, the global, regional and modular metrics showed a tendency of progressive decrease at the acute and early delayed stages, and a partial/full recovery at the late-delayed stage. This changed pattern illustrated that the radiation-induced brain damage began at the acute reaction stage and were aggravated at the early-delayed stage, and then partially recovered at the late-delayed stage. Furthermore, the spearman’s correlations between the abnormal nodal metrics and temporal dose were calculated and high correlations were found at the temporal (MTG.R and HES.L), subcortical (INS.R), prefrontal (ORBinf.L and ACG.L), and parietal (IPL.R) indicating that these regions were more sensitive to dose and should be mainly considered in radiotherapy treatment plan.


Introduction
Nasopharyngeal carcinoma (NPC) is a malignant tumor, and it is mostly found in Southern China and Southeast Asia (Chan, 2010;Tabuchi et al., 2011). Radiotherapy (RT) with or without adjuvant chemotherapy is the primary treatment for patients with NPC. However, the normal brain tissues surrounding the tumor are inevitably irradiated during cranial irradiation, causing brain abnormalities and cognitive decline. These abnormalities may compromise the quality of life of patients with NPC. Based on the pathophysiology of the side effects of RT, the time following RT can be classified into acute reaction period (days-weeks) (post-RT-AC), early-delayed period (1-6 months) (post-RT-ED), and late-delayed period (> 6 months) (post-RT-LD) (Lell, 2015). The RT-related brain changes are different during different periods but how the RT-related brain damage evolves over time is still unclear. Therefore, it is essential to further explore the temporal brain changes after completing RT which may facilitate clinical diagnosis and early intervention.
Recently, few cross-sectional or longitudinal studies have demonstrated that normal-appearing brain tissues underwent different changes at different post-RT periods in patients with NPC using various magnetic resonance imaging (MRI) analysis techniques Guo et al., 2018;Lv et al., 2019;Wu et al., 2020;Qiu et al., 2021). Specifically, our previous longitudinal studies found that the volumes of the gray matter  and white matter (WM) in bilateral temporal subfields  and bilateral hippocampal subfields (Lv et al., 2019) decreased over time after RT. In addition, the cross-sectional or longitudinal studies on cortical brain morphology revealed progressive RT-induced reduction in cortical volume, cortical thickness, and cortical surface area, mainly in the temporal, basal occipital, and basal frontal lobes Zhang et al., 2018). Aside brain morphological alteration, the WM microstructure changed after RT in patients with NPC (Wang et al., 2012;Xiong et al., 2013;Chen et al., 2015Chen et al., , 2020Duan et al., 2016;Leng et al., 2017Leng et al., , 2019Ding et al., 2018). Diffusion tensor imaging (DTI) is the only non-invasive MRI technique to assess brain white matter microstructure in vivo (Le Bihan et al., 2001). Most DTI studies adopted the regions of interest (ROI)-based analysis strategy to detect the microstructural changes in the temporal lobe of patients with NPC (Wang et al., 2012;Xiong et al., 2013;Chen et al., 2015). They found that diffusion metrics, such as FA and ADC in the temporal lobe, exhibited dose-related dynamic alterations over time after RT. Nevertheless, the ROI-based analysis is limited to specific regions and cannot reflect whole-brain changes. Recently, some studies investigated the changes in whole-brain WM at different post-RT periods by voxel-based analysis (Duan et al., 2016;Leng et al., 2017Leng et al., , 2019Ding et al., 2018). They found that RT-induced brain alterations were dynamic and extensive, and were not limited to the temporal lobe.
However, the voxel-based analysis cannot reflect the dynamic interaction of distinct brain regions. The graph theory analysis models brain connectivity as a network to assess the structural and functional brain organization (Sporns, 2011), offering an opportunity to better understand how the brain changes from a network perspective. The structural connectivity (SC) network is usually considered to be the physical substrate of the functional connectivity (FC) network. In patients with NPC, functional and structural brain network topology change after RT (Ma et al., 2016;Tian and Zhao, 2017;Qiu et al., 2018;Leng et al., 2019;Chen et al., 2020). For structural brain networks, a longitudinal DTI study reported that both global and local efficiencies, as well as the nodal topology, were altered in post-RT patients (Tian and Zhao, 2017). This study only investigated the difference between pre-RT and post-RT, but did not consider the different patterns of brain changes at different post-RT periods. Subsequently, a cross-sectional DTI study on three points (baseline, post-RT-ED, and post-RT-LD) found that structural topological properties were altered in the post-RT-ED but began recovering in the post-RT-LD (Chen et al., 2020). Nevertheless, in this cross-sectional study, the data with different post-RT durations were not from the same group of patients with NPC; the cohort effect could compromise the ability to detect the RT-induced brain alteration; the study did not investigate the acute reaction period which exhibits different side effect of RT when compared to the post-RT-ED and post-RT-LD periods. Inclusion of three post-RT periods will facilitate better understand the RT-related brain changed patterns over time.
Therefore, this work will adopt a longitudinal study with four time points (baseline, post-RT-AC, post-RT-ED, and post-RT-LD) to investigate the dynamic changes in structural brain network. Our cohort group included 35 patients with NPC, and each patient was followed up with four repeated scans: prior to RT, 0∼3, 6, and 12 months follow-up after the completion of RT. The topological properties of the structural network at the global, regional, and modular levels were calculated. Based on the analysis of these topological properties, the dynamic brain changes after RT and the relationship between these brain alterations and radiation dose were assessed.

Materials and methods Patients
Forty-three newly diagnosed treatment-naïve patients with NPC (aged 18-60) were initially enrolled. The inclusion criteria were as follows: right-handedness, no alcoholism or substance dependence, no high blood pressure, no diabetes, no brain tumors, no visible brain lesions, no history of cranial trauma, no history of any psychiatric or neurological disease, no current medications that may affect cognitive function, and

Treatment
All patients were treated with intensity-modulated radiotherapy (IMRT) (n = 32) or tomotherapy (TOMO) (n = 3), the details of which have been reported by previous studies (Sun et al., 2013;Tang et al., 2015). The prescribed regimen included a total dose of 68-70 Gy in 30-33 fractions at 2.12-2.33 Gy/fraction to the planning target volume (PTV) of the primary gross tumor volume (GTVnx), 60-70 Gy to the PTV of GTV of involved lymph nodes (GTVnd), 60-64 Gy to the PTV of the high-risk clinical target volume (CTV1), and 54-58 Gy to the PTV of the low-risk clinical target volume (CTV2). All patients received one fraction daily over a period of about 45 days, five consecutive days per week. Based on the guidelines defined by the 7th edition of the AJCC staging system for NPC, the patients with stage I to IIa disease received no chemotherapy, those with stage IIb received concurrent chemotherapy, and those with stages III to IVa-b received concurrent chemotherapy with/without neoadjuvant/adjuvant chemotherapy (Edge et al., 2010).

Follow-up procedure
To assess the dynamic alterations in structural brain network topology after RT, we repeatedly performed MRI scanning at the following stages for each patient: before initiation of RT (baseline), 0∼3 months (post-RT-AC), 6 months (post-RT-ED), and 12 months (post-RT-LD) after RT. Since the MRI data at each stage were acquired from the same group (35 patients), a longitudinal comparison strategy was performed to avoid potential bias due to cohort effect.

Data preprocessing and tractography
The data preprocessing included the following steps: (1) denoising the DW images using Marchenko-Pastur PCA (Veraart et al., 2016); (2) correcting the eddy current and head motion-induced distortion with an affine transformation; (3) skull stripping for the T1-weighted images and non-DW images (b = 0 s/mm 2 ) with FSL-Brain Extraction Tool (BET).
Whole-brain fiber reconstruction was performed for each diffusion data in native space using probabilistic tracking. Anatomically constrained tractography , seeding from the interface between grey matter and white matter, was used to achieve an anatomically plausible trajectory. A total of 10 million (M) seeding streamlines were initially generated and tracked. Finally, Spherical-deconvolution Informed Filtering of Tractograms (SIFT) (Smith et al., 2013) was performed to filter the streamlines from 10 to 1 M for improving the quantitative nature of whole-brain streamline reconstruction.
All these preprocessing steps and fiber tracking were accomplished within MRtrix3, 1 which is an open-source software package and includes scripts that interface with external packages, such as FSL 2 (Jenkinson et al., 2012). Figure 1 shows the flow chart of structural network construction, which is also accomplished within MRtrix3. First, 90 brain regions (nodes) were created for each participant with the automated anatomical labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002). Particularly, for each participant, the non-diffusion images (b = 0 s/mm 2 ) were co-registered to the corresponding T1-weighted images with an affine transformation. Meanwhile, the T1-weighted images were nonlinearly transformed to the Montreal Neurological Institute (MNI) space using the ICBM-152 brain template. Thereafter, these two transformations were inversed and combined into one transformation, which was applied to wrap the AAL from The flowchart of white matter structural network construction. (1) The non-diffusion image b0 was co-registered to the corresponding T1-weighted image with an affine transformation T b0→T1 for each participant.

Structural network construction
(2) The T1-weighted image was non-linearly transformed to the Montreal Neurological Institute (MNI) space using the ICBM-152 brain template, resulting in a non-linear transformation T T1→MNI . (3) These two transformations (T b0→T1 and T T1→MNI ) were inversed and combined into one transformation (T MNI→b0 ), which was applied to wrap the AAL from the MNI space to the native diffusion space of each participant. Thereafter, 90 participant-specific brain regions were created. (4) A whole-brain fiber reconstruction was performed for each subject in native space using the probabilistic tracking in conjunction with DTI. (5) Finally, the 90 × 90 symmetric connectivity matrix for each participant was constructed by calculating the mean FA values of streamlines that connect each region pair.
MNI space to the native diffusion space of each participant. Finally, the 90 × 90 symmetric connectivity matrix for each participant was constructed by calculating the mean FA values of streamlines that connect each node pair.

Structural network analyses
The global and regional network metrics, as well as the modular metrics, were calculated to characterize the topological properties of altered structural networks. All the following network metrics were calculated using GRETNA. 3 The global and regional network metrics The global metrics calculated in our study consisted of global efficiency (E glob ), local efficiency (E loc ), cluster coefficient (Cp), shortest path length (Lp), normalized cluster coefficient (γ), normalized characteristic path length (λ), and smallworldness (σ). For regional properties, the following two nodal 3 http://www.nitrc.org/projects/gretna/ metrics were considered: nodal efficiency (NE) and degree centrality (DC). The definition and interpretation of these network metrics can be referred to Rubinov and Sporns (2010).
To avoid both spurious connections and bias of a single sparse threshold, the area under the curve (AUC) under sparsity, ranging from 27 to 40% with an interval of 0.5% for each global and regional measures, was calculated for the following statistical analysis.

The modular metrics
With regards to the modularity analysis, the total number of modules and the associated module membership of nodes were optimized by maximizing modularity Q, the detailed definition and interpretation of which can be referred to Newman, 2006. Particularly, the Louvain algorithm in the Brain Connectivity toolbox 4 was used to optimize the Q-value under varying sparsity, ranging from 0.05 to 0.3. Generally, a Q-value > 0.3 indicated a strong modular structure (Fortunato and Barthelemy, 2007;Hilger et al., 2017). Given that the number of modules and their membership varied across different sparsity, sparsity was fixed at 0.1 (Q = 0.39) to obtain a quite reasonable and consensus modularity partition, as shown in Figure 2. Thereafter, to assess the modular segregation, the within-module connectivity and between-module connectivity, which are defined as the strengths of edges within a single module and between a pair of modules, respectively, were calculated.

Statistical analysis
The one-way repeated measures analysis of variance (ANOVA) was used to compare four groups at the baseline, and three follow-up stages (post-RT-AC, post-RT-ED, and post-RT-LD) after RT for all of these metrics. All the measured data satisfied the assumptions of normality and homogeneous variance. The assumption of sphericity was violated for the regional network metrics in several brain regions, where the Greenhouse-Geisser method was used to correct the sphericity. Thereafter, a post-hoc analysis (multiple comparisons) was performed by using paired t-test to compare each pair within the four groups. When the differences between the paired observations did not follow a normal probability distribution, Wilcoxon Signed-Rank test, which is a nonparametric equivalent of the paired t-test, was used instead. Finally, false discovery rate (FDR) correction was used for multiple comparisons.
In addition, dose-response analysis was performed by calculating the Spearman's rank correlation coefficient (r-value) of the association between the abnormal nodal metrics and the radiation dose of ipsilateral temporal lobe.

Results
The global analysis    Table 1.
The regional analysis  Table 2, and the relevant information of 90 regions from the AAL atlas together with corresponding abbreviations are listed in Supplementary Table 1.
For nodal efficiency, 18 regions had significant differences among the four groups. They were found in the prefrontal

The modularity analysis
Four modules were identified according to the mean network matrix of whole patients at baseline. The detailed information of four modular networks are summarized in Table 3. The modularity analysis results are shown in Figure 2.
For the within-module connectivity strengths, only module 2 showed a statistically significant difference among the four stages. The module 2 comprises 21 regions, including the prefrontal lobe, frontal lobe, and parts of the parietal lobe. It has a high overlap with the default mode network (DMN), and is related to the normal cognitive and emotional functions. Specifically, the connectivity strength within module 2 at the  post-RT-ED was significantly lower than those at baseline and post-RT-AC. In addition, the within-module connectivity strengths for each module showed a tendency of progressive decrease at post-RT-AC and post-RT-ED and then exhibited a recovering trend at the post-RT-LD, when compared to baseline, but no significant difference was found. For between-module connectivity strengths, no significant differences were found for each pair of modules. However, all the between-module connection strengths, except for that between modules 1 and 3, showed a tendency of sustained and progressive decrease at post-RT-AC and post-RT-ED and then exhibited a recovering trend at post-RT-LD, when compared to baseline. Figure 5 shows the correlations between the abnormal nodal parameter metrics (NE and DC) and the radiation dose of ipsilateral temporal lobe. The mean and/or maximum doses of ipsilateral temporal lobe were correlated with the changed NE and DC in several regions which were distributed in temporal, subcortical, prefrontal and parietal. In brief, the changed NE and DC were positively correlated with mean and/or maximum dose at ACG.L, INS.R, HES.L, and IPL.R, and negatively correlated at ORBinf.L and MTG.R. In addition, more brain regions exhibited a significant dose correlation with NE and DC at late-delayed period; more brain regions were correlated to the mean dose than the maximum dose. The spearman's correlations (r-values) with significant differences (p-values) for the mean dose and the maximum dose were shown in Supplementary Tables 2, 3, respectively.

Discussion
To our knowledge, this study is the first longitudinal cohort study to monitor the RT-induced alterations of brain structural network in patients with NPC after RT. The DTI probabilistic tractography and graph theoretical approach were used to assess RT-related brain changes at the global, local, and modular levels; the following findings were obtained: (1) E loc shows a significant difference among four stages. Both E glob and E loc show a tendency of progressive decrease at post-RT-AC and post-RT-ED and a partial recovery at post-RT-LD. (2) Except for the ORRinf.L and DCG.L, all other regions exhibited significant  reductions in the nodal efficiency and degree centrality at post-RT-AC and post-RT-ED, and most of these regions showed a partial or full recovery at post-RT-LD.
(3) The within-module connectivity strength of modular 2 exhibited significant and progressive decrease at post-RT-AC and post-RT-ED, compared to baseline, and showed a partially recovering trend at post-RT-LD. All the within-and between-module connectivity strengths, except for that between modules 1 and 3, showed a tendency of sustained and progressive decrease at post-RT-AC and post-RT-ED. Thereafter, a recovering trend at post-RT-LD was exhibited. All these findings imply that the brain injures begin at post-RT-AC, are aggravated at post-RT-ED, and undergo brain reorganization at the post-RT-LD. (4) The temporal irradiation dose was significantly correlated to the altered nodal parameters at the temporal (MTG.R and HES.L), subcortical (INS.R), prefrontal (ORBinf.L and ACG.L) and parietal (IPL.R), which suggests that these regions were more sensitive to dose and should be paid more attention during RT treatment plans. The global network analysis revealed that the structural brain network possessed small-world properties (λ ≈ 1, γ > 1, and σ > 1), at baseline and all three follow-up stages (post-RT-AC, post-RT-ED, and post-RT-LD). These results illustrate that the small-world networks are relatively robust to the changes of brain white matter (He et al., 2009;Colombo, 2013;Xu et al., 2017). For the presented global measures, only E loc had statistically significant difference among the four groups. E loc represents the efficiency of information exchange within a local subnetwork or among adjacent regions (Jiang et al., 2020). Reduced E loc in a structural brain network may arise from the RT-associated injures of the fiber tracks (e.g., demyelination and axonal damage) (Nazem-Zadeh et al., 2012;Qiu et al., 2021). The results of significant decrease in E loc at post-RT-ED and post-RT-LD were compatible with findings of prior fMRI studies reporting lower efficiency of information transfer after RT (Ding et al., 2018;Leng et al., 2021). Notably a significant decrease in E loc firstly occurred 6 months after RT in our structural network study, later than the significant abnormalities in global properties of functional networks (<6 months) (Leng et al., 2021). These findings were plausible because brain function might be more vulnerable or sensitive to attack (Karim et al., 2017). In addition, both E glob and E loc showed a tendency of progressive decrease at post-RT-AC and post-RT-ED and partial recovery at post-RT-LD, although this trend was not statistically significant. These inconspicuous changes in trend of E glob and E loc may explain the inconsistent and unstable results from previous studies. Some DTI studies found a gradual and irreversible white matter damage (Nagesh et al., 2008;Welzel et al., 2008;Ding et al., 2018), whereas other groups found that the DTI metrics decreased in the early stage but partially recovered later (Wang et al., 2012;Xiong et al., 2013;Chen et al., 2015).
The significant alteration of the nodal parameters (nodal efficiency and degree centrality) among the four stages was mainly located in the temporal, frontal, prefrontal, parietal, and subcortical regions. Most of these regions showed a progressive decrease during 0-6 months post-RT and a partial or full recovery 12 months post-RT. This result may indicate that the structural brain reorganization mainly occurred in the late-delay stage, which is generally consistent with findings of previous studies (Wang et al., 2012;Xiong et al., 2013;Duan et al., 2016;Chen et al., 2020). However, some regions, including MTG.L&R, HES.L, PreCG.L&R, SFGdor.L, and SMA.L, exhibited a sustained decrease without recovering tendency within 1 year after RT, which may be due to two reasons: vulnerability of these regions to radiation causing an irreversible damage and the need of these regions for a longer recovery time (>12 months), which cloud not be observed in this 1-year longitudinal study after RT. The bilateral temporal lobes, including MTG.L&R and HES.L, exhibited decreased nodal parameters without recovering tendency over the time after RT. This observation was not surprising because the temporal lobe is often inside the target volume and inevitably receives high-dose radiation and may suffer from severe injury. Late-delayed temporal injuries have been well documented as irreversible, and sometimes presented as necrosis of temporal lobes on routine medical imaging examinations (Mao et al., 2014;Lv et al., 2019). Additionally, the nodal parameters showed significant changes in the prefrontal, frontal, and parietal regions, which were outside the irradiation field. Previous TBSS analysis (Duan et al., 2016) revealed that the fractional anisotropy values were significantly lower in the frontal, parietal, and occipital WM after RT. A previous VBM study found a reduced GM volume in the frontal and parietal cortices (Lv et al., 2014). Altogether, the changes in nodal parameters in the prefrontal, frontal, and parietal regions may arise from the degeneration of associated white matter fibers or radiationinduced disruption of the blood brain barrier (BBB) (van Vulpen et al., 2002). Notably, the increased nodal parameters in the ORBinf.L and DCG.L might act as a compensatory mechanism that maintains normal cognitive function. The subcortical regions, including the THA.L and INS.R, exhibited a different changing pattern, when compared with baseline. Specifically, the INS.R shows a "decrease-decrease-recover" pattern after RT for both nodal efficiency and degree centrality, whereas THA.L begins to decrease 12 months post-RT (post-RT-LD) for nodal efficiency. The alteration of structural brain network in the insular and thalamus is probable, given that both regions are parts of the paralimbic system that are sensitive to irradiation. In addition, these findings are consistent with those of previous studies (Ding et al., 2018;Qiu et al., 2018;Yang et al., 2019;Zhang et al., 2020;Nan et al., 2022), which reported functional and/or morphological changes in the thalamus and insula.
The dose-correlation analysis shows the nodal parameters (NE and DC) had a positive correlation with temporal dose at ACG.L, INS.R, HES.L, and IPL.R, which may be due to the compensatory change in structural brain network that interconnects these regions. Whereas the nodal parameters had a negative correlation with temporal dose at ORBinf.L and MTG.R, indicating that a higher dose reduces the information transfer efficiency to these regions. In addition, through acute reaction stage to late-delayed stage, the number of significant dose-correlation brain regions increased. This finding suggests that the dose effect on brain change is more notable at the latedelayed stage. Furthermore, some brain regions were correlated to the mean dose and/or maximum dose which illustrates that both the mean dose and the maximum dose should be considered for the protection of normal organs.
This study explored the changed patterns of structural modularity over time after RT in patients with NPC. We found that the connectivity strength within module 2 at the post-RT-ED were significantly weaker than those at baseline and post-RT-AC, indicating radiation-induced disruption of topological organization of module 2. The module 2 mainly includes the prefrontal lobe, frontal lobe, and parts of the parietal lobe. The areas of module 2 and the DMN have a large overlap, and the DMN is associated with normal cognition and emotion (Alves et al., 2019). Moreover, the module 2 includes the medial prefrontal lobe, whereas module 4 includes the left side of the temporal lobe and parietal lobe. Several fiber bundles run between the medial prefrontal lobe and temporal lobe, which is highly related to memory processing (Vertes et al., 2007). The decreased connectivity strengths within module 2 and between modules 2 and 4 at the acute-and early delayed stages may be due to the damage of axonal fiber tracts between the medial prefrontal lobe and temporal lobe. These findings support the psychological disorders, cognitive dysfunction, and mood disorders commonly found in patients with NPC after RT (Tang et al., 2012;Mo et al., 2014;Wu et al., 2014). In addition, a "decrease-decrease-partially recovery" pattern was observed for the connectivity strengths within each module and between each pair of modules, although no significant alterations were found except for connectivity strengths within module 2. These observed results were roughly consistent with the findings on nodal parameters, further implying that the brain undergoes recovery and reorganization of structure to a certain extent at the late-delayed stage.
Despite the merits of this longitudinal study, several limitations were identified. First, the 1-year follow-up was insufficient to monitor all the dynamic changes in structural network properties after RT over time. A longer period, ranging over several years, should be considered to understand whether the injured structural network topology will eventually recover to "baseline" with time. Second, this study included 35 patients with NPC; this sample size was not large enough. A larger cohort size of patients with NPC is needed to provide more reliable statistical results and to accurately reveal the dynamic changing pattern of structural brain network after RT. Third, the relationship between the alterations in structural brain network and cognitive decline were not explored because of incomplete neurocognitive outcomes.

Conclusion
The follow-up data were used to track the dynamic changes in structural brain network after RT in patients with NPC. Our study found that the radiation-induced alterations in topological properties mainly began at the acute reaction stage, were aggravated at the early delayed stage, and then partially recovered at the late-delayed stage. The dynamic change patterns of topological properties facilitate to better understand how the radiation-induced brain injuries evolves over time and the early detection of radiation-induced changes in normal-appearing brain tissue to improve patient survival. In addition, a dose-correlation alteration was found in the temporal (MTG.R and HES.L), subcortical (INS.R), prefrontal (ORBinf.L and ACG.L), and parietal (IPL.R), indicating that these regions were more sensitive to dose and should be mainly considered in radiotherapy treatment plan.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement
The studies involving human participants were reviewed and approved by the Sun Yat-sen University Cancer Center. The patients/participants provided their written informed consent to participate in this study.

Author contributions
YY, YF, GF, and JL contributed to design of the study and data collection. XZ, JPa, PX, and YL were responsible for experimental implementation. YL, PX, CY, JPe, and XL performed the data analysis. XZ, JPa, GF, YY, and YF contributed to the manuscript writing. All authors read and approved the published version of the manuscript.