Mapping QTL for spike fertility and related traits in two doubled haploid wheat (Triticum aestivum L.) populations

Background In breeding programs, the selection of cultivars with the highest yield potential consisted in the selection of the yield per se, which resulted in cultivars with higher grains per spike (GN) and occasionally increased grain weight (GW) (main numerical components of the yield). In this study, quantitative trait loci (QTL) for GW, GN and spike fertility traits related to GN determination were mapped using two doubled haploid (DH) populations (Baguette Premium 11 × BioINTA 2002 and Baguette 19 × BioINTA 2002). Results In total 305 QTL were identified for 14 traits, out of which 12 QTL were identified in more than three environments and explained more than 10% of the phenotypic variation in at least one environment. Eight hotspot regions were detected on chromosomes 1A, 2B, 3A, 5A, 5B, 7A and 7B in which at least two major and stable QTL sheared confidence intervals. QTL on two of these regions (R5A.1 and R5A.2) have previously been described, but the other six regions are novel. Conclusions Based on the pleiotropic analysis within a robust physiological model we conclude that two hotspot genomic regions (R5A.1 and R5A.2) together with the QGW.perg-6B are of high relevance to be used in marker assisted selection in order to improve the spike yield potential. All the QTL identified for the spike related traits are the first step to search for their candidate genes, which will allow their better manipulation in the future. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-03061-y.


