lmage-based root phenotyping for field-grown crops:An example under maize/soybean intercropping

HUI Fang ,XlE Zi-wen ,LI Hai-gang ,GUO Yan ,LI Bao-guo ,LlU Yun-ling ,MA Yun-tao

1 College of Land Science and Technology,China Agricultural University,Beijing 100193,P.R.China

2 College of Grassland,Resource and Environment,Inner Mongolia Agricultural University,Hohhot 010011,P.R.China

Abstract Root architecture,which determines the water and nutrient uptake ability of crops,is highly plastic in response to soil environmental changes and different cultivation patterns. Root phenotyping for field-grown crops,especially topological trait extraction,is rarely performed. In this study,an image-based semi-automatic root phenotyping method for field-grown crops was developed. The method consisted of image acquisition,image denoising and segmentation,trait extraction and data analysis. Five global traits and 40 local traits were extracted with this method. A good consistency in 1st-order lateral root branching was observed between the visually counted values and the values extracted using the developed method,with R2=0.97. Using the method,we found that the interspecific advantages for maize mainly occurred within 5 cm from the root base in the nodal roots of the 5th-7th nodes,and that the obvious inhibition of soybean was mostly reflected within 20 cm from the root base. Our study provides a novel approach with high-throughput and high-accuracy for field research on root morphology and branching features. It could be applied to the 3D reconstruction of field-grown root system architecture to improve the inputs to data-driven models (e.g.,OpenSimRoot) that simulate root growth,solute transport and water uptake.

Keywords:root phenotyping,high-throughput,image analysis,intercropping,maize (Zea mays L.),soybean (Glycine max L.)

1.lntroduction

Root systems determine plant water and nutrient uptake,and ultimately affect plant growth and yield (Lynch 1995;Blum 2005;De Smetet al.2012;Chenet al.2019;Lobetet al.2019). Root phenotypes are highly plastic in response to soil environmental changes (Unger and Kaspar 1994;Dexter 2004;Wuet al.2021;Yanget al.2021) and cultivation patterns (Muonekeet al.2007;Duet al.2018).The container size and growth medium also have strong influences on root development and morphology,and have led to nonnegligible differences in root phenotypes between pot-and field-grown plants (Tianet al.2017;De la Fuenteet al.2018;Sunget al.2019). Therefore,field phenotyping remains a key focus of root research.

Root phenotyping in the field is a major challenge because soil is heterogeneous and opaque. Numerous attempts have been made to obtain information about root development and morphology in the field. Among these methods,auger sampling is commonly used to obtain partial biomass and phenotypic traits of the root system to characterize the overall root extension and distribution by root length density (RLD) (Gajriet al.1994;Liet al.2006;Gaoet al.2010;Liet al.2016). Studies have shown that the number,position and size of sampling holes strongly influence RLD estimates (Rossi and Nuutinen 2004;Wu Qet al.2018). Minirhizotrons are transparent interfaces designed to observe local root growthin situand have been applied to investigate the influences of drought,irrigation and fertilizer treatments on the RLD distribution in the soil profile (Liedgens and Richner 2001;Hanet al.2016;Chenet al.2018). Shovelomics is a global root measurement method for rapidly observing the number,angle and simple branching pattern of the root crown in maize (Trachselet al.2011),legumes (Burridgeet al.2016) and wheat (Slacket al.2018;Yorket al.2018). Recently,an integrated method that includes field sampling system for larger root cores and phenotypic analysis has been proposed for field-grown maize (Wu and Guo 2014). The topological structure of lateral roots (LRs) in a relatively complete root system can be quantified and analysed with this method. However,manual correction of image processing results from WinRHIZO is extremely time-consuming and limits its potential application in root phenotypic research.

Because of its low cost and high efficiency,image-based analysis provides more possibilities for high-throughput root phenotyping.It has been widely used in root phenotyping platforms in a wide variety of root growth systems,including soil-filled (Joshiet al.2017;Bodneret al.2018;Gaoet al.2019),hydroponic/semi-hydroponic (Jeudyet al.2016;Mohamedet al.2017;Wu Jet al.2018) and gel/agar-based systems (Clarket al.2013;Yasrabet al.2019). These platforms are non-destructive and high-throughput,but they cannot be used to observe a whole mature root system.Along with the root phenotyping platforms,a series of image processing software packages have been developed,such as GiA Roots (Galkovskyiet al.2012),AutoRoot (Poundet al.2017) and GT-Roots (Borianneet al.2018). These tools mainly focus on global root traits (Seethepalliet al.2020) and are not good choices for analysing images from field experiments (Betegón-Putzeet al.2019;Yasrabet al.2019). Therefore,for the field root research,constructing a method of image processing and trait extraction is a current bottleneck.

