Hui You,Sundus Zafar,Fan Zhang,Shuanging Zhu,Kai Chen,Congong Shen,Xiuqin Zhao,Wenzhong Zhang,Jianlong Xu,,d,*

a Rice Research Institute,Shenyang Agricultural University,Shenyang 110866,Liaoning,China

b The National Key Facility for Crop Gene Resources and Genetic Improvement,Institute of Crop Sciences,Chinese Academy of Agricultural Sciences,Beijing 100081,China

c Shenzhen Branch,Guangdong Laboratory for Lingnan Modern Agriculture,Agricultural Genomics Institute at Shenzhen,Chinese Academy of Agricultural Sciences,Shenzhen 518120,Guangdong,China

d Hainan Yazhou Bay Seed Lab/National Nanfan Research Institute(Sanya),Chinese Academy of Agricultural Sciences,Sanya 572024,Hainan,China

Keywords:Milling and appearance quality Quantitative trait locus/loci(QTL)Heterosis Hybrid rice Pyramiding breeding

ABSTRACT Development of hybrid rice with high yield and grain quality is a goal of rice breeding.To investigate the genetic mechanism of heterosis for rice milling and appearance quality in indica/xian rice,QTL mapping was conducted using 1061 recombinant inbred lines(RILs)derived from a cross of the xian rice cultivars Quan9311B(Q9311B)and Wu-shan-si-miao(WSSM),and a backcross F1(BC1F1)population developed by crossing the RILs with Quan9311A(Q9311A),combined with phenotyping in two environments.The F1 hybrid(Q9311A×WSSM)showed various degrees of heterosis for milling and appearance quality.A total of 142 main-effect QTL(M-QTL)and 407 pairs of epistatic QTL(E-QTL)were identified for five milling and appearance quality traits and grain yield per plant(GYP)in the RIL,BC1F1 and mid-parental heterosis(HMP)populations.Differential detection of QTL in three populations revealed that most additive loci detected in the RILs did not show heterotic effects,but some of them did contribute to BC1F1 trait performance.Unlike heterosis of GYP,single-locus overdominance and epistasis were the main contributors to heterosis for milling and appearance quality.Epistasis contributed more to the heterosis for milling quality than to that for appearance quality.Three(four)QTL regions harboring opposite(consistent)directions of favorable allele effects for GYP and grain quality were identified,indicating the presence of partial genetic overlaps between GYP and grain quality.Three strategies are proposed to develop hybrid rice with high yield and good grain quality:1)pyramiding favorable alleles with consistent directions of gene effects for GYP and grain quality at the M-QTL on different chromosomes;2)introgressing favorable alleles for GYP and grain quality into the parents and then pyramiding and fixing these additive effects in hybrids;and 3)pyramiding overdominant and dominant loci and minimizing or eliminating underdominant loci from the parents.

1.Introduction

Rice(Oryza sativa)is the staple food of more than half of the world’s population.With improvement in living standards,the main goal of rice breeding has gradually expanded from increasing yield potential to including grain quality[1].

Rice quality varies from place to place and from person to person,and can be assessed by many criteria.It can be divided into four components:milling,appearance,cooking,and nutritional quality.Milled rice appearance strongly influences market value[2].Milling quality includes brown rice rate,milled rice rate,and head milled rice rate.Appearance quality is usually evaluated by the degree of endosperm chalkiness and proportion of chalky grain.Because most rice quality traits are quantitatively inherited and strongly influenced by environment[3],it is challenging to improve quality by conventional breeding[4].

The milling and appearance quality traits of rice are controlled by quantitative trait loci(QTL)[5,6].Many genes controlling rice quality traits have been cloned in recent decades,including GS3,GW2,GW5,GS5,GW6a,and GL7[7-12],which regulate grain size and shape.Only a few QTL,including GIF1,Chalk5,cyPPDK,OsRab5a,PDIL1-1,and SSG4[13-18],that affect grain appearance by regulating the endosperm filling and accumulation of starch have been cloned.Starch,the principal storage substance in the endosperm,is affected by the expression of starch synthesis genes and is thus closely associated with chalkiness[19].Some QTL affecting appearance quality have been fine-mapped.qPGWC-7 was mapped to a 44-kb DNA interval on chromosome 7 containing 13 predicted genes[20],and qPGWC-8 was narrowed to a 142-kb region on chromosome 8 between InDel markers 8G-7 and 8G-9[21].No gene affecting milling quality has yet been cloned.But recently,qBRR-10 controlling brown rice rate was localized to a 39.5-kb region on chromosome 10 containing two candidate genes[22].

Heterosis is a common phenomenon in which an F1hybrid performs better than either inbred parent[23].In crop genetic research,the mechanism of heterosis has always been a hot topic,and several hypotheses have been proposed:the dominance hypothesis,which proposes the masking of deleterious recessive parental alleles in the hybrid;the overdominance hypothesis,which attributes heterosis to the superiority of heterozygotes over parental homozygotes at individual loci;and the epistatic hypothesis,which postulates the contribution of positive epistatic interactions between non-allelic genes[24,25].In recent years,with the development of DNA markers and genome sequencing technology,researchers have acquired a better understanding of the genetic mechanisms of heterosis of rice yield traits using F2,BCF1,and testcross F1populations derived from distantly related parents or parents of well-known hybrid cultivars[26-29].However,there are no reports explaining the genetic mechanisms responsible for heterosis of rice grain quality traits.

