Mapping QTL for Phenological and Grain-Related Traits in a Mapping Population Derived from High-Zinc-Biofortified Wheat

Genomic regions governing days to heading (DH), days to maturity (DM), plant height (PH), thousand-kernel weight (TKW), and test weight (TW) were investigated in a set of 190 RILs derived from a cross between a widely cultivated wheat-variety, Kachu (DPW-621-50), and a high-zinc variety, Zinc-Shakti. The RIL population was genotyped using 909 DArTseq markers and phenotyped in three environments. The constructed genetic map had a total genetic length of 4665 cM, with an average marker density of 5.13 cM. A total of thirty-seven novel quantitative trait loci (QTL), including twelve for PH, six for DH, five for DM, eight for TKW and six for TW were identified. A set of 20 stable QTLs associated with the expression of DH, DM, PH, TKW, and TW were identified in two or more environments. Three novel pleiotropic genomic-regions harboring co-localized QTLs governing two or more traits were also identified. In silico analysis revealed that the DArTseq markers were located on important putative candidate genes such as MLO-like protein, Phytochrome, Zinc finger and RING-type, Cytochrome P450 and pentatricopeptide repeat, involved in the regulation of pollen maturity, the photoperiodic modulation of flowering-time, abiotic-stress tolerance, grain-filling duration, thousand-kernel weight, seed morphology, and plant growth and development. The identified novel QTLs, particularly stable and co-localized QTLs, will be validated to estimate their effects in different genetic backgrounds for subsequent use in marker-assisted selection (MAS).


Introduction
Micronutrient deficiencies affect two billion people globally. Annually, they accounts for 45% of all child death under the age of five [1]. Furthermore, 52 million children are emaciated and 155 million children are stunted due to micronutrient deficiencies [2]. The iron (Fe) and zinc (Zn) deficiencies are more pronounced and affect one-third of the population in developing nations, particularly children and pregnant women [3,4]. Various interventions have been proposed to combat micronutrient malnutrition, namely, pharmaceutical supplementation, industrial fortification, dietary diversification and biofortification [5,6]. Among them, genetic biofortification employing conventional, molecular or transgenic methods to produce nutrient-rich crop varieties is considered as a sustainable, cost-effective long-term strategy for addressing nutritional needs [7]. With this concern, various biofortified wheat varieties were developed and deployed for production in India, such as 'WB02 , 'PBW01Zn HUW711 and 'Zinc Shakti,' which have a 20-40% higher zinc content than local controls [8].
The adaptation of biofortified cultivars by farmers depends mainly on their yield potential. Grain yield is a complex trait and is an outcome of the combined effect of several yield-contributing agro-morphological and physiological traits [9,10]. Important yield component traits such as DH and DM play a key role in determining the crop's

