A single nucleotide deletion in the third exon of FT-D1 increases the spikelet number and delays heading date in wheat (Triticum aestivum L.).

The spikelet number and heading date are two crucial and correlated traits for yield in wheat. Here, a quantitative trait locus (QTL) analysis was conducted in F8 recombinant inbred lines (RILs) derived from crossing two common wheats with different spikelet numbers. A total of 15 stable QTL influencing total spikelet number (TSN) and heading date (HD) were detected. Notably, FT-D1, a well-known flowering time gene in wheat, was located within the finely mapped interval of a major QTL on 7DS (QTsn/Hd.cau-7D). A causal indel of one G in the third exon of FT-D1 was significantly associated with total spikelet number and heading date. Consistently, CRISPR/Cas9 mutant lines with homozygous mutations in FT-D1 displayed an increase in total spikelet number and heading date when compared with wild type. Moreover, one simple and robust marker developed according to the polymorphic site of FT-D1 revealed that this one G indel had been preferentially selected to adapt to different environments. Collectively, these data provide further insights into the genetic basis of spikelet number and heading date, and the diagnostic marker of FT-D1 will be useful for marker-assisted pyramiding in wheat breeding.

Introduction

Wheat is one of the largest grain crops in the world today. Continuing population growth is increasing the demand for the development of novel ways to enhance its yield potential (Boden et al., 2015). Wheat yield can be divided into three components, including thousand grain weight (TGW), grain number per spike (GNS) and spike number per hectare (Simmonds et al., 2014). Generally, GNS can be broken down into the spikelet number and grains per spikelet (Quarrie et al., 2006). An increased number of fertile florets per spikelet could probably increase in yield potential in wheat (Sakuma et al., 2019; Sakuma and Schnurbusch, 2020). Thereby, modifying the spikelet number could be a strategy to improve yield potential in wheat, since it is associated with a higher GNS (Lewis et al., 2008).

Spikelet number is quantitative in nature and controlled by numbers of genes/QTL (quantitative trait loci). Identification of DNA markers associated with spikelet number would allow marker‐assisted selection (MAS) and increase the efficiency for improving yield potential during breeding. Over the past two decades, the successful application of quantitative‐genetic methodology facilitated the identification of QTL for spikelet number (Cui et al., 2012). Many QTL across multiple genetic backgrounds have been identified on almost all 21 chromosomes in wheat (Chen et al., 2020b; Deng et al., 2011; Ma et al., 2007, 2019; Yao et al., 2019; Zhai et al., 2016; Zhang et al., 2015; Zhou et al., 2017). Only few genes involved in spikelet number have been characterized through map‐ or homology‐based cloning, such as TaMOC1, TaAPO1, TB1 and WFZP (Dixon et al., 2018b; Du et al., 2021; Muqaddasi et al., 2019; Zhang et al., 2015). Despite these efforts, cloning of the QTL controlling spikelet number will provide an important entry for studying the gene network involved in the spike development in wheat (Lewis et al., 2008).

Heading date is usually determined by vernalization requirement, photoperiodic response and earliness per se, which is crucial for wide adaptation to diverse environments and affecting crop yield (Gawroński et al., 2014; Lewis et al., 2008; Shimada et al., 2009; Yasuda and Shimoyama, 1965; Zhang et al., 2018). Vernalization is defined as the acquisition of the ability to flower by a low‐temperature treatment (Chouard, 1960; Fu et al., 2005). Genetic studies demonstrated that vernalization process is controlled by at least four critical genes designated Vrn‐1, Vrn‐2, Vrn‐3 and Vrn‐4 (Kippes et al., 2015; Yan et al., 2003, 2004, 2006). Allelic variations in Ppd‐1 influence sensitivity to daylength in temperate cereals (Boden et al., 2015; Gauley and Boden, 2021). Wheat varieties that are photoperiod sensitive required a period of long days to permit initiation of the floral primordia without undue delay (Worland, 1996). Earliness per se is described as the variation in flowering observed once photoperiod and vernalization requirements were fully satisfied (Gomez et al., 2014). The previous studies have shown a positive correlation between spikelet number and heading date (Boden et al., 2015; Lewis et al., 2008). Thus, identification and comparison of genomic regions governing the spikelet number and heading date will provide us useful information for genetic improvement of these two traits.

FLOWERING LOCUS T (FT) encoding a PEBP protein that could be transported from the leaf phloem to shoot apical meristem to initiate floral transition is a well‐conserved florigen gene across flowering plants (Tsuji and Taoka, 2014). The ortholog gene of Arabidopsis FT in wheat is FT‐1 (also named as Vrn‐3) (Yan et al., 2006). Growing evidence suggests that beyond flowering, FT‐1 also plays an important role in wheat spikelet development (Sakuma and Schnurbusch, 2020). For example, an absence of FT‐B1 was found an increase in spikelet number when grown under different temperature regimes (Dixon et al., 2018a; Finnegan et al., 2018). Brassac et al. (2021) identified a non‐synonymous mutation in FT‐B1 that mainly influenced TSN with a minor effect on heading date. FT‐A1 with variations in the promoter or exon might influence the spikelet number and heading date (Chen et al., 2020b; Yu et al., 2017). Although some studies have identified QTL associated with spikelet number and heading date in the flanking region of FT‐D1 (Chen et al., 2020b; Isham et al., 2021), further characterization of the genetic control of FT‐D1 in spikelet number using transgenic technology is still an area to be elucidated.

In this study, we reported the mapping of stable QTL for total spikelet number (TSN) and heading date (HD) using recombinant inbred lines (RILs) derived from a cross between Hesheng2hao (HS2) and Nongda4322 (4332). Of ten genomic regions harbouring 15 stable QTL, three had pleiotropic effects for total spikelet number and heading date. One major QTL for TSN and HD on chromosome 7DS was further validated. A causal G indel in the third exon of FT‐D1 located in the fine‐mapping interval was highly associated with these two traits. Furthermore, analysis of mutant lines developed by CRISPR/Cas9 provided genetic evidence for the function of FT‐D1 in regulating spikelet number and heading date.

