Consensus genomic regions for grain quality traits in wheat revealed by Meta‐QTL analysis and in silico transcriptome integration

Grain quality traits are the key factors that determine the economic value of wheat and are largely influenced by genetics and the environment. In this study, using a meta‐analysis of quantitative trait loci (QTLs) and a comprehensive in silico transcriptome assessment, we identified key genomic regions and putative candidate genes for the grain quality traits protein content, gluten content, and test weight. A total of 508 original QTLs were collected from 41 articles on QTL mapping for the three quality traits in wheat published from 2003 to 2021. When these original QTLs were projected onto a high‐density consensus map consisting of 14,548 markers, 313 QTLs resulted in the identification of 64 MQTLs distributed across 17 of the 21 chromosomes. Most of the meta‐QTLs (MQTLs) were distributed on sub‐genomes A and B. Compared with the original QTLs, the confidence interval (CI) of the MQTLs was smaller, with an average CI of 4.47 cM, while the projected QTLs CI was 11.13 cM (2.49‐fold lower). The corresponding physical length of the MQTL ranged from 0.45 to 239.01 Mb. Thirty‐one of these 64 MQTLs were validated in at least one genome‐wide association study. In addition, five of the 64 MQTLs were selected and designated as core MQTLs. The 211 quality‐related genes from rice were used to identify wheat homologs in MQTLs. In combination with transcriptional and omics analyses, 135 putative candidate genes were identified from 64 MQTL regions. The findings should contribute to a better understanding of the molecular genetic mechanisms underlying grain quality and the improvement of these traits in wheat breeding.