Phenotyping and Data Analysis
All the genotypes of a mapping population were phenotyped for three agro-morphological (DM, DH, PH) and two grain-related traits (TKW and TW). The RILs, along with the parents, were grown in a randomized-complete-block design with two replications, under the condition of assured irrigation. The data on DH was measured by counting the number of days from germination to 50% of the plants heading in a plot. DM was measured by counting the number of days from germination to physiological maturity, when more than 50% of spikes were ripe and had turned yellow. PH was measured from the ground to the tip of the spike excluding awns, at the late grain-filling stage. TKW was measured with the Seed Count digital imaging system (model SC5000; Next Instruments Pty Ltd., Condell Park, NSW, Australia), which was standardized to measure TKW. The Seed Count system can rapidly and accurately measure wheat grain samples, determining the grain number and the grain physical characteristics based on flatbed scanner technology.
The Meta-R (Multi Environment Trial Analysis with R) version 6.0 software was used for phenotypic analysis (https://excellenceinbreeding.org/toolbox/tools/multi-environmenttrail-analysis-r-meta-r (accessed on 20 March 2022)). Best linear unbiased predictors (BLUPs) of each RIL were obtained for an individual year and across years, and used in QTL analysis. The broad-sense heritability (h 2 BS ), pairwise correlation co-efficient (r g ) between traits, and the coefficient of variance (CV) were also estimated using Meta-R. The statistical equations used are presented below.
BLUPs (individual environment): Y ik = µ + Repi+ Gen k + ε ik BLUPs (across environments): Y ijkl = µ + Env i + Repj (Env i ) + Gen l + Env i × Gen l + ε ijkl where, Y ik -trait of interest, µ-mean effect, Repi-effect of the ith replicate, Gen k -effect of the kth genotype, ε ijk -error associated with the ith replication and the kth genotype, Env i -effects of the ith environment, Env i × Gen l -environment by genotype interaction, σ 2 g -genotype variance, σ 2 e -error variance, nreps-number of replicates, σ 2 ge -genotypeby-environment-interaction variance component, nEnvs-number of environments in the analysis, and ASED-average standard error of the differences between pairs of means.

Genotyping, Linkage Map Construction, and QTL Analysis
The RILs were genotyped with 40,059 DArTseq markers [39]. The DArTseq markers were generated by following genotyping-by-sequencing (GBS) techniques, which employs a combination of restriction enzymes and reduce the genome complexity, to obtain a representation of the whole genome. The quality filtering of the FASTQ files was carried out using a Phred quality score of 30. More stringent filtering was also performed on the barcode sequences, using a Phred quality score of 10, which represents 99.9% of basecall accuracy for at least 75% of the bases. Allele calls for SNP were generated using a proprietary analytical pipeline developed by DArT P/L. The markers were quality filtered using the bin function in the QTL IciMappingv4.2 software (http://www.isbreeding.net) (accessed on 1 November 2022), to remove parental non-polymorphic markers, markers with >30% missing data, and the redundant markers. Both linkage-map construction and QTL mapping were carried out with the QTL IciMappingv4.2 software. The linkage map was constructed with the LOD threshold of 3.0 between adjacent markers [40], the markers were ordered with the nnTwoOpt algorithm and rippled using 5 cM window size. The QTLs were mapped by applying the LOD threshold of 2.5, following the inclusivecomposite-interval-mapping (ICIM) method.

In Silico Analysis of QTL and QTL Additive-Effects Estimation
The in silico analysis was carried out in the Ensemble Plants database (http://plants. ensembl.org/index.html (accessed on 1 November 2022)) of the bread-wheat genome with the basic local alignment search tool (BLAST). The sequence information of the markers falling within the confidence interval of the QTL was used to identify putative candidate genes within the vicinity of the markers.
The RILs were grouped based on genotypic classes of multiple stable and major-effect QTLs, using flanking markers, and they were compared with the actual phenotype traits across environments, to identify the QTL additive-effect and best combinations of QTLs for each trait.

Genetic Variability and Correlation Analysis
Individual environments and the pooled mean across environments, including phenotype ranges of RILs and their parents are presented as box plots in Figure 1. The lowest mean was recorded in Y1 for all three agro-morphological traits, whereas the highest mean was observed for both grain-related traits in Y1. The frequency distribution of the population is presented in Figure S1. A wide range of variation was observed among the RILs, and few transgressive segregants surpassing both the parents appeared for all the traits. The population exhibited a continuous and near-normal distribution for all the traits, suggesting a polygenic nature. The h 2 BS, CV and genotypic variance are provided in Table 1    The correlation of agro-morphological and grain-related traits is presented in Table 2. The correlation among DH, DM, and PH was highly significant and positive (p < 0.001). However, the correlation between TKW and all three agro-morphological traits was significant and highly negative except in Y3 for PH (p < 0.01-0.001). Furthermore, the correlation between TW and DH, and TW and DM were highly significant and negative in Y2 and for the pooled mean (p < 0.05-0.001). The correlation between TW and PH was significant and positive for the pooled mean (p < 0.01). The correlation between grain-related traits was significant and negative in Y2 (p < 0.01). We also observed no correlation or negative correlation for GZnC and GFeC with TW, DH, DM (p < 0.05-0.001), but positive correlation with TKW (p < 0.001). Interestingly, PH showed positive correlation with GFeC (p < 0.05-0.001) and exhibited no association with GZnC [32].

Linkage-Map Construction
The linkage map was constructed using 909 DArTseq markers. The parental nonpolymorphic markers, markers with <0.05 and >0.95 allele frequency, >30% missing data and >20% heterozygosity were discarded, to obtain high-quality informative markers. The linkage map spanned a genetic distance of 4665 cM, with an average marker density of 5.13 cM. The subgenome and chromosome-wise distribution of markers are given in Table 3. The highest number of markers were mapped on subgenome A (412) followed by subgenome B (370), and the lowest number of markers (127) were mapped on subgenome D. A detailed description of the framework linkage-map is provided in [32].  Table 4 and Figure 2. For DH, the six QTLs were detected on chromosomes 2B, 5A, 5B, and 7D, with the phenotypic variance explained (PVE) ranging from 3.20 to 37.98%. All the identified QTLs were stable, as they were identified in two environments (Y1 and Y2) along with the pooled mean, except for QDH-5A.3, which was identified in Y2 and the pooled mean. Two major QTLs, i.e., QDH-5A.1, located at 122 cM on 5A, and QDH-7D, located at 104 cm on 7D, recorded the highest PVE, ranging from 26.21 to 37.98% and 12.52 to 16.49%, respectively. All the identified QTLs had positive alleles from the Kachu parent except QDH-7D, which had positive alleles from the Zinc-Shakti parent.    For DM, five QTLs were mapped on chromosomes 2B, 5A, 5B, 6A, and 7D, with a PVE ranging from 3.33 to 38.23%. Three (QDM-2B, QDM-5A, and QDM-7D) out of five QTLs were consistently expressed in two environments, along with the pooled mean. A major and consistent QTL, i.e., QDM-5A, flanked between 5411867 and 1141498 was identified on chromosome 5A at 122 cM, with PVE ranging from 31.54 to 38.23%. Similarly, another major and consistent QTL, i.e., QDM-7D, flanked between 2249010 and 100024878 was mapped on chromosome 7D at 105-106 cM, with a PVE ranging from 11.65 to 11.81%. However, the third consistent QTL, i.e., QDM-2B, flanked between 3570063 and 3935335 on chromosome 2B and located at 289-290 cM was recorded with a low PVE, ranging from 3.39 to 6.29%. The remaining two QTLs, i.e., QDM-5B and QDM-6A were identified in one environment each on chromosomes 5B and 6A, with a PVE of 4.32% and 3.33%, respectively. All the QTLs of DM except QTL on chromosome 7D had positive alleles from Kachu, whereas QDM-7D had positive alleles from Zinc-Shakti.

Quantitative-Trait-Locus Additive-Effects
The additive-effects of the stable and pleiotropic QTLs were investigated for agromorphological and grain-related traits ( Table 5). The combination of two or more QTLs has an additive-effect for agro-morphological traits (PH, DM, and DH); however, there is no significant additive-effect for grain-related traits (TKW and TW). Three QTL (QDH-5A.1, QDH-5A.2, and QDH-7D) combinations had the highest average across the environments for DH, and this combination was found in six RILs. For DM, the highest average was observed in the QTL combination of QDM-2B, QDM-5A, and QDM-7D, and this combination was identified in seven RILs. The three QTL combinations including QPH-3A, QPH-3D, and QPH-5A yielded the highest average for PH.

Discussion
The essential micronutrients play a pivotal role in the proper functioning of metabolic pathways and physiological activities of both plants and humans. The micronutrient deficiency hinders both physical and mental development, and causes various health-related complications such as anemia, infection, infant mortality, diarrhea, stunting, wasting, etc. Efforts have been made to develop biofortified cultivars to address the micronutrient deficiency, and some have been released for commercial cultivation in India, namely, WB-02, PBW01Zn, HUW711 and Zinc-Shakti. However, the adaptation of these varieties by farmers hugely depends on their yield potential. Hence, we tried the genetic dissection of yield-contributing phenological and grain-related traits, using a population developed with high-zinc-biofortified wheat, "Zinc Shakti" and the widely adopted cultivar "Kachu". The phenotype-based selection in conventional wheat-breeding has improved wheat yield for several decades across the major wheat-growing regions; however, yield plateaus have been observed in many crops including wheat, and the genotype-based modern breeding tools may complement the conventional breeding approaches to break the yield plateaus. Recent efforts have led to the sequencing of the complex wheat genome, which could be valuable in the varietal improvement programs through the identification of marker systems, putative candidate genes, etc. Several QTLs/MTAs have been identified for yield and contributing traits; however, additional genomic studies using varied genetic materials may be useful, as the saturation point has not been reached [41]. Furthermore, phenotyping plays a critical role in the identification of QTLs, as many of the identified QTLs may be genetic-material-specific or environment-specific QTLs. Additionally, the wheat genome is highly complex, and there is always the possibility of detecting novel genomic regions for yield and component traits.

Mapping Population
The RILs were developed by crossing the high-Zn parent 'Zn-Shakti', derived from a synthetic hexaploid wheat (Pedigree: CROC1_/ AE.SQUARROSA(210) // INQALAB91*2 / KUKUNA /3 / PBW343*2 / KUKUNA) and the high yielding and highly adopted line in South Asia, 'Kachu', popularly known as 'DPW-621-50 (Pedigree: KAUZ//ALTAR-84/(AOS)AWNED-ONAS/3/MILAN/KAUZ/4/HUITES). The parents were distantly related (coefficient of parentage (COP = 0.17) and SNP-based diversity (π = 0.26)). This is highly likely to be the reason for obtaining stable QTLs for both nutritional [32] and yield-related traits in our study. The significant effect of environment and genotypeenvironment interactions (GEI) was observed in the expression of yield and component traits (Figure 1). The RIL population has been extensively tested consecutively for three years, as the magnitude of GEI is an important factor in the detection of environmentspecific QTL(s) as well as consistent QTL(s).

Correlation Studies
The highly positive and significant correlation among DH, DM, and PH suggests that some common genetic mechanism existed in the expression of these traits. The strong positive association of agro-morphological traits was further supported by the identification of a common genomic region flanked between 5411867 and 100024127, which harbors QTLs (QPH-5A, QDH-5A.1, and QDM-5A) for the same three agro-morphological traits (PH, DH, and DM). However, the correlation between TKW and all the three agro-morphological traits was significant and highly negative, except in Y3 for PH. Significant correlations found in this study have also been reported in earlier studies [23,42]. The positively associated traits (DH, DM, and PH) can be improved simultaneously; however, for the improvement of negatively associated traits, the adoption of suitable breeding-methods to break the undesirable linkages is envisaged. In our previous studies [32,38], we observed no correlation or a negative correlation for GZnC and GFeC with TW, DH, and DM, but a positive correlation with TKW. PH had a positive correlation with GFeC, and exhibited no association with GZnC. Interestingly, we found co-localized QTLs for GFeC with DM, TKW and TW, and for GZnC with DH, DM, PH, and TKW.

Linkage Map
The linkage map was constructed with 909 DArTseq high-quality informative markers in a set of 190 RILs spanning a map length of 4665 cM, with an average marker density of 5.13 cM, and was utilized for the QTL analysis. A total of thirty-seven QTLs, including twelve for PH, six for DH, five for DM, eight for TKW, and six for TW were identified. The highest number of QTLs were identified for subgenome A (18 QTLs), followed by subgenome B (12 QTLs) and subgenome D (7 QTLs). The subgenome-wise distribution of QTLs is similar to the marker distribution pattern of subgenomes, as the highest number of markers were mapped on subgenome A (412), followed by subgenome B (370), and the lowest number of markers (127) were mapped on subgenome D. The type and size of a mapping population and frequency and distribution of markers on the linkage map are considered to be key determinants for QTL identification. In previous reports also, the pattern of QTL distribution for grain-quality traits was similar [43]; however, the enrichment of subgenome D with additional markers substantially increased the power of QTL identification [44]. Furthermore, a total of fifteen stable and three pleiotropic QTLs were identified.

Pleiotropic QTL Regions
Three common genomic regions flanked between 5411867 and 100024127 on 5A, 1206846 and 1082014 on 6A, and 1698406 and 100027274 on 6A, harboring co-localized QTLs associated with two or more traits, were identified. The significant associations of traits in the desired direction are always beneficial for the simultaneous improvement of the associated traits. One of the parents of the mapping population is a high-zinc parent (Zinc-Shakti) and the other parent is a high yielding, widely adopted parent (Kachu). Therefore, there is always a possibility of identifying novel QTLs for both grain quality and agro-morphological traits. In our earlier study, several QTLs were reported for grain iron-and zinc-concentration [32]. A total of eight genomic regions identified for various agro-morphological and grain-related traits in the present study correspond to our previous study [32] for grain iron-and zinc-concentration.  (QFeC-6A). Therefore, effective utilization of these stable and co-localized QTLs in varietal-improvement programs through MAS may simultaneously improve the grain yield and grain micronutrient-concentration in wheat.

QTL Utilization Strategy
The exploitation of the QTLs of agro-morphological traits depends on the breeding objective to develop lines suitable for different growing environments. The correlation studies indicated that the best ideotype would be reduced DH and DM, and increased PH, for enhancing grain-micronutrient (GFeC and GZnC) and yield (TKW) traits. Therefore, combining PH QTLs QPH-3A and QPH-3D from Zinc Shakti and QPH-5A from Kachu, and individual QTLs, QDH-7D and QDM-7D from Zinc Shakti, would be of great value for attaining the objective. Combining QTLs caused only a slight increase in TW, while TKW further stresses the multigenic nature and breeding complexity of these yield-determining traits. The pleiotropic regions found on 6A between 106.5 and 112.5 cM for GFeC, DM, TKW and TW indicate the role of this chromosome in governing multiple complex traits.

Putative Candidate Genes
The putative candidate genes underlying stable and pleiotropic QTLs for agro-morphological and grain-related traits were identified through a BLAST search ( Table 6). The QTLs identified on various chromosomes were located in gene coding-regions related to proteinbinding factors, photoreceptors, transporters, etc. For example, QDH-2B and QDM-2B associated with days to heading and maturity encode MLO-like protein (TraesCS2B02G097900) genes, which were found to have a functional role in the regulation of agro-morphological traits in rice. For instance, OsMLO12, OsMLO4, OsMLO10, and OsMLO8 showed preferential expression in mature pollen, in the root tips, throughout the roots except at the tips, and in the leaves and trinucleate pollen, respectively [59]. Furthermore, in osdxr mutants that exhibited defects in the light response, OsMLO1, OsMLO3, OsMLO8, and four calmodulin genes were down-regulated.
The three QTLs (QDH-5A.1, QDM-5A and QPH-5A) associated with days to heading, maturity, and plant height, encode Phytochrome (TraesCS5A02G391300). Phytochromes play an important role in light-signaling and the photoperiodic control of flowering-time in plants. A candidate gene, i.e., HvPHYTOCHROME C, was identified for the early-maturity 5 locus modulating the circadian clock and photoperiodic flowering in barley [60]. Similarly, the ectopic expression of a phytochrome B gene from Chinese cabbage in Arabidopsis promotes seedling de-etiolation, dwarfing in mature plants, and delayed flowering [61].
Plant height is one of the important agro-morphological traits, and significantly contributed to yield gains during the Green Revolution. Putative candidate genes such as Zinc finger, RING-type, and Cytochrome P450 regulate plant growth in many crop plants. One QTL, i.e., QPH-3D, associated with plant-height, encodes Zinc finger, and RING-type (TraesCS3D02G095600), found to have a role in flowering-time modulation, abiotic-stress tolerance and plant growth and development. The characterization of the zinc finger transcription-factor gene (Cm-BBX24) suggests the dual role of flowering-time modulation and abiotic-stress tolerance in the chrysanthemum, partly by the influence of GA biosynthesis [62]. The other zinc finger transcription-factor gene, i.e., TANDEM ZINC-FINGER PLUS3 (TZP) controls the expression of growth-promoting transcriptional regulators via the direct association with light-regulated promoter elements in Arabidopsis. Elucidating how such novel protein-complexes coordinate gene expression may help breeders to optimize plant growth and development [63]. The Really Interesting New Gene (RING) finger proteins characterized by the RING domain, which contains around 40-60 residues, are thought to be E3 ubiquitin ligases. The RING-finger proteins play significant roles in plant growth, stress resistance, and signal transduction in crop plants [64]. The RING-type E3-ubiqitin-ligase barley genes (HvYrg1-2) control the characteristics of both vegetative organs (leaf width and weight, plant weight and height, flowering-time, grain-filling duration, root-growth) and seed components (thousand-grain weight and seed morphology) in barley [65]. A novel RING finger protein, i.e., PLANT ARCHITECTURE and GRAIN NUMBER 1 (PAGN1) regulates the number of grains per panicle, the plant height, and the number of tillers in rice [66].
The same QTL, i.e., QPH-3D, which encodes Zinc finger, RING-type also encodes Cytochrome P450 (TraesCS3D02G077400). The cytochrome P450 gene, i.e., CsCYP85A1, is a putative candidate for super compact-1 (scp-1), which controls the plant height in cucumber [67]. Another cytochrome P450 family member, OsCYP96B4, reduces plant height in a transcript dosage-dependent manner in rice [68]. Amino-acid transporters play a key role in resource-allocation processes that support plant growth and development; an amino-acid transporter, i.e., LHT1, plays a role in plant-growth in rice, and its disruption leads to growth inhibition and low yields [69]. The grain-weight QTL QTKW-3A.1 encodes pentatricopeptide repeat (TraesCS3A02G232800), which affect grain-weight in rice and maize. The novel QTL (qKW9) related to kernel-size encodes a pentatricopeptide repeat protein that affects photosynthesis and grain-filling in maize [70]. The peptide transporter OsNPF7.3 enhances nitrogen allocation and increases grain-yield in rice through enhanced tiller numbers, panicles per plant, filled-grain numbers per panicle, and grain nitrogen-content [71].

Conclusions
The study with 190 RILs has shown that DH, DM, PH, TKW, and TW were quantitatively inherited traits. The strong positive correlation among DH, DM, and PH suggested the possibility of improving these traits simultaneously. Furthermore, biofortification traits such as GFeC and GZnC could also be increased simultaneously, as they exhibited a significant and positive correlation with TKW; this may be due to the elimination of the dilution effect in the Zinc-Shakti parent, through careful selection for high GFeC and GZnC, along with high TKW. A total of thirty-seven QTLs, including twelve for PH, six for DH, five for DM, eight for TKW, and six for TW were identified. A total of 20 stable QTLs were identified in more than one environment and associated with the expression of DH, DM, PH, and TKW. Similarly, three novel pleiotropic genomic-regions harboring co-localized QTLs governing two or more traits were also identified. Three QTLs (QPH-4A, QPH-3B, QTKW-5B) identified in the present study were also identified in previous reports, and therefore could be potential candidates for MAS. Several putative candidate genes encode important molecular functions such as tissue-specific expression, light-signaling and photoperiodic modulation of flowering-time, abiotic-stress tolerance, and plant growth and development, which have a role in controlling both agro-morphological and grain-related traits. Further validation and functional characterization of the candidate genes to elucidate the role of these genes in wheat is envisaged.

Conflicts of Interest:
The authors declare no conflict of interest.