Results

Phenotypic performance of the two parents and RIL population

Based on the phenotypic data of the two investigated traits collected from two parents in six environments, we found that 4332 had higher total spikelet number (TSN) than that of HS2 across all environments (Figure 1a,b). In addition, 4332 showed significantly delayed heading date (HD) in E4 and E5 environments when compared with HS2 (Figure 1c). Transgressive segregation was observed for each trait in the RIL population (Figure S1). TSN displayed normal distribution, whereas HD exhibited asymmetrical distribution (Figure S1). Both traits showed high broad‐sense heritability (; TSN, 0.95; HD, 0.96). Pairwise correlation analysis showed a positive correlation between TSN and HD (r = 0.472).

Figure 1

Spike morphologies and phenotypic data of HS2 and 4332. (a) Spikelet number (left) and rachis (right) of HS2 and 4332 grown in Beijing. The bars represent 1 cm. (b, c) Statistical analysis of heading date (b) and spikelet number (c) in HS2 and 4332. Values are means ± standard deviations (SD). Only one repeat was measured for HD at E1, E2 and E3 environments. E1, 2014–2015 Hebei; E2, 2014–2015 Shandong; E3, 2014–2015 Shanxi; E4, 2015–2016 Hebei; E5, 2015–2016 Shandong; E6, 2015–2016 Shanxi. Significant differences are indicated by * (P < 0.05), ** (P < 0.01), *** (P < 0.001), ns (no significant difference) (Student's t‐test).

Identification of genomic regions harbouring stable QTL for spikelet number and heading date

Quantitative trait locus that could be detected in ≥3 individual environments and the combined analysis (BLUP) are regarded as ‘stable QTL’. According to the criterion, 15 stable QTL for TSN and HD were identified within 10 genomic regions of nine chromosomes (1B, 2A, 3B, 3D, 5B, 5D, 6B, 7A and 7D) (Figure 2, Table S1). ‘Putative QTL’ that could not meet this criterion are also listed in Table S1 and shown in Figure 2.

Figure 2

Chromosomal and physical locations of QTL regions associated with total spikelet number per spike (TSN) and heading date (HD). Three centimorgan (cM) scales are shown on the left. Black ellipses on chromosome bars indicate the approximate positions of the centromeres. Vertical bars represent the confidence interval for the location of each QTL. Black and red triangles indicated that positive alleles of stable QTL were contributed by 4332 and HS2, respectively. Grey and pink arrows indicate that positive alleles of putative QTL were contributed by 4332 and HS2, respectively. Red bars on the physical map represent the approximate positions of the QTL identified in the present study. The known positions of Ppd‐A1 and Vrn‐B1 are presented in blue arrows.

Five stable QTL on chromosome arms 1BL, 3BL, 5DS and 7AL were found to have effects on TSN but with no significant effect on HD (Figure 2, Table S1). The QTL region on chromosome arm 1BL harboured a QTL (QTsn.cau‐1B.1), with 4332 contributing the allele associated with increased TSN. QTsn.cau‐3B.1 was flanked by the markers BS00037536_51 and Excalibur_c30527_559 on the chromosome arm 3BL, with the favourable allele contributed by 4332. In addition, HS2 conferred effects for increasing TSN at the QTsn.cau‐5D.1 locus. Two adjacent TSN QTL (QTsn.cau‐7A.1 and QTsn.cau‐7A.2) were identified on chromosome arm 7AL, with the superior alleles from 4332.

Four stable QTL regions on chromosome arms 5BL, 6BL and 7AS mainly influenced HD (Figure 2, Table S1). A stable QTL on chromosome arm 5BL was designated as QHd.cau‐5B, with the positive allele from 4332. QHd.cau‐6B was flanked by markers Ex_c17379_1431 and BS00034339_51 on 6BL, and HS2 conferred the allele for delayed HD. Two stable QTL, QHd.cau‐7A.1 and QHd.cau‐7A.2, located in an adjacent interval on chromosome 7AS, only contributed 1.75%–6.11% and 2.91%–5.29% of the phenotypic variation, respectively, in the detected environments.

The stable QTL regions on chromosome arms 2AL, 3DS and 7DS could simultaneously influence TSN and HD (Figure 2, Table S1). A major QTL for TSN (QTsn.cau‐2A.1) was detected on chromosome arm 2AL, which accounted for 12.50%–25.48% of the observed variation. In addition, QHd.cau‐2A was detected in this region flanked by markers CAP8_c607_659 and BS00022301_51. Two QTL (QTsn.cau‐3D and QHd.cau‐3D) on chromosome arm 3DS shared the uniform confidence interval, and the flanking markers were wsnp_Ex_c2258_4232538 and Ra_c4231_465. A QTL controlling TSN on chromosome 7DS (QTsn.cau‐7D) which explained 15.19%–30.93% of the phenotypic variation was identified over all environments. In addition, one major QTL for HD (QHd.cau‐7D) was co‐located with QTsn.cau‐7D, which explained 39.13%–54.32% of the phenotypic variation. For both regions on chromosome arms 2AL and 7DS, 4332 contributed the alleles for increased TSN but delayed HD. Different from that of chromosome arms 2AL and 7DS, 4332 provided the allele that was associated with delayed heading on chromosome 3DS, whereas HS2 conferred the increasing alleles for TSN.

Fine mapping of the QTL on chromosome 7DS

To delimit the genomic interval of QTL on the short arm of chromosome 7D (QTsn/Hd.cau‐7D), we used the linkage map of 7DS consisting of 23 markers that had been reported by Chen et al. (2020a) to conduct the first step of fine mapping. QTsn/Hd.cau‐7D was further mapped to the genomic region flanked by markers Xcau.7D‐3 and Xcau.7D‐5, and the corresponding physical interval was from 56.9 Mb to 72.5 Mb according to IWGSC (International Wheat Genome Sequencing Consortium) RefSeq v1.0 (IWGSC, 2018) (Figure 3a). Next, eight SSR markers were used to screen the segregation populations of RIL120, a residual heterozygous line that showed heterozygosity within the interval of QTsn/Hd.cau‐7D (Chen et al., 2020a). Consequently, eight recombinants, designated as NF1~NF8 (NF, used for deriving NIL families), were detected that had crossover breakpoints within the QTsn/Hd.cau‐7D region (Figure 3b). NF1~NF8 were self‐pollinated to produce corresponding segregation families, and the values for TSN and HD were compared between two homozygous groups in each NF segregation family (Figure 3b). No significant difference for TSN and HD was detected in NF4, whereas NF4332 (with 4332 haplotype) had 2.59%–6.25% higher TSN (P < 0.01) and headed later 0.7–2.4 days (P < 0.05) than those of NFHS2 (with HS2 haplotype) in NF1~3 and NF5~8 (Figure 3b). Collectively, QTsn/Hd.cau‐7D was delimited into the interval of approximately 6 Mb flanked by the markers Xcau.7D‐4 and Xcau.7D‐7 (Figure 3b).

Figure 3

Fine‐mapping of QTsn/Hd.cau‐7D. (a) The results of QTL mapping using genetic linkage map with 23 markers. ‘a’ and ‘R2’ represent additive effect and explained phenotypic variation rate, respectively. (b) Left side is the graphical genotypes of the recombinants (NF1~NF8). On the right side, the means of TSN and HD (mean ± SD) of the homozygous plants in NF families are presented. a‘No.’, number of plants in each NF family. b‘A’, Plants with HS2 alleles, c‘B’, Plants with 4332 alleles. (c) Further fine‐mapping of QTsn/Hd.cau‐7D using four new markers. Left side is the graphical genotypes of the recombinants. Right side is the comparisons of TSN and HD between NILHS2 and NIL4332 derived from corresponding recombinants at three environments. White, grey and black rectangles represent HS2, heterozygous and 4332 alleles, respectively. Significant differences are indicated by * (P < 0.05), ** (P < 0.01), *** (P < 0.001), **** (P < 0.0001) (Student's t‐test).

To further map the gene(s) associated with TSN and HD, BC5F2 population was developed by backcrossing RIL222 to 4332 (as the recurrent parent). Four polymorphic InDel markers (7D‐ID‐6, 7D‐ID‐33, 7D‐ID‐8 and 7D‐ID‐12) were developed and anchored to the interval between Xcau.7D‐4 and Xcau.7D‐7 according to the physical positions (Figure 3c, Tables S2 and S3). Using these new markers, BC5F2 plants showing heterozygosity within the interval Xcau.7D‐4~7D‐ID‐12 were self‐pollinated to conduct the next fine‐mapping step. Ultimately, nine recombinants were screened from different populations, including four from BC5F3 (BF3‐1 to BF3‐4), four from BC5F4 (BF5‐1 to BF5‐4) and one from NF3 (NF3‐1) (Figure 3c). Subsequently, nine corresponding sets of NILs derived from these nine recombinants were planted in three locations (Table S4). Significant differences in TSN and HD were observed between NILHS2 and NIL4332 derived from BF3‐1~4 and BF5‐1~3 across all environments (P < 0.001) (Figure 3c). However, there was no significant difference for TSN and HD between NILHS2 and NIL4332 derived from BF5‐4 and NF3‐1 in at least two locations (Figure 3c). Based on these results, QTsn/Hd.cau‐7D was narrowed down to 2.73 Mb physical region flanked by markers 7D‐ID‐6 and 7D‐ID‐9 (Figure 3).

To estimate the additive and dominance effects of QTsn/Hd.cau‐7D, the TSN and HD of NF6 family were selected for further analysis, including 145 heterozygous plants and 160 homozygous plants with alternative haplotypes across the interval from Xcau.7D‐2 to Xcau.7D‐7. Additive effects were observed for both TSN (d/a = −0.05) and HD (d/a = 0.18) (Table S5).

Cloning of the candidate gene FT‐D1

Based on the wheat gene annotation database (IWGSC, 2018), the finely mapped genomic interval of QTsn/Hd.cau‐7D encompassed 34 high‐confidence genes, including TraesCS7D01G111600, a well‐known FT gene controlling flowering (Yan et al., 2006) (Table S6). To anchor the candidate genes, we firstly analysed the expression patterns of the 34 genes on the Wheat eFP Browser (http://bar.utoronto.ca/efp_wheat/cgi‐bin/efpWeb.cgi; Ramirez‐Gonzalez et al., 2018). A total of 11 genes were found to be mainly expressed in shoot apical meristem, leaf and/or spike, which will most probably influence heading date and/or spikelet development (Table S7). Moreover, we investigated the sequence variations of these 11 genes in two parents according to the resequencing data. An indel of one G was detected between HS2 and 4332 within the third exon (+840 bp downstream of the translation start codon ATG) of TraesCS7D01G111600 (FT‐D1) (Figure 4a, Table S6). There was no other variation between two parents by amplifying a ~5.5 kb genomic sequence of FT‐D1 including 3.3 kb‐upstream sequence, exons, introns and 1.2 kb‐downstream sequence from HS2 and 4332 (Figure S2). The genotypes of HS2 and 4332 were named FT‐D1(G) and FT‐D1(ΔG), respectively (Figure 4a). As a consequence, the deletion of one G in 4332 which occurred outside the PEBP domain led to a frameshift mutation of amino acid sequence (Figure S3), revealing a new lipoprotein attachment site with two disulphide bonding sites (Bonnin et al., 2008). Based on the polymorphism of FT‐D1 between HS2 and 4332, a semi‐thermal asymmetric reverse PCR (STARP) marker (SFT‐D1) was developed to genotype the RILs (Figure 4b). Notably, the homozygous and heterozygous genotypes can be easily identified by the SFT‐D1 marker (Figure 4c). As expected, SFT‐D1 was located to the finely mapped interval of QTsn/Hd.cau‐7D (Table S3). Moreover, SFT‐D1 was used to genotype the seven NF families (NF1~3 and NF5~8). The result of t‐test indicated that FT‐D1 was associated with delayed HD and increased TSN (Figure 4d).

Figure 4

Cloning of FT‐D1 and development of STARP marker. (a) Schematic diagram of nucleotide polymorphism for FT‐D1. The polymorphic site and relative positions are indicated on the genomic sequence of FT‐D1. ‘−’ and ‘Δ’ represent deletion. Exons are indicated by black boxes, flanking regions and introns are indicated by black solid lines. FT‐D1(ΔG) and FT‐D1(G) represent 4332 and HS2 alleles, respectively. (b) STARP marker was designed according to the nucleotide polymorphism of FT‐D1. The red hexagons and blue solid rectangles indicate nucleotides (SNP) and substituted nucleotides, respectively. SFT‐D1.F1 and SFT‐D1.F2 specially amplify FT‐D1(G) and FT‐D1(ΔG), respectively. TGCTGACGAC in red indicates the introduced insertion in SFT‐D1.F2. ‘‐’ means deletion. (c) PCR products of HS2, 4332, Ae. tauschii (DD) and hexaploid (AABBDD) wheat accessions and progenies of recombinants amplified by STARP marker. (d) Effects of FT‐D1 in populations derived from NF1~NF3 and NF5~NF8 using SFT‐D1. The number of plants used for test was shown on the corresponding boxes. Significant differences are indicated by * (P < 0.05), ** (P < 0.01), *** (P < 0.001), **** (P < 0.0001) (Student's t‐test).

Phenotypic analyses of homozygous and transgene‐free ft‐D1 mutants in wheat

To determine whether FT‐D1 was associated with the TSN and HD in wheat, CRISPR/Cas9‐mediated gene editing was performed to knock out FT‐D1 in a photoperiod insensitive cultivar ‘CB037’ background, in which the genotype of FT‐D1 was the same as that in HS2 (Figures S4 and S5). Sequence analysis showed that all mutations occurred in the sgRNA‐2 target site (target‐2) in the first exon (Figure 5a). After that, three homozygous ft‐D1 mutant lines containing mutations only in FT‐D1 were obtained and subsequently designated as ko‐1, ko‐2 and ko‐3 (Figure 5a, Figure S6). Specifically, ko‐1 was mutated by the deletion of a guanine nucleotide (G), ko‐2 was disrupted by the deletion of two G, and ko‐3 carried seven nucleotide deletions in D genome (Figure 5a, Table S8). Protein sequence analysis revealed that those mutations resulted in amino acid changes and premature termination (Figure S7). Furthermore, these three lines were confirmed to be transgene‐free by genomic PCR with a primer set recognizing the sequence of the CRISPR/Cas9 construct (Figure S8). Field phenotypic analysis showed that the heading date of all three mutant lines was 2–3 days later than that of CB037, and there was no significant difference between CB037 and negative control (NC) (Figure 5b,e). Further analysis revealed that the spikelet number of the three ft‐D1 mutant lines was significantly higher than that of wild type (WT, CB037) (P < 0.01) (Figure 5c,f). In addition, we observed that the plant height and spike length in ft‐D1 mutants were all significantly elevated compared with that of CB037 (Figure 5c,d,g,h), whereas only ko‐3 showed significantly higher grain number per spike (GNS) compared with that of CB037 (Figure 5i). Taken together, these results revealed that FT‐D1 gene played an important role in the regulation of spikelet number and heading date.

Figure 5

Phenotype of FT‐D1 knock‐out mutant lines induced by CRISPR/Cas9. (a) The gene structure of FT‐1, locations of the two target sites of FT‐D1 and mutant types of three ft‐D1 mutant lines induced by CRISPR/Cas9. Deleted nucleotides are shown by ‘.’. The numbers represent the number of nucleotides involved in the deletion. (b–d) Phenotypes of heading (b), spike (c) and plant height (d) of CB037, negative control (NC) and ft‐D1 mutant lines. Scale bars in (b) and (d) represent 10 cm; bars in (c) represent 1 cm. (e–i) Statistical analysis of phenotype between CB037 and ft‐D1 mutant lines. (e) heading date; (f) spikelet number; (g) plant height; (h) spike length; (i) grain number per spike. Data were means ± SD, and t‐test was used to compare mutant lines, NC and CB037. *P < 0.05, **P < 0.01, ***P < 0.001; ****P < 0.0001; ns, no significant difference.

Association analysis and implication in breeding of FT‐D1

In further exploration of the evolution and domestication of FT‐D1, the marker SFT‐D1 was used to genotype a diversity panel of wheat accessions, including 1430 common wheat (BBAADD) and 88 Aegilops tauschii (DD) from different countries (Tables S9–S11). All Ae. tauschii presented the FT‐D1(G) allele, indicating that the deletion of one G in FT‐D1 might occur after the formation of allohexaploid species (Figure 4c, Table S9). To determine the distribution of the FT‐D1(G) or FT‐D1(ΔG) allele in relation to the origin of geographic area, a total of 1412 hexaploid wheat accessions with known origin were used for analysis (Tables S10 and S11). In China, FT‐D1(G) was the prevalent allele in Yellow and Huai River Valleys Facultative Wheat Zone (92.87%) and Northern Winter Wheat Zone (86.1%) (Figure 6a, Table S10). By contrast, up to 84.0% of the common wheat accessions from foreign countries carried the FT‐D1(ΔG) allele (Figure 6a, Table S11). In addition, the FT‐D1(ΔG) allele occurred at a higher frequency (62.3%–67.7%) in Xinjiang Winter‐Spring Wheat Zone, Southwestern Autumn‐Sown Spring Wheat Zone, and Middle and Lower Yangtze Valleys Autumn‐Sown Spring Wheat Zone (Figure 6a, Table S10).

Figure 6

Geographic distributions and contributions of FT‐D1 haplotypes based on the functional STARP marker SFT‐D1. (a) Geographic distributions of FT‐D1 haplotypes in the worldwide and Chinese accessions. I, Northern Winter Wheat Zone; II, Yellow and Huai River Valleys Facultative Wheat Zone; III, Middle and Lower Yangtze Valleys Autumn‐Sown Spring Wheat Zone; IV, Southwestern Autumn‐Sown Spring Wheat Zone; V, Southern Autumn‐Sown Spring Wheat Zone; VI, Northeastern Spring Wheat Zone; VII, Northern Spring Wheat Zone; VIII, Northwestern Spring Wheat Zone; IX, Qinghai‐Tibetan Plateau Spring‐Winter Wheat Zone; X, Xinjiang Winter‐Spring Wheat Zone. (b) Contributions of FT‐D1 to TSN and HD in 150 wheat accessions grown in three different environments using SFT‐D1 marker. 15LF, 2014–2015 Shanxi Linfen; 16SY, 2015–2016 Shaanxi Sanyuan; 16HB, 2015–2016 Hebei Shijiazhuang. *P < 0.05, **P < 0.01, ***P < 0.001; ****P < 0.0001.

To investigate the implications of our findings on wheat breeding, we analysed the effects of the FT‐D1(ΔG) and FT‐D1(G) allele on TSN and HD in a diversity panel of 150 accessions in three environments (Figure 6b, Table S12). As expected, the accessions with FT‐D1(ΔG) allele had more TSN, but with a delayed heading date than those with the FT‐D1(G) allele (P < 0.01) across all three environments (Figure 6b, Table S12). These results indicated that FT‐D1 was conserved in the hexaploid wheat and FT‐D1(ΔG) allele was effective for TSN and HD in different genetic backgrounds.

Discussion

Pleiotropic QTL for spikelet number and heading date on chromosome arms 2AL, 3DS and 7DS

Spike morphogenesis in wheat is subdivided into spike initiation and growth phases (Sreenivasulu and Schnurbusch, 2012). Theoretically, the more spikelet number could be attributed to a longer thermal duration for spikelet primordia production and/or a faster rate of spikelet production, whereas the former was more closely associated with the increased spikelet production than the latter (Gaju et al., 2009; Guo et al., 2018a). Consistent with this, genes associated with heading date are always employed in manipulation of spikelet number in wheat. For example, photoperiod insensitive alleles of Ppd‐1 brought forward the time of terminal spikelet, and hence reduced the spikelet number and GNS (Guo et al., 2018b; Snape et al., 2001). VRN1 and/or FUL2 mutations delayed flowering initiation with highly significant effects on spikelet number (Li et al., 2019). The effect of Eps‐Am1 is to shorten the vegetative and spike initiation phases, thereby producing fewer spikelets (Lewis et al., 2008).

Here, we also found that heading date was significantly positively correlated with TSN (r = 0.472, P < 0.01). Moreover, three of seven stable QTL for heading date (QHd.cau‐2A, QHd.cau‐3D and QHd.cau‐7D) were associated with TSN (Figure 2, Table S1), providing further genetic evidence that the heading date co‐regulated with spikelet differentiation and the maximum spikelet number per spike determination (Abeledo et al., 2002). The superior alleles for increasing TSN and HD in the QTL regions on chromosome arms 2AL and 7DS were provided by 4332. Therefore, selection of these QTL for a higher spikelet number would be inevitably accompanied by a delay of heading date and vice versa. However, in the QTL region on chromosome 3D, HS2 carried the positive allele for increasing TSN as well as advancing HD. These results indicated that the mechanism of the QTL region on chromosome 3D controlling HD and spikelet number is different from that of the QTL regions on chromosomes 2A and 7D. Collectively, we proposed that these QTL could be attractive targets for marker‐assisted selection to develop high spikelet number or early flowering varieties.

A single nucleotide deletion of FT‐D1 contributed to increased spikelet number and delayed heading date

Here, we demonstrated that a single nucleotide deletion of FT‐D1 would be the causal mutation of increased spikelet number but delayed heading date. This single nucleotide (G) deletion in the third exon of FT‐D1 was detected in 4332, which led to a frameshift mutation of amino acid sequence (Figure 4a, Figure S3). Compared with previous studies, we found that this variation was relatively conservative in different wheat materials. For example, Bonnin et al. (2008) analysed the nucleotide polymorphisms of the D copy FT in wheat and observed an insertion–deletion of one G in the third exon. The 1‐bp deletion of G was also found in 2174 but not in Jagger by Chen et al. (2010). In molecular breeding, functional markers are very useful for enhancing the precision and accuracy in marker‐assisted selection (MAS) of the target gene (Collard and Mackill, 2008). Hence, the development of markers based on the causal single nucleotide deletion of FT‐D1 can directly differentiate the alleles conferring increased spikelet number but delayed heading date. Recently, a new genotyping method named semi‐thermal asymmetric reverse PCR (STARP) with advantages of high accuracy, simple assay design, low operational costs and compatibility for a variety of platforms, was developed by Long et al. (2017). Based on 1‐bp polymorphism, the co‐dominant STARP marker SFT‐D1 was developed and mapped to the interval between Xcau.7D‐8 and 7D‐ID‐8 on chromosome 7D (Table S3). Collectively, this simple, robust and economical STARP marker will simplify and streamline MAS for spikelet number and heading date in wheat breeding.

Diverse functions of FT and FT‐like genes

FT is known as a florigen gene promoting flowering in flowering plants. In the present study, we demonstrated the function of FT‐D1 involved in spikelet number. Recent studies showed that FT and FT‐like genes had diverse roles in plant development beyond flowering and spikelet number. The result of Kinoshita et al. (2011) defined a new cell‐autonomous role of FT in regulating stomatal opening. Tsuji et al. (2015) showed that Hd3a (FT homolog) protein accumulates in axillary meristems to promote branching in rice. AcFT1 promotes bulb formation of onion, which adds to the growing evidence that FT genes play a wide role in controlling development decisions (Lee et al., 2013). The ft2aft5a mutants produced significantly increased numbers of pods and seeds per plant compared with WT in soya bean (Cai et al., 2020). In wheat, Shaw et al. (2019) reported that FT2 (the closest paralogue of FT‐1) played an important role in spike development and fertility. A QTL cluster on 7DS harbouring FT‐D1 was identified to be associated with grain yield, fertile spikelet number per spike and thousand kernel weight (Isham et al., 2021). Here, we demonstrated the multiplex roles of FT‐D1 not only on TSN and HD but also on PHT and SL using the ft‐D1 mutants, which were probably associated with the function of FT‐1 in regulating gibberellic acid (GA) genes during the spike development (Pearce et al., 2013). Previous studies also reported that the delayed heading allele of RFT1 (RICE FLOWERING LOCUS T1) in rice increased plant height (Zhu et al., 2017). It is not clear whether FT is transported to organs outside the shoot apex (Lee et al., 2013); therefore, the new molecular mechanism of FT‐D1 controlling different traits is still worth to elucidation.

The domestication and potential implication value of FT‐D1 in wheat breeding

The spikelet is the basal unit of inflorescence and is crucial for reproduction and final yield (Cai et al., 2014). Considering the influence of FT‐D1 on TSN and HD, the domestication and geographic distribution of FT‐D1 alleles will be of great interest for genetic improvement in wheat breeding. FT‐D1(ΔG) allele was not detected in the diploid progenitors (Ae. tauschii, DD) of allohexaploid wheat (BBAADD), suggesting that this allele may arise during the hexaploidization event that yielded hexaploid wheat or the process of domestication (Table S9). In China, FT‐D1(G) was the prevalent allele in Yellow and Huai River Valleys Facultative Wheat Zone (92.87%) and Northern Winter Wheat Zone (86.1%), which may be partially attributed to the cropping system that entails growing three crops in 2 years (He, 2001); that is, wheat varieties with early maturity were required to ensure planting of summer crops. By contrast, up to 84.0% common wheat accessions from foreign countries carried FT‐D1(ΔG) allele. We speculate that this variant is selected due to the its favoured potential and positive impact on the grain yield (Brassac et al., 2021). Taken together, these data indicated that the two FT‐D1 alleles have been widely used in wheat breeding programmes to enable adaption of wheat to special environmental conditions.

The rational design of yield and quality traits is a powerful strategy for meeting the challenges of future crop breeding (Zeng et al., 2017). Therefore, further isolation of the genes related to grain yield potential is needed to improve the wheat yield (Sakuma and Schnurbusch, 2020). Here, we used the CRISPR/Cas9 technology to explore the role of FT‐D1 on TSN and HD. Compared with WT, ft‐D1 mutants showed higher spikelet number with minor delayed heading date. Notably, significant GNS difference was observed between ko‐3 and WT. GNS is largely determined by the spikelet number and fertile florets within a spikelet (Quarrie et al., 2006; Sreenivasulu and Schnurbusch, 2012), and thus, the unexpected sterile spikelets of ko‐1 and ko‐2 possibly explained the mild significant difference in GNS from WT. Even so, the average of GNS of all mutants was higher than that of WT and the negative control (NC) (Figure 5i). These results indicated that ft‐D1 might have a potential in increasing grain yield. Thus, pyramiding ft‐D1 with other genes controlling different complex traits involved in yield could be a possible approach in improving wheat yield.

Experimental procedures

Plant materials and field trials

The mapping population containing 271 RILs was bred from a cross between Hesheng2hao (HS2) and Nongda4332 (4332) (Chen et al., 2020a). 4332 is a common germplasm resource characterized by high spikelet number, whereas HS2 is a high‐yield line with relatively fewer spikelet numbers. The RIL population and parents were planted in a randomized complete block design with three replications at three locations (Shandong, Hebei, Shanxi) during the 2014–2015 and 2015–2016 seasons, providing data for six environments (Table S4). Each line of RILs, along with the two parents, was evenly sown in two‐row plots (1.5 m long and 0.3 m apart) at a sowing rate of 30 seeds per row.

In addition, a total of 150 winter wheat accessions were selected to plant in single‐row plots with three replications at the following three environments: Shanxi Linfen (2014–2015 season), Hebei Shijiazhuang (2015–2016 season) and Shaanxi Sanyuan (2015–2016 season). Segregation families derived from residual heterozygous line (RIL120) were sown in rows 1.5 m long and 0.3 m apart at a sowing rate of 20 seeds per row at Hebei during the 2016–2017 growing season. NILs derived from BC5F3 were planted during 2019–2020 season in Hebei Jize, Shanxi Linfen and Shandong Qingdao; NILs derived from BC5F4 and NF3 were planted during 2020–2021 season in Hebei Jize, Shanxi Linfen and Henan Xinxiang (Table S4). The planting pattern of all NILs was same as 150 winter wheat accessions.

The wheat accessions, which included 1430 hexaploid wheat (967 Chinese accessions and 463 foreign accessions) (BBAADD), and 88 Ae. tauschii (DD) from different countries were gathered to detect the allele distribution for the target gene (Tables S9–S11).

Phenotypic evaluation in field conditions

For the RIL population, 150 winter wheat accessions and NILs, the heading date (HD) of each field plot was recorded as the developmental stage, by counting days from the sowing date to the date when approximately 50% of the spikes had fully emerged from the flag leaf sheath in each genotype. For RIL population planted at Shandong, Hebei and Shanxi during the 2014–2015, only one repeat was measured for HD. For 150 winter wheat accessions and NIL populations planted in one‐row plots, the spikelet number data of each genotype were collected from five spikes of main tillers in each replicate at maturity; for RIL population planted in two‐row plots, ten spikes of each genotype were measured. The resultant TSN data were the average ± SD (standard deviations) of three replicates for each environment.

For the fine‐mapping families NF1~NF8, the HD of the segregation families was scored when the first head of single plant was fully visible. At maturity, the main spike of each single plant was sampled, and the TSN was assessed.

Statistical analysis

The best linear unbiased prediction (BLUP) for TSN and HD across all environments was computed using SAS9.2 (SAS Institute Inc., Cary, NC) with the PROC MIXED procedure to reduce the influence of environmental effects on phenotypic variation. The descriptive statistics of RILs from each environment including the means and Pearson's correlation coefficient were analysed using software SPSS version 20.0 (SPSS, Chicago, IL). The Shapiro–Wilk test was performed with R software (V. 3.2.2; R Core Team, 2019) to test departures from normal distribution. The broad‐sense heritability () for each trait was estimated using the lme4 package in R software (Bates et al., 2015), and the formula is , where is genotypic effect, σ 2 is the residual error, n is the number of environments (Xu et al., 2017). Significance analysis was calculated using Student's t‐test.

QTL analysis

The genetic linkage map of HS2/4332 RILs used in the present study had been described in Chen et al. (2020a). The means of each trait under six environments and the adjusted mean values of BLUP were collected for QTL analysis. QTL analysis was performed using the Windows QTL Cartographer version 2.5 software (Wang et al., 2012), with the composite interval mapping (CIM) method (Zeng, 1994). The parameter exploited for the QTL calculations was set according to Chen et al. (2020a). QTL with threshold LOD values ≥2.5 and similar confidence intervals (±2 LOD away from the peaks of likelihood ratios) were considered identical QTL and were named according to McIntosh et al. (2017).

InDel markers' development

The resequencing data of HS2 and 4332 had been submitted by Cheng et al. (2020) to the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA533588. Insertion/deletions within the interval between markers Xcau.7D‐4 and Xcau.7D‐7 were searched according to the resequencing data of HS2 and 4332. Then, the corresponding sequences were blasted and used to design primers in Primer3 v.0.4.0 (http://bioinfo.ut.ee/primer3‐0.4.0/). The sequences of new InDel markers are listed in Table S2. The 10 μL PCR system followed as 5 μL 2×Taq PCR StarMix, 2 μL primer, 2 μL 50–100 ng/μL DNA template and 1 μL H2O.

Fine mapping of QTsn/Hd.cau‐7D

We adopted a recombinant‐derived progeny testing strategy to fine map QTsn/Hd.cau‐7D. Briefly, we used segregation families derived from two generations after self‐pollinated of RIL120, a residual heterozygous line that has been identified and used in our previous study (Chen et al., 2020a), to screen recombinants. Then, each recombinant was self‐pollinated and the derived‐homozygous lines were phenotypically evaluated to narrow down the target interval. In addition, RIL222 harboured HS2 alleles within the QTsn/Hd.cau‐7D was backcrossed to 4332 (4332 as recurrent parent) to obtain BC5F2 population. Using four new markers, four, four and one recombinants were screened from BC5F3, BC5F4 and progenies of NF3, respectively. Subsequently, 10–15 lines with homozygous HS2 alleles and 10–15 lines with homozygous 4332 alleles (NILs) were randomly selected from corresponding generations of nine recombinants, which were planted in three locations to evaluate phenotype (Table S4).

For the QTL controlling TSN and HD on chromosome arm 7DS, NF6 derived from RIL120 were divided into three genotypic classes (homozygous families and heterozygous families) based on their corresponding marker genotypes. The inheritance mode of an individual QTL can be classified into four categories: additive (d/a ≤ 0.20), partial dominance (0.20 < d/a < 0.80), dominance (0.80 ≤ d/a < 1.20) and overdominance (d/a ≥ 1.20), as described in the previous studies (Li et al., 2017; Zhai et al., 2018).

Cloning the genomic DNA of FT‐D1

Based on the sequence alignment, three set of primers for FT‐D1 (TraesCS7D01G111600) were designed to amplify HS2 and 4332 (Table S2). PCR products were gel‐purified and cloned into the pGEM‐T Easy Vector (Promega) for sequencing. Nucleotide sequence alignments of parents were performed using DNAMAN 8 (8.0.8.789).

STARP marker development

For distinguishing an insertion–deletion of one G between HS2 and 4332, we developed semi‐thermal asymmetric reverse PCR markers (STARP), as described by Long et al. (2017) with minor modifications. STARP marker included two asymmetrically modified AMAS primers and their same reverse primer. AMAS‐primer 1 (SFT‐D1.F1) was designed to uniquely amplify the FT‐D1(G) allele. Meanwhile, the AMAS‐primer 2 (SFT‐D1.F2), with its introduced 10 bp insertion (TGCTGACGAC) at 5′ terminus, could specially amplify the FT‐D1(ΔG) allele. Nucleotide substitution principle was followed as depicted in Long et al. (2017). Two AMAS primers and their reverse primer were mixed in a proportion of 1:1:2 and diluted. The PCR program was 95 °C for 5 min, followed by 9 cycles of 3‐step touchdown PCR program: 95 °C for 30 s, 68 °C for 30s and 72 °C for 30 s, with the annealing temperature being decreased by 1 °C per cycle and then 30 cycles of the 3‐step PCR program as usual. To increase resolution, a 10% non‐denatured polyacrylamide gel electrophoresis (PAGE) was adopted to separate the length polymorphism.

Construction of CRISPR/Cas9 transgenic vector

Two target sites were designed to recognize regions in the promoter sequence and coding sequence of FT‐1 on the website E‐CRISP Design (http://www.e‐crisp.org/E‐CRISP/), respectively (Table S2). We blasted these two target sequences on the Ensembl Plants database (Howe et al., 2019; http://plants.ensembl.org/index.html) and found that the target‐1 matched FT‐D1 perfectly but mismatched with FT‐A1 and FT‐B1, whereas target‐2 targeted three homoeologous genes of FT‐1 (Table S2). Fragments were amplified from the pCBC‐MT1T2 plasmid using four primers containing two target sites of sgRNA with two BsaI restriction site. The PCR products were digested with BsaI and inserted into BsaI‐digested pBUE414 to yield the fused expression vector. The resulting CRISPR/Cas9 construct was bombarded into the cultured immature embryo of CB037, a spring, photoperiod insensitive (carried a 2089 bp deletion in Ppd‐D1 promoter) wheat cultivar with high transformability (Figure S4) (Beales et al., 2007; Wang et al., 2017).

Homozygous mutant line development and phenotypic evaluation

T0 mutants were identified via amplifying the genomic DNA using a set of conserved primer, which were planted in a greenhouse under LD photoperiod and controlled temperature (24 °C during 16‐h light period and 20 °C during 8‐h dark period). Because target‐2 could target three homoeologous genes of FT‐1, primer pairs c2‐A1, c2‐B1 and c2‐D1 specific for FT‐A1, FT‐B1 and FT‐D1, respectively, were developed to analyse the mutant type of each transgenic plants. For T0 and T1 generations, the PCR products of each transgenic plant were cloned into pEASY‐T1 simple cloning vector (TransGen Biotech, Beijing, China), and the number of randomly sequenced clones was at least 30 (T0) and 6 (T1). The PCR products of T2 and T3 generations were sequenced directly. In addition, a marker (CR‐VT) was designed for detecting CRISPR/Cas9 (Table S2). Homozygous and transgene‐free lines (ft‐D1) were obtained by self‐pollinating of T0 twice. The detailed information for the mutant genotypes of ft‐1 in T0 to T3 generations is listed in Table S8. CB037 (WT), ft‐D1 (T4) and negative control (NC) were planted in the field of Beijing. A randomized complete block experiment with three replications was set, each replication of each material contained two rows. Before harvest, we measured plant height from eight plants in each row. Similarly, the data of spike length and spikelet number were collected from eight main spikes in each row. Ultimately, the mean of the heading date, plant height, spike length and spikelet numbers were collected from three replications. Grain number per spike was calculated by counting that of ten spikes divided by 10 after harvest.

Conflict of interest

The authors declare no conflict of interests.

Author contributions

ZN conceived the project; XC developed the RIL population; ZC, LC, XC and DD performed field trail of RILs and QTL analysis; ZC, WK, FH, HX and YZ participated in fine‐mapping and transgenic experiment; ZC, XC, LC and XW collected data for 150 wheat accessions used to analyse the contribution of FT‐D1; QS, XC, JX, MX, WG, ZH, ZS, JL, HP and YY assisted in revising the manuscript; ZC and WK analysed experimental results; ZC and ZN wrote the manuscript.

Figure S1. Histograms of the HS2/4332 recombinant inbred population for total spikelet number (TSN) and heading date (HD) based on the adjusted mean values (BLUP) collected from all environments.

Figure S2. Comparation analysis of FT‐D1 sequence between HS2 and 4332.

Figure S3. FT‐D1 peptide alignment between HS2 and 4332.

Figure S4. Comparation analysis of Ppd‐D1 promoter between Chinese Spring (CS) and CB037.

Figure S5. The gene sequence of FT‐D1 in CB037.

Figure S6. Comparison of whole gene sequences of FT‐A1 and FT‐B1 among CS (Chinese Spring), CB037 and three ft‐D1 mutants.

Figure S7. Peptides alignment between CB037 (7D) and ft‐D1 mutant lines.

Figure S8. Three transgene‐free lines were identified by used a pair of primer recognizing the sequence of the CRISPR/Cas9 construct.

Click here for additional data file.

Table S1. Effects of stable and putative QTL for TSN and HD in individual environments.

Table S2. Sequence of primers used in this study.

Table S3. Physical location of markers on the linkage map of 7DS.

Table S4. Planting environments and planting date of materials.

Table S5. Estimation of the additive and dominance effects of QTsn/Hd.cau‐7D on TSN and HD using NF6.

Table S6. Putative genes harbored in the interval between markers 7D‐ID‐6 and 7D‐ID‐9.

Table S7. Expression of 11 high confidence genes according to the gene expression atlas on the Wheat eFP Browser.

Table S8. The genotypes of ft‐1 mutant lines in T0 to T3 generations.

Table S9. FT‐D1 genotype in 88 Ae. Tauschii.

Table S10. Origins and genotypes of FT‐D1 in Chinese wheat accessions.

Table S11. FT‐D1 genotype in foreign wheat accessions.

Table S12. Genotype and phenotype of 150 wheat accessions used to analyze the contribution of FT‐D1.

Click here for additional data file.