Fine mapping of qPH9, a major quantitative trait locus, responsible for plant height in foxtail millet [Setaria italica (L.) P. Beauv.]

Plant height is vital for crop yield by influencing plant architecture and resistance to lodging. Although lots of quantitative trait loci (QTLs) controlling plant height had been mapped in foxtail millet, their contributions to phenotypic variation were generally small and mapping regions were relatively large, indicating the difficult application in molecular breeding using marker-assisted selection (MAS). In the present paper, a total of 23 QTLs involving in 15 traits were identified via a high-density Bin map containing 3024 Bin markers with an average distance of 0.48 cM through an F2 population. Among them, qPH9, with a large phenotypic variation explained (51.6%) related to plant height, was one of the major QTLs. Furthermore, qPH9 was repeatedly detected in multi-environments under field conditions using two new developed F2 populations from the same F1 plant, and was narrowed down to a smaller interval of 281 kb using 1024 recessive F2 individuals from the same F1 plant. Finally, we found that there was an extremely significant correlation between marker MRI1016 and plant height, and further speculated that Seita.9G088900 and Seita.9G089700 could be key candidates of qPH9. This study laid an important foundation for the cloning of qPH9 and molecular breeding of dwarf varieties via MAS.

RAD-seq Restriction site-associated DNA sequencing FPKM Fragments per kilobase million SNP Single nucleotide polymorphism

InDel Insertion-deletionIntroduction
Plant height (PH), mainly determined by length of stem internode and number of stem node, is one of important crop yield traits by directly or indirectly affecting heading time, plant structure, lodging resistance, etc. Numerous studies have shown that phytohormones, such as gibberellins (GAs) (Sun 2011), brassinosteroids (BRs) (Clouse et al. 1996), auxin (Multani et al. 2003), and strigolactones (SLs) (Jiang et al. 2013), participate in PH regulation. Mutations in plantohormone biosynthesis or signal transduction may lead to phenotypic variation of PH. For instance, rice mutants, semidwarf1 (OsGA20ox2) and dwarf18 (OsGA30ox2) involving in GAs biosynthesis (Spielmeyer et al. 2002;Itoh et al. 2001), Dwarf1 (OsGID1) and Gibberellin Insensitive Dwarf2 (OsGID2) involving in GAs signal transduction (Ueguchi-Tanaka et al. 2005;Sasaki et al. 2003), exhibited diverse variation in PH. Additionally, PH is also controlled by co-located QTLs or pleiotropic genes. In rice, DTH7(Ghd7), encoding a CCT domain protein, played a pleiotropic role on plant height, heading date and grains per panicle (Xue et al. 2008;Gao et al. 2014); Ghd8, encoding a CCAAT boxbinding protein and belonging to the HAP3 subfamily, was thought a pleiotropic locus affecting grain yield, heading date, and plant height (Yan et al. 2011). Foxtail millet is prone to lodging, which may lead to a decrease both in yield and quality (Tian et al. 2010). The utilization of dwarf lines was thought to be a main strategy to improve crop yield by increasing the resistance to lodging (Gooding et al. 2012). However, in foxtail millet, there were few applications of dwarf lines in production for lack of findings on efficient genes or QTLs. Although some corresponding QTLs had been identified, such as qPH1.1 (Wang et al. 2017), qph1, qph4, qph5, qph6, qph7 and qph9 (Zhang et al. 2017), qPH1-1, qPH1-2, qPH5-2, qPH5-3, qPH8 and qPH9 (Wang et al. 2019), and several candidates including SiD1, SiD2, and Seita.1G242300 were cloned and characterized (Zhao et al. 2019;Xue et al. 2016;He et al. 2021), their contributions to phenotypic variation were generally minor or mapping regions were relatively large, indicating the difficult application in molecular breeding using MAS.
In the present study, QTL mapping was carried out using a high-density Bin map constructed through resequencing strategy in an F 2 population derived from a cross of two foxtail millet cultivars with low and high plant height, respectively. A total of 23 QTLs for 15 traits were identified, in which qPH9, a major QTL for PH, had been repeatedly validated in two environments. Furthermore, qPH9 was fine mapped onto an interval of 281 kb through a large recessive F 2 population. Finally, the InDel marker MRI1016 tightly linked to PH was validated through an analysis of variance (ANOVA), and key candidates of qPH9 were identified. Our findings provided valuable information for cloning of qPH9 and molecular breeding of dwarf varieties through MAS in foxtail millet.

