Bo Zhng,Feixing Qi,Gng Hu,Yiki Yng,Li Zhng,Jinghu Meng,Zhongmin Hn,Xingchun Zhou,Hiyng Liu,b,Mohmmed Ayd,c,Yongzhong Xing,*

a National Key Laboratory of Crop Genetic Improvement and National Center of Plant Gene Research,Huazhong Agricultural University,Wuhan 430070,Hubei,China

b Hubei Collaborative Innovation Center for Grain Industry,Yangtze University,Jingzhou 434100,Hubei,China

c Plant Research Department,Nuclear Research Center,Atomic Energy Authority,Abo-Zaabal 13759,Egypt

Keywords:

A B S T R A C T Bulked-segregant analysis is a time-and cost-saving strategy for identifying major quantitative trait loci(QTL)in a mapping population.Bulked-segregant analysis combined with whole-genome sequencing(BSA-seq)was performed to rapidly identify QTL for heading date,plant height,and panicle length in a large F2 population derived from two landraces:Chuan 7(C7)and Haoboka(HBK).Twenty plants with extremely low or high phenotypic values for the target traits were selected from an F2 population of 940 plants to construct low-and high-value bulks.Three pairs of bulks for the three traits were constructed,resulting in six DNA pools.BSA-seq revealed nine QTL:four for heading date,three for plant height,and two for panicle length.These QTL were validated in a random F2 population or BC4F2 populations.The major novel plant height QTL,qPH8,acting additively with an effect equivalent to that of semi-dwarf 1(sd1),is potentially valuable for hybrid rice breeding.qPH8 controls mainly the elongation of basal internodes.The C7 allele of qPH8 reduces plant height and increases lodging resistance without yield penalty,suggesting a potential role for qPH8 in improving rice plant architecture.

1.Introduction

Rice is an essential component of the diets and livelihoods of more than 3.5 billion people[1].During the past 60 years,rice grain yields have increased largely owing to dwarfing breeding and heterosis utilization and have further increased since the 1990s owing to the development of‘‘green super rice”and molecular breeding combined with advanced cultivation techniques[2–4].However,to meet the demands of increasing population,improving rice production remains urgent.

Mining new QTL for yield traits is the key to further increasing rice production.QTL mapping is an efficient way to identify new genes.Conventional genetic mapping usually uses all individuals collected from a population to perform genotype analysis using an array of molecular markers[5].This technique requires genotyping and phenotyping a large number of progeny from a biparental mapping population to ensure enough statistical power,a time-consuming and labor-intensive task.Bulked-segregant analysis(BSA)has been used to overcome this issue by genotyping only lines with extreme phenotypes instead of a large number of plants in the mapping population[6,7].The initial application of BSA in rice involved searching for markers linked to QTL associated with fertility,disease resistance,or photoperiod sensitivity and expected to be controlled by major genes[8–10].With the swift development of sequencing technology,genome-wide association studies(GWAS)have accelerated the identification of candidate genes for a given trait owing to the genome-wide availability of super-high-density single nucleotide polymorphism(SNP)markers[11–13].Nevertheless,whole-genome sequencing of a large number of plants or lines in a mapping population is still expensive.

BSA combined with whole-genome sequencing(BSA-seq)can dramatically reduce genotyping costs by using selective sampling,and the statistical power in QTL-mapping is comparable to that of full-population analysis[5,14].In BSA-seq,plants of extreme phenotypes for a target trait from a segregating population are pooled separately and sequenced using a next generation sequencing(NGS)platform.The allele frequencies in the pools are then calculated and compared to identify genomic regions with frequencies differing significantly between the bulks and thus associated with the target trait.The confidence interval for the associated genomic region is determined by a simulation approach[14,15].The BSAseq approach has been successfully used to map QTL for agronomic traits in multiple plant populations,such as F2populations,recombinant inbred lines(RILs)and backcross populations[16–18].The availability of analytical tools in the form of standalone software packages such as QTL-seq[14]and QTLseqr[19]has greatly simplified BSA-seq.

