Quantitative Trait Locus (QTLs) Mapping for Quality Traits of Wheat Based on High Density Genetic Map Combined With Bulked Segregant Analysis RNA-seq (BSR-Seq) Indicates That the Basic 7S Globulin Gene Is Related to Falling Number

Numerous quantitative trait loci (QTLs) have been identified for wheat quality; however, most are confined to low-density genetic maps. In this study, based on specific-locus amplified fragment sequencing (SLAF-seq), a high-density genetic map was constructed with 193 recombinant inbred lines derived from Chuanmai 42 and Chuanmai 39. In total, 30 QTLs with phenotypic variance explained (PVE) up to 47.99% were identified for falling number (FN), grain protein content (GPC), grain hardness (GH), and starch pasting properties across three environments. Five NAM genes closely adjacent to QGPC.cib-4A probably have effects on GPC. QGH.cib-5D was the only one detected for GH with high PVE of 33.31–47.99% across the three environments and was assumed to be related to the nearest pina-D1 and pinb-D1genes. Three QTLs were identified for FN in at least two environments, of which QFN.cib-3D had relatively higher PVE of 16.58–25.74%. The positive effect of QFN.cib-3D for high FN was verified in a double-haploid population derived from Chuanmai 42 × Kechengmai 4. The combination of these QTLs has a considerable effect on increasing FN. The transcript levels of Basic 7S globulin and Basic 7S globulin 2 in QFN.cib-3D were significantly different between low FN and high FN bulks, as observed through bulk segregant RNA-seq (BSR). These QTLs and candidate genes based on the high-density genetic map would be beneficial for further understanding of the genetic mechanism of quality traits and molecular breeding of wheat.


INTRODUCTION
Wheat (Triticum aestivum L.) is an important cereal crop, being one of the main sources of food for approximately 35% of the world's population (FAOSTAT, 2013;Curtis and Halford, 2014;Liu J. et al., 2017). Currently, wheat is widely consumed and processed into bread, noodles, cakes, pasta, beer, and other products. Enhancement in the quality and yield of wheat has significant impacts on food security and human health. Improving end-use quality to meet the increasing market demand is becoming a critical topic in crop sciences (Kong et al., 2013;Amiri et al., 2018). The quality traits of wheat are controlled by many genes and are easily influenced by the environment. Quantitative trait locus (QTL) analysis and genome-wide association studies are regarded as useful tools for understanding the genetic mechanism and for identifying markers for marker-assisted selection (MAS) of quantitative traits of wheat.
Grain hardness (GH), which defined endosperm texture, is the most important determinant for wheat quality and classification (Pasha et al., 2010). Difference in GH is suggested to be due to the continuity of the protein matrix and the strength with which it physically entraps starch granules (Pomeranz and Williams, 1990). Changes in GH affect milling, baking, and other end-use applications of wheat. Hard wheat is suited for fermented products such as bread, because starch granules were liable to be broken under stronger grinding forces, and these broken starch granules with higher water-absorbing capacity and enzymatic digestion rate are beneficial for yeast growth and fermentation (Henry, 1996;Pasha et al., 2010). In contrast, for soft wheat, more intact starch granules remain during the milling process; soft wheat is usually used in the production of biscuits, cakes, pastries, and confectionaries. The major QTL for GH is located at the Ha loci of chromosome 5DS (Mattern, 1973;Law et al., 1978). Pina-D1 and Pinb-D1 genes are close to the Ha loci, encoding "friabilin, " which is an endospermspecific lipid-binding protein that has been proven to largely determine the GH of wheat (Morris, 2002). However, up to 40% of the variation in GH was due to other unknown factors (Weightman et al., 2008).
Grain protein content (GPC) is another key factor determining the nutritional value and end-use quality of wheat. For example, wheat with higher GPC is suitable for making bread, pasta, and yellow alkaline Chinese noodles (Ramen); however, for white-salted Japanese noodles (Udon), wheat flour with low or moderate GPC is preferred (Distelfeld et al., 2008). In general, GPC is a complex trait controlled by multiple quantitative QTL/genes (Tsilo et al., 2010;Simons et al., 2012;Kumar et al., 2018) and largely influenced by the environment (Chope et al., 2014;Fatiukha et al., 2019). Approximately 367 QTLs identified for GPC have been mapped on all chromosomes of wheat, explaining the phenotypic variance of 0.6-66%. Most major QTLs have been detected on chromosomes 2A, 2B, 3A, 4A, 6B, 7A, and 7B (Kumar et al., 2018). Gpc-B1 on chromosome 6BS, explaining up to 66% of the phenotypic variation for GPC (Olmos et al., 2003), was further fine mapped to an NAC transcription factor (NAM-B1), which affects nutrient remobilization (Uauy et al., 2006), is an important QTL, and is deployed in several breeding programs (Kumar et al., 2018).
The Hagberg falling number (FN), which reflects α-amylase activity in mature grains and simulates the rheological properties of starch in cooking, is an important grading index for wheat quality. Wheat flour with excessively low FN (LFN) related to higher α-amylase activity usually produces a sticky dough and discolored and poorly structured loaves (Zhang et al., 2014;Newberry et al., 2018). Elevated levels of α-amylase activity in mature grains are usually associated with preharvest sprouting (PHS) and late maturity α-amylase (LMA), which downgrades the grain (Munkvold et al., 2009). PHS refers to grain germination in the spikes of the mother plant before harvest especially during rainy days, which reduces not only the yield, but also nutrition and industry quality (Gubler et al., 2005;Newberry et al., 2018). LMA is a genetic defect with the characteristic of a single α-amylase isoform 1 (α-Amy-1) during later stages of grain development (Mares and Mrva, 2008;Newberry et al., 2018). QTLs for PHS and LMA rated by FN were mapped to almost all wheat chromosomes (Zanetti et al., 2000;Groos et al., 2002;Osa et al., 2003;Kulwal et al., 2004Kulwal et al., , 2011Mares et al., 2005;Mori et al., 2005;Kottearachchi et al., 2006;Fofana et al., 2008;Mohan et al., 2009;Munkvold et al., 2009;Jaiswal et al., 2012;Cabral et al., 2014), of which most QTLs were distributed on chromosome 3A, 3B, 3D, 6B, 7B, 4A, and 4B and were related to plant height, grain color, and grain dormancy. However, the mechanism for PHS and LMA remains unclear, with few genes cloned.
Starch pasting properties usually receive far less attention than other quality traits, but pasting parameters measured by rapidviscosity analysis (RVA) can simulate flour changes in the cooking process and are closely associated with Asian noodle and steamed bread quality Jin et al., 2016). The Wx gene located on chromosomes 7A, 4A, and 7D responsible for amylose synthesis has large effects on starch pasting properties (Nakamura et al., 1993;Hung et al., 2007;Blazek and Copeland, 2008;Jin et al., 2016). Although many QTLs for pasting properties were detected on all chromosomes (Bao et al., 2000;Kuchel et al., 2006;Zhao et al., 2009;Deng et al., 2015;Wang et al., 2017), few key genes have been cloned and characterized.
Thus, although many QTLs have been identified to determine the genetic basis of quality traits in wheat, the average distances of most genetic maps are too large, or markers are mainly distributed on partial regions of chromosomes. Therefore, rare genes were cloned based on QTL analysis, and few molecular markers can be developed for MAS in wheat breeding. Therefore, high-density genetic maps are necessary for the genetic dissection of complex quality traits of wheat.
High-throughput single-nucleotide polymorphism (SNP) genotyping plays an important role in constructing high-density genetic maps of wheat (Wang et al., 2014;Maccaferri et al., 2015). Specific-locus amplified fragment sequencing (SLAF-seq) based on restriction-site associated DNA tag sequencing (RAD-seq) was developed for economic and efficient high-throughput SNP discovery (Sun et al., 2013;Qi et al., 2014). SLAF-seq can provide abundant InDel and SNP markers and has been applied in highdensity genetic map construction and candidate functional gene identification for many species in recent years Han et al., 2019;Ma et al., 2019;Sun et al., 2019;Zhuang et al., 2019;Wang et al., 2020). In addition, for wheat, as the genome sequencing has been completed, a predesign experiment in SLAF-seq can be performed to evaluate restriction enzymes and sizes of restriction fragments, which have improved the efficiency of SLAF-seq. SLAF-seq has been mainly used for the identification of desirable genes on alien chromosomes or for yield traits such as thousand seed weight of wheat; however, there has been no report on the quality traits of wheat based on SLAFseq yet (Hu et al., 2016;Li G. et al., 2016;Yin J. L. et al., 2018).
In this study, a high-density genetic map was constructed based on recombinant inbred lines (RILs) derived from Chuanmai 42 and Chuanmai 39 through SLAF-seq, and QTLs for quality traits were identified. Chuanmai 42, with high yield and good resistance, is a backbone parent used in wheat breeding in southwest China. Chuanmai 39 has high protein content and strong gluten. The QTL and candidate genes would be beneficial for MAS in wheat breeding for quality.

