Stability evaluation of the PROSPECT model for leaf chlorophyll content retrieval

The radiative transfer model, PROSPECT, has been widely applied for retrieving leaf biochemical traits. However, little work has been conducted to evaluate the stability of the PROSPECT model with consideration of multiple factors (i.e., spectral resolution, signal-to-noise ratio, plant growth stages, and treatments). This study aims to investigate the stability of the PROSPECT model for retrieving leaf chlorophyll (Chl) content (Cab). Leaf hemispherical reflectance and transmittance of oilseed rape were acquired at different spectral resolutions, noise levels, growth stages, and nitrogen treatments. The Chl content was also measured destructively by using a microplate spectrophotometer. The performance of the PROSPECT model was compared with a commonly used random forest (RF) model. The results showed that the prediction accuracy of PROSPECT and RF models for Cab did not produce significant differences under varied spectral resolutions ranging from 1 to 20 nm. The ranges of the relative root mean square errors (rRMSE) of the PROSPECT and RF models were 12%-13% and 11.70%-12.86%, respectively. However, the performance of both models for leaf Chl retrieval was strongly influenced by the noise level with the rRMSE of 13-15.37% and 12.04%-15.80% for PROSPECT and RF, respectively. For different growth stages, the PROSPECT model had similar prediction accuracies (rRMSE = 9.26%-12.41%) to the RF model (rRMSE = 9.17%-12.70%). Furthermore, the superiority of the PROSPECT model (rRMSE = 10.10%-12.82%) over the RF model (rRMSE = 11.81%-15.47%) was prominently observed when tested with plants growth at different nitrogen treatment levels. The results demonstrated that the PROSPECT model has a more stable performance than the RF model for all datasets in this study.


