*Article* **Phytoremediation of Secondary Salinity in Greenhouse Soil with** *Astragalus sinicus***,** *Spinacea oleracea* **and** *Lolium perenne*

**Shumei Cai, Sixin Xu, Deshan Zhang, Zishi Fu, Hanlin Zhang \* and Haitao Zhu**

Institute of Eco-Environment and Plant Protection, Shanghai Academy of Agricultural Sciences, Shanghai 201403, China; caishumei@saas.sh.cn (S.C.); xsxofsaas@163.com (S.X.); zds234@163.com (D.Z.); fzs@foxmail.com (Z.F.); htzhu123@163.com (H.Z.)

**\*** Correspondence: zhanghanlinchick@163.com

**Abstract:** Phytoremediation is an effective and ecological method used to control soil secondary salinization in greenhouses. However, the plant–soil interactions for phytoremediation have not been studied sufficiently. In this study, three crop species (*Astragalus sinicus* (CM), *Spinacea oleracea* (SP) and *Lolium perenne* (RY)) were compared in a greenhouse experiment. The results showed that all three crops increased the soil microbial biomass, the abundance of beneficial microorganisms, available phosphorus and soil pH, and reduced the soil salt content. The crop nutrient accumulation was positively correlated with the relative abundance of bacterial 16S rRNA sequences in the soil. CM and RY respectively increased the relative abundances of *norank\_f\_Gemmatimonadaceae* and *norank\_f\_Anaerolineaceae* within the soil bacterial community, while SP increased the relative abundances of *Gibellulopsis* within the fungal community. Correlation analysis revealed that pH and total dissolved salts were the vital factors affecting soil microbial communities in the secondary salinized soil. Our results suggest that phytoremediation could effectively alleviate secondary salinization by regulating the balance of soil microbial community composition and promoting crop nutrient accumulation.

**Keywords:** phytoremediation; secondary salinization; salt tolerance; microbial diversity; nutrient accumulation

#### **1. Introduction**

Secondary soil salinization, which is mainly caused by intensified anthropogenic agricultural production, has been recognized as an extensive form of land degradation [1–3]. High inputs of agrochemicals, high evaporation, and mineral leaching in the intensified production system significantly intensify secondary soil salinization, as well as high Na+ accumulation in surface soil, which restricts agricultural productivity worldwide [4,5]. The salinization causes soil compaction and an imbalance in nutrient supply, which directly harms the normal growth of crops. Furthermore, the salinization alters the status of soil microorganisms, thereby indirectly affecting the entire ecological environment, and thus hindering the sustainable development of agricultural production [6].

Phytoremediation can alleviate secondary salinization in facility cultivation soils and reduce the dependence on mineral fertilizers [7]. In previous studies, it has been suggested that soil microbes are susceptible to farming practices, and that selecting an effective crop has positive effects on microbial communities and functions [8]. For example, applying a green manure crop has been shown to significantly change the soil microbial community composition and function [9]. Fungi, bacteria, and actinomycetes have been found to be active in the rhizosphere of Italian ryegrass [10]. The symbioses of these microorganisms accelerated nutrient cycling processes [11]. Therefore, understanding the structure of the soil microbial community and its responses to applications of different types of crops is necessary to elucidate the effects of microbes on secondary salinization.

**Citation:** Cai, S.; Xu, S.; Zhang, D.; Fu, Z.; Zhang, H.; Zhu, H. Phytoremediation of Secondary Salinity in Greenhouse Soil with *Astragalus sinicus*, *Spinacea oleracea* and *Lolium perenne*. *Agriculture* **2022**, *12*, 212. https://doi.org/10.3390/ agriculture12020212

Academic Editors: Chengfang Li and Lijin Guo

Received: 21 December 2021 Accepted: 28 January 2022 Published: 1 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In this study, it was hypothesized that the saline soil biochemical properties and microbial communities would change consistently with the type of crop species planted. Further, it was postulated that there would be considerable linkage between crop nutrient accumulation, soil salinization indicators and soil biochemical properties during the process of phytoremediation. The objectives of this study were to clarify the impact of different crop species on soil biochemical properties and microbial communities in cultivation facility soil, and to explore the linkage between crop nutrient accumulation and the composition of soil microbial communities.

#### **2. Materials and Methods**

#### *2.1. Field Site Setup, Management and Sampling*

A greenhouse experiment was performed at the Zhuanghang Comprehensive Experimental Station of the Shanghai Academy of Agricultural Sciences, Fengxian, Shanghai, China (30◦53 20.0 N, 121◦23 06.4 E). The study site was flat and the soil type was calcareous alluvium. Three replicates of four treatments were arranged in a randomized block design using 30 m × 2 m plots constructed in January 2015 (Figure 1). Nylon screen fabric was erected around every plot to avoid runoff effects, and it was buried beneath the soil surface with a height of 40 cm. Four treatments were set up, including the fallow control (CK), Chinese milk vetch (*Astragalus sinicus* L., CM), Spinach (*Spinacia oleracea* L., SP) and Ryegrass (*Lolium perenne* L., RY). CM, SP and RY were selected because they are the major native winter cover crop species that are easily accessible and widely applied to ameliorate soil salinization. Seeds of CM, SP and RY were obtained from Shanghai Nongle Planting Co. Ltd. (Shanghai, China). After 3 years of continuous planting, soil samples were collected on 30 January, 2018. In each separated plot, soil samples from the 0−20 cm surface layer were collected from 8 points to form a mixed composite soil sample, which was then divided into two parts, with one part air dried prior to the determination of basic physicochemical properties, and the other stored at −80 ◦C prior to the DNA extraction.

**Figure 1.** The schematic for field site setup, management and sampling. CK, control; CM, *Astragalus sinicus*; SP, *Spinacia oleracea*; RY, *Lolium perenne*.

#### *2.2. Crop Yield and Nutrient Accumulation*