Intercropping is an ancient practice that has been used for centuries because of its higher temporal and spatial resource utilization efficiency (Romeroet al.2013) and grain yield advantages (Fanget al.2013;Yinet al.2015) compared to monocropping. Intercropping can also improve biodiversity and soil fertility,and reduce diseases,insects and weeds(Zhuet al.2000;Liet al.2007;Wanget al.2014;Conget al.2015). Intensive studies have analysed biomass,yield and morphological responses to aboveground interactions between species (Muonekeet al.2007;Liuet al.2017;Fanet al.2018). Furthermore,image-based aboveground phenotyping has also been attempted in the field for the 3D quantification of intercropping crops (Burgesset al.2017;Zhuet al.2018). In contrast,field root phenotyping is rarely reported,especially in intercropping systems.Consequently,the aims of this study are to:(1) develop an image-based root phenotyping method for field-grown crops;and (2) assess and analyse the global and local root phenotypic traits,including morphology and branching features,of sole crops and intercrops.

2.Materials and methods

2.1.Field experiments

Field experiments were carried out from May to October 2017 at the Lishu Experimental Station of China Agricultural University (43°16′N,124°26′E,196 m a.s.l.),Jilin Province,China. The soil texture is a silty clay loam (USDA soil classification scheme),and the soil layer is more than 2 m.The experiments were set up in a completely randomized block design with three treatments and three replicates.The three treatments were sole maize,sole soybean and intercropped maize/soybean. The plot size was 24 m×10 m.The maize (Zea maysL.) cultivar was Xianyu 335(XY335),and the soybean (Glycine maxL.) cultivar was Jiyu 47 (JY47). The intercropping consisted of two rows of maize and two rows of soybean. The interplant distance was 20 cm for maize and 10 cm for soybean. The interrow distance was 50 cm for both intercrops and sole crops.The date of sowing was May 11. Before sowing,chemical fertilizers were applied by mechanical rotary tillage at a rate of 80 kg N ha-1(urea),120 kg P2O5ha-1(potassium sulfate)and 100 kg K2O ha-1(calcium superphosphate). Urea was top-dressed for the maize stands at the maize jointing and belling stages,with 80 kg N ha-1for sole maize and 40 kg N ha-1for intercropped maize at each stage.

Roots were sampled at 17,35,56,and 84 days after emergence (DAE). Each plot represented one sampling.Three samplings were performed for each treatment. The root sampling system designed by Wu and Guo (2014)was used to obtain the root systems in the field. The system comprised a hammer module (Fig.1-A),a sampling cylinder (Fig.1-B) and a lift module (Fig.1-C). The diameter(Fig.1-B) and height of sampling cylinder were both 50 cm.After the shoot of a plant to be sampled was cut off,the sampling cylinder was placed on the soil surface with the plant in the centre of the sampling cylinder. The hammer was struck on the cross bearing to gradually drive the cylinder into the soil (Fig.1-A). When the cylinder was completely in the soil (Fig.1-B),the lift module was used to remove the cylinder from the soil (Fig.1-C and D). The sampling time of one root core was approximately 20 min. After removing the cylinder (Fig.1-E),the roots were cleaned with a group of adjustable water jets (Fig.1-F and G). About 2 h was required for washing out each sample. The root system was then stored in a refrigerator at 4°C.

Fig.1 Image-based root phenotyping method for field-grown crops.The method includes image acquisition by two strategies (A),image border and background denoising (B),root segmentation (C),global (D),and local (E) root trait extraction and evaluation,and data processing and analysis (F).

At 17,35,26,41,56,84,and 113 DAE,the shoot and root biomasses sampled for the three treatments were also used to describe the plant growth dynamics. All samples were dried for 24-48 h at 75°C to a constant weight and then their dry biomass was measured.

2.2.Method frame