Plant Materials and Field Trials
The RILs, including 193 lines derived from Chuanmai 42 (♀) × Chuanmai 39 (♂), were applied to construct genetic maps. Chuanmai 42 is a soft and red-grained wheat with high yield and resistance. Chuanmai 39 is a hard and white-grained wheat with higher gluten strength. The population was grown in Shuangliu, Sichuan Province, China, in 2016(E1) and 2017(E2), and in Shifang, Sichuan Province, China, in 2017. The double-haploid (DH) population, including 376 lines derived from Kechengmai 4 (♀) × Chuanmai 42 (♂), was applied for QTL validation and was planted in Shifang, Sichuan Province, China, in 2018-2019. Each line was planted in 2 blocks; each block was 2 m long, with three rows 30 cm apart. Crop management was implemented according to local cultivation practices. Grains were collected separately by each block and dried naturally. Grains were milled using FOSS Cyclotec CT1093.

Phenotyping
Grain protein content was detected by near-infrared reflectance spectroscopy on a Perten DA-7200 instrument (Perten Instruments, Huddinge, Sweden). GH was determined using a single-kernel characterization system 4100 (Perten Instruments, Springfield, IL, United States). Grains were ground to whole meal using a 1-mm-sieve Cyclotec mill (Foss Tecator AB, Höganäs, Sweden). FN was measured using Falling Number FN 1000 (Perten, Sweden). Pasting properties parameters were analyzed by Micro Visco-Amylo-Graph (A. W. Brabender Instruments, South Hackensack, NJ, United States) as follows: whole meal (15 g) was suspended in distilled water (98 mL) and 10% (m/v) AgNO 3 (2 mL). The suspension was then heated from 50 • C to 92 • C at a rate of 7.5 • C/min, held at 92 • C for 5 min, cooled to 50 • C at 7.5 • C/min, and held at 50 • C for 1 min; the rotation speed was 250 rpm. α-Amylase activity was determined according to Whan et al. (2014) and expressed in the Ceralpha unit per gram four as determined by Bradford assays (Bradford, 1976) on the CERALPHA extracts. All traits were measured for two replicates of each block of one line.

