Xioxing Zhn,Hui Sho,Win Zhng,Wig Huo,Willim Dvi Bthlor,Png Hou,Enli Wng,Guohu Mi,Yuxin Mio,Higng Li,*,Fusuo Zhng

aCollege of Resources and Environmental Sciences,China Agricultural University,Beijing 100193,China

bBiosystems Engineering Department,Auburn University,Auburn,AL 36849,USA

cInstitute of Crop Sciences,Chinese Academy of Agricultural Sciences,Key Laboratory of Crop Physiology and Ecology,Ministry of Agriculture,Beijing 100081,China

dCSIRO Agriculture Flagship,Canberra,ACT 2601,Australia

eDepartment of Soil,Water and Climate,Minnesota University,St.Paul,MN 55108,USA

Keywords:Hybrids Leaf area Leaf number Plant density Environmental conditions

A B S T R A C T Accurate leaf area simulation is critical for the performance of crop growth models.Area of fully expanded individual leaves of maize hybrids released before 1995(defined as old hybrids)has been simulated using a bell-shaped function (BSF)and the relationship between its parameters and total leaf number(TLNO).However,modern high-yielding maize hybrids show different canopy architectures.The function parameters calibrated for old hybrids will not accurately represent modern hybrids.In this study,we evaluated these functions using a dataset including old and modern hybrids that have been widely planted in China in recent years.Maximum individual leaf area(Y0)and corresponding leaf position(X0)were not predicted well by TLNO(R2=0.56 and R2=0.70)for modern hybrids.Using recalibrated shape parameters a and b with values of Y0 and X0 for modern hybrids,the BSF accurately predicted individual leaf area(R2=0.95–0.99)and total leaf area of modern hybrids(R2=0.98).The results show that the BSF is still a robust way to predict the fully expanded leaf area of maize when parameters a and b are modified and Y0 and X0 are fitted.Breeding programs have led to increases in TLNO of maize but have not altered Y0 and X0,reducing the correlation between Y0,X0,and TLNO.For modern hybrids,the values of Y0 and X0 are hybrid-specific.Modern hybrids tend to have less-negative values of parameter a and more-positive values of parameter b in the leaf profile.Growth conditions,such as plant density and environmental conditions,also affect the fully expanded leaf area but were not considered in the original published equations. Thus, further research is needed to accurately estimate values of Y0 and X0 of individual modern hybrids to improve simulation of maize leaf area in crop growth models.

1.Introduction

Accurate simulation of leaf area is critical for crop growth models[1,2].Two approaches based on leaf number have been used to simulate the area of fully expanded individual leaves of maize in previous studies [1, 3–5]. In the first approach,fully expanded leaf area is simulated using four discrete equations for leaves 1–3,4–11,12 to TLNO-4,and TLNO-3 to TLNO[3].However,the area of leaves above the 12th leaf is often underestimated by this function,especially for high-yielding modern maize hybrids[6,7].

The second approach is independently described by three equations(Eqs.(1)–(3),Table 1).The parameters in these equations are typically fitted to measured field data and have shown in previous studies to give good predictions of leaf area for a range of hybrids growing in different soils[1, 5, 8]. Keating and Wafula found that the values of parameters Y0,X0,a,and b can be estimated from TLNO(Eqs.(4)–(7),Table 1).The BSF(Eq.(1),Table 1)used with parameters estimated by Eqs.(4)–(7)gave good predictions of individual leaf area and total leaf area of maize grown in moderate temperature regions with a low plant population of 2.2 plants m?2[9].

Birch et al.[10]modified the relationships between TLNO and the four parameters using five maize hybrids grown in subtropical conditions with a plant density of 7 plants m?2(southeastern Queensland,Australia)(Eqs.(8)–(11),Table 1).The model simulated leaf area well under the optimal temperature for leaf expansion.However,validation using an independent data set was not performed[10].Elings[11]calibrated Eqs.(4)–(7)(Table 1)for tropical cultivars growing under tropical conditions at various water and nitrogen levels at a plant density of 5.33 plants m?2.However,coefficients of determination (R2) were too low to justify estimation of function parameters from TLNO.It is likely that leaf expansion was reduced by limitations imposed by water,N and temperature[11–13].

