Mapping dynamic QTL dissects the genetic architecture of grain size and grain filling rate at different grain-filling stages in barley

Grain filling is an important growth process in formation of yield and quality for barley final yield determination. To explore the grain development behavior during grain filling period in barley, a high-density genetic map with 1962 markers deriving from a doubled haploid (DH) population of 122 lines was used to identify dynamic quantitative trait locus (QTL) for grain filling rate (GFR) and five grain size traits: grain area (GA), grain perimeter (GP), grain length (GL), grain width (GW) and grain diameter (GD). Unconditional QTL mapping is to detect the cumulative effect of genetic factors on a phenotype from development to a certain stage. Conditional QTL mapping is to detect a net effect of genetic factors on the phenotype at adjacent time intervals. Using unconditional, conditional and covariate QTL mapping methods, we successfully detected 34 major consensus QTLs. Moreover, certain candidate genes related to grain size, plant height, yield, and starch synthesis were identified in six QTL clusters, and individual gene was specifically expressed in different grain filling stages. These findings provide useful information for understanding the genetic basis of the grain filling dynamic process and will be useful for molecular marker-assisted selection in barley breeding.

NC). Correlation analysis between the traits was performed using the "PROC CORR" program. ANOVA of each component was performed using PROC GLM procedure. The broad-sense heritability was calculated using the formula: h B 2 = σ g 2 /(σ g 2 + σ ge 2 /n + σ e 2 /rn), where σ g 2 was the genetic variance, σ ge 2 was the genetic-by environment interaction variance, σ e 2 was the error variance, and r and n were the number of repetitions for each genotype and environment, respectively.
Dynamic QTL analysis. The genetic map consisting of the 1962 markers cover a total of 1375.8 cM genomic regions was used to screen QTLs 43 . Unconditional phenotypic values were measured at different stages of I, II, III, IV, V, VI and VII after flowering time. Conditional phenotypic values were the incremental grain filling-related traits values in adjacent stages ΔT2 (II-I), ΔT3 (III-II), ΔT4 (IV-III), ΔT5 (V-IV), ΔT6 (VI-V) and ΔT7 (VII-VI), since the grain size traits were recorded from stage II, so the conditional phenotypic values were calculated from ΔT3, and the data collection method was used according to Zhu 33 . The conditional phenotypic value ΔT1 of GFR from the flowering time to the stage I was the phenotypic value of stage I, so the effects were considered to be the unconditional genetic effects of stage I. QTL IciMapping 4.1 with inclusive composite interval mapping (ICIM) model was used to analyze the locus and effects of unconditional QTLs 44,45 . Conditional variable analysis method 33 combined with ICIM was used to perform conditional QTL analysis of GA, GP, GL, GW, GD and GFR at each stage. Furthermore, in order to eliminate the interference of row type (Rt) and caryopsis type (Ct), we used QTL.gCIMapping.GUI software 46 for covariate QTL analysis using Rt and Ct as covariates. The scan step and probability in stepwise regression (PIN) were set to 1 cM and 0.001, respectively. The logarithm of the odds (LOD) threshold was set to 3.0 after 1000 permutations on a 0.05 Type 1 error, so the QTL was declared based on the LOD threshold of 3.0, and these QTLs were termed 'identified QTLs' . QtL integration. Goffinet and Gerber 47 first used the QTL meta-analysis method to integrate QTLs from independent experiments and determine the corresponding confidence intervals for consensus QTLs, this