Background
Wheat (Triticum aestivum L.) is one of the most cultivated and consumed worldwide cereals. Its production has to increase to supply the growing world population demand [1][2][3]. Improving the cultivar's yield potential by breeding (i.e., the yield of a cultivar adapted to the environment, which is growing without water or nutrient deficits and with no biotic stress, [4]) is a sustainable alternative to guarantee increases in world production [5,6]. Wheat breeding of yield potential has been based on empirical selection of yield per se due to the complexity of the character, the scarcity of knowledge and the lack of useful tools with real applicability in breeding programs [7]. This selection generally resulted in more grains per spike (GN), and hence, increased grains per unit area (no consistent trend in spikes per unit area has been reported) [8][9][10][11][12][13][14][15]. The grain weight (GW) has not shown any change with breeding; except for some recent reports in which yield potential was positively associated with its increment [16][17][18]. The selection process could be more efficient using molecular markers. The implementation of single nucleotide polymorphism (SNP) markers in plant breeding has increased the pace and precision of plant genetic analysis, which in turn, facilitated the implementation of crop molecular improvement [19]. SNP markers have been increasingly used for QTL mapping studies as they are the most frequent variations in the genome, and they provide a high map resolution [19,20,21]. Therefore, the identification, understanding and incorporation of yield related QTL could be a useful selection tool for a breeding program.
The most common approach, looking for genetic bases to further improving yield potential, is based on the numerical component analyses, GN and GW (see references quoted in Table S1). The GN is understood as the result of the total spikelets per spike (TS) and the grains per spikelet, being the former associated with the spike length (SL) and compactness of the spike (CN). Several QTL have been reported for the GW and the GN itself and their numerical components during the last years. Many studies identified stables QTL for these traits widespread in the genome (see Table S1). However, considering the IWGSC Ref. Seq. v1.0 wheat genome assembly [22], we identified QTL reported for the same trait that are located at the same position (Table 1). For example, a QTL for GW was detected in 6 studies on chromosome 7A-within 664.3-719. 6 Mb [23][24][25][26][27][28] (Table 1). Two important QTL for SL were detected on chromosomes 2D and 5A [26,[29][30][31][32][33][34][35][36][37] ( Table 1). Additionally, two QTL for TS were detected on chromosomes 5A and 7A in several studies [27, 30-32, 36, 38-41] (Table 1).
From the crop physiology approach, the GN depends on the florets that reach the fertile stage at anthesis (fertile florets per spike, FF) and on the proportion of them that sets grains (grain set, GST, grains per fertile floret). Both depend on the assimilate availability, the first one for the growing spike and developing florets during the 20 days before anthesis [57][58][59], and the second one during the 10 days after anthesis [57,60]. This would explain the high importance in GN and FF determination of: (i) the spike dry weight achieved at anthesis (SDW) [10,12,61,62]; and (ii) the dry matter partitioned within the spike between florets/grains and spike structure parts, i.e. the fertile floret efficiency (FFE, fertile florets per g of SDW) [63] and the fruiting efficiency (FE, grains per g of SDW, or FEm grains per g of chaff at maturity) [15,[63][64][65][66][67][68]. It has been reported that in modern elite cultivars the SDW was less important to explain GN variation than the fruiting efficiency [15,[63][64][65][66][67][68][69]. The GST is considered to be high in relative modern cultivars (i.e., > 80% of fertile florets set grains) [10,67,70], but it has been recently shown that it can be as low as 60% [45,63]. Then, the amount of assimilates partitioned within the spike to grains, to chaff (CH) or to its structures (glume, lemma, palea, awns GLPA-, and rachis R-), as well as the GST and SDW, are worthy of study. Only few reports looked for the genetic bases of these traits [26,45,48,71,72] (Table  S1).
In a previous work [72], dealing with one of the DH populations used in the present study, we identified and validated the novel QFEm.perg-3A for FEm on chromosome 3A, and the first known QTL for FFE and FE, QFFE.perg-5A, located on chromosome 5A. This last QTL was also detected when the FEm was measured, agreeing with Basile et al. [71] who detected two regions within this QTL associated with FEm. Despite studying the pleiotropic effect of both QFFE.perg-5A and QFEm. perg-3A on the associated traits of spike fertility previously mentioned, we did not look for new major and stable QTL for those traits themselves.
The aim of this study is to identify stable and major QTL for the spike fertility and related traits (numerical and physiological components) and to discuss the possible pleiotropic effects among them and with the previously reported QFEm.perg-3A and QFFE.perg-5A. Two doubled haploid populations were used, Baguette Premium 11 × BioINTA 2002 (BP11xB2002) and Baguette 19 × BioINTA2002 (B19xB2002) derived from the crosses of elite cultivars adapted to the central region of the wheat-producing area of Argentina.
In the present study, we report 305 QTL for different spike fertility and related traits distributed throughout the wheat genome. Nevertheless, only 28 QTL are considered major and stable. Eight genomic regions that group some of the significant and stable QTL have been identified and their pleiotropic effects on other related traits have been also analysed. Finally, we have found two regions (R5A.1 and R5A.2) and a QTL (QGW.perg-6B) that resulted in a final higher spike yield.

Genetic linkage map construction
The linkage map of BP11xB2002 consisted of 7,323 SNPs and two functional markers for the vernalization genes Vrn-A1 [73] and Vrn-B1 [74] (Tables S2, S3). All the SNPs represented 723 loci across the 21 wheat chromosomes. The map covered 2605.3 cM in length with an average locus spacing of 4.7 cM (Tables S2, S3). The linkage map of B19xB2002 was previously described in Pretini et al. [72]. Briefly, the B19xB2002 map consisted of 10,936 SNPs and the Vrn-A1 and Vrn-B1 markers. All the SNPs represented 739 loci across the 21 wheat chromosomes (Tables S4, S5). The map covered 2,221.7 cM in length with an average locus spacing of 4.3 cM (Tables S4, S5). Although the genome length of each population was similar, distribution of the markers in the three genomes was not uniform. In BP11xB2002, genomes A and B were similar with at least three times the number of polymorphic markers than in genome D, with 2,955, 3,513 and  [18] 857 markers, respectively (Table S3). In B19xB2002, the marker uneven distribution was higher with 4,126, 5,448 and 1,364 markers in genomes A, B and D, respectively (Table S5). The number of loci in genomes A and B were three times the number of loci in genome D in both populations. For BP11xB2002 there were 311 and 300 loci in genomes A and B and 113 loci in genome D whereas for B19xB2002 there were 324 and 317 loci in genomes A and B and 98 loci in genome D (Tables S3, S5).

Phenotypic analysis
The means and ranges of the 14 studied traits across the five environments (E1 to E5) for the three parents and both DH populations are detailed in Table S6. According to BLUE values, B19 and BP11 parents had higher FF, FFTS (fertile florets per total spikelet), FFFS (fertile florets per fertile spikelet) and GN, whereas B2002 had higher SDW, SL, TS, CH, R, GLPA and GW (Table 2). BP11 showed the highest and B19 the lowest FS value, while B2002 was in-between ( Table 2). All traits presented a normal distribution across each environment and BLUE values, with a transgressive segregation from both parent lines in both populations ( Table 2, Table S6). The narrow-sense heritability values ranged from 0.31 to 0.86, depending on the trait and the DH population ( Table 2). The phenotypic performances of the GN, FF, SDW, CH and GST in both populations were already described in Pretini et al. [63]. Briefly, the range of mean values for the BP11xB2002 and B19xB2002 populations based on the BLUE values, were: (i) 29.4 to 53.4 and 27.2 to 50.0 grains per spike for GN; (ii) 37.1 to 65.5 and 34.7 to 52.9 florets per spike for FF; (iii) 305 to 502 and 279 to 480 mg per spike for SDW; (iv) 307 to 535 and 323 to 564 mg per spike for CH; and finally (v) 0.6 to 1.1 and from 0.5 to 1.3 for GST (Table 2).
In relation to the traits determining spike structure at anthesis, the highest SL, considering the agronomic environments (E1 to E3) and both populations, was observed in E1 (~ 107 mm), while the lowest SL was observed in E2 (~ 93) ( Table S6). The E5 (a non-agronomic summer season) showed even a lower SL than E2 (~ 73 mm)  (Table S6). The TS ranged similarly for both DH populations within the agronomic environments (from ~ 21 to 24 spikelets per spike in E2 to E1); while the lowest TS was observed in the non-agronomic environment (E5, ~ 16 spikelets per spike, Table S6). The FS ranges were also similar for both populations (from ~ 17 to 20 fertile spikelets per spike in E2 to E1) and the lowest FS was observed in the E5 (~ 11 fertile spikelets per spike, Table S6). The CN was the lowest for both populations in the E2 (4.4 mm per spikelet, Table S6); while the highest (4.6 mm) was detected in E1 for BP11xB2002 and in E5 for B19xB2002 (Table S6). For both populations the FFTS ranged from ~ 2.0 to 2.4 fertile florets per total spikelet (E2 to E1, Table S6), while the lowest FFTS was detected in the E5 with 1.7 fertile florets per total spikelet (Table S6). Even though the range explored by the FFFS was similar for both populations, ~ 2.4 to 2.7 fertile florets per fertile spikelet, the maximum and minimum values were associated with different environments in each population (E5 to E3 in BP11xB2002 and E2 to E1 in B19xB2002, Table S6).
The CH partitioning between R and GLPA varied from 14 to 22% for the R, and from 78 to 86% for the GLPA, depending on DH population and environment. Similarly to the CH, within the agronomic environments the highest R was detected in the E3 (~ 98 mg per spike) and the lowest in the E2 (~ 60 mg per spike) for both populations (Table S6). The R measured in the E5 was even lower than the one of E2 (48 mg per spike, Table S6). The GLPA varied from ~ 216 to 564 mg per spike within the agronomic environments, reaching 234 mg per spike in the E5 (Table S6). The GW ranged from ~ 32 to 47 mg (E4-E3) for BP11xB2002 and from ~ 29 to 41 mg (E2-E5) for B19xB2002 including the non-agronomic environment (Table S6).
For both populations, when the spikes were longer (higher SL), they bore more total and fertile spikelets per spike (TS and FS) distributed in a laxly way (higher mm per node or higher CN) (Figures S1 and S2); this resulted in increased fertile florets per spike (FF). In both populations, the longer spikes were heavier (higher SDW and CH), which showed positive correlation with the FF but reduced efficiency to set fertile florets or grains per unit of spike growth (negative correlations SDW vs FFE and SDW vs FE in Figures S1 and S2). On the other hand, theses efficiencies to set florets and grains were positively correlated with the ability of a spikelet to bear florets (positive correlation and FFE vs FFFS), increasing the fertile florets and grains per spike (positive correlation FFE vs FF and FE vs GN) (Figures S1 and S2). Although the higher GN was positively correlated with the spike yield (YLD), it was negatively associated with the GW in both populations. Meanwhile, the GW contributed to YLD only in one population (BP11xB2002) not showing Table 2 Means, ranges, heritability and Shapiro-Wilk test for all traits based on the BLUE values a SL spike length (mm), TS total spikelets per spike (n° spike −1 ), CN compactness of the spike (mm node −1 ), FF fertile florets per spike (n° spike −1 ), FS fertile spikelets per spike (n° spike −1 ), FFTS fertile florets per total spikelet (n° spikelet −1 ), FFFS fertile florets per fertile spikelet (n° spikelet −1 ), SDW spike dry weight at anthesis (mg spike −1 ), R rachis (mg spike −1 ), GLPA: glume + lemma + palea + awns (mg spike −1 ), CH: chaff (no-grain spike dry weight at harvest, mg spike −1 ), GN: grain number per spike (n° spike −1 ), GW: grain weight (mg), GST: grain set b SD standard deviation c W: A modification of the test of Shapiro-Wilks for normality. Mahibbur S2). It is interesting that the higher efficiencies to set fertile florets (FFE) and grains (FE) resulted in an increased efficiency to set yield per unit of spike growth at anthesis in both populations (r = 0.48 or 0.19 for FFE vs YLD/ SDW p < 0.05; and r = 0.81 or 0.68 for FE vs YLD/SDW p < 0.0001, in B19xB2002 and BP11x B2002, respectively). These efficiencies are fundamental considering the limitation of assimilates for spike growth during the preanthesis period (for a thorough discussion see Pretini et al. [72]).

QTL mapping analysis
A total of 305 QTL were identified across 5 environments and BLUE distributed on the 21 chromosomes (Table S7). Nevertheless, only 28 QTL were stable, i.e. present in at least 3 individual environments or BLUE analysis with a LOD > 2.5 considering a single population or a combination of both populations but with the contribution of the same germplasm (Baguette or B2002), and major, i.e. the R 2 > 10% in one environment at least (Table 3). Those stable and major QTL were distributed on the 1A, 2A, 2B, 2D, 3A, 3B, 5A, 5B, 6A, 6B, 7A and 7B chromosomes ( Table 3). For SL 19 QTL were detected (Table S7); however, three of them on chromosomes 2B (QSL.perg-2B), 5A (QSL.perg-5A) and 7A (QSL.perg-7A) were considered major and stable across environments. Both QSL.perg-2B and QSL.perg-7A were detected on BP11xB2002 while the QSL.perg-5A was present in B19xB2002 ( Table 3). The increasing allele was always contributed by B2002 with an additive effect that ranged from 3.1 to 3.8 mm for QSL.perg-2B, from 2.0 to 4.0 mm for QSL.perg-5A and from 2.2 to 4.3 mm for QSL.perg-7A (Table 3). Significant epistatic interaction between QSL.perg-2B and QSL.perg-7A was detected (P = 0.042). The QSL.perg-7A allele from BP11 produced a greater reduction in SL in the presence of the allele from the same parent for QSL. perg-2B ( Figure S3a).
For CN, two QTL out of 22 (Table S7), on chromosomes 2A (QCN.perg-2A) and 5A (QCN.perg-5A) were considered major and stable across environments. The QTL were identified in both populations. The QCN.perg-2A was detected in one environment in BP11xB2002 (E3) and in three environments in B19xB2002 (E1, E5 and BLUE). The QCN.perg-5A was detected in two environments in BP11xB2002 (E2 and BLUE) and in one environment in B19xB2002 (E3, Table 3). The increasing allele of both QTL was contributed by B2002 and had an additive effect ranging from 0.13 to 0.20 mm per node (Table 3).
From the 21 QTL identified for FFTS (Table S7), two QTL on chromosomes 5A (QFFTS.perg-5A) and 5B (QFFTS.perg-5B) were considered major and stable. The QFFTS.perg-5A was detected in one environment in BP11xB2002 (E3) and in three environments in B19xB2002 (E1, E3 and BLUE) while the QFFTS.perg-5B was detected in four environments in BP11xB2002 (E1, E2, E3 and BLUE) ( Table 3). For both QTL, the increasing allele was contributed by the Baguette parents with an additive effect that ranged from 0.09 to 0.13 and 0.09 to 0.16 fertile florets per total spikelet per spike for QFFTS.perg-5A and QFFTS.perg-5B, respectively (Table 3).  Only one QTL for FFFS on chromosome 7B (QFFFS. perg-7B) was major and stable from the 22 QTL identified (Table S7). It was detected in B19xB2002 (Table 3), in two environments (E1 and E3) and the BLUE ( Table 3). The increasing allele was contributed by B19 with an additive effect that ranged from 0.09 to 0.15 fertile florets per fertile spikelet (Table 3).
Three QTL for R on chromosomes 2A (QR.perg-2A), 3B (QR.perg-3B), and 6A (QR.perg-6A) were major and stable; out of the 31 QTL identified for R (Table S7). The QR. perg-3B was only detected in BP11xB2002 while the QR. perg-2A and QR.perg-6A were identified in both populations. The QR.perg-2A was present in two environments in BP11xB2002 (E1, BLUE) and in one environment in B19xB2002 (E5), whereas the QR.perg-6A was detected in two environments for each DH population ( Table 3). The increasing allele of QR.perg-2A and QR.perg-3B was contributed by the Baguette parents with an additive effect that ranged from 3.3 to 8.0 and 4.4 to 9.5 mg per spike, respectively. In contrast, the increasing allele of QR.perg-6A was contributed by B2002 with an additive effect varying from 3.0 to 8.4 mg per spike (Table 3).
Despite being detected 23 QTL for GLPA, ( Table  S7) Table 3). The increasing allele for QGLPA.perg-1A, QGLPA.perg-3A and QGLPA.perg-7A was contributed by B2002, with an additive effect ranging from 27 to 33, 11 to 24 and 13 to 33 mg per spike, respectively (Table 3). On the other hand, the increasing allele for QGLPA.perg-5B was contributed by the Baguette parents with an additive effect that ranged from 18 to 41 mg per spike (Table 3). Significant epistatic interaction between QGLPA.perg-1A and QGLPA.perg-3A was detected (P = 0.021). The GLPA.perg-1A allele from B19 produced a greater reduction in GLPA in the presence of the allele from both parents (B19 and B2002) for GLPA.perg-3A, but the effect was grater in the presence of the B19 ( Figure S3b).
(See figure on next page.) Fig. 1 Genomic regions represented on chromosomes with markers and their reference position. *The corresponding physical distances (Mb) of the QTL regions were obtained by blasting the flanking SNP markers (± 1 LOD) of the most separated QTL in the region to the Chinese Spring RefSeq v1.0 sequence ** The "a" indicates that BP11 or B19 allele increases the corresponding trait and the "b" indicates that B2002 allele increases the corresponding trait. Red letters are for the BP11xB2002 population and blue letters are for the B19xB2002 population *** Solid and dot lines indicates the markers for each hotspot **** In GN BLUE the-1 LOD SNP (Table S9) was located in a map space without markers, the closest one is 341.3 cM apart, for this reason this environment was ruled out to determine the interval of the R5A.1   For GN only one QTL on chromosome 5A (QGN. perg-5A) was major and stable out of the 18 QTL identified (Table S7). The QTL was detected in B19xB2002 (Table 3), in two environments (E2 and E5) and in the BLUE ( Table 3). The increasing allele was contributed by B19 with an additive effect ranging from 1.8 to 2.3 grains per spike (Table 3).
From the 21 QTL detected in the analysis for GW, only two QTL on chromosomes 5A (QGW.perg-5A) and 6B (QGW.perg-6B) were major and stable. The QGW.perg-6B was present in B19xB2002 while the QGW.perg-5A was detected in two environments in BP11xB2002 (E3 and BLUE) and in one environment in B19xB2002 (Table 3). Baguette parents with an additive effect that ranged from 1.2 to 3.0 mg contributed to the increasing allele for QGW.perg-5A (Table 3). In contrast, for QGW.perg-6B, the increasing allele was contributed by B2002 ranging the additive effect from 1.6 to 1.8 mg (Table 3).
No QTL was considered major and stable for SDW and GST despite the 24 and 14 QTL respectively identified in the analysis (Table S7).

Discussion
Most of the breeding progress in wheat yield potential has been achieved by selection of yield per se due to the lack of reliable secondary traits and molecular information available to be used in marker assisted selection (MAS) [7]. The yield potential improvement was, in most cases, consequence of increased GN [8][9][10][11][12][13][14][15], though effects of GW have been reported recently [16][17][18]. In the present paper we identified one stable and major QTL for GN on chromosome 5A (QGN.perg-5A) mapping in the same position that the one reported by Guan et al. [28]. Nevertheless, we recently reported this position as primarily controlling the fertile floret efficiency (FFE, fertile florets per g SDW) when the QFFE.perg-5A was identified and validated [72]. Then, as the GN is the result of the FFE and GST (both defining FE) together with the SDW [76,77], the QFFE.perg-5A can be detected as a QTL associated with GN, highlighting the relevance of FFE and the QTL validated to define GN. This result exemplifies the importance of dissecting the traits into simpler and more heritable components because it enables a better search for the actual candidate gen in further research.
In relation to GW, we detected two QTL, one on chromosome 5A and other on chromosome 6B ( Table 3). The first one has already been reported [24,26,47], but the second one, the QGW.perg-6B is novel. This QTL is located 157.7 Mb apart from the homeologous gene of GW2 in B genome (GW2-B1), associated with grain size [78], suggesting that it would not be a candidate gene to explain the phenotypic variations observed.
The GN is a complex trait itself, being the result of many numerical and physiological spike fertility related traits. In the present study, 25 major and stable QTL for spike fertility and related traits were detected (without considering the one for GN and the two of GW already mentioned in the previous paragraph). There were only two traits, SDW and GST, for which no stable and major QTL were detected. This agrees with the low narrowsense heritability observed (see Table 2) and highlights the high impact of environment on those traits (see Table 3 in Pretini et al. [63]). Considering the three QTL detected for SL, the QSL.perg-2B is 13.4 Mb apart from the one already described by Cui et al. [41] (Table S1), and the QSL.perg-5A is located in the same region as a previously reported QTL [26,29,32,35,36] (Table 1).
In contrast, no equivalent regions have been detected in previous studies for the QSL.perg-7A. Regarding the three QTL identified for TS, the QTS.perg-2D partially overlaps with one previously reported [54] (Table S1). Meanwhile, the QTS.perg-7A is in the same region as a previously identified QTL [27,31,32,36,38,39,41], (Table 1), and co-localizes with the recently described WAPO-A1 gene (674.07 Mb) that modifies the total number of spikelet per spike [79]. Finally, the QTS.perg-3A detected in our work has not been previously described. In relation to the CN, the QCN.perg-5A is in the same region as a previous detected QTL [32,36], but the QCN.perg-2A is a novel one.
Interestingly, the two QTL detected for FF are novel. The QFF.perg-2B is approx. 539 Mb apart from the QTL for FF detected by Guo et al. [45] discarding that it is the same region, while for QFF.perg-7B, no equivalent regions have been reported previously. For the FS, only the QFS. perg-2B is 540 Mb apart from another QTL detected previously [27,34] (Table 1). The remaining two QTL, QFS. perg-3A and QFS.perg-5B, do not share their regions with other previous works.
For the rest of the traits analysed in this study (FFTS, FFFS, R, GLPA and CH), no previous reports are available to the best of our knowledge (Table S1) [80], known to increase the number of grains through higher fertile florets per spikelet, has been identified.
As many of the spike fertility traits detected in this study had similar positions, we identified eight genomic regions that share 17 major and stable QTL for the different traits (R1A, R2B, R3A, R5A.1, R5A.2, R5B, R7A and R7B). Only in two of these regions (R5A.1 and R5A.2) other QTL for the same trait have been previously described. The remaining six regions are identified for the first time as important hot spots for spike fertility traits (Fig. 1). Interestingly, the R5A.1 region, which contains QTL for SL, CN and GN, is close to the QFFE.perg-5A identified and validated for fertile floret efficiency in Pretini et al. [72]. The allele of B2002 parent increases the SL and CN while the allele from B19 parent increases the GN via QFFE.perg-5A. These results agree with the performance of the parental lines described in the present study ( Table 2).
The region R5A.2, which includes QFFTS.perg-5A and QGW.perg-5A, coincides with the location of the vernalization response gene Vrn-A1; while the R5B region, which includes QFS.perg-5B and QGLPA.perg-5B, coincides with the location of the other vernalization response gene Vrn-B1. The three parental lines of the two DH populations used in the present study are spring wheats (Vrn-A1b/vrn-B1 /vrn-D1 for B19 and BP11 and vrn-A1 /Vrn-B1 /vrn-D1 for B2002); and mostly insensitive to photoperiod (Ppd-D1a allele). This agrees with the very close anthesis dates of the lines within each population described in Pretini et al. [72], except for the summer sowing (E5) of B19xB2002 population, in which the range was higher. Furthermore, to test the effect of the two genes, we made an ANOVA using the functional markers for Vrn-A1 and Vrn-B1 as source of variation for time to anthesis and small differences were detected between the anthesis dates for both populations. In BP11xB2002, only five to seven days difference in time to anthesis was associated to the allelic constitution of Vrn-A1 or Vrn-B1, respectively (Table S10). Similarly, for B19xB2002, three days difference in anthesis date was detected depending on the allelic constitution for both genes (Table S10). No epistatic interaction was observed between Vrn-A1 and Vrn-B1 in BP11xB2002, while a difference of up to 7 days to anthesis was observed depending on the allelic constitution in the B19xB2002. Based not only on those results but also on the fact that most of the QTL included in the R5A.2 and R5B regions were not expressed in the summer sowing of the environment E5 (except for QGLPA. perg-5B in B19xB2002), these QTL are not considered to be masking an important phenology effect. In contrast, it could be indicating that the Vrn-A1 and Vrn-B1 allelic variation in the population might have a pleiotropic effect on the spike traits located in those regions, with little impact on phenology in the tested conditions.
The spike fertility and related traits are correlated, positively or negatively, depending on the trait (Figures S1 and S2) [81]. In addition, a negative correlation is usually observed between GN and GW [11,16,82], which was also present in our data set. Then, we enquired about the possible pleiotropic effects of each of the eight regions detected over the other spike related traits, GW and final yield per spike (YLD), following the Fig. 2. For this reason, we performed an ANOVA for each of the evaluated traits using the QTL peak marker as fixed and the environments as random class variables in the model. Four regions had a significant effect on GN (R2B, R3A, R5A.1 and R5A.2), six on GW (R1A, R2B, R5A.1, R5A.2, R7A.1, R7A.2, and R7B), but only two in spike YLD (R5A.1, R5A.2). For the R5A.1 region (QSL.perg-5A, QCN.perg-5A and QGN.perg-5A) when the QTL from B19 is present, it results in a shorter spike (-6% SL) with similar TS (-2%) or FS (ns), due to a reduction in the distance between spikelets (-5% CN). The FF increases 4% due to higher FFE (+ 10%), despite a reduction in the SDW (-3%) which is accompanied by a 3% increment of the FFFS. The FF increment together with the higher GST (+ 8%) results in an increment in the GN (+ 7%), which translates into higher yield (+ 3%) despite a significant reduction in GW (Fig. 2). As we previously mentioned, this region includes the QFFE.perg-5A identified and validated for fertile floret efficiency in Pretini et al. [72], also within the B19xB2002 population and showed similar pleiotropic effects to the R5A.1 region. The other region that resulted in a final higher YLD was the R5A.2, which contained the QGW.perg-5A and QFFTS.perg-5A. When this region from B19 is present, the SL is not affected, but the distance between spikelets is increased (+ 3% CN), reducing the TS (-2%). The FFTS and the FFFS increase 3 and 2%, respectively and the FFE is higher (+ 6%). Nevertheless, the GN is not significantly improved. The YLD improvement of R5A.2 (+ 5%), when B19 alleles are present, is consequence of the increased GW (+ 6%) (Fig. 2). As far as we know, the pleiotropic effect of these regions had not been previously reported, except for Pretini et al. [72] for the QFFE.perg-5A, which is within the R5A.1. We made a similar analysis of pleiotropic effect for each QTL identified (data not shown), being the QGW.perg-6B the only one that has a pleiotropic effect in YLD. When the B2002 alleles are present, the spikes are longer (+ 2% SL), but the TS and CN are not significantly modified. Nevertheless, higher FS were detected (+ 2%), which was counterbalanced by a reduction in the FFFS (-2%) resulting in no impact on GN. The YLD increment (+ 5%) was consequence of the increased GW (+ 10%). This is an interesting result highlighting the relevance of this QTL for the first time.

Conclusion
From the 14 analysed traits, only two of them did not show major and stables QTL (SDW and GST). For the rest of the 12 traits, there were up to 28 significant and stable QTL and 8 hotspot regions detected. Based on the complex pleiotropic analysis preformed it is concluded that the R5A1 and R5A.2 regions together with the QGW. perg-6B are of high relevance to be used in MAS to improve a set of traits related with spike yield potential. All the QTL identified for the spike related traits are the first step to search for their candidate genes, which will allow their better manipulation in the future.

Plant materials
Two doubled haploid (DH) populations were developed from the crosses between Baguette Premium Cycles to anthesis in optimal sowing dates are similar for the three parents [59]. The GW of B2002 was higher compared to the one of BP11 and B19, whereas the GN of B19, followed by BP11, was higher than that of B2002 [68].

Experiments and phenotyping
The DH populations were grown in two experimental sites: EEA Pergamino (33 • [72]. When plants reached anthesis stage (Z6.1, [83]) five to three median spikes were selected from a larger sample (0.5 m from central row) in E1, E2 and E3, while three main stem spikes were chosen in E5. No samples were taken in E4 at anthesis. The spike length (SL, mm), the number of total spikelets per spike (TS), the number of fertile spikelet (FS), and the number of fertile florets per spike (FF) were measured following the methodology described in Pretini et al. [63]. The ratio between SL and TS was used to determine the spike compactness (CN, mm spikelet node −1 ), whereas the ratios FF/TS and FF/ FS were used to estimate the number of fertile florets per total spikelet per spike (FFTS) and the number of fertile florets per fertile spikelet per spike (FFFS), respectively. After drying in an oven at 70 • C for 76 h the spike dry weight at anthesis (SDW) was estimated.
When plants reached maturity (Z9, [83]) a second spike samples were performed in a similar procedure as the one described for anthesis (including E4 environment).
Before threshing by hand, all the spikes were dried in an oven and weighted. For E1, E3, E4 and E5, the rachis (R) and the rest of the no-grain parts (glume + lemma + palea + awns, GLPA) were separated when threshing and weighted. The chaff (no-grain spike dry weight at maturity) was calculated as the sum of R + GLPA. For E2 the chaff was estimated by subtracting the weight of all the grains from the dry weight of the spike before threshing, because no chaff dissection was performed. The grain number (GN) of each spike was counted in E2, E3, E4 and E5 using an automatic counter. The grains from E1 were discarded because they were severely affected by Fusarium head blight. The grain weight (GW) was estimated as the ratio between the weight of all grains and the GN. The grain set (GST) was estimated by the ratio between GN/FF. All the phenotypic data used in this work for both populations is available in Tables S11 and S12.

Data analyses
For each DH line, the mean value of each trait was calculated across the two replicates for E1 to E4 and the six replicates for E5. The Shapiro-Wilk test and the quantile-quantile (q-q) plot was performed to test for normal distribution. The analysis of variance (ANOVA) was performed using the Infostat/p software [84]. In addition, Best Linear Unbiased Estimator (BLUE) was estimated for each DH line including all tested environments; as random variable using R v3.3.2 and the Pearson's correlations with the BLUE values were made to determine the relationship between all traits. The narrow-sense heritability of the traits was calculated as: where σ 2 G is the genotypic (additive) variance, σ 2 GE is the G × E interaction variance, E is the number of environments, R is the number of replications, and σ 2 RES is the error variance [85].

Linkage map construction and QTL analysis
The DH populations and the three parents were screened with the iSelect 90 K array containing 90,000 wheat SNP markers [86]. Additionally, Vrn-A1 [73] and Vrn-B1 [74] markers were added to the DH genetic map. The SNP markers with a high number of missing/ heterozygous data (> 20%) were discarded for the construction of the linkage map. SNPs with a 1:1 segregation distortion greater than 20% were also eliminated. Then, the dataset was reduced by merging SNPs markers with identical segregation with the Python script, merger.py 1 [72]. Finally, the publicly available R package "Rqtl" [87] was used for the linkage map development. The physical position of SNPs associated with phenotypic traits was established by BLAST against the IWGSC Ref. Seq. v1.0 wheat RES ER genome assembly [22]. Complete linkage maps developed  for both populations are available at the Tables S2 and S4.
The mean value of the trait in each environment and the BLUE values (which were treated like an additional environment) were used in the QTL mapping. The QTL analyses was performed with QTL Cartographer 2.5 software [88] through composite interval mapping (CIM) with the standard model. For the standard model we used a control marker number of 5, a window size of 10 cM and a forward and backward regression method with 500 permutations at α = 0.05. A LOD value of 2.5 was selected as a uniform threshold for all analyses. Detected QTL for a given trait with overlapping support intervals (< 50 Mb) were considered as equivalents. The QTL were considered "stable" if they were detected in a minimum of three environments and were defined as "major stable" if they present a R 2 > 10% in one environment at least. For all evaluated traits in each individual DH population, we performed a factorial ANOVA using the peak marker for each major and stable QTL as class variables in the model, along with all possible two-way interactions in the case that more than one QTL was detected in order to determine significant epistatic interactions.