SLAF Library Preparation and Genotyping
Total DNA was extracted from seedlings according to the CTAB procedure (Wang, 1997), quantified using NanoDrop 2000, and stored at −80 • C. SLAF library was prepared and genotyped according to Sun et al. (2013). The DNA digestion sites and the length and distribution of the resulting fragments were investigated using the wheat reference genome (IWGSC RefSeq v1.0 1 ), and RsaI was chosen as the appropriate restriction enzyme. Oryza sativa L. was used to control SLAF library preparation. Fragments of 464-484 bp were selected to generate paired-end reads on an Illumina HiSeq-2000 sequencing platform (Illumina, San Diego, CA, United States) at Biomarker Technologies Corporation in Beijing 2 . Then, low-quality reads (quality score B20e) were filtered out and sorted to each line according to duplex barcode sequences using SLAF_Poly.pl software (Biomarker, Beijing, China). High-quality reads were barcoded, terminal 5 bp were removed, and clean reads were mapped to the wheat reference genome (IWGSC RefSeq v1.0 1 ) using SOAP software (Li et al., 2008).
Genotyping was performed using the Bayesian approach to ensure quality (Sun et al., 2013). First, a posteriori conditional probability was calculated using the coverage of each allele and the number of SNPs. Then, the genotyping quality score translated from the probability was used to select qualified markers for subsequent analysis. Low-quality markers for each marker and each individual were counted and the worse markers or individuals were deleted during the dynamic process. When the average genotype quality scores of all SLAF markers reached the cutoff value, the process was stopped. Only polymorphic SLAF tags having full parental homozygosity (aa × bb) were chosen. Then, high-quality SLAF markers were obtained for high-density linkage map construction referring to the following filter standard. First, average sequence depths should be > 2-fold for each line and > 4-fold for the parents. Second, markers with more than 90% (A and B genome) or 75% (D genome) missing data were removed. Third, markers with significant segregation distortion (P < 0.01) were excluded.

Construction of a High Genetic Linkage Map
Based on the genotyping results of the 193 RILs and the wheat reference genome (IWGSC RefSeq v1.0 1 ), a genetic linkage map including 21 linkage groups was constructed using HighMap software . The genetic distances of marker and genotypes of each RIL are listed in Supplementary Table 1.
The physical positions were obtained by referring to the wheat reference genome, and the collinearity between the genetic and physical positions was measured using the Spearman correlation coefficient.

QTL Mapping
Additive QTL mapping (individual environment) was performed with the package "R/qtl" in R using Haley-Knott regression. The significance of the mapped QTLs was determined at an experimental probability of error P < 0.05, using genome-wide LOD thresholds ≥ 3. The best linear unbiased predictors (BLUPs) were calculated for quality traits over different environments using SAS. The walking step was set to 1 cM. Genetic maps of chromosomes with significant QTL were drawn using MapChart version 2.1.

QTL Validation
Quantitative trait loci QFN.cib-3D was selected for validation in the DH derived from Kechengmai 4 × Chuanmai 42 containing 376 lines. As shown in Supplementary Table 2, four pairs of primers (3D55982, 3D55956, 3D56098, and 3D56055) were designed for QFN.cib-3D based on the four flanking SLAF markers (SLAFs) (M56098, M56055, M56098, and M56055). Parents and 310 lines from DH were randomly chosen for amplification with the four pairs of primers. At least three amplifications were applied to each line. The amplified bands were sequenced by TSINGKE Biological Technology Company (Beijing, China). Then, the SNP differences in the amplified sequence were compared among parents and DH lines, to identify the effect of QFN.cib-3D on FN.

Annotation of Genes Within QTL Region and Comparison With Previous Studies
Flanking SLAF markers of QTLs stable in at least two environments were used to blast against IWGSC RefSeq v1.0 3 to obtain the physical locations and were compared to previously reported QTLs. The sequences of previously reported QTLs were obtained from GrainGenes 4 or "iwgsc_refseqv1.0_Marker_mapping_summary_2017Mar13/" downloaded from https://urgi.versailles.inra.fr/download/iwgsc/. Genes within the QTL region were retrieved from CDS sequences in IWGSC_RefSeq_Annotations_v1.0 5 . Gene function was analyzed using UniProt, COG, GO, Swissprot, and KEGG.

Bulked Segregant Analysis by RNA-Seq
Referring to FN in three environments, 10 lines of RILs derived from Chuanmai 42 and Chuanmai 39 with extremely high FN (HFN) and LFN from RILs were chosen for bulk segregant RNAseq (BSR). Total RNA was extracted from grains at 35 days after flowering using the TRIzol method for each chosen line. Then, equal amounts of RNA from the 10 lines with extremely high FN were mixed as HFN bulk, and equal amounts of RNA from 10 lines with extremely low FN were mixed as LFN bulk. Two libraries including the two bulks were constructed using the NEBNext Ultra RNA Library Preparation Kit (New England Biolabs, United States) and sequenced through Illumina HiSeq TM 2500 platform (Illumina, United States). Reads with adaptors or more than 10% unknown nucleotides and lowquality reads (> 50% bases with a quality score ≤ 20) were removed. The filtered reads were aligned to IWGSC RefSeq v1.0 using HISAT2.0.4.
Gene expression differential display was analyzed using the DEGseq package (version 1.18.0). The statistical significance of differentially expressed genes (DEGs) was determined using a combination of multiple tests and false discovery rate (FDR). Genes with FDR < 0.05 were classified as significant DEGs. Gene function was analyzed using UniProt, GO, KEGG, and Mapman. GO and KEGG pathway enrichment analyses were performed using GO Seq 2.12 and KOBAS v2.0, with FDR < 0.05.
Single nucleotide polymorphisms were obtained initially using the GATK package (Genome Analysis Toolkit, v3.2-2; Broad Institute, United States). These SNPs were then filtered based on the following criteria: (1) sequencing depth for each SNP ≥ 7; (2) SNP index of HFN and LFN is both lower than 0.3. The SNP index value was calculated using the MutMap method (Abe et al., 2012). Then, the SNP index for each SNP was calculated through the following formula: SNP index = |(SNP index of HFN bulk) -(SNP index of LFN bulk)|. The average value for the SNP index in the corresponding window was calculated using a sliding window with a window size of 1 Mb and slides at a size of 1 kb. SNPs with SNP index > 0.80 in candidate regions were considered as candidate loci related to FN.

