Genome-Wide QTL Mapping for Wheat Processing Quality Parameters in a Gaocheng 8901/Zhoumai 16 Recombinant Inbred Line Population

Dough rheological and starch pasting properties play an important role in determining processing quality in bread wheat (Triticum aestivum L.). In the present study, a recombinant inbred line (RIL) population derived from a Gaocheng 8901/Zhoumai 16 cross grown in three environments was used to identify quantitative trait loci (QTLs) for dough rheological and starch pasting properties evaluated by Mixograph, Rapid Visco-Analyzer (RVA), and Mixolab parameters using the wheat 90 and 660 K single nucleotide polymorphism (SNP) chip assays. A high-density linkage map constructed with 46,961 polymorphic SNP markers from the wheat 90 and 660 K SNP assays spanned a total length of 4121 cM, with an average chromosome length of 196.2 cM and marker density of 0.09 cM/marker; 6596 new SNP markers were anchored to the bread wheat linkage map, with 1046 and 5550 markers from the 90 and 660 K SNP assays, respectively. Composite interval mapping identified 119 additive QTLs on 20 chromosomes except 4D; among them, 15 accounted for more than 10% of the phenotypic variation across two or three environments. Twelve QTLs for Mixograph parameters, 17 for RVA parameters and 55 for Mixolab parameters were new. Eleven QTL clusters were identified. The closely linked SNP markers can be used in marker-assisted wheat breeding in combination with the Kompetitive Allele Specific PCR (KASP) technique for improvement of processing quality in bread wheat.

Dough strength and starch pasting characteristics are quantitative traits controlled by multiple genes. Several QTLs for Mixograph (Huang et al., 2006;McCartney et al., 2006;Patil et al., 2009;Tsilo et al., 2011;Li et al., 2012;Deng et al., 2013;Zheng et al., 2013;Prashant et al., 2015) and RVA parameters (McCartney et al., 2006;Sun et al., 2008;Zhao et al., 2009;Deng et al., 2015) were reported using different populations. However, these reports were primarily based on simple sequence repeat (SSR) markers, resulting in relatively large genetic distances defining QTL and markers due to the limited numbers of such markers, which leads to difficulty for gene cloning based on QTL analysis and marker-assisted selection (MAS) in wheat breeding. Therefore, it is necessary to develop higher density linkage maps for dissecting the genetic factors associated with these complex traits. It is widely accepted that the major effect QTLs for Mixograph and RVA parameters are associated with HMW-GSs and Wx genes, respectively. However, the associations of many other genes involved in dough rheological and starch pasting properties, and the genetic basis of these complex traits have not been well-documented. In addition, no QTLs for Mixolab parameters have been reported so far.
Gaocheng 8901, a bread wheat variety developed by the Gaocheng Agricultural Research Institute in Hebei province, has excellent dough strength and exhibits outstanding breadmaking quality. Zhoumai 16, developed by Zhoukou Academy of Agricultural Sciences in Henan province, is a high yielding, multiple disease resistant variety. Therefore, it was interesting to dissect the genetic components of yield and quality traits in these varieties. The aims of the present study were to: (1) construct a high-density linkage map with the 90 and 660 K SNP assays, (2) provide a comprehensive insight into genetic loci for dough rheological and starch pasting properties using a genome-wide QTL mapping approach, and (3) identify SNP markers closely linked to QTLs associated with quality traits for MAS in wheat breeding.

Evaluation of Flour Properties with Mixograph, Rapid Visco-Analyzer, and Mixolab
A 10-g Mixograph (National Mfg. Co., Lincoln, NE) was used to assess dough mixing characteristics and MPT, MPV, MPW, and MTxW were measured according to AACC (2000) method 54-40A. A RVA (Newport Scientific, Australia) was employed to evaluate the starch pasting properties of flour samples and PV, trough viscosity (TV), BD, final viscosity (FV), setback (SB), and peak time (PTI) were scored following AACC (2000) method 76-21 with minor modifications, viz., the reaction solution of water was replaced by 170 mg/L AgNO 3 to eliminate the effect of α-amylase activity in flour on starch pasting properties.
A Mixolab (Chopin Technologies, France) was used to determine dough mixing and pasting properties of wheat flour simultaneously during dough mixing. About 50-g of flour was put into the Mixolab bowl and an appropriate amount of water was added to ensure that the torque of the dough was in the 1.1 ± 0.07 Nm range. Processing was divided into five stages based on the "Chopin 12heat" protocol as follows: establishing equilibrium at 30 • C for 8 min, then heating to 90 • C at a rate of 4 • C/min for 15 min, holding at 90 • C for 12 min, cooling to 50 • C at a rate of 4 • C/min for 10 min, and finally holding at 50 • C for 5 min. The mixing speed was kept constant at 80 rpm. The parameters water absorption (WA), development time (DT), stability time (ST), C1 (the torque of maximum point in the first mixing stage), C2-C5 (the torque of end points in the corresponding mixing stages) were recorded during the procedure.

