Higher Radiation Use E ﬃ ciency Produces Greater Biomass Before Heading and Grain Yield in Super Hybrid Rice

: To reveal the physiological mechanism underlying the yield advantage of super hybrid rice compared with inbred super rice, a super hybrid rice cultivar Yliangyou 3218 (YLY) and an inbred super rice cultivar Zhendao 11 (ZD) were ﬁeld grown under ﬁve nitrogen (N) fertilizer rates in 2016 and 2017. The average grain yield of YLY across nitrogen fertilizer rates was 10.1 t ha − 1 in 2016 and 9.7 t ha − 1 in 2017, 29.6% and 21.3% higher than that of ZD in 2016 and 2017, respectively. YLY showed higher above-ground biomass accumulation, especially growth before heading, which was mainly due to its faster green leaf area index (GLAI) formation and greater maximum GLAI (GLAI max ). The daily radiation interception (RI daily ) was 15.0% higher in YLY than ZD, but the accumulated radiation interception (RI acc ) before heading showed little di ﬀ erence between them because ZD had a longer growth duration. The radiation use e ﬃ ciency (RUE) of YLY before heading was 54.7% higher than that of ZD (YLY, 2.12 g MJ − 1 ; ZD, 1.37 g MJ − 1 ). Our result demonstrated that the yield advantage of YLY was due to its higher above-ground biomass before heading, which was mainly achieved by its improvement in RUE rather than radiation interception.


Introduction
Rice is the main staple food in the world. To meet the substantial increases in food demand due to the rapidly increasing population, rice yield needs to improve further. China, one of the largest rice producers, has successfully increased rice yield through genetic improvement and improved crop management practices [1][2][3]. However, compared with a 3.7% yearly yield increase in 1980s, the increment dropped to 0.9% in 1990s [4]. In order to break the yield ceiling, a "super rice" mega-project was established in 1996 in China [5]. Today, 131 "super" rice varieties have been released. Moreover, a yield more than 12 t ha −1 was found in the field experiments which adopted these varieties [6].
Many studies have been carried out to explore the mechanism underlying the high yield potential of "super" rice [7][8][9]. Compared with ordinary rice varieties, "super" rice varieties have a larger panicle size [10,11], larger green leaf area, longer leaf area duration [4], higher leaf photosynthesis rate [12], greater above-ground biomass [9] and improved root physiological traits [7]. However, it is argued that enhanced properties at the organ level would lead to increased crop yield. For example, Alex et al. [13] found that the enhancement of single leaf photosynthesis had no impact on canopy

Experimental Setup and Design
Field experiments were carried out in 2016 and 2017 in Rugao County (32 • 26 24"N, 120 • 29 24"E), Jiangsu Province, China. The chemical properties of the upper soil (0-20 cm) were as follows: pH 7.54, organic matter 20.2 g kg −1 , total N 1.97 g kg −1 , Olsen-P 10.8 mg kg −1 , NH 4 OAc-K 92.0 mg kg −1 . A split-plot design with three replicates was used in the experiments: N rates were the main plots and cultivars were the subplots. The five N treatments were as follows: N0 (0 kg N ha −1 ), N90 (90 kg N ha −1 ), N180 (180 kg N ha −1 ), N270 (270 kg N ha −1 ), N360 (360 kg N ha −1 ). N fertilizer supplied as urea was used on four application dates: 40% as basal fertilizer, 20% at mid-tillering stage (MT), 20% at panicle initiation stage (PI), 20% at spikelet differentiation stage. Phosphorus (75 kg P 2 O 5 ha −1 ) and potassium (90 kg N ha −1 ) were supplied equally for all treatments. Phosphorus fertilizer supplied as calcium superphosphate was used as basal fertilizer and potassium fertilizer supplied as potassium chloride Agronomy 2020, 10, 209 3 of 13 was used at the basal and panicle initiation stage at a ratio of 1:1. Yliangyou 3218 (YLY, super hybrid rice) and Zhendao 11 (ZD, inbred super rice) were used in this study. YLY (Y58S × Xianghui 3218) is an indica hybrid cultivar developed in 2005. ZD (ZD88 × Wuyunjing 8) is a conventional japonica cultivar developed in 2004. The individual plot size was 25 m 2 . Pre-germinated seeds were sown in the seedbed one month prior to transplanting. Uniform seedlings were selected for transplanting at a hill spacing of 20 × 14 cm with two seedlings per hill on 22 June 2016 and 20 June 2017. Weed, pest and disease were controlled according to local agricultural practices for rice production.