Statistical Analysis
Analysis of variance was conducted by the general linear model using SPSS Statistics 20. The phenotypic variance of each trait included nγσ 2 G + γσ 2 GE + σ 2 e , where n is the number of environments, r is the number of blocks, σ 2 G is the genetic variance, σ 2 GE is the genotype × environment variance, and σ 2 e is the error variance. The entry-based broad sense heritability (H 2 ) was measured by H 2 = σ 2 G /(σ 2 G + σ 2 GE /n + σ 2 e /nγ). The significance test and Pearson correlation were analyzed by SPSS Statistics 20.

Phenotypic Variation of Quality Traits
Grain protein content, GH, FN, and starch pasting parameters exhibited continuous variation in each environment, of which GH showed a bimodal distribution (Supplementary Figure 1). GPC, GH, and FN of the male parent Chuanmai 39 were consistently and significantly higher than those of the female parent Chuanmai 42 in all three environments ( Table 1). Chuanmai 42 was characterized by higher peak viscosity (PV), through viscosity (TV), and final viscosity (FV), but no significant difference in starch pasting parameters was found between parents ( Table 1). The RILs exhibited high variations among these quality traits, of which GH showed the highest coefficient of variation (CV) across the three environments, from 40.62% (E1) to 42.37% (E3), and peak time (PT) showed the lowest CV ranging from 4.36% (E2) to 5.21% (E3) ( Table 1).
The broad-sense heritability (H 2 ) values of FN, GH, and GPC were 0.78, 0.71, and 0.75, respectively, which indicated that phenotypic variations were mainly due to genetic differences, although genotype, environment, and genotype × environment (G × E) interaction had significant effects on all traits (Supplementary Table 3).
The Pearson correlation was analyzed among quality traits (Supplementary Figure 2). Higher correlation was detected among pasting parameters especially PV, TV, and FV (r = 0.482 * * -0.975 * * ), in the three environments. In terms of the correlation between pasting parameters and GPC, FN, and GH, higher and significant correlations were found in E1 as follows: r = 0.401 * * -0.697 * * , r = 0.519 * * -0.722 * * , and r = 0.254 * * -0.559 * * . The correlation between pasting properties and GPC and FN was slightly higher than that between pasting properties and GH. Among these pasting properties, PT and FV were more relevant to GPC, FN, and GH with higher correlation coefficients (r = 0.551 * * -0.697 * * ). In terms of the correlation among GH, FN, and GPC, a higher and significant correlation was found in E1. For E1, a relatively higher correlation was detected between GH and FN (r = 0.722 * * ) than between FN and GPC (r = 0.501 * * ) and between GPC and GH (r = 0.431 * * ).

SLAF-seq and Genotyping of RILs Population
In total, 233.90 Gb of raw data, including 1,169.49 M pair-end reads, were obtained by high-throughput sequencing of the SLAF libraries. The average Q30 ratio and guanine-cytosine content for parents and progenies were 90.61% and 43.70%, respectively (Supplementary Table 4). A total of 397,199 high-quality SLAFs harboring 4,448,250 SNPs were detected. After removing markers lacking polymorphisms between parents and low-quality markers with average sequence depths less than 4-fold, 114,140 SNPs were retained and classified into eight segregation patterns (aa × bb, ab × cc, cc × ab, ef × eg, hk × hk, lm × ll, nn × np). Only aa × bb pattern including 88,234 SNPs was used for construction of the genetic map.
The A genome included 3,759 SLAFs (29.66%), covering 859.89 cM with an average marker density of 0.30 cM. The B genome contained 6,887 SLAFs (53.96%) and spanned 1,022.87 cM with an average interval of 0.20 cM, and the D genome included 2,386 SLAFs (15.89%) covering 977.18 cM with an average interval of 0.66 cM (Supplementary Table 5). SLAFs are concentratedly distributed on the B genome, and unfortunately, the polymorphisms detected in the D genome were still relatively low.
Collinearity analysis results of the SLAFs between the genetic map and the physical map of wheat are shown in  Table 6). This indicated that the genetic map constructed by SLAFs had a sufficient coverage over the wheat genome, and the majority of SLAFs on the linkage map were of the same order as those on the corresponding chromosomes of the physical map of the wheat genome.

Additive QTL Analysis for Quality Traits
A total of 30 QTLs were identified for FN, GPC, GH, and starch pasting properties ( Table 2, Supplementary Table 7 Table 2,  Supplementary Table 7, and Figure 2). Five of the 30 QTLs were stable in at least two environments also detected by the BLUPs (Table 2 and Figure 2). The average phenotypic variance explained (PVE) by individuals of the 5 QTLs varied from 4.32% to 25.74% ( Table 2).
For FN, the number of QTLs (8) and major or stable QTLs (3) detected were much higher than other quality traits ( Table 2 and Supplementary Table 7). All QTLs conditioned higher FN value through Chuanmai 39 alleles except QFN.cib-3D, with PVE ranging from 2.85 to 25.74% ( Table 2 and Supplementary  Table 7). Three QTLs on chromosomes 3A, 2B, and 3D could be detected in two environments and by BLUP, of which QTL qFN-3D.1 and qFN-3D.2 had the highest PVE of 16.58%-25.74% with additive effects ranging from 25.57 to 56.63 s ( Table 2).
Five QTLs were detected for GPC, of which only QGPC.cib-4A located on chromosome 4A was detected in all three environments and by BLUP, explaining 5.11-12.24% phenotypic variations ( Table 2 and Supplementary Table 7). All these QTLs resulting in higher GPC content were contributed by Chuanmai 39 alleles ( Table 2 and Supplementary Table 7).
Only one QTL located on the ends of chromosome 5D (QGH.cib-5D) was detected for GH with PVE reaching 33.31-47.99% (Table 2). This QTL contributing to harder grains with a hardness index increase by 13.17-15.62 was derived from Chuanmai 39 alleles and was stable in all three environments and detected by BLUP ( Table 2).
In total, 16 QTLs with a PVE of 1.89-17.38% were detected for starch pasting parameters except FV, including 5, 2, 1, 6, and 2 QTLs for PV, TV, PT, BD, and SB, respectively, but no one stable QTL was detected in at least two environments (Supplementary Table 7).