Plant materials
To map QTLs for PH, we developed an F 2 population derived from a cross between Henggu12 and Chang-nong35. Henggu12, the female parent, was characterized by lower biomass, earlier heading date, lower plant height, smaller panicle, and multiple tillers (Li et al. 2015); in contrast, Changnong35, the male parent, exhibited higher biomass, later heading date, higher plant height, larger panicle and no tiller (Guo et al. 2008). Henggu12 was crossed with Chang-nong35 to obtain F 1 in 2015, and the hybrid F 1 was selected to self-pollinate to generate F 2 populations: population 1 included 182 plants for preliminary mapping of QTLs for 15 traits in 2017; populations 2 and 3 included 210 plants and 297 plants for verification of qPH9 in Changzhi and Yuanyang, respectively; population 4 was consisted of 5000 plants, in which 1024 recessive individuals were used for fine mapping of qPH9 in 2018.
The F 2 individuals and the biparent were grown in row plot with a row length of 2 m, and resonable planting densities of 30 seeds per row for the F 2 individuals, and 60 seeds per row for the biparent. In 2017, population 1 to map QTLs for 15 traits were sown at Millet 1 3 Research Institute of Shanxi Agricultural University (MRI), Changzhi (113° 08′ E, 36° 18′ N), China, with an average temperature of 28℃/16℃(day/night) in crop season from May to September. In 2018, population 2 and 3 for qPH9 verification were sown at MRI and the Modern Agricultural Science and Technology Experiment Base, Yuanyang (113° 96′ E, 35° 05′ N), China, with an average temperature of 32℃/ 21℃ (day/night) in crop season from May to September, while population 4 were sown at MRI. Additionally, the natural population (Supplementary Table S1) for correlation analysis in 2018, and the biparent for expression analysis in 2021 were planted with the same sowing method to biparent described above.

Phenotypic investigation and statistical analysis
A total of 15 traits were observed in the present study. During maturity, leaf number of the main tiller (LNT), tiller number (TN), and branch number per panicle (BNP) were counted by manual; plant height (PH, cm), main panicle length (MPL, cm), and neck length (NL, cm) were measured by a ruler; main panicle diameter (MPD, cm) and main stem diameter (MSD, cm) were measured using a vernier caliper; flag leaf length (FLL, cm) and flag leaf width (FLW, cm) were measured using a hand-held laser leaf area meter (YMJ-D, Zhejiang Tuopuyunnong technology Co., Ltd). After harvest, main panicle weight (MPW, g), main grain weight (MGW, g), fresh weight per plant (FWP, g), straw weight per plant (SWP, g), and thousand-grain weight (TGW, g) were measured by electronic balance with an accuracy for 0.01 g. For the natural population, 5 plants for every traits were measured, and the mean value was regarded as the corresponding trait phenotype value. All the raw data for 15 traits were recorded using Microsoft Excel 2007. Descriptive analysis (mean, maximum, minimum, skewness, and kutosis) and bivariate correlation analysis (Pearson correlation coefficient and bilateral detection) were conducted using SPSS 17.0 software (IBM, Chicago, IL, USA). Frequency analysis was done using the R software (version 3.6.2).