Wu-shan-si-miao(WSSM)and 9311 are elite inbred xian rice varieties with high yield,good grain quality and wide adaptability,and have become backbone parents in two-line and three-line hybrid rice breeding in the past 20 years in China.It is estimated that 31 two-line hybrid combinations and 17 three-line hybrid combinations have been developed using WSSM as a direct paternal line,and six two-line hybrid combinations have been developed using 9311 as a direct paternal line and dozens of two-line hybrid combinations developed using the paternal lines derived from 9311 as a main lineage so far in China(https://www.ricedata.cn).Quan-you-si-miao(QYSM),produced by Q9311A(derived from 9311,and near-isogenic to Q9311B)and WSSM,a superior xian hybrid cultivar with high yield,good grain quality,and wide adaptability,has become one of the predominant three-line xian hybrid cultivars widely grown in the whole south China rice region covering the upper,middle,and lower areas of the Yangtze River.In view of QYSM’s high yield and good grain quality,it was desirable to dissect the genetic mechanisms underlying the heterosis of rice grain quality and yield and reveal genetic relationships between them.

The objective of this study was to dissect the heterotic genetic basis of milling and appearance quality in a high-yield and goodquality rice hybrid,using deep next-generation sequencing and phenotyping of a set of recombinant inbred lines(RILs)and backcrosses developed from them.

2.Materials and methods

2.1.Plant materials

QTL mapping was performed using 1061 F2:8RILs derived from the cross of Quan9311B(Q9311B)and Wu-shan-si-miao(WSSM)by single-seed descent,and the corresponding 1061 BC1F1progeny developed by crossing the RILs with Quan9311A(Q9311A)(Fig.S1).Q9311A is a sterile line near-isogenic to Q9311B.Thus,Q9311B was used to cross with WSSM for maintaining normal fertility in the RILs and Q9311A was used to cross with the RILs for easy development of a testcross population.Other materials included an F1hybrid,Quan-you-si-miao(QYSM),which was developed from Q9311A and WSSM.Q9311A and WSSM are considered backbone parents in hybrid rice breeding programs in China for their superior general combining ability and grain quality.

2.2.Field experiment and trait evaluation

The RIL and BC1F1populations were generated in the experimental field of Chinese Academy of Agricultural Sciences(18.3°N,109.3°E),Sanya,Hainan province,in the 2018 winter season.The RILs,BC1F1s,parents(Q9311B and WSSM),and QYSM were sown in the seedling nursery on May 17 and July 12,2019 in Hefei(32.3°N,117.6°E)of Anhui province and Nanning(22.1°N,107.5°E)of Guangxi province,respectively.Owing to its infertility,Q9311A was not included in the field experiment for trait evaluation,and its trait values were evaluated using the near-isogenic line Q9311B.The 25-day-old seedlings were transplanted into a four-row plot with six plants in each row at 13×30 cm spacing for two replications in a randomized complete block design.Field water,weed,and fertilizer management followed local practice.

At the mature stage,seeds from six uniform plants in the middle of each plot were harvested and air-dried for three months in drying houses to a 13.5% grain moisture content.After removal of unfilled grains,the grain was weighed to calculate grain yield per plant(GYP,g).Samples of 30 g whole rice seeds were milled into brown rice with a rice huller JLGJ-45(Taizhou Food Instrument,Taizhou,Zhejiang,China)and then milled into polished rice with a rice milling machine JNMJ 6(Taizhou Food Instrument).Three milling quality traits were measured:brown rice rate(BRR,%),milled rice rate(MRR,%),and head milled rice rate(HMRR,%).Two appearance quality traits:proportion of grain with chalkiness(PGWC,%),and degree of endosperm chalkiness(DEC,%),were measured with a scanning machine(SC-E,Wanshen Technology Company,Hangzhou,Zhejiang,China).All measurements were performed with two replicate samples and the mean was used for data analyses.

2.3.Genotyping

Genomic DNA for SNP genotyping was isolated from approximately 100-mg fresh leaf samples of 5-week-old seedlings for each of the 1061 RILs,Q9311B,and WSSM using a modified cetyltrimethylammonium bromide method[30].Whole-genome shotgun sequencing of the parental lines(Q9311B and WSSM)and RILs generated 150 bp paired-end reads consisting of～2.4 Tb of sequence.The RILs were sequenced with～10×genome coverage,and the two parental lines were sequenced with～40×genome coverage using the Illumina sequencing platform NovaSeq 6000(CapitalBio Technology Inc.,Beijing,China).The sequence reads of the parental lines and RILs were aligned against the Nipponbare RefSeq(IRGSP 1.0)rice reference genome sequence[31]with BWA 0.7.1[32]using default parameters,and PCR duplicates were removed by the MarkDuplicates module in Picard Tools 1.119(https://broadinstitute.github.io/picard/).The raw reads were also realigned to identify highly polymorphic regions using the IndelRealigner function in GenomeAnalysisTK 3.4.0[33].Sequence variants between the parental lines were called with UnifiedGenotyper in GenomeAnalysisTK.Only uniquely mapped reads were used for subsequent single-nucleotide polymorphism(SNP)calling.Genotype calling of each RIL was performed based on the SNP alleles between the parents.A total of 741,928 homozygous SNPs polymorphic between Q9311B and WSSM were used for SNP calling in the 1061 RILs.The SNPs were further filtered by removal of those with minor allele frequencies＜0.01 or missing rates＞0.20,leaving 156,373 high-quality and evenly distributed SNPs,located within the genic regions of 15,043 genes based on the gene annotation of the Nipponbare RefSeq from RAP-DB(released on June 26,2019)[34].A set of 1196 genes was removed because one of the two parental alleles was rare or showed a high frequency of the heterozygote in the RIL population.Then,an allelic genotyping dataset containing two genotypes(Q9311B type and WSSM type)at each of 13,847 genes was generated to be used for bin map construction for the 1061 RILs.

2.4.Bin map construction and QTL mapping

A bin map is a type of genetic map constructed using bin(binning of redundant markers)genotypes,which can save running time of linkage and QTL mapping by greatly reducing number of redundant markers.The bin map of the RIL population was constructed based on the two parental genotypes at 13,847 genes.To remove redundancy,only a single gene was retained to represent each bin,either one gene with a minimum missing rate or a random gene when the missing rates were equal.Genes with unique genotypes(those that did not belong to a redundant bin)were excluded from bin-map construction.As a result,855 bins were identified with mean size～440 kb and used for construction of bin map using the BIN function in IciMapping QTL 4.2[35].Allelic genotyping datasets containing two bin genotypes(Q9311B type and WSSM type)at each of 855 bins for the 1061 RILs,and the genotype for each BC1F1hybrid were deduced from the gene-based allelic genotypes of its parental RIL and Q9311B.Specifically,for each of the 13,847 genes,if the parents(RIL and Q9311B)have the same allelic genotype,their BC1F1hybrid should be the homozygous genotype of Q9311B;if the parents had different homozygous genotypes,the allelic genotype of their BC1F1hybrid was deduced as the heterozygote;and if RIL parent had the heterozygous genotype,the BC1F1was treated as missing.The bin genotypes were used for QTL mapping.

QTL mapping for six traits(three milling quality traits,two appearance quality traits,and GYP)was performed separately for the RIL,BC1F1,and the corresponding HMPpopulations evaluated in both environments.The HMPpopulation was produced by estimation of mid-parental heterosis for each BC1F1as follows.The mid-parental heterosis value,HMP(%)=(F1-MP)/MP×100,where F1is the trait value of a BC1F1and MP is the mean value of the corresponding paternal RIL and Q9311B[36].All traits were analyzed using the BIP(biparental populations)function in IciMapping QTL.The population settings of RILs and BC1F1s were F1RIL and P2BC1F1,respectively.The inclusive composite interval mapping of additive(ICIM-ADD)QTL method was performed to identify main-effect QTL(M-QTL)using the default settings in which P values for entering a variable(PIN)were set at 0.001 and the scanning step was set at 1.0 cM.The inclusive composite interval mapping of the digenic epistatic(ICIM-EPI)QTL method detected possible digenic epistatic QTL(E-QTL)using the default settings.The corresponding scan step and PIN for E-QTL mapping were set at 5 cM and 0.0001,respectively.The(logarithm of odds)LOD threshold values 3.0 and 5.0 were determined according to 1000 permutations at a 95% confidence level to detect M-QTL and E-QTL,respectively[37].The physical positions of a QTL were retrieved based on the left and right markers of the detected interval.

2.5.Inference of QTL gene action

The gene effects of heterosis can be dissected through QTL mapping using segregation populations such as the RILs derived from QYSM and the BC1F1s between Q9311A and the RILs,where the additive effect(a)of M-QTL can be detected in the RILs while the sum of additive effect and dominance effect(a+d)is detected in BC1F1s.When HMPs are used for QTL analysis of heterosis,the inherited effect should be a dominance effect(d).The detected M-QTL can be divided into four types in RILs and BC1F1s following Mei et al.[38]and Kim et al.[39].It was called the additive effect when a QTL was detected only in the RILs or BC1F1s(A).When QTL were simultaneously detected in the RILs and BC1F1s,or in the HMPof BC1F1s,QTL with|2d/(a+d)|＜1 or|2a/(a+d)|＞1 were designated as incompletely dominant QTL(ID).QTL with|2d/(a+d)|＞1 or|2a/(a+d)|＜1,or those detected only in HMPdatasets were designated as overdominant QTL(OD),and QTL showing negative values were defined as underdominant QTL(UD).Only ID-type,OD-type,and UD-type M-QTL were used for subsequent heterosis analysis of the detected M-QTL.

2.6.Identification of candidate genes

QTL with heterotic effects simultaneously detected in both environments for milling and appearance quality traits were used for candidate gene analysis.All genes in each QTL interval were identified for further candidate-gene analysis using haplotype analysis of 732 accessions(386 xian,219 japonica/geng,46 of intermediate type,67 aus/boro,and 14 basmati/sadri accessions)from the 3K Rice Genome Project(3K RGP)(Table S1).The five milling and appearance quality as well as GYP phenotypic data were collected in Sanya.All available high-quality SNPs with a minor-allele frequency of more than 0.05 and/or a missing rate of less than 20% located in these genes were retrieved from 18 M SNP data generated from 3K RGP in the Rice SNP-Seek Database[40].Haplotype analysis was performed for each of the candidate genes in each QTL region using all nonsynonymous SNPs located inside the gene CDS region.Haplotypes shared by more than 10 accessions were compared.For a single candidate gene,the haplotype with the best phenotype of the corresponding trait was assigned as a favorable haplotype.

2.7.Statistical analysis

Differences in the mean phenotypic values among the haplotypes were evaluated by one-way(ANOVA)using the agricolae package in R[41].Phenotypic differences among the check parents and the relative hybrids were assessed by Duncan’s multiple comparison test,and among RILs and the BC1F1s by Student’s t-test.HMPwas tested with a Student’s t-test based on the contrast between the F1hybrid mean and the mean performance of the corresponding parental lines[36].Phenotypic correlation analyses of the six traits were computed using the corrplot package in R.Broad-sense heritabilities(H2)were calculated using the AOV module implemented in IciMapping QTL.

3.Results

3.1.Trait performance and correlation

The field performance of all experimental materials is shown in Table 1.Q9311B showed lower MRR and HMRR than WSSM in HF.Heterosis was detected in DEC in the two environments and PGWC in HF,and various degrees of heterosis were also found for milling quality and GYP in the two environments and PGWC in NN.The RIL and BC1F1populations showed wide segregation in the six traits,especially for HMRR,PGWC,DEC,and GYP in the two environments.The BC1F1population showed lower heterosis for milling and appearance quality than QYSM,having HMPof 0.26(0.81)%for BRR,5.62(0.19)% for MRR,2.97(7.43)% for HMRR,-14.69(-1.38)% for PGWC,-16.41(-9.45)% for DEC,and-4.46(-6.05)% for GYP in HF(NN),indicating that hybrid breakdown occurred in the RILs both in GYP and appearance quality due to inbreeding depression,as shown by the finding that PGWC,DEC,and GYP of the RILs accounted for 137.5(157.6)%,145.4(157.4)%,and 88.3(96.3)% of QYSM in HF(NN).The broad-sense heritabilities(H2)of RILs(BC1F1s)were 0.45(0.43)for BRR,0.41(0.42)for MRR,0.51(0.56)for HMRR,0.69(0.71)for PGWC,0.68(0.70)for DEC,0.24(0.48)for GYP,respectively(Table 1).Thus,quality and yield traits were strongly affected by environment,especially the milling quality and yield traits.

There were respectively 21(2.0%),28(2.6%),21(2.0%),16(1.5%),26(2.5%)and 30(2.8%)BC1F1lines having significantly higher BRR,MRR,HMRR,and GYP and lower PGWC and DEC than QYSM in HF,and respectively 19(1.8%),14(1.3%),15(1.4%),18(1.7%),15(1.4%)and 21(2.0%)BC1F1lines having significantly higher BRR,MRR,HMRR,GYP,and lower PGWC and DEC than QYSM in NN.In both environments,there were 14 BC1F1lines with significantly higher quality traits than QYSM and similar GYP(Table S2).

As expected,significantly positive correlations were found between BRR and MRR,MRR and HMRR,and DEC and PGWC in RILs and BC1F1s under the two environments(Table 2).Significant positive correlations were found between MRR and PGWC,and MRR and DEC in RILs under the two environments.Similarly,significant positive correlations were found between PGWC and GYP,and DEC and GYP,while negative correlations were found between HMRR and PGWC,HMRR and DEC,and HMRR and GYP in BC1F1s under the two environments,indicating that some correlations were associated with genetic background.Some correlations varied with environment:negative correlations between BRR and PGWC,and BRR and DEC,and positive correlations between PGWC and GYP,and DEC and GYP were detected in RILs only in HF and not in NN;and negative correlations between MRR and PGWC,and MRR and DEC were detected in BC1F1s only in HF while positive correlations between BRR and GYP,and MRR and GYP were detected in BC1F1s only in NN.HMRR,PGWC,and DEC showed significantly positive correlations between the two environments in both RILs and BC1F1s,and the same was true for GYP in BC1F1s,indicating that HMRR,PGWC,and DEC showed consistent trends in performance in the two environments and two backgrounds.

3.2.Identification of QTL affecting grain quality and yield traits

The 855-bin linkage map constructed for the RIL population spanned 2268 cM,with chromosome lengths varying from 106 cM(chromosome 7)to 334 cM(chromosome 2)(Table S3).Using genotypic datasets at 855 bins for the 1061 RILs,BC1F1s and HMPs,142 M-QTL for grain milling and appearance quality and GYP of the lines(RILs and BC1F1s)and corresponding HMPs were identified in HF and NN.Of the 142,30 were identified in both environments(Fig.1;Table S4).

For three milling quality traits,eight(nine)M-QTL for BRR were identified in HF(NN),explaining 23.2(3.4)%,8.8(17.5)%,and 3.6(18.1)% of total phenotypic variance in the RILs,BC1F1s,and HMPs,respectively(Table S4).These 17 M-QTL included 8 additive,4 overdominance-type(OD-type),and 5 underdominance-type(UD-type)(Fig.2).qBRR6.1 was detected in RILs in both environments,and qBRR7 was detected in BC1F1s in HF and in HMPs in NN with different gene actions.For MRR,eight(ten)M-QTL were identified in HF(NN),explaining 11.2(8.0)%,18.4(14.0)%,and 13.9(10.1)% of total phenotypic variance in the RILs,BC1F1s,and HMPs,respectively(Table S4).These 18 M-QTL included 12 additive,4 OD-type,and 2 UD-type(Fig.2).qMRR2 was detected in BC1F1s in both HF and NN,and qMRR6.3 was detected in RILs in both HF and NN.A total of 12(14)M-QTL for HMRR were identified in HF(NN),explaining 8.1(12.5)%,29.3(29.0)%,and 25.0(17.8)% of total phenotypic variance in the RILs,BC1F1s,and HMPs,respectively(Table S4).These 26 M-QTL included 16 additive,2 OD-type,and 8 UD-type(Fig.2).Six M-QTL were identified in both environments,including two(qHMRR1.3 and qHMRR7.2)in RILs,two(qHMRR1.1 and qHMRR12.1)in BC1F1s,and two(qHMRR1.2 and qHMRR12.2)in HMPs.

For two appearance quality traits,23(16)M-QTL for PGWC were identified in HF(NN),explaining 19.1(32.5)%,50.7(44.1)%,and 33.8(19.7)% of total phenotypic variance in the RILs,BC1F1s,and HMPs,respectively(Table S4).These 39 M-QTL included 19 additive,9 OD-type,10 UD-type,and 1 incomplete-dominance type(ID-type)(Fig.2).Six M-QTL were simultaneously detected in the two environments,including two(qPGWC2.4 and qPGWC5.3)in RILs,two(qPGWC1.2 and qPGWC5.1)in BC1F1s,and two(qPGWC8 and qPGWC10.4)in HMPs.qPGWC6.3 was detected in BC1F1s in HF and HMPs in NN with differing gene actions.qPGWC6.5 and qPGWC6.6 were separate with opposite additive effects in the same region of 25,182,993-26,442,731 bp on chromosome 6 in NN.For DEC,15(18)M-QTL were identified in HF(NN),explaining 33.9(28.7)%,54.8(47.5)%,and 19.1(11.4)% of the total phenotypic variances in the RILs,BC1F1s,and HMPs,respectively(Table S4).qDEC5.3 and qDEC10.1 showed largest phenotypic variations of respectively 40.0% in BC1F1s in HF and 24.7% in RILs in HF.These 33 M-QTL included 24 additive,4 OD-type,4 UD-type,and 1 ID-type(Fig.2).Two M-QTL(qDEC2.3 and qDEC8.3)in RILs,three(qDEC1.2,qDEC5.1,and qDEC5.3)in BC1F1s were all detected in both HF and NN.qDEC6.5 was detected in BC1F1s in NN and in HMPs in HF,and qDEC10.1 was detected in RILs in HF and in BC1F1s and HMPs in NN.

For GYP,19(20)M-QTL were identified in HF(NN),explaining 6.8(15.3)%,42.2(44.2)%,and 37.4(47.1)%of total phenotypic variance in the RILs,BC1F1s,and HMPs,respectively (Table S4).qGYP10.3 showed largest phenotypic variations of 10.0% and 21.0% detected in BC1F1s in HF and NN,and 7.3%and 19.4%in HMPs in HF and NN,respectively.These 39 M-QTL included 17 additive,19 OD-type,and 3 UD-type(Fig.2).Six M-QTL were identified in both environments,including two QTL(qGYP2.3 and qGYP4.2)in BC1F1s,two(qGYP1.3 and qGYP4.1)in HMPs,and two(qGYP3.3 and qGYP10.3)in both BC1F1s and HMPs.

Among 49 M-QTL with additive effects detected in the RILs,92 with additive and dominant effects detected in BC1F1s,and 69 dominant effects seen in HMPs,there was only one M-QTL(qPGWC1.2)commonly found between the RILs and HMPs in NN,and eight(qBBR6.3,qMRR6.3,qPGWC1.2,qPGWC3.3,qPGWC5.3,qDEC1.2,qGYP1.1,and qGYP2.3)commonly detected between the RILs and BC1F1s in NN,suggesting that most additive loci detected in the RILs lacked heterotic effects because they were not detected in HMPs.And the fixed favorable alleles from Q9311 or WSSM at above eight M-QTL did contribute to BC1F1’s performance but without significant heterotic effect.A total of 30 M-QTL was generally detected between BC1F1s and HMPs,indicating that these parts of loci not only contributed to BC1F1performances,but also had contributions to heterosis.

3.3.QTL epistasis for grain quality and yield traits

A total of 198,119,and 90 E-QTL pairs were identified in RILs,BC1F1s,and HMPs,respectively(Tables S5,S6).Of these,11,8,and

2 pairs involved one M-QTL and the other a random locus in RILs,BC1F1s,and HMPs,respectively.No E-QTL was observed between two M-QTL.The total phenotypic variation explained(PVE)of EQTL was 33.3(11.4)%,17.5(27.3)%,and 4.6(20.6)% for BRR,27.3(19.2)%,16.5(22.2)%,and 18.4(14.3)% for MRR,39.2(18.0)%,33.1(29.6)%,and 21.7(17.4)% for HMRR,45.8(46.4)%,60.9(44.8)%,and 24.6(15.3)% for PGWC,47.0(38.9)%,35.2(47.1)%,and 21.8(15.3)% for DEC,and 14.9(15.3)%,31.4(37.4)%,and 27.6(0)% for GYP in HF(NN)in RILs,BC1F1s,and HMPs,respectively,indicating that epistasis played a very important role in heterosis.Around 37.8%,41.9%,and 34.4% of common epistasis were detected for BRR,MRR,and HMRR between BC1F1s and HMPs in the two environments,significantly higher than those of 4.0%,3.3%,and 7.7%of PGWC,DEC,and GYP,respectively,indicating that epitasis contributed more to the three milling quality traits with heterotic effect than to the appearance quality and yield traits in BC1F1s.

Table 1Phenotypic performance and broad-sense heritability of grain milling and appearance quality as well as grain yield in the parent lines,QYSM(Q9311A×WSSM),RILs,BC1F1s(Q9311A×RILs),and their HMPs in Hefei(HF)and Nanning(NN).

Table 2Correlations between quality and yield traits estimated in RIL population and BC1F1(Q9311A×RILs)population in Hefei(upper diagonal)and Nanning(lower diagonal).

Fig.1.Genomic distributions of 142 main-effect QTL for quality and yield traits identified in the RILs,BC1F1s,and HMPs in Hefei(HF)and Nanning(NN).QTL represented in the figure are those stably identified in both environments.BRR,brown rice rate;MRR,milled rice rate;HMRR,head milled rice rate;PGWC,proportion of grain with chalkiness;DEC,degree of endosperm chalkiness;GYP,grain yield per plant.

Fig.2.Numbers of M-QTL classified by gene action for quality and yield traits detected in RILs,BC1F1s,and HMPs.BRR,brown rice rate;MRR,milled rice rate;HMRR,head milled rice rate;PGWC,proportion of grain with chalkiness;DEC,degree of endosperm chalkiness;GYP,grain yield per plant.

Phenotypic performances of interaction genes in quality and yield traits were further dissected using previously cloned genes.Among the known cloned genes within the detected M-QTL,we found a pair of genes Wx(located together with qBRR6.1,qMRR6.1,qPGWC6.1,and qDEC6.1 in the region of 486,085-825,521 bp on chromosome 6)×Du1(located together with qPGWC10.4 and qDEC10.1 in the region of 19,112,236-19,801,017 bp on chromosome 10)interacting with each other.Du1 functions in starch biosynthesis by influencing the splicing of Wxband the expression of other genes involved in the rice starch biosynthesis pathway[42].If genotypes of Q9311B and WSSM are abbreviated as A and B,respectively,there were four different homozygous genotypes(A×A,B×B,A×B,and B×A)at the two interacting genes can be classified in RILs,and three heterozygous genotypes(A×H,H×A,and H×H)classified in BC1F1s.Comparing the trait performance between homozygous genotype and heterozygous genotype,we could identify favorable combinations balancing grain quality and yield traits.When the genotype combination of Wx×Du1 was A×H,the quality traits were highest,but yield was lowest.While both genes were heterozygous(H×H),yield and quality were in a well-balanced state(Fig.3A).We identified two other genes,one being YDA1(located with qPGWC4.1 in the region of 19,962,993-30,642,129 bp on chromosome 4)×Gn1a(located with qBRR1.1,qHMRR1.1,and qGYP1.1 in the region of 3,958,628-5,799,607 bp on chromosome 1),both affecting yield,as indicated that YDA1 had been shown to positively regulate rice grain size and weight of rice through the YDA1-OsMKK4-OsMAPK6 signaling pathway[43],and Gn1a affected yield by controlling the number of spikelets[44].Another pair of genes GW5(located with qHMRR5.1,qPGWC5.1,qDEC5.1,and qGYP5 in the region of 4,578,959-5,986,015 bp on chromosome 5)×GL7(located with qHMRR7.2 in the region of 22,056,455-26,044,575 bp on chromosome 7),which affected both quality and yield.GW5 played a role in the brassinosteroid signaling pathway,regulating rice grain width and weight[9],and OsSPL16-GL7 regulatory module determined the grain shape and improved the yield and quality of rice at the same time[45].The genotype combinations of above two pairs of interactions(YDA1×Gn1a and GW5×GL7)with A×A reached the tradeoff of high-yield and good-quality(Fig.3B,C).

3.4.QTL regions affecting GYP and grain quality traits

Seven QTL regions were identified on chromosomes 1,2,4,5,8,and 10 that simultaneously affected GYP and milling and appearance quality traits in the same environment(Fig.4;Table S7).Three regions were detected with opposite directions of favorable allele effects for GYP and grain quality traits,including the region of 6,172,143-8,331,441 bp on chromosome 1,harboring qHMRR1.2 and qGYP1.3 detected in HMPs in HF,qGYP1.3,qPGWC1.2,and qDEC1.2 detected in BC1F1s in HF,and qHMRR1.2,qGYP1.3,and qPGWC1.2 detected in HMPs in NN;the region of 26,663,354-27,3 16,157 bp on chromosomes 2 harboring qDEC2.2,qPGWC2.3,qHMRR2.1,and qGYP2.3 detected in BC1F1s in NN;and the region of 24,367,369-25,270,058 bp on chromosome 8 harboring qDEC8.3 and qGYP8 detected in RILs in HF.Four regions were detected with consistent directions of favorable allele effects for GYP and grain quality traits,including the region of 8,135,122-8,389,827 bp on chromosome 2 harboring qPGWC2.1 and qGYP2.2 detected in BC1-F1s in HF,the region of 5,111,957-6,692,156 bp on chromosome 4 harboring qMRR4.1 and qGYP4.2 detected in BC1F1s in NN,the region of 4,578,959-5,986,015 bp on chromosome 5 harboring qHMRR5.1,qGYP5,qPGWC5.1,and qDEC5.1 detected in BC1F1s in NN,and the region of 19,543,841-19,801,017 bp on chromosome 10 harboring qHMRR10.3,qGYP10.4,and qPGWC10.4 detected in HMPs in NN.It appeared that there is partial genetic overlap between GYP and grain quality with either consistent or opposite directions of favorable alleles for each M-QTL.

3.5.Candidate genes for milling and appearance quality traits

Among the 24 M-QTL affecting five milling and appearance quality detected in both environments(Table S8),seven heterotic M-QTL were selected for candidate gene analysis.Fourteen candidate genes for the seven M-QTL were identified based on haplotype analyses,and phenotypes of the haplotypes were characterized in xian and geng populations(Table S9).

Fig.3.Phenotypic performances of interacting genes for quality and yield traits in Hefei.(A)Wx×Du1,(B)YDA1×Gn1a,(C)GW5×GL7.A×A,B×B,A×B,and B×A indicate four homozygous genotypes in RILs.A×H,H×A,and H×H indicate three heterozygous genotypes in BC1F1s.Different letters represent significant differences at P＜0.05 based on Duncan’s multiple-range test.BRR,brown rice rate;MRR,milled rice rate;HMRR,head milled rice rate;PGWC,proportion of grain with chalkiness;DEC,degree of endosperm chalkiness;GYP,grain yield per plant.

Fig.4.M-QTL affecting grain quality and GYP traits detected in RILs,BC1F1s and the corresponding HMPs in Hefei(HF)and Nanning(NN).Number on the left side of the chromosome indicates the genetic position of the QTL in cM.The linkage map spans 2268 cM with a mean of 2.65 cM between neighboring bins.

For qPGWC1.2 and qDEC1.2,four candidate genes(Os01g0217050,Os01g0217400,Os01g0221500,and Os01g0221600)were identified(Fig.5;Table S9).For Os01g0217050 and Os01g0217400,five haplotypes were identified in the xian population.Only three haplotypes for Os01g0217050 and two haplotypes for Os01g0217400 were identified in geng population,and the favorable Hap 2 in xian population and Hap 1 in the geng population all significantly improved PGWC and DEC as well as GYP.Hap 1 in the geng population also significantly improved three milling quality traits.Two haplotypes were identified for each of Os01g0221500 and Os01g0221600,and Hap 2 both in xian and geng populations significantly improved PGWC and DEC but had no significant effect on GYP,while Hap 2 both in geng populations also significantly increased MRR and HMRR.

For qPGWC6.3,two candidate genes(Os06g0209600 and Os06g0211300)were identified.Two haplotypes were identified for Os06g0209600 in the xian or geng population,and Hap 2 in the xian population and Hap 1 in the geng population significantly improved PGWC but had no significant effect on GYP.Two haplotypes were identified for Os06g0211300,Hap 1 significantly improved PGWC and DEC but reduced GYP in the xian population.The Os06g0211300 had a significant effect on GYP but no significant effects on milling and appearance quality in the geng population(Table S9).

For qDEC6.5,three candidate genes(Os06g0267500,Os06g0268201,and Os06g0269300)were identified.Three haplotypes were identified for Os06g0267500 in both xian or geng populations,and Hap 2 significantly improved DEC,PGWC,GYP,MRR,and HMRR in the xian population,while Hap 4 significantly improved all milling and appearance quality and GYP in the geng population.For Os06g0268201,three and two haplotypes were identified in the xian and geng populations,respectively,Hap 2 had significantly improved DEC and PGWC but had no effect on GYP in xian population.The gene Os06g0268201 had significant effect on GYP but no effect on all milling and appearance quality in geng population.Four and three haplotypes were identified for Os06g0269300 in xian and geng population,respectively,and Hap 1 significantly improved DEC,PGWC,GYP,MRR,and HMRR in xian population,while Hap 5 significantly improved all milling and appearance quality and GYP in geng population(Table S9).

Fig.5.Exon-intron structure and haplotype analysis of four candidate genes in whole,xian,and geng populations for qPGWC1.2 and qDEC1.2.Different letters indicate differences in trait values among haplotypes by Duncan’s test at P＜0.05.

For qBBR7, three candidate genes (Os07g0636200,Os07g0636900,and Os07g0637200)were identified.Four and two haplotypes were identified for Os07g0636200 in the xian and geng populations,respectively,and Hap 2 significantly increased BRR but decreased GYP,PGWC,and DEC in the xian population.The gene Os07g0636200 had significant effects on HMRR and GYP but no effect on BRR in the geng population.Two and four haplotypes were identified for Os07g0636900 in the xian and geng populations,respectively,and Hap 2 had significantly improved all milling and appearance quality and GYP in the geng population.The gene Os07g0636900 had significant effect GYP but no effects on BRR in xian population.Five haplotypes were identified for Os07g0637200,and Hap 1 significantly increased BRR but decreased GYP and the two appearance quality traits in the xian population.The gene Os07g0637200 had no effect on milling and appearance quality and GYP in the geng population(Table S9).

One candidate gene Os08g0446700 was identified for qPGWC8,two haplotypes were identified in the geng population,and Hap 2 significantly improved PGWC,DEC,BRR,and MRR.One candidate gene Os12g0171600 was identified for qHMRR12.1,five and three haplotypes were identified in the xian and geng populations,respectively.Hap 1 significantly improved HMRR,MRR,GYP,PGWC,and DEC in the geng population.The gene Os12g0171600 had no effect on HMRR but had significant effects on BRR,PGWC,DEC,and GYP in xian population(Table S9).

Based on the above phenotypic performances of candidate genes for grain quality M-QTL,many candidate genes showed pleiotropic effects with different favorable haplotypes between the xian and geng populations,suggesting that these genes were associated with variation in grain quality between the two subspecies.

4.Discussion

4.1.Genetic mechanisms of heterosis for milling and appearance quality

Despite the extensive study of the genetic mechanisms of heterosis for rice yield and its components[46-49],few studies have been conducted on heterosis of grain quality traits.In the present study,single-locus overdominance was the primary genetic role in rice yield(GYP)heterosis,in agreement with Li et al.[50]and Melchinger et al.[51].In contrast,single-locus overdominance and epistasis contributed mainly to the heterosis of rice milling and appearance quality.Three OD-type M-QTL(qBRR6.3,qMRR6.3,and qBRR6.4)and a UD-type M-QTL(qMRR6.1)were identified as heterotic loci for milling quality and mapped together with three known genes associated with grain shape and a cloned gene related to starch synthesis,i.e.,GW6a[11],GL6[52],GS6[53],and Wx[54].GW6a and GL6 were both positively involved in regulating grain size and grain weight[11,52],and GS6 negatively regulated grain size in rice[55].Generally,milled rice rate is positively correlated with grain weight[56,57].It is well explained that GW6a and GL6 were identified with overdominant effects on heterosis of BRR and MRR,while GS6 was identified with an underdominant effect on MRR heterosis.Rice Wx gene,encoding granule-bound starch synthase,showed incomplete dominance to waxy allele(wx)and gene dose effect,which is the primary gene controlling amylose synthesis and directly affected amylose content in rice endosperm and pollen[58].Moreover,there was a negative correlation between amylose content and MRR[49].Thus,this study showed the Wx gene to have an underdominant effect on MRR heterosis.

Similarly,two UD-type M-QTL(qPGWC10.4 and qDEC10.1),an OD-type M-QTL(qPGWC4.1),and an ID-type M-QTL(qPGWC3.3)were identified for heterosis of appearance quality and mapped together with known genes such as starch synthase gene Du1[59],grain-filling gene GIF1[13],and three cloned genes related to grain size and grain weight:YDA1[43],qTGW3[60],and SRL2[61].Du1 may function as a regulator of starch biosynthesis by affecting the splicing of Wxbpre-mRNAs and the expression of other genes involved in rice starch biosynthetic pathways[42].GIF1 and YDA1 are both positively involved in regulating grain size and grain weight[13,43],while qTGW3 negatively regulates rice grain size and grain weight[60].Grain weight is positively correlated with rice appearance quality[62,63],while amylose content has a significant positive correlation with chalky grain rate[64].Thus,the biological functions of these genes affecting grain weight and starch synthesis may determine the genetic effects on heterosis of appearance quality.

The number of instances of epistasis with heterotic effect for milling quality was much greater than that for appearance quality.Comparing epistasis observed in BC1F1s and HMPs in the two environments,more epistasis was found between BC1F1s and HMPs in the three milling quality than in the two appearance quality traits.Specifically,respectively 48.3%,56.5%,and 47.8% of epistasis was estimated for BRR,MRR,and HMRR in BC1F1s and HMPs in both environments,much higher proportions than the 5.3% for PGWC and 6.7% for DEC.In view of the number of M-QTL detected for milling and appearance quality,we conclude that heterosis of BRR,MRR,and HMRR was controlled mainly by epistasis while that of PGWC,DEC,and GYP was controlled mainly by OD or UD,suggesting that the genetic mechanism underlying milling quality differed from that underlying appearance quality.

4.2.Influence of environment on heterosis for milling and appearance quality

Crop yield heterosis is affected by environmental conditions[65-67].We previously reported[68]that the role of heterotic loci for yield was affected by genetic background and environment.However,it is unclear how environmental factors affect the heterosis of rice quality traits.The present study revealed a similar proportion(80%)of environment-specific M-QTL for milling heterosis compared to 82% of environment-specific M-QTL for yield heterosis.However,a slightly lower proportion(76%)of environment-specific M-QTL was seen for appearance heterosis.The environment-specific M-QTL for milling(appearance)quality heterosis with PVE＞3% were identified mainly in NN 56.5(66.7)%.The loci qBRR7,qMRR6.3,qHMRR12.1,qPGWC1.2,qPGWC5.3,qPGWC6.3,qDEC1.2,and qDEC10.1 were identified with additive effects in HF but as the heterotic loci in NN,and the opposite was true for qDEC6.5.In addition,numerous E-QTL were detected for heterosis of quality traits that were environmentspecific.These results suggest that the heterotic loci of milling and appearance quality in hybrid rice were differently affected by environmental conditions.

4.3.Comparisons of stably detected grain quality M-QTL with previously cloned genes and their candidate gene inference

M-QTL that affected the same grain quality simultaneously identified in two environments were compared with genes previously reported to lie at the same or nearby physical positions.Some identified M-QTL were located near previously reported genes.For instance,qBRR6.1 located around 825,521 bp on chromosome 6,which affected BRR in RILs,was mapped in the same region with GS6,which negatively regulates grain size[55]and Wx,which affects amylose content[58].qMRR6.3 in the 26,159,547-26,303,837 bp region on chromosome 6,which affected MRR in the RILs,was mapped together with GW6a,which is positively involved in regulating grain size and weight[11].qHMRR1.1 in the 3,958,628-4,980,670 bp region on chromosome 1,which affected HMRR in BC1F1s,was mapped together with D2 for grain size,grain number and yield[69],Gn1a for grains per panicle[44],and GW5L negatively regulates grain width and weight[70].qHMRR7.2 in the 22,056,455-26,044,575 bp region on chromosome 7,which affected HMRR in RILs,was mapped together with OsGASR9,which positively regulates grain size and yield[71],and GL7 for grain size[12].qPGWC2.4 in 33,034,282-35,705,713 bp region on chromosome 2,which affected PGWC in RILs,was mapped together with OsbHLH107 for grain size[72].qPGWC5.1 and qDEC5.1 in 3,981,395-5,986,015 bp region on chromosome 5,which affected PGWC and DEC in BC1F1s,respectively,were mapped together with GW5 for grain width and weight[9].qPGWC10.4 and qDEC10.1 in 19,307,415-19,801,017 bp region on chromosome 10,which affected PGWC in HMPs in the two environments,and DEC in RILs in HF and BC1F1s and HMPs in NN,respectively,were mapped together with Du1,which regulates starch biosynthesis[59].GW6a and Du1 were mapped together with qMRR6.3 and qPGWC10.4 with OD and UD heterotic effects,respectively.These two genes probably had heterotic effects on grain milling and appearance quality.Allelisms between M-QTL for grain quality mentioned above and the previously reported genes await confirmation via further fine-mapping and gene cloning.

For the seven stable M-QTL newly identified with heterotic effects(Table S9),candidate genes were further inferred based on other supporting evidence such as gene expression database.For instance,the region 6.39-6.77 Mb on chromosome 1 harboring qPGWC1.2 and qDEC1.2 contain four candidate genes,only Os01g0221500 was expressed in seeds(https://rice.uga.edu/cgibin/ORF_infopage.cgi?orf=LOC_Os01g12200),and encodes an uncharacterized conserved protein,so it is the most likely candidate gene.Of the two candidate genes for qPGWC6.3,Os06g0211300 is seemingly a candidate gene that encodes a sterol carrier protein,and participates in lipid metabolic process.The gene has higher expression level in seeds(http://rice.uga.edu/cgibin/ORF_infopage.cgi?orf=LOC_Os06g10890),suggesting that it may function in seed development,thereby affecting milling and appearance quality.Of the three candidate genes of qDEC6.5,Os06g0269300,encoding an NHL(abbreviation of the proteins NCL1,HT2A,and LIN-41)repeat-containing protein in rice,is probably a candidate gene.A gene(FUWA)encoding an NHL domain-containing protein has been reported[73].Compared with wild-type plants,fuwa plants are slightly shorter,and the kernels are wider,thicker and shorter with increase in endosperm weight,filling rate,and thousand-kernel weight.Of the three candidate genes for qBRR7,Os07g0636900 encodes a pentatricopeptide repeat(PPR)domain protein in rice.Previous studies[74]have demonstrated that PPR proteins are essential for maize kernel formation.Loss of function of specific PPR proteins results in empty pericarps and tiny,malformed kernels in various genetic backgrounds.Thus,Os07g0636900 is probably a candidate gene for qBRR7.Os12g0171600,a candidate gene for qHMRR12.1,encodes an expressed protein of unknown function with a higher expression level in seeds(http://rice.uga.edu/cgi-bin/ORF_infopage.cgi?orf=LOC_Os12g07380).Thus,it is most likely a candidate gene for qHMRR12.1.

4.4.Breeding strategy aimed at high yield and quality in hybrid rice

In this study,we identified some M-QTL for different traits in the same genomic intervals,such as a region on chromosome 5(4.58-5.99 Mb)containing qHMRR5.1,qPGWC5.1,qDEC5.1,and qGYP5(Table S7).This region contains the well-known grain size gene GW5[9].The WSSM alleles at the identified QTL increased GYP and HMRR and reduced PGWC and DEC,suggesting that GW5 may be used for improving both yield and grain quality of hybrid rice simultaneously.Similarly,it would be possible to increase grain yield and grain quality together by pyramiding favorable alleles with consistent directions of gene effects for GYP and grain quality at the QTL regions of 8,135,122-8,389,827 bp on chromosome 2,5,111,957-6,692,156 bp on chromosome 4,and 19,543,841-19,801,017 bp on chromosome 10,to achieve ideal combinations of M-QTL associated with yield and quality.Some pairwise genetic interactions with their favorable genotypic combinations,such as Wx×Du1 with A×H,YDA1×Gn1a,and GW5×GL7 with A×A could be used to improve high yield and quality in hybrid rice breeding through marker-assisted pyramiding of favorable alleles at each locus that are dispersed among multiple parents.

In this study,eight M-QTL(qBRR6.3,qMRR6.3,qPGWC1.2,qPGWC3.3,qPGWC5.3,qDEC1.2,qGYP1.1,and qGYP2.3)were identified to contribute to BC1F1’s performance but with negligible heterotic effect,given that they were detected only in BC1F1s and not in HMPs.The same fixed favorable alleles with additive effects in Q9311B and WSSM at these eight M-QTL influenced the performance of BC1F1s.Thus,introgression of favorable alleles at selected loci into two parents of a hybrid can be performed,and pyramiding of these favorable alleles in the parents will fix the additive effects and improve the performance of the hybrid F1for both GYP and grain quality.This is why elite inbred lines are the basis for successful hybrid breeding.The finding that 30 M-QTL were detected in both BC1F1s and HMPs suggests that these loci contributed to BC1F1performances,and contributed to heterosis:that is,the 30 M-QTL contributed to BC1F1performance with their heterotic effects.Thus,further improvement of yield and quality in hybrid rice may be realized by pyramiding OD/D loci that are polymorphic between the parents and minimizing or eliminating UD loci that are polymorphic between the parents.

There are marked variations in grain appearance and grain shape between the two subspecies.For instance,the geng-type alleles of qSW5/GW5,OsSPL16/GW8,and GS6 increase the width of rice grain[9,55,75-78],while the xian-type ones of GS5 and OsSPL13/GW7 increase grain length[10,79].Rice milling and appearance quality is highly correlated with grain shape[80,81].In the present study,many candidate genes showed differing phenotypic effects on grain quality and GYP between the two subspecies.So,in rice breeding for grain quality and GYP,we can transfer favorable geng alleles (haplotypes) of Os01g0217050,Os01g0217400,and Os07g0636900 into a xian background to improve milling quality and GYP of a xian cultivar,or transfer favorable xian alleles of Os06g0211300,Os06g0268201,Os07g0636200,and Os07g0637200 into a geng background to improve appearance quality and GYP of a geng cultivar,or pyramid these non-allelic favorable alleles to simultaneously improve milling and appearance quality as well as GYP in rice breeding to exploit xian-geng inter-subspecific heterosis.

5.Conclusions

Single-locus overdominance and epistasis were the main contributors to heterosis for rice milling and appearance quality,in distinction from yield heterosis,which was contributed primarily by single-locus overdominance.Epistasis contributed more to heterosis of milling quality than to that of appearance quality,and some candidate genes of grain quality M-QTL showed phenotypic effects on milling and appearance quality traits that differed between the two subspecies.Three strategies are proposed for increasing heterosis of grain yield and grain quality:to introgress favorable allele at QTL or genes underlying grain quality and yield such as GW5 or to pyramid favorable alleles with consistent directions of gene effects for GYP and grain quality at the QTL regions on different chromosomes;to introgress numerous favorable alleles at target loci for grain yield and quality into two parents of a hybrid,and pyramid these favorable additive alleles in both parents to fix additive effects and improve the performance of the F1hybrid;and to pyramid OD/D loci that are polymorphic between the parents and minimize or eliminate UD loci that are polymorphic between the parents.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement

Hui You:Investigation,Writing-original draft.Sundus Zafar:Writing-review&editing.Fan Zhang:Data curation,Formal analysis.Shuangbing Zhu:Investigation.Kai Chen:Investigation.Congcong Shen:Investigation.Xiuqin Zhao:Investigation.Wenzhong Zhang:Conceptualization.Jianlong Xu:Funding acquisition,Project administration,Conceptualization,Writing-review & editing.

Acknowledgments

This work was funded by the Key Research and Development Project of Hainan Province(ZDYF2021XDNY128),the Hainan Yazhou Bay Seed Lab Project(B21HJ0216),and the Agricultural Science and Technology Innovation Program and the Cooperation and Innovation Mission(CAAS-ZDXT202001).The authors thank Win-All Hi-Tech Seed Co.,Ltd.for their assistance with planting.

Appendix A.Supplementary data

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