Image-based semi-automatic root phenotyping was developed with MATLAB R2018a Language (MathWorks,USA) and WinRHIZO Pro 2009b Software (Regent,Canada). The method consisted of four submodules:image acquisition,image processing,trait extraction,and data analysis (Fig.2). The image acquisition submodule involved root scanning with two strategies (Fig.2-A). The image processing submodule involved image denoising(Fig.2-B) and root segmentation (Fig.2-C). The trait extraction submodule included the extraction of global root traits with the batch processing function of WinRHIZO(Fig.2-D) and local trait extraction and evaluation (Fig.2-E).The data analysis submodule integrated the output files from MATLAB and WinRHIZO to analyse and visualize the root trait data (Fig.2-F).

Fig.2 Root sampling in the field. A,the hammer strikes the cross bearing. B,the sampling cylinder completely in the soil.C,lifting the sampling cylinder out of the soil. D,intact root core. E,removing the sampling cylinder. F,root washing. G,cleaned maize (left) and soybean roots (right).

2.3.lmage acquisition

Before root scanning,individual maize root systems were separated manually into nodal roots according to 1st-7th nodes,which were named N1 to N7 according to their appearance time (Fig.3-A). N6 indicates the nodal root in the 6th-node. A nodal root consisted of the main root and LRs (Fig.3-B). Three representative roots of each node were selected from three maize replicates for detailed scanning.

To reduce root overlap,the maize nodal roots and soybean roots were both cut into 5 cm root segments along the main root (Fig.3-D). Two strategies were adopted for detailed scanning depending on the number of LRs in each 5 cm segment. For root segments with a high LRs number,1st-order LRs were cut from the main root and then scanned together with main root (Fig.4-A). For root segments with a low LRs number,the root segments,including its LRs,were scanned directly (Fig.4-H). All roots were scanned at a resolution of 600 dots per inch (dpi) (ScanMaker i800 Plus,Microtek,China). For one person,about 0.5-3 days were required to scan a root depending on the complexity of LRs. After scanning,the roots were oven-dried for 24-48 h at 75°C (to constant weight) and weighed.

Fig.3 Diagrams of maize nodal root rank (A),root types in maize (B) and soybean (C),and each 5 cm root segment (D).

2.4.Image denoising and segmentation

A program based on the MATLAB Language was compiled for image denoising and segmentation. The inputs included two types of original image because of the two scanning strategies. Type one was a scanned image with a high LRs number (Fig.4-A);type two was a scanned image with a low LRs number (Fig.4-H).

Image denoising was performed on each original image to remove noise from the border and background (soil impurities and air bubbles) to obtain a denoised image. The border pixels were cut out based on the border width of the scanning plate (Fig.4-B and I). Background denoising was performed according to the following three steps:(1) Image binarization. The image was binarized based on the adapting thresholding algorithm with the ‘imbinarize’ function (Bradley and Roth 2007). (2) Connected component labelling. A connected component in a binary image is a set of pixels that form a connected group. Connected component labelling is the process of identifying the connected components in a binary image and assigning each one a unique label. The connected components were labelled based on a two-pass algorithm with the ‘bwlabel’ function (Haralick and Shapiro 1992). (3) Background denoising. Background noise was removed based on the parameters of Area (actual number of pixels in the region) and MajorAxisLength (length in pixels of the major axis of the ellipse that has the same normalized second central moments as the region) of the connected components (Fig.4-C and J). The parameters of Area and MajorAxisLength were set to 90 pixels (about 4 mm) and 30 pixels (about 1 mm),respectively.

Fig.4 Image denoising and segmentation. A,the original scanned image with a high lateral roots (LRs) number. B-G,the results of border denoising,background denoising,image erosion,the main root extracted with dilation and median filter algorithms,the lateral roots and segmentation visualization of the scanned image with a high LRs number. H,the original scanned image with a low LRs number. I-N,the results of border denoising,background denoising,image erosion,the main root,the lateral roots and segmentation visualization of the scanned image with a low LRs number. The main root pixels are cyan in the G and N images.

Root segmentation was performed according to the following two steps:(1) Main root segmentation (Fig.4-E and L). The main root pixels were separated from the denoised image based on the pixel coordinates of the maximum area of connected components,which were extracted based on the erosion (Fig.4-D and K),dilation (structural element was ones (3)) and median filter (filter size was five pixels;Fig.4-E and L) algorithms with the ‘bwmorph’ function.(2) LR segmentation. The LR pixels were obtained by subtracting the separated main root pixels from the denoised image (Fig.4-F and M). Then,LRs were separated from each other based on the labels of connected components.The segmentation result was visualized with the ‘label2rgb’function (Fig.4-G and N).