The experimental crops were planted by sowing the equivalent of 75 kg ha−<sup>1</sup> of seed per plot in early October, and the crops were grown and harvested until the end of January annually. Equal amounts of irrigation water were supplied to each plot and no fertilization was used during the period of phytoremediation. The former crop for all treatments was pakchoi (*Brassica chinensis* L.), and urea (N 46%); compound fertilizer (17-17-17) and potassium sulphate (K2O 52%) were applied as the N, P and K fertilizer, with application rates of N 120 kg ha<sup>−</sup>1, P2O5 45 kg ha−<sup>1</sup> and K2O 90 kg ha<sup>−</sup>1. Crops were harvested at the same time as the soil samples were taken. The selected uniformly growing plants were taken to the laboratory immediately and were dried to a constant weight in preparation for nutrient determination. The total nitrogen (N), phosphorus (P) and potassium (K) contents of the dry matter were quantified using the Kjeldahl nitrogen determination method, the vanadium-molybdenum-yellow photometric method, and the flame photometry approach, respectively [12].

#### *2.3. Determination of Soil Physicochemical Properties*

Soil chemical properties, including the pH, total N, soil organic matter (SOM), alkalihydrolyzable nitrogen (AMN), available P (AP), available K (AK), and total dissolved salts (TDS), were tested according to the methods of Bao (2000). The soil pH was tested using a soil-to-water ratio of 1:2.5. The soil total N was determined via the Kjeldahl method and SOM was determined via the potassium dichromate oxidation method. The soil AMN content was measured using the alkaline hydrolysis diffusion method. The AP and AK were measured using the molybdenum blue colorimetric method and the flame photometry method, respectively. The TDS in the soil were determined using the gravimetric method. The soil microbial biomass C (MBC) and N (MBN) were measured using the chloroform fumigation method [13].

#### *2.4. Soil DNA Extraction and Microbial Community Analysis*

Bacterial and fungal DNA were extracted as three replicates from each soil sample using a FastDNA Spin Kit for Soil and were stored at −80 ◦C. The bacterial V3–V4 region of the 16S rRNA gene was amplified using the primers 338F and 806R [6]. The internal transcribed spacer (ITS) region of the fungal rRNA gene was amplified using ITS1F and ITS2R [14]. All PCR reactions were performed according to the methods described by Cai et al. [15]. Pyrosequencing was carried out by Majorbio Bio-Pharm Biotechnology Co., Ltd., Shanghai, China, using the Illumina Miseq PE250 platform. After high-throughput sequencing and optimization, 566,160 bacterial 16S rRNA sequences with 235,957,809 bp were obtained from the four treatments (N = 12), and the average sequence length was 416.8 bp. Meanwhile, 727,233 fungal ITS sequences with 173,088,882 bp were obtained, and the average sequence length was 238.0 bp. Sequence data were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under the accession number SRP273207.

#### *2.5. Real-Time Quantitative Polymerase Chain Reaction (qPCR)*

qPCR was used to examine the effects of bioremediation on soil microbial abundance. Standard reactions were performed for all samples in triplicate with an ABI7500 Real-time PCR System (Applied Biosystems INC, Foster City, CA, USA) using the SYBR green qPCR method. The standard curves and amplification curves are shown in Figures S1–S4. The qPCR mixture (20 μL) contained 10 μL of Maxima SYBR green/ROX qPCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA), 0.8 μL of each primer, 1.0 μL of template DNA and 7.4 μL of dd H2O. The amplification conditions of 16S rRNA comprised predenaturation for 5 min at 95 ◦C, followed by 30 amplification cycles of denaturation at 95 ◦C for 30 s, annealing for 30 s at 55 ◦C and extension at 72 ◦C for 1 min. The amplification conditions of ITS were nearly the same despite the difference of annealing temperature at 62 ◦C. The gene abundances of each reaction were calculated based on the constructed standard curves and then converted to copies per gram of soil, assuming 100% DNA extraction efficiency.

#### *2.6. Data Analysis*

The effects of the different crop treatments on the physicochemical soil properties and crop biomass values were tested using one-way ANOVA in SPSS 17.0 (SPSS Inc., Chicago, IL, USA). QIIME (1.7.0) software was used to calculate the alpha and beta diversities of the soil bacterial and fungal communities. The OTUs were used to characterize the alpha diversity. The Chao1, ACE, Shannon and Simpson indices were calculated. Principal coordinates analysis (PCoA) of the unweighted UniFrac distances between the samples was used to characterize the similarities (beta diversity) in the bacterial and fungal communities among the treatments [16]. The vegan data package in R was used for redundancy analysis (RDA), which was used to identify factors that affected microbial community structure.

#### **3. Results**

#### *3.1. Response of Soil Biochemical Properties to Green Manure Crops*

The application of green manure in the form of the three different crops improved the biochemical properties associated with soil fertility (Table 1). Compared to CK, all cultivation treatments displayed lower TDS contents, but higher soil pH, MBC and MBN contents (*p* < 0.05). The AP and AK contents were also significantly higher than in the control, irrespective of the type of crop applied (*p* < 0.05). The AMN and MBN contents were significantly higher in the CM treatment than in the other treatments (*p* < 0.05). The AP and AK contents were significantly higher in RY than in the other treatments (*p* < 0.05).


**Table 1.** Biochemical properties of the soil in each bioremediation treatment.

Note: TDS, total dissolved salts; SOC, soil organic carbon; AMN, alkali-hydrolyzable nitrogen; AP, available phosphorus; AK, available potassium; MBC, microbial biomass carbon; MBN, microbial biomass nitrogen. CK, control; CM, *Astragalus sinicus*; SP, *Spinacia oleracea*; RY, *Lolium perenne*. Different letters in the same column indicate a significant difference between treatments at the 0.05 level (*n* = 3).

#### *3.2. The Yield and Nutrient Uptake and Accumulation of the Different Crop Species*

Crop yield significantly differed both in terms of shoot biomass and root biomass (Figure 2). RY had a significantly higher yield than SP and CM (*p* < 0.01), with the whole fresh biomass of RY reaching 83.8 kg ha<sup>−</sup>1. As shown in Figure 3, the three crops displayed the highest cumulative uptake for K, followed by N, and then P. The cumulative K uptake by RY was significantly greater than that exhibited by SP and CM (*p* < 0.01). It was also observed that the root-to-shoot ratios of dry weight for SP and CM were significantly higher than the ratio of RY (Table S1). Further, significant differences were observed in the nutrient contents of the three crops (Figure S5). The N content of RP was significantly lower than that of CM and SP. The P content was higher in SP than in CM and RY. The average K content of the three crops was in the following order: SP > RY > CM.