Meteorological Conditions
A meteorological station located next to the field was used to record air temperature (°C), rainfall (mm), and photosynthetic active radiation (MJ m −2 ). Thermal times (TT, • Cd) were estimated by daily mean temperature n days after transplanting (T an ): where T b was the base temperature when crop growth could be negligible. For rice, T b was 10 • Cd according to previous studies [22,23].

Green Leaf Area Index and Above-Ground Biomass
Green leaf area index was determined by a plant canopy analyzer (PCA, LAI-2200C, LI-COR, Lincoln, NE, USA) at MT, PI and heading stages (HD). Meanwhile, six hills of plants were harvested from each plot, and then these plants were separated into two or three parts: leaf, stem and panicles. Samples were dried with an oven until a stable dry weight was attained.

Estimation of the Growth Pattern of Green Leaf Area
As suggested by Yin et al. [24], a beta growth function was used to fit the dynamics of green leaf area index: where t is thermal days ( • Cd), t e ( • Cd) represents the time when the growth of green leaf area ended, and t m ( • Cd) is the time at which a maximum green leaf area index was obtained. GLAI max (m 2 m −2 ) is the green leaf area index at t e . The maximal growth rate of green leaf area (C m , m 2 m −2°C d −1 ) was calculated as follows:

Fractional PAR Interception (f) and Light Extinction Coefficient (k)
When GLAI determination was completed, measurements of photosynthetic active radiation (PAR) above and below the canopy were conducted with a 1m line quantum sensor (MQ-306, APOGEE, North Logan, UT, USA) from 11:00 to 13:00 in order to estimate fractional PAR interception (f ). PAR was recorded at least 10 locations in each plot. The f was estimated as follows: where PAR 1 is PAR above the canopy and PAR 2 is PAR below the canopy. The relationship between GALI and f can be described with an asymptotic equation: where k is the light extinction coefficient, an indicator of canopy structure, which can be solved with the known value of f and GLAI. In the current study, there was no difference in k values among treatments, therefore an averaged k (e.g., 0.428) of all the treatments was used.

Light Interception and Radiation Use Efficiency (RUE)
As mentioned above, photosynthetic active radiation was recorded and daily radiation interception (PAR a ) was calculated based on it: The reflection of PAR by rice canopy was neglected. The sum of PAR a was the accumulated PAR interception (RI acc ) at given growth period. As shown in Equation (6), PAR a calculation requires daily values of GLAI, in the present study, GLAI was measured at given periods as previous studies did. Thus, daily values of GLAI was estimated from transplanting to HD by the Equation (2) in 2016 and a linear decay of GLAI from HD to PM was assumed ( Figure S1). There was not enough sampling data for non-linear fitting in 2017 as 2016 did, instead, a linear behavior of GLAI from MT to HD was assumed according to Cosentino et al. [25]. Further, the slope of linear relationship between RI acc and AM was RUE.

Yield and Yield Components
At the maturity stage (PM), the grain yield was estimated from a 5 m 2 area plants in the central part of each plot and adjusted to a 14% water content. Then, 12 hill plants in each plot were collected to determine the above-ground biomass and yield components (panicles per hectare, spikelets per panicle, grain filling percentage and 1000-grain weight). The filled grains and unfilled grains were distinguished by winnowing, and the grain filling percentage was the ratio of filled grains to the sum of spikelets.

Statistical Analysis
Analysis of variance (ANOVA) was used to analyze the data by SPSS 18.0 (SPSS, Chicago, IL, USA), and the significant difference between cultivars and nitrogen treatments within a cultivar was accessed by the least significant difference (LSD) test (p < 0.05). Model fitting to determine parameters associated with the dynamics of green leaf area index was conducted by OriginPro 8.0 software (OriginLab Corporation, Northampton, MA, USA).

Climatic Conditions
As shown in Figure 1