Heading date,plant height,and panicle length are agronomic traits of rice that contribute to yield by determining regional adaptability,plant architecture,and panicle size,respectively.A series of genes associated with rice photoperiodic flowering have been extensively studied and can be roughly summarized in two regulatory pathways:theOsGI-Hd1-Hd3a/RFT1conserved pathway(comparable to that inArabidopsis)and theEhd1-Hd3a/RFT1ricespecific pathway[20,21].Among these genes,several that encode flowering repressors with natural variation,such asGhd7[22],Ghd8/DTH8[23,24],Ghd7.1/OsPRR37[25,26]andHd16[27],delay rice heading by downregulatingEhd1expression under long-day conditions,and combinations of these genes largely define the ecogeographical adaptations and yield potential of cultivated rice[28,29].Some of them,such asGhd7[30],Ghd8[31],Ehd1[32],andHd1[33],have also been applied to the directional genetic improvement of rice heading date.Rice plant height is controlled by numerous genes that act in a complex regulatory network and that are involved mainly in the biosynthesis or signal transduction of phytohormones such as gibberellins(GA),brassinosteroids,and strigolactones[34].The single recessive genesd1,which is known as a‘‘green revolution”gene,encodes a key enzyme(GA20ox2)involved in the penultimate step of the GA biosynthetic pathway,and its defective alleles have been used extensively in modern rice breeding[35,36].Hundreds of QTL for rice panicle length have been identified,and most of them simultaneously control other traits,such as heading date,plant height,and tiller number(https://archive.gramene.org/qtl).The functions of genes controlling panicle length are diverse,including peptide transporters such asSP1[37]andOsNPF7.7[38],transcriptional activators such asLAX1[39]andSP3[40],and growth-regulatory factors such asOsGRF4[41].

In this study,BSA-seq was performed in a large F2population derived from two landraces to identify further QTL for heading date,plant height,and panicle length.Single-marker analysis was performed for QTL validation using F2and BC4F2populations,and polymorphism analysis and fine-mapping were performed to identify new QTL for these three traits.

2.Methods and materials

2.1.Plant materials and field experiments

An F2population was generated from a cross between maternal Haoboka(HBK)and paternal Chuan 7(C7).Single sequence repeats(SSR)and insertion/deletion(InDel)markers with polymorphism between C7 and HBK and evenly distributed on the 12 chromosomes were used to construct chromosome-segment substitution lines(CSSLs)by marker-assisted selection(MAS)to ensure that the genome fragments of the donor parent in the CSSLs covered the entire genome.A total of 100 BC4F1plants with HBK as the recurrent and C7 as the donor parent were selected by MAS and self-pollinated to produce 100 BC4F2populations.Seeds were sown in a seedling bed at the experimental stations in Wuhan(Hubei,China)or Lingshui(Hainan,China),and 25-day-old seedlings were transplanted into the field.A total of 940 F2plants(F2_940)were grown at Wuhan in the summer of 2016.The first 200 plants of the F2_940 population were selected as a random subpopulation(F2_200)for subsequent QTL validation.The BC4F2populations for QTL validation were grown at Wuhan in the summer of 2017.A BC4F3population consisting of 2000 plants was grown at Wuhan in the summer of 2018 for fine mapping ofqPH8,and the progeny of the recombinants were planted in Lingshui in the winter of 2018 for progeny tests.

2.2.Phenotype measurement and bulk construction

F2plants were measured in the field for three traits:heading date,plant height,and panicle length.Heading date was scored as the number of days from sowing to the emergence of the first panicle on the plant.Plant height and panicle length were measured 20 days after rice heading.Plant height was recorded as the distance from the soil surface to the tip of the tallest panicle,and panicle length was scored as the mean length of the three tallest panicles.Among the NILs forqPH8,thousand-grain weight,grain yield per plant and total number of spikelets per plant were measured with the Yield Traits Scorer[42].Spikelets per panicle were recorded as the total number of spikelets divided by the number of panicles.Lodging index was measured following the previously reported method[43]with minor modifications.Twelve main culms per line were sampled 20 days after heading,and then the center site of the fourth internode(breaking point)of each culm was selected for measurement of bending strength(BS)using a plant stem strength-testing machine YF-1200(Henan Yunfei Technology Development Co.,Ltd.,Zhengzhou,Henan,China).The distance between the two fulcrums was constant at 6 cm.The lodging index of the fourth internode was calculated by the following formula:

The length is the distance from the breaking point to the top of the panicle.Weight is the total fresh weight of the main culm above the breaking point,including the panicle.

Twenty early-heading and 20 late-heading plants were collected from F2_940 to construct the low-value(HD-E)and highvalue(HD-L)bulked pools for heading date.Twenty of the shortest and 20 of the tallest plants formed the low-value(PH-S)and highvalue(PH-T)bulked pools for plant height.Likewise,two sharply contrasting bulks(PL-S and PL-L for shortest and longest panicles)for panicle length were constructed.

2.3.DNA extraction and genotyping

High-quality genomic DNA of 200 random plants(F2_200),six bulks,100 BC4F1plants,targeted BC4F2populations,BC4F3population and progenies of recombinants were extracted with a modified cetyl-trimethyl ammonium bromide(CTAB)method[44].BC4F1plants were genotyped with the RICE6K SNP array(China National Seed Group Co.,Ltd.,Beijing,China).Genotyping by sequencing was used for each bulk,with 25-fold genome coverage.Genomic DNA of each individual in the F2_200,BC4F2and BC4F3populations was amplified with primers using rTaq polymerase from TaKaRa in buffer I according to the manufacturer’s indications.For each PCR,DNA was initially incubated for 4 min at 95 °C,followed by 35 cycles of amplification(95 °C for 30 s,58°C for 30 s and 72 °C for 30 s).SSR markers were selected from the Gramene Markers Database(https://archive.gramene.org/markers)based on their physical position.InDel markers were selected from a LInDel[45]marker set designed for amplifying large fragment insertions/deletions between twoindicarice lines,Zhensan 97 and Minghui 63,and onejaponicarice line,Nipponbare.All markers used in this study are listed in Table S1.

2.4.Whole-genome sequencing,read mapping and SNP calling

Whole-genome sequencing was performed on the Illumina HiSeq Xten platform(Tianjin,China).Raw reads were processed with Trimmomatic[46]to remove low-quality reads and any adapter contamination.Clean reads containing 150 bp paired-end reads were aligned to the Nipponbare reference genome(release 7)using Burrows-Wheeler Aligner(BWA)software[47,48].SAMtools[49]was used to convert sam files to bam files,and Picard[49]software was employed to fix matches,sort read groups,and remove PCR duplicates from bam files.SNPs were called with the Genome Analysis Toolkit(GATK)[50],generating a VCF file for both paired bulks for each trait.SNPs between parents were initially used to calculate SNP indices for paired bulks.All the software and commands were run in a Linux server environment.

2.5.BSA-seq analysis

NGS-based BSA was performed with QTLseqr[19],an R package for bulked-segregant analysis combined with next-generation sequencing,derived from the original analysis pipeline of QTLseq[14].TheVariantsToTabletool[51]was used to extract the necessary fields from the VCF file containing SNP variants.SNP data were imported into QTLseqr with theimportFromGATK()function,and the total reference allele frequency,per-bulk SNP-index,and ΔSNP-index of each SNP were calculated.ThefilterSNPs()function offered options for further filtering SNPs to reduce noise and improve results based on reference allele frequency,total read depth,per-bulk read depth,and genotype quality score.After SNP filtering,therunQTLseqAnalysis()function was performed to count the SNPs within the set window bandwidth(1 Mb)and calculate a tricube-smoothedΔSNP-index.TheΔSNP-index was simulated for more than 10,000 replicates for each bulk based on the population type and size,and the quantiles(95%and 99%)from the simulations were used to estimate confidence intervals.An actualΔSNP-index averaged over a sliding window that lies outside the confidence interval indicates a putative QTL.TheplotQTLStats()andgetQTLTable()functions were called to extract the analysis results and significant QTL regions.The tricubesmoothedΔSNP-index values for heading date,plant height,and panicle length were calculated in genomic windows of 1 Mb and plotted against all 12 rice chromosomes(Figs.S1–S3).QTL were identified at the 95% confidence level based on theΔSNP-index peaks.All commands and codes are available at https://github.com/bmansfeld/QTLseqr.

2.6.Validation of QTL in F2_200 and BC4F2 populations

The QTL identified by BSA-seq were validated in the F2_200 and BC4F2populations by single-marker analysis of variance(ANOVA).QTL flanking markers were selected to genotype the plants of the mapping populations.Plants in a population were classified into three groups according to their marker genotype,and then ANOVA was used to test for significant differences among group means.AP-value＜0.05 was used to define a significant marker–QTL association.TheR-squared value from ANOVA was interpreted as the proportion of phenotypic variation explained(PVE)by the QTL.

2.7.Fine mapping of qPH8

The flanking markers ofOsSPY,a negative regulator of gibberellin signaling located inqPH8region[52],were used to identify recombinants between the markers,andOsSPYwas identified in the large BC4F3population comprising 2000 plants to excludeOsSPYas the candidate gene ofqPH8.The phenotypes of plant height in selected recombinants were confirmed by self-progeny tests.The plant height of each homozygous genotype in the progeny of recombinants was also used to identify significant differences among group means.Details of the markers used for fine-scale mapping are presented in Table S1.

3.Results

3.1.Phenotypic variation of three traits in the F2 population and bulked pools

Compared with HBK,C7 showed lower values of heading date and plant height but greater panicle length(Table S2).All three traits showed continuous distributions with transgressive segregation in both directions in the F2_940 population(Fig.1A–C).Panicle length was highly correlated with plant height but not with heading date,and no significant correlation was detected between plant height and heading date(Table S3).The variation in heading date,plant height,and panicle length in the population was large,ranging from 80 to 143 days,90 to 218 cm and 18 to 42 cm,respectively.The F2_200 population showed similar frequency distributions for these traits(Fig.1A–C).Twenty extremephenotype plants from each tail of the trait distribution in the F2_940 population were selected to form a bulked pool(Fig.1A–C).The bulks of heading date were mutually exclusive with bulks of plant height and panicle length,but bulks of plant height and panicle length were not mutually exclusive.Five percent of the plants were common to PH-S and PL-S,and 10% of plants were common to PH-T and PL-L.The mean heading dates of HD-E and HD-L were 84.2 and 140.6 days,respectively.Similarly,the mean plant heights of PH-S and PH-T were 107.2 and 202.1 cm and the mean panicle lengths of PL-S and PL-L were 22.0 and 36.7 cm(Fig.1D–F).All the bulks showed higher or lower phenotypic values than both parents,indicating that QTL with positive effects were present in both parents(Fig.1D–F).The mean phenotypic values of other two traits in each pair of bulks were also recorded.There was no significant difference in either plant height or panicle length between HD-E and HD-L(Table S4).There was no significant difference in heading date between PH-S and PH-T,but the panicle length of PH-S was much shorter than that of PH-T(Table S4).Similarly,there was a significant difference in plant height but no significant difference in heading date between PL-S and PL-L(Table S4).

3.2.Whole-genome sequencing of bulks

A total of 577.4 million raw paired-end reads were generated for the six bulks and two parents by whole-genome sequencing(Table S5).After removal of low-quality reads and adapter contaminants,a total of 575.8 million clean reads and 86.4 Gb of clean bases were retained.The mean proportion of reads with Q30 scores was 95.7%,and the mean GC content was 44.3%(Table S5).Of these sequencing data,76.8 and 67.5 million clean reads were generated from HD-L and HD-E,respectively,accounting for respectively 30.1×and 26.2×coverage of the rice genome.Respectively 98.9%and 98.5% of the reads from HD-L and HD-E were mapped to the reference genome of Nipponbare(Table S5).Respectively 76.5 and 51.4 million clean reads were generated from PH-T and PH-S,accounting for 29.5×and 19.9×genome coverage.The mapping proportions of reads from PH-T and PH-S were 97.7% and 97.8%,respectively(Table S5).Approximately 70.1 and 72.2 million clean reads were obtained from PL-L and PL-S,respectively,accounting for 27×and 28×genome coverage.Respectively 97.4%and 97.8% of reads from PL-L and PL-S were mapped to the reference genome(Table S5).Variant calling resulted in 511,393 SNPs for heading-date bulks,543,319 SNPs for plant-height bulks,and 455,074 SNPs for panicle-length bulks.After further SNP filtering with QTLseqr,98,301,79,877 and 70,444 SNPs for heading date,plant height,and panicle length bulks,respectively,were retained for QTL mapping(Table S5).

Fig.1.Frequency distributions of three traits in F2_940 and F2_200 populations and bulk selections.(A–C)Frequency distributions of heading date(A),plant height(B)and panicle length(C)in the F2_940(gray)and F2_200(black)populations.Black and white arrows indicate the mean phenotypic values of C7 and HBK,respectively.(D–F)Comparisons of heading date(D),plant height(E)and panicle length(F)among the paired bulks(high bulk and low bulk)and parents(C7 and HBK).Paired bulks of the target traits were selected from each tail based on the frequency distributions and are shown by black boxes.

Fig.2.QTL for heading date,plant height,and panicle length identified by QTLseqr.(A–C)Distribution ofΔSNP-index for heading date(A),plant height,(B)and panicle length(C)calculated in a 1 Mb sliding window using a tricube-smoothing kernel.Red lines and blue lines indicate 95%and 99%confidence intervals,respectively.QTL were detected at a 95% confidence level.

3.3.QTL identification by QTLseqr

Four QTL for heading date were identified(Fig.2A;Table 1).The genomic interval ofqHD3was 6,144,761–7,012,695 bp on chromosome(Chr.)3,and theΔSNP-index peak was-0.44 at 6,975,882 bp.qHD6was located in the interval of 472,611–3,563,230 bp on Chr.6,with the highestΔSNP-index of 0.57 at 2,440,541 bp.qHD9was located in the interval of 18,800,356–20,362,313 bp on Chr.9,with a peak of-0.48 at 19,391,881 bp.qHD10was located in the interval of 17,330,748–19,251,928 bp on Chr.10,with a peak of 0.46 at 18,107,096 bp.The known heading daterelated genes with natural variation were located in or near these QTL regions(Table 1).OsMADS1is located 84.4 kb fromqHD3[53],qHD6colocalizes withHd3aandRFT1[54],andEhd1is 254.7 kb fromqHD10[55].However,no known heading-date genes were identified nearqHD9(Table 1).

Table 1QTL identified for heading date,plant height,and panicle length.

Three QTL for plant height were detected(Fig.2B;Table 1).qPH1.1andqPH1.2were located in the intervals of 27,377,900–30,224,989 bp and 35,310,855–39,301,634 bp,respectively,on Chr.1.TheΔSNP-index peak ofqPH1.1was 0.52 at 28,386,060 bp,and that ofqPH1.2was 0.65 at 37,859,039 bp.qPH8was located in the interval of 27,109,783–28,439,167 bp on Chr.8,with aΔSNPindex peak of-0.56 at 27,991,461 bp.The confidence intervals ofqPH1.2andqPH8contained the known plant height genesSD1andOsSPY,respectively[56,57],whereas no plant height genes with natural variation were reported in the interval ofqPH1.1(Table 1).

Two QTL were identified for panicle length on chromosomes 1 and 5(Fig.2C;Table 1).qPL1was located in the interval of 37,565,115–41,595,779 bp,with aΔSNP-index peak of 0.46 at 38,533,219 bp.The confidence intervals ofqPL1andqPH1.2overlapped.qPL5was located in the interval of 22,545,745–25,532,28 1 bp,with aΔSNP-index peak of 0.49 at 24,331,365 bp(Table 1).No panicle length-associated genes with natural variation have been identified in this QTL region.

3.4.QTL verification in the F2 and BC4F2 populations

All QTL detected by BSA-seq were validated by single-marker analysis in either the F2_200 population or the BC4F2populations(Table 2).For heading date,all QTL exceptqHD9were detected in the F2_200 population.The additive effects ofqHD3,qHD6andqHD10were-5.8,-7.4,and 2.4 days,respectively,indicating that C7 alleles atqHD3andqHD6promoted rice heading,while the C7 allele atqHD10delayed rice heading(Table 2).qHD6is a major QTL that explained 38.6% of heading date variation in the F2progeny(Table 2).qHD3,qHD9andqHD10were also identified in the BC4F2populations,and their additive effects were,-2.2,2.2,and 7.4 days,respectively(Table 2).No BC4F2population targetingqHD6was identified for its validation.

qPH1.2andqPH8were detected in the F2_200 population.They explained 25.1%and 27.3%of the variation in plant height,respectively,with additive effects of 13.3 and-14.3 cm,respectively(Table 2).qPH1.1,qPH1.2andqPH8were also identified in the BC4F2population,and their additive effects were 4.4,16.7 and-14.6 cm,respectively(Table 2).C7 alleles atqPH1.1andqPH1.2increased plant height but those atqPH8reduced plant height.The segregation ratio of plant height atqPH1.2was 3:1(high:dwarf)in the BC4F2population,while the ratio was 1:2:1(high:medium:dwarf)atqPH8(Fig.S4).The dominance effects ofqPH1.2andqPH8were 15.3 and-1.1 cm,respectively,in the BC4F2population,indicating thatqPH1.2was dominant but thatqPH8was additive(Table 2).

qPL1andqPL5were detected in both the F2_200 and BC4F2populations(Table 2).Their additive effects were 1.8 and 1.6 cm in the F2population,respectively,and 2.7 and 1.4 cm in BC4F2population,respectively.C7 alleles at both QTL increase panicle length.qPL1andqPL5individually explained 11.3% and 12.2%,respectively,of the variation in the F2population(Table 2).

3.5.Polymorphism analysis of known genes in the QTL regions

To confirm whether the known genes were the causative genes underlying the corresponding QTL,polymorphism analysisbetween C7 and HBK was performed based on their coding sequences(Table 3).Parental-allele sequencing showed that C7 carried the same strong functionalSD1allele as did Kasalath containing SD1-GR,which significantly increased plant height[58],but HBK contained a weak allele similar to Nipponbare’s SD1-EQ.Thus,SD1was likely the gene underlyingqPH1.2/qPL1(Table 3).Sequence alignment indicated that an insertion/deletion polymorphism in the splice site ofOsMADS1resulted in an alternatively spliced protein in HBK(Table 3),but this variation did not lead to a change in heading date[59].Thus,OsMADS1was not the candidate gene ofqHD3.Hd3aandRFT1,two closely linked florigen-associated genes,were in the confidence interval ofqHD6.C7 and HBK carried type 3 and type 1 alleles atHd3a,respectively[60];atRFT1,both parents carried the alleles of previously reported Group_IV and Group_I[61],respectively(Table 3).However,there were no significant differences in function between type 3 and type 1 alleles or between Group_IV and Group_I[60,61].Further comparisons of the promoter(1.5 kb)ofHd3abetween the parents showed that C7 belonged to type 3,which was expressed at a significantly higher level than was the type 1 promoter of HBK[60](Table S6).Accordingly,Hd3awas considered the candidate gene ofqHD6.OsSPYwas in the interval ofqPH8and separately identified as Hap_III and Hap_I in C7 and HBK based on combinations of the two SNPs at+26 and+2498 from the translation start site(Table 3).However,the effects of Hap_I and Hap_III on plant height showed no significant differences in a germplasm collection[57].Ehd1was nearqHD10and identified as type 7 and type 1 in C7 and HBK(Table 3),respectively,but there were no functional differences between the two alleles[60].Sequence alignment of the promoter ofEhd1indicated that the haplotypes ofEhd1in C7 was identified as type 6 with a large deletion(49 bp),and type 1 was identified in HBK[60](Table S6).Previous study[32]indicated that the 49 bp deletion in the promoter ofEhd1resulted in down-regulation ofEhd1expression,and that introgression of this variation into Zhenshan 97 significantly delayed rice heading.Accordingly,Ehd1was assigned as the candidate gene ofqHD10.

Table 3Polymorphism analysis of coding regions of known genes associated with the corresponding QTL between C7 and HBK.

Thus,SD1,Hd3a,andEhd1were assigned as candidate genes ofqPH1.2/qPL1,qHD6,andqHD10,respectively.The other five QTL were novel.

3.6.qPH8 is a novel major QTL for plant height

The additive effect ofqPH8on plant height was comparable to that ofqPH1.2in both the F2and BC4F2populations,but the C7 alleles showed opposite effects at the two QTL(Table 2),attracting our attention.To further excludeOsSPYas the candidate ofqPH8,two flanking SSR markers,RM264 and RM477,were used to select recombinants from a BC4F3population(Fig.3A).A total of 14 recombinants were identified.Among them,six recombinants whose heterozygous regions overlapped and covered the full target region were selected forqPH8fine mapping(Fig.3B).Progeny tests confirmed that the recombinants R1,R2,R3,and R4 were heterozygous atqPH8,given large variation in plant height in their own progenies and the observation that R5 and R6 were homozygous atqPH8,with little variation in plant height.Moreover,OsSPYdid not cosegregate withqPH8in the progenies,whereas RM264 and B2(an SSR marker)did cosegregate withqPH8(Fig.3B).Accordingly,qPH8was assigned as a new QTL for plant height,located upstream of B1(Fig.3B).

The genetic effect ofqPH8was further estimated by comparison of the plant heights of two homozygous NILs.qPH8-C7andqPH8-HBKshowed large differences in plant heights,measuring 140.9 and 169.6 cm,respectively(Fig.4A–C;Table 4).Differences in plant height between the NILs became significant after the rice culms had rapidly elongated beginning at 8 weeks after sowing,suggesting that the control of plant height byqPH8was closely associated with culm elongation(Fig.4D).Although each internode ofqPH8-C7was shortened,the basal internodes,including the third to sixth internodes(counting from the panicle to the stem base),showed a higher percentage of length reduction byqPH8-C7(Fig.4E,F).These results indicated thatqPH8is a major QTL for plant height and controls internode elongation mainly of the basal rather than the uppermost internode.

A comprehensive comparison of the phenotypic performance between these two NILs was performed(Table 4).They showed no significant differences in heading date,tiller number,or seed setting rate.However,the panicles ofqPH8-C7were 3.2 cm shorter than those ofqPH8-HBK,andqPH8-C7had 9.0 more spikelets per plant thanqPH8-HBK(Table 4).The thousand-grain weight ofqPH8-C7andqPH8-HBKwere 37.8 g and 40.5 g,respectively.The larger thousand-grain weight ofqPH8-HBKbut with less spikelets per plant ultimately resulted in no difference in yield per plant in comparison withqPH8-C7(Table 4).The lodging index(0.36±0.02)of the fourth internode ofqPH8-C7was significantly smaller than that ofqPH8-HBK(0.45±0.03)(Table 4).

4.Discussion

4.1.Stringent selection for bulks in a large F2 population is important for reliable QTL detection

Fig.3.qPH8 is a novel QTL for plant height.(A)Graphical genotype of a BC4F3 individual in which a heterozygous genomic fragment containing qPH8 was introgressed into the end of the long arm of Chr.8.(B)qPH8 was mapped upstream of SSR marker B1,and the cloned gene OsSPY was not a candidate gene for qPH8.Black and white bars represent homozygous genomic fragments of C7 and HBK,respectively.Gray bars indicate heterozygotes.HH,homozygous genotype of HBK;CC,homozygous genotype of C7.Numbers in parentheses are numbers of plants.The P-value indicates the significance of plant height between HH and CC,which was calculated based on Student’s t-test.

Fig.4.qPH8 controls plant height mainly by determining basal internode elongation.(A,B)Visualization of the genome composition of two NILs,qPH8-HBK and qPH8-C7.Red lines indicate homozygous genomic regions introgressed from C7 based on SNP markers.(C)Plant height of the two NILs.(D)Plant height growth curve of the two NILs.(E,F)Comparisons of internode lengths between the two NILs.The highest internode was recorded as the first internode.Scale bars,20 cm for(C)and(E).

In this study,we selected 20 plants with extremely phenotypic values from 940 plants to construct a bulk and nine QTL for heading date,plant height,and panicle length were identified by BSA-seq(Table 1).All detected QTL were validated in F2or BC4F2populations,confirming reliable QTL detection.Although the ΔSNP-index peaks of these detected QTL exceeded the established threshold,the numerical values of peaks were not very large,especially these of major QTL(Table 2).Unlike conventional QTL mapping,the BSA strategy is based on allele-frequency differences between two contrasting pools consisting of plants with high and low phenotypic values[6,8].It is assumed that plants with extremely high phenotypic values carry homozygous alleles with additive effects,and vice versa.In general,only a few QTL for a target trait can be detected using the BSA-seq strategy in an F2population[14,62,63].If four QTL control trait variation(assuming no interaction),in principle,one of 256 plants is homozygous at all four QTL.If a pool is composed of 20 plants,at least 5120 plants should be grown in the field to select the expected genotypes.If the size of the population is inadequate,some selected plants in the pools will carry heterozygous alleles at one or two loci,greatly influencingthe difference in allele frequency between pools.Strict screening of extreme plants to construct bulks from a large segregating population is conducive to the aggregation of alleles with the same directional effect.Thus,an increased difference in allele frequency would be observed in the QTL regions between two contrasting bulks;theΔSNP-index in turn increases,which is conducive to QTL detection.

Table 4Comparison of agronomic traits between two homozygous qPH8-NILs.

qPH8was fine-mapped and its causal gene awaits future identification.BothqPH1.2andqPH8were major QTL for plant height and did not affect heading date,accounting for the low correlation between plant height and heading date in the F2population.qPH8not only controlled plant height but also influenced panicle length(Table 4).However,qPH8was not identified in PL bulks because theΔSNP-index peak in the region did not exceed the threshold(Fig.S3).The PVE of all these detected QTL was greater than 10% in the F2population(Table 2),a finding consistent with previous simulation results[14]indicating that the power of QTL detection was high when the QTL relative contribution to total phenotype variation was greater than 10%.qPH8explained only 5.9% of the variation in panicle length in the F2population(Table S7),a level too low for detection in PL bulks by QTLseqr.But its identity was validated in the BC4F2population and NILs(Table 4;Table S7).

4.2.Increasing the bulk size and sequencing depth improves the resolution of QTL mapping

The QTL intervals(95% confidence interval)mapped with QTLseqr averaged 2.5 Mb(Table 1),the approximate resolution of the QTL-seq method for QTL mapping[14].Advanced approaches of BSA-seq such as GradedPool-seq and QGT-seq can be used to narrow a QTL region and even to fine-map a QTL rapidly by improving the algorithm,adding additional pools,or optimizing the mapping population[63,64].GradedPool-seq allows high resolution(～400 kb)by classifying the mapping population into several graded bulks and then using a modified BSA with the Ridit analysis[63].The QGT-seq strategy combines QTL partitioning to convert a quantitative trait into a near-qualitative trait,sequencing of bulked-segregant pools from a large segregating population,and the use of the maximum-likelihood method for identifying candidate genes[64].Both methods require many bulked samples(20%–30% of individuals from the total population for each bulk)and sufficient sequencing depth to ensure high mapping resolution,which is largely influenced by recombination frequency and genotyping accuracy[63,64].However,in the present study,only 20 plants(2.1%)were sampled for each bulk,and the limited number of plants of each bulk were not sufficient for fine-mapping the target QTL.It is recommended to select as many individual plants with extreme phenotypes as possible for bulk construction and increase the sequencing depth,substantially increasing the resolution of BSA-seq for QTL mapping.

4.3.qPH8 is conducive to improving hybrid rice plant architecture

Application ofsd1in rice breeding for semidwarf stature triggered the first‘‘green revolution”,which led to substantial increases in rice production owing to increased lodging resistance and planting density[56].Most elite rice cultivars carry nonfunctional or weak alleles ofsd1[36].However,overexploitation of this single gene can reduce the genetic diversity of rice cultivars and increases the risk of pest and disease pressure[65].To address this problem,many genes associated with plant height have been cloned,but most of them are not suitable for breeding,owing to their undesirable effects.The effect ofqPH8on plant height was equivalent to that ofsd1in both the F2and BC4F2populations,and the C7 allele ofqPH8reduced plant height by 28.7 cm in the HBK background(Table 4).The decrease in plant height was due mainly to the shortening of the basal internodes,resulting in strong resistance to lodging(Fig.4E and F;Table 4).The PH was reduced byqPH8,whereas the yield per plant was not significantly reduced,suggesting the potential value of this QTL for rice breeding.In the process of semidwarf hybrid rice breeding,recessive dwarfing genes are introgressed into both parents.qPH8is an additive gene,and one parent carrying the C7 allele ofqPH8may confer semidwarf stature in hybrid rice with a midparent plant height value.

CRediT authorship contribution statement

Bo Zhangperformed most of the research and drafted the manuscript.Feixiang Qi and Jianghu Mengwere responsible for the QTL validation.Yikai Yang and Li Zhangperformed the phenotypic data collection and bulked DNA extractions.Gang Huconducted the SNP calling.Zhongmin Han,Xiangchun Zhou,Haiyang Liu and Mohammed Ayaadprovided technical assistance.Yongzhong Xingdesigned the experiment,supervised the study,and revised the manuscript.

Declaration of competing interest

Authors declare that there are no conflicts of interest.

Acknowledgments

We thank Mr.Jianbo Wang for his excellent work in the paddy fields.This study was partially supported by the National Natural Science Foundation of China(31701391)and the National Key Laboratory of Crop Genetic Improvement Self-Research Program(ZW18B0101).

Appendix A.Supplementary data

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