**Figure 2.** Shoot and root biomass (fresh weight [FW]) of *Lolium perenne* (RY), *Spinacia oleracea* (SP) and *Astragalus sinicus* (CM). Vertical bars represent the standard error of the mean. \*\* denotes statistically significant differences between crop varieties (*p* < 0.01, Duncan's test).

**Figure 3.** Nutrient accumulation of the three different crop species, *Astragalus sinicus* (CM), *Spinacia oleracea* (SP) and *Lolium perenne* (RY). Vertical bars represent the standard error of the mean.

#### *3.3. Abundance and Diversity of Doil Microbial Communities*

Compared to the CK samples, significantly more 16S rRNA sequences but fewer ITS copies (*p* < 0.05) were found in the CM and SP samples (Table 2). However, RY had no significant effects on the number of ITS sequences in the soil compared with CK. Quantitative PCR (qPCR) data also indicated that the bacteria-to-fungi (B/F) ratio declined in the following order: SP > CM > RY > CK.


**Table 2.** Bacterial (16S rRNA) and fungal (ITS) gene copy numbers in soil samples.

Note: CK, control; CM, *Astragalus sinicus*; SP, *Spinacia oleracea*; RY, *Lolium perenne*. Different letters in the same column indicate a significant difference between treatments at the 0.05 level (*n* = 3).

The α-diversity analysis showed that the bacterial and fungal community richness (Chao1 and ACE) and diversity (Shannon and Simpson) indices varied markedly among the treatments (Table S2). Crop application increased the bacterial and fungal richness indices (*p* < 0.05) and the bacterial Shannon indices (*p* < 0.05) when compared to the control.

The crop treatments were related to an increase in the relative abundance of Proteobacteria and Bacteroidetes, and a decrease in the relative abundance of Actinobacteria for soil bacteria at the phylum level (Figure 4). The fungal community at the phylum level was comparable among all soil samples except for SP samples, which had a high abundance of Basidiomycota and Unclassified\_k\_Fungi (Figure 5). With respect to the bacterial community at the genus level, CM increased the relative abundances of *norank\_f\_Gemmatimonadaceae*. With respect to the fungal community, CM and RY both increased the relative abundances of *Chaetomium* and *Humicola*. A higher relative abundance of *Gibellulopsis* was observed in the SP samples.

**Figure 4.** Relative abundance maps of the dominant bacterial taxa in the soils of the different crop treatments based on 16S rRNA sequences at the phylum (**A**) and genus (**B**) levels. CK, control; CM, *Astragalus sinicus*; SP, *Spinacia oleracea*; RY, *Lolium perenne*.

#### *3.4. Relationships between the Relative Abundance of Soil Microorganisms and Crop Nutrient Accumulation*

The relative abundances of 16S rRNA and ITS sequences in the soils were significantly correlated with crop nutrient accumulation. The accumulation of all the nutrients by the crops was positively correlated with the relative abundance of soil bacterial and fungal sequences (Figure S6).

**Figure 5.** Relative abundance maps of the dominant fungal taxa in the soils of the different crop treatments based on internal transcribed spacer sequences at the phylum (**A**) and genus (**B**) levels. CK, control; CM, *Astragalus sinicus*; SP, *Spinacia oleracea*; RY, *Lolium perenne*.

RDA based on the soil biochemical properties explained 85.13% of the variation in the first two components of the 16S rRNA community diversity (Figure 6A). The first component (RDA1) represented 58.31% of the variability and was dominated by pH and MBN. The second component (RDA2) represented 26.82% of the variability and was dominated by AK and TDS. With respect to the ITS community diversity, the first two trait axes of the RDA accounted for 50.86% and 43.03% of the total variation, respectively; AMN scored high on the first axis, and MBN, SOC, AP, TDS, pH and MBC scored high on the second axis (Figure 6B). For the bacterial 16S rRNA community, most bacterial genera were clustered and scattered in the directions of pH, TDS and AK. Meanwhile, most fungal genera were clustered and scattered in the directions of TDS, AMN, pH and AK.

**Figure 6.** Redundancy analysis (RDA) of soil microbial community diversity and soil biochemical properties using the five most dominant genera according to bacterial 16S rRNA (**A**) and fungal internal transcribed spacer (**B**) sequences. TN, total nitrogen; SOC, soil organic carbon; AMN, alkalihydrolyzable nitrogen; AP, available phosphorus; AK, available potassium; TDS, total dissolved salts; MBC, microbial biomass carbon; MBN, microbial biomass nitrogen. CK, control; CM, *Astragalus sinicus*; SP, *Spinacia oleracea*; RY, *Lolium perenne*.

#### **4. Discussion and Conclusions**

#### *4.1. Crop Biomass Accumulation and Nutrient Absorption in Secondary Salinized Soil*

Depending on the ability of a crop to adapt to the stress of secondary salinization, crop growth and nutrient absorption differs. In this study, consistent with the yields of the crops, the order of the total amount of nutrient absorption and accumulation of the three crops was as follows: RY > SP > CM (Figures 1 and 2). Previous studies have shown that when the saline conditions of the soil are aggravated, the growth of legumes is inhibited, and the amount of biomass and nutrient absorption and accumulation decreases [17]. CM, which is a legume used as green manure, can obtain the nutrients required for crop growth through biological nitrogen fixation even in soils with low fertility. However, milk vetch is sensitive to the soil pH and salt content, which restricts its growth and nitrogen-fixing ability in the face of saline-alkali adversity [9,18]. In the present study, the fresh biomass of CM was 30.8 kg ha<sup>−</sup>1, which was only 1/3 of the average yield of RY.

In this study, it was confirmed that crop yield and nutrient content together determine the amount of nutrient accumulation. For example, the yield of RY was significantly higher than that of the other crops (Figure 1), meanwhile it also performed better with respect to nutrient accumulation. Overall, compared to CK, all cultivation treatments increased soil available nutrients except for soil AMN in RY and SP (Table 1), which may be due to high N accumulation and low N fertilizer application. The results also showed that the root-to-shoot ratio was another key factor affecting crop nutrient accumulation. Previous studies on soil salinization in cultivation facilities have noted that aboveground plant parts were less sensitive to environmental changes than belowground parts [19,20]. Similar variation trends were observed for the root-to-shoot ratios of dry weight and the aboveground N contents of the three crops in the present study. This suggests that the variation in root-to-shoot ratio of dry weight could be a good indicator of aboveground plant N uptake status in cultivation facility soils subjected to secondary salinization.

#### *4.2. Plant–Soil Feedback in Soil Subjected to Secondary Salinization*

The mechanisms underlying plant–soil feedback in agrosystems are complex. Previous studies have reported that when exposed to salt stress, plants actively change their strategy for the absorption of inorganic ions, and synthesize proline and other substances to osmotically adjust the cytoplasmic microenvironment [21,22]. Through these changes, the plants can resist the damage caused by saline-alkali stress. The results of the present study show that when crops are grown on secondary salinized soil, the absorption of K<sup>+</sup> by

crops increases. In addition, the content of AP in the soil solution increased, which was possibly related to the pH value and K+ saturation in the soil solution.

Furthermore, the results of the present study confirmed that phytoremediation increased the B/F ratio in the secondary salinized soil of the facility, especially for SP (Table 2). Previous studies have reported that SP had strong salt tolerant and antifungal ability, compared with other vegetable species on the aspects of phytoremediation and food safety [23,24]. It was found in the present study that leguminous green manure (CM) increased the relative abundances of the bacterial groups *norank\_f\_Gemmatimonadaceae*, and of the fungal genera *Chaetomium* and *Humicola*. *Gemmatimonadaceae* has been reported for the capacity of accumulating polyphosphate [25]. *Chaetomium* and *Humicola* have been found to be major groups of biological control agents, which not only reduce the incidence of soil-borne pathogens and plant disease, but also degrade a wide range of recalcitrant compounds [26]. These fungi possess a variety of genes that produce high-value enzymes, including chitinase and glucanase. The present study revealed that the nutrient accumulation of the crops was positively correlated with soil microbial communities, and soil pH, MBN, AK and TDS play important roles in maintaining microbial flora balance. However, we recommend future studies using dependency analysis of accuracy methods to create a holistic view of soil microbial succession and crop nutrient accumulation in the cultivation facility soils subjected to secondary salinization.

The utilization of phytoremediation in the form of planting salt-tolerant crops can alleviate the secondary salinization of soils in cultivation facilities. Such bioremediation can optimize the structure of the soil microbial community by increasing the soil microbial biomass, AP, AK and soil pH, and by reducing the soil salt content. The bacterial and fungal community compositions differed among the soils planted with the different salt-tolerant crop species. This study stresses the importance of phytoremediation for soils subjected to secondary salinization, and confirms that the crop species influences the correlations between crop nutrient accumulation and soil microbial community compositions.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/agriculture12020212/s1, Figure S1: qPCR standard curve for the 16S rRNA gene; Figure S2: qPCR amplification curve for the 16S rRNA gene; Figure S3: qPCR standard curve for the ITS gene; Figure S4: qPCR amplification curve for the ITS gene; Figure S5: Shoot and root nutrient contents of the three different crop species; Figure S6: Correlations (Pearson's *P*-value) between crop nutrient accumulation and soil microbial abundance for Chinese milk vetch (A,B); spinach (C,D) and ryegrass (E,F); Table S1: Biomass accumulation and root-to-shoot ratio of the bioremediation crops; Table S2: Bacterial (16S rRNA) and fungal (ITS) α-diversity in soil samples.

**Author Contributions:** Data curation and writing-original draft preparation, S.C.; investigation, S.X. and Z.F.; formal analysis, D.Z.; visualization and writing—review and editing, H.Z. (Hanlin Zhang); supervision, H.Z. (Haitao Zhu). All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Domestic Cooperation Program of Shanghai Science and Technology Commission [Grant No. 20025800500], the Agricultural Achievement Transformation and Demonstration Application Program of Shanghai Science and Technology Commission [Grant No. 20392003400] and the Outstanding Team Program of Shanghai Academy of Agricultural Science [Grant No. Hu-Nong-Ke-Zhuo 2022 (008)].

**Institutional Review Board Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest in this manuscript.

#### **References**


## *Article* **Root System Architecture Differences of Maize Cultivars Affect Yield and Nitrogen Accumulation in Southwest China**

**Song Guo 1, Zhigang Liu 2, Zijun Zhou 1, Tingqi Lu 3, Shanghong Chen 1, Mingjiang He 1, Xiangzhong Zeng 1, Kun Chen 1, Hua Yu 1, Yuxian Shangguan 1, Yujiao Dong 1, Fanjun Chen 2, Yonghong Liu <sup>4</sup> and Yusheng Qin 1,5,\***


**Abstract:** Root system architecture (RSA) plays a critical role in the acquisition of water and mineral nutrients. In order to understand the root characteristics that contribute to enhanced crop yield and N accumulation high-yielding and N efficient cultivars under N-stressed conditions. Here, grain yield, N accumulation and RSA traits of six dominant maize cultivars (CD30, ZH311, ZHg505, CD189, QY9 and RY1210) grown in the Southwestern part of China were investigated in field experiment under three different N regimes in 2019–2020; N300 (300 kg N ha<sup>−</sup>1), N150 (150 kg N ha−1) and N0 (no N supplied). Using Root Estimator for Shovelomics Traits (REST) for the quantitative analysis of maize root image obtained in the field, RSA traits including total root length (RL), root surface area (RA), root angle opening (RO), and root maximal width (RMW) were quantified in this study. The results showed that Yield, N accumulation and RSA were significantly affected by N rates, cultivars and their interactions. Grain yield, N accumulation and root weight showed a similar trend under N300 and N150 conditions compared to N0 conditions. With the input of N fertilizer, the root length, surface area, and angle increase, but root width does not increase. Under the N300 and N150 condition, RL, RA, RO and RMW increased by 17.96%, 17.74%, 18.27%, 9.22%, and 20.39%, 18.58%, 19.92%, 16.79%, respectively, compared to N0 condition. CD30, ZH505 and RY1210 have similar RO and RMW, larger than other cultivars. However, ZH505 and RY1210 have 13.22% and 19.99% longer RL, and 11.41% and 5.17% larger RA than CD30. Additionally, the grain yield of ZH505 and RY1210 is 17.57% and 13.97% higher compared with CD30. The N accumulation of ZH505 and RY1210 also shows 4.55% and 9.60% higher than CD30. Correlation analysis shows that RL, RA, RO and RMW have a significant positive correlation with grain yield while RO and RMW have a significant positive correlation with N accumulation. Linear plus plateau model analysis revealed that when the RO reaches 99.53◦, and the RMW reaches 15.18 cm, the N accumulation reaches its maximum value under 0–300 kg N ha−<sup>1</sup> conditions. Therefore, selecting maize cultivars with efficient RSA suitable for different soil N inputs can achieve higher grain yield and N use efficiency.

**Keywords:** maize; root system architecture; nitrogen rates; cultivars; yield

#### **1. Introduction**

The root is an essential organ in plants, and it plays an important role in nutrient uptake, growth and yield formation [1,2]. Root system architecture (RSA), including root

**Citation:** Guo, S.; Liu, Z.; Zhou, Z.; Lu, T.; Chen, S.; He, M.; Zeng, X.; Chen, K.; Yu, H.; Shangguan, Y.; et al. Root System Architecture Differences of Maize Cultivars Affect Yield and Nitrogen Accumulation in Southwest China. *Agriculture* **2022**, *12*, 209. https://doi.org/10.3390/ agriculture12020209

Academic Editors: Chengfang Li and Lijin Guo

Received: 26 December 2021 Accepted: 30 January 2022 Published: 1 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

length, root numbers, root surface area, root angle and root width, is an important trait in crops for the acquisition of underground resources [3,4]. The RSA of maize is influenced by genotype and environmental factors such as water, nutrients and temperature [5,6]. Nitrogen (N) is the key limiting nutrient in crop production. At the same time, Ninputinagricultural systems are also an important factor affecting environmental degradation and climate change [7]. Enhancing N use efficiency (NUE) is one of the most effective ways in sustainable agriculture to meet the 2050 global food demand projected [7,8]. Understanding the relationship between N uptake and utilization efficiency and RSA in maize is an important step towards improving maize productivity. Breeding new varieties based on RSA differences will improve N use efficiency (NUE) in maize production [3,9].

Modern varieties with deeper root distribution can increase yield under low N conditions [10]. A strong root system is an important factor for high yields [11]. There is a positive correlation between root weight and above-ground biomass and ultimately yield [12]. Under high planting density, a medium root system with more root distribution is more likely to result in high yield [6]. The interaction between RSA and soil N absorption determines yield largely [6,13]. Lynch considered that steep, cheap and deep are ideal RSA for obtaining N fertilizer and water in a low N input system in maize [5]. The maize root length and surface of different eras showed an increasing trend followed by a decreasing trend in China [14]. It has been reported that the effect of root horizontal distribution on grain yield is greater than that of root horizontal distribution [15]. Under low N conditions, the root horizontal expansion decreased [16]. The RSA of modern maize varieties in China is characterized by "horizontal contraction and vertical extension", which is more suitable for planting at a higher N level [17]. Response to N fertilizer varies with genotype, the yield and root biomass of maize varieties with high N efficiency are higher than those with low N efficiency [10,18].

High temperature, drought, poor soil fertility and nutrient leaching are persistent agronomic challenges in spring maize production in the central Sichuan Basin [19]. Thus, it is crucial to identify maize varieties with an ideal RSA suitable for this environment. Root Estimator for Shovelomics Traits (REST) is a simple, rapid and effective method for the quantitative analysis of plant root images obtained in the field [20–22]. In this study, we used this high-throughput root phenotype analysis method to study the genotypic differences and N response of RSA among maize varieties grown in this area.

#### **2. Material and Methods**

#### *2.1. Plant Materials*

A field experiment was conducted at Sichuan Agricultural Research Institute Modern Agriculture Experimental Station, Deyang City, Sichuan Province (30◦36.784 N, 105◦01.322 E) in 2019–2020. A total of 6 maize varieties were tested (Table 1), all of which are currently dominant high-yield spring maize cultivars in the hilly region of central Sichuan, namely Chengdan 30 (CD30), Zhenghong505 (ZH505), Zhenghong 311 (ZH311), Chuandan 189 (CD189), Quanyu 9 (QY9) and Rongyu 1210 (RY1210).

**Table 1.** Characteristics of cultivars used in this study.


#### *2.2. Experimental Design*

The experimental design was a split-plot with three replicates, with N fertilizer treatments in the main plots and the cultivars in the subplots. The variable between plots was three application rates of N fertilizer, namely 300 kg N ha−<sup>1</sup> (N300), 150 kg N ha−<sup>1</sup> (N150), and no N supplied (N0). The plots were fertilized with 90 kg ha−<sup>1</sup> P2O5 and 150 kg ha−<sup>1</sup> K2O. Phosphorus and potassium fertilizers were applied before sowing, 50% of the N fertilizer was applied as a base dressing, and the remaining 50% was applied at the stem jointing stage in N300 and N150 treatments. The subplots area was 20 m<sup>2</sup> (5-m-long × 4-m-width). Maize was over-seeded on 1 April in 2019 and 7 April in 2020. At the V3 stage, seedlings were thinned to a final density of 50,000 plants ha−1, which is the optimum population for maize cultivars grown in local areas. Cultivars were harvested on 11 August in 2019 and 15 August in 2020. The field was irrigated before sowing. Plots were kept free of weeds, insects and diseases with chemicals based on standard practices. No irrigation was applied during the growing season. Rainfall throughout the growing season was as shown in Figure 1.

**Figure 1.** Temperature and precipitation of Zhongjiang County during the study period in 2019 and 2020. Note: Prec, Total monthly precipitation (mm); Temp, Average monthly temperature (◦C).

Soil physical and chemical characteristics were evaluated at the beginning of the experiment for each treatment by analyzing three soil samples. The topsoil layer (0–30 cm) contained organic matter 11.5 g kg<sup>−</sup>1, total N 0.99 g kg−1, alkali-hydrolyzable N 74.0 mg kg<sup>−</sup>1, available phosphorus (Olsen-P) 14.5 mg kg−1, ammonium acetate extractable potassium (K) 172 mg kg−<sup>1</sup> and pH 7.94 (1:1.25 *g*/*v*) in N300 treatment. The chemical characteristics in N150 treatment were as follows: organic matter 10.8 g kg−1, total N 0.81 g kg−1, alkali-hydrolyzable N 32.5 mg kg<sup>−</sup>1, Olsen-P 17.5 mg kg−1, ammonium acetate extractable potassium (K) 161.8 mg kg−<sup>1</sup> and pH 7.98. In No treatment, chemical characteristics were organic matter 9.6 g kg<sup>−</sup>1, total N 0.69 g kg−1, alkali-hydrolyzable N 8.7 mg kg−1, Olsen-P 11.9 mg kg−1, ammonium acetate extractable potassium (K) 143.7 mg kg−<sup>1</sup> and pH 7.85. The soil type is brown and is classified as Cambisols with sandy loam according to the FAO classification system (IUSS Working Group WRB, 2015).

#### *2.3. Agronomic Trait Measurements*

At silking and physiological maturity stage, three uniform plants from each plot were cut at the soil surface and separated into leaves, stem and grain (only at maturity). At the silking stage, roots were excavated within a soil volume of 30 cm (length) × 30 cm (width) × 25 cm (depth) for each plant and were then shaken off a large fraction of the soil adhering to the root system. Afterward, the roots were washed under low pressure using a water hose and nozzle. Root imaging and processing were as described by Colombi et al. [21]. Briefly, root images were captured in the imaging tent using a digital camera (Canon EOS 400D, Canon, Tokyo, Japan) with a 35 mm fixed focal lens (Canon EF 35 mm f/2.0, Canon, Tokyo, Japan). The image size was 35 × 52.5 cm resulting in a pixel size of 0.13 mm. Root images analyses were processed using REST software (Figure 2). RSA traits, including total root length, surface area, angle opening, and maximal width, were quantified in this study.

**Figure 2.** Root image processing with Root Estimator for Shovelomics Traits (REST).

At both sampling dates, all samples were heat-treated at 105 ◦C for 30 min, dried at 65 ◦C to constant weight. After obtaining dry matter weight, the samples were ground into fine powder for N measurement. N concentration at silking (leaves and stem) and maturity (leaves, stem, and grain) were determined by the semi-micro Kjeldahl procedure [23]. At maturity, the ears in each plot were harvested. Grain was oven-dry to determine grain moisture content. The grain yield was determined and then standardized to 14% moisture.

#### *2.4. Statistical Analysis*

Statistical analysis was performed using SPSS Statistics 17.0 (SPSS, Inc., Chicago, IL, USA). Three-way analysis of variance (ANOVA) was used to test for significant differences among N treatment, cultivar, year and N treatment × cultivar × year interaction. N level and cultivar were treated as a fixed effect. The least significant difference test (LSD) was used to evaluate significant differences between means when a significant effect was detected by ANOVA. Means for each cultivar were used for correlation analysis and linear platform model fitting.

#### **3. Results**

#### *3.1. Grain Yield and N Accumulation Properties*

Across the two years, significant differences in grain yield and N accumulation were found among N treatments (Tables 2 and 3). ANOVA showed significant effects of N levels (0, 150 and 300 kg ha−1) (N), cultivar (C), and years (Y) on grain yield, ED, KPR and N accumulation at silking and maturity. The interaction of N × C had a significant effect on N accumulation at the silking stage, and grain yield and N accumulation at maturity. Due to the less rain from April to June 2020 (Figure 1), the grain yield and N accumulation were significantly lower than those in 2019.

The yield in N300 was similar to that of N150, with both of them having higher yields than N0, except that there was no difference in 2019. N accumulation showed a similar trend under N300 and N150 conditions compared to N0 conditions. The grain yield and N accumulation among cultivars were significantly different. The grain yield per plant ranged from 126.16 to 148.33 g, and the maximum value was 17.58% in ZH505, higher than the cultivar with the minimum value, CD30 (Table 2). There was no significant difference in grain yield betweenZH505, ZH311, CD189 and RY1210, although they were significantly higher than CD30. The grain weight of RY1210 between 2019 and 2020 had no difference at the N300 level. At the silking stage, the difference of N accumulation between varieties was mainly on the leaf. ZH311 showed the highest value in N accumulation, and 10.95% higher than CD30 which has the minimum value. However, the difference in N accumulation between varieties was mainly on leaf and grain at maturity. The minimum value of total N accumulation was observed in QY9 and 10.57% lower when compared with RY1210 (Table 3).

**Table 2.** Analysis of variance in GW, HKW, EL, ED, RPE and KPR on maize of six cultivars under three N conditions.


Within N or cultivar or year, numbers followed by different letters indicate significant difference (*p* < 0.05). \* significant at *p* < 0.05, \*\* significant at *p* < 0.01, ns: not significant (*p* > 0.05). GW, Grain weight; HKW, hundredkernel weight; EL, Ear length; BHL, Bald head length; ED, Ear diameter; RPE, Rows per ear; KPR, Kernels per row.

#### *3.2. Root System Architecture Traits Evaluation*

Variance analysis of root traits showed that N treatments had a significant effect on root traits (Table 4). The interaction of the N × cultivar had a significant effect on all root traits. Therefore, the N × cultivar was further analyzed. Root dry weight, total length, surface area, and angle opening in N150 were similar to that of N300, with both of them having higher values than the N0 condition. In addition, the maximal width in N150 was higher than that of N300 and N0 treatments. Under N150 condition, the root angle, root width, root length and root surface area increased by 19.93%, 16.79%, 20.39% and18.58% compared with no fertilizer treatment. The RW among cultivars was significantly different. The maximum value of RW was 20.36 g in RY1210 higher than the cultivar with the minimum value, QY9. The maximum value of total length was 1651.29 cm in RY1210 higher than

the cultivar with the minimum value, CD30. The maximum value of the surface area was 2851.67 cm<sup>2</sup> in ZH505 higher than the cultivar with the minimum value, CD189. The root angle opening of RY1210 and CD30 were larger than others. The root maximal width of ZH 505, RY1210 and CD30 were larger than others.


**Table 3.** Evaluation of N accumulation traits of six cultivars under three nitrogen (N) conditions.

Within N or cultivar or year, numbers followed by different letters indicate significant difference (*p* < 0.05). \* significant at *p* < 0.05, \*\* significant at *p* < 0.01, ns: not significant (*p* > 0.05).



Within N or cultivar or year, numbers followed by different letters indicate significant difference (*p* < 0.05).

\* significant at *p* < 0.05, \*\* significant at *p* < 0.01, ns: not significant (*p* > 0.05).

The root dry weights under N150 and N300 were increased by 62.24% and 58.80%, respectively, compared with N0 (Figure 3). It indicated that N application might increase maize root weight, while excessive N will reduce it (Figure 4). The root dry weight was also significantly different among cultivars. Under N0 treatment, RY1210 showed higher root weights, while QY9 had the lowest values in two years. Under N150 treatment, RY1210, and CD189 had higher root weights, while QY9 had the lowest values. The root weight under N300 treatment, of which CD189 had the least root weight compared to other varieties. CD189 were very sensitive to N stress, and the root weight decreased significantly under N deficiency or excess.

**Figure 3.** Root biomass of six maize cultivars under three N treatments in 2019 and 2020. Bars indicate standard error. Different lower case letters indicate significant differences at *p* < 0.05 among the cultivars in the same N treatment. Different capital letters indicate significant differences among the N treatments.

**Figure 4.** The distribution of the root system of six maize cultivars in three N treatments.

The total root length and root surface area were regulated by N fertilizer (Figure 5). There was no significant difference in root total length and surface area between N150 and

N300, which were higher than that of N0. Compared with N0, the root length and root surface area of N150 were increased by 20.39% and 18.58%, respectively. The root length and root surface area of N300 were increased by 17.96% and 17.74%, respectively (Table 4). This shows that N can promote the growth of roots. However, excessive application of N fertilizer inhibits root elongation and root surface area increase. The root length and root surface area of most maize varieties showed a trend of increasing at first and then decreasing at the three N levels, while RY1210 and CD189 showed an increasing trend all the time, indicating that RY1210 and CD189 were not sensitive to high N. The root length of RY1210 under N0 treatment is higher than that under N300 treatment. Root surface area under N0 is similar to that of N300. It indicates that this genotype is not sensitive to low N stress and belongs to low N tolerance varieties.

The root angle and width are regulated by N as well (Figure 5); they increased at first and then decreased with the increase in N supply, while the difference between N150 and N300 was significant in root width. The root angles of N150 and N300 were increased by 19.92% and 18.27% compared to the N0 treatment, respectively. The root maximal width of N150 and N300 was increased by 16.79% and 9.22% compared with the N0 treatment, respectively. These results indicate that N might promote the growth of maize roots, and the angle and width of the root system are significantly increased. However, excessive N may inhibit the increase in root angle and width.

The root angle opening and width between genotypes were significantly different. Under three N treatments, CD30 and RY1210 had larger root angles, CD189 and QY9 showed smaller ones, while ZH505 had a larger root width at N150 and N300, respectively. However, the change in root width was not completely consistent with the root angle. Under three N treatments, CD30, ZH505 and RY1210 had larger root widths, while QY9 had a smaller root width.

#### *3.3. Relationship between RSA, Grain Yield, and N Accumulation*

Significant correlations were found between RSA, grain yield and N accumulation (Figure 6). The grain yield increased logarithmically with increasing root length and surface area, and the regression multiple R<sup>2</sup> values were 0.84 and 0.69 (*p* < 0.01). Furthermore, N accumulation increased logarithmically with increasing root length and surface area (*p* < 0.01). With the increase of root length of maize roots, the grain yield and N accumulation continued to increase, while after reaching a certain point, the yield and N accumulation did not continue to increase, but showed a downward trend (Figure 7).

Yield and N accumulation are significantly related to root angle and root width (Figure 8). With the increase of root angle and root width, the yield and N accumulation continue to increase, while their continuous increase cannot further increase yield. On the contrary, too large root width will reduce maize yield and N accumulation. The trends in yield and N accumulation at the three N levels are consistent with the linear + plateau model (*p* < 0.05). Results from the model showed that when the root angle reaches 99.53◦ and 97.39◦, the N uptake will reach the plateau value of 2.56 g plant−<sup>1</sup> and 2.11 g plant−<sup>1</sup> in 2019 and 2020, respectively; when the root width reaches 15.18 cm and 14.83 cm, the y N accumulation plateau value will be 2.34 g plant−<sup>1</sup> and 1.90 g plant−<sup>1</sup> in 2019 and 2020, respectively. Therefore, when the root angle of cultivars reaches 99.53◦, and the root width reaches 15.18 cm, higher yield and N accumulation can be obtained.

**Figure 5.** Root total length, surface area, angel opening, and maximal width of six maize cultivars under three N treatments in 2019 and 2020: (**a**) root total length in 2019, (**b**) root total length in 2020, (**c**) root surface area in 2019, (**d**) root surface area in 2020, (**e**) root angel opening in 2019, (**f**) root angel opening in 2020, (**g**) root maximal width in 2019, (**h**) root maximal width in 2020. Bars indicate standard error. Different lower case letters indicate significant differences at *p* < 0.05 among the cultivars in the same N treatment. Different capital letters indicate significant differences among the N treatments.

**Figure 6.** Correlation coefficients between grain yield, N accumulation, and root system traits of six cultivars under three nitrogen (N) conditions. \* significant at the 0.05 probability level; \*\* significant at the 0.01 probability level. GW, Grain weight; HKW, hundred-kernel weight; EL, Ear length; BHL, Bald head length; ED, Ear diameter; RPE, Rows per ear; KPR, Kernels per row; TNS, total N accumulation at silking; TNM, total N accumulation at maturity; RW, root weight; RO, root angle opening; RA, root surface area; RMW, root maximal width; RL, root length.

**Figure 7.** Relationship among total root length, surface area, grain yield, and N accumulation in six maize cultivars under three different N treatments. (**a**) Relationship between total root length and grain yield, (**b**) Relationship between root surface area and grain yield, (**c**) Relationship between total root length and N accumulation, (**d**) Relationship between root surface area and N accumulation. \* significant at *p* < 0.05, \*\* significant at *p* < 0.01.

**Figure 8.** Relationship among angle opening, maximal width, and N accumulation in six maize cultivars under three different N treatments. (**a**) Relationship between root angle opening and N accumulation, (**b**) Relationship between root maximal width and N accumulation. The white and black cicles represent the 2019 and 2020 values, respectively.

#### **4. Discussion**

#### *4.1. Influence of RSA Traits on Grain Yield and N Accumulation*

Under field conditions, 80–90% of maize roots are distributed in 0–30 cm soil layers [24]. The root system of maize continues to grow and develop during the growth period, reaching the maximum at the silking stage and then begins to senesce [25]. Therefore, the silking stage is a critical period for analyzing root traits in fields [6]. The growth and development of roots are affected by genetics and the environment [26]. In this study, ANOVA showed significant effects of N levels, cultivar and years on grain yield, N accumulation and RSA (Tables 2–4). Their interactions (N × cultivar) suggested the existence of strong environmental effects on grain yield, NUE and RSA of the different maize cultivars. The differences in root biomass, root total length, surface area, root angle, and width among maize cultivars were significant. These root differences are the key factors that cause differences in N uptake and grain yield among these cultivars (Figures 6 and 7).

The function of the root system depends on the biomass and spatial distribution of the root system. Larger root systems are often closely related to higher yield, biomass and N uptake [27,28]. In this study, ZH505 and RY1210 have a larger root system and higher yield and N accumulation (Tables 2 and 3). The root length of RY1210 is longer, and it has a higher N uptake during the silking period because the root length of maize is directly related to the uptake and utilization of nitrate-nitrogen in the soil [29]. The root surface area of ZH505 is larger, and it increases the ability to obtain N in the soil through the root surface by diffusion in the rhizosphere [30]. The root system of QY9 is small, its yield and N uptake are low because the smaller root system will restrict the aboveground access to water and nutrients and will also lead to a decline in yield [5]. The angle between the maize root system and the ground is 10–80◦, and the angle of the root system affects the depth of root penetration [5]. The root expansion width can reflect the horizontal distribution of the root system, ultimately affecting nutrient uptake, and gradually decreasing as the root system extends downward. In this study, RY1210 had a larger root angle and width, which can expand the root growth space and promote the uptake of nutrients and water, while the root angle of ZH505 is smaller, but the N uptake is higher, which may be because the deeper root can get more N in deep soil [6]. Although CD30 had a larger root angle and width, smaller root length and root surface area affected N absorption, while modern varieties need to have higher N acquisition capacity in deep soil [18].

Adverse weather conditions, such as droughts, elevated temperatures can reduce maize yields. In this study, the grain weight of RY1210 between 2019 and 2020 was similar at the N300 level. RY1210 has larger root width and root opening angle, which may be important for yield in arid environments (Figure 5). Drought stress inhibited the growth of shoots and reduced the number of lateral roots, but the growth of roots usually continued [3]. The deeper the root angle, the greater the root growth angle has better water and N absorption capacity [31]. Under high planting density, competition between N and water among the plants, medium root systems with more root distribution are more likely to result in high yield [32]. RY1210 and ZH505 have the same RMW, but RO is smaller, and RA is larger in ZH505, making it more suitable for high-density planting (Table 4).

#### *4.2. Correlation between RSA, Grain Yield and N Accumulation*

N fertilizer application affects maize yield and N uptake by affecting root growth and distribution in the soil [33] (Tables 2–4). The growth of roots is usually improved by increasing N application, but excessive N application can inhibit the growth of lateral roots and the elongation of roots, while N deficiency can promote the increase in root biomass [16,34]. In this study, root architecture significantly correlated with yield and N accumulation, while there was no significant correlation between root angle and root width and yield (Figure 6). It may be because, under excessive N, the yield has nothing to do with root angle and root width. Studies have shown that on the clay soil of northeast China, the total root length of maize continues to increase with the increase of N application, and the total root length reaches the maximum at 168–240 kg N ha−1, then decreases with the increase of N application [10]. In this study, excessive N application not only led to a decrease in yield but also reduced the root width (Figure 5, Table 4). Under low N conditions, the root angle and width are significantly reduced (Figure 5). This is because N deficiency forces the root system to obtain N fertilizer in deep soil, and the growth becomes steeper [35].

A robust root system is an essential feature of maize N-efficient varieties. Studies have shown that N efficient varieties have higher root length and root biomass under low and medium N conditions [3]. In this study, QY9 and CD189 had significantly reduced root weight under N deficiency, while CD30, ZH505 and RY1210 had higher root biomass (Table 3). We found that applying an appropriate amount of N fertilizer (150 kg N ha−1) can significantly increase root angle by 19.93%, root width by 16.79%, root length by 20.39%, and root surface area by 18.58% compared with no fertilizer treatment. Therefore, the rational application of N fertilizer can promote the development of the maize root system and improve the absorption of nutrients in the soil, resulting in an efficient N absorption and high crop yield. The study found that biochar application can increase root angle by about 14%, root width about 20%, and root surface area about 54%, thereby increasing yield by 45%. The linear + plateau model revealed that when the RO reaches 97.39◦, and the RMW reaches 15.18 cm, the N uptake will reach a plateau (Figure 8). However, an excessively large angle will affect the depth of root penetration, which is not efficient in the uptake and utilization of nutrients and water from deep soil and reduces the plant's ability to resist adverse environmental conditions [3]. Under high-density conditions, horizontally distributed roots reduce competition among plant roots and improve yield than vertically distributed roots [15]. Therefore, it is more suitable to plant maize varieties with medium root size at high density, because too large a root system will lead to competition between root systems, resulting in decreased yield [6].

#### **5. Conclusions**

Root architecture can be an important index to evaluate high-yielding and N efficient cultivars. In this study, the RSA of RL and RS displayed a significant positive correlation with grain yield and N accumulation. The RSA of RO and RMW also showed a significant positive correlation with N accumulation, no grain yield. Oversize RO and RMW cannot further improve the N accumulation. Therefore, although CD30 has larger RO and RMW, smaller RL and RS result in lower yield and N accumulation. The root system of the cultivar ZH505 and RY1210 are moderate in size, reducing the excessive carbon consumption and penetrating down to absorb N in deeper soil. The root architecture of RY1210 has good

adaptability to the arid climate. In comparison, ZH505 is suitable for high-density planting. It can be expected, selecting maize cultivars with an ideal root system architecture and applying the appropriate N fertilizer inputs in the hilly area of Southwest China can improve N efficiency and crop yield in a sustainable way.

**Author Contributions:** S.G., X.Z. and Y.Q. conceived of and designed the study; S.G., Z.L., M.H. and Z.Z. analyzed the data; S.G. wrote the manuscript; S.G., K.C., H.Y., Y.S. and Y.D. carried out the field measurements and soil analysis; F.C. and Y.L. assisted with manuscript writing and editing; T.L. and S.C. provided the material. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (41907080), the National Key R&D Program of China (2018YFD0200700), the Sichuan Science and Technology Program (2020YFH0171), and the Research Foundation of Sichuan Academy of Agricultural Sciences (2019QYXK022, 2020BJRC005).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data can be achieved by contacting with the first author and respondense author.

**Acknowledgments:** We sincerely thank the reviewers for their helpful comments and supplementary proposal. We thank Yuanzhong Bianji for editing the English text of a draft of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