Yield and Yield components
Across two cultivars, nitrogen fertilizer demonstrated a significant enhancement of grain yield in both 2016 and 2017. Meanwhile, YLY and ZD achieved maximum yield at different N rates, i.e., for ZD it was achieved at about 180 kg ha −1 N, but the grain yield of YLY did not increase significantly when N rate was above 90 kg ha −1 ( Figure S2). When compared the yield between these two cultivars, it was found that YLY had a yield advantage than ZD. In 2016, the grain yield of YLY was 16.0% to 46.7% higher than that of ZD at different N rates. Similarly, the yield increment was 4.3% to 41.1% in 2017 (Table 1).
When considering the yield components, the higher grain yield of YLY compared with ZD was mainly due to significant increase in panicles ha −1 and spikelets panicle −1 . As for the grain filling percentage and 1000-grain weight, they were higher in ZD than YLY in 2016. However, no significant difference was observed in 2017(( Table 1).

Yield and Yield components
Across two cultivars, nitrogen fertilizer demonstrated a significant enhancement of grain yield in both 2016 and 2017. Meanwhile, YLY and ZD achieved maximum yield at different N rates, i.e., for ZD it was achieved at about 180 kg ha −1 N, but the grain yield of YLY did not increase significantly when N rate was above 90 kg ha −1 ( Figure S2). When compared the yield between these two cultivars, it was found that YLY had a yield advantage than ZD. In 2016, the grain yield of YLY was 16.0% to 46.7% higher than that of ZD at different N rates. Similarly, the yield increment was 4.3% to 41.1% in 2017 (Table 1). When considering the yield components, the higher grain yield of YLY compared with ZD was mainly due to significant increase in panicles ha −1 and spikelets panicle −1 . As for the grain filling percentage and 1000-grain weight, they were higher in ZD than YLY in 2016. However, no significant difference was observed in 2017 (Table 1).

Above-Ground Biomass (AM)
For each rice cultivar, above-ground biomass increased with the increasing application of N fertilizer either in 2016 or 2017. After averaging above-ground biomass across different growth stages and all N treatments, the annual above-ground biomass of YLY increased by 37.5% relative to that of ZD. As for the growth before heading, the above-ground biomass increment was 47.4%. Moreover, no significant difference was found between 2016 and 2017 in the above-ground biomass of each rice cultivar for all N rates ( Figure 2; Table S1).

Above-Ground Biomass (AM)
For each rice cultivar, above-ground biomass increased with the increasing application of N fertilizer either in 2016 or 2017. After averaging above-ground biomass across different growth stages and all N treatments, the annual above-ground biomass of YLY increased by 37.5% relative to that of ZD. As for the growth before heading, the above-ground biomass increment was 47.4%. Moreover, no significant difference was found between 2016 and 2017 in the above-ground biomass of each rice cultivar for all N rates (Figure 2; Table S1).

Green Leaf Area Index (GLAI) and Fractional PAR Interception (f)
As shown in Figure 3, GLAI increased significantly due to increasing N application rate in 2016 and 2017 for each cultivar. When considering the GLAI difference between ZD and YLY, it was found that GLAI of YLY averaged across different N rates and sampling dates was 35.2% higher than that

Green Leaf Area Index (GLAI) and Fractional PAR Interception (f)
As shown in Figure 3, GLAI increased significantly due to increasing N application rate in 2016 and 2017 for each cultivar. When considering the GLAI difference between ZD and YLY, it was found that GLAI of YLY averaged across different N rates and sampling dates was 35.2% higher than that of ZD in 2016, and the GLAI enhancement was 26.2% in 2017 (Table S2). Then GLAI dynamics in 2016 was fitted by the beta function, as shown in Table 2 and at least 91% variation could be explained by this equation. The maximum GLAI (GLAI max ) decreased due to limited N supply. The GLAI max averaged across nitrogen rates in YLY was 10.7% higher than that of ZD. The time when maximum of green leaf area index is reached (t e ) was delayed because of increasing N rate in ZD, however, the opposite trend was found in YLY. The t e averaged across nitrogen rates in YLY and ZD was 1226 and 1403°Cd, respectively. Increased nitrogen supply delayed the time at which the maximum growth rate is achieved t m in both cultivars. The t m averaged across nitrogen rates in YLY and ZD was 256 and 346°Cd, respectively. Irrespective of cultivars, enhanced nitrogen supply improved the maximum growth rate of green leaf area C m . Moreover, C m in YLY was consistently higher than that in ZD (Table S3).
Agronomy 2020, 10, 209 7 of 13 opposite trend was found in YLY. The te averaged across nitrogen rates in YLY and ZD was 1226 and 1403 ℃d, respectively. Increased nitrogen supply delayed the time at which the maximum growth rate is achieved tm in both cultivars. The tm averaged across nitrogen rates in YLY and ZD was 256 and 346 ℃d, respectively. Irrespective of cultivars, enhanced nitrogen supply improved the maximum growth rate of green leaf area Cm. Moreover, Cm in YLY was consistently higher than that in ZD (Table S3).  GLAImax, the maximum value of green leaf area index (m 2 m −2 ); tm, the time at which the maximum growth rate is achieved ( d); t ℃ e, the time when maximum of green leaf area index is reached ( d); C ℃ m, the maximum growth  GLAI max , the maximum value of green leaf area index (m 2 m −2 ); t m , the time at which the maximum growth rate is achieved (°Cd); t e , the time when maximum of green leaf area index is reached (°Cd); C m , the maximum growth rate of green leaf area (m 2 m −2°C d −1 ). Data with different letters represents significant difference (p < 0.05) between N treatments within cultivar; * represents significant difference (p < 0.05) between YLY and ZD under the same N treatment.
Similar to GLAI, f was restricted due to N deprivation for each year and both cultivars (Figure 4). After averaging data across different sampling dates and N rates, the average f of YLY increased by 19.3% and 10.4% compared with those of ZD in 2016 and 2017, respectively (Table S4). rate of green leaf area (m 2 m −2 d ℃ -1 ). Data with different letters represents significant difference (p < 0.05) between N treatments within cultivar; * represents significant difference (p < 0.05) between YLY and ZD under the same N treatment.
Similar to GLAI, f was restricted due to N deprivation for each year and both cultivars ( Figure  4). After averaging data across different sampling dates and N rates, the average f of YLY increased by 19.3% and 10.4% compared with those of ZD in 2016 and 2017, respectively (Table S4).

Radiation Interception and Radiation-Use Efficiency (RUE)
For both cultivars within each year, an increasing trend was found for accumulated radiation interception (RIacc) with elevated nitrogen supply. RUE was estimated from the slope of the linear relationship between RIacc and AM, which shown that the average RUE across N rates and experiment years of YLY and ZD was 2.12 and 1.37 g MJ −1 , respectively. The RUE averaged across cultivars and N rates in 2017 was slightly higher than that in 2016 ( Figure 5).

Radiation Interception and Radiation-Use Efficiency (RUE)
For both cultivars within each year, an increasing trend was found for accumulated radiation interception (RI acc ) with elevated nitrogen supply. RUE was estimated from the slope of the linear relationship between RI acc and AM, which shown that the average RUE across N rates and experiment years of YLY and ZD was 2.12 and 1.37 g MJ −1 , respectively. The RUE averaged across cultivars and N rates in 2017 was slightly higher than that in 2016 ( Figure 5).

Leaf Net Photosynthesis Rate (A)
Leaf net photosynthesis rate (A) was enhanced by increased N supply in most cases ( Figure 6). During tillering stage, A of YLY was significantly higher than that of ZD under N0 treatment. However, there was no difference between these two cultivars under N180 treatment. During the panicle initiation stage, irrespective of N rates, A of YLY was consistently superior to that of ZD.

Leaf Net Photosynthesis Rate (A)
Leaf net photosynthesis rate (A) was enhanced by increased N supply in most cases ( Figure 6). During tillering stage, A of YLY was significantly higher than that of ZD under N0 treatment. However, there was no difference between these two cultivars under N180 treatment. During the panicle initiation stage, irrespective of N rates, A of YLY was consistently superior to that of ZD.

Leaf Net Photosynthesis Rate (A)
Leaf net photosynthesis rate (A) was enhanced by increased N supply in most cases ( Figure 6). During tillering stage, A of YLY was significantly higher than that of ZD under N0 treatment. However, there was no difference between these two cultivars under N180 treatment. During the panicle initiation stage, irrespective of N rates, A of YLY was consistently superior to that of ZD.

Discussion
RI acc or RUE: what contributed more to greater biomass accumulation and yield of YLY?
In our study, YLY achieved 4.1-51.4% higher grain yield than ZD across nitrogen rates and years (Table 1). This result confirmed superior yield potential of super hybrid rice in previous studies [26,27]. In terms of yield components (number of panicles, spikelet per panicle, weight of 1000 spikelet and spikelet sterility), it was found that the yield superiority of YLY compared with ZD mainly resulted from the spikelet per panicle, which was consistent with previous studies [26,27]. Besides sink size, source-related characteristics, such as greater biomass accumulation and higher leaf area index after heading are also important for yield potential of super rice achievement [10,28,29]. Our data also confirmed these with respect to source traits as mentioned above (Figures 3 and 4). And averaged across nitrogen rates (N0, N90 and N180) the HI of YLY were 0.564 and 0.567 in 2016 and 2017, and the corresponding values of ZD were 0.513 and 0.520 ( Figure S3). This finding was in line with Mae et al. [30], who reported that a large-grain rice cultivar, Akita 63 characterized as high sink capacity (the number of spikelets per unit land area × grain size) exhibited a high HI. It is demonstrated that temporal development of biomass gain (e.g., growth rate) fluctuates largely due to variation in canopy photosynthesis, therefore biomass increment during different stages varies [31], indicating that the contribution of rice growth during different stages to final grain yield is not equal. Our analysis found that rice yield was more dependent on biomass growth before heading rather than after heading in both years (Figure 7), and a similar result was observed by Takai et al. [32]. demonstrated that temporal development of biomass gain (e.g., growth rate) fluctuates largely due to variation in canopy photosynthesis, therefore biomass increment during different stages varies [31], indicating that the contribution of rice growth during different stages to final grain yield is not equal. Our analysis found that rice yield was more dependent on biomass growth before heading rather than after heading in both years (Figure 7), and a similar result was observed by Takai et al. [32]. The yield difference between YLY and ZD was determined by biomass growth before heading, and the biomass improvement in YLY is related to its larger N accumulation and higher nitrogen utilization efficiency (Table S5), which was consistent with Huang et al. [33]. Biomass production is dependent on RIacc and RUE. Then we further analyze RIacc and RUE of these two cultivars during rice growth stages. The results found that daily radiation interception (RIdaily) and RUE of YLY was both enhanced compared with that of ZD (Figures 5 and Figure S3). To better explain the relationship between radiation component (RIacc and RUE) and yield, the beta growth function was used to fitting the green leaf area growth, and the function explained green leaf area growth perfectly ( Table 2). According to the results of simulation, the enhanced RIdaily was mainly achieved by faster growth rate of green leaf area (e.g., smaller tm and larger Cm) given that no significant difference exists in the values of light extinction coefficient (k) between cultivars or N rates. Similar values of k demonstrate that the canopy structure of these two cultivars showed little difference, and previous studies also reported that prevailing rice cultivars are close to the optimum canopy structure due to the ideotype The yield difference between YLY and ZD was determined by biomass growth before heading, and the biomass improvement in YLY is related to its larger N accumulation and higher nitrogen utilization efficiency (Table S5), which was consistent with Huang et al. [33]. Biomass production is dependent on RI acc and RUE. Then we further analyze RI acc and RUE of these two cultivars during rice growth stages. The results found that daily radiation interception (RI daily ) and RUE of YLY was both enhanced compared with that of ZD ( Figure 5 and Figure S3). To better explain the relationship between radiation component (RI acc and RUE) and yield, the beta growth function was used to fitting the green leaf area growth, and the function explained green leaf area growth perfectly ( Table 2). According to the results of simulation, the enhanced RI daily was mainly achieved by faster growth rate of green leaf area (e.g., smaller t m and larger C m ) given that no significant difference exists in the values of light extinction coefficient (k) between cultivars or N rates. Similar values of k demonstrate that the canopy structure of these two cultivars showed little difference, and previous studies also reported that prevailing rice cultivars are close to the optimum canopy structure due to the ideotype approach [6,19,34]. Furtrher, many studies have shown that k was not influenced by nutrition status, such as nitrogen and phosphorous deficiency [15,35]. Because GLAI was observed discontinuously in the present study, simulated GLAI from the beta growth function was used to recalculate RI acc and RUE, and the simulated RUE predicted RUE well ( Figure S5). The results found that although the RI daily was improved in YLY compared with ZD, the RI acc of the former was equal to the latter before heading, because the growth duration of ZD was longer than YLY (Tables S6 and S7). Because of the delayed growth of green leaf area (e.g., larger t e , Table 2), the RI acc in ZD after heading was larger than that in YLY. The growth duration may be explained by genotype difference and environmental condition of the experiment site. As suggested by Yoshida [36], indica genotypes prefer warmer conditions compared with japonica genotypes. Thus, the low temperature during the end of rice growth season induced earlier senescence in YLY, and a similar result was found in Wei et al. [27]. As for RUE, the averaged value across N rates in ZD before heading was 1.33 g MJ −1 and in YLY was 2.12 g MJ −1 , which was close to the values reported by previous studies (Table S6). The increased RUE in YLY may be attributed to its larger leaf photosynthesis rate during panicle initiation stage ( Figure 6). Considering the complexity of canopy structure, the improved leaf photosynthesis may not necessarily result in increased canopy photosynthesis. To gain a comprehensive understanding of the mechanism of RUE, canopy photosynthesis during the whole rice growth season needs to be investigated. Observed and simulated data both showed that the enhancement of RUE was consistently larger than that of RI acc ( Figure 5 and Figure S4). Thus, it is reasonable to draw the conclusion that the biomass enhancement of YLY compared to ZD was mainly achieved by RUE improvement, indicating that yield superiority of YLY was more closely related with RUE. A cross location analysis stated that the rice yield over a wide range of regions was more dependent on RI acc than RUE although they admit that RUE improvement is a feasible way to yield promotion based on a nearly a twofold variation in RUE (e.g., 0.99-1.88 g MJ −1 ) [37]. As discussed above, there exists large variation in RUE during the whole rice growth period, and using an integrated RUE to represent the rice canopy RUE may miss some important information. Take our study for example, the RUE value of the whole season in YLY and ZD showed no difference according to our simulation (Table S6), but it is not reasonable to draw the conclusion that rice yield is more dependent on RI acc since RUE values before heading between these two cultivars indeed had significant difference. In the present study, we stated that RUE improvement was more related to rice yield potential achievement, but we have to admit that the results obtained here may be limited by the specific environmental conditions and cultivars. Cross location studies should be implemented to better understand the most crucial determinants of rice yield potential.

Conclusions
Super hybrid rice YLY showed a higher grain yield than inbred super rice ZD. The yield advantage of YLY was mainly achieved by greater above-ground biomass production before heading, larger sink size as showing more spikelets per panicle. Considering the larger increment of RUE than RI acc when comparing the difference between YLY and ZD, we draw the conclusion that the RUE contributed more to yield improvement than RI acc .
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/2/209/s1, Figure S1. Schematic representation of the calculation of simulated GLAI for rice from transplanting to maturity in 2016; Figure S2. Estimation of the minimum amount of nitrogen fertilizer for the highest yield of YLY and ZD using the model of linear plus platform.; Figure S3. The harvest index (HI) of YLY and ZD averaged across N0, N90 and N180 in 2016 and 2017; Figure S4. Normalized value of growth rate (GR), daily radiation interception (RI daily ) and radiation use efficiency (RUE) of YLY relative to ZD for the different nitrogen treatments in 2016 and 2017; Figure S5. The relationship between predicted RUE derived from simulated GLAI and observed RUE calculated with field-observed GLAI in 2016. The solid line is the 1:1 line and the dashed line is the fitted linear regression. Table S1 ANOVA analyses of above-ground biomass (AM) in 2016 and 2017.; Table S2 ANOVA analyses of green leaf area index (GLAI) in 2016 and 2017; Table S3 ANOVA analyses of parameters of green leaf area growth in 2016; Table S4 ANOVA analyses of fractional PAR interception (f ) in 2016 and 2017; Table S5 Nitrogen accumulation (N acc ), aboveground biomass (AM) and nitrogen utilization efficiency (NUtE) from transplanting to heading in 2016 and 2017.; Table S6 Accumulated radiation interception (RI acc ), daily radiation interception (RI daily ) and radiation use efficiency (RUE) estimated by green leaf area index (GLAI) derived from leaf growth function in 2016; Table S7 The growth stages for YLY and ZD at Rugao County, Jiangsu Province, China in 2016 and 2017.
Author Contributions: Y.P. undertook the field sampling, analysis and manuscript preparation; S.G. assisted with filed sampling and analysis of 2016; X.M. and S.W. assisted with sampling analysis of 2017; K.X., J.L. and Z.L. provided manuscript editorial support. S.G. conducted experiment design and data analysis. All authors have read and agreed to the published version of the manuscript.