Results
Phenotypic variation. The phenotypic data of the parents and DH line were listed in Table 1. The value of each trait in Huadamai6 was significantly higher than Huai11 (P < 0.01) except for GW at stage II. For the two parents, GFR increased slowly during the initial two stages, after which GFR increased rapidly and reached its maximum at stage IV, and then decreased continuously in two years (Fig. 1B). In the DH population, the tendency of GFR in the seven sampling stages was basically the same as that of the parents (Fig. 1A). The frequency distribution of the five grain size traits was shown in Supplementary Fig. S2. Each trait was continuous variation in all sampling stages, and the absolute values of skewness and kurtosis of most traits were less than 1 except for GFR at stage VII, indicating that these traits followed a normal distribution and were suitable for QTL analysis (Table 1). ANOVA showed statistically significant effects for genotype, environment and G × E interaction with all traits (P < 0.01) (Supplementary Table S1). The broad-sense heritability (h B 2 ) was estimated between 58.9% and 98.5% (Supplementary Table S1).
The Pearson correlation coefficients (r) between GA, GP, GL, GW, GD, TGW, GFR mean , GFR max and GFR traits at seven sampling stages over two years were given in Fig. 2. Grain size traits were significantly correlated with each other at each stage (not listed), and grain size traits at all stages were significantly positively correlated with GFR at stage II, III, IV and V, GFR mean , GFR max and TGW, whereas GFR at stages I, VI and VII were mostly unrelated or negatively correlated with grain size traits at each stage. TGW showed high positive correlation coefficients with GFR mean and GFR max (r > 0. 9), while GFR at stage III, IV and V were also positive correlation coefficients with TGW (r > 0.7). The initial and final stages of GFR was not correlated with GFR mean and TGW, but negatively with GFR max (Fig. 2).