2.5.Extraction and evaluation of root traits at global and local scales

Root traits,e.g.,the root length,surface area,volume,average diameter and number of tips of all roots,main root,and LRs,were automatically extracted with the batch processing function of WinRHIZO Software. Root volume and surface area were both calculated from the length and diameter. The number of tips contained all levels of LRs(i.e.,1st-order,2nd-order,…). Global and local root traits were calculated for the entire root system,a nodal root and a 5-cm segment (Table 1).

WinRHIZO can only extract the Num1LRs by manual operation. We automatically calculated the number of 1storder LRs and number of 2nd-order and above LRs in the current study. The number of 1st-order LRs (Num1LRs) was counted based on the results of LR segmentation mentioned above with MATLAB (i.e.,counting the number of connected components of LRs in each root segment,Table 2). The number of 2nd-order and above LRs (Num2LRs) was calculated by subtracting Num1LRs from the number of root tips.

Table 1 Root traits at global and local scales extracted by our method

Table 2 Outputs of MATLAB

Based on 125 randomly selected scanning images of maize and soybean (at 56,73 and 89 DAE),comparisons between the artificially counted Num1LRs in each image and Num1LRs extracted from the method were made by using linear fitting. The slope,R2and root mean square error (RMSE) of the fitting equation were used to evaluate the accuracy of Num1LRs extracted with this method.

2.6.Data processing and visualisation

An R language program based on the package ‘tidyverse’(Wickhamet al.2019) was compiled for the output files from MATLAB and WinRHIZO to integrate and analyse the root trait data. The ‘ggplot2’ (Wickham 2016) and‘cowplot’ (Wilke 2019) packages were used for data visualization.t-test analysis (P＜0.05) was used to analyse all data with the package ‘STAT’. MATLAB and R codes used in this paper are available upon request to the corresponding author.

3.Results

3.1.Evaluation of Num1LRs extracted with the method

Comparing the segmentation results with the visually counted values for Num1LRs in the two crops in threegrowth stages,there was a good consistency,withR2=0.97,slope=0.99 and RMSE=2.56 (Fig.5-K). A visualization of the root segmentation for the 1st to 10th root segments of a nodal root is provided in Fig.5-A-J,in which the cyan pixels are the main root and the pixels of other colours are LRs.

Fig.5 Visualization (A-J) and evaluation (K) of root segmentation. A-J,a selected sample nodal root was segmented into 10 parts,and A-J are the 1st to 10th root segments of the nodal root in maize. Cyan pixels are the main root,and the remaining coloured pixels are lateral roots. K is the validation of the number of 1st-order lateral roots (Num1LRs) between the visually counted value and the segmentation result from the method,n is the number of scanning images. DAE,days after emergence.

3.2.Global root traits of maize and soybean during the whole growth period

The total root length (Fig.6-A and F),surface area (Fig.6-B and G),volume (Fig.6-C and H) and tips (Fig.6-D and I)increased rapidly for maize and soybean under the two planting patterns from 17 to 56 DAE. The values of these traits in intercropped maize were greater starting at 56 DAE,and 10-25% higher at 84 DAE,than those in sole maize(Fig.6-A-D). In contrast,intercropped soybean had lower values for those traits than sole soybean,and an almost 50% reduction occurred at 84 DAE (Fig.6-F-I). In both intercrops and sole crops,the maize root average diameter gradually decreased to approximately 0.5 mm (Fig.6-E) and soybean root average diameter increased initially and then decreased to approximately 0.4 mm (Fig.6-J). The total root average diameter in intercropping was slightly smaller than that in the sole crop for both maize and soybean at the late growth stage,but there were no significant differences.The total root average diameter represents the overall root thickness. In general,the thickness of roots is positively correlated with their ability to penetrate the soil (Materecheraet al.1992;Jinet al.2013).

Fig.6 Dynamic total root length (A and F),surface area (B and G),volume (C and H),tips (D and I),and average diameter (E and J) of maize and soybean root systems.

3.3.Local root traits of maize