In previous studies,calibrations of the BSF and its parameters have used older maize hybrids released before 1995 and grown at low plant densities[9–11].Because newer high-yielding maize hybrids have different canopy architecture that is sensitive to plant density and environmental condition,it is expected that function parameters calibrated for old hybrids will not accurately represent modern hybrids.The objective of this study was to calibrate and evaluate the BSF and its parameters for modern hybrids, and to develop a generalized method to estimate individual leaf area and total leaf area for a wide range of hybrids differing in total leaf number grown in different locations and environmental conditions without nutrient or water stress.

Table 1–Principal equations used for estimation of fully expanded leaf area in maize.Y0 is the area of the largest leaf,X0 is the position of the largest leaf,a and b are dimensionless empirical constants.

2.Materials and methods

2.1.Experimental site and design

Two experiments were conducted to collect leaf area data for modern maize hybrids for model calibration(Table 2).The first experiment (experiment 1) was conducted at the Quzhou Experiment Station,Hebei province,China(36°52′N,115°02′E).Maize(cv.Zhengdan 958)was sown at a density of 7.5 plants m?2.There were three phosphorus(P)levels with four replicates in a randomized complete block,including P0(no P application),P75(75 kg P2O5ha?1),and P300(300 kg P2O5ha?1)as superphosphate,but all leaf area data were collected from treatments P75and P300,which provided sufficient P to the crop.Overall fertilizer applications were high enough so that N,P,or K deficiencies would not limit maize growth, based on previous results.Irrigation was applied to ensure no water stress during the growing season.

The second experiment(experiment 2)was conducted at Lishu county,Jilin province,China(43°21′N,124°04′E).Six hybrids commonly grown in northeast China were evaluated to measure leaf expansion.All treatments were replicated three times in a randomized complete block design.Rainfall was 736 mm during the growing season in 2016,an amount much higher than the evapotranspiration of maize production in this area and ensuring that water demand was met. All field experiments were rigorously controlled for pests and diseases.

Table 2–Site and design of experiment 1 at Quzhou,China and experiment 2 at Lishu,China.Leaf area data from experiments 1 and 2 were used for model calibration.Maize was supplied with sufficient nutrients and water to prevent stress in both experiments.

2.2.Field sampling and data collection

In experiment 1(2016,Quzhou),three plants were randomly selected and destructively sampled in each plot at four growth stages(V4,V9,silking,and grain filling).Each leaf position was recorded as the leaf number counted from the bottom.In experiment 2 (2016, Lishu), five plants were tagged after emergence in each plot.Leaf area and leaf position were measured and recorded at V5,V7,silking,and grain filling.Length and maximum width of fully expanded leaves were measured.Length was determined from the collar to the tip of fully expanded leaves and width was measured at the widest part of the leaves.It was assumed that leaf length and width did not change after full expansion. The fully expanded individual leaf area was calculated as length×maximum width×0.75[10,14,15].The value of 0.75 was considered anacceptable average of reported values:0.72[9],0.73[1],0.75[15], 0.79 [10]. Total leaf number per plant (TLNO) was recorded after the flag leaf had appeared.

Table 3–Leaf area data assembled from the literature for old and new hybrids from different sites around the world.Data from study 1 were used for model calibration and data from studies 2–9 were used for model evaluation.Data from studies 10 and 11 were used to analyze canopy characteristics under different densities and environmental conditions.Maize was supplied with sufficient nutrients and water to eliminate stresses in all experiments.

Table 4–Error in simulated total leaf area(TLA,cm2)using the bell-shaped function(BSF)(Eq.(1),Table 1)with three different methods of estimating parameters.

2.3.Data collection from published studies

Fig.1–Area of fully expanded individual leaves(LAn,cm2)observed in experiment 1 at Quzhou and experiment 2 at Lishu(symbols),and simulated(lines)using the bell-shaped function(BSF)(Eq.(1),Table 1)in combination with the relationship coefficients developed by Keating(a)(b)(c).Simulated total leaf area per plant(TLA,cm2)was compared with observed total leaf area per plant(TLA,cm2)in experiments 1 and 2(d).