INTRODUCTION
Wheat (Triticum aestivum L.) is one of the most important food crops worldwide, providing about 20% of the calories, protein, and fiber for the global human diet (Asseng et al., 2011;Ling et al., 2013;Shiferaw et al., 2013;Yadav et al., 2019), and its economic value has also been determined by the intrinsic grain quality factors affecting the final products (Ammiraju et al., 2001). Therefore, improving quality traits has long been an important goal in wheat breeding Jin et al., 2016;Nelson et al., 2006). Due to the importance of protein and gluten content, test weight (TW) for bread production, nutritional value, and economic impact of end product quality, wheat quality breeding has also mainly focused on improving these essential components Ruan et al., 2021;Sun et al., 2009).
Grain quality is a complex quantitative trait that is largely influenced by genetic traits, ecological factors, and cultivation practices (Guo & Wang, 2006;Li et al., 2005;Quraishi et al., 2017;Souza et al., 2004;Zhang et al., 2004). Stem length, plant height, spikelet length, flowering date, and several yield traits have been reported to affect quality traits (Jiang et al., 2019). Number of spikelets, grains per spike, thousand grain weight, and grain yield are significantly and negatively correlated with some quality traits of wheat (Yang et al., 2020). Identification of important quantitative trait loci (QTLs) and further utilization of superior genes are essential for genetic improvement of grain quality in modern breeding programs. Over the past two decades, numerous studies using traditional biparental linkage mapping (Groos et al., 2003;Krishnappa et al., 2021;Ma et al., 2012) and genome-wide association studies (GWAS) (Fiedler et al., 2017;Gao et al., 2021;Lou et al., 2021;Sallam et al., 2020) have been conducted for grain quality traits in wheat. QTLs for protein content (Blanco et al., 2002;Huang et al., 2006;Prasad et al., 2003;Yang et al., 2013) and TW (Mccartney et al., 2005;Sun et al., 2009Sun et al., , 2010Yang et al., 2013) have been reported on all chromosomes, while QTLs for gluten content have been reported on all chromosomes except chromosome 3D (El-feki et al., 2013;Li et al., 2009;Ma et al., 2012). Most of these QTLs contribute nearly 50% to explaining phenotypic variations. However, due to differences in genetic background of populations, genetic map, marker types, and environmental conditions, many of the previously identified QTLs are of low accuracy and have no direct application value in wheat molecular breeding (Kumar et al., 2020;Liu, Mullan et al., 2020;Yang et al., 2007).
Meta-analysis is a method of integrating QTLs obtained in previous studies to perform more comprehensive and accurate experiments for improvement of complex quantitative traits (Goffinet et al., 2000;Hanocq et al., 2007;Sosnowski et al., 2012;Veyrieras et al., 2007). This method has been successfully and widely used in plant breeding (Chardon et al., 2004;

Core Ideas
• We collected 508 quantitative trait loci (QTLs) associated with quality traits in wheat and performed meta-QTL (MQTL) analysis. • We identified 64 MQTLs, 31 of which were validated by marker-trait association in genome-wide association studies, and five were selected as core MQTLs. • Through the comprehensive platform, 1219 important candidate genes were excavated. • Through transcriptomics and gene function analysis, 135 highly expressed genes were identified. • This study aims to integrate QTL for quality traits in wheat and identify important genomic regions and key genes. Wang et al., 2013), and has achieved good results in QTL integration of various quantitative traits in several crops, for example, yield-related traits of maize (Zea mays L.) , insect-resistant combinations in maize , yield-related traits in rice (Oryza sativa L.) (Bhatia et al., 2017), and quality traits of cotton (Gossypium hirsutum L.) (Adhikari et al., 2017). In addition, useful results have been reported from meta-QTL (MQTL) studies on several related traits in wheat, such as yield-related  traits, grain quality traits (Jernigan et al., 2017), root-related traits (Soriano et al., 2019), and disease resistance (Liu, Salsman et al., 2020;Soriano et al., 2015;Venske et al., 2019). However, most of these studies did not examine the candidate genes behind the MQTLs, because the reference sequence of the wheat genome was not available in the past. With the advent of the high-quality Chinese spring wheat reference genome reported by the International Wheat Genome Sequencing Consortium (International Wheat Genome Sequencing Consortium [IWGSC] et al., 2018) and a large amount of wheat transcriptome data on a user-friendly platform (Borrill et al., 2016;Ramírez-González et al., 2018), meta-analysis has recently been used to identify key candidate genes, opening unexpected opportunities for breeding. Many key genes cloned by homology-based methods have been confirmed in yield-related MQTL regions , such as TaVrn1 (Yan et al., 2003), TaVrn2 (Yan et al., 2004), TaVrn3 (Yan et al., 2006), TARHT-B1 (Rht1) (Peng et al., 1999), two copies of TaPpd (Beales et al., 2007), Tags-D1 (Zhang et al., 2014), TaCKX2 (Zhang et al., 2011), and TaTGW6 (Hanif et al., 2016). However, in wheat, only a few major quality genes have been isolated by map-based cloning, so understanding the molecular basis of QTL/genes controlling grain quality still lacks essential knowledge.
In this study, the major genomic regions and putative candidate genes for grain quality traits in wheat were identified through meta-analysis. The objective of this study was to integrate the genomic regions associated with previously published QTLs for quality traits to obtain MQTLs with smaller confidence intervals (CIs) and transcriptome evidence for candidate gene identification, which will help improve our understanding of the genetic basis controlling grain quality. These reliable QTLs and candidate genes can then be used in further breeding research to improve wheat quality.

Initial QTL projection and MQTL analysis
After initial QTL data collection, all individual QTLs were projected onto a reference map by using BioMercator v4.2.3 software (Sosnowski et al., 2012). The high-density reference map developed by Bilgrami et al. (2020) was used, which is an integration of genetic reference maps from two previously published genetic maps (Maccaferri et al., 2015;Soriano et al., 2019). This map contains 14,548 markers, including simple sequence repeat (SSR), diversity arrays tech-nology (DART), and single nucleotide polymorphism (SNP), and other types of markers, with a total length of 4813.72 cM and a marker density of 3.02/cM. The QTLs that could not be found on the consensus map were eliminated. On each chromosome, the MQTL were calculated using the following two methods. In the first method, the meta-analysis proposed by Goffinet et al. (2000) was applied when the number of original QTLs on a chromosome was less than 10. Based on this approach, the best MQTL model with the lowest Akaike information criterion (AIC) values for QTL integration and identification of consensus MQTL positions in BioMercator was selected. In the second method, if the number of QTLs in a chromosome was at least 10, a two-step approach proposed by Veyrieras et al. (2007) was used. In the first step, the collected QTLs on individual chromosomes were clustered using standard parameters. The number of potential MQTLs per chromosome were then estimated based on the following five selection criteria, including AIC, Akaike information criterion correction (AICc), Akaike information criterion 3 candidate models (AIC3), Bayesian information criterion, and approximate weight of evidence. A QTL model that had the lowest values of the selection criteria was considered the best optimal model for the next step of meta-analysis. In the second step, the 95% CI and the positions of each MQTL were determined according to the optimal model selected in the previous step. The QTLs were integrated so that the peak position of the original QTLs was in the MQTL CI (Daryani et al., 2022), whereas the QTLs without the minimum AIC values were to be discarded.

