PLANT DENSITY, IRRIGATION AND NITROGEN MANAGEMENT: THREE MAJOR PRACTICES IN CLOSING YIELD GAPS FOR AGRICULTURAL SUSTAINABILITY IN NORTH-WEST CHINA

Agriculture faces the dual challenges of food security and environmental sustainability. Here, we investigate current maize production at the field scale, analyze the yield gaps and impacting factors, and recommend measures for sustainably closing yield gaps. An experiment was conducted on a 3.9-ha maize seed production field in arid north-western China, managed with border and drip irrigation, respectively, in 2015 and 2016. The relative yield reached 70% in both years. However, drip irrigation saved 227 mm irrigation water during a drier growing season compared with traditional border irrigation, accounting for 44% of the maize evapotranspiration (ET). Yield variability under drip irrigation was 12.1%, lower than the 18.8% under border irrigation. Boundary line analysis indicates that a relative yield increase of 8 to 10% might be obtained by optimizing the yield-limiting factors. Plant density and soil available water content and available nitrogen were the three major factors involved. In conclusion, closing yield gaps with agricultural sustainability may be realized by optimizing agronomic, irrigation and fertilizer management, using water-saving irrigation methods and using site-specific management. mixture (H 2 SO 4 -HClO 4 ) and continuous flow analysis. NO 3– -N and NH 4+ -N were determined by extracting samples with 1 mol·L − 1 KCL solution and continuous flow analysis. AP was determined by extracting samples with 0.5 mol·L − 1 NaHCO 3 solution and the molybdate colorimetric method. AK was determined by extracting samples with 1 M NH 4 OAc solution and flame photometry. Soil gravimetric water content (  ') and bulk density (BD) were determined by the oven drying method. Soil texture was determined with a laser particle size analyzer (Mastersizer2000,


INTRODUCTION
Two main challenges facing agriculture are food security and environmental sustainability [1,2] . A 60% to 110% increase in global crop production (from the 2005 to 2007 levels) may be needed by the year 2050 because of the increasing population and energy needs [3,4] . Slowing down of rate of increase in cereal yield or even yield stagnation has also been a global trend in recent years [5] . Increasing crop production without expanding the land area used involves increasing productivity on the existing arable land. One way is to increase potential yields and another is to close the yield gap between those achieved by farmers and potential yields [1,6] . However, increasing potential yields is challenging because crop harvest index and radiation interception efficiency have been greatly improved since the first green revolution. Great room for improvement exists in increasing the conversion efficiency of intercepted radiation into biomass energy, but there is a two-to three-decade gap between the conceptual design at the experimental level and availability to farmers [7] . It is vital to close existing yield gaps for food security.
Crop yields at the field scale are the comprehensive result of multiple factors. For example, climatic conditions and cultivar features define the potential yield, water and nutrients limit the potential yield, and weeds and pests serve as reducing factors [8] . Irrigation and nutrient management are essential parts of crop production. Global agriculture consumes about 70% of the freshwater withdrawals [1,9] . Global water consumption is projected to increase by 21% by 2050 [9] and domestic, industrial, and livestock water demands are still growing rapidly. Thus, shrinkage of the share for irrigation may be inevitable in the face of limited water resources. Climate change, ecosystem degradation, unsustainable use of groundwater, and pollution exacerbate water scarcity [9] . Global fertilizer use also faces serious problems as in some regions insufficient nutrients are a major factor limiting crop yields. However, excessive fertilization has caused environmental problems in other regions, especially in East Asia [1,2,10] . Increasing crop yields with sustainability requires greater efforts to increase water and fertilizer use efficiency.
Numerous studies have investigated the interrelationships among crop yield, resource use efficiency and multiple variables at different scales. Agronomists have screened minimum data sets (MDS) to assess correlations between soil quality and yield at the farm or regional scales [11][12][13] . The MDS may include soil physical variables (e.g., bulk density and clay, silt and sand contents), soil chemical properties (e.g., pH in water and organic carbon, total nitrogen, ammonium nitrogen, exchangeable calcium and available phosphorus contents), and soil microbial variables (e.g., sucrase activity, earthworm numbers, relative amounts of bacterivorous nematodes, microbial biomass C and biomass N). Mueller et al. [2] examined yield gaps at the global scale and concluded that 75% of attainable yields could be realized by increasing irrigation area and nutrient application on all underachieving areas. Furthermore, minimal changes in total worldwide nitrogen and phosphate use would be possible by reducing nutrient imbalances and inefficiencies. Chen et al. [14] evaluated factors affecting summer maize yield and nitrogen use efficiency on small farms on the North China Plain and concluded that plant density was more important than other management practices in closing yield gaps. Li et al. [15,16] analyzed the spatiotemporal distribution of irrigation water productivity and its driving factors in the Hexi Corridor, north-west China. They concluded that increasing irrigation water productivity should rely on the use of water-saving technology, optimum fertilizer rates, agricultural mulching film, and pesticides. Although a large amount of research has been conducted, sustainable agricultural development in China continues to face many challenges such as water shortage, low resource use efficiency, and environmental contamination [17] . Moreover, few system analyses and interdisciplinary studies have been conducted.
Further information is needed on current crop production conditions to increase crop yield, reduce agricultural resource inputs or increase resource use efficiency, and reduce the impact of agriculture on the environment. For example, the extent of current yield gaps, the major causes of these, and how much they can be closed and by what means are all poorly understood. Boundary line analysis (BLA) has been widely used to determine critical levels of nutrients [18,19] and analyze yield gaps [20] . It simplifies a multivariate model into analysis of the relationships between the dependent variable and each independent variable to make interpretation easier [21] . The objectives of the present study were to investigate current yields and yield-limiting factors in maize production at the field scale, determine the major yield-limiting factors and quantify the yield gaps using BLA, and compare the feasibility of and measures for closing yield gaps and maintaining agricultural sustainability.

Experimental site and field management
The experiment was conducted on a 3.9-ha maize seed production field during the growing seasons of 2015 and 2016 at Shiyanghe Experimental Station of China Agricultural University (37°52′ N, 102°50′ E, 1581 m), located in Wuwei City, Gansu Province, north-west China. This region is located in an arid inland temperate climatic zone with mean annual pan evaporation of 2000 mm, mean annual precipitation of 164 mm, and an average groundwater table lower than 25 m below the ground surface [22] . Precipitation during the maize growing season was 149 and 115 mm in 2015 and 2016, respectively. Field management was operated by a farmer according to local practices. The cultivars (the female cultivar was FN1390 in 2015 and WH01 in 2016) and agronomic recommendations were provided by local seed companies. The field was partly covered by plastic film mulch, with every 1.2 m width covered strip across the whole field was separated by 0.46 m width uncovered strip. In 2015, every four rows of the female parents with one row of the male parents were sown on each mulched strip at a sowing rate of 9.34 × 10 4 female plants ha −1 with 2.34 × 10 4 male plants ha −1 . In 2016, every seven rows of the female parents with one row of the male parents were planted on two adjacent mulched strips at a sowing rate of 8.17 × 10 4 female plants ha −1 with 1.17 × 10 4 male plants ha −1 . Each mulched strip was controlled by two laterals of the drip irrigation system (the inner flat emitter drip tape was used with an emitter spacing of 30 cm and a flow rate of 3 L· h −1 ) in 2016. The average spacing within rows was 25.8 cm in both years. Female parents were sown on April 13-14, 2015 and April 23-24, 2016. The pollination period was extended by sowing the male parents twice, about 5 and 10 days later than the sowing date of the female parents. Female parents were detasseling on July 7 and harvested on September 16-18 in both years.
Basal fertilizers were spread on the soil surface and plowed evenly into the soil before soil mulching with plastic film. In 2015 the amounts of basal fertilizers were 216 kg· ha −1 N, 105 kg· ha −1 P, 62 kg· ha −1 K and 37.5 kg· ha −1 magnesium-zinc sulfate (Zn + MgO 3 17%). A further 207 kg· ha −1 N as urea was top-dressed and then dissolved and infiltrated into the soil with irrigation water at the first irrigation event. Also, foliar spraying with KH2PO4 was conducted at the jointing stage. In 2016, only 205 kg· ha −1 N and 101 kg· ha −1 P were applied as basal fertilizers, and two top-dressings were fertigated with 139 kg· ha −1 N on June 23-24 and 70 kg· ha −1 N on July 15-16.
In 2015, border irrigation (irrigation controlled or directed by short dikes around the areas treated) was used. Five irrigation events were used with the irrigation units marked with a red dashed line ( Fig. 1(a)). Outside the red dashed line, only the first four irrigation events were used. Also, winter irrigation with more than 150 mm was applied in the last winter to control salinity and provide adequate soil moisture for seedling emergence in 2015. Winter irrigation in 2015 was canceled because drip irrigation (surface drip irrigation) was planned in 2016. Instead, irrigation was applied with a drip irrigation system just after the female parents were sown to assure seedling emergence in 2016. Also, another seven irrigation events were applied during the growing season. The total amount of irrigation was 679 mm (winter irrigation was not included) under border irrigation in 2015 and 452 mm under drip irrigation in 2016. Irrigation schedules are listed in Table 1.  Note: Irrigation amount is the area-weighted average for the irrigated units. The total is the area-weighted average of total amount of irrigation during the growing season, and winter irrigation is omitted.

Data collection and preprocessing
In 2015, 185 sampling locations were used, comprising 169 sampling points on an approximate 15-by-15m grid and another 16 sampling points placed randomly with the distance to the grid points less than 15 m.
In 2016, 184 sampling points were used in the same field, comprising 156 sampling points on an approximate 15-by-15-m grid and another 28 sampling points placed randomly with the distance to the grid points less than 15 m. A further 18 sampling points, labeled with ID numbers 185-202 ( Fig. 1(b)), were added on an adjacent field under the same agronomic management in 2016. These additional sampling points were not included in the statistical analysis as there are no corresponding samples from 2015; however, they were used in searching the maximum yield, clustering soil texture zones and fitting boundary lines. Sampling points were adjusted slightly when too close to the field ridge. The coordinates of each sampling point were obtained using a GPS tracker (Trimble Recon, Sunnyvale, CA) and then projected into a Gauss-Kruger projected coordinate system (Beijing 1954, 3 Degree GK CM 102E), and the elevation was obtained with an optical level (Fig. 1).
Each sampling site was a 3.3 m (two mulched strips) by 3.1 m (12 plants per row) area. Emergence was counted to calculate plant density. Averaged emergence was 74% under border irrigation and 81% under drip irrigation. All maize ears were harvested in the representative area then weighed for fresh mass (Mfresh), ten ears were oven-dried at 80°C to determine ear water content (', according to Eq. (1)) and the mass ratio of grain to ear (p) was calculated. The grain yield was calculated on dry matter basis (Eq. (2)). Then the relative yield (RY, expressed as actual crop yield divided by the observed maximum yield in each year. Equation (3) was used in the analysis to avoid the influences of meteorological and cultivar differences between the two production years.
' = (mfreshmdry)/mdry (1) where ' is the ear water content. mfresh and mdry are the fresh and dry mass of ten ears, respectively, g. Yield is the grain yield in dry mass, Mg· ha − . Mdry and Mfresh are the dry and fresh mass of all the ears harvested on the representative area, kg. p is the mass ratio of grain to ear. A is the harvested area, m 2 . 10 is the coefficient converting grain yield from kg· m − to Mg· ha − . RY is the relative yield. Yieldmax is the maximum of all the observed grain yield in each year.
Soil samples were collected four times each year to a depth of 100 cm at 20 cm depth intervals in the representative area at the start and end of the growing season, and on June 8 and July 10, 2015 (about one week after irrigation), and June 30 and August 6, 2016 (one day before irrigation). Four soil samples at the same soil depth were mixed and processed to determine soil properties. All samples were used to determine soil gravimetric water content (') and nitrate-and ammonium-nitrogen (NO3 --N and NH4 + -N).
All samples collected at the start of each growing season were used to determine soil texture (clay, silt and sand fractions) and the volume-weighted mean diameter (Dvol). Also, samples from 0 to 40 cm depth at the start and the end of each growing season were used to determine pH in water, electrical conductivity (EC), organic matter (OM) content, total nitrogen (TN), total phosphorus (TP), soil available phosphorus (AP), and soil available potassium (AK) contents. Bulk density (BD) in the top 100 cm of the soil profile was determined on the 130 subsample sites with drip irrigation. These sampling sites included ID numbers 1-155 (odd numbers only), all ID numbers 157-202, and other ID numbers 2, 32, 34, 48, 76 and 110.
Soil pH was determined in a 1:2.5 soil-water suspension with a pH meter after shaking for 2 min and standing for 30 min. Soil electrical conductivity was determined in a 1:5 soil-water extract using a conductivity meter after shaking for 3 min. Soil OM was determined by the Walkley-Black rapid dichromate oxidation method [23] . TN was determined by Kjeldahl digestion and continuous flow analysis (AutoAnalyzer 3, Bran + Luebbe, SEAL Analytical GmbH, Norderstedt, Germany). TP was determined by digesting soil samples with an acid mixture (H2SO4-HClO4) and continuous flow analysis. NO3 --N and NH4 + -N were determined by extracting samples with 1 mol· L −1 KCL solution and continuous flow analysis. AP was determined by extracting samples with 0.5 mol· L −1 NaHCO3 solution and the molybdate colorimetric method. AK was determined by extracting samples with 1 M NH4OAc solution and flame photometry. Soil gravimetric water content (') and bulk density (BD) were determined by the oven drying method. Soil texture was determined with a laser particle size analyzer (Mastersizer2000, Malvern Instruments Ltd., Malvern, UK).
Soil available nitrogen (AN0-100 cm) is expressed as the average of available nitrogen (nitrate-and ammonium-nitrogen summed) in top 100 cm of the soil profile. Other soil properties (pH, EC, OM, TN, TP, AP, AK) were the averages in the top 40 cm of the soil profile.
Soil available water content in the top 100 cm of the soil profile (PropAWC) is expressed as the fraction of the total available water content (the moisture held between field capacity and permanent wilting point).
PropAWC was calculated using Eq. (4) to Eq. (8). The permanent wilting point (WP) was estimated using the pedotransfer functions (PTFs) described by van den Berg et al. [24] (Eq. (6)). Field capacity (FC) and BD were estimated using the PTFs when no observed data were available (Eq. (7) and Eq. (8)). (6) FC = pr1 + pr2 × Clay + pr3 × Clay 2 + pr4 × Silt + pr5 × Silt 2 + pr6 × Clay × Silt (7) BD = pr1 + pr2 × Clay + pr3 × Clay 2 + pr4 × Silt + pr5 × Silt 2 + pr6 × Clay × Silt (8) where SWCActual, SWCWP, and SWCFC are the soil water storage in the top 100 cm of the soil profile when soil was sampled, soil potential reached −1.5 MPa, and field capacity reached, respectively, mm. i (i = 1, 2, ⋯, 5) is the number of the soil layer 0-20, 20-40, 40-60, 60-80 and 80-100 cm. 'Actual and Actual are the soil gravimetric and volumetric water contents when the soil was sampled. FC and WP are the volumetric water contents when field capacity reached, and soil potential reached −1.5 MPa, respectively. Clay and silt are the fractions of soil particle size ˂ 2 m and between 2 and 50 m (when clay was 30%, then 30 would be used in Eq. (6) to Eq. (8)), respectively. BD is the soil bulk density, Mg· m −3 . pr1 to pr6 are the parameters of the PTFs estimated by nonlinear regression with XLSTAT version 2014.5.03 (Table 2). When estimating the parameters of FC (Eq. (7)), data was collected from other research at the same study area (sample size: n = 26; the measured FC ranged from 18.4% to 38.2%). For the parameters of BD at different soil depths (Eq. (8)), both the measured data set in this study (bulk density of the 0-100 cm depth soil profile with 20 cm depth intervals from 130 sampling sites) and collected data set from other studies (5 sampling sites) were used (sample size: n = 135). . For BD at different soil depths, both the measured data set from the present study (bulk density of 0-100 cm depth of soil profile with 20-cm depth intervals from 130 sampling sites) and a data set from other studies (5 sampling sites) were used (sample size, n = 135).

Data analysis 2.3.1 Boundary line analysis and quantification of yield gaps
The yield gap is commonly defined as the difference between yield potential and actual yield achieved in the field. Yield potential, also called potential yield, is defined as the yield of a cultivar when grown in an environment to which it is adapted, with nutrients and water non-limiting and other stresses effectively controlled [8] . However, it is difficult to measure the yield potential as it can be difficult to provide an entirely stress-free testing environment. Yield potential is therefore usually estimated by the maximum yield of sizable samples under optimum management. As the main concern of this study was closing yield gap as determined by optimizing field management, and RY (each season) was used to avoid the influence of meteorological difference between the two growing seasons, we define the yield gap (Total) as the difference between the maximum RY (RYmax) and the RY achieved: BLA has been described in previous studies [18][19][20] . The boundary line is usually determined by fitting a line through the upper boundary points of data scatter (yield as the dependent variable and yield-limiting factors as the independent variables). It is assumed that the boundary line represents the dependent variable influenced only by the single independent variable in the absence of any other limiting factors if the data set is sufficient [21] . Here, a large area was affected by plant density deficiency. Thus the boundary lines of other yield-limiting factors (except for plant density) were inevitably influenced by two or more yieldlimiting factors in most cases (at least one was the independent variable, and another was plant density). This conflicts with the assumption of the BLA. We therefore gave priority to determining the boundary line of plant density. We used RY as the dependent variable for the boundary line of plant density, but adjusted RY as the dependent variable for boundary lines of other yield-limiting factors by eliminating the reduction effect of all the yield-limiting factor(s) that were analyzed (Fig. 2). The procedure is described below. Supposing crop yield (we used RY) was the outcome that potential yield (as the potential yield was difficult to measure, we used the observed maximum yield instead, RYmax = 1) under specific climatic conditions constrained by yield-limiting factors (Xi, i = 1, 2, ..., n) and unknown factors (, including the reduction effect of others unknown yield-limiting factors, random error and interaction effects between factors Xi if they exist). The reduction effect of each yield-limiting factor was expressed as a yield correction coefficient, which was a function of each factor f(X1), f(X2), ..., f(Xn), f(). This may be expressed as follows.
As the assumption of BLA, the boundary line for the factor X1 represents the dependent variable influenced only by a single independent variable X1 in the absence of any other limiting factors when the data are sufficient. Then Eq. (10) can be simplified as follows, as all the reduction coefficient function of factors (X2, X3, ..., Xn, ) are equal to 1: (11) where RYmax is the maximum RY, BLA(X1) is the boundary line function of factor X1. According to Eq. (10) and Eq. (11), the reduction effect of other factors (X2, X3, ..., Xn, ) can be expressed as: When factor X1 is managed to the optimum level, the reduction effect of X1 can be eliminated, i.e., the yield correction coefficient function f(X'1) ≡ 1 (the symbol f(X'i) is used to represent the reduction function after optimizing factor Xi distinguished from the yield correction coefficient function f(Xi) before optimizing factor Xi; RY'Xi is used to represent the predicted RY when the first i yield-limiting factor(s) X1, X2, ..., Xi is/are managed to the optimum level(s). Then the predicted RY (RY X 1 ′ ) after optimizing factor X1 can be expressed as: The difference (∆RY X 1 ) between RY and RY X 1 ′ can be explained as the yield gap determined by single factor X1 when factor X1 was managed to the optimum level: Then the RY X 1 ′ can be used as the new dependent variable to determine the boundary line of the following factor X2 free from the yield reduction effect of the previous factor X1. Similarly, the boundary lines of other factors can be determined in sequence. The general expressions for factors Xi (i = 2, 3, ..., n) are as follows: where RY X −1 ′ is the predicted RY with the reduction effect of previous factors (X1, X2, ..., Xi-1) eliminated by managing them to the optimum levels, which was also the dependent variable when fitting the boundary line of factors Xi; RY X ′ is the predicted RY with the reduction effect of the first i factors (X1, X2, ..., Xi-1, Xi) eliminated by managing them to the optimum levels; BLA(Xi) s the boundary line function of factor Xi; ∆RY X is the yield gap determined by a single factor Xi.
When all the boundary lines of the yield-limiting factors considered had been determined, the yield gap of unknown factors  was estimated by Eq. (17): where RY was the relative yield achieved in the field, RYmax was the maximum RY; ∆RY X 1 , ∆RY X 2 , ..., ∆RY X and ∆RY ε were the yield gaps determined by single yield-limiting factors Xi (i = 1, 2, ..., n) and unknown factors (), respectively.
The main steps for quantifying the yield gaps by multiple yield-limiting factors are as follows: (1) Determining the sequence of fitting the boundary lines for different yield-limiting factors. The factor that was more likely to be the dominant yield-limiting factor and with less uncertainty should be given priority (plant density in our study). However, the order should be regulated if the upper boundary points, when fitting the boundary line of the current factor, are commonly outside the optimum levels of the latter factors (RY is constrained by the latter factors).
(2) Fitting the boundary line of every single factor (Xi). We used yield-limiting factor Xi as the independent variable and the RY as the dependent variable for the first factor X1, but we adjusted relative yield (RY X 1 ′ ) as the dependent variable for the subsequent factors Xi + 1 (i = 1, 2, ..., n). Also, if a point of the x-value of any yield-limiting factors (Xi + 1, ..., Xn) (after the current analysis factor Xi) is not within its optimum range, it is also deleted as an outlier. Finally, fitting an appropriate model through the upper boundary points for every single factor. More details about the searching procedure are given by Schnug et al. [19] and we conducted the procedure using the R version 3.5.2 statistical package [25] .
To describe the reduction intensity of different yield-limiting factors we defined the optimum range and acceptable range. The optimum range is defined as the independent variable interval with RY of the boundary line (maximum attainable RY) equal to 100%, and the acceptable range as the independent variable interval with RY of the boundary line ˃ 95%. The independent variable between the acceptable range was supposed to be acceptable stress and outside the acceptable range was regarded as non-ignorable stress.
The results of the BLA are shown in Fig. 2 and Table 3. Soil properties, comprising soil organic matter, total nitrogen, total phosphorus and available phosphorus, were not included in the BLA, as the boundary lines of these properties were approaching the horizontal line RY = 1.    I  38  43  6  12  1  44  56  0  0  0   II  2  6  92  --3  9  87  --III  0  18  68  12  2  0  0  56  39  5   IV-1  5  5  89  -------IV-2  -----11  16  73  --V  -----2  7  91  --VI  --86  12  2  --78  19  3   VII  --87  12  1  --100  0  0 Note: ID, order of fitting the boundary line; Num, female plant density; AK, soil available potassium (at 0-40 cm) at the end of the growing season; AN, soil available nitrogen (at 0-100 cm) at the end of the growing season; PropAWC, soil available water content (in the fraction of total available water content); EC, soil electrical conductivity at the start of the growing season; pH, soil pH at the start of the growing season; optimum range, the independent variable interval that relative yield of the boundary line equaled 100%; acceptable range, the independent variable interval that relative yield of the boundary line was ˃ 95%. When determining the parameters of boundary lines with R version 3.5.2 [25] , lm() function was used for the factor Num (ID: I), and both the t-test and F-test were conducted; nls() function was used for other factors (ID: II-VII), and only the t-test was conducted.
※ The value of t1 or t2 was set manually to the minimum or maximum of the independent variable with the adjusted RY approaching 1.000 (> 0.999) if the fitting program failed for too many parameters or too small a sample size. ¶ Soil sampling dates were June 8, 2015, August 6, 2016 and September 13, 2016. § Nutrient convertion coefficient: 1 mg· kg −1 = 6.02 kg· ha −1 (at 0-40 cm soil), 1 mg· kg −1 = 14.1 kg· ha −1 (at 0-100 cm soil); † Yield reduction level is expressed as the maximum attainable relative yield (relative yield of the boundary line) that was reached, and 0.95 or 1 represents the relative yield of the boundary line that reached 0.95 or 1. This is consistent with the lower or upper limits of the acceptable range and optimum range.

Hierarchical cluster analysis
To analyze the relationship among soil moisture, nutrients and crop yield at the same soil zones, hierarchical cluster analysis was conducted using XLSTAT version 2014.5.03. Soil Dvol in the top 100 cm of the soil profile (20-cm-depth intervals) was selected as the clustering variable. Ward's method was chosen as the agglomeration method, and Euclidian distance was used to measure the dissimilarity coefficient. As a result, the whole field was divided into four soil texture zones (Fig. 3). Also, the box plots of the remaining soil NPK in different soil texture zones are shown in Fig. 4.

Interpolation method
The inverse distance weighted method was used for interpolation when plotting the spatial distribution of yield and yield-limiting factors [26] .
] ⁄ (18) where, ( 0 ) is the interpolated value at point 0 , ( )is the measured value at sampling points (i = 1, 2, ..., m), di is the distance between 0 and , p (ranging from 1 to 100) is the power value of the distance di. For the parameters settings when interpolating RY and yield-limiting factors, both the searching major and minor semiaxes were 15 m, minimum and maximum neighbors of each sector were 1 and 2, four-sector with 45° offset was selected, and an optimum power value p was used according to cross-validation by minimizing the root mean square prediction error. The value of p was 1.75, 1.96, 3.94, 1.58, 1.00, 1.00, 1.00 and 1.00 for RY, RY', Num, AK, AN, PropAWC20150608, pH and EC in 2015, respectively, and p was 1.00, 1.00, 1.00, 1.00, 1.24, 1.00, 1.00, 1.00 and 1.40 for RY, RY', Num, AK, AN, PropAWC20160806, PropAWC20160913, pH and EC in 2016. However, to maintain the soil texture of the predicted point consistent with the nearest observed point when interpolating soil texture zone, one sector and one neighbor for each sector were selected, and p could be any value. Interpolation was conducted using ArcGIS 10.4 and the interception plots are shown in Fig. 5 and Fig. 6. (a)

Hierarchical cluster analysis and zonal statistics of soil NPK
The hierarchical cluster analysis results are shown in Fig. 3. The whole field was divided into four soil zones. The soil texture according to the USDA soil texture classification system of the cluster centroid in each soil zone is listed in Table 4. The area percentages of Zones 1, 2, 3 and 4 were 41%, 30%, 16% and 13%, respectively. Note: Soil texture, the soil texture of the class centroid, calculated by averaging the sand, silt and clay contents at 0-60 cm or 60-100 cm depth in each soil zone, according to the USDA soil texture classification system; AN0-100 cm, AP0-40 cm and AK0-40 cm, the remaining soil available nitrogen at 0-100 cm, available phosphorus at 0-40 cm and available potassium at 0-40 cm depths at the end of the growing season. ¶ Significance test of difference by independent-samples t-test using the SPSS 20.0 software package (P < 0.05).
Soil nutrients in the four different zones under border and drip irrigation at the end of the growing season are summarized in Table 4. Throughout the field the remaining soil available nitrogen at 0-100 cm depth increased from 16.3 mg· kg −1 N under border irrigation to 30.4 mg· kg −1 N under drip irrigation, and available phosphorus at 0-40 cm depth increased from 12.1 mg· kg −1 P under border irrigation to 14.8 mg· kg −1 P under drip irrigation. Both the soil available nitrogen and available phosphorus under drip irrigation were significantly higher than in the same soil texture zones under border irrigation (P < 0.05). Accumulation of soil available nitrogen and available phosphorus is also shown in box plots (Fig. 4) and spatial interception plots (Fig. 5). However, throughout the field the difference in remaining soil available potassium at 0-40 cm depth was not significant under the border and drip irrigation, with an average of 104.3 mg· kg −1 K under border irrigation and 101.9 mg· kg −1 K under drip irrigation. Only in Zone 2 was the difference in remaining available potassium significant (P < 0.05) with average values of 103.5 mg· kg −1 K under border irrigation and 95.4 mg· kg −1 K under drip irrigation (Table 4).

Grain yield and water consumption under both irrigation methods
Maize grain yield and water consumption are summarized in Table 5. The maximum yields were 6.  Table 5). The difference in RY between the two irrigation methods was not significant in the same soil zone (P < 0.05). In the same growing seasons Zone 4 was the zone with the lowest RY under both irrigation methods. However, the only significant difference was between Zones 2 and 4 under border irrigation (P < 0.05). Drip irrigation showed a marked advantage in increasing RY at area A (mainly belonging to Zones 3 and 4 with coarser-textured soil) and area B (belonging to Zone 1 with finer-textured subsoil) where soil nitrogen and salinity were both at higher levels. However, drip irrigation decreased RY at area C ( Fig. 6(a)). Drip irrigation decreased the yield variability with a coefficient of variation of 12.1% throughout the field compared with 18.8% under border irrigation. Note: RY, the relative yield; CV, the coefficient of variation; Yieldmax, the observed maximum yield (dry grain); P and I, the total precipitation and irrigation water during the growing season, respectively; ET, crop evapotranspiration by the eddy covariance method. ¶ The significance test of difference by independent-samples t-test using the SPSS 20.0 software package (P < 0.05). † ET was the eddy covariance result in both years, according to other literature on the same field. ET from April 15 to September 16, 2015 was the eddy covariance results from Qin et al. [27] ; ET from April 25 to September 16, 2016 was adjusted from the result of He et al. [28] by plus or minus the ET of missing or extra days.
Irrigation water input during the growing season was 679 and 452 mm under border and drip irrigation, respectively. Drip irrigation saved 227 mm irrigation water during a drier growing season compared with border irrigation and this was 44% of the average ET throughout the growing season (ET, the eddy covariance measurements of the same field by Qin et al. [27] and He et al. [28] ). Also, 91% of the total water input (precipitation plus irrigation) was converted into crop ET under drip irrigation but the corresponding percentage was only 61% under border irrigation (Table 5).

Yield-limiting factor analysis and yield gap quantification
The parameters of the BLA method, the optimum range and the acceptable range, are listed in Table 3. We recommend the optimum range as the thresholds for field management and yield-limiting factors outside the acceptable range should be given priority in closing yield gaps. Also, the area proportion of different stress intensity by different yield-limiting factors is also summarized in Table 3. The area in which yield was constrained by each yield-limiting factor (deficiency or excessive) is delineated in the interpolation maps (Fig. 6).
Yield gaps determined by yield-limiting factors are summarized for different soil texture zones in Table 6. Throughout the field the yield gap determined by all the currently considered yield-limiting factors () was 7.9% under border irrigation and 9.7% under drip irrigation (Table 6). However, yield gaps differed with soil zone and irrigation method. Zone 1 with finer-textured soil had the minimum yield gaps determined by currently considered factors and Zone 4 with coarser-textured soil had the maximum yield gaps under both irrigation methods. The yield gaps determined by currently considered factors () ranged from 5.9% at Zone 1 to 12.4% at Zone 4 under border irrigation, and 7.0% at Zone 1 to 13.8% at Zone 4 under drip irrigation. The RY was expected to increase from ~70% at present to ~80% (range 76% to 84% in different soil zones under both irrigation methods) if the currently considered yield-limiting factors were managed to the optimum ranges. The yield gaps determined by unknown factors () throughout the field were 21.6% under border irrigation and 19.2% under drip irrigation. Zone 1 had the maximum unexplained yield gap with values of 24.2% under border irrigation and 21.8% under drip irrigation.  Note: RY, the observed relative yield; RY', the expected relative yield when all the factors considered (unknown factors  not included) were managed to the optimum ranges; RY' = RY + ; , the total yield gap of the current factors considered (unknown factors  not included);  = Num + PropAWC + AN + EC + AK + pH; Num, PropAWC, AN, EC, AK, pH, , yield gap determined by female plant density, soil available water content, soil available nitrogen, soil electrical conductivity, soil available potassium, soil pH, and unknown factors, respectively; RY' +  = 1.
Irrespective of irrigation method the first three yield-limiting factors were plant density (Num), soil available water content (PropAWC), and soil available nitrogen (AN) throughout the field ( Table 6). The yield gaps determined by plant density (Num) were 4.1% and 4.9% in 2015 and 2016, respectively. The slightly larger yield gap determined by plant density in 2016 may be explained by the larger areas constrained by female plant deficiency that occurred almost throughout the study area compared with 81% of the study area in 2015 (Table 3). Throughout the field the yield gaps determined by soil available water content (PropAWC) were 1.4% under border irrigation and 2.5% under drip irrigation ( Table 6). The larger yield gap determined by soil available water content under drip irrigation was the result of the larger area affected by water stress. Yield gap determined by soil available water content was considerable in all soil textures except Zone 1 under drip irrigation with values of 3.1%, 4.5% and 7.5% in Zones 2, 3 and 4, respectively. However, only Zone 4 had a large yield gap of 6.0% determined by soil available water content under border irrigation (Table 6). In Zone 4, soil available water content had become the predominant yield-limiting factor under both irrigation methods ( Table 6). The yield gap determined by soil available nitrogen (AN) was similar throughout the field with values of 1.1% and 1.0% under the border and drip irrigation, respectively. However, the cause of the yield gap was different. Both soil nitrogen deficit and excess occurred under border irrigation, with soil nitrogen deficit occurring mainly in the coarser-textured Zones 3 and 4, and soil nitrogen surplus occurred mainly in finer-textured Zone 1 (Fig. 6(c)). However, only soil nitrogen surplus occurred under drip irrigation, mainly in Zone 1 ( Fig.  6(h)). Accumulation of soil available nitrogen led to the area proportion affected by excessive nitrogen increasing from 14% of the total area under border irrigation to 44% under drip irrigation (Table 3).
Although other yield-limiting factors considered had little influence on yield gaps throughout the field, they need to be considered in specific soil texture zones. For example, a large area at Zone 1 was affected by soil salinity at the early growth stages, accompanied by soil nitrogen surplus (Fig. 6((c-d),(h-i))), and Zones 3 and 4 with coarser-textured topsoil were also affected by available potassium deficit ( Table 6, Fig.  6(b,g)).

Thresholds of yield-limiting factors for field management
Studies have been conducted to recommend the thresholds for field management with different methods such as BLA [18,19] , the Cate-Nelson graphical procedure [29] , or least square regression between crop yield and yield-limiting factors. The results may differ slightly from each other. As these ranges are the results of analyzing the relationship between the yield of a specific crop and yield-limiting factors at specific sampling time, they would be misleading if the sampling time or crop species were ignored.
Here, the optimum range of female parents was (7.6-7.8) × 10 4 plants ha −1 , a lower value than the results of Jiang et al. [22] who found that maximum yield was achieved at a female plant density of (9.4-10.6) × 10 4 plants ha −1 (a further 20% male parents were not included). The optimum range was reached earlier in our study because the BLA method is concerned with the upper boundary points of scatter, and factors other than plant density were at optimum levels. Also, the narrow optimum range was due to few data being distributed at higher plant density in our study. However, this range is consistent with the optimum plant density of (9.7-16.4) × 10 4 plants ha −1 for hybrids released between 1965 and 2010 [30] if the male parents are taken into account. A slightly higher plant density is recommended for increased yield stability as modern hybrids are tolerant to high plant density [30] .
According to the BLA of soil available water content in 2016 we recommend 37% to 51% of soil available water content as the lower limit for irrigation management. This is consistent with the recommendation of 35% to 55% (the maximum allowable water depletion is 45% to 65% of total soil available water) to meet full water requirements [31,32] . Soil available water content exceeded 100% at some sampling sites (Fig. 2), especially under border irrigation, so that soil water content exceeded field capacity. This is related to the common occurrence of finer-textured subsoil in the field, and this may lead to a longer time required when downward flux becomes negligible after large rates of irrigation [33] . The linear-plateau model was used in fitting the relationship of yield response to soil available water content, as there was no apparent yield reduction by excessive irrigation in our study. However, large rates of irrigation under border irrigation should be avoided because of the risk of hypoxia [34] .
Soil nutrients at the end of the growing season were used to indicate thresholds for fertilizer management, as soil nutrients during this period tended to be lower than that of the earlier growth stages. Another consideration was that residual soil nitrogen was more stable than at other growth stages, as nitrogen was applied more than once and complex movement with soil water. The optimum range of average soil available nitrogen in the 100 cm was 6.04 to 26.1 mg· kg −1 (Table 3). However, high nitrogen inputs should be avoided because of the risk of nitrate leaching [35] and the negative influence of salinity on germination and early seedling establishment [36] . The lower limit of the optimum range in our study is comparable to the critical level of 8 mg NO3 --N kg −1 (bulk density assumed to be 1.49 Mg· m −3 ) [29] . Soil phosphorus limitation did not occur according to the BLA method, with soil available phosphorus ranging from 5 (at the end of the growing season) to 51 mg· kg −1 P (at the start of the growing season when maize was sown). However, soil available phosphorus outside this range should be considered for phosphorus deficiency or excess, which may lead to deficits of micronutrients such as Zn or Fe [37] . We recommend the lower limit of 81 mg· kg −1 K for soil available potassium which is comparable to the lower limit of 90 mg· kg −1 K according to other studies [38] .

The current situation and field management problems
Here, irrespective of whether border or drip irrigation was used, current field management achieved a RY of 70% for the whole field level, which was at a high level given the yield plateau 75% to 85% of the potential yield in the literature [8] . However, it was achieved at the cost of high resource consumption.
Under border irrigation only 61% of the total water input (precipitation plus irrigation) was converted to crop ET (Table 5), so that that the remaining water was wasted in the form of deep percolation or water storage in the root zone at the end of the growing season. Drip irrigation showed a substantial advantage over border irrigation in saving irrigation water. Other studies also report that drip irrigation required 30% to56% less water with an equal or higher yield compared with the sprinkler or border irrigation [39,40] . The poor performance of border irrigation is well known, as irrigation efficiency is significantly affected by field microtopography even with laser leveling [41] . Also, application efficiency (Ea, the percentage of irrigation water contributing to target irrigation depth) was only 49% to 82% under border irrigation [42,43] , but reached ~90% under both sprinkle and drip irrigation [42] . Expanding drip irrigation may ease the water scarcity situation in north-west China.
Using drip irrigation calls for corresponding changes in fertilizer management because current high fertilizer inputs greatly threaten resources and the environment. Fertilizer overuse in China has been widely reported in recent years [1,2,10,17] . Here, despite the application of ˃ 400 kg· ha −1 N, nitrogen deficit occurred at some locations in the coarser-textured Zones 3 and 4 under border irrigation (Fig. 6c). However, drip irrigation succeeded in maintaining soil nitrogen within the root zone, as top-dressing nitrogen was fertigated with greater frequency but lower quantity at each event and there was less nitrate leaching due to reduced irrigation water volume. As a result, nitrogen deficit was eliminated. However, only one growing season after using drip irrigation the area in which remaining available nitrogen at 0-100 cm depth exceeded the upper limit of the optimum range (26.1 mg· kg −1 or 369 kg· ha −1 N) increased to 44% under drip irrigation compared to 14% under border irrigation (Table 3). This upper limit reached 1.8 times the maize demand of 200 kg· ha −1 N throughout the growing season [31] , and crop yield began to be constrained by nitrogen surplus (Fig. 2(c)). In these areas, reducing nitrogen application would save nitrogen, reduce the energy used by the fertilizer industry [1] , decrease agricultural impact on drinking water, atmosphere or ecosystem [44] , increase the partial factor productivity of nitrogen [14,45] , and also increase crop yields.

Measures and suggestions for closing yield gaps and sustainable agriculture
We recommend the optimum range as the threshold for field management as discussed in Section 4.1. Any practices that can maintain the yield-limiting factors within the optimum ranges should be taken into consideration. Here, plant density was the predominant factor limiting yield throughout the field irrespective of border or drip irrigation. This is consistent with the conclusion of Chen et al. that plant density makes a slightly higher contribution to the yield gap than other management practices [14] . Both the emergence and sowing rate determined plant density. Emergence was influenced by numerous variables such as soil temperature, soil moisture, soil air quality, soil burial, compaction, soil fertility and salinity [36,46] , which can be manipulated indirectly or directly through field management [46] . Any practices that can guarantee emergence should be considered, including delaying sowing date to avoid low temperatures at the germination-emergence stage, ensuring suitable soil moisture content, lowering soil salinity and even increasing the sowing rate if reduced emergence is unavoidable.
Soil available water content and available nitrogen were the other two important yield-limiting factors following plant density, and they can be manipulated through irrigation and nitrogen management. Modern irrigation methods such as drip irrigation should be given special priority, because of the advantages in irrigation management even in soils with small ranges of non-limiting water [47] . Other advantages include lower consumption of irrigation water [40] , improved soil physical properties [48] , higher soil aeration [34] , and reduced weed growth [40] compared to border irrigation. Also, balanced fertilizer management was needed to maintain a high yield level with lower fertilizer input. However, it seems to be difficult or impossible to achieve the optimum range for the whole field with current uniform irrigation and fertilizer management, as these factors are closely related to soil texture and spatial variability commonly occurs in soils. Hence, site-specific management is needed. Frequent irrigation with a smaller volume each time would be favorable in the coarser-textured soil zones as irrigation frequency had a large influence on grain yield even using the same amount of total irrigation [49] . In the finer-textured soil zones in which nitrogen and salt accumulation occurred, low fertilizer rates are recommended. Utilizing site-specific management can improve nitrogen management and increase nitrogen use efficiency, and reduce nitrogen losses to the environment [50] , which was also true for irrigation management. This supports the use of site-specific management on larger farms and at regional scale.
According to the yield gap analysis of different yield-limiting factors, multiple goals of 8% to 10% RY increase at the whole field scale, reducing agricultural inputs and promoting agricultural sustainability may be realized by managing the important factors within the optimum ranges. However, this relies on optimizing agronomic management [45] , balanced fertilizer management, modern irrigation methods, sitespecific management [50] , and other technical support of precision agriculture, such as the timely acquisition of field information, variable-rate fertilization and irrigation techniques [51] . However, further yield increases appear not to be economical or even possible when a RY of 78% to 81% was attained throughout the field but a further yield gap of 19% to 22% RY occurred due to unknown factors ( Table 6). Unknown factors may include weeds and pests, variation in seed vigor and plant health status, spatiotemporal variability of soil water and nutrients, and random error. Sustainable agriculture must be practiced to ensure future food, water resources and environmental security.

Conclusions
Using a field investigation and yield-limiting factor analysis, the main conclusions of this study are as follows.
(1) Whatever the irrigation method, current field management can achieve a RY of 70%. However, drip irrigation consumed less irrigation water and decreased yield variability compared to traditional border irrigation. Drip irrigation is recommended to minimize water scarcity in arid north-west China. (2) The optimum ranges of different yield-limiting factors can be used as the thresholds for field management. About 8% to10% RY can be increased by managing the considered yield-limiting factors within the optimum ranges by relying on the theoretical and technical support of precision agriculture.
When an RY of 78% to 81% is reached, further yield increase appears to be uneconomic or impossible. (3) Plant density, followed by soil available water content and remaining available nitrogen, explained the major causes of yield gaps by considered yield-limiting factors (not including the unknown factors). Also, in coarser-textured soil zones soil available water content explained a greater yield gap than the other two factors. Hence, plant density, irrigation and nitrogen management must be given special emphasis in closing yield gaps and developing sustainable agricultural systems.