Resequencing and SNP genotyping
Genomic DNA of leaves was isolated using CTAB method (Chen and Ronald 1999). Resequencing was carried out according to the Illumina protocol.
Genomic DNA was sheared into ~ 350-base pairs (bp) fragments to construct the library. Paired-end reads (150 bp) were sequenced using the HiSeq 2500 system (Illunima, Inc., San Diego, CA, USA). Lowquality reads (quality score < 20e) were filtered out, and raw reads were sorted to each sample according to barcode sequences. After trimming, clean reads from each sample were mapped according to Yugu1 genome sequence (Setaria italica v2.2) using the Burrows-Wheeler Aligner (BWA) software ). The Samtools software ) was used for mark duplicates, and the Genome Analysis ToolKit (GATK) software (McKenna et al. 2010) was used for local realignment and base recalibration. A SNP set was formed by combining GATK and samtools, and SNP calling was analyzed with default parameters. Finally, the SNPs identified between biparent were considered polymorphic, and the polymorphic SNPs with aa × bb pattern were selected for a subsequent analysis.
For the Bin calling, sliding window method (15 SNPs for each window sliding a site) was used to obtain the genotype. Windows containing more than 13 "aa" or "bb" types were genotyped as "aa" or "bb," respectively. Fifteen adjacent SNP intervals with the same genotype across the entire F 2 population were combined into a recombination Bin. The Bin markers, with significant segregation distortion (χ 2 test, P < 0.0001), were excluded from the subsequent Bin map construction. Recombination breakpoints were determined by the physical locations where the genotype of window changed. Briefly, when the entire F 2 population did not have the recombination breakpoints in 100 kb interval, these regions were combined to one bin; while there are more than two recombination breakpoints in 100 kb interval, they were determined as one breakpoint. Additionally, when short heterozygous regions were located in the middle of consecutive "aa" or "bb" genotype regions and these heterozygous regions were shorter than window size, the genotype was manually corrected to match with the adjacent genotype.
Secondly, the modified logarithm of odds (MLOD) scores between markers were calculated to further confirm the robustness of markers for each LGs. Markers with MLOD scores < 5 were filtered prior to ordering. Thirdly, HighMap strategy described by Liu et al (2014) was utilized to order the Bin markers and correct genotyping errors within LGs. The specific process was as follows: (1) Bin markers were selected using spatial sampling. One marker was taken randomly in a priority order of test cross, and markers with a recombination frequency smaller than a given sampling radius were excluded from the marker set; (2) simulated annealing was applied to searching for the best map order. Blocked Gibbs sampling was employed to estimate multipoint recombination frequencies of the parents after the optimal map order of sample markers was obtained; (3) the error correction strategy of SMOOTH was then conducted according to parental contribution of genotypes ( Van et al. 2005), and a k-nearest neighbor algorithm was applied to impute missing genotypes (Huang et al. 2012). Skewed markers were then added into this map by applying a multipoint method of maximum likelihood. Map distances were estimated using the Kosambi mapping function (Kosambi 1943).

QTL analysis and validation of qPH9
QTL analysis was conducted through the Bin map and phenotype (raw data of 15 traits) of 182 F 2 individuals from population 1, while validation of qPH9 was carried out by genotyping of new developed polymorphic InDel markers from the preliminary mapping region of qPH9 and phenotype of F 2 individuals from populations 2 and 3. Interval mapping (IM) analysis with default parameters of Map QTL ® 5 (Ooijen 2004) was conducted to identify initial candidate QTLs and validate qPH9. The threshold corresponding to 0.995 confidence was firstly considered with a threshold set of 1000 times permutation test; when no QTL was detected under the threshold, that corresponding to 0.95 confidence was considered, and otherwise that corresponding to 0.90 confidence was considered.

Extraction of total RNA and reverse transcription
Total RNA was extracted from the stem-top-second of the biparent at booting stage using RNAiso Plus (TaKaRa Bio Inc., Shiga, Japan) according to the manufacturer's instruction. The amount and quality of the total RNA were checked through electrophoresis in 1.0% agarose gel, and the concentration of RNA was measured by eppendorf biophotometer (Hamburg, Germany). Equal amount of 2.0 μg total RNA was reversely transcribed to cDNA using PrimeScript II 1st strand cDNA synthesis kit (TaKaRa Bio Inc., Shiga, Japan) with oligo dT primers.
qRT-PCR was conducted in a 10.0 μL reaction volume containing 1.0 μL template cDNA, 5.0 μL TB Green Premix Ex Tap II (TaKaRa Bio Inc., Shiga, Japan), 2.0 μL primer (2 μmol L −1 ), and 2.0 μL ddH 2 O. The reaction was cycled as follows: 95℃ for 3 min, followed by 40 cycles of 95℃ for 10 s, 58℃ for 30 s, 72℃ for 30 s, followed by a melt curve at 65-95℃ by increments of 0.5℃/s. Differences between the cycle threshold (Ct) values of target genes (Seita.9G088900 and Seita.9G089700) and Siactin-7 were calculated as ΔCt = Ct target -Ct actin . Expression level of target gene relative to actin was determined as 2 −ΔCt , and the average value of 2 −ΔCt was used to determine the expression of target gene. Gene expression data were analyzed with ANOVA using the SPSS 17.0 software (IBM, Chicago, IL, USA).