Data of individual leaf area and total leaf number from published studies were assembled.These experiments represented old and new hybrids from around the world(Table 3)[9,10,13,16–23].These experiments were selected because they had been managed to eliminate nutrient and water stresses.The data were reported as graphs in the literature and the data from these graphs were extracted using GetData Graph Digitizer(version 2.22,developer:S.Fedorov).To avoid variation from calculation,leaf area taken from the literature was recalculated as length×maximum width×0.75 if the value of the coefficient was not 0.75.The data from study 1(1985,Katumani,Kenya)were used for calibration of Eqs.(4)–(7)(Table 1),and data from studies 2–9(Table 3)were used for evaluation.Data from study 10(1996–1998,Pergamino and Salto,Argentina)and study 11(2012,Luohe and Jiamusi,China)were used to analyze canopy characteristics under different plant densities and environmental conditions.

Fig.2–Area of fully expanded individual leaves(LAn,cm2)observed in experiment 1 at Quzhou and experiment 2 at Lishu(symbols)and simulated(lines)using the bell-shaped function(BSF)(Eq.(1),Table 1)in combination with the relationship coefficients for parameters a,b developed by Keating,but Y0,X0 fitted to measured data(a)(b)(c).Simulated total leaf area per plant(TLA,cm2)was compared with observed total leaf area per plant(TLA,cm2)in experiments 1 and 2(d).

2.4.Data analysis

The observed individual leaf areas from study 1(1985,Katumani)and experiments 1(2016,Quzhou)and 2(2016,Lishu),were combined to fit BSFs(Eqs.(1)–(3),Table 1)to determine one set of parameters Y0,X0,a,and b using the nonlinear least squares algorithm.The relationships between the four parameters of the BSF and TLNO were also determined by nonlinear least squares regression using Microsoft Excel(version 2010,Microsoft Corporation,Redmond,Washington,USA).

The root mean square deviation(RMSD)and coefficient of determination(R2)were used to assess the degree of fit of the functions between observed(O)and simulated(S)values.The RMSD indicates an average weighted difference between observed and simulated values.where n is the number of O ?S pairs and Oais the mean of the observed values.

Coefficient of variation (CV, %) was calculated using procedures in Microsoft Excel(version 2010,Microsoft Corporation,Redmond,Washington,USA).

3.Results

3.1. Evaluation of Keating and Wafula's coefficients for modern hybrids

The prediction accuracies of the three BSFs(Eqs.(1)–(3),Table 1)were compared in simulating leaf area using the datasets from Quzhou and Lishu in 2016(Table 4,Tables S1,and S2).The TLNO values of the modern hybrids were 19 for Zhengdan 958 at Quzhou,and 21 for Xianyu 335 and 23 for Liangyu 99 at Lishu. Simulated individual leaf area for modern hybrids using Keating and Wafula's coefficients and equations(Eqs.(4)–(7),Table 1)combined with the BSF(Eq.(1),Table 1)did not match observed values well(Fig.1-a,b,c).The R2ranged from 0.62 to 0.80 across the three hybrids.This poor fit was due primarily to overestimation of Y0and X0.Simulated Y0was 943 cm2for Zhengdan 958,1067 cm2for Xianyu 335,and 1191 cm2for Liangyu 99,values much higher than the observed Y0values,by 30%,32%,and 47%,respectively.The simulated average of X0was higher by 2.1 leaves compared to observed values for the three hybrids. The discrepancy between simulated and observed X0was smallest for Zhengdan 958,at 1.3 leaves(Fig.1-a).Thus,the total leaf area(TLA)of all modern hybrids was overestimated by 14%–26%(Fig.1-d,R2=0.96)using the coefficients developed by Keating and Wafula for old hybrids.

Fig.3–Relationships between parameters Y0(a),X0(b),a(c),and b(d)of the bell-shaped function(BSF)(Eq.(1),Table 1)and total leaf number(TLNO)based on hybrids from study 1(Katumani)and experiments 1(Quzhou)and 2(Lishu)in the absence of water and nitrogen deficits.

In the next step we used observed values of Y0and X0along with Eqs.(6)and(7)(Table 1)for the same three modern hybrids.This new set of equations gave much better predictions of individual leaf area for the three modern hybrids(Fig.2-a,b,c).The R2values for simulated and observed individual leaf area were 0.97 for Zhengdan 958,0.97 for Xianyu 335,and 0.96 for Liangyu 99.However,total leaf area was underestimated,especially for Liangyu 99(by 16%)because of the narrow shape of the bell-shaped curve in simulation caused by small values of parameters a and b(Fig.2-d).