Effects of QTL Combination on FN
To determine the combination of three major or stable QTLs (QFN.cib-3A, QFN.cib-2B, and QFN.cib-3D) on FN, the RILs were grouped into eight genotypes by using the flanking markers for each QTL (Figure 3A). Each genotype represents one QTL combination, which refers to a group of alleles from different loci that are inherited from parents and expressed in the progeny. Genotypes containing fewer than three lines were not analyzed.

QTL Effect Confirmation and a PCR-CTPP Marker Development
The effect of QFN.cib-3D was also verified in the DH population containing 376 lines derived from Kechengmai 4 × Chuanmai 42. The parent lines were polymorphic for four pairs of primers (Supplementary Table 2) designed for the four flanking SLAFs of these QTLs, and the RILs were classified into two groups through amplification and sequencing. A total of 310 lines were randomly chosen for QFN.cib-3D validation, and the four primers (3D56098, 3D56055, 3D55982, 3D55956) were closely linked in RILs derived from Chuanmai42 and Chuanmai39. Of these, 147 lines with the same alleles as Kechengmai 4 or Chuanmai 39 displayed significantly lower FN (277.03 ± 59.62 s) than the other lines (300.57 ± 45.94 s) with the same allele  Frontiers in Plant Science | www.frontiersin.org FIGURE 2 | Distribution of QTLs on chromosomes. Orange-marked QTLs were detected in E1, blue ones were detected in E2, purple ones were detected in E3, and black ones were detected by BLUP. Underlined QTLs denote those that could be detected in at least two environments or by BLUP.
as Chuanmai 42 ( Figure 3B). This indicated that for QFN.cib-3D, the alleles from Chuanmai 42 had positive effects on maintaining HFN. A polymerase chain reaction with a confronting two-pair primers (PCR-CTPP) marker was designed for QFN.cib-3D to detect its presence or absence. The same allele as Chuanmai 42 was amplified with two DNA fragments of 483 and 787 bp, whereas the same allele as Kechengmai 4 or Chuanmai 39 was amplified with 360 and 787 bp (Figures 3C,D). The marker profile was in accordance with the sequencing results.

Identification of SNPs and Candidate Genes for FN Based on BSR
The RNA-seq result was successfully submitted to the SRA database of NCBI (accession no. PRJNA661989). In total, 80.17 and 93.94 million clean reads were obtained from the LFN and HFN bulks. The total mapped rates to IWGSC RefSeq v1.0 for the LFN and HFN bulk were 84.16 and 80.29%, respectively (Supplementary Table 9). In total, 282 significant DEGs were screened, including 170 downregulated genes and 112 upregulated genes by LFN/HFN bulk ( Figure 4A). Among these genes, 10 genes were enriched in four pathways, including protein processing in endoplasmic reticulum, sphingolipid metabolism, starch and sucrose metabolism, and galactose metabolism based on KEGG enrichment analysis ( Figure 4B). Through GO enrichment, 130 of these genes were enriched in 18 biological processes, 7 cellular components, and 27 molecular functions ( Figure 4C). These DEGs comprised α-Amy-3 and α-amylase inhibitor (Figure 4A), which is consistent with the variation in α-amylase activity between HFN and LFN bulk ( Figure 4D). A total of 117,664 SNPs and 5,993 InDels were obtained between the HFN and LFN bulks (Supplementary Figure 4). Under the threshold of SNP/InDel index > 0.80, 27 candidate regions distributed on chromosomes 1B, 2A, 2B, 3A, 3B, 3D, 4A, 5B, 7A, and 7B, including 37 SNPs and InDels related to FN, were screened (Table 3 and Supplementary Figure 4). Most of the candidate regions were located on 3D (13), followed by 3A (9) ( Table 3 and Supplementary Figure 4). Compared to the QTL results, two candidate regions partially overlapped with QFN.cib-3A and QFN.cib-3D, respectively, suggesting that there were genes controlling FN in the two regions. They were located on chromosome 3A (656618001-660197000, 668527001-669767000) including three SNPs and 3D (565077001-566487000) including 1 SNP. By searching QFN.cib-3A and QFN.cib-3D, two genes annotated as Basic 7S globulin and Basic 7S globulin 2 annotated in QFN.cib-3D were detected in 282 significant DEGs. The gene expression of Basic 7S globulin and Basic 7S globulin 2 was significantly higher (log 2 FC LFN/HFN = 1.39) and lower (log 2 FC LFN/HFN = −3.57) in the LFN than HFN bulk ( Figure 4A). They are probably candidate genes related to FN.