To analyse the growth differences between sole maize and intercropped maize for different ranks of nodal root,we calculated the proportions of nodal root length (Fig.7-A)and nodal root tips (Fig.7-B) to those of total root traits in the entire maize root system. Regarding N1-N3,there were no obvious differences between sole crops and intercrops at 17 and 35 DAE. Intercropped maize had a higher proportion of nodal root length and nodal root tips than sole maize in N4-N5 at 35 DAE (5-15%,P＞0.05),in N5-N7 at 56 DAE(5-25%,P＞0.05) and 84 DAE (5-20%,P＜0.05). A 5-25%smaller proportion in intercropped maize than in sole maize was found in N1-N3 at either 56 DAE (P＞0.05)or 84 DAE(P＜0.05).

Fig.7 Proportions of nodal root length (A) and nodal root tips (B) to those of total root traits in the entire maize root system at 17,35,56 and 84 days after emergence. The nodal root ranks are from 1 to 7 (N1-N7). The red dotted lines show the intersegment point of variation in the proportions.

The root length and root tips of N4-N7 in each root segment were calculated to investigate the growth differences between sole maize and intercropped maize from the root base to the tip (Fig.8). The results showed that length and number of LRs were mainly distributed in the first 20 cm of the N4-N7. Overall,a relatively greater root length was found in N4-N7 of intercropped maize than of sole maize at 56 and 84 DAE,especially within 5 cm from the nodal root base (Fig.8-A). The results of root tips between intercropped maize and sole maize were basically consistent with those of root length,except for N5 at 84 DAE (Fig.8-B).

Fig.8 Root lengths and root tips of the 4th to 7th nodal roots (N4-N7) in each nodal root in each root segment from maize root base to tip at 56 and 84 days after emergence. Black points are averages.

3.4.Local root traits of soybean

The root length,Num1LRs and Num2LRs of each 5 cm root segment in 0-60 cm from the soybean root base at 56 and 84 DAE were shown in Fig.9-A-C. In the soybean root system,at least 70 and 90% of root traits’ values were located within 5 and 10 cm from the root base,respectively.No obvious difference was found in root segments at more than 20 cm. During the whole growth period,the dynamics of root length (Fig.9-D),Num1LRs (Fig.9-E) and Num2LRs(Fig.9-F) in the 0-5 cm root segment were consistent with those of the entire soybean root system,as shown in Fig.6-F and I. In the 5-20 cm root segment,the growth trend of LRs number was similar to root length during the whole growth period. Compared to the values of root length,Num1LRs and Num2LRs of intercropped soybean,these values were greater for sole soybean within the 10 cm root segment at 56 and 84 DAE and in the 10-20 cm root segment at 84 DAE.

Fig.9 Root length,number of 1st-order lateral roots (Num1LRs) and number of 2nd-order and above lateral roots (Num2LRs) in each 5 cm root segment of soybean. A-C,the results of the 0-60 cm root segment from the root base at 56 and 84 days after emergence. D-F,the results of each 5 cm root segment in 0-20 cm from the root base at 17,35,56 and 84 days after emergence.

4.Discussion

4.1.lmage-based root phenotyping method for fieldgrown crops

In this study,we developed a novel semi-automatic method of image-based root phenotyping for field-grown crops.The data acquisition process depended on human labor and the data processing process were improved.The automatic image processing and data extraction of the root branches in our method are complementary to the batch image analysis in WinRHIZO. Branching features can only be extracted manually in WinRHIZO Pro (Table 3),as WinRHIZO sometimes confuses the main root with LRs.Researchers must therefore repeatedly check for errors in each image to extract the root topological information (Wu and Guo 2014). If we analysed the root branching features with WinRHIZO in our study,manual correction would be needed for more than 10 000 scanned images. In contrast,our method had high throughput and high accuracy for the extraction of Num1RLs (Fig.5). It potentially enhances the efficiency of trait extraction for under 2nd-order LRs,which constitute 85% of the total LRs (Wu and Guo 2014).

Table 3 Time requirements of data extraction and processing by our method and WinRHIZO

Segmentation is a key step in image processing and analysis. The root scanning strategies adopted in our method provided high-quality original images with limited root overlap. Thus,the erosion and dilation algorithms and connected component analysis were simple but effective for this study (Siddiqiet al.2009;Moseset al.2016). For the analysis of early seedling roots with similar diameters,machine-learning approaches have been applied to segmentation and feature extraction and have shown a several-fold increase in speed compared with semiautomatic analysis (Yasrabet al.2019;Falket al.2020). However,this speed comes at the cost of a high dependence on image quality,a large number of image samples and a substantial amount of marking and training time to achieve high prediction accuracy. Unfortunately,even by machine learning,it is still difficult to achieve satisfactory accuracy at present on relatively complete mature roots with complex branching architecture analysis (Yasrabet al.2019;Falket al.2020). In contrast,improving the quality of the original scanned images is more feasible.