SNP Genotyping and Functional Marker Analyses
A 90 K iSelect SNP array containing 81,587 markers and a 660 K SNP array with 630,517 markers were used to genotype all 176 RILs and 59 randomly selected RILs, respectively, at CapitalBio Corporation (Beijing, China; http:// www.capitalbio.com). Genotypic clusters involving each SNP were confirmed using the polyploidy version of GenomeStudio software (Illumina, http://www.illumina.com).
Twenty-one functional markers identifying glutenin subunits, the 1B.1R translocation, kernel hardness and waxy alleles were used to screen the parents, and polymorphic markers were subsequently used to test the RILs. Sequences of PCR primers and fragment sizes are listed in Table S1.
Construction of a High-Density Linkage Map, QTL Analysis, and Identification of Candidate Genes SNPs of poor quality, or with more than 10% of missing data, or segregation distortion of more than 0.35 were removed. Three procedures were followed in constructing the high-density linkage map. Firstly, among 81,587 SNPs (90 K) used in screening all 176 RILs, 12,205 were polymorphic between the parents. Of those, 2080 had more than 10% missing data, and 828 were not anchored to the linkage map. Finally, 9297 SNPs were used to construct a linkage map (defined as Map 1) using Joinmap 4.0 software (Stam, 1993; http://www.kyazma.nl). However, many SNPs were mapped at the same loci or within 0.01 cM among them. Therefore, BIN-Mapping function from IciMapping 4.0 (Li et al., 2007; http://www.isbreeding.net) was used to construct a skeleton map containing 2375 markers for subsequent QTL analysis. Secondly, among 90 and 660 K SNP analyses of 59 random RILs,12,205 and 109,209,respectively, were polymorphic between the parents. Of those, 4809 and 20,498, respectively, had more than 10% missing data, and 2081 and 31,371, respectively, were not anchored to the linkage map. Thus, 5315 and 57,340 from two arrays, respectively, were used to construct the linkage map (defined as Map 2) with Joinmap 4.0 software; this was used to enrich the markers in several QTL regions with larger marker intervals from Map 1. Thirdly, Maps 1 and 2 were integrated into a high-density linkage map with MergeMap Online (Wu et al., 2011;http://www.mergemap.org).
The short and long chromosome arms of each linkage group were confirmed according to the wheat 90 and 660 K consensus SNP maps. The linkage map was constructed using MapChart 2.2 (Voorrips, 2002;http://www.earthatlas.mapchart. com). QTL analysis was carried out with QTL Cartographer 2.5 using composite interval mapping  http://statgen.ncsu.edu/qtlcart/WQTLCart.htm) based on Map 1. Logarithm of odds (LOD) scores ranged from 1.8 to 2.6 for all traits tested according to 2000 permutations tests at a probability of 0.01. Therefore, a LOD score of 2.6 was used for declaring significant QTL. The QTL × Environment (QE) interaction was analyzed by IciMapping 4.0 using the multienvironment trials (MET) (Li et al., 2007; http://www.isbreeding. net), with a LOD threshold based on 1000 permutation tests at a probability of 0.01. For the 90 K iSelect SNP genotyping assay, candidate genes were confirmed following Zhai et al. (2016). For the 660 K iSelect SNP genotyping assay, the sequences of flanking SNP markers tightly linked to QTL were blasted in the wheat genomes database (http://plants.ensembl.org/Triticum_ aestivum/Info/Index) to determine SNP markers corresponding to the original genes, and these gene sequences were used as queries to blast the NCBI database (http://www.ncbi.nlm.nih. gov/) to identify putative gene functions. BLAST hits were filtered to an e-value threshold of 10 −5 with an identity higher than 75%. Collinearity analysis was conducted by Ensembl Plants database (http://plants.ensembl.org/index.html) and blast hits were filtered with an e-value threshold of 10 −16 .

Statistical Analysis
Analyses of variance (ANOVA) and correlation were performed by SAS 9.0 (SAS Institute Inc., Cary, NC, USA). ANOVA was conducted using the PROC MIXED procedure, where environments were treated as fixed effects, and lines, line × environment interactions and replicates nested in environments were all treated as random effects. The broad-sense heritabilities (h 2 B ) of all the traits were calculated by h 2 B = σ 2 g /(σ 2 g + σ 2 ge /r + σ 2 ε /re), where σ 2 g , σ 2 ge and σ 2 ε were estimates of line, line × environment interaction and residual error variances, respectively, and e and r represented the numbers of environments and replicates, respectively. Pearson's correlation coefficients were calculated between traits with the PROC CORR procedure based on mean values across three environments.

Phenotypic Description
Gaocheng 8901 performed much better than Zhoumai 16 in Mixograph parameters and Mixolab dough rheological scores based on data in three environments (Table 1), whereas Zhoumai 16 showed better performance than Gaocheng 8901 for RVA parameters and Mixolab starch pasting properties. Wide variations were observed for all traits in the RIL population ( Table 1). The frequency distribution for all traits indicated continuous variation, with transgressive segregation (Figure S1), indicating polygenic inheritance.
ANOVA showed that genotype, environment and G × E interaction had significant effects on all traits except for MPW and PTI. Genotype and environment showed significant effects on MPW and PTI. Broad-sense heritabilities based on plot means ranged from 0.56 to 0.94 ( Table 2).

High-Density Linkage Map
The high-density linkage map comprised 8067 (90 K) and 38,894 (660 K) SNP markers, and spanned a total length of 4121 cM involving all 21 chromosomes, with an average chromosome length of 196.2 cM, ranging from 78.0 cM (3D) to 387.5 cM (5B) (Tables S2, S3). The A genome included 20,012 SNPs (42.6%) covering a length of 1439.8 cM and an average marker density of 0.07 cM; the B genome had 22,142 SNPs (47.1%) covering 1783.6 cM and an average marker density of 0.08 cM; the D genome included 4807 SNPs (10.2%), a length of 897.7 cM, and an average marker density of 0.19 cM. The number of SNP markers in each chromosome ranged from 43 (5D) to 6042 (4A). The marker intervals ranged from 0.02 (4A) to 2.14 (5D) cM, with a means of 0.09 cM.

QTLs for Quality Traits
A total of 119 additive QTLs were identified for 17 quality parameters using the 90 K iSelect SNP genotyping assay linkage map; 15 stable QTLs explained more than 10% of the phenotypic variation across environments (Table 4). In addition, 39 QTLs exhibited significant interaction with environments ( Table 4).

QTL Clusters
Some QTLs for different quality parameters were mapped in the same or nearby marker intervals, possibly due to a pleiotropic effect of a single gene or tightly linked genes. Eleven QTL clusters were identified on chromosomes 1AL, 1BL, 1DL, 2BS, 3B, 4B (2), 5AS, 5BL, 6BL, and 7DL based on Map 1 ( represented a gap of 10.0 cM in the linkage map from the 90 K assay. Additional six markers from the 660 K chip were newly mapped to this region and delineated this QTL to a 2.0 cM interval flanked by markers CAP8_c2839_118 and AX-110915310 ( Figure S2). The distance of the closest marker to the LOD contour peak of QMPV.caas.5AL was shortened from 1.1 to 0.7 cM by the additional marker from the 660 K assay. An additional two markers from the 660 K assay were mapped to the QTL cluster on chromosome 4B in the region of 47.4-49.7 cM, and the interval of QMTxW.caas.4B, QDT.caas.4B and QST.caas.4BS was reduced from 2.3 to 1.6 cM, flanked by markers AX-109563308 and AX-111689026. An additional two markers from the 660 K assay were newly mapped to the QTL cluster on chromosome 5AS in the region of 30.8-35.4 cM, and the interval for QFV.caas.5AS was shortened from 4.6 to 0.3 cM, flanked by markers RFL_Contig2251_434 and AX-108730091.

Phenotypic Evaluation
Mixograph and RVA are most commonly used to evaluate the dough rheological characteristics and starch pasting properties in wheat breeding, respectively. Recently, a new device Mixolab has been introduced to assess dough rheological characteristics, starch pasting properties and flour enzyme activity in one test sample, reducing labor requirements. The instrument also provide information on effects of different ingredients (Bonet et al., 2006;Rosell et al., 2007;Marco and Rosell, 2008), indicating that it could be used in areas of the bakery products industry. Since, it is a relatively new device, the correlations between Mixolab parameters and other traditional instruments such as Mixograph and RVA have been not completely established. The  significant correlations between MPT and Mixolab parameters for dough properties such as DT (r = 0.80), ST (r = 0.85), and C2 (r = 0.60) were found in the present study. FV also exhibited significant correlations with C3 (r = 0.31), C4 (r = 0.44), and C5 (r = 0.38). These results suggested that wheat quality can be effectively assessed by Mixolab parameters. However, it needs about 50-g flour in one Mixolab test, whereas only 10 and 3.5-g of flour were needed in one Mixograph and RVA analysis, respectively. As a result, Mixograph and RVA are more suitable for analyzing early generational material in wheat quality breeding while Mixolab could be used in bakery factories. Although the present study showed genotype, environment and QE interaction had significant effects on all traits except MPW and PTI, the broad-sense heritabilities assessed for all traits ranged from 0.56 to 0.94, suggesting that genotype had the largest contribution to source of variation, in agreement with Patil et al. (2009).

High-Density Linkage Map
Linkage maps play a crucial role in identifying QTL, cloning genes, MAS and genome structure analysis (Maccaferri et al., 2014). The QTL mapping was based on the genotyping data with the 90 K array. Because a low coverage of SNP markers on some chromosomal regions, several QTLs were mapped in larger marker intervals. Thus, the 660 K SNP array was used to genotype two parents and 59 RILs randomly chosen from the population to saturate the chromosomal regions with larger marker intervals. In the present study, we firstly obtained a high-density linkage map based on the 90 and 660 K iSelect SNP assays, including 46,961 polymorphic SNP markers. Of these markers, 6596 (1046 from the 90 K SNP assay and 5550 from the 660 K SNP assay) were newly anchored to the bread wheat linkage map (Table S4), by comparing with the 90 and 660 K consensus maps (Wang et al., 2014;Prof. Jizeng Jia, pers. comm.). The average density of the map was 0.09 cM/marker, indicating a higher markerdensity than developed previously with DArT (Echeverry-Solarte et al., 2015) or SNP from the 90 K assay (Gao et al., 2015). The marker-densities in regions surrounding important QTLs (QC3.caas.3AS.2 and QMPV.caas.5AL) and QTL clusters (4B, 5AS, 5BL, and 7DL) were significantly increased by the high-density linkage map ( Figure S2); consequently increasing the resolution of QTL mapping. It is interesting that the original marker orders of these regions, except for chromosome 7DL from the 90 K SNP assay, are unchanged in the high-density linkage map, indicating that the regions are conserved. Although a good coverage of the genomes was obtained, polymorphism in the D genome was still relatively low. In summary, this highdensity linkage map is valuable for fine mapping, candidate gene discovery and MAS in wheat breeding.

Comparison with Previous Studies
It was proven that HMW-GSs showed larger contribution to dough properties (Branlard et al., 2001;He et al., 2005;Liu et al., 2005). In the present study, the major effect QTLs for Mixograph parameters on chromosomes 1AL, 1BL, and 1DL should be conferred by HMW-GSs flanked by BS00030036_51 and Excalibur_c44668_382, Glu-B1 and wsnp_Ra_c7527_12935330, BS00018250, and Glu-D1, respectively, in agreement with previous reports (Nelson et al., 2006;Mann et al., 2009).
The stable QTL on chromosome 2BS with Gaocheng 8901 allele increased MPT (QMPT.caas.2BS.1) and MTxW(QMTxW.caas.2BS). Zhang et al. (2009) indicated that an appropriate ratio of quantity of glutenin to gliadin had larger contribution to dough Mixograph properties. This QTL may be corresponding to the region related to the quantity of gluten or gliadin protein fractions in grain protein using a genome-wide association study of a bread wheat world core collection (Plessis et al., 2013). QMPV.caas.3BL, with allele from Zhoumai 16 increasing MPV, is in a similar position on chromosome 3BL that influenced glutenin macropolymer content, Zeleny sedimentation volume and grain protein content (Sun et al., 2008). This coincides with previous study that there was significant correlation between protein content and MPV (Tronsmo et al., 2003).
Gaocheng 8901 (HI = 66) with alleles Pina-D1b/Pinb-D1a/Pinb-2v2 showed slightly harder than Zhoumai 16 (HI = 53) with Pina-D1a/Pinb-D1b/Pinb-2v3. Grain hardness mainly affected by Ha locus on chromosome 5DS exhibited significant FIGURE 1 | Genetic maps and QTLs for quality traits identified in the Gaocheng 8901/Zhoumai 16 population from the 90 K iSelect SNP genotyping assay. Single chromosome was in red and blue region on the chromosome indicated the confidence interval of QTLs detected. QTL names were on the right with different colors for different traits. See footnote to Table 1 for abbreviations.
influence on milling quality. It has been reported that hard wheat has much more the amount of damaged starch than soft wheat (Barak et al., 2012). The association analysis of genotypic and phenotypic data using T-test indicated that both Pina and Pinb conferred significant effect on WA (P < 0.01). Therefore, QTL for WA on chromosome 5DS was contributed by both Pina and Pinb genes at Ha locus, in agreement with Ma et al. (2007).
One stable QMPT.caas.5AL positioned at 87 cM on chromosome 5AL was not previously reported. A stable QTL QMPT.caas.7DL at 93 cM on chromosome 7DL is different from one reported by Tsilo et al. (2011)  Functional genes, transporters and transcription factors associated with starch metabolism were summarized by Singh et al. (2015). The sequences of these genes were used as queries against the T3 marker database (Zhai et al., 2016;http:// triticeaetoolbox.org/wheat/viroblast/viroblast.php) to identify SNP markers corresponding to the original genes, and then the genetic map constructed in the present study was inspected for the presence of the same markers. Kukri_rep_c101946_496 derived from an isoamylase 2 gene, was mapped on chromosome 1AL at a distance of 2.45 cM from the LOD contour peak marker for QSB.caas.1AL. Wsnp_Ex_rep_c66900_65313836 derived from an isoamylase 3 gene, was mapped on chromosome 5AL at a distance of 9.50 cM from the LOD contour peak marker of QFV.caas.5AL (Table S5). The QTL on chromosome 7D is not a Wx gene since both parents were wild type at all three Wx loci.
Although 119 additive QTLs were detected for dough rheological and starch pasting properties, 50 of them were grouped into 11 clusters. QTL clusters for dough strength on chromosomes 1AL, 1BL, and 1DL detected in the present study were associated with HMW-GSs, consistent with previous reports (Huang et al., 2006;Deng et al., 2013). The interval 24.9-34.2 cM on chromosome 2BS is a QTL cluster impacting MPT, MTxW, WA, ST and C2, which coincides with their significant phenotypic correlations. A similar co-localized region for dough development time, stability and softening was also located on chromosome 2B (Maphosa et al., 2015). A QTL cluster for MPV, MPW and PTI between BS00072151_51 and GENE-1617_131, positioned in the interval of 61.3-81.9 cM on chromosome 3B, is likely to be in the same region for MPT, MPI, and MTxW reported by Deng et al. (2013). Another QTL cluster for MPT, MTxW, DT, ST, C2, and C3 identified in the interval 35.8-63.1 cM on chromosome 4B, is likely the same as a QTL cluster for Mixograph midline peak time, midline tail width, envelope peak energy and envelope peak time (Patil et al., 2009). The interval 83.5-89.7 cM on chromosome 4BL is a QTL cluster influencing MPT, DT and ST, in agreement with their phenotypic correlations. This QTL cluster is reported for the first time.
The QTL cluster for BD, FV, ST, and C2 between Kukri_c18268_79 and Kukri_c60091_331 positioned in the interval 20.1-35.8 cM on chromosome 5AS is reported for the first time. Another QTL cluster for PV, FV, SB, C3, C4, and C5 in the interval 167.9-177.8 cM on chromosome 5BL is also new. The QTL cluster for MPV, PV, TV, WA, and C4 was in the interval 55.1-73.7 cM on chromosome 6BL, is in a similar pleiotropic region for Farinorgraph dough stability, PV, TV, FV, and PTI between Xgwm644 and Xgwm608b reported previously (Sun et al., 2008). The QTL cluster for MPT, BD, and FV in interval 86.3-103.5 cM on chromosome 7DL is new. The presence of QTL cluster may be attributed to either pleiotropic effects of a single QTL or a few loci closely linked. Most of traits involved in a QTL cluster were significantly correlated ( Table 3).

Potential Implication in Wheat Breeding
KASP is used to detect InDels or SNPs. It is suitable for genotyping several SNP markers for a number of samples. SNP markers closely linked to traits identified by genomewide association study or QTL mapping will be successfully converted into KASP assay for marker-assisted breeding in the future. In this study, the stable QTLs with larger contribution such as QMPT. The development of a high density linkage map with the wheat 90 K array and comparative genomics provide a powerful tool in searching for potential candidate genes in wheat. Bioinformatics analysis of the mapped SNP markers in the important QTL regions for dough rheological and starch pasting properties identified eight candidate genes. Of them, four genes were involved in biosynthesis of amino acids, two related to starch and sucrose metabolism, and two associated with fatty acid biosynthesis (Table S6). However, since a number of biological processes are associated with these candidate genes, more detailed experimental analyses will be needed to confirm their roles in dough strength and starch pasting properties.
Although there are many reports on QTL for quality traits in wheat (Huang et al., 2006;McCartney et al., 2006;Sun et al., 2008;Patil et al., 2009;Zhao et al., 2009;Tsilo et al., 2011;Li et al., 2012;Deng et al., 2013Deng et al., , 2015Zheng et al., 2013;Prashant et al., 2015), few genes are identified in single QTL analysis. The main reason is that SSR markers previously used in QTL analysis have insufficient density and are often in noncoding regions. The 90 and 660 K iSelect SNP genotyping assays are developed from transcriptome and genome sequencing, respectively (Wang et al., 2014;Prof. Jizeng Jia, pers. comm.), and therefore provide the better approach in searching for candidate genes.
With rapid development of next-generation sequencing, genomic sequence of many species such as rice (Oryza sativa L.) (Goff et al., 2002), brachypodium (Brachypodium distachyum L.), (The International Brachypodium Initiative, 2010), wheat (Triticum aestivum L.) (Brenchley et al., 2012), barley (Hordeum vulgare L.) (The International Barley Genome Sequencing Consortium, 2012) were released on web, which provides possibility for comparative genomics analysis. Syntenic relationships among grass species are served as a tool for fine mapping and predicting gene function (Sorrells et al., 2003;Ouyang et al., 2014). As shown in Figure S3, collinearity analysis on QTL confidence interval indicated that the collinearity between wheat and barley was better than Brachypodium and rice, and the homologies among chromosomes were consistent with previous reports (Sorrells et al., 2003;The International Brachypodium Initiative, 2010 , and QWA.caas.7AS) showed the best conservation of gene order among these species, suggesting that these regions may facilitate further fine mapping and discovery of candidate genes.

ETHICS STATEMENT
We declare that these experiments comply with the ethical standards in China.

AUTHOR CONTRIBUTIONS
HJ carried out the experimental and wrote the paper. WW, JL, YZ, and JY participated in field trials. SZ contributed to flour milling. ZL, XX, and ZH designed the experiment and assisted in writing the paper.

ACKNOWLEDGMENTS
The authors are grateful to Prof. R. A. McIntosh, at Plant Breeding Institute, University of Sydney, for review of this manuscript. This study was supported by the National Natural Science Foundation of China (31461143021, 31371623), Beijing Municipal Science and Technology Project (D151100004415003), and China Agriculture Research System (CARS-3-1-3).
Table S1 | PCR primers for gene-specific markers used for genotyping the parents and RILs.