High-Density Genetic Map for MAS in Wheat Quality Breeding
Numerous QTL analyses have been conducted to clarify the genetic mechanism controlling complex traits of yield, quality, and resistance in wheat. However, because of the large genome size and limited genome sequence information in wheat, the one common problem is that the interval of QTL mapping is too large to further clone genes, especially for complex quality traits. Therefore, it is necessary to construct high-density genetic maps for gene identification and marker development in MAS for wheat breeding. Recently, high-density linkage map construction was facilitated by high-throughput SNP genotyping such as 55K, 90K, and 660K SNP arrays (Cui et al., 2017;Liu J. et al., 2017; Liu et al., 2018;Xu et al., 2019;Guo et al., 2020;Xin et al., 2020;Yang et al., 2020), thus increasing the accuracy and shortening the confidence interval of QTL mapping even to 0.09 cM/marker compared to those based on simple sequence repeats (SSRs) or DArTs (Jin et al., 2016).With the rapid development of nextgeneration high-throughput DNA sequencing, some procedures including restriction site-associated DNA tag sequencing (RADseq) and genotyping-by-sequencing are used in sequence-based marker development for genetic map construction, which have some advantages such as abundance, uniform distribution, and cost-effectiveness (Ganal et al., 2009;Liu H. et al., 2017). SLAF-seq based on RAD-seq for SNP discovery and genotyping in large populations has been widely applied for highdensity genetic map construction, QTL analysis, and gene cloning of extensive species such as rice (Li D. et al., 2016;Li et al., 2018;Quan et al., 2018;Zhu et al., 2018), soybean (Chen et al., 2016;Cao et al., 2017;Li et al., 2017;Kong et al., 2018;Pei et al., 2018;Zhang et al., 2018), cotton Palanga et al., 2017;Ali et al., 2018;Keerio et al., 2018), and peanut Wang et al., 2018). Some studies have investigated chromosomal localization (Li G. et al., 2016; and gene cloning for agronomic traits and resistance based on SLAF-seq of wheat (Hu et al., 2016;Yin J. L. et al., 2018). In our study, high-density genetic map construction and QTL analysis based on SLAF-seq were first applied for quality traits of wheat. Compared with most maps constructed previously by dozens or hundreds of markers such as SSRs or DArTs, the map in this study contained up to 12,674 SLAF markers; the average density of markers increased to 0.23 cM, and the region of gap > 5 cM in each chromosome was almost lower than 5%. Based on this genetic map, we identified a QTL QGH.cib-5D on the short arm of chromosome 5D, which is near the major QTL Ha loci for GH detected in previous studies (Mattern, 1973;Law et al., 1978). This suggests that the genetic map constructed by SLAFs is effective and reliable. In particular, genome sequencing of wheat has already been completed, which makes this genetic map more beneficial for discovering candidate genes for wheat quality.

Starch Pasting Properties Were Controlled by Many Loci With Minor Effects and Greatly Influenced by the Environment
Starch pasting properties predicted by RVA parameters were significantly associated with Asian noodle quality Wang et al., 2017). However, fewer QTL analyses for starch pasting properties have been conducted compared with those for other quality traits, and hardly any genes have been cloned, partly due to the complexity of the trait. In this study, 16 QTLs for starch pasting properties were detected, but no one was stable in at least two environments, even if SB (CV = 19.49-23.93%) and BD (CV = 26.94-37.08%) presented large variations among parents and RILs. This is consistent with the results of some previous studies that showed QTLs for pasting parameters were not stable in different environments (Mattern, 1973;Deng et al., 2015).
It also suggested that starch pasting properties were complex, controlled by many loci with minor effects, and greatly influenced by the environment. Although previous studies indicated that the Wx gene plays an important role in starch properties (Law et al., 1978;Giroux and Morris, 1998;Deng et al., 2015), no QTLs around the Wx gene were detected in this study. In addition, some QTLs including Ha loci on chromosome 5D identified for GH, and QTLs for gluten content correlated with Glu-B3 gene on chromosome 1B, were also detected for RVA parameters such as PV, FV, and SB . These results suggest that starch pasting properties are associated with other grain quality traits and are regulated by complex genetic factors. For some QTLs explaining higher phenotypic variance, such as QBD.cib-1A on chromosome 1A with PVE of 17.27%-17.38%, further studies are required in more environments.

Grain Protein Content
Quantitative trait locis for GPC were detected on all wheat chromosomes (Kumar et al., 2018). In this study, only one stable QTL, QGPC.cib-4A on chromosome 4A, was identified across all three environments (Table 2 and Figure 2). Previous studies also found GPC QTLs on chromosome 4A (Groos et al., 2003;León et al., 2006;Kunert et al., 2007;Peleg et al., 2009;Raman et al., 2009;Blanco et al., 2012;Chen et al., 2019), but they were far away from QGPC.cib-4A (Figure 5A). Considering the low density of makers in this area, further verification is required to confirm whether QGPC.cib-4A is a new QTL. The QGPC.cib-4A interval contained no genes, but closely adjacent to this loci, there were five consecutive genes classified into NAC domain-containing and no apical meristem (NAM) protein. NAC domain-containing proteins are a class of plant-specific transcription factors that play important roles in the growth and stress response. Through RNA interference, Uauy et al. (2006) showed that an NAC transcription factor, classified into NAM protein and encoded by TtNAM-B1 gene, accelerates senescence and increases nutrient remobilization from leaves to developing grains, resulting in an increase in grain protein, Zn, and Fe content by more than 30% in hexaploid wheat. This suggested that these five NAM genes closely adjacent to QGPC.cib-4A have effects on GPC.
In addition, QGPC.cib-4A is near the gene encoding the coatomer subunitβ'-1 (COPI-β'), which is primarily responsible for intra-Golgi and Golgi to endoplasmic reticulum transport of proteins (International Rice Genome Sequencing Project, 2005). To date, there have been no reports about the relationship between COPI-β'and GPC, but the efficiency of protein transport may affect GPC.
In this study, a major and stable QTL QGH.cib-5D was also identified at the end of chromosome 5DS with a PVE of 33.31%-47.99% in all three environments (Table 2 and Figure 2). Pina-D1 and Pinb-D1 were located within this QTL detected by BLUP ( Figure 5B). To analyze the relationship between QGH.cib-5D, Pina-D1, and Pinb-D1, Pina-D1, and Pinb-D1 were amplified by the marker designed by Chen et al. (2013) and compared between Chuanmai 39 and Chuanmai 42. There was no difference in Pinb-D1; Pina-D1 was missing in Chuanmai 39, which exhibited the hard phenotype. This suggested that the variance in GH was mainly due to Pina-D1. Pina-D1 was further amplified and compared to the polymorphism of flanking markers (M85128 and M85140) of QGH.cib-5D in 147 random lines of RILs derived from Chuanmai 39 and Chuanmai 42. Specifically, if QGH.cib-5D was controlled by Pina-D1, null Pina-D1 should correspond to "a" genotype of M85128 and M85140, the same as hard graintype Chuanmai 39, whereas the presence of Pina-D1 should correspond to the "b" genotype of M85128 and M85140, the same as soft grain-type Chuanmai 42. All these consistent types are shown by a transparent circle in Supplementary Figure 5A. However, the polymorphism of Pina-D1 was not consistent with that detected by the flanking markers of QGH.cib-5D for 36 lines; of them, 25 lines shown in green circles did not contain Pina-D1 but corresponded to the "b" genotype. Eleven lines shown in red circles contain Pina-D1 but corresponded to the "a" genotype (Supplementary Figure 5A). M85128 and M85140 had advantages in screening hard grains (hardness index > 70) compared to Pina-D1 (Supplementary Figure 5A). Probably genes other than Pina-D1 affect GH.
There are some reports about the molecular mechanism of GH difference excluding the effect of Pina-D1 and Pinb-D1. The main difference between soft and hard wheat was in adhesive strength among starch granules and the protein matrix surrounding starch granules, and proteins in soft grains are easier to separate from starch granules. Previous studies indicated that there was no relationship (Miller et al., 1984;Pomeranz et al., 1985) or some positive relationship (r = 0.26-0.43) between GPC and GH (Groos et al., 2004). One QTL identified for GH was close to the GluA3 locus on chromosome 1A, which is a candidate for GPC (Arbelbide and Bernardo, 2006). However, in this study, no QTL was identified for GPC and GH simultaneously, and the correlation coefficients between GH and GPC was lower than 0.15 in the two environments (Supplementary Figure 2). This means that the difference in GH is not mainly due to GPC variance and other carbohydrates such as lipids (Morrison et al., 1989;Panozzo et al., 1993;Pasha et al., 2010) and pentosans (Hong et al., 1989;Bettge and Morris, 2000;Pasha et al., 2010).
In addition to semidwarfing genes Rht-B1 and Rht-D1 associated with plant height, some agronomic traits and disease resistance were found to have an effect on GH, which was reflected by major QTL on chromosome 4BS and 4DS close to the semidwarfing (gibberellin-insensitive) genes; Reduced height (Rht)-B1 and Rht-D1 explained up to 20 and 34% of GH variance of RILs derived from the soft × soft type (Wang et al., 2012). Pina-D1 and Pinb-D1 were also found to play a role in plant defense mechanisms against pathogens (Bhave and Morris, 2008;Gautier et al., 1994). This suggested that difference in GH was regulated by gibberellin and other plant hormones, and it was also related to plant defense.
In our study, the physical localization of QGH.cib-5D revealed 33 genes were in this region. Only 11 genes express in the grain based on gene expression data of 22 genes from WheatExp (Supplementary Figure 5B) 6 . The 11 genes include Ethylene receptor, Carotenoid cleavage dioxygenase participating abscisic acid (ABA) synthesis, Sugar phosphatase and transport genes, Fantastic four-like protein, Protein kinase, Trypsin family protein, and Protein yippee-like. NBS-LRR disease resistance protein-like gene was also found in QGH.cib-5D. These genes related to plant hormone and disease defense will also be analyzed in further studies.