Unconditional QTL analysis.
A total of 204 unconditionally identified QTLs for GFR mean , GFR max , GA, GP, GL, GW, GD and GFR were detected on 7 chromosomes, individually accounted for 1.12-78.03% phenotypic variation (Supplementary Table S2). Up to 152 identified QTLs with overlapping CIs were detected and integrated into 38 consensus QTLs. The other 52 non-overlapping identified QTLs were also considered to be consensus QTLs, of which 8 QTLs for GFR mean and GFR max were regarded as the unconditional consensus QTLs (Table 2,  Supplementary Table S2). In total, 90 unconditional consensus QTLs were detected (Fig. 3, Table 2).
For the GFR max and GFR mean , eight QTLs were detected on chromosome 1H, 2H, and 7H over two years, individually explaining 2.44-74.83% of PVE (Table 2). Among these QTLs, two main-effect QTLs were detected on 2H, ucqGFR mean 2 was located at 125 cM close to the marker 2HL_22930294, and ucqGFR max 2-2 was located at 127 cM near the marker 2HL_22930005. These alleles from Huadamai6 contributed PVE 68.26-74.83%.
A total of 30 identified QTLs on 5 chromosomes (excluding 4H and 6H) were detected for GL. Of these QTLs, integrating 24 QTLs with overlapping CIs into 6 consensus QTLs and 6 non-overlapping QTLs were also considered to be consensus QTLs ( Table 2, Supplementary Table S2). At stages III-VII, the major QTL ucqGL2-1 was 125.94 cM from 2_527241334 and explained up to 15.52-51.03% of PVE. At stages IV-VII in Y2, the main-effect ucqGL7-1 was closely linked to the marker 7HL_13143105 of 65.25 cM, accounting for 25.01-45.19% PVE. ucqGL7-3 tightly linked to the 7HS_29196961 was detected at stage VI and VII in Y2, explaining 23.52-27.20% of PVE. ucqGL7-5 closely linked to the marker Bmac31 of 150.92 cM was expressed at stages V-VII in Y1 and at stage II in Y2, explaining 17.10-43.21% of PVE. In addition, two QTLs ucqGL2-2 and ucqGL3-1 were detected at one or two stages in both years.
For the GW, 33 identified QTLs distributed on 7 chromosomes except 5H were detected over two years. Of which, 28 QTLs with overlapping CIs and 5 non-overlapping QTLs were integrated into 10 consensus QTLs ( Table 2, Supplementary Table S2). ucqGW1-1, and ucqGW3-2 were expressed at one or two stages across two years. ucqGW4 was detected at stage IV and VII in Y2 and exhibited negative additive effects. In addition, ucqGW7-2 was 65.13 cM from 7HL_13143105 and expressed at stage IV, V and VII in Y1 and at stages IV-VI in Y2, explaining 7.63-46.68 of PVE. ucqGW3-2 and ucqGW7-2 expressed additive effects in opposite directions at different stages. Only ucqGW2 closely linked to the 2HL_22930294 of 124.69 cM was identified steadily at all stages in both years, explaining up to 32.27-75.51% of PVE.
For the GD, 39 identified QTLs located on 7 chromosomes were detected across two years. Of which, 26 overlapping were integrated into 6 consensus QTLs, the other 13 non-overlapping QTLs were also considered as consensus QTLs (Table 2, Supplementary Table S2). ucqGD1-1, ucqGD3-3 and ucqGD5-2 were detected at one or two stages over two years with −0.14 to −0.05 additive effects. The major QTL ucqGD2 closely linked to the 2_527241334 of 125.93 cM was stably expressed at all stages in both years, accounting for 18.90-78.03% PVE. In addition, ucqGD7-1 tightly linked to the 7HL_37199773 of 65.88 cM was expressed at stage VI in Y1 and at stages IV-VII in Y2, accounting for 10.08-31.65% PVE.
For the GFR, 36 identified QTLs distributed on 7 chromosomes except 6H and 4H were detected in both years. The 36 identified QTLs were integrated into 12 consensus QTLs, including 4 non-overlapping identified QTLs (Table 3, Supplementary Table S3). The main-effect QTL cqGFR2-2 was stably expressed at ΔT2, ΔT3, ΔT4, ΔT6 and ΔT7 in Y1 and at all intervals in Y2, explaining 14.92-65.72% of PVE and the additive effects had opposite directions in different intervals. Another major QTL cqGFR3-2 was detected at ΔT3 and ΔT5 in Y1, and ΔT2, ΔT3, ΔT5 and ΔT7 in Y2, accounting for 1.40-22.86% PVE with additive effects ranging from −0.07 to 0.16. Similarly, cqGFR2-1 and cqGFR3-1 had opposite additive effects at different time intervals.
Sixteen identified QTLs, including three QTLs each at ΔT3, ΔT4 and ΔT6, six at ΔT5 and one at ΔT7 were detected for the GA. Nine overlapping and 7 non-overlapping QTLs were integrated into 11 conditional consensus QTLs (Table 3, Supplementary Table S3). cqGA3-1, cqGA3-2 and cqGA7-4 were repeatedly detected at a specific time interval in both years, and the alleles increasing GA were contributed by Huadamai6. Another QTL cqGA2-2 closely linked to the 2HL_17075593 of 127.54 cM was consistently expressed at ΔT4 in Y1 and at ΔT3 and ΔT4 in Y2, explaining 27.71-40.97% of PVE. The remaining 7 conditional consensus QTLs were only detected at a single interval in a particular year. www.nature.com/scientificreports www.nature.com/scientificreports/ For the GP, five identified QTLs consisting of three QTLs at ΔT3 and one each at ΔT6 and ΔT7 were detected on 7H. The five identified QTLs were integrated into 4 conditional consensus QTLs (Table 3, Supplementary  Table S3). Only cqGP7-2 tightly linked to the 7HL_3360534 of 68.99 cM was expressed simultaneously at two intervals (ΔT6 and ΔT7) in Y2, with explaining 17.05-41.39% of PVE. Among these consensus QTLs, three of them increased GP contributing from the Huadamai6 allele, and one from Huaai11.
For the GW, 16 identified QTLs were detected, including four each at ΔT3, ΔT5 and ΔT6, three at ΔT4 and one at ΔT7. Eight identified QTLs with overlapping CIs and 8 non-overlapping identified QTLs were integrated into 11 conditional consensus QTLs (Table 3, Supplementary Table S3). cqGW2-4 was expressed at ΔT3 and ΔT4 in both years, explaining 13.33-55.68% of PVE. Another major QTL cqGW3-2 closely linked to the  www.nature.com/scientificreports www.nature.com/scientificreports/ 3HL_32354397 of 32.19 cM was repeatedly detected at ΔT3 in Y1 and at ΔT6 in Y2, and explained 16.11-20.59% of PVE. However, the additive effects of cqGW3-2 at different intervals were opposite. In addition, cqGW3-1 was detected at ΔT5 in Y1 and at ΔT3 in Y2, with opposite additive effects. The remaining 8 consensus QTLs were detected at a single interval in particular year, 6 of them increased the GW were from Huaai11 alleles.
For the GD, 3, 2, 2 and 4 identified QTLs were detected at ΔT3, ΔT4, ΔT5 and ΔT7, respectively. The 11 identified QTLs were integrated into 9 conditional consensus QTLs, the 7 non-overlapping identified QTLs were also regarded as consensus QTLs (Table 3, Supplementary Table S3). cqGD2-2 and cqGD7-1 were stably detected at the same interval over two years. The major QTL cqGD2-2 was 129 cM from the 2HL_34260490, accounting for 18.55-20.41% PVE. The other 7 consensus QTLs were only detected at one interval in a particular year. None of the consensus QTLs were expressed at all intervals, and some QTLs showed opposite effects in different intervals.
In conclusion, most of the QTLs detected by the unconditional QTL mapping method were detected in the V, VI and VII stages, while the QTLs identified by the conditional QTL mapping method were expressed in different periods (Fig. 4).
In addition, to reduce the interference of Rt and Ct on the grain shape and grain filling rate, we used the genome-wide composite interval mapping (GCIM) to perform covariate QTL analysis with Rt, Ct, Rt + Ct as covariates, and detected 118, 109 and 84 QTLs on all seven chromosomes, respectively (Supplementary Table S4). Among these 311 covariate QTLs, 140 covariate QTLs were consistent with unconditional or conditional consensus QTLs, and the remaining 171 were new QTLs. Comparing the consensus QTLs for grain filling rate and grain size traits detected by the five QTL mapping methods, we identified 34 major consensus QTLs (repeatedly detected in at least two QTL mapping methods), including 2, 4, 7, 9, 5, 5, 1 and 1 consensus QTLs for GFR, GA, GP, GL, GW, GD, GFR max and GFR mean , respectively (Table 4). Of these major consensus QTLs, most of them were identified in two QTL mapping methods, and only six were simultaneously identified in at least four QTL mapping methods.