Eq.(1)gave the best simulations(Fig.3)compared to Eqs.(2)and(3)(Figs.S1 and S2)and yielded less error in simulated total leaf area(TLA)(Table 4)than Eqs.(2)and(3)(Tables S1 and S2).Accordingly,evaluation was conducted only for Eq.(1)in this study.

3.2.Modified relationship coefficients

New coefficients were developed for Eqs.(4)–(7)(Table 1)using the three datasets(study 1,experiments 1 and 2).To evaluate the robustness of the new estimates,data were selected from study 1(1985,Katumani;Table 3)and field experiments at Quzhou and Lishu that provided a range of TLNO from 12 to 23.The coefficients in Eqs.(4)–(7)were estimated using the nonlinear least squares algorithm in Microsoft Excel 2010,which were then used with the BSF(Eq.(1),Table 1)to develop new estimates for Y0and X0,a and b and ultimately individual leaf area and total leaf area per plant.The new equations were(Fig.3-a,b,c,d):

The best linear fit was found for the relationships between Y0and X0,vs TLNO(Fig.3-a,b),while the best-fitting curve was found for the relationships between a and b and TLNO(Fig.3-c,d).The R2was only 0.56 for Y0and 0.70 for X0.The Y0and X0values did not linearly increase when TLNO was ＞16.Thus,one line does not represent the relationships between Y0,X0,and TLNO well.In contrast,the calibrated model gave high correlations between a,b,and TLNO,with R2values of 0.99 for a and 0.97 for b.

Fig.4–Area of fully expanded individual leaves(LAn)observed in studies 2–9(symbols)and simulated(lines)from the bellshaped function(BSF)(Eq.(1),Table 1)with parameters a and b estimated according to Fig.3-c,d,but Y0,X0 fitted to observed values(a–e).Simulated total leaf area per plant(TLA,cm2)was compared with observed total leaf area per plant(TLA,cm2)in Studies 2–9(f).

3.3.Evaluation of the new coefficients

Eqs.(14)–(17)were evaluated using the data from studies 2–9 using three methods (Table 4). In method 1, all four parameters,including Y0,X0,a,and b of the BSF,were derived from the new calibrations in Eqs. (14)–(17). In method 2,observed Y0was used while other parameters were derived from estimations using Eqs.(15)–(17).In method 3,observed Y0and X0were used,while parameters a and b were estimated by Eqs.(16)and(17).

RMSD of individual leaf area(RMSDL)ranged from 64 cm2to 232 cm2with an average of 138 cm2between observed and simulated individual leaf area using method 1(Table 4).The average RMSDL(aveRMSDL) decreased as more observed values of parameters were used.The aveRMSDLdecreased to 113 cm2for method 2 and 36 cm2for method 3.All RMSDLvalues were lower than 48 cm2across the experiments for method 3.RMSD of total leaf area(RMSDT)showed the same trend.The RMSDTwas 1216 cm2for method 1,332 cm2for method 2,and 260 cm2for method 3.

The individual leaf areas were simulated well by the BSF(Eq.(1),Table 1)using method 3 for the seven hybrids grown at different sites(Table 4),including Lincoln,GL420,Volga,DK636,Hycorn 83,Nongda 486 and Dekalb XL82(Fig.4-a–e).The R2values for the regression of simulated on observed individual leaf area were 0.95 to 0.99(Fig.4-a–e),which were much higher and more stable than the R2values of 0.20–0.99 for method 1 and 0.20–0.99 for method 2(data not shown).The total leaf areas(TLAs)in all eight studies shown in Table 3 were also simulated well by method 3,giving an R2of 0.98(Fig.4-f).

Table 5–Coefficient of variation(CV,%)for maximum individual leaf area(Y0),leaf position of maximum leaf area(X0)and total leaf number(TLNO)for six maize hybrids from experiment 2 at Lishu.

3.4.Variation of parameters