Falling Number
Falling Number is a quick evaluation method for α-amylase activity in wheat grains. Extremely high α-amylase activity with lower FN in mature wheat grain often results in downgrading quality and a sharp fall in wheat prices. PHS and LMA are two major causes of elevated α-amylase activity, but they differ in their formation mechanisms. LMA refers to high-pI α-amylase mainly synthesized during latter grain development, which seems to have little effect on grain germination and quality (Wrigley, 2006;Mares and Mrva, 2014;Newberry et al., 2018). However, the mechanism of LMA synthesis is not clear, and most QTLs for LMA were distributed on the long arm of chromosome 7B and were located near α-Amy-2 (Mrva and Mares, 2001;Mrva et al., 2009;Barrero et al., 2013;Mares and Mrva, 2014;Mohler et al., 2014). α-Amy-1 on chromosome 6 and α-Amy-4 genes were also involved in LMA synthesis (Barrero et al., 2013;Mieog et al., 2017). In addition, Rht-B1 and Rht-D1 genes also seem to be related to LMA, because some LMA QTLs were near the Rht-B1 and Rht-D1 genes on chromosomes 4B and 4D, and lower LMA expression genotypes were usually semidwarf (Kunert et al., 2007;Kulwal et al., 2012;Mohler et al., 2014;Börner et al., 2017). This is probably because Rht and α-amylase activities were both regulated by gibberellin. In this study, no stable QTLs were identified in the above regions. However, during late grain development, it was found that the transcript level of α-Amy-3 on chromosomes 5B and 5D of the LFN bulk was relatively higher than that of the HFN bulk through BSR results. In addition, the expression of 12 α-amylase inhibitor genes was significantly lower in the LFN than in the HFN bulk. It seems that α-Amy-3 was involved in α-amylase activity difference among RILs, possibly resulting from an α-amylase inhibitor, although the role of α-Amy-3 during grain development process is not clear yet.
In our study, eight QTLs on seven chromosomes were identified for FN (Table 2, Figure 2, and Supplementary Table 7). This indicated that FN is a complex trait and is regulated by many genes. We paid particular attention to three major or stable QTLs, including QFN.cib-3A, QFN.cib-2B, and QFN.cib-3D. QFN.cib-2B on chromosome 2B was detected in two environments, which partially overlapped with previous QTL QPhs.cnl-2B.1, which was simultaneously identified for PHS, grain dormancy, and rate of germination with the closest marker WMC474 (Table 2 and Figure 5C) (Munkvold et al., 2009). Subsequently, some genes related to grain dormancy, including an ABA receptor and other genes involved in calcium signaling, were speculated to be candidates for QPhs.cnl-2B.1 through comparative genetic analysis of rice and Brachypodium (Somyong et al., 2014).
We also identified a stable QTL QFN.cib-3A spanning 59.90-65.16 cM on chromosome 3A with a PVE of 10.95-18.30% (Table 2 and Figure 2). It could be detected in all three environments, also by BLUP over three environments and BSR (Tables 2, 3). This QTL is close to the downstream region (about 153 kb) of the Viviparous-1 (Vp-1) gene (Figure 5D), which is a transcription factor that is positively correlated to grain dormancy and embryo sensitivity to ABA and inhibit α-Amy genes expression (Hoecker et al., 1995;Bailey et al., 1999;Yang et al., 2007). Therefore, the effect of QFN.cib-3A on FN may be attributed to the function of Vp-1 gene. In addition, 95 genes were located in this QTL; three genes including nitrate transporter NRT1-2, NRT1/PTR family protein 2.2, and cysteine proteinase inhibitor were potential candidate genes for FN.
QFN.cib-3D spanning 140.72-145.13 cM on chromosome 3D was detected in two environments and by BLUP, which also explained the relatively higher phenotypic variation (16.58-25.74%) (Table 2 and Figure 2). Its effect was also validated by RILs derived from Kechengmai 4 × Chuanmai 42 ( Figure 3B). QFN.cib-3D was located between the previous QTL for PHS and grain color, flanked by marker Xgwm314 and Xcfd9 (Groos et al., 2002). All these QTLs were located in another previous QTL for germination index, sprouting index and FN, flanking by Xwmc552 and Xwmc533 (Fofana et al., 2008; Figure 5E). This indicates the presence of important genes in this region controlling FN. A total of 149 genes were annotated in QFN.cib-3D. The Tamyb10 gene, which encodes an MYB transcription factor and is a candidate for R loci of grain color, grain dormancy, and PHS, was included in QFN.cib-3D (Himi et al., 2005(Himi et al., , 2011. Previous studies showed that redgrained wheat (R-A1b/R-B1b/R-D1b, dominant wild-type alleles) is usually more resistant to PHS than white-grained varieties (R-A1a/R-B1a/R-D1a, recessive mutant alleles) (Gfeller and Svejda, 1960;Debeaujon et al., 2000;Groos et al., 2002;Fofana et al., 2008;Seshu and Sorrells, 2008;Mares and Mrva, 2014) and that the R loci upregulate flavonoid biosynthesis involved in antioxidation, such as grain dormancy in cereal plants (Himi et al., 2005) or enhanced grain dormancy through increasing the sensitivity of embryos to ABA (Himi et al., 2002). Tamyb10 may be candidate gene for QFN.cib-3D, based on that the effect of QFN.cib-3D on HFN was derived from the red-grained Chuanmai 42 allele.
However, no significant difference was found in transcript levels of either Tamyb10 or Vp-1 through BSR result (data not shown). This suggested that other genes in QFN.cib-3A and QFN.cib-3D affected FN. The transcript levels of Basic 7S globulin and Basic 7S globulin 2 annotated in QFN.cib-3D were significantly higher and lower in the LFN than in the HFN bulk, respectively. In addition, eight other continuous genes that belong to Basic 7S globulin 2 were in this region. They encode proteins with xylanase inhibitor domain and play roles in inhibiting cell wall degradation and radicle extension during germination. Therefore, the series of Basic 7S globulin 2 genes could be related to α-amylase activity and could be candidates for QFN.cib-3D on chromosome 3D.

QTL Clusters for Quality of Wheat
In this study, the location of QFN.cib-4D detected for FN and QPV.cib-4D detected for PV is partially overlapping (Figure 2 and Supplementary Table 7), which corresponds to the phenomenon that the measurement of FN and starch pasting properties both go through starch swelling and pasting. It suggests that some genes involved in starch synthesis and metabolism on chromosome 4D also affect FN. However, these genes have a relatively small impact on FN (PVE = 2.85-3.45%) and PV (PVE = 6.39-6.84%) and are easily disturbed by environmental deviation (Supplementary Table 7).
In addition, QGH.cib-5D on chromosome 5D identified for GH with PVE of 33.31-47.99% is also detected for FN (QFN.cib-5D.1) and PV (QPV.cib-5D) with PVE of 14.31-22.93% and 6.78% (Figure 2, Table 2, and Supplementary Table 7). Although QFN.cib-5D.1 and QPV.cib-5D for FN and PV were detected only in single environment, it may be a major QTL that was also identified by BLUP (Supplementary Table 7). It is interesting that the expression of Pina-D1 and Pinb-D1 gene in LFN was relatively higher and lower compared to HFN bulk (data not shown). This indicated that variation of GH has large effect in FN in special environment.
A previous study also found that QGh.caas-5D for GH flanked by marker Xcfd18 and Ha was detected in PV with a PVE of 10.4%, and Pinb-D1b at the Ha locus explained a large portion of the phenotypic variances for PV, FV, and SB, especially for pasting temperature (71.5%) . Some studies have confirmed the obvious effects of GH on starch properties. Soft wheat has better starch properties than hard wheat because more intact starch granules obtained from soft grains during the mill process resulted in higher flour viscosity compared to hard grains (Rogers et al., 1993;Xu et al., 2005;León et al., 2006).
A relatively higher PVE of QFN.cib-5D.1 means GH also has a large effect on FN in a special environment. This was also reflected by the fact that GH and FN are both related to Rht gene regulation by gibberellin (Tan et al., 2010;Gooding et al., 2012;Wang et al., 2012;Mohler et al., 2014;Zhang et al., 2014); therefore, some genes involved in the gibberellin response are probably candidates for QFN.cib-5D.1. As mentioned above, there is a pleiotropic gene or genes in QFN.cib-5D.1 controlling FN, grain texture, starch pasting properties, and mixograph parameter (Fofana et al., 2008;Zhang et al., 2009;Tan et al., 2010;Börner et al., 2017).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
QL designed the experiment, constructed the high density genetic map, identified the QTLs and candidate genes and wrote the manuscript. ZP was involved in experimental design, instruction, and improve the manuscript. YG was involved in wheat quality determinations. TL, JL, ZZ and HZ were involved in QTLs identification. GD, HL and MY were involved in population construction and growing and improve the manuscript. All authors contributed to the article and approved the submitted version.