Phenotypic and correlation analysis
In this study, a total of 15 traits were evaluated (Supplementary Table S2) among 182 F 2 individuals and biparent. The phenotypic variations for each trait had a wide range with a large variation coefficient, and the absolute value of the skewness and kurtosis of most traits was < 1, indicating that these traits had great potentials for genetic improvement. Additionally, apart from NL with a normal distribution, the other traits exhibited a dramatic distribution in F 2 population: PH, LNT, MPL, MPD, MSD, FLL, FLW, TGW, and BNP had a distinct bimodal distribution, and FWP, SWP, MPW, MGW, and TN deviated from a normal distribution (Fig. 1).
The Pearson correlation coefficient with bilateral detection was performed across 15 traits (Supplementary Table S3). TN with the other traits showed negative correlation; NL with PH, MPL, MPD, MSD, BNP, FLL, FLW, FWP, SWP, MPW, MGW, and TGW showed negative correlation; positive correlation was observed among the other traits. Even some traits exhibited extremely significant correlation. For instance, the correlation coefficients between PH and MPL, MSD, LNT, FLL and MSD were 0.960, 0.941, 0.928, and 0.901, respectively; LNT and MPL, MSD and FLL were 0.937, 0.925, and 0.918, respectively; MPL and FLL were 0.943; MPD and MSD were 0.916, indicating these traits might be affected by a pleiotropic gene or closely linked loci, respectively.
Construction of the high-density genetic map Biparent and 182 F 2 individuals from population 1 were sequenced in the present study. Totally, we obtained 8.08 Gbp clean data for Henggu12 with an average coverage of 17 × , 7.38 Gbp for Changnong35 with an average coverage of 16 × , and 278.38 Gbp for offspring with an average coverage of 3.23 × . The genetic map was constructed using 3024 Bin markers originated from 651,047 SNPs, and the distance of the map was 1457.40 cM with an average distance of 0.48 cM between adjacent Bin markers (Table 1; Supplementary Fig. S1; Supplementary Table S4). Interestingly, compared with other chromosomes, Chr. IX shared the largest number of markers, the maximum length of total distance, the minimum length of average distance, and the minimum gap. Additionally, colinearity analysis showed that the genetic map had a good collinearity with the reference genome (Fig. S2).

QTL analysis for 15 traits
MapQTL was carried out to detect QTLs controlling 15 traits. A total of 23 QTLs were identified on Chr. IX (Table 2). Collectively, only one QTL was identified for PH, LNT, MSD, FLL, FLW, FWP, SWP, MGW, TGW, and BNP, with a PVE from 31.2 to 54.5%, an additive effect from − 0.22 to − 56.92, and a dominant effect from 0.38 to 116.76, respectively; two QTLs were identified for MPL, MPD, MPW, and NL, respectively, in which qMPL9-2 for MPL and qMPD9-2 for MPD were two major effect QTLs and accounted for 47.9% and 45.4% of the PVE, respectively; for TN, 5 QTLs (qTN9-1, qTN9-2, qTN9-3, qTN9-4, and qTN9-5) with a PVE from 12.2 to 20.7%, an additive effect from 0.42 to 0.59, and a dominant effect from − 0.44 to − 0.65 were detected, respectively. Notably, qPH9 with a large PVE (51.6%), a major QTL for PH, was observed and shared the same position with qMSD9 and qFLL9, in charge of MSD and FLL, respectively.