In experiment 2,the observed range of TLNO across replications was 21–23 for Liangyu 918,20–21 for Jinqing 202,20–22 for Nonghua 101,and 21–22 for Shengrui 999(Table 5).Only Liangyu 99(23 leaves)and Xianyu 335(21 leaves)appeared to produce the same number of leaves across replications.Parameter Y0appeared to be a hybrid trait,and was different for hybrids that produced the same total leaf number.Hybrid Liangyu 99 had a maximum Y016%larger than that of Liangyu 918 but both had 23 leaves.Jingqing 202 with 21 leaves had a maximum Y0of 698 cm2,which was only 83%that of Shengrui 999.No increase in Y0was observed with increasing TLNO for Liangyu 918,Jingqing 202,and Shengrui 999.The range of X0was 9.0–12.3 for the hybrids in Table 5.Shengrui 999 always had a larger X0value than other hybrids with the same total leaf numbers except Jingqing 202 with 21 leaves.No significant positive correlation was found between X0and TLNO for Liangyu 918 and Nonghua 101. The X0of Nonghua 101 decreased from 11.5 to 9.0 when TLNO increased from 21 to 22.The averages coefficients of variation(CV,%)were 3.99%for Y0and 3.24%for X0.

Plant density did not significantly change TLNO for DK4F37 and DK752 from the Pergamino,Argentina experiment and DK969 and DK757 from the Salto, Argentina experiment(Table 3;study 10)even when plant density ranged from 3 to 12 plants m?2(Fig.5-a,b).In contrast,Y0decreased by 15%at Pergamino and 22%at Salto with increasing plant density from 3 to 12 plants m?2(Fig.5-c,d).The slope of Y0was much steeper at Salto than at Pergamino. Increasing plant density also reduced X0at both Pergamino(7%)and Salto(7%)(Fig.5-e,f).The decrease of Y0due to increase of plant density was much higher than X0at both Pergamino and Salto(Fig.5-c–f).

Hybrid Zhengdan 958 was sowed in Luohe(Henan province,China,33°35′N,114°01′E)and Jiamusi(Heilongjiang province,China,46°48′N,130°19′E)at a plant density of 6 plants m?2in 2012(Table 3;Study 11).A lower TLNO was observed at Luohe(18 leaves)than Jiamusi(21 leaves)(Fig.6-a).The average value of Y0in Luohe was smaller than that in Jiamusi by 33%(Fig.6-b).The average value of X0in Luohe was 13.3 but in Jiamusi was 12.0,lower than that at Luohe(Fig.6-c).

4.Discussion

4.1.Fitting individual leaf area of modern hybrids with old equations