Integration of unconditional and conditional QTLs. Unconditional and conditional consensus QTLs
were further integrated into unique QTLs (Details of unique QTLs were shown in Supplementary Table S5), and listed in Table 5. For the GFR, unique QTL (uqGFR2-3) was repeatedly detected at three different stages and multiple time intervals over two years. For the GA, unique QTL uqGA2-2 was stably detected at five different periods and multiple intervals across two years. For the GL, uqGL7-1 was steadily expressed in several different periods and time intervals. For the GW, unique QTL uqGW2-4 was repeatedly expressed at all sampling stages. For unconditional and conditional QTL integration, some unique QTLs were detected at multiple stages, but not detected at the final stage, these unique QTLs might play an important role in the grain filling stage.
QTL clusters in genome. QTL clusters were defined as a region containing multiple QTLs of various traits within approximately 20 cM 53 . Total of 11 QTL clusters were detected on 1H (two clusters), 2H (one cluster), 3H (three clusters), 5H (one cluster), and 7H (four clusters) ( Table 6). Among these QTL clusters, three QTL clusters (C1, C6 and C8) affected all grain size traits, and the cluster C3 affected all traits. The remaining seven QTL clusters affected at least four traits and three of them were associated with three different grain size traits. The consensus QTLs associated with all traits in C1 and C6 showed negative additive effects, and the alleles that increased GFR and grain size were from Huaai11. Conversely, the alleles of consensus QTLs that increased all grain size traits in C3 were from Huadamai6. In addition, QTLs in these three QTL clusters (C2, C7 and C10) increased GFR and grain size traits alleles were from different parents.