4.2.Global and local root phenotypic traits extracted using the method

Although image analysis approaches (Buckschet al.2014)and platforms (Daset al.2015) have been developed for field root phenotyping,these tools were specially designed for field-imaging protocols such as shovelomics,which mainly focus on global root traits. In this study,we quantified five global root traits and 40 local root traits using our method.Under maize/soybean intercropping,the global root traits(Fig.6) and biomass (Appendix A) extracted by our method indicated that interspecific interaction promoted intercropped maize growth (Fig.6-A-E) and inhibited intercropped soybean growth (Fig.6-F-J) in both shoots and roots(Appendix A). Similar findings were also reported in previous studies of maize/soybean intercropping under different planting conditions (Duet al.2018;Iqbalet al.2019),e.g.,planting density (Muonekeet al.2007),wide-narrow strips(Liu and Song 2012;Liuet al.2017;Renet al.2017),and nutrient treatments (Liuet al.2014;Lvet al.2014). This kind of interaction was defined as asymmetric interspecific facilitation (Liet al.2006).

The asymmetric details of the root phenotypic responses in different ranks of nodal roots (Fig.7) and from the root base to the tip (Figs.8-9) were also observable in the local root traits. Based on the local root traits extracted using the method,we further determined that the considerable advantages for intercropped maize roots were mainly reflected within 5 cm from the root base (Fig.8) in N5-N7(Fig.7),and that the inhibition of intercropped soybean roots mostly occurred within 20 cm from the root base (Fig.9).These results indicated that corn roots could obtain more water and nutrients from soil than soybean roots in the surface of nutrient-rich soil. This ability of the maize root system,resulting in asymmetric interspecific competition for land (Gaoet al.2010),water (Renet al.2017) and nutrients (Liuet al.2017) between the intercropped maize and soybean,also contributed to the increased grain yield and biomass of intercropping (Appendix A).

4.3.Potential applications of the method

The root phenotypic traits extracted with this method may improve the inputs to data-driven simulation models such as R-SWMS (Pohlmeieret al.2010),OpenSimRoot (Postmaet al.2017) and CRootBox (Schnepfet al.2018). Root growth parameters,such as the elongation rate,can be calculated based on our method and used as model inputs to simulate root growth (Morandageet al.2019;Pagèset al.2020),solute transport (Devauet al.2011;Postmaet al.2014;Gerardet al.2017) and water uptake (Javauxet al.2013;Lobetet al.2014;Meunieret al.2020).

The three-dimensional architecture of the field-grown root system can be reconstructed by combining the extracted values from the method with the skeleton of the root system measured in the field (Wuet al.2015;Huiet al.2018). Using reconstructed 3D root models,Wu Qet al.(2018) proposed an improved two-core sampling strategy based on an areaweighted algorithm to optimize soil-coring strategies for field-grown maize.

5.Conclusion

The image-based root phenotyping method developed in this study provides a novel approach with high-throughput and high-accuracy for field research on root morphology and branching features. Using this method,we found that the interspecific advantages for maize mainly occurred within 5 cm from the root base in the nodal roots of the 5th-7th nodes,and that the obvious inhibition of soybean was mostly reflected within 20 cm from the root base. This method could be applied to the 3D reconstruction of field-grown root system architecture to improve inputs to data-driven models that simulate root growth,solute transport and water uptake.

Acknowledgements

We would like to thank Li Shilin and Liu Fusang,Zhu Binglin,Wu Wenfeng and Che Yingpu,from China Agricultural University for their help in conducting field sampling and image scanning. This work was supported by the National Key Research and Development Program of China(2016YFD0300202),the Science and Technology Project of Yunna,China (2017YN07) and the Science and Technology Major Project of Inner Mongolia,China (2019ZD024 and 2020GG0038).

Declaration of competing interest

The authors declare that they have no conflict of interest.

Appendixassociated with this paper was available on http://www.ChinaAgriSci.com/V2/En/appendix.htm