Introduction 
The photosynthetic capacity of crops is often used as an essential indicator to evaluate the potential crop yield [1,2] . Chlorophyll (Chl) is one of the primary biochemical traits related to the process of photosynthesis [3] . Thus, accurate measurement of leaf Chl content (C ab ) is critical for monitoring crop growth and health status in precision agriculture and crop breeding [4] . Since Chl mainly absorbs red and blue-violet light in the visible spectrum (VIS, 400-750 nm), it shows troughs in the reflectance curve of the red and blue-violet light regions which makes possible the retrieval of C ab from the spectrum of the leaf [5] . Some researchers have developed various non-destructive methods to estimate C ab from leaf reflectance or transmittance properties, including physically-based and empirically-based methods [6] . Compared with empirical methods, the physically-based methods evolved from radiative transfer models (RTMs) are versatile in predicting leaf biochemical components [7] . PROSPECT, one of the most popular RTMs, has been widely used to retrieve C ab and is based on rigorous physical and mathematical theories [8,9] . The refractive index and the maximum incident angle of the leaf in the PROSPECT model are applied to establish the relationship between leaf biochemical traits (Chl content, dry matter content, equivalent water content, carotenoid content) and leaf surface reflectance and transmittance [10,11] . The PROSPECT model has also been significantly improved and applied to leaf biochemical component retrieval for numerous plant species [6,7,12,13] . Li et al. [12] optimized this model by coupling the PROSPECT model with the continuous wavelet transform (CWT) to suppress the surface reflectance effect and improve the accuracy of chemical composition retrieval. Additionally, this model has also been used to quantify the carotenoid, dry matter, and anthocyanins [14][15][16] . In view of the difficulty of accurately and non-destructively estimating the leaf protein and nitrogen content in a timely manner, researchers have further proposed the assessment of protein and lignin content of fresh and dry leaves based on the PROSPECT model [17][18][19] . Further, several studies have empirically correlated the C ab with spectral features using machine learning algorithms, such as random forest (RF), partial least squares (PLS), support vector machine (SVM) and their feasibility and accuracy have also been demonstrated [20][21][22][23][24][25][26][27] . Spectral indices (e.g., vegetation indices (VI), and normalized difference vegetation indices (NDVI)) and derivative spectra were found to be the most preferable features for the estimation of the C ab of various crops [3,25,[28][29][30][31][32][33] . Most of the aforementioned, empirically based methods are only used for estimating C ab at a single growth stage, and their performance is highly influenced by the spectral resolutions and signal-to-noise ratio (SNR) of the sensors. Unlike machine learning algorithms, which require retraining for each plant type or growing conditions, the PROSPECT model has shown superior performance in estimating C ab and can be used to fit experimental data of a wide range of plant types. The stability of these models in different growth stages or using datasets with varying spectral resolutions has not been investigated. Therefore, the RF model, a typical empirically based model, was chosen as an example to evaluate the stability of the PROSPECT model by comparison.
Generally, the quality of spectral data could be easily affected by environmental factors as well as sensor performance. At the same time, the inversion of the PROSPECT model is sensitive to the noise level of the dependent variable (e.g., reflectance). Therefore, it is worthwhile to investigate the performance of the PROSPECT model to retrieve C ab with the reflectance at different spectral resolutions and noises levels acquired at different crop growth stages.
This study aims to investigate the performance of the PROSPECT model for estimating C ab of oilseed rape leaves under different growing conditions and to compare its performance with the RF model. Varied spectral resolutions were used and noise was added to the dataset to simulate real-world problems and to ensure that the model can generalize well for the dataset collected using different sensors under varied environmental conditions. The performance of both models under varied growth stages and nitrogen treatment levels was also explored in this study in order to comprehensively evaluate the stability of the PROSPECT model.

Experimental design
Two-year oilseed rape experiments were conducted in 2017/2018 and 2018/2019 during the winter seasons at the Agricultural Research Station of Zhejiang University, Hangzhou, China (30°18′26′′N, 120°4′29′′E). In the two experiments, the oilseed rape of cultivar "ZD630" was first sown in the prepared seedbeds using high fertility soil in mid-October, and the seedlings at the three-leaf stage were transplanted to plastic pots with a volume of 3 L. During the first experiment, three nitrogen (N) levels, including a control treatment (no added nitrogen: N0), 311.74 g/pot (N1) and 623.48 g/pot (N2) were used to obtain various leaf C ab content samples. In the second experiment, three nitrogen levels, including a control treatment (no added nitrogen: (N0), 311.74 g/pot (N1), and 623.48 g/pot (N2) were also applied. For two experiments, the nitrogen was applied in three batches. The first batch was a basal fertilizer (50%), which was added into the soil before transplanting, the second batch (30%) was added at the seeding stage, and the third (20%) at the bolting stage. The fertilizer was applied in the form of urea (CH 4 N 2 O). Calcium magnesium phosphate (P 2 O 5 = 478.13 g/pot) and potassium chloride (K 2 O = 239.06 g/pot) were added as supplementary nutrients.

Data collection and processing 2.2.1 Measurement of hemispherical reflectance and transmittance
The hemispherical reflectance (R) and transmittance (T) were measured using an ASD FieldSpec4 Hi-Res Spectroradiometer (Analytical Spectral Devices, Boulder, CO, USA), with a spectral range of 350-2500 nm, equipped with an RTS-3ZC integrating sphere (Analytical Spectral Devices, Boulder, CO, USA). The sampling interval of this spectrometer is 1.4 nm at 350-1000 nm and 1.1 nm at 1000-2500 nm with a spectral resolution of 3 nm in the visible and near infrared (VIS-NIR) spectrum and 8 nm in the shortwave-infrared (SWIR) spectrum.
The white plug (polytetrafluoroethylene) and light trap (metal coated with black paint) with a reflectivity of about 99% and 0, respectively, and 10 W collimated ASD CL-10 halogen light source were assembled in this integrating sphere.
The measurement principle of hemispherical reflectance and transmittance is shown in Figure 1. The radiance of the sample was measured immediately after dark current and white plug radiance acquisition. The hemispherical reflectance and transmittance of the leaf were determined using Equation (1) and Equation (2). R and T of each sample were averaged with three repeated measurements. In this study, 240 and 316 leaf samples were collected in 2017/2018 and 2018/2019, respectively.
where, R s is the hemispherical reflectance of the leaf; I s1 is the reflected radiance of the leaf; I d1 is the radiance of the light trap; I r1 is the radiance of the white plug, and R r is the reflectance of the white plug (99%).
where, T s is the hemispherical transmittance of the leaf; I s2 is the transmitted radiance of the leaf; I d2 is the radiance of the light trap; I r2 is the radiance of the white plug and R r is the reflectance of the white plug (99%).
a. Reflectance measurement b. Transmittance measurement Figure 1 The experimental setup for directional-hemispherical reflectance (a) and transmittance (b) measurements of leaf samples Generally, the R and T at the spectral region of 400-2500 nm were used to retrieve various leaf biochemical traits [34] . Since this study aims to quantify the leaf C ab , the spectral range of 400-1000 nm was taken into account [35] . Public dataset (ANGER 2003) was used to prove our claim of considering this spectral region for the leaf C ab estimate. First, the resolution of raw spectrum was reduced by binning (3 nm, 5 nm, 10 nm, 15 nm and 20 nm), which is a common method for changing spectral resolutions in order to study the effect on both models [36] . Then, by adding varied proportions of Gaussian white noise (1%, 2%, 3%, 5%, and 10%) into the raw spectrum, a stability evaluation of the noise addition on the two models was also carried out ( Figure 2). Finally, datasets of varied growth stages and nitrogen treatment levels were also considered for the stability assessment of the two models.

Measurement of leaf C ab
After measurements of hemispherical reflectance and transmittance, each fresh leaf sample was cut into 12 round disks of a diameter of 0.85 cm with three replicates. Four round disks were then put into 9 mL of ethanol solution (95%, v/v) for more than 24 hours for Chl extraction [37] . The absorbance of Chl at specific wavelengths (A 649 and A 665 ) was measured with a microplate spectrophotometer (Epoch 2, BioTek Instruments, Winooski, VT, United States). The leaf C ab was calculated according to the following formulas [38] : Chl a = 13.95×A 665 -6.88×A 649 (3) Chl b = 24.96×A 649 -7.32×A 665 (4) Chl (a+b) = Chl a + Chl b (5) The content (μg/cm 2 ) of Chl a and Chl b was derived from the following equation: where, V is the volume of the ethanol solution and S is the area of the round disks per leaf sample, C ab is leaf Chl content. Datasets under varied growth stages and nitrogen treatment levels were collected in order to obtain a wide range of leaf C ab . One-way Analysis of Variance (one-way ANOVA) was then used to test all measured samples to evaluate the effectiveness of nitrogen treatment levels (N0, N1, and N2) on leaf C ab . The effectiveness of nitrogen treatment was considered significant at p < 0.05. Before one-way ANOVA, the measured leaf C ab value was examined by the Shapino-Wilk normality test method to confirm to Gaussian distribution. Multiple comparisons and the Tukey method were selected to correct error after passing the normal distribution test. GraphPad Prism 8 (GraphPad Software, Inc., San Diego, CA, USA) was used for the statistical analysis, and PROSPECT and RF models were performed with MATLAB 2018 (MathWorks, Inc., Natick, Ma, USA).

Models description 2.3.1 PROSPECT model
The PROSPECT model describes the structural properties of the leaf in detail and establishes the relationship between the biochemical components of the leaf and the reflectance and transmittance of the leaf based on strict physical and optical assumptions [10,34,39] .
The PROSPECT model was initially developed based on the plate model in which the leaf was assumed to be composed of homogeneous elementary layers with the same refractive index and absorption coefficient everywhere [40] . In fact, the interior of leaf is not uniform. With continuous advancements of research in this field, the leaf internal structure was described as N homogeneous elementary layers and N-1 air layers which have a specific absorption coefficient and refractive index for each layer, thus forming the preliminary PROSPECT model that has gradually increasing in use [10] . The advantages of the PROSPECT model include its ability to provide accurate spectral estimates with limited input parameters.
In the PROSPECT model, the reflectance and transmittance that are calculated involve leaf Chl content (C ab , expressed in µg/cm 2 ), leaf dry matter content (C m , expressed in g/cm 2 ), equivalent water content (C w , expressed in g/cm 2 ), carotenoid content (C xc , expressed in µg/cm 2 ), and structure parameters (Nj, number of layers). Minimizing the difference between the measured spectrum and the modeled spectrum is key to optimizing the estimates of these input parameters. A program could be used to minimize the fitting between the PROSPECT model predicted reflectance and transmittance (R mod and T mod ) and the measured reflectance and transmittance (R meas and T meas ), regardless of the value of the initial input parameters [9,39,41] . Therefore, merit function J was defined to describe the difference between the predicted reflectance and transmittance and the measured reflectance and transmittance at each wavelength: where k (λ) is the total absorption coefficient of all biochemical components at the wavelength λ; K spe ,i (λ) is the absorption coefficient of biochemical component I; R mod and T mod are the predicted reflectance and transmittance; R meas and T meas are the measured reflectance and transmittance; C i,j is the content of biochemical component i in leaf j; N j is the structure parameter of leaf j; n is the number of biochemical components, and p represents the number of wavelength intervals between 400 and 1000 nm. As the most popular version, the PROSPECT-5 model was selected to evaluate stability for leaf C ab retrieval in this study.

Random forest model
The RF model is a common supervised machine-learning model used for regression and classification analysis [42,43] . It is characterized by the use of a resampling technique to randomly extract samples from the original training dataset to generate a new training set [44] . The regression trees were established with randomly selected samples and the predicted results were obtained by voting. Finally, the result of each regression tree was integrated and generated simultaneously to complete the final prediction. The RF model has been widely used to predict Chl fluorescence, biomass, leaf area index (LAI), and leaf C ab in the literature [45][46][47][48][49][50] . Therefore, the RF model was selected as a typical empirical model to predict leaf C ab . The dataset was divided into the training dataset (70%) and the test dataset (30%) for the establishment of the RF regression model. The training dataset and the test dataset were selected by a systematic sampling method in this study, which was easy to operate and can avoid random errors.

Models evaluation
In order to evaluate the performance and predictive ability of the PROSPECT and RF models, the prediction accuracy of the two models was tested. The coefficient of determination (r 2 ), root mean square error (RMSE), and relative root mean square error (rRMSE) were employed to estimate the performance of each model. The determination coefficient indicates the fitting of the predicted values and measured values. The RMSE quantizes the deviation between the predicted values and the measured values. Correspondingly, the rRMSE is obtained by taking the RMSE divided by the mean value of all samples. The smaller the rRMSE, the better the performance of the model. The calculation of these three performance evaluation indicators is as follows:  (11) where, ŷ is the predicted leaf C ab value calculated by the prediction models; y i is the measured leaf C ab value; y is the average value of all measured leaf C ab , and n is number of samples measured in this study.

Wavelength selection and validation
The accuracies of leaf Chl retrieval using public dataset ANGER 2003 (OPTICLEAF-database (ipgp.fr)) at two different spectral regions are shown in Figure 3. It is obvious that using the small spectral region (400-1000 nm) can provide the same accuracy as full spectral regions (400-2500 nm).
To reduce the computational complexity and eliminate unnecessary information, the spectral range of 400-1000 nm was considered for subsequent analysis and verified by PROSPECT model.

Distribution of leaf C ab
According to the results of leaf C ab distribution which are shown in Figure 4, it can be demonstrated that with an increase in nitrogen treatment levels, average leaf C ab increased proportionally. The average leaf C ab increased at the beginning and then decreased as growth stages continued. This is because four-leaf growth a. 400-1000 nm b. 400-2500 nm Figure 3 Results of leaf C ab retrieval using ANGER 2003 dataset. Leaf C ab retrieval result at 400-1000 nm (a); leaf C ab retrieval result at 400-2500 nm (b).
stage fertilizer was applied after the first measurement and before the second measurement. Therefore, leaf C ab decreased synchronously as growth stages lasted and fertilizer decreased in the soil, so that leaf C ab increased at the beginning of the four-leaf growth stage.
These findings hold important potential for regulating crop growth and field management.
The results of a one-way ANOVA analysis under varied nitrogen treatment levels and growth stages are shown in Figure 5. With increases in growth stages, significant differences among varied nitrogen treatment levels become more apparent, especially for N0 and N2. These results indicated that varied nitrogen treatment levels showed an obvious difference in leaf C ab and confirmed that nitrogen treatments were effective in this experiment.

Effects of spectral resolutions on model stability
Leaf C ab retrieval results of different spectral resolutions are summarized in Figure 6. The results showed that the PROSPECT model was more stable compared to the RF model, with an r 2 range of 0.55-0.56. In contrast, the RF model obtained a larger variation in outcomes when the spectral resolution was changed and r 2 ranged between 0.37 and 0.68, indicating that there was weak linear relationship between the predicted and measured values for the RF model. The variation in the rRMSE was relatively small for the PROSPECT model (12%-13%) and the RF model (11.69%-12.86%), which confirmed the robustness of both models for the retrieval of leaf C ab from leaf R and T with varied spectral resolutions.
With decreases of the spectral resolution, the prediction errors of PROSPECT and RF models showed minor changes. These findings indicated that the performance of both models was not sensitive to varied spectral resolutions.

Effects of noise addition on model stability
The performance of the PROSPECT and RF models was also compared under the effects of different noise ratios. The results (Figure 7) showed that the prediction accuracy of the two models decreased significantly with the increase of noise ratio. For the PROSPECT model, the rRMSE increased from 13.00% to 15.37%, and the rRMSE of the RF model increased from 12.04% to 15.80%, as shown in Figure 7. It can be concluded that both models were severely affected by noise interference, while the PROSPECT model was less affected than the RF model. According to the linear correlation analyses between the predicted and measured values which were presented in Figure 7, the r 2 obtained from the PROSPECT model (r 2 = 0.46-0.56) was relatively more consistent than the RF model (r 2 = 0.68-0.01). A stronger linear relationship between the predicted and the measured values was found in the PROSPECT model under varied noise addition.  Figure 7 Results of leaf C ab retrieval using the PROSPECT (a, c, e, g, i, k) and RF (b, d, f, h, j, l) models with different noise ratios: no noise (a, b); 1% noise addition (c, d); 2% noise ratio (e, f); 3% noise ratio (g, h); 5% noise addition (i, j) and 10% noise ratio (k, l).

Effects of varied growth stages on model stability
Based on the results of one-way ANOVA analysis, it was confirmed that different growth stages and nitrogen treatment levels showed significant differences on leaf C ab . The effects of these two factors on the PROSPECT and the RF models were explored to evaluate the prediction accuracy of leaf C ab .
As for varied growth stages, the results showed that the PROSPECT (rRMSE = 9.26%, Figure 8a) and the RF (rRMSE = 9.17%, Figure 8d) models exhibited the best performance at the three-leaf growth stage.

Effects of nitrogen treatment levels on model stability
As shown in Figure 9, the prediction accuracy of the PROSPECT and RF models showed different patterns under three nitrogen treatment levels. The performance of the PROSPECT (rRMSE = 11.34%, Figure 9a) and RF (rRMSE = 11.81%, Figure 9d) models was similar at N0. As for N1 and N2, the PROSPECT model (N1: rRMSE = 10.10%, Figure 9b; N2: rRMSE = 12.82%, Figure 9c) had better performance than the RF model (N1: rRMSE = 12.78%, Figure 9e; N2: rRMSE = 15.47%, Figure 9f). The PROSPECT model had higher prediction accuracy and stability under three nitrogen treatment levels, while there was a larger variation in the prediction accuracy of the RF model. In particular, the prediction accuracy for N2 was very low. All the above results showed that the performance of both models was affected by varied nitrogen treatment levels, but the degree of impact for them was different. The prediction capability of the PROSPECT model under different nitrogen treatment levels was more consistent indicating a robust prediction capability compared to the RF model.

Discussion
This study evaluates the stability of PROSPECT and RF models for leaf C ab estimation of oilseed rape grown in a two-year experiment. It also confirmed the advantages of the PROSPECT model for leaf C ab retrieval with complex data sources. Leaf C ab is an essential pigment for monitoring the growth and health status of a plant. In field management, different growth stages and fertilizer treatments can affect leaf C ab as well. According to the results of leaf C ab distribution which are shown in Figure 4, it can be demonstrated that with an increase in nitrogen treatment levels, average leaf C ab increased proportionally. The average leaf C ab increased at the beginning and then decreased as growth stages continued. This is because four-leaf growth stage fertilizer was applied after the first measurement and before the second measurement. Therefore, leaf C ab decreased synchronously as growth stages lasted and fertilizer decreased in the soil, so that leaf C ab increased at the beginning of the four-leaf growth stage. These findings hold important potential for regulating crop growth and field management.
Although there are some empirical methods related to C ab retrieval, the empirical models still show some disadvantages, e.g., poor transferability, the need for large sample inputs, and difficulty of interpretation.
With the development of physically-based models, the PROSPECT model has become a popular means of leaf biochemical component retrieval.
This study focused mainly on the stability evaluation of the PROSPECT model for leaf C ab retrieval of oilseed rape. With the rapid development of optical technology in the field of plant phenotypes, hyperspectral technology plays an important role in modern farmland management. Various hyperspectral devices with different spectral resolutions obtain plant spectral signals for further analysis of the growth status of crops. Differences in spectral resolution are one of the most important factors affecting the quality of the raw spectrum data. Thus, this article first evaluates the stability of the PROSPECT model under varied spectral resolutions. As the results showed in Figure 6, sensors with lower spectral resolution can also be developed for the retrieval of leaf C ab and field management. The potential of reducing spectral resolution to retrieve oilseed rape leaf C ab for both models was confirmed, which affirmed the feasibility of lower spectral resolution sensors to retrieve leaf C ab and provide a cost-effective method for future field management. As for noise addition, both models showed an increased tendency in prediction error, which suggested that noise addition had a great impact on the stability of both models (Figure 7). In view of this, more attention should be given to minimizing the impact of noise in experiments. In addition, it was demonstrated that the PROSPECT and RF models showed similar estimation accuracy in varied growth stages ( Figure  8).
In addition, for varied nitrogen treatment levels, the PROSPECT model showed more stable performance than the RF model ( Figure 9). According to the distribution of leaf C ab shown in Figure 4 and the inversion results in Figure 8 and figure 9, it can be seen that the PROSPECT model has stable performance for both discrete data points and concentrated data points. However, for the RF model, the prediction accuracy of concentrated data points was better than that of discrete data points. This allows a significant reduction of the field management sensors and equipment cost when compared with high spectral resolution multispectral sensors.
In addition, devices with high signal-to-noise ratio must be a priority, which can help ensure the quality of the datasets. For varied growth stages and different soil nutrient conditions, the advantages of the PROSPECT model will be more obvious. With reference to these findings, it is possible to find a more suitable leaf biochemical components prediction method after analyzing the composition of the acquired datasets. Investigating the stability of two representative models provides support for researchers to select appropriate models for leaf C ab retrieval and inspires new ideas for the development of portable field sensors.
An advantage of this study is that the stability of the PROSPECT model was evaluated using two-year experiments compared with the RF model under multiple factors such as varied spectral resolutions, noise addition, growth stages, and nitrogen treatment levels. These results reveal that it not only contributes well to the selection, application and development of phenotypic sensors in oilseed rape, but also is in support of the selection of leaf biochemical components prediction methods.
Overall, the stability assessment of leaf biochemical components prediction models with consideration for multiple scenarios is promising and necessary for plant phenotypic research. However, the dataset used in this study was still limited, and more oilseed rape cultivars and crops should be considered in future studies to examine the applicability and robustness of the results in this research. Several versions of the PROSPECT model and various empirical models have been developed so far [12,16,39,51,52] . This article only took two representative models as examples, but it is a potential idea in future crop science research to carry out stability analyses of more prediction models.
The estimation of biochemical components at leaf scale was taken into account in this article, and larger scales (canopy scale or ecosystem scale) estimation applications should also be considered to meet the requirements of smart agriculture.

Conclusions
In this study, stability of the PROSPECT model was assessed in order to estimate leaf C ab of oilseed rape. The results demonstrated that the prediction accuracy of the PROSPECT and RF models for leaf C ab retrieval was not highly influenced by reducing the spectral resolution. However, when random noise was added into the raw spectrum, the prediction accuracy of both models decreased significantly, especially for the RF model. This finding showed that the PROSPECT model had slightly higher stability against noise compared to the RF model. For varied growth stages and nitrogen treatment levels, the PROSPECT model also showed more robust prediction capability. All these findings showed that the PROSPECT model can be of potentially significant value when dealing with low spectral resolution, varied growth stages, different nitrogen treatment levels, and a small number of samples. It can be concluded that the PROSPECT model could be adopted universally in the prediction of the biochemical content of plants due to its great advantages in dealing with data collected by sensors with varied internal specification, and they also provide reference for further improvements and the promotion of biochemical component retrieval models.

Background work
All authors made significant contributions to this manuscript. Haiyan Cen, Li Zhai and Liang Wan designed the experiment; Li Zhai, Liang Wan and Xiaoran Li performed field data collection; Li Zhai, Liang Wan, Dawei Sun, Alwaseela Abdalla and Yueming Zhu contributed to the data analysis; Li Zhai, Haiyan Cen, Dawei Sun and Alwaseela Abdalla wrote the manuscript; Haiyan Cen and Yong He provided suggestions on the experiment design and discussion sections.