Discussion
Many studies have reported that genetic difference in grain yield was related to difference in GFR in barley 7,9,10 . In addition, GFR mean and GFR max have been reported in maize and wheat as important factors regulating grain weight 42,54 . The GFR mean and GFR max were significantly associated with TGW (r > 0.9) in our experiments (Fig. 2), indicating that GFR mean and GFR max also promoted grain weight in barley.
Grain development is a dynamic process that is regulated by three physiological stages: (1) grain formation period, mainly the division of endosperm cells and the formation of basic structure of seeds, during which there www.nature.com/scientificreports www.nature.com/scientificreports/ is almost no accumulation of dry matter; (2) the linear growth period of dry matter, the accumulation of dry matter in the most vigorous period of the grain, the grain weight during the period almost increased linearly; (3) maturity, the grain weight increased slowly during this period [18][19][20][21]55 . The substance accumulation mainly occurs in the linear growth period of dry matter. In this study, seven sampling from grain formation period to maturity were used to evaluate GFR. The GFR at stage III, IV and V were significantly positively correlated with TGW (r > 0.7), while, there was no correlation between the initial and final stages of GFR and TGW (Fig. 2). This results were basically consistent with the previous studies on the grain weight of linear dry matter accumulation period 20 .
Grain size can be divided into components such as GA, GP, GL, GW and GD 56,57 . In barley, previous studies on grain shape were mainly at maturity, and only major QTLs controlling GL traits were found [57][58][59][60] . However, QTLs detected at maturity may not observe their genetic effects during specific periods of crop development, and dynamic QTL analysis can better understand the developmental behavior of quantitative traits. We found that some of the QTLs for GFR and grain size traits identified on the 2H and 7H chromosomes co-localized with QTLs for yield-related traits, seedling traits and dwarf gene btwd1 detected in previous studies. For example, certain major consensus QTLs tightly linked to 2_527241334 at 126 cM for grain filling rate and grain size detected here, are likely the same to seedling traits QTL qSH2-191 and qSFW2-191 identified by Wang et al. 61 and also likely same to the qSms2-7 and qTgw2-1 for yield-related traits reported by Wang et al. 29 . In addition, a major QTL ucqGL7-5 for grain length closely linked to the Bmac31 of 151 cM identified on 7H chromosome, which may be in the same locus as the dwarfing gene btwd1 reported by Ren et al. 43 . The findings indicated that GFR and grain size traits are closely related to yield and yield-related traits.  In previous studies, QTLs for grain area were detected on chromosomes 1H, 2H, 3H, 4H, 5H and 6H 57,62,63 , and QTLs associated with grain perimeter were identified on chromosomes 1H and 3H 57 . QTLs associated with grain length and width were previously detected on all seven chromosomes 57-60,62-67 , and QTLs for grain diameter were detected on chromosomes 2H, 3H, 4H, 6H, and 7H 67 . Ayoub et al. 57 and Sharma et al. 62 reported a major QTL on chromosome 2H, nearby the locus vrs1 and affected all grain size traits. Most of QTLs for grain sizes identified on 2H in this study were also distributed near the morphological marker vrs1, which are likely same to the QTL reported by Ayoub et al. 57 . Walker et al. 66 detected two QTLs for grain length on 3H, located near markers 2_0662 and 2_1272, respectively. The physical location of the two genetic markers was queried using the Barleymap website (http://floresta.eead.csic.es/barleymap/), we found that the marker 3HL_42780152 in this experiment was close to the physical position of 2_1272, indicating that the QTL cqGL3 adjacent to 3HL_42780152 is likely the same locus reported by Walker et al. 66 . Many new QTLs controlling grain size were detected in our experiment, in which some major QTLs (ucqGA7-1, ucqGP7-2, cqGP7-2, ucqGL7-1 and ucqGD7-1) located near 65 cM of the 7H chromosome were repeatedly detected at different stages or at the same stage of different environments, indicating that these regions might be an important novel locus affecting grain size traits.
Dynamic QTL mapping identified 196 unconditional QTLs and 95 conditional QTLs (Supplementary  Table S2, Table S3). Most QTLs were detected at V, VI and VII stages (Fig. 4). By integrating the unconditional QTL for GFR and grain size traits, some unconditional consensus QTLs were detected simultaneously at several stages, and the majority of the QTLs were identified at stages IV-VII. These results indicated that QTLs associated with grain size traits are selectively expressed during grain filling period. Additionally, we found that some unconditional and conditional consensus QTLs, such as ucqGW3-2, cqGFR2-2, and cqGW3-1, had a combination of identified QTLs with opposite additive effects, indicating that certain QTLs had different expression patterns at different developmental stages or environments. This expression pattern has also been reported for plant height in rapeseed 68 and for grain filling rate in corn 69 .
In this study, 11 clustered QTL regions controlling grain filling rate and grain size traits were found on chromosomes 1H, 2H, 3H, 5H and 7H (Fig. 3, Table 6), and these co-localized QTLs were mainly concentrated on the chromosome 2H, 3H and 7H. Anchoring the SNP markers located in important QTL regions to the Morex genome via the barleymap website (http://floresta.eead.csic.es/barleymap/) identified eight related candidate genes 70 , of which, five associated with plant height and grain size, and three genes were involved in the biosynthesis and metabolism of starch and (1,3;1,4)-β-glucan. Notably, the most important QTL cluster region at the 2HL_16222169 -M_149956_1482 was located on chromosome 2H, containing 18 QTLs controlling all grain size and grain filling rate traits (Fig. 3, Table 6). The vrs1 gene was detected in this region, which was reported in previous studies to control row-number phenotype, affecting grain size and thousand grain weight traits 29,57,61,62 . The Nud gene controlling the hulled/naked caryopsis was detected within QTL cluster affecting GFR, GFR max ,  GA, GP, GL, GW and GD in interval 56.5-70 cM on chromosome 7H, and QTLs affected TGW were located at this locus 71 . Since these two genes have large effects on the population used here, some of the minor-effect QTLs may be affected and difficult to detect. Therefore, to eliminate the interference of these two genes, we performed a covariate QTL analysis to find more new QTLs. Through covariate QTL analysis, we found that 171 (55%) of the new QTLs were undetected by either unconditional or conditional QTL mapping methods (Supplementary  Table S4). The effects of these QTLs on grain filling rate and grain size traits are not as obvious as the vrs1 and Nud genes, but they had an important role underlying these traits. Gibberellin 20-oxidase gene (Hv20ox 2 ) was found in a QTL cluster for GFR, GA, GL, GW and GD between 3HL_32425970 and 3_504106156, located in the interval of 17.5-34.13 cM on chromosome 3H, which was a functional gene regulating barley sdw1/denso, and certain QTLs for yield, grain size and plumpness were co-localized with this gene 72,73 . Within another QTL cluster for GA, GP, GL and GW between 6HL_29119759 and 3 HL_24669457 on chromosome 3H, the flowering time gene HvFT1 was detected, which was the dominant plant transition from vegetative state to reproductive state, affecting the flowering and maturity of barley 74,75 . Furthermore, within the QTL cluster associated with GA, GP, GL and GD in the interval of 150.5-153.5 cM on chromosome 7 H, the novel dwarf gene btwd1 was identified, which not only affect plant height, but affect grain yield at the btwd1 locus 29,43 . The content of starch in barley is 62% to 77% of the grain dry weight, and the grain filling process is mainly the accumulation process of starch. The biosynthesis of starch mainly involves ADP glucose pyrophosphorylase, starch synthase, starch branching enzyme and starch debranching enzyme [76][77][78] . Interestingly, two genes related to starch biosynthesis and metabolism, as well as a gene involved in (1,3;1,4)-β-glucan synthesis, were found in three QTL clusters in this study. The GBSS1b gene, the key gene involved in amylose synthesis, was detected in the C3 region of chromosome 2H, which was close to the 2HL_22930294 marker. The ucqGFR2-3 closely linked to the marker 2HL_22930294 was detected simultaneously at stage II over two years. The GBSS1b transcripts were abundant in the pericarp of flowering and initial grouting 79 , so we considered this to be a candidate gene for the QTL locus. The amy2 gene, which was the most important gene for starch degradation during malting and saccharification 80 , was located in the C8 region of chromosome 7H, close to the 7HL_37199773 marker at 65 cM. Previous studies have reported QTLs underlying malt activity and amylopectin content at the amy2 locus 81 . The HvCslF6, a key gene regulating (1,3;1,4)-β-D-glucans biosynthesis, was found in the C9 region of chromosome 7H 82 . This gene was expressed low in the early stage of grain development and then rapidly up-regulated as the activity of synthetase increased 83 .
In summary, we identified 90, 55 and 311 consensus QTLs using unconditional, conditional and covariate QTL mapping methods, respectively, and detected 34 main-effect QTLs that were simultaneously expressed at multiple stages. The results indicated that these major QTLs were not only expressed at maturity but were also in the early stages of grain development. In addition, eight predicted candidate genes involved in grain yield and starch synthesis pathways were identified in the six clustered QTL regions, which might play an important role in controlling GFR and grain size traits. These findings enhance our understanding of the genetic mechanism of barley grain filling process.