Validation of qPH9
To verify qPH9, QTL analysis was conducted using the F 2 individuals (population 2: 210 plants, population 3: 297 plants) from the same F 1 plant in two environments (Changzhi and Yuanyang, China). A total of 12 InDel markers (Supplementary Table S5), showing polymorphisms between two parental lines in the candidate region of qPH9, were used to validate qPH9. The results showed that the distinct bimodal distribution for PH was repeatedly observed in two environments, demonstrating that PH was controlled by a major QTL or gene, and that qPH9 was repeatedly detected and shared the same chromosome position in two environments, although there were some differences in LOD score and PVE, which were 66.69 and 79.0% in Changzhi ( Fig. 2a

Fine mapping and candidate analysis of qPH9
To narrow down the mapping interval of qPH9, 13 polymorphism InDel markers (Supplementary  Table S5) were developed from its preliminary mapping physical region (Fig. 3a), and 1024 recessive F 2 individuals (dwarf plants) from population 4 were screened for recombinants. The results showed that qPH9 was located on the interval of 281 kb between MRI1014 (5,314,541) and MRI018 (5,595,508) (Fig. 3b). Furthermore, ANOVA analysis revealed that there was an extremely significant correlation between MRI1016 and PH not only in F 2 population in two environments, but also in the natural population in Changzhi (Supplementary Table S8).

Discussion
As cost of genome sequencing became decreased, QTL mapping for complex traits, and also for the quality-quantity traits, has become extremely effective by resequencing strategy. Previously, in Capsella rubella, CrFLC controlling flowering time, as a quality-quantity gene was identified and cloned using two distinct F 2 population (Yang et al. 2018), and qCOK2 with a PVE of 31.9% for color of kernel was detected and fine mapped through an F 2 population in maize . In this study, we preliminarily mapped, repeatedly validated and fine mapped the qPH9 using four different F 2 populations Fig. 1 Distribution of 15 traits in the F 2 population derived from a cross between Henggu12 and Changnong35 ◂ from the same F 1 plant derived from a cross between Henggu12 and Changnong35. Inspiringly, qPH9 could be repeatedly detected in multi-environments, and finally it was fine mapped onto a smaller interval of 281 kb, indicating that it was effective to fine map QTLs for quality-quantity trait using several different F 2 mapping populations from an F 1 individual via resequencing strategy. These results provided a further insight on fine mapping of QTLs associated with quality-quantity trait. 1 3 In the present study, extreme QTL cluster for all QTLs except qNL9-1 and qNL9-2 was observed. QTL cluster was also found and thought to be associated with ecological significance in wild annuals of rice (Onishi et al., 2007). In this study, biparent had distinct features in multi-traits associated with ecological significance including heading date, sensitivity to light and temperature, and plant height. Hengu12 could be grown across China and the performance could be almost normal due to its shorter heading date and insensitive to light and temperature (Li et al. 2015), but Changnong35 was only planted in north China, indicating that ecological adaptability might led to the QTL cluster on Chr. IX in foxtail millet. Also, we found all QTLs for 15 traits were mapped onto Chr. IX and none was present on the other chromosomes. This phenomenon might be resulted from those causes as follows: (1) extremely significant correlation: this was found almost among all traits (Supplementary Table S3), indicating major QTLs responsible for those traits might be tightly linked together each other; (2) quality-quantity trait: almost all traits exhibited a distinct bimodal distribution (Fig. 1), indicating that these traits could be of quality-quantity in this population. Due to the nature of quality-quantity traits, major QTLs were more easily detected and the minor QTLs were difficultly detected in an F 2 mapping population; (3) ecological adaptability illustrated above. To this end, it was necessary to develop NIL populations to further map minor QTLs.
QTL co-location is a widespread phenomenon in plant genome. Co-located QTL was usually used to assess the selective genetic effects, or to understand pathways of genetic activity for breeders and geneticists (Hill and Zhang 2012). Generally, co-located QTLs were caused by pleiotropy, such as Q gene responsible for free-threshing, spike rachis fragility, plant height, spike shape, and heading date (Simons et al. 2006), and Rht-B1 gene responsible for plant height, grain size, and weight (Guan et al. 2020) in wheat, or close linkage of different genes affecting multiple traits such as Dw2 and Ma1 in sorghum in charge of plant height and flowering time, respectively (Klein et al. 2008). In foxtail millet, previous studies showed that this phenomenon also existed (Zhang et al. 2017;Wang et al. 2019;Liu et al. 2020;Mauro-Herrera and Doust 2016;Fang et al. 2016). In the present study, QTLs controlling PH, MSD, and FLL with large effects were also co-located onto the same region on Chr. IX, which was consistent with significant positive phenotypic correlation between these traits, indicating that qPH9 might be a pleiotropic QTL.
In this study, qPH9 was fine mapped onto the interval of 281 kb (5,314,595,508) between InDel marker MRI1014 and MRI1018 on Chr. IX. Previous studies showed that 9 QTLs for PH on Chr. IX, consisting of H9a: 3,452,395-9,964,754 (Mauro-Herrera and Doust 2016), qph9: 1,559,957-1,859,997 (Zhang et al. 2017), qPH9: 8,207,836-8,386,583 (Wang et al. 2019 Fine mapping of qPH9. a qPH9 was preliminarily mapped on Chr. IX with a physical interval of 9.2 Mb. Scale of Y-axis on the left showed LOD value, and the right showed PVE, scale of X-axis showed the chromosomes; The gray line represented the threshold value of QTL mapping, the blue line represented LOD value, and the red line represented PVE. b qPH9 was narrowed down to an interval of 281 kb between MRI1014 and MRI1018.The marker name, physical position, and number of recombinants were indicated 42,767,054-43,133,505, qPH9.5: 43,430,778-43,849,397, and qPH9.6: 46,893,323-47,598,99 (He et al. 2021), had been detected in foxtail millet. Distinctly, except for H9a, physical regions of the other QTLs were different from ours (5,314,595,508). Notably, H9a had been repeatedly detected not only at different growth stages but also in various environments (Mauro-Herrera and Doust 2016), in spite of the physical region was distinctly larger than qPH9 in this study. Interestingly, qPH9 was also repeatedly detected in multienvironments, and had the similar PVE from 76.4 to 79% in different environments, indicating that qPH9 was a stable and high hereditary locus. Indubitably, there have been one or several stably major QTLs on Chr. IX responsible for PH, and qPH9 in the present study was the finest one mapped so far.
To analyze the candidate genes for qPH9, we searched for the putative functions of genes across the physical region of qPH9 by Phytozome database. We found that there were 36 annotated genes within the physical region of qPH9, of which eight genes, including Seita.9G088400, Seita.9G088500, Seita.9G088900, Seita.9G089000, Seita.9G090100, Seita.9G091200, Seita.9G091400, and Seita.9G092100 exhibited higher expression level in stem either in Xiaomi or in JG21, and one, Seita.9G089700, exhibited higher expression both in Xiaomi and in JG21 (Supplementary Table S9). Bioinformatics analysis showed that Seita.9G088900 was homologous to OsMADS14 and AtAGL7 regulating the development of inflorescence meristem and flowering in rice and Arabidopsis, respectively (Jeon et al. 2000;Han et al. 2014), and Seita.9G089700 was homologous to ZmPHYC involving in flowering time and plant height in maize (Li et al. 2020). Moreover, analysis of recombinant plants near the markers MRI1014 (5,314,541) and MRI1016 (5,490,668) showed that MRI1016 has the smaller number of recombinant plants. Thus, we speculated that Seita.9G089700 and Seita.9G088900 could be better candidates for qPH9. Furthermore, we analyzed their expression level and sequences between Henggu12 and Changnong35. qRT-PCR analysis showed that, for Seita.9G089700, the relative expression in Henggu12 was 26 times lower than that in Chang-nong35, and for Seita.9G088900, the relative expression in Henggu12 was 2.8 times higher than that in Changnong35 (Fig. S3), indicating that there were of significant expression differences of Seita.9G088900 and Seita.9G089700 in stems between Henggu12 and Changnong35. In addition, we investigated the genomic sequence and complete coding sequence of Seita.9G088900 and Seita.9G089700 between Henggu12 and Changnong35. The results showed that there was no difference not only on genomic sequence but also on complete coding sequence for Seita.9G088900 between Henggu12 and Chang-nong35. Meanwhile, for Seita.9G089700, the complete coding sequence had a 28-bp Deletion and a 6-bp Insertion (+ 2299-bp to + 2327-bp downstream from ATG) in the second exon in Henggu12 compared to Changnong35 (Fig. S4), and furthermore, we obtained its whole genomic sequence in Chang-nong35 but only incomplete sequence in Henggu12, indicating that a complex genomic structural variation would be likely to occur in Henggu12, which remained to further verify.

Conclusion
In the present study, we identified a stable major locus (qPH9) controlling plant height in multi-environments, and further fine mapped it onto an interval of 281 kb. Moreover, we developed an InDel marker MRI1016 tightly linked to PH and found that Seita.9G088900 and Seita.9G089700 might be key candidates of qPH9. These results were conducive to cloning of qPH9 and to molecular breeding of dwarf varieties via MAS in foxtail millet.
Author contribution JW conceived and supervised the complete study; XFD and ZLW performed the experiments and carried out the bioinformatics work; EHG, SCL, KNH, YXL, and LYZ were responsible for the field trial; XFD and JW wrote and revised the manuscript. All authors have read and approved the final manuscript.