The original equations of Keating and Wafula(Eqs.(4)–(7)as empirical functions fail to simulate Y0,X0and the bell-shaped curve of modern hybrids(Figs.1 and 2)[9].This is likely because the TLNO of modern hybrids is beyond the scope of the original leaf area model calibration.Keating and Wafula found that Y0and X0showed a linear relationship with TLNO with a range of 12 to 17 leaves for maize under non-limiting conditions[9].However,the TLNO of modern hybrids widely planted in China since 1995 has increased to ＞18 leaves(Fig.1,Tables 2 and 3).The measured Y0and X0values did not continue increasing linearly with increasing TLNO over 16 leaves as the original equations described(Fig.3-a,b)for these new hybrids(2016,Quzhou,China;2016,Lishu,China).

Using the data of Keating and Wafula[9],Birch et al.developed a nonlinear regression of Y0and TLNO(Eq.(8)).This equation could simulate maximum individual leaf area for plants with ＞16 leaves in Birch's study[10].However,it did not simulate Y0well for modern hybrids in this study because it overestimated Y0of modern hybrids when TLNO is above 16.In the same study,X0was found to have a correlation with TLNO(Eq.(9)),which did not fit modern hybrids well because X0of modern hybrids did not increase when TLNO was above 16(2016,Quzhou;2016,Lishu)(Fig.3-b)[10].

4.2.Correlations between Y0,X0,and TLNO for modern hybrids

Fig.5–Relationships between total leaf number(TLNO),maximum individual leaf area(Y0),leaf position of maximum leaf area(X0)and density(D)based on hybrids in study 10.Two hybrids,DK4F37 and DK752,were sown at three plant densities during 1996–1997 at Pergamino(a)(c)(e)and the two hybrids DK969 and DK757 were planted at three plant densities during 1997–1998 at Salto(b)(d)(f).Values are the means of two hybrids.

The parameters Y0and X0of new hybrids did not linearly increase with TLNO as in old hybrids,but tended to plateau(Fig.3-a,b).Modern hybrids have a lower ear and ear leaf position than old hybrids [24, 25]. As low ear position decreases the likelihood of lodging. Breeders tend to select hybrids with a low ear position.Thus,the parameter X0of modern hybrids is more hybrid-specific than old ones,and may have a poor correlation with TLNO(Fig.3-b;Table 5).

High plant density is a strategy for maximizing yield in high-yield management.TLNO in maize is determined by thermal time and photoperiod between the fourth and sevenleaf stage[26],and is not sensitive to plant density(Fig.5-a,b).However,high plant density significantly reduces area of leaves in the vicinity of the ear including Y0due to competition,but has little effect on the early leaves in many hybrids(Fig.5-c,d)[22,27].High plant density reduces radiation quality and quantity for each leaf and intensifies competition for radiation,leading to a decrease in Y0due to inadequate photosynthate supply[6,10,26].

Temperature also appears to affect the Y0of hybrids.The hybrid Zhengdan 958 was sown at different sites from 32°52′N to 48°08′N in China from 2010 to 2012[23].Average daily temperatures decreased greatly from south to north [23].Environmental variables,particularly temperature in these sites,led to parameter Y0variations for the hybrid Zhengdan 958[11,17].Higher temperatures in low latitude regions gave smaller Y0than those in high latitude regions(Fig.6-b),the same result reported by Reid et al.[28].

Fig.6–Observed total leaf number(TLNO),maximum individual leaf area(Y0),and leaf position of maximum leaf area(X0)for hybrid Zhengdan 958 in study 11(a)(b)(c).Zhengdan 958 was sown at a plant density of 6 plants m?2 in 2012 at Luohe(Henan province,China)and Jiamusi(Heilongjiang province,China).Error bars represent standard error(n=3).Asterisk indicates significance at**P ＜0.01,and*P ＜0.05.

Increasing plant density also enhances competition for assimilation between stem and leaves,leading to small values of Y0at low X0values(Fig.5-c–f).The ideal canopy structure for high plant density is low leaf area index above the ear to alleviate shading of middle leaves by upper leaves[5,25,27].Selection for high-yielding hybrids under high plant density by plant breeders may have led to changes in the relationship between parameter Y0and TLNO(Fig.3-a;Table 5).

As Keating and Wafula[9]mentioned,growing conditions could change the relationship between the parameters of the BSF and TLNO,which can be simulated under non-limiting conditions,because TLNO has a major influence on these parameters.In the present study,we found that parameters Y0and X0cannot be accurately predicted by TLNO for modern hybrids,likely because breeding has led to changes in these relationships(Fig.3-a,b;Table 5).The failure of Y0and X0as predicted by Eqs.(14)and(15)(Fig.3-a,b)may be due to the limited variation of Y0and X0with TLNO across diverse hybrids and sites.In addition,parameters Y0and X0are sensitive to genetic improvement and environmental factors(Table 5;Fig.4-a–e;Fig.6-b,c).To simulate individual leaf area of a specific hybrid more precisely, we suggest that the parameters Y0and X0be fitted using known values from experiments(Table 4,method 3).

4.3.Correlations between a,b,and TLNO for modern hybrids

In recent years,maize breeding programs for high grain yield have resulted in hybrids that maintain high individual leaf radiation interception at high plant density.We improved the original equations to estimate parameters a and b using our dataset from study 1(1985,Katumani)and experiments 1 (2016, Quzhou) and 2 (2016, Lishu),and found the new calibrations gave strong correlations with TLNO(Fig.3-c,d).The parameters a and b increase progressively to a maximum value with increasing values of TLNO. Modern hybrids showed a wider bell-shaped curve with a lower slope than older hybrids (Fig. 2-c).This finding indicates that leaf area at lower and upper positions showed a greater increase over time than leaves at middle positions in newer hybrids.

5.Conclusions

The BSF appears to be robust for simulating the fully expanded area of maize leaves with modified parameters a,b and fitted Y0and X0across older and more modern hybrids.We found that modern hybrids have broader(less negative values of a)and more skewed(more positive values of b)leaf profiles.The parameters Y0and X0are hybrid-specific for modern hybrids and are affected by plant densities and environmental conditions.Further research should quantify the impacts of these factors on Y0and X0of modern hybrids.

Acknowledgments

This work was supported by the National Basic Research Program of China(973-2015CB150400)and by the National Institute of Food and Agriculture(ALA014-1-16016),U.S.Department of Agriculture,Hatch project under ALA014-1-16016.

Appendix A.Supplementary material

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2018.03.008.