Delineation of MQTL physical positions and verification by GWAS
To determine the physical positions of all MQTLs, the flanking markers within the target MQTLs were first found through a manual search based on the 95% CIs and the location of each MQTL. Then, flanking target MQTL markers were searched using the "Marker Information" function of the Triticeae Wheat Multi-Omics Center (http://wheatomics. sdau.edu.cn) to find the physical positions. If the physical positions of the flanking marker were not found, the GrainGenes database (https://wheat.pw.usda.gov/GG3) or the DART database (https://www.diversityarrays.com) were used to search for flanking marker sequences. The most likely physical positions were identified using the BLASTN program based on the Chinese spring chromosome RefSeq V1.0 from the Chinese Triticeae Multi-omics Center.
Quality trait data from 13 GWAS published from 2017 to 2021 were collected and used to verify the accuracy of these MQTL regions. Phenotypic data from these studies were collected with population sizes ranging from 163 to 2038, including spring wheat populations, winter wheat populations, and mixed populations of spring and winter wheat. Similarly, the physical position of marker-trait associations (MTAs) in these studies was determined using the source papers or BLASTN of the flanking marker sequences (Gudi et al., 2022;Yang et al., 2021) and compared with the physical coordinates of MQTLs (chromosome-wise). A single MTA falling within the genomic region of one or more MQTLs was considered co-localized .

Search for putative genes in the MQTL regions
Putative genes in the MQTL region were identified by the following three methods: (1) To identify the major putative candidate genes in the MQTL region using the strategy of orthologous wheat-rice comparison, the basic information on all functionally verified putative genes related to rice grain quality was downloaded from China Rice Data Center (https:// www.ricedata.cn/). Based on IWGSC RefSeqv1.1 (http:// wheat.cau.edu.cn/TGT), the homologous genes in wheat were found using Triticeae-Gene Tribe (http://wheat.cau.edu.cn/ TGT). The putative genes located in the MQTL regions were considered as important candidate genes affecting grain quality . (2) To obtain the core MQTL, the selection criteria (Venske et al., 2019) were as follows: (a) at least two overlapping QTL projections must be present in the resulting MQTL; (b) physical distance of less than 20 Mb; and (c) genetic distance of less than 1.0cM for each MQTL interval. To this end, high-probability genes within each highly refined MQTL were listed and then designated as putative candidate genes using the Triticeae Multi-omics Center (http://wheatomics.sdau.edu.cn) based on genomic features annotated by IWGSC_v1.1_HC_gene. (3) For the remaining MQTLs with longer CI, only a genomic region of 2 Mb (1 Mb region on both sides of the MQTL peak) was analyzed for the occurrence of putative candidate genes . The physical position of the peak of the MQTLs was calculated according to the method proposed by Saini et al. (2022).
The "Genes" tool of the WheatGmap database (https:// www.wheatgmap.org/tools/gene/information) was used to obtain detailed information (locus ID information and functional description) about the gene models in the MQTL regions detected by the last two methods. In addition, known wheat genes associated with the quality trait were searched for, and their physical positions were determined by the Wheat Multiomics Center (http://wheatomics.sdau. edu.cn) using their available IDs from previous corresponding studies. The physical positions of these genes were then compared to the genomic MQTL regions to determine their co-localizations (Saini et al., 2021).

Expression analysis of putative genes in the MQTL regions
Functional analysis of all genes from the MQTL region was performed using gene ontology (GO) enrichment on the GENEDENOVO cloud platform (https://www.omicshare.com/tools/) Then, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (https://www.omicshare.com/tools/) was used to filter out some metabolic pathways of all genes, such as starch and sucrose metabolism, protein processing, and so on. Finally, all genes available from the MQTL regions were subjected to further expression analysis using an Expression Visualization and Integration Platform (expVIP, http://www.wheat-expression.com) (Borrill et al., 2016), which includes expression data from 18 tissues throughout the wheat growth period (Ramírez-González et al., 2018). The expression levels of the putative candidate genes were assessed using at least two transcripts per million (TPM) and displayed using the TBtools (https://github.com/CJ-Chen/TBtools/releases) of TPM based on the normalized scale method. The transcriptomic data (http://www.wheat-expression.com) in this study included five tissues at different growth stages, such as grain at 2 dpa (day post anthesis), 14 dpa, and 30 dpa stages; spike at two nodes detectable, flag leaf stage and anthesis; leaf (excluding flag leaf) at tillering stage, 2 dpa and sowing stage; stem at 1 cm spike, two nodes detectable, and anthesis stage; roots at seedling stage, three leaf stage, and flag leaf stage.

Grain quality trait characteristics obtained from QTL studies
Reported QTLs for grain quality in wheat were from 41 previous studies published between 2003 and 2021 (Table  S1). They were identified in 39 biparental mapping populations, including 32 RILs, seven DH, and two F 2 populations. In total, we found 508 original QTLs for grain quality of wheat. Among these 508 QTLs, 46 QTLs were for PC, 241 QTLs were for GPC, 35 QTLs were for FPC, 70 QTLs were for WGC, and 116 QTLs were for TW. These initial 508 QTLs were distributed across all 21 chromosomes belonging to seven homeologous groups (1-7) and three sub-genomes (A, B, and D), ranging from nine QTLs on 3D to 55 QTLs on 1B. The number of QTLs in the three sub-genomes also varied, with 164 QTLs (32.28%) in sub-genome A, 225 QTLs (44.29%) in sub-genome B, and 119 QTLs (23.43%) in subgenome D. Furthermore, the initial CI of these QTLs ranged The Plant Genome

F I G U R E 1 The information of QTL and MQTL used in this study. (a)The distribution of projected QTLs and MQTLs on chromosomes. (b)
The reduction degree of QTL confidence interval (CI, 95%) after meta-QTL analysis. The pink and blue bars represent the average CI length (cM) of MQTL and projected QTL on chromosome, respectively, and the broken line represents the reduction folds of the QTL CI length. MQTL, meta-QTL; QTL, quantitative trait loci. from 0.01 to 105.73 cM, with an average of 16.545 cM. Of the 508 QTL, 42.13% (214) had an initial CI of less than 10 cM, and 67.72% (344) had an initial CI of less than 20 cM. Accordingly, individual PVE ranged from 0.07% to 46.27%, with an average of 9.59%. Only 33.86% (172) of loci had PVE values greater than 10%, indicating that most loci were secondary QTL (Table S2). Of the 508 QTLs initially identified, only 313 (61.61%) loci were successfully mapped to the QTL consensus map ( Figure S1). The remaining 195 QTLs were eliminated, due to the unavailability of common markers between the consensus and original linkage maps and the large CI of the original QTL. The 313 QTLs were distributed on all 21 chromosomes except 4D, 6B, and 7B, and the number of QTLs on each chromosome was also significantly different. More than 25 QTLs were found on chromosomes 2B (29), 3B (30), and 7A (26), whereas fewer than 10 QTLs were found on chromosomes 1D (9), 3D (4), and 5D (8). These results showed that QTLs related to grain quality traits in wheat were more distributed on A (142/313, 45.37%) and B (112/313, 35.78%) sub-genomes, but less distributed on D (59/313, 18.85%) sub-genome ( Figure 1a). The initial CI of these QTLs ranged from 0.01 to 44.51 cM, with an average of 11.13 cM. Of the 313 QTLs, 54.63% (171) had an initial CI of less than 10 cM, and 82.11% (257) had an initial CI of less than 20 cM ( Figure 1b). Accordingly, individual PVE ranged from 1.00% to 33%, with an average of 9.6%.

Identification of MQTL for grain quality traits
Of the 313 projected QTLs, 64 MQTLs were finally identified based on the minimum model value and at least two overlapping original QTLs ( Figure 2). Most MQTLs were associated with at least two different traits, reflecting the multi-trait effects of grain quality formation. Among the 64 MQTLs, 10, 55, and 12 MQTLs contained QTLs for PC, GPC, and FPC, respectively; 22 and 24 MQTLs contained QTLs for WGC and TW, respectively; and 16 and nine MQTLs contained QTLs for GPC, WGC and GPC, FPC, respectively ( Table 1). The 95% CI of the obtained MQTL ranged from 0.07 to 25.72 cM, with an average of 4.47 cM, which was 2.49 times less than that of the projected QTLs (11.13 cM) (Figure 1b), indicating that the mapping of these MQTLs was more accurate. The distribution of MQTLs on chromosomes was also different. For example, there were more than five MQTLs on chromosomes 2B, 4A, and 5B, whereas only one or two MQTLs were mapped on chromosomes 3D, 6D, and 7D ( Figure 1a). The physical locations of these 64 MQTLs were determined by querying the multi-omics database for wheat, and their physical intervals ranged from 0.45 to 239.01 Mb. Five core MQTLs, such as MQTL1A.1, MQTL1A.2, MQTL2A. 5, MQTL2D. 2, and MQTL5B.6 were positioned with a closer genetic distance of less than 1.0 cM, physical intervals of less than 20 Mb, and at least two overlapping QTL projections (Table 1). Thus, the selection criteria for further search for putative candidate genes were met.

Verification of MQTL by previously published GWAS
The GWAS results on grain quality traits of wheat published in recent years were used to further verify the authenticity of the above MQTLs (Table S3). Considering the relatively long linkage disequilibrium decay distance of wheat (∼5 Mb), the MTAs obtained from GWAS near MQTLs in the 5 Mb physical region were associated with MQTLs. Therefore, by comparing the physical coordinates of MQTLs with the physical coordinates of MTAs for each trait, 31 of the 64 MQTLs (48.44%) were co-located in at least one MTA (Figure 3). Some MQTLs were co-located with multiple MTAs. Among them, 12 MQTLs that co-localized with at least two MTAs, including MQTL1B.3 with 14 MTAs, followed by MQTL1A.1 and MQTL7A.2 with five MTAs. In addition, the physical coordinates of MQTL2D.4 (16.5-49.9 Mb) was also overlapped with the physical coordinates of MQTTL2D.3 (16.3-28.1Mb) identified in previous study .

Search for candidate genes within MQTL regions
Several strategies were used to identify candidate genes affecting wheat grain quality. First, the cloned genes controlling rice quality traits were searched in detail and 211 functional genes were identified. Then, their orthologous genes affecting quality traits in wheat were identified. Among them, there were 65 genes in 34 MQTL regions (Table S4), which were previously reported to regulate rice quality traits through multiple pathways, such as protein disulfide isomerase, starch synthase I, fragrance genes, seed disorder, amino acid transporter, histidine transporter, lysine-histidine transporter and others. The orthologous wheat genes related to these rice genes also were identified by BLASTP. Of these, 79 candidate genes were found in 34 MQTL regions (with an average of two per MQTL). Several candidate genes in the MQTL regions had similar effects on quality-related traits in wheat and rice. For instance, TraesCS4A02G214200 (MQTL4A.1) and PDIL1 affected PC, TraesCS7A02G120300 (MQTL7A.2) and OsSSI affected starch content, and TraesCS7D02G309400 (MQTL7D.1) and OsHT affected amino acid content. This result also indicated that the functions of the putative candidate genes are relatively conserved in rice and wheat. Moreover, these putative candidate genes have a high degree of confidence because their orthologous genes have been thoroughly investigated for their functions in rice quality traits. Second, 516 putative candidate genes were found among the five core MQTLs (MQTL1A.1, MQTL1A.2, MQTL2A. 5, MQTL2D. 2, and MQTL5B.6) (Table S5), and 625 candidate genes were identified by calculating the peak value of the remaining 25 MQTLs with longer CI (Table  S6). A total of 1217 putative candidate genes (excluding duplicate genes) were identified in 64 MQTL regions, with an average of 19.02 candidate genes per region. Among these putative candidate genes, 66 were genes for proteins with F-box family domains (such as TraesCS3D02G465200, TraesCS3D02G464100, and TraesCS3D02G465400 on MQTL3D.1), 27 were genes for cytochrome P450 proteins T A B L E 2 Summary of 135 putative candidate genes with significant expression (TPM > 2) within MQTLs TraesCS5B02G529300 Shaggy-related protein kinase delta Os03g0841800
Lastly, GO and KEGG pathway enrichment analysis were used to obtain the functional classification of these genes, and they were mainly associated with biological processes (20 subfunctions), molecular functions (13 subfunctions), and cellular components (13 subfunctions) (Figure 4). The most enriched GO terms related to biological processes were metabolic processes (433/1217, 35.58%) and cellular processes (391/1217, 32.13%). The most enriched GO terms related to molecular function involved binding (588/1217, 48.32%) and catalytic activity (429/1217, 35.25%). As for the cellular component, the genes were mainly enriched in the cell (197/1217, 16.19%) and in the cell part (196/1217, 16.11%). In the KEGG metabolic pathways, some of these genes were strongly involved in quality-related metabolic pathways, such as starch and sucrose metabolism, protein export, and citrate cycle (TCA cycle) ( Figure 5).

In silico transcriptome assessments of candidate genes
In silico expression analysis showed that 135 putative candidate genes were highly and specifically expressed in grain, spike, leaf, and root (TPM > 2) ( Figure 6, Table 2). The expression patterns of these putative candidate genes could be further divided into three classes. The putative candidate genes of class I were mainly expressed in the grain at the 14 dpa and 30 dpa stages. The putative candidate genes of class II were mainly expressed in the leaf at the 2 dpa stage and in the roots at the three leaf and flag leaf stages. The putative class III candidate genes were strongly expressed in the spike at the flag leaf stage and in two nodes detectable stages. Moreover, the expression of some genes such as TraesCS5B02G216000, TraesCS7A02G145200, and TraesCS4A02G267000 varied in the same tissue at different growth stages.

DISCUSSION
The discovery of molecular markers and advances in QTL mapping strategies have facilitated the identification of many QTLs controlling quality traits in wheat. However, many of these QTLs have low PVE values and large CI, so they were not suitable for marker-assisted plant breeding . In the present study, we performed a meta-analysis of 508 initial QTLs associated with three quality traits from 41 independent studies and found that these 508 initial QTLs were unevenly distributed across three genomes, with most QTLs on genomes A and B and the fewest on the D genome. The 95% CI of 64 MQTLs obtained in this study was 3.70 times lower than that of the original QTL studies. In a similar study, 141 MQTLs were detected with an average CI of 1.4 Mb, which was 8.8 times less than the CI (>12.1cM) of the original QTLs . Based on the released reference sequence of the Chinese Spring wheat genome, the 64 MQTLs detected in this study have unique physical positions, with physical intervals ranging from 0.45 to 239.01 Mb. Moreover, the distribution of MQTL was unequal across the different chromosomes, with the A and B genome carrying significantly more MQTLs than the D genome. The results were consistent with previous MQTL analyses of grain quality, yield, and yield-related traits Wang et al., 2015;Yang et al., 2021). It is possible that there are fewer QTLs for our target traits on the D sub-genome. However, it is also possible that QTL detection is affected by the lower level of DNA polymorphisms present on the D sub-genome. Because the D sub-genome is a recent evolutionary addition to the hexaploid wheat genome, gene flow from Aegilops tauschii to cultivated wheat is limited, resulting in relatively low genetic variation (Kumar et al., 2012).
To validate the MQTL results based on previously GWAS reports, we found that 48.44% (31/64) of the 64 MQTLs identified in this study were co-localized with MTAs identified by GWAS, which is consistent with previous studies (Aduragbemi et al., 2022;Yang et al., 2021) comparing MQTLs with MTAs. When MQTL were verified by more than one GWAS-MTAs study and included many original QTLs from independent experiments, the MQTL examined in our study may be more stable and accurate. In this case, nine MQTLs were detected at least twice in 13 GWAS reports, of which MQTL1A.1, MQTL4A.1, MQTL5A.3, and MQTL7A.2 were detected three times. This further confirmed the reliability of the present MQTLs. The identification of these MQTLs provided a basis for the accurate search of quality-affecting candidate genes that can facilitate molecular cloning and functional identification of genes related to grain quality in wheat.
Among the 64 MQTLs found in this study, five MQTLs were designated as core MQTLs. The interval of these core MQTLs had a relatively small genetic distance (n < 1.0cM), a smaller physical distance (n < 20Mb), and contained more original QTLs (n > 2) (Venske et al., 2019). These five core MQTL had a smaller CI (0.26cM), which was 17.19-fold smaller than all MQTL (4.47cM) in this study. Moreover, four of the five core MQTLs were controlling two traits, for example, MQTL1A.1 (GPC, WGC), MQTL1A.2 (GPC, FPC), MQTL2A.5 (GPC, WGC), and MQTL5B.6 (PC, GPC). In addition, studies show that GPC and WGC, and GPC and FPC are highly positively correlated (p < 0.01) (Huang et al., 2006;Liu et al., 2017;Raman et al., 2009) with correlation coefficients of 0.74 and 0.92, respectively. Three of the five core MQTLs contained initial QTLs detected in quite diverse and varied segregating populations, such as MQTL5B.6, which consisted of four initial QTLs from three different populations, with an average PVE of 20.02% (Gonzalez-Hernandez et al., 2004;Krishnappa et al., 2021;Yang et al., 2013). MQTL2A.5 was composed of three original QTLs from two different populations, with an average PVE of 7.84% (Cui et al., 2016;Peleg et al., 2009). MQTL1A.2 was composed of two initial QTLs from two different populations, with an average PVE of 6.58% (Raman et al., 2009;Sun et al., 2016). Interestingly, all five core MQTLs controlled for the F I G U R E 3 Validation of MQTL by MTAs for three quality traits in wheat from 13 different GWAS reports. The first circle of rectangles represents the initial QTL within MQTL, with the color from red to blue indicating the number of initial QTL from high to low. The remaining rectangles represent the colocation between MQTL and MTA, with the color from red to blue indicating the number of colocations from high to low. GWAS, genome-wide association study; MTA, marker-trait association; MQTL, meta-QTL; QTL, quantitative trait loci. same trait (GPC). The physical positions of the genes identified by previous GWAS, and meta-analyses overlapped with the MQTL positions in the current study (Gudi et al., 2022;Yang et al., 2021). This also proves the reliability of the currently determined MQTLs to be used for map-based cloning of quality traits in further studies.
Previous studies have shown that polyphenol oxidase (PPO) is the major cause of time-dependent discoloration in raw wheat dough (Beecher et al., 2011). Its homologs PPO-A1, PP0-A2, PPO-D1, and PPO-D2 were strongly expressed in grains and their physical location was found within the regions of MQTL2A.2 and MQTL2D.1 in the current study. The two genes Traescs2D02G033900 and Traescs2D02G468600 encoding PPO, were also identified in the interval of our MQTL2D.2 and MQTL2D.1. In addition, the early study showed that OsmtSSB1-mediated mitochondrial function plays a critical role in determining subaleurone cell fate in rice , and its related gene TraesCS1A02G327700 was found in the region of MQTL5B.2 in this study.
As shown in the previous studies, MQTL regions correlate strongly with gene density in the genome (Kumar et al., 2020;Maccaferri et al., 2015;Saini et al., 2022). For F I G U R E 4 Level 2 gene ontology (GO) terms for all putative candidate genes from MQTL regions. MQTL, meta-QTL.

F I G U R E 5
Twenty of all KEGG enrichment pathways for all genes from MQTL regions. Enrichment factors are shown on the horizontal axis and paths are shown on the vertical axis. Different colors represent different p-value. From magenta to blue, the p-value range from large to small, and the enrichment degree is more and more significant. The size of the circle represents the number of genes enriched in this pathway. KEGG, Kyoto Encyclopedia of Genes and Genomes; MQTL, meta-QTL.
instance, they identified 15,772 gene models for yield, baking quality, and GPC (Quraishi et al., 2017) and 2533 gene models for rheological dough characteristics, nutritional properties, and processing quality (Gudi et al., 2022). In our study, 1217 genes were obtained in the interval of 64 MQTLs and selected for further analysis. Some of our putative candidate genes were strongly implicated in the pathways of starch and sucrose metabolism, protein export, and the TCA cycle ( Figure 5). The Plant Genome F I G U R E 6 Expression characteristics of 135 candidate genes in five tissues. From blue to magenta, the expression of the gene in the relevant tissue ranges from low to high. dpa, day post anthesis.
Starch and sucrose metabolisms play an important role in carbohydrate partitioning in the plants and have a major impact on wheat quality. For example, granule-bound starch synthase I, which controls amylose content , and key enzymes that control amylopectin content, namely soluble starch synthase (Li et al., 2003), de-branched enzyme, and branched enzyme (Crofts et al., 2017), have an impact on final quality by affecting processing quality.
In our study, the expVIP platform was used to identify 1217 putative candidate genes in the MQTL region. Only the putative candidate genes exhibiting more than two TPM were considered (Venske et al., 2019;Wagner et al., 2013), combined with some genes previously reported to affect grain quality traits in wheat. For instance, the GPC, Glu-1, Glu-3, and Gli genes, which encode NAD transcription factor, high and low molecular weight glutenin subunit (HMW-GS/LMW-GS), and gliadin, respectively, and are directly involved in grain protein formation (Krystkowiak et al., 2017;Plessis et al., 2013;Uauy et al., 2006). Grain protein is one of the three nutrients in food, which is of great significance to the color, aroma, and taste of food. We found a total of 135 putative genes with high and/or specific expression patterns in grains, ears, or leaves (TPM > 2) ( Figure 6). These putative candidate genes showed tissue and developmentdependent expression patterns that can strongly influence grain quality traits. For example, TraesCS1D02G317300, encoding the high molecular weight glutenin subunit, was specifically expressed in grain 14 dpa. The content of HMW-GS in wheat grains is low, accounting for only about 10% of the traditional storage protein but may explain about 60% of the variation in bread baking quality of wheat (He et al., 2005;Singh et al., 2016). TraesCS2A02G035700 and TraesCS2D02G034900, encoding the arginase family, were specifically expressed in cereals. In rice, the homologous gene OsARG is the key enzyme of arginine metabolism. Loss of arginase activity in the nglf-1 mutant resulted in accumulation of arginine in panicles and increased seed formation rate, grain yield, and number of full grains (Ma et al., 2013). Studies have demonstrated that histone deacetylases (HDACs) play a key role in plant growth and development, including flower development , seed development (Wu et al., 2000), root development (Xu et al., 2005), and cell proliferation (Nelissen et al., 2005) and cell death (Bourque et al., 2011;Huang et al., 2007) during organ growth (Dong et al., 2016). TraesCS2D02G075800 encodes NAD-dependent protein deacetylase and conserved HDACs family (Imai et al., 2000). These genes could be further investigated to reveal their possible role in wheat quality and their application in breeding programs. Although many genes related to wheat quality are rarely studied, their homologous genes in rice and Arabidopsis have been reported to potentially play key roles in plant growth and development. For example, some homologous genes in rice have been implicated in the regulation of rice quality traits, such as OsmtSSB1, WTG1, and OspPLAIIIα (Huang et al., 2017;Li et al., 2021;Liu et al., 2015). Other identified genes related to protein and sugar transport were also strongly expressed in the corresponding tissues, suggesting that these 135 putative candidate genes are most likely involved in the regulation of grain quality traits in wheat.

CONCLUSION
This study confirms the validity and reliability of the metaanalysis method for understanding the genetic mechanisms underlying grain quality traits in wheat. The candidate genes we identified all have a direct or indirect effect on grain development and ultimately grain quality. The critical genes were involved in three metabolic pathways involving starch and sucrose metabolism, protein export, and the TCA cycle. Some of the genes had functions as previously reported for grain development and quality formation. The discovered MQTL regions and candidate genes can be used in marker-assisted selection to improve grain quality in future breeding research programs.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare no conflicts of interest.

D A T A AVA I L A B I L I T Y S T A T E M E N T
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. Some of the data are shown in Supporting Information.