Variables influencing yield-scaled Global Warming Potential and yield of winter wheat production

Field Crops Research 227 (2018) 19–29

Contents lists available at ScienceDirect

Field Crops Research journal homepage: www.elsevier.com/locate/fcr

Variables influencing yield-scaled Global Warming Potential and yield of winter wheat production

T

Elżbieta Wójcik-Gront Department of Experimental Design and Bioinformatics, Warsaw University of Life Sciences, Poland Warsaw University of Life Sciences – SGGW, Nowoursynowska 159, 02-776, Warszawa, Poland

A R T I C LE I N FO

A B S T R A C T

Keywords: GHG emission Carbon footprint Yield variability Multivariate analysis Regression trees

The main goal of this study was to determine the major drivers of variation for the yield-scaled Global Warming Potential (GWP) and yield of winter wheat in Poland focusing on environmental, genetic and management variables. The yield-scaled GWP is the GWP calculated per grain unit expressed in kg CO2 equivalent kg−1 yield. The analysis was performed using multivariate statistical methods: CART (classification and regression tree) and RF (random forest). This is the first study in Poland to focus on variables besides those used in GWP calculations and those influencing the yield-scaled GWP variability of winter wheat production. Identification of these variables contributes to the creation of more environmentally friendly wheat cropping systems. In this study, regression-tree based analysis revealed that soil quality and water availability during crucial stages of plant growth are the most influential input variables of the yield-scaled GWP and yield variability in winter wheat production. Not surprisingly, environmentally favorable conditions for wheat growth contribute to its high yields yet require less intensified agronomic management. The strong influence of water availability in June and July, at the end of plant growth, applies to the undesirable effect of excess water leading to plant diseases which result in lower yield. N fertilizer has a strong effect on GWP of winter wheat production. However, this study also shows that nitrogen is not one of the most influential variables of wheat yield variability. Thus, increasing its average use in Polish environmental conditions from 107.3 kg ha−1 to 147.3 kg ha−1 might not increase yield sufficiently for its use to be justified. More important variables of yield variability were the use of fungicides and growth regulators, which are applied at much smaller rates (1.7 and 0.8 kg ha−1, respectively) than N fertilizer and positively influence efficient winter wheat production.

1. Introduction

emissions during all off-farm and on-farm practices, a Life Cycle Assessment (LCA) (Consoli et al., 1993) should be applied to provide reliable information on the Global Warming Potential (GWP) per unit of crop production area. However, for planning more environmental friendly cereal cropping systems an assessment of both GHG emissions and crop yield per unit of crop production area is essential. Hence, GWP can also be related to the amount of a crop and be expressed per unit of grain yield. This is yield-scaled GWP (van Groenigen et al., 2010). One of the most important cereals in world food production is wheat (881 million tons in 2016, FAOSTAT, 2018). Poland is a substantial contributor to worldwide wheat production. Especially popular is winter wheat (Triticum aestivum L). The Statistics Poland (GUS) reports that from 1999 to 2015 annual winter wheat production increased by 28%, from 9.0 to 11.6 million tons, with average grain yield increase from 3.5 to 4.8 t ha−1, while the wheat production area was reduced by 9%. Still, this relatively low average yield is mainly due to the economical limitations of many Polish farms. When nutrient supplies in the

Cereal-based cropping systems contribute to different greenhouse gas (GHG) emissions. These emissions come from manufacturing, processing and applying fertilizers and pesticides, from fossil fuels required in field operations for soil cultivation, spreading agrochemicals, harvesting and machinery production. The majority of these GHG emissions come from nitrogen fertilizer usage. There are not only GHG emissions associated with N fertilizer manufacture but also direct and indirect GHG emissions from N application to soil. Direct emissions of N2O result when excess N fertilizer undergoes incomplete denitrification. In addition, indirect N2O emissions can result outside of the farm boundary when nitrate lost from the field through leaching is subsequently incompletely denitrified (Charles et al., 2006; Snyder et al., 2009; Williams et al., 2010; Linquist et al., 2012). GHGs cause chemical changes to the atmosphere and therefore contribute to Global Warming (GW). To establish a comprehensive assessment of overall GHG

E-mail address: [email protected]. https://doi.org/10.1016/j.fcr.2018.07.015 Received 5 February 2018; Received in revised form 26 July 2018; Accepted 30 July 2018 0378-4290/ © 2018 Elsevier B.V. All rights reserved.

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

fertilization, fungicides and growth regulators which were not applied in the moderate intensity system. In the experiments, mean application levels for the fertilizers were 107.7 kg for N ha−1 for the moderate intensity system and 147.4 for the high intensity one, and 53.5 kg for P2O5 ha−1 and 88.4 kg for K2O ha−1 for both management intensities. In the high intensity management system average foliar fertilization was 6 kg ha−1. Average herbicide and insecticide use was 1.3 and 0.2 kg ha−1, respectively. Fungicides were used only in the high intensity system at the mean amount of 1.7 kg ha−1. The average level of chemical seed protection was 0.3 kg ha−1 and growth regulators were used only in the high intensity system and were applied on average at the level of 0.8 kg ha−1. Winter wheat was sown at the end of September or at the beginning of October and harvested in August. Fertilizers and crop protection were applied following standard recommendations for winter wheat production.

form of fertilizers, pests, weeds, and diseases are effectively managed, cultivars are carefully selected, crop growth is determined mostly by plant available water during the growing season, soil properties influencing the available-water holding capacity, and by solar radiation and temperature (Edreira et al., 2017). This is the case in intensive crop production where high wheat yields (9.0 t ha−1) in Polish agricultural conditions are achieved in high-yielding treatments and are comparable to yield in areas such as Germany (9.5 t ha−1), which employs the world’s most intensive agronomic management. Consequently, on the one hand the process of winter wheat grain production contributes to GW and on the other, wheat growth is vulnerable to negative impacts of climate change such as increasing temperatures, more extreme weather and more variable rainfall patterns (Pachauri et al., 2014; Mase et al., 2017; Mäkinen et al., 2018). A number of studies have been conducted using LCA to calculate the yield-scaled GWP of winter wheat production (Hillier et al., 2009; Chen et al., 2014; Wójcik-Gront and Bloch-Michalik, 2016). This study, however, focuses on drivers of variation in yield-scaled GWP in intensive wheat production based on supplementary variables to those included in GWP calculations (Wójcik-Gront and Bloch-Michalik, 2016). The variables studied relate to weather, soil, winter wheat cultivar and crop management and come from long term experiments across the whole of Poland. Identifying those variables benefits sustainable wheat cropping systems by providing management practices that maximize crop productivity, while simultaneously minimizing negative environmental impacts (Grassini and Cassman, 2012). Of the major cereal crops, wheat accounts for one of the largest global consumptions of N fertilizer (Heffer, 2013). Higher N inputs offer the potential for increases in yield. On the other hand, nitrogen fertilizer application has been recognized as the key factor for increasing N2O emissions from the agricultural sector. Wójcik-Gront and BlochMichalik (2016) report that the largest impact on GWP in Polish cereal production came from applying N fertilizer (avg. 68%) and when N fertilizer use was increased by 40 kg ha−1 that led to higher values of yield-scaled GWP. This might suggest inefficient N use. On that account, the analysis presented in this study investigates if nitrogen was indeed an important driver of yield improvement in comparison to other environment and management related variables. Both quantitative and qualitative variables were used to study the variability of yield-scaled GWP and yield. In the whole analysis, 6 qualitative and 14 quantitative variables were used. To analyze such a mixture of variables and the large, unbalanced amount of observation units classic statistical methods might fail to find meaningful agricultural patterns. Thus, to select the most influential variables, the classification and regression tree (CART) and random forest (RF) (Breiman et al., 1984; Tulbure et al., 2012) methods were used. As these methods are relatively new in this research field, this manuscript evaluates the potential for CART and RF to inform agronomic management decisions along with determination if environmental conditions are indeed so important in intensive winter wheat production.

2.2. Data used for characterizing the environment, genes and management in winter wheat yield-scaled GWP and yield analysis For each yield data point, weather, environmental, genetic and management variables were described (Table1). To describe the major drivers of variation in the yield-scaled GWP of winter wheat production in Poland the following variables were used: - agro-region of Poland (Studnicki et al., 2018), from 1 to 6 (Fig. 1); These are regions with similar relative performances of cultivars. There are differences between regions in agronomic conditions and presumably in agro-technical culture differentiating the yield-scaled GWP. In the cultivar recommendation, cultivar performance is evaluated in a specific trial location. However, in many cases the trial location does not overlap with the specific farmland and this information might not be useful for wheat producers. This is why the recommendation is given for the whole agro-region where locations with similar environmental conditions are clustered (Tapley et al., 2013). - climatic water balance – KBW for April and May, KBW for May and June and KBW for June-July during each growing season; This is water available to a plant during its crucial development phases, expressed as a difference between precipitation and evaporation, dependent on temperature, humidity, wind and solar radiation (Doroszewski et al., 2012; The Institute of Soil Science and Plant Cultivation – IUNG, 2018); this was partitioned into 16 bins from sufficient water availability > −50 mm, -50- −59 mm, …, to extreme water deficiency -190- −199 mm. - genetic variables: cultivar utilization type (based on grain protein content, amylolytic and proteolytic enzymes activity or grain density), country of origin of cultivar, cultivar frost resistance; The cultivars tested in the Polish Post-Registration Variety Testing System are selected based on the assumption that they give high yield in Polish environmental conditions. The Polish climate is moderate and varies from maritime to continental (Graczyk and Kundzewicz, 2016) and as such is similar to that in other CentralEastern European countries. Hence, the yield is highly differentiated regarding cultivar but it is over-exaggerated when an input variable has that many levels (in this case a genotype with 127 cultivars). Then, CART can be biased toward choosing this variable (Hill and Lewicki, 2006). Thus, instead of using cultivars, their characteristics were used. - soil quality (valuation classes according to the soil quality evaluation system in Poland compatible with regulations of the Council of Ministers – CM, 2012); In Poland, the class reflects the agricultural value of soils, the lower the class the more fertile the soils. The main criteria for classifying soil into a given valuation class are: grain-size composition, thickness of the humus horizon and humus contents, hydrological properties of the soil, reaction and calcium carbonate contents, structure and usage of the area.

2. Materials and methods 2.1. Yield experimental data Yield data were gathered from the Polish Post-Registration Variety Testing System (PRVTS, 2018), which evaluates the yield and other related traits of newly released cultivars in multi-environmental trials. Winter wheat yield data from 17,870 observations were used. This was associated with 75 field trial locations across the whole of Poland (Fig. 1), 7 growing seasons (from 2009/2010 to 2015/2016), 127 modern, high yielding winter wheat cultivars and two agronomical intensities. The data used in this study contained observations from moderate input intensity with mineral fertilization, seed preparation, and herbicide and insecticide use and from high input intensity with additional (40 kg ha−1 yr−1) nitrogen fertilization, use of foliar 20

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

Fig. 1. Points show 75 locations of the experiments which took place between 2009/2010 and 2015/2016 in Poland. Lines display the 16 Polish voivodeships and colors indicate the 6 regions used in the multivariate analysis. Grid straight lines represent the geographic coordinate system for Poland.

measure of the average amount of a specific GHG discharged into the atmosphere by a given source relative to units of activity (UNFCCC, 2018). EF for each agricultural chemical and fuel use is provided in Table A.2. Then, the calculated GWP was divided by the experimental yield to calculate the yield-scaled GWP. The value of a gas emission is expressed as carbon dioxide equivalents (CO2 eq.). More recent values of the GWP 100 with inclusion of climate-carbon feedbacks CO2 equivalents were used to calculate the yield-scaled GWP, i.e. the GWP of 1 kg of CH4 and 1 kg of N2O is the equivalent to the GWPs of 34 kg and 289 kg of CO2, respectively (IPCC, 2013). However, a few emission factors employed from other studies used old GWP values and those remained unchanged. The differences between the old GWPs of CH4 and N2O and the newest ones are several percent and should not influence the outcome of this study.

- date of sowing; - the number of days from sowing to harvesting; - prior crop. The variables, their frequencies and types (quantitative or qualitative) are presented in Table 1. Additionally, the amount of seeds, the rate of nitrogen, phosphorus, potassium and foliar fertilization, and the amount of fungicides, herbicides, insecticides and growth regulators per hectare were used to describe the winter wheat yield variability (Fig. 2). 2.3. LCA approach and assumptions while calculating the yield-scaled GWP A detailed description of the methodology used to calculate the yield-scaled GWP is provided in Appendix A based on the work of Wójcik-Gront and Bloch-Michalik (2016). The total GWP of winter wheat production is the sum of GHG emissions from fuel use in field operations (Table A.1.), seed preparation, production and application of mineral fertilizers, non-organic pesticides and growth regulators, and agricultural machinery usage. This is expressed as the amount of CO2 equivalents. In this study, the total emission from one hectare of the crop cultivation was calculated as the sum of emissions from each experimental input, i.e. the application rate multiplied by the appropriate emission factor (EF). EF is a

2.4. Statistical analysis The impact of cultivar and environment on yield-scaled GWP and additionally, management practices on yield in winter wheat production was determined using decision tree algorithm based methods: CART and random forest. CART is a non-parametric data mining method which predicts the value of a target variable by learning decision rules inferred from the data. When the target variable is quantitative, this is a non-parametric regression method that recursively 21

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

of yield-scaled GWP and yield of winter wheat. However, due to the main disadvantage of a simple tree model, its instability resulting from small changes in the data, the analysis in this study was compared with the output from Random Forest (RF) method conducted on the same dataset. This method consists of an ensemble of unpruned regression trees which are created by bootstrapping, i.e. sampling with replacement from the original dataset and creating bootstrap samples with the same number of observation units as in the original dataset. For each bootstrap sample, a regression tree is constructed. Then, the process is repeated. Finally, a prediction is made based on the weighted average of the ensemble of trees (Breiman, 1996). RF is a very good method for prediction accuracy. However, RF does not produce one tree as an output, so one CART model is much easier to interpret. The method has also been successfully used in agronomy (Tittonel et al., 2008; Tulbure et al., 2012; Fukuda et al., 2013; Everingham et al., 2016). The CART and RF models for the analysis of variables influencing yield-scaled GWP from life cycle assessment of winter wheat production were created with STATISTICA ver. 12 (StatSoft Inc., 2014). The total number of observations used in the models was 17,870. To avoid overfitting, when the CART model overreacts to random variation in the data lowering its predictive performance, the classic strategy of “pruning” overgrowth trees was used, i.e. the 10-fold cross validation method. Other parameters set in the CART analysis were: the stop rule – trimming to variance, stopping parameters – the minimum number of observations in a shared or final node was 1787 (10% of all observations), the maximum number of shared nodes was 1000, the starting value of the random number generator – 1, the rule of one standard error for cross validation procedure. Another output of the CART procedure is the importance ranking of the input variables in terms of their impact on the target variable variability. To judge the relative importance of each input variable its improvement scores are calculated at each split. This is performed for the actual split variables and for variables that provide competing splits and that may never appear in the tree. The variable significance is expressed on a scale from 0 to 1. A variable can be significant, even though it might not have been used for any split in the final CART tree. When using RF, the randomly selected, subsample proportion – with replacement – from an original set of observation units was set at 50%. The proportion of randomly chosen observations that will serve as a test sample for the created model was determined as 30%. In the RF model, each tree can be grown with a randomized set of input variables selected at each split while building a single tree. The number must be smaller than the total number of variables in the dataset. In this study, this number was log2M+1, M – the number of input variables to minimize the correlation between the RF trees. In the case of yieldscaled GWP M = 7, so log2M+1 = 3 and for yield M = 16, so log2M +1 = 5. The RF method produces an ensemble of regression trees; thus, its results lack interpretability. The number of trees created was 500 to stabilize the errors. The stopping parameters in the RF analysis were as follows: minimum number of cases 446, minimum number of levels 10, minimum number in child node 5, and maximum number of nodes 100. The stopping condition number of cycles to calculate mean error was – 10 and the percentage decrease in training error – 5. This states that if the training error does not improve by at least 5% over 10 cycles, then training should stop. Before performing CART, the descriptive analysis of target variable – yield was conducted to distinguish outliers. During the CART analysis outliers can be grouped into a separate subset but they might increase the relative error of the model so they should be identified and, if possible removed a priori. The reasoning behind outliers might vary, so CART was performed on data without outliers to enable more flexibility in translation to other similar studies. From the dataset, 4.9% of all observations were removed based on identified outliers. An observation unit was considered to be an outlier based on the 25th and 75th percentile, i.e. if the value was greater than the 75th percentile + 1.5•(the 75th percentile – the 25th percentile) or lower than the 25th percentile

Table 1 The qualitative and quantitative variables used in the analysis of variability of yield-scaled GWP and yield in winter wheat production. Variable

Values

Frequency

Variable type (qualitative or quantitative)

region

1 2 3 4 5 6 partitioned into 16 bins from wet to dry conditions > −50 mm to -190−199 mm I – best quality II – very good quality IIIa, IIIb – good quality, IVa, IVb – medium quality soils cereal legume rapeseed root crop 20 September – 11 October, 13, 17-20 October 281-335

2887 1979 2643 4145 3637 2579

qualitative

superior variety good quality variety bread variety non-baking and forage cookie variety spelt Australia Austria Belgium Czech Republic Denmark France Germany Lithuania Poland Sweden Switzerland UK 1-9 scale (1- low, 9 high resistance to frost)

589 9450 6970 825 32 4 4 375 35 4 14 2296 8499 16 6442 74 66 45

KBW April-May KBW May-June KBW June-July

soil valuation class

prior crop

sowing date

nr of days from sowing to harvesting cultivar type

cultivars country of origin

cultivar frost resistance

quantitative

687 3553 10453 3177

qualitative

2893 5977 8918 82

qualitative

quantitative

quantitative

qualitative

qualitative

quantitative

partitions the multidimensional space spanned by input variables into subsets including observations with target values as homogeneous as possible according to variables producing the highest impurity (variance) reduction (Strobl et al., 2009). In this study, the regression tree partitions the observations of all winter wheat genotypes and growth environments dichotomously splitting them into subsets to minimize the variability of a target variable (the yield-scaled GWP or yield) within the subset while maximizing it between the subsets. The method works best on a considerable number of observations. It has many advantages compared to standard statistical techniques: it does not need data transformation because there is no statistical distribution assumption regarding input and target variables; it performs with unbalanced and missing data, works on quantitative and qualitative data and combinations of the two, and deals with non-linearity and highorder interactions. For quantitative variables, CART only makes splits and does not pick an arbitrary value to produce a better model fit. CART has so far been used in several agricultural studies (Ferraro et al., 2009; Krupnik et al., 2015). The flexibility of CART makes it suitable for the problem in this study of detecting quantitative or qualitative variables

22

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

Fig. 2. Histograms for all input variables (in kg) and the target variables – yield (t) and yield-scaled emission (kg CO2 eq. kg−1 grain). Next to histograms are corresponding box-plots with medians in red. The first and third quartiles are outlined by the boxes. Cases below the lower fence or above the upper fence are considered outliers (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

– 1.5•(the 75th percentile – the 25th percentile). This led to 17,870 observations being taken into the CART classification tree method and further analysis.

variables used in the analysis: prior crop, soil valuation class, location’s agro-region and cultivar characteristics. There were also quantitative ones: KBW, sowing date, number of days from sowing to harvest (Table 1). Before performing the analysis, the Spearman correlations (as the variable distribution is not always normal) between variables were verified (Table 2). In most cases, there was a weak correlation. A positive, strong correlation is observed for variables of inputs used in more intensive agronomical management system (e.g., fungicides and growth regulators, 0.87). Also, potassium and phosphorus were characterized with strong positive correlations. In locations where higher doses of one fertilizer were used, the other fertilizer and pesticide rates were probably also higher. Another significant correlation coefficient was observed between nitrogen application rate and the yield-scaled GWP (0.44).

3. Results 3.1. Descriptive statistics of the analyzed data The histograms and box-plots of the quantitative variables used in the analysis of variability in target variables, yield-scaled GWP and yield, in winter wheat production are presented in Fig. 2. These include the amount of seeds, nitrogen, phosphorus, potassium, and foliar fertilization, and herbicides, insecticides, fungicides and growth regulator used per hectare of the field. The histograms for these variables are highly skewed with high numbers for zero values. Distribution for N fertilizer is very wide. This is caused by using additional (around 40 kg ha−1) amounts at high input intensity; thus, it might be two distributions combined with peaks separated by 40 kg. In many cases, the data were collected across more than one season in the same location, so previous treatments might influence subsequent amounts of fertilizers applied. Moreover, there were additional qualitative

3.2. Yield-scaled GWP A single CART with cross validation performed on 17,870 observations having yield-scaled GWP as the target variable identified soil quality as the main variable explaining yield-scaled GWP variability (Fig. 3). Unsurprisingly, observations with higher average yield-scaled 23

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

Table 2 Correlation matrix (with Spearman coefficients) between quantitative variables used in the analysis of variability in yield-scaled GWP and yield in winter wheat production.

seeds nitrogen phosphorus potassium foliar fertilization herbicides insecticides fungicides growth regulator yield emission/yield

seeds

nitrogen

phosphorus

potassium

foliar fertilization

herbicides

insecticides

fungicides

growth regulator

yield

emission/yield

1.00

0.10 1.00

−0.09 −0.05 1.00

−0.01 0.08 0.65 1.00

−0.09 0.59 −0.01 0.00 1.00

−0.13 −0.10 0.04 0.04 0.04 1.00

−0.07 −0.04 −0.05 −0.03 0.02 −0.08 1.00

−0.02 0.68 −0.01 −0.01 0.75 0.01 0.04 1.00

−0.03 0.66 0.02 0.05 0.70 −0.03 0.03 0.87 1.00

−0.19 0.27 0.07 0.15 0.23 0.00 −0.02 0.34 0.36 1.00

0.24 0.44 0.04 0.07 0.23 −0.05 −0.01 0.17 0.15 −0.69 1.00

GWP 0.33 kg CO2 eq.kg−1 yield (group 2) came from less fertile soils (groups IVa, IVb). These were further split by KBW calculated for June and July. The same variable split observations with lower average yield-scaled GWP 0.28 kg CO2 eq.kg−1 yield (group 3) obtained on better soils for wheat cultivations. However, this time the split is made for a different value of the variable. This results from the way the yieldscaled GWP is calculated. The GWP obtained for one hectare of winter wheat field is divided by its yield; thus, the higher the yield the lower yield-scaled GWP. Soils types I, II and IIIa, IIIb are good quality soils allowing higher cropping. They are favorable for wheat growth, with sufficient or high nutrient reserves, and adequate water holding capacity, with structure suitable for root growth. Thus, these soils do not

need external nutrient inputs in the same rates as less fertile soils do. High values of KBW indicate high precipitation or lack of sun radiation and this leads to higher yield-scaled emissions (group 4 and group 8). In the further splits, other variables explained yield-scaled GWP variability. These were region, prior crop, length of growing period (nr days), sowing and KBW calculated for April and May and KBW for May and June. The cultivar describing variables were unimportant while creating the splits. Following the various splits that resulted in the final subsets (nodes), the interpretation of the way a subset was created might be derived. The data were classified into 17 final subsets. CART can identify the conditions that give the lowest yield-scaled GWP. The lowest yield-scaled GWP (0.20 kg CO2 eq.kg−1 yield) was obtained for

Fig. 3. CART for yield-scaled GWP (kg CO2 eq.kg−1 grain). In white are marked splits in grey final subsets. 24

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

many experiments where winter wheat growing season was the longest. In regions 3 and 5, the average use of fungicides was the highest and in region 5 the use of growth regulators was the highest. The next important variable was prior crop. In CART, lower yieldscaled GWP was observed for prior crops other than cereals. Previous studies using results from a long-term experiment showed that the yield of wheat is lower in a monoculture than in a crop rotation system (Berzsenyi et al., 2000). The length of the growing season was identified as a moderately important variable. Lower yield-scaled GWP was rather obtained for shorter growing seasons. Sometimes harvest can be delayed due to high precipitation, which might be the reason for yield loss. Cultivar describing variables were not important in the tree construction. Wheat quality class should not influence the yield-scaled GWP variability which is an independent trait and this was confirmed by the CART analysis results (importance – 0.08). When the cultivar frost resistance was considered, it also turned out to be unimportant while creating the trees (0.10). Eventually, when the cultivar country of origin was considered its importance was evaluated as 0.26. This suggests that cultivar is indeed not important in influencing the yieldscaled GWP variability or that there are other characteristics differentiating target variables which are not available in this study like, for example drought or water-logging resistance.

Fig. 4. Variable importance (from 0 – unimportant to 1 – most important) in CART and Random Forest models for yield-scaled emission.

the following variable combinations: good quality soils (soil types I, II, IIIa, IIIb), KBW calculated for June and July lower than −99 mm (drier conditions), KBW for May and June lower than −89 mm, and number of growing days lower than 299. The Pearson’s correlation coefficient calculated between predicted and experimental yield-scaled GWP for the tree with cross validation pruning was 0.57. In the case of RF, the Pearson’s correlation coefficient calculated between predicted and experimental yield-scaled GWP was 0.66. RF provides a hierarchy of variables regarding their impact on yield-scaled GWP based on an ensemble of regression threes, which in this analysis was different than in CART (Fig. 4). CART analysis identified KBW calculated for months June-July as the most influential variable (importance value 1.00), while RF selected soil class as having the highest importance value. Other variables with importance higher than 0.5 were as follows: in CART – soil class (0.98), region (0.78), and number of days from sowing to harvesting (0.61) and prior crop (0.60); in RF – KBW June-July (0.84) and region (0.54). The importance of the remaining variables was estimated to be between 0.08 to 0.49 in both models (Fig. 4). The results suggest a relative homogeneity of the environmental conditions within an agro-region and variability of yieldscaled GWP between the agro-regions. This is similar to the results of a Polish winter wheat grain yield analysis based on linear mixed models conducted by Derejko et al. (2016). In the analysis, the most important variables turned out to be water availability and soil class, which are environmental conditions of winter wheat production. During the final months of plant growth drier weather is more favorable for obtaining higher yields. This is related to yield variation which depends heavily on soil quality, precipitation and water evaporation, although not in a linear manner. In the CART tree, the variable soil class was shown to give the most reduction in variability. However, in the overall variable importance calculations the variable KBW June-July was the most important. The next important variable was agro-region. The first split in CART was due to soil class. On soils with lower agricultural productivity (IVa and IVb) in regions 2, 3, 6 the yield-scaled GWP was lower than in regions 1, 4, 5. In region 2, average nitrogen use for soils IVa and IVb across the analyzed years was medium (123.4 kg ha−1) while average yield was quite high (8.4 t ha−1). In region 3, medium N use (132.2 kg ha−1) gave high yield (9.0 t ha−1) and in region 6 yield was low (7.3 t ha−1), but N rates were also relatively low N (105.5 kg ha−1). In regions 1, 4, 5, one of the following two situations existed: the N dose was very high (region 1 average N rate 160 kg ha−1) or yield was low (region 4 – yield 7.4 t ha−1). The differences between Polish agro-regions are caused by many related variables including soil fertility, weather, the length of growing season and agro-technical intensity. As far as the soil class is concerned there are some differences between Polish agro-regions regarding its distribution. There were no experiments on soil class I in regions 2, 3 and 5 and there were no experiments on soil class IVb in regions 1 and 3. In regions 1 and 3, the average yield across all years and soil classes was the highest. In region 1, there were

3.3. Yield To examine winter wheat yield variability, besides the variables used in the yield-scaled GWP analysis additional variables were used: amount of nitrogen, phosphorus, potassium and foliar fertilizers, seeds, insecticides, herbicides, fungicides and growth regulators. In the calculations regarding the amount of emissions from one hectare of winter wheat production, pre-sowing treatment and lime were included in the analysis. These two variables were not used in this analysis though, as the pre-sowing treatment was dependent on the amount of seeds used on the field, so there would be a correlation, and lime is applied every four years and affects the wheat yield especially on acidic soils but there were no detailed data on the date of application. In the single CART analysis, the data were classified into 16 final subsets. The CART analysis revealed growth regulator as the variable where the split gave the most reduction in yield variability (Fig. 5). Observations where the amount of growth regulator applied to one hectare of field was lower than 0.22 (group 2) had a lower average yield (equal to 8.48 t), and observations with a higher rate had an average yield of 9.76 t (group 3). Both groups were further split by KBW calculated for the months JuneJuly in such a way that in each group observations for KBW lower than −99 mm (drier conditions) had a higher yield (group 5, yield = 9.19 t and group 21, yield = 10.42 t). In the further splits, other variables explained yield variability. These were agro-region, soil class, KBW May-June, seeds, and fungicides. The highest yield (11.17 t) was obtained for the following variable combinations: growth regulator greater than 0.224, KBW June–July lower than −99 mm, soil class I, II or IIIa, IIIb and KBW May-June lower than −59 mm. The Pearson’s correlation coefficient calculated between the yield predicted by CART and the experimental one was 0.59. In the case of RF, the correlation coefficient was 0.75 and the hierarchy of variables regarding their impact on yield was slightly different (Fig. 6). Both regression tree based methods identified environmental conditions as the most influential variables (importance value close to 1.00). Variables with importance higher than 0.5 were as follows: in CART – region (1.00), growth regulator (0.90), soil class (0.89), KBW June-July (0.88), fungicides (0.78), KBW May-June (0.67), and seeds (0.61); in RF – KBW June-July (1.00), growth regulator (0.87), fungicides (0.83), region (0.74), soil class (0.72), nitrogen (0.61) and KBW May-June (0.61). The importance of the remaining variables was estimated to be between 0.08 to 0.48 in both models, i.e. from unimportant to moderately important (Fig. 6). 25

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

Fig. 5. CART for yield (t).

analysis. There are yields for both intensity levels: lower and higher ones with fungicides, growth regulators, foliar fertilization and an additional 40 kg of nitrogen applied in the production. A very weak yield dependence on nitrogen can be observed (the points form almost horizontal lines for most years). The yield increase between moderate and high agricultural management intensity is mostly caused by using growth regulators and fungicides. The results also indicate no significant influence of insecticides or herbicides, although the literature reports that yields are greatest when these are applied (Oerke, 2006). This is mainly because their use was similar in each of the two studied management intensities. In this study, seed density was negatively correlated with the wheat yield and this was moderately important in explaining yield variability. Fig. 6. Variable importance (from 0 - unimportant to 1 – most important) in CART and Random Forest models for yield.

4. Discussion A very important outcome of this analysis is that N application rate used in excess is not a driver of yield improvement. Yield dependence on nitrogen application rates is non-monotonic (Hawkesford, 2014; Zhou and Butterbach-Bahl, 2014). If nitrogen fertilizer is applied in excess of uptake by crops, with attempts to increase yields simply by increasing inputs, this can contribute to agricultural inefficiency and environmental damage, for example atmospheric pollution and groundwater contamination or surface water eutrophication which can threaten human health (Basso and Ritchie, 2005; Li et al., 2007). Based on the results from the moderate and high intensity agronomic management, there might be a potential to improve site-specific N management aimed at increasing yield and reducing soil emissions. This could be achieved by further research with the use of a gradient of N application rates in each specific environmental condition to specify N rates according to local climate and soil (Cui et al., 2008). Overall, the results support previous research indicating that increasing N fertilizer application may result in increased nitrous oxide (N2O) emissions from

While creating the optimal CART tree for yield, the variable growth regulator was shown as the one splitting the dataset in such a way that most variability reduction was obtained. In the first CART split, higher values of yield were observed for growth regulators used in amounts higher than 0.2 kg ha−1. As in CART, only one variable is used in the model at any given split, correlated variables may perform in almost the same way. This variable is highly correlated with another main variable in the analysis, i.e. fungicides (correlation coefficient 0.87). Thus, the observed grain yield increase is mostly due to application of additional protection, i.e. growth regulators and fungicides, which were shown to be very important in the regression tree based analysis applied in this study. Berry et al. (2008) stated that across elite wheat cultivars, fungicides increased yields by an average of 1.78 t ha−1 and hence reduced GHG emissions per ton. In Fig. 7, there are data on winter wheat yield dependence on nitrogen application rates in each of the 7 years taken into account in the 26

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

Fig. 7. Yield (kg) dependence on nitrogen application rates (kg) in each year of experiments included in the analysis. Dots represent data for lower intensity management and squares data for high management intensity with fungicides, growth regulators, foliar fertilization and an additional 40 kg of nitrogen applied in winter wheat production. Lines represent linear regressions and determination coefficients for fits. The highest N rates (above 200 kg ha−1 in 2011, 2014 and 2016) were observed in 4 locations in west Poland. In general, regions 1 and 3 had higher average N application rates.

a lower impact from soil class might be observed in the yield-scaled GWP regression tree analysis. The study by Bouwman et al. (2002) also summarizes that emissions from poorly drained soils exceed emissions from well drained ones. When oxygen availability is restricted, it enhances denitrification processes. That, combined with lower yield due to high soil moisture, contributes to high yield-scaled GWP. The KBW impact in creating the CART and RF trees is already very strong, so increased N2O emissions due to high soil water content might not change the tree structure but rather only emphasize it. While calculating yield-scaled GWP for winter wheat production, uncertainty over the result was estimated based on assumed uncertainty of the emission factors and application rates used. The biggest input to the estimated total yield-scaled GWP uncertainty came from the emission factor for using N fertilizers (Wójcik-Gront and Bloch-Michalik, 2016). This is due to the nature of the phenomena. Thus, in the sensitivity analysis the input from nitrogen was lowered and increased in comparison to the actual input used in the calculations and analysis by 9.5 kg CO2 eq.kg−1 N applied. When calculating yield-scaled GWP, the denominator (yield) was not changed but the numerator was lower or higher. However, in both models the environment was still the most important variable for lower and higher yield-scaled GWP, thus indicating a strong dependence on winter wheat yield performance in specific environmental conditions. The most important variables contributing to winter wheat yield variability were weather and soil growing conditions. The same observation was made by Rozbicki et al. (2015) based on an analysis using linear mixed models. Winter wheat production is subject to year-to-year fluctuations influenced by seasonal and inter-annual climatic conditions in many regions of the world (Tapley et al., 2013; Derejko et al., 2016). Variation in grain yield is related to grain number per area unit and therefore depends on the availability of resources between the onset of stem elongation and shortly after flowering during grain filling (in climates similar to that in Poland, this would be May and June). In this study, KBW May-June was an important variable in explaining yield variability and this relates to sufficient water availability during that period. Soil class was also a very important variable; where higher quality classes contributed to higher yields. Wheat production requires soils with high agronomical productivity. In intensive wheat

soil and higher yield-scaled GWP (Galloway et al., 2003; Tonitto et al., 2006; Gruver and Galloway, 2008; Bouwman et al., 2013; Chen et al., 2014; Wang and Dalal, 2015). The LCA-based GWP used does not take into account temporal variation in N2O emissions from soils. That, however, needs long-term, year-round monitoring to obtain representative estimates of nitrous oxide emissions from agricultural systems (Zhou et al., 2017), which was not available at the time when these experiments took place. Zhou et al. (2017) noticed that in a rice-wheat rotation system substantial N2O emissions only occurred following events of N fertilization and drainage practices. Huang et al. (2017) proved that optimized N application, where the fertilizer rates were estimated based on crop nitrogen demand, reduced cumulative N2O emissions by 18–37%, which in turn decreased yield-scaled N2O emissions by 38–42%. A promising way to reduce fertilizer-induced nitrous oxide emissions is to use nitrification inhibitors (Liu et al., 2013) which have been shown to have no negative effect on grain yield in rain-fed areas (Villar and Guillaumes, 2010). Soil N2O emissions are known to be highly dependent on mineral nitrogen availability and soil moisture and are highly sensitive to the inter-annual variability of climatic conditions. Bouwman et al. (2002) report that organic soils show higher emissions compared to mineral ones. Based on 139 field studies with 896 measurements, N2O emissions increased with increasing soil organic content. Additionally, there are also emissions due to changes in soil carbon. A large quantity of carbon is stored as soil organic carbon (SOC). However, there is still no detailed information on the size of SOC stocks, their spatial distribution, and the CO2 emissions from soils due to changes in land use and land cover (Scharlemann et al., 2014). On mineral soils, conversion of uncultivated land to cropland can result in the loss of 20–40% of the original soil carbon stocks (Davidson and Ackerman, 1993). Better quality soils have higher amounts of SOC in comparison to lower quality soils. Thus, GHG emissions due to conversion of uncultivated land to cropland might be higher in the case of good quality soils (in this case, soil classes I and II). Moreover, soils with high organic carbon content are prone to high N2O emissions. This could influence yield-scaled GWP, which would increase more in cases of better quality soils in comparison to low quality soils when taking into account additional emissions from soils with high SOC. As a result, 27

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

significant impact on GHG emissions in winter wheat production. However, this variable’s importance in creating regression trees regarding target variable winter wheat grain yield was moderate (importance around 0.55). Intensively managed croplands are the most important source of leached nitrate. An imbalance between applied nitrogen inputs and crop demand indicates the potential to further improve environmental performance through better management of these inputs. This study demonstrated that it is possible to reduce the negative environmental impact of N2O emissions without compromising crop productivity. More important variables of yield variability were use of fungicides and growth regulators. When planning to obtain the highest possible yield one should consider using growth regulators and fungicides, an average dose of nitrogen (around 110 kg ha−1), a prior crop other than cereals and cultivars recommended by PRVTS. Combining assessment of both GHG emissions and wheat yield is necessary in optimizing cropping practices. This study is applicable to defining goals for sustainable intensification of cropping practices and to defining policy trade-offs between GHG emission mitigation and food security. This study also shows that CART and Random Forest are capable of detecting the actual variables influencing the target variable, especially when a large data set is used.

production, the best available soils are chosen. Environmental effects have been shown to be responsible for approximately 80% or more of total experimental variation in wheat grain yield (Mohammadi et al., 2010; Tapley et al., 2013; Slafer et al., 2014; Derejko et al., 2016). In many cases, water has an important influence on yield stability. In Poland, most winter wheat production is rain-fed. Surprisingly, a more yield restricting factor is water overflow in the final winter wheat growth phase, in June and July. In the yield-scaled GWP tree, soils with lower agricultural fertility (IVa and IVb) with KBW June–July ≥ −50 give the highest yield-scaled GWP. These soils are strongly influenced by high precipitation. They are characterized by low soil permeability which can negatively influence plant conditions. Weather conditions in Poland are characterized by high changeability. Sometimes during harvest unfavorable days with continuous rainy, wet, humid weather occur. Optimizing irrigation or drainage might increase yield as well as lead to lowering yield-scaled GWP (Lv et al., 2011). Too much water during final development stages of winter wheat growth can cause lower yield. Likewise, artificially drained fields might achieve higher yields. To some extent, the variability of yield in winter wheat production in the climatic and soil conditions of Poland can be neutralized by appropriate agronomic management intensity and selection of cultivars. Thus, yield-scaled GWP can be maintained at the lowest possible level. An interesting finding is that fungicides and growth regulators are very important variables in explaining winter wheat yield variability. This finding is in agreement with the observation that when there are high amounts of precipitation combined with low solar radiation causing water excess in the soil foliar diseases can attack plants, particularly in winter wheat that has not been treated with fungicide. Other undesirable effects of excess water available to plants are proliferation of pests, leakage of nutrients, inhibition of oxygen uptake by roots, and water-logging during harvest (Zampieri et al., 2017) which leads to yield loss. Thus, the yield is positively correlated with the amount of growth regulators used. Soil affected by excess water also favors a loss of nitrogen, lowering its availability to the plants (Huang et al., 1995). As the relationship between precipitation, fertilization and crop production is not linear, it is essential that this is further studied and understood to enable sustainable food production. Pesticides are not neutral to the environment. One of the harmful effects of pesticides is their acute lethal toxicity to bees (Huntzinger et al., 2008; Stanley et al., 2015). As has been proved, the use of fungicides is crucial in yield improvement, especially in cold, moist climates, but farmers have to choose the least environmentally-risky chemicals and use them in concentrations not exceeding legislatively determined safe levels. The situation encountered in the data set used in the presented CART analysis is unbalanced data. There are agro-regions used in the analysis with no observations from soil class I. When compared to other soil classes, this class has a lower number of examples. However, in regression trees the presence of balanced or unbalanced data sets is not critical, since as long as there are enough data the method performs reasonably well. In this study, the data set consisted of a considerable number (17,870) of observations. Moreover, 10-fold cross validation was incorporated into the analysis and the results were compared with the outcome of bootstrapping based method RF. Therefore, these two methods can be successfully used in a multivariate analysis.

References Basso, B., Ritchie, J.T., 2005. Impact of compost, manure and inorganic fertilizer on nitrate leaching and yield for a 6-year maize-alfalfa rotation in Michigan. Agric. Ecosyst. Environ. Appl. Soil Ecol. 108, 329–341. Berry, P.M., Kindred, D.R., Paveley, N.D., 2008. Quantifying the effects of fungicides and disease resistance on greenhouse gas emissions associated with wheat production. Plant Pathol. 57, 1000–1008. Berzsenyi, Z., Győrffy, B., Lap, D., 2000. Effect of crop rotation and fertilisation on maize and wheat yields and yield stability in a long-term experiment. Eur. J. Agron. 13 (23), 225–244. Bouwman, A.F., Boumans, L.J.M., Batjes, N.H., 2002. Emissions of N2O and NO from fertilized fields. Summary of available measurement data. Global Biogeoch. Cy. 16 (4), 1058. Bouwman, A.F., Beusen, A.H.W., Griffioen, J., Van Groenigen, J.W., Hefting, M.M., Oenema, O., Van Puijenbroek, P.J.T.M., Seitzinger, S., Slomp, C.P., Stehfest, E., 2013. Global trends and uncertainties in terrestrial denitrification and N2O emissions. Philos. Trans. R. Soc. B -Biol. Sci. 368. Breiman, L., 1996. Bagging predictors. Mach. Learn. 26 (2), 123–140. Breiman, L., Friedman, J.H., Olshen, R.A., Stone, C.G., 1984. Classification and Regression Trees. Wadsworth International Group, Belmont, California, USA. Charles, R., Jolliet, O., Gaillard, G., Pellet, D., 2006. Environmental analysis of intensity level in wheat crop production using life cycle assessment. Agric. Ecosyst. Environ. 113, 216–225. Chen, X., Cui, Z., Fan, M., Vitousek, P., Zhao, M., Ma, W., Wang, Z., Zhang, W., Yan, X., Yang, J., 2014. Producing more grain with lower environmental costs. Nature 514, 486–489. CM, 2012. Regulation of the Council of Ministers on Soil Quality Classes. Dz.U.2012.0.1246 – Rozporządzenie Rady Ministrów Z Dnia 12 Września 2012 r. w Sprawie Gleboznawczej Klasyfikacji Gruntów. Consoli, F., Allen, D., Boustead, I., Fava, J., Franklin, W., Jensen, A.A., de Oude, N., Parrish, R., Perriman, R., Postlethwaite, D., Quay, B., Séguin, J., Vignon, B. (Eds.), 1993. Guidelines for Life-Cycle Assessment: A ‘Code of Practice’. Society of Environmental Toxicology and Chemistry (SETAC), Brussels. Cui, Z., Zhang, F., Chen, X., Miao, Y., Li, J., Shi, L., Xu, J., Ye, Y., Liu, C., Yang, Z., Zhange, Q., Huang, S., Bao, D., 2008. On-farm evaluation of an in-season nitrogen management strategy based on soil Nmin test. Field Crop. Res. 105 (1-2), 48–55. Davidson, E.A., Ackerman, I.L., 1993. Changes in soil carbon inventories following cultivation of previously untilled soils. Biogeochemistry 20, 161–193. Derejko, A., Studnicki, M., Mądry, W., Gacek, E., 2016. A comparison of winter wheat cultivar rankings in groups of polish locations. Cereal Res. Comm. 44 (4), 628–638. Doroszewski, A., Jadaczyszyn, J., Kozyra, J., Pudełko, R., Stuczyński, T., Mizak, K., Łopatka, A., Koza, P., Górski, T., Wróblewska, E., 2012. Podstawy Systemu Monitoringu Suszy Rolniczej. (eng. Fundamentals of the agricultural drought monitoring system) Woda Środ. Obsz. Wiej 12 (38), 77–91 2. Edreira, J.I.R., Mourtzinis, S., Conley, S.P., Roth, A.C., Ciampitti, I.A., Licht, M.A., Kandel, H., Kyveryga, P.M., Lindsey, L.E., Mueller, D.S., Naeve, S.L., Nafziger, E., Specht, J.E., Stanley, J., Staton, M.J., Grassini, P., 2017. Assessing causes of yield gaps in agricultural areas with diversity in climate and soils. Agr. Forest Meteorol. 247, 170–180. Everingham, Y., Sexton, J., Skocaj, D., Inman-Bamberet, G., 2016. Accurate prediction of sugarcane yield using a random forest algorithm. Agron. Sustain. Dev. 36, 27. FAOSTAT, 2018. FAOSTAT Crop Production Database. Available at: http:// www.fao.org/faostat/en/#data/QC (Accessed 26 April 2018). . Ferraro, D.O., Rivero, D.E., Ghersa, C.M., 2009. An analysis of the factors that influence sugarcane yield in Northern Argentina using classification and regression trees. Field Crops Res. 112 (2-3), 149–157.

5. Conclusions In this study, the regression tree based analysis revealed that soil quality and weather are the most influential variables in terms of the yield-scaled GWP variability and yield of winter wheat. This means that environmental favorable conditions in particular can contribute to high yields close to potential grain yield of winter wheat. Adequate precipitation distribution during intensive plant growth is key to obtaining high yields. Based on positive correlations between nitrogen rate and yield-scaled GWP (0.44), these fertilizer application rates had a 28

Field Crops Research 227 (2018) 19–29

E. Wójcik-Gront

S., Hegerl, G., Howden, M., Jiang, K., Jimenez Cisneroz, B., Kattsov, V., Lee, H., Mach, K.J., Marotzke, J., Mastrandrea, M.D., Meyer, L., Minx, J., Mulugetta, Y., O'Brien, K., Oppenheimer, M., Pereira, J.J., Pichs-Madruga, R., Plattner, G.K., Pörtner, H.O., Power, S.B., Preston, B., Ravindranath, N.H., Reisinger, A., Riahi, K., Rusticucci, M., Scholes, R., Seyboth, K., Sokona, Y., Stavins, R., Stocker, T.F., Tschakert, P., van Vuuren, D., van Ypserle, J.P., 2014. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. In: Pachauri, R., Meyer, L. (Eds.), IPCC, Geneva, Switzerland. PRVTS –– the Polish Post-Registration Variety Testing System, available at: http://www. coboru.pl/English/index_eng.aspx. Rozbicki, J., Ceglinska, A., Gozdowski, D., Jakubczak, M., Cacak-Pietrzak, G., Madry, W., Golba, J., Piechocinski, M., Sobczynski, G., Studnicki, M., Drzazga, T., 2015. Influence of the cultivar, environment and management on the grain yield and breadmaking quality in winter wheat. J. Cereal Sci. 61, 126–132. Scharlemann, J.P.W., Tanner, E.V.J., Hiederer, R., Kapos, V., 2014. Global soil carbon: understanding and managing the largest terrestrial carbon pool. Carbon Manag. 5 (1), 81–91. Slafer, G.A., Savin, R., Sadras, V.O., 2014. Coarse and fine regulation of wheat yield components in response to genotype and environment. Field Crops Res. 157, 71–83. Snyder, C.S., Bruulsema, T.W., Jensen, T.L., Fixen, P.E., 2009. Review of greenhouse gas emissions from crop production systems and fertilizer management effects. Agric. Ecosyst. Environ. 133, 247–266. Stanley, J., Sah, K., Jain, S.K., Bhatt, J.C., Sushil, S.N., 2015. Evaluation of pesticide toxicity at their field recommended doses to honeybees, Apis cerana and A-mellifera through laboratory, semi-field and field studies. Chemosphere 119, 668–674. StatSoft Inc, 2014. STATISTICA (Data Analysis Software System), Version 12. www. statsoft.com. Strobl, C., Malley, J., Tutz, G., 2009. An introduction to recursive partitioning: rationale, application, and characteristics of classification and regression trees, bagging, and random forests. Psychol. Methods 14, 323–348. Studnicki, M., Derejko, A., Wójcik-Gront, E., Kosma, M., 2018. Adaptation patterns of winter wheat cultivars in agro-ecological regions. Sci. Agric in press. Tapley, M., Ortiz, B.V., van Santen, E., Balkcom, K.S., Mask, P., Weaver, D.B., 2013. Location, seeding date, and variety interactions on winter wheat yield in Southeastern United States. Agron. J. 105, 509–518. Tittonel, P., Shepherd, K.D., Vanlauwe, B., Giller, K.E., 2008. Unravelling the effects of soil and crop management on maize productivity in smallholder agricultural systems of western Kenya-An application of classification and regression tree analysis. Agric. Ecosyst. Environ. 123, 137–150. Tonitto, C., David, M.B., Drinkwater, L.E., 2006. Replacing bare fallows with cover crops in fertilizer-intensive cropping systems: a meta-analysis of crop yield and N dynamics. Agric. Ecosyst. Environ. 112, 58–72. Tulbure, M.G., Wimberly, M.C., Boe, A., Owens, V.N., 2012. Climatic and genetic controls of yields of switchgrass, a model bioenergy species. Agric. Ecosyst. Environ. 146 (1), 121–129. UNFCCC - United Nations Framework Convention on Climate Change, available at: https://unfccc.int/process/transparency-and-reporting/greenhouse-gas-data/greenhouse-gas-data-unfccc/definitions. van Groenigen, J.W., Velthof, G.L., Oenema, O., van Groenigen, K.J., van Kessel, C., 2010. Towards an agronomic assessment of N2O emissions: a case study for arable crops. Eur. J. Soil Sci. 61, 903–913. Villar, J.M., Guillaumes, E., 2010. Use of nitrification inhibitor DMPP to improve nitrogen recovery in irrigated wheat on calcareous soil. Span. J. Agric. Res. 8, 1218–1230. Wang, W., Dalal, R.C., 2015. Nitrogen management is the key for low-emission wheat production in Australia: a life cycle perspective. Eur. J. Agron. 66, 78–82. Williams, A.G., Audsley, E., Sandars, D.L., 2010. Environmental burdens of producing bread wheat, oilseed rape and potatoes in England and Wales using simulation and system modelling. Int. J. Life Cycle Ass. 15, 855–868. Wójcik-Gront, E., Bloch-Michalik, M., 2016. Assessment of greenhouse gas emission from life cycle of basic cereals production in Poland. Zemdirbyste-Agric. 103, 259–266. Zampieri, M., Ceglar, A., Dentener, F., Toreti, A., 2017. Wheat yield loss attributable to heat waves, drought and water excess at the global, national and subnational scales. Environ. Res. Lett. 12, 064008. Zhou, M., Butterbach-Bahl, K., 2014. Assessment of nitrate leaching loss on a yield-scaled basis from maize and wheat cropping systems. Plant Soil 374, 977. Zhou, M., Zhu, B., Wang, X., Wang, Y., 2017. Long-term field measurements of annual methane and nitrous oxide emissions from a Chinese subtropical wheat-rice rotation system. Soil Biol. Biochem. 115, 21–34.

Fukuda, S., Spreer, W., Yasunaga, E., Yuge, K., Sardsud, V., Müller, J., 2013. Random Forests modelling for the estimation of mango (Mangifera indica L. cv. Chok Anan) fruit yields under different irrigation regimes. Agric. Water Manage. 116, 142–150. Galloway, J.N., Aber, J.D., Erisman, J.W., Seitzinger, S.P., Howarth, R.W., Cowling, E.B., Cosby, J., 2003. The nitrogen cascade. Bioscience 53, 341–356. Graczyk, D., Kundzewicz, Z.W., 2016. Changes of temperature-related agroclimatic indices in Poland. Theor. Appl. Climatol. 124 (1-2), 401–410. Grassini, P., Cassman, K.G., 2012. High-yield maize with large net energy yield and small global warming intensity. PNAS 109 (4), 1074–1079. Gruver, N., Galloway, J.N., 2008. An Earth-system perspective of the global nitrogen cycle. Nature 451 (7176), 293–296 GUS – Statistics Poland (Główny Urząd Statystyczny), available at: http://stat.gov.pl/en/. Hawkesford, M.J., 2014. Reducing the reliance on nitrogen fertilizer for wheat production. J. Cereal Sci. 59, 276–283. Heffer, P., 2013. Assessment of Fertilizer Use by Crop at the Global Level. International Fertilizer Industry Association (IFA), Paris, France. Hill, T., Lewicki, P., 2006. Statistics: Methods and Applications a Comprehensive Reference for Science, Industry, and Data Mining. StatSoft, Inc., USA. Hillier, J., Hawes, C., Squire, G., Hilton, A., Wale, S., Smith, P., 2009. The carbon footprints of food crop production. Int. J. Agric. Sustain. 7, 107–118. Huang, B., Johnson, J.W., NeSmith, D.S., Bridges, D.C., 1995. Nutrient accumulation and distribution of wheat genotypes in response to waterlogging and nutrient supply. Plant Soil 173 (1), 47–54. Huang, S., Lv, W.S., Bloszies, S., Shi, Q.H., Pan, X.H., Zeng, Y.J., 2017. Effects of fertilizer management practices on yield-scaled ammonia emissions from croplands in China: a meta-analysis. Field Crops Res. 192, 118–125. Huntzinger, C.I., James, R.R., Bosch, J., Kemp, W.P., 2008. Fungicide tests on adult alfalfa leafcutting bees (Hymenoptera: Megachilidae). J. Econ. Entomol. 101 (4), 1088–1094. Myhre, G., Shindell, D.,Bréon, F.‐M., Collins, W. Fuglestvedt, J., Huang, J., Koch, D., Lamarque, J.‐F., Lee, D., Mendoza, B., Nakajima, T., Robock, A., Stephens, G., Takemura, T., Zhang, H. 2013. Anthropogenic and Natural Radiative Forcing. In: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [Stocker, T.F., D. Qin, G.‐K., Plattner, M. Tignor, S.K. Allen, J. Boschung, A. Nauels, Y. Xia, V. Bex and P.M. Midgley (eds.)].Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. IUNG, the Institute of Soil Science and Plant Cultivation, data on the Climatic Water Balance – KBW: http://www.susza.iung.pulawy.pl/KBW/. Krupnik, T.J., Ahmed, Z.U., Timsina, J., Yasmin, S., Hossain, F., Al Mamun, A., Mridha, A.I., McDonald, A.J., 2015. Untangling crop management and environmental influences on wheat yield variability in Bangladesh: an application of non-parametric approaches. Agric. Syst. 139, 166–179. Li, X., Hu, C., Delgado, J.A., Zhang, Y., Ouyang, Z., 2007. Increased nitrogen use efficiencies as a key mitigation alternative to reduce nitrate leaching in north China plain. Agr. Water Manage. 89, 137–147. Linquist, B., van Groenigen, K.J., Adviento-Borbe, M.A., Pittelkow, C., van Kessel, C., 2012. An agronomic assessment of greenhouse gas emissions from major cereal crops. Glob. Change Biol. Bioenergy 18, 194–209. Liu, C., Wang, K., Zheng, X., 2013. Effects of nitrification inhibitors (DCD and DMPP) on nitrous oxide emission, crop yield and nitrogen uptake in a wheat-maize cropping system. Biogeosciences 10, 2427–2437. Lv, L., Wang, H., Jia, X., Wang, Z., 2011. Analysis on water requirement and water-saving amount of wheat and corn in typical regions of the North China Plain. Front. Agric. China 5 (4), 556–562. Mäkinen, H., Kaseva, J., Trnka, M., Balek, J., Kersebaum, K.C., Nendel, C., Gobin, A., Olesen, J.E., Bindi, M., Ferrise, R., Moriondo, M., Rodríguez, A., Ruiz-Ramos, M., Takáč, J., Bezák, P., Ventrella, D., Ruget, F., Capellades, G., Kahiluoto, H., 2018. Sensitivity of European wheat to extreme weather. Field Crop. Res. 222, 209–2017. Mase, A.S., Gramig, B.M., Prokopy, L.S., 2017. Climate change beliefs, risk perceptions, and adaptation behavior among Midwestern U.S. Crop farmers. Clim. Risk Manag. 15, 8–17. Mohammadi, R., Roustaii, M., Haghparast, R., Roohi, E., Solimani, K., Ahmadi, M., Abedi, R., Amri, A., 2010. Genotype × environment interactions for grain yield in rainfed winter multi-environment trials in Iran. Agron. J. 102, 1500–1510. Oerke, E., 2006. Crop losses to pests. J. Agr. Sci. 144 (1), 31–43. Pachauri, R.K., Allen, M.R., Barros, V.R., Broome, J., Cramer, W., Christ, R., Church, J.A., Clarke, L., Dahe, Q., Dasgupta, P., Dubash, N.K., Edenhofer, O., Elgizouli, I., Field, C.B., Forster, P., Friedlingstein, P., Fuglestvedt, J., Gomez-Echeverri, L., Hallegatte,

29