Novel Loci for Kernel Hardness Appeared as a Response to Heat and Combined Heat-Drought Conditions in Wheat Harboring Aegilops tauschii Diversity

: Kernel hardness inﬂuences the milling and baking quality of wheat. Stress environments such as heat and combined heat-drought can produce harder kernels, thereby affecting the overall wheat quality. Beside puroindoline genes that are known to determine hardness, other QTLs contribute to the hardness. These QTLs, especially under stress conditions, need extensive research. Moreover, understanding the modiﬁcation or stabilization of hardness under stress condition and the relationship with stress tolerance will facilitate the selection of superior lines that maintain both high yield and quality even under the stress environment. Therefore, in the current work, we aimed to identify the genetic loci and marker trait associations (MTAs) that contributes for hardness under optimum conditions in Japan, and heat and combined heat-drought (HD) conditions in Sudan. We used a panel of multiple synthetic derivatives (MSD) having diverse Aegilops tauschii genome segments and investigated the association between hardness stabilization and stress tolerance. Under stress conditions, we observed that less reduction of kernel weight is associated with either low change or stable kernel hardness. We identiﬁed 47 markers associated with hardness under all conditions; the D genome was the main contributor. For the ﬁrst time, we found a signiﬁcant association with hardness under stress conditions on chromosome 4D. We dissected several candidate genes associated with the change of hardness under stress conditions. Our results will improve the understanding of the genetic factors that affect wheat hardness stability.


Introduction
Kernel hardness is an important quality trait that greatly influences the milling and baking quality of wheat. The world wheat trade is largely carried out based on hardness grades. Wheat hardness is a quantitative trait with classes ranging from soft to hard [1]. Two genes mainly determine hardness, puroindoline a and b (Pina and Pinb), which are located on the short arm of chromosome 5D (5DS) and form the molecular basis of wheat hardness [2]. The presence of the wild-type form of these genes is associated with the soft kernel phenotype in hexaploid wheat. The absence or mutation in either gene results in a hard kernel phenotype [3].
Although puroindoline genes are considered major determinants of wheat hardness, many studies have demonstrated the complex nature of this trait [4,5], and suggest that hardness is affected by several factors. These factors include abiotic stresses such as heat and drought. Heat and drought stresses pose a serious risk to agriculture [6,7]. Several studies have demonstrated that stress from high temperatures and drought can accelerate kernel filling. These stresses compress the timing of key events during wheat kernel development, such as increasing the production of storage proteins and starch synthesis in the endosperm under stress, which together affect kernel hardness [8]. Furthermore, kernel hardness has been reported to be negatively associated with the thousand-kernel weight [9,10], as well as with shape traits such as kernel length, width, and diameter [9]. Proteins are considered important components of wheat grains since they govern wheat's end-use quality. Alterations of the protein fraction composition due to drought and heat stress are primarily due to changes in the total nitrogen quantity accumulated during the seed-filling phase [11]. In turn, they increase hardness because the protein content is positively correlated with hardness [12].
In addition to the puroindoline genes, several studies have identified alleles that contribute to hardness; these were located on most of the 21 wheat chromosomes, but especially on chromosomes 1B, 2A, 4B, 5A, 5B, 5D, and 7D [1,[13][14][15]. Several researchers studied the hardness of plants grown under stress condition [16,17]. However, the genetic factors that affected the change in hardness remained unclear. Therefore, better understanding of the change or stability in hardness under stress environment is essential because the factors or QTLs contributes to hardness stabilization or modification under stress environment will have a great value for wheat quality breeding. Therefore, we conducted this study to dissect the genetic factors that contribute to hardness under optimum, heat, and combined heat-drought (HD) conditions and to investigate the association between hardness stabilization or modification and stress tolerance.
We used a panel of multiple synthetic derivatives (MSD) lines that harbor a wide diversity from a wild species, Aegilops tauschii. Our results revealed a significant marker trait association on chromosome 4D under heat and combined heat-drought (HD) environments. The tolerant MSD lines, with low reduction in their kernel weight, have more stable hardness than sensitive ones. The germplasm source identified and characterized in this study will be an excellent source to breed wheat stress-tolerant cultivars with stable hardness.

Plant Materials
In this study, we used a multiple synthetic derivative (MSD) population of 400 lines developed by crossing and backcrossing 43 synthetic wheat lines with the Japanese wheat cultivar 'Norin 61' [18]. Based on heading time and vernalization requirements, 140 MSD lines were selected out of the 400 original lines and tested in Sudan under heat and HD conditions. All the MSD lines were genotyped using the DArT-seq platform [18].

Field Experiment
The field experiments were conducted in Japan and Sudan. We choose Japan as the optimum condition, because it is considered as having favorable conditions for wheat, while we selected Sudan for the heat and HD experiment because it has been recognized as the global platform for heat tolerance research [19]. In Japan, the 400 MSD lines were grown in a field at the Arid Land Research Center, Tottori University for two seasons (2015/2016 and 2018/2019), using an augmented randomized complete block design with eight blocks. We used four replicated checks in each block; "Norin 61" (the MSD parent), "Imam" and "Tagana" (Sudanese heat-tolerant cultivars), and "Safedak Ishkashim" (a Tajikistan landrace) [19]. In Sudan, we selected 140 lines because they do not require vernalization treatment and adapt to the Sudanese conditions. Those lines were evaluated under heat and HD conditions using alfa-lattice design, with two replications. The drought was imposed by withholding the irrigation after heading as described in our previous study [19].

Hardness and Protein Content Measurements and SEM Observation
The measurement of kernel hardness (index) of the MSD lines under optimum, heat, and HD conditions was performed with a bulk sample of 100 clean, unbroken wheat kernels using the single-kernel characterization system (SKCS 4100, Perten Instruments, Waltham, MA, USA) at the NARO Western Region Agricultural Research Center. This machine measures the force required to crush individual grains of a sample between two surfaces using (index) unit. The MSD lines were classified into hard and soft based on significant differences from the value for the background cultivar 'Norin 61'.
We measured the nitrogen content of about 50 mg of wheat kernel powder using a CN Corder (Model MT-700; Yanaco, Inc., Kyoto, Japan). We recorded the nitrogen to carbon ratio, and then calculated the nitrogen content as a percentage and converted that value to the crude protein content by multiplying the N content by the conversion factor 5.95 [20].
To detect physical changes in response to stress, we examined the internal structure of the wheat kernels using a transverse section of the kernel under a scanning electron microscope (JSM-6610LV, JEOL, Peabody, MA, USA).

Statistical Analysis and Hardness Stress Indexes
We performed an analysis of variance (ANOVA) for the data from the Japanese experiment to detect differences between accessions using version 1.4 of the Plant Breeding Tools software (PBTools, http://bbi.Irri.org/products), Laguna, Philippines. We estimated the best linear unbiased predictions (BLUPs) across the two seasons and used the results for GWAS. For the Sudan experiments, we performed ANOVA using GenStat 18 (VSN International, Rothamsted Research, Harpenden, Hertfordshire, UK). Broad-sense heritability (H 2 ) was defined as H 2 = V G /(V G +V E ), where V G and V E are the genetic and environment estimation, respectively [21]. We calculated relative performance for hardness under heat (HI) and HD conditions (HDI). Also, we calculated heat susceptibility index (HSI) and combined heat-drought susceptibility index (HDSI). All calculations were conducted using the following formulas: HI = (performance in heat environment/performance in optimum environment) × 100% (1) HDI = (performance in HD environment/performance in optimum environment) × 100% (2) where Y h is the phenotypic means for each genotype under heat or HD condition, Y is the phenotypic means for each genotype under optimum condition, X h is the means for all lines under heat or HD condition, and X is the means for all lines under optimum conditions.

Association Analysis
We performed genome-wide association analysis for kernel hardness using a mixed linear model (MLM) with 14,355 DarT-seq markers by version 5 of the TASSEL software [22]. Tassel MLM products were used to generate Manhattan plot using CM plot package in R. We created Manhattan plots using p < 3 × 10 −3 as the threshold to identify significant associations.

Bioinformatics Analysis
To identify candidate genes for hardness and the related proteins, we selected the marker trait associations (MTAs) that were identified as significant under all conditions, and we BLAST them against the International Wheat Genome Sequencing Consortium (IWGS) RefSeq V.1 chromosomes, using URGI with BLAST option (https://urgi.versailles. inra.fr/blast/), accessed on 25 May 2021. Then, we searched for the candidate genes with high confidence in distance (±500 kbp) for the genome region. We used version 1.0 of IWGSC_Ref_seq to search for genes with high confidence. To identify the protein function, we used version 1.1 of IWGSC_Ref_Seq_Annotations along with EnsemblPlant (https://plants.ensembl.org), accessed on 25 May 2021. We investigated the expression levels of all candidate genes that highly contribute to hardness and compared them to the expression of the puroindoline genes using the Wheat Expression Browser expVIP [23]. This lead to understand the association between the candidate genes and hardness. This browser allows comparisons across several studies by taking an input of RNA-seq that collected from these studies. The output is viewable as browser interference with interactive filtering, sorting and export option. This allows an easy access for the researchers [23]. The visualization of the expression was realized using GENEVESTIGATOR software ( https://genevestigator.com/gv/start/start.jsp), Tokyo, Japan, accessed on 25 May 2021.

Phenotypic Variation and Diversity in Hardness among MSD Lines under Optimum and Stress Conditions
Under an optimum environment, significant variation (p < 0.001) in hardness was detected among the MSD lines. No significant differences were detected between the two seasons (Table 1), while the interaction between the environment and seasons (G × S) was significant (p < 0.05). The MSD lines exhibited a wide range of variation for hardness (18.5 to 47.8), whereas, the backcross parent (Norin 61) was 29.6. We selected 140 MSD lines out of the original 400 lines based on their suitability to grow under Sudanese conditions to understand the effect of the heat and HD environments on the hardness and to detect the different factors that influence the hardness under stress conditions. Therefore, the results of the 140 lines will be presented and discussed.
Under stress environments, significant differences were observed among the genotypes (p < 0.05) ( Table 1). The environmental effect was significant, whereas the G × E effect was not significant. High heritability (>90%) was observed under the optimum as well as heat and HD environments (Table 1).
Under heat and HD environments, the genotypes shifted to became harder. This appeared on the frequency distribution, which ranged from 19.6 to 56.9 under heat, and from 25.  Hardness is one of the traits reported to correlate with kernel weight and other kernel-shape traits [9]. In our previous study [19], we discussed the effect of the stress environments on the kernel shape traits, including the kernel weight, and we identified some tolerant genotypes under heat or HD environment based on their kernel weight reduction. The correlation between the hardness and the other kernel-shape traits showed that kernel weight was significantly negatively correlated with hardness under optimum and stress conditions. However, this correlation became stronger under stress environments (Table S1). All the kernel-shape traits except kernel weight correlated with hardness only under stress environments (Table S1).
One of our objectives in the present study was to investigate whether the stabilization of kernel hardness under stress condition was associated with stress tolerance. Since kernel weight is correlated with hardness, the kernel weight susceptibility indexes, which assess the reduction under stress compared with optimum conditions. Thus, it can be used as tolerance indicators. We investigated the relationship between hardness and kernel weight for heat and heat-drought susceptibility index (HSI and HDSI) for the MSD lines ( Figure 2). We observed a relationship between stress tolerance (a low kernel weight reduction) and hardness stabilization under stress condition. The kernel weight HSI ranged from 0.44 to 1.55, and the kernel weight HDSI ranged from 0.67 to 1.36. The hardness HSI ranged from −2.5 to 9.2, and the hardness HDSI ranged from −2.9 to 8.4. Genotypes with HSI and HDSI values less than 0.75 and 0.79, respectively (i.e., kernel weight reductions less than 20 and 27%, respectively) were considered to be tolerant, whereas the genotypes with HSI and HDSI values greater than 1.3 and 1.2, respectively (i.e., kernel-weight reductions greater than 35 and 41%, respectively) were considered stress-susceptible.
There was a clear association between having stable hardness and tolerance to heat or HD stress. The kernel weight HSI showed that 17 genotypes had a low kernel weight reduction, of which 12 had only a slight change in hardness, which remained stable or became softer (as shown in red color in Figure 2a). Based on HDSI, ten genotypes had a low kernel weight reduction, of which six had a low change in hardness (as shown in red color in Figure 2b). In contrast, among 14 genotypes with high kernel weight reduction, 8 were largely affected in their hardness as shown in yellow color in (Figure 2). Of the 10 genotypes that had a high reduction in kernel weight under HD, 6 of them had a high change in hardness (with yellow color Figure 2). Among these genotypes, we found five genotypes shared for both HSI and HDSI ( Figure 3). The tolerant line MSD187, which we identified in our previous study [19], had low change in hardness under both heat and HD (Figure 3a), whereas the sensitive line MSD259 showed a large reduction in hardness under both heat and HD (Figure 3b). The cultivars 'Norin 61', which is the parent of the MSD lines, and 'Imam', had relatively high reductions in kernel weight and moderate changes in hardness under heat conditions compared with MSD187 ( Figure 3a). We further tested these lines for their internal structure and protein content to investigate the relationship between the kernel weight reduction and stabilization or modification of hardness.

Internal Structure and Protein Content
We used a scanning electron microscope (SEM) to examine the internal structure of the mature kernels of the selected tolerant (MSD187), sensitive (MSD259), 'Norin 61' accessions and of the 'Imam' (Figure 4a). The starch granules of 'Norin 61' were loosely packed under optimum conditions as well as heat and HD condition, which indicates a soft kernel. 'Imam', which was hard under optimum as well as heat and HD, showed embedded starch granules with no gaps, which indicates hard kernels (Figure 4a). We observed that the surface of starch granules of the tolerant line MSD187 were separated from the protein matrix under optimum, heat and even under HD with no change due to the stresses (Figure 4a). In contrast, the starch granules of the sensitive line MSD259 were separated under optimum; however, under heat and especially HD, starch granule were tightly attached to the protein matrix, and were embedded in the protein matrix, indicating an increase in hardness as response to the stresses. We observed an increase in protein content under heat conditions in the four tested lines comparing to optimum condition (Figure 4b), which could be due to the big variation between optimum in Japan and stress conditions in Sudan. However, in the Sudanese conditions, when HD and heat were compared, the protein content of the tolerant line MSD187 decreased compared to the sensitive line MSD259 and 'Norin 61' (Figure 4b). This indicated that a decrease in protein content is associated with a low change in hardness.

Marker Traits Association of Hardness and Hardness Indexes under Optimum and Stress Conditions
In the current study, we performed a genome-wide association study (GWAS) to identify genetic loci that contribute to hardness under optimum, heat, and HD conditions. We found a total of 47 statistically significant markers ( Table 2). The Manhattan plots revealed several significant associations for hardness under all conditions ( Figure 5). The A genome contributed to hardness only under stress conditions, in which chromosome 1A under HD, and 4A under both heat and HD environments showed significant MTAs ( Figure 5, Table 2). In the B genome, only chromosome 5B contributed to hardness under optimum as well as stress condition. The D genome contributed more strongly to hardness than the A and B genomes, with the D genome explaining 91.4% of the total contribution under optimum, heat, and HD environments.
We observed MTAs on chromosome 5D under all three conditions ( Figure 5, Table  2). The marker with the highest significance was 1,127,970 with 22.6% of the phenotypic variation explained (PVE) under optimum environment. Also, significant MTAs were observed on chromosome 4D, under heat and HD ( Figure 5). Interestingly, these association appeared only under the stress environments. The marker with the highest significance was 1062681|F|0-26 on chromosome 4D was associated with heat and HD with a PVE of 20.5 and 17.3%, respectively. In addition, we obtained significant associations on chromosomes 6D and 7D under stress condition, but with a few significant markers under optimum conditions. Marker 1074408|F|0-12 on chromosome 6D was associated with heat and HD with PVE of 13.7 and 13.1%, respectively. Marker 3222372|F|0-27 on chromosome 7D was associated with heat conditions and had a PVE of 19.8%.
Hardness indexes revealed 45 MTAs for hardness heat index and hardness heat drought index (HI, HDI) distributed in the A, B, and D genomes. The highest contribution was from the D genome (64%, of which 62% was located on chromosome 4D), versus 22 and 14% for the A and B genomes, respectively (Table S2, Figure S1). Marker 1062681|F|0-26 on chromosome 4D was highly significant under heat and HD environments and was associated with both hardness and the hardness indexes, indicating that it contributed only under stress condition (Table S2). This marker explained 20.5, 17.3, and 17.5% of the PVE of heat, HD and HI, respectively ( Table 2 and Table S2).

Common and Specific MTAs Associated with Kernel Hardness under Optimum and Stress Conditions
We identified stable markers for hardness across the three conditions, mainly in the D genome on chromosomes 5D and 6D (Figure 6a). We also detected significant markers in the D genome associated only with stress conditions on chromosomes 4D, 5D, 6D, and 7D ( Figure 6b). Meanwhile, we identified markers from the A, B, and D genomes associated only with the two hardness indexes (Figure 6c, Table S2), indicating that these markers were associated with the modification of hardness in response to stress condition (Figure 6c). We identified markers associated with stress for hardness and the hardness indexes (Figure 6d). Most of the markers that we detected in the D genome were common across all conditions, whereas markers in the A and B genomes were specific to the HD condition (Figure 6e). Under the optimum condition, there were significant associations with hardness only on the short arm of chromosome 5D, however, under the stress environments of heat and HD, there was another association on the long arm ( Figure 5). This indicates that there are gene(s) on the chromosome 5D other than puroindoline contributing to the hardness under stress conditions. We noticed that all the markers associated with kernel weight and shape traits on chromosome 5D under the stress environments located on the long arm [19]. This result led us to investigate the relationships between markers for grain hardness and the kernel weight or shape traits under stress on this chromosome arm. However, we did not identify pleiotropic markers for both hardness and the kernel weight or shape traits ( Figure S2).

Candidate Genes for Hardness and Gene Expression
We searched for candidate genes associated with the significant markers. We targeted the markers with a high PVE combined with a function related to hardness. The resulting candidates are listed in Table 3. Marker 1,127,970 on chromosome 5D was stable under all conditions and it encodes a puroindoline gene. Marker 3,532,985 on chromosome 5D detected under optimum and heat environments encoded for sucrose transporter gene SUT. We further investigated the SUT allele's contribution to hardness under optimum and heat environments (Figure 7a-c). We found that the "C" allele and "N" allele were associated with decreased hardness, whereas, the "A" allele was associated with increased hardness (Figure 7a,b). However, when we checked a stress-tolerant line (MSD187), a sensitive line (MSD259), and their parent 'Norin 61', we found that both the sensitive line and 'Norin 61' harbored the "C" allele that decreased hardness, whereas the tolerant line harbored the "N" allele ( Figure 7c). Marker 3,947,128 on chromosome 1A was associated with hardness under HD and encodes for antigen that causes celiac disease in human, which is an autoimmune disease that triggered by wheat gluten. We investigated the allele's contributions to celiac disease (Figure 7d-f) and found that under optimum condition, the allele's contribution was almost the same (Figure 7d). However, under HD, allele "C" contributed to increase hardness. In contrast, the allele "A" contributed to decrease hardness. Markers 3944483, 1073587|F|0-24, and 5,357,783 were common and stable markers under all environments and encodes for NB-ARC protein. Marker 998809|F|0-7 on chromosome 4D encodes an ATP/ADP transporter. Moreover, we found that a marker associated with a candidate no-apical-meristem (NAM) gene was associated with marker 1,099,828 on chromosome 3B. We detected the expression for the candidate genes using the expression data from expVIP databases [23] and compared their expression to puroindoline genes (Pina and Pinb) expression. The expression of puroindoline was high on seed part such as endosperm, starchy endosperm, seed coat and aleurone ( Figure S3). The candidate gene TraesCS5D02G001200 on the chromosome 5D showed expression on seed parts besides shoot, leaves and root. The candidate gene TraesCS4D02G047400 on the chromosome 4D showed expression on seed part and spike ( Figure S3). The candidate gene TraesCS1A02G00 8000 on the chromosome 1A that encodes for the antigen that cause celiac disease in human showed a high expression level same as puroindoline genes ( Figure S3).

Phenotypic Variation and Diversity for Hardness
Hardness is an important milling quality trait that is known to be controlled by two puroindoline genes (Pina and Pinb) on the short arm of chromosome 5DS [2]. In the current study, we detected the effect of heat and HD environments on hardness in the MSD population. We aimed to identify MTAs, and possible candidate genes that contributed to the stability or modification of kernel hardness under heat and combined heat-drought conditions, and to investigate whether the stabilization of kernel hardness under stress was associated with stress-tolerance genes.
The wide variation observed among the MSD genotypes indicates high diversity in this population, which includes genes derived from 43 Aegilops tauschii accessions [18]. However, most of the genotypes had soft kernels and resembled the backcross parent, 'Norin 61'. Gedye et al. [24] characterized several synthetic wheat lines and found that most of the population had soft kernels.
MSD population shifted to become harder under heat and HD conditions. This could be attributed to the increased protein content under stress condition [25]. In addition, stress affected the plants, causing smaller kernels; this will lead to reduced kernel weight and overall yield [26]. This is due to the difference of the conversion ratio of carbohydrate and protein, i.e., the conversion ratio of carbohydrate is highly affected by the stress but not much in protein, resulting in smaller grains with a higher protein content. These seeds become harder than plump seeds with much starch [26].
Moreover, under heat and drought stress, the short kernel-filling duration results in less starch accumulation in wheat [27]. This, in turn, could also lead to increased hardness since starch content is negatively associated with kernel hardness [1].
Climate change is adversely affecting wheat production around the world, leading to yield losses and decreased quality [6]. Therefore, identifying varieties that have both high yield and high quality is critical for food security in the context of global climate change. These varieties will also be crucial for breeding programs. The resulting varieties will increase satisfaction for all stakeholders in the wheat value chain [28]. We observed a relationship between stress tolerance and stable kernel hardness or being softer under heat and HD environments. This could be attributed to the starch accumulation that occurs under stress conditions. Prathap and Tyagi [29] found that stress-tolerant rice lines can accumulate starch under drought stress due to the action of starch synthase enzymes and that this could lead to softer kernels under stress.
We confirmed our results by measuring the kernel protein content and examining the endosperm's internal structure using an SEM. One important difference between soft and hard kernels was the degree of adhesion between starch granules and the protein matrix; starch granules in soft kernels were separated from the protein matrix [30]. Hard kernels instead had tight adhesion between starch granules and the protein matrix, and the starch granules were embedded in the protein matrix [31]. The stress-tolerant line MSD187 showed separated starch granules under heat and HD conditions. In contrast, the sensitive line MSD259 showed embedded starch granules with no gaps under stress condition. This indicates that the sensitive line's kernels became harder under stress condition.
Under HD conditions, the protein content increased in the sensitive line, whereas the protein content for the tolerant line decreased. This result suggests that a smaller change in hardness is associated either with less change in the protein content or with decreasing protein content. Protein content is known to increase under stress condition. Peterson et al. [32] reported that one of the factors associated with increasing kernel hardness was an increase of the protein content. All these findings suggest that the stress-tolerant MSD lines could maintain their kernel hardness under stress conditions. These lines therefore represent an excellent resource for breeding stress-tolerant wheat genotypes with high yield and a range of end-use qualities.

Marker Traits Association for Hardness under Optimum and Stress Conditions
Our results indicated that the D genome contributed strongly to kernel hardness under all conditions, with a diverse range of D-genome markers associated with hardness. Among the MTAs in the D genome, we detected a significant association for marker 1,127,970 on chromosome 5D, whose associated gene encodes a puroindoline protein. This marker was stable across all conditions. We observed a high expression level for puroindoline on seed parts such as endosperm, starchy endosperm, and aleurone, whereas it showed low expression levels in shoots and roots. Similar results were reported by Kiseleva et al. [5].
In addition to the puroindoline genes, other markers associated with hardness were found on chromosome 5D under optimum as well as stress environments. Among these markers, marker 3,532,985 far with (1.52) Mb from Pinb contributed to hardness with 16.1 and 12.1% PVE under optimum and heat environments, respectively. The candidate gene TraesCS5D02G001200 for this marker encodes the sucrose transporter TaSUT. Starch formation in wheat kernels requires importation of sugar in the form of sucrose via a sucrose transporter [33]. TaSUT has five homologs, including TaSUT2-5D, which is located on chromosome 5D [34]. A recent study by Al-Sheikh Ahmed et al. [35] showed that TaSUS1 is associated with increased kernel weight under drought stress. This may explain the association between kernel weight and hardness. However, further investigation of the alleles in the stress-tolerant and sensitive lines indicate that both had alleles that decreased the hardness. This finding demonstrates that the activity of the sucrose transporter is not sufficient on its own to stabilize kernel hardness under stress. We found high expression of this gene in all wheat tissues, and this is logical since those genes are associated with sink-source relationships and starch formation.
We detected a marker 3,944,483 on chromosome 5D that was stable under all conditions. This marker is associated with the candidate gene TraesCS5D02G037000, which encodes an NB-ARC protein. The NB-ARC is a stress-resistance protein that contains a central nucleotide-binding domain [36]. It works as a dynamic signaling molecules performing reversible interaction which confer resistance to a wide variety of microbes [37]. Previously, Giroux et al. [38] reviewed the role of puroindoline genes in the plant defense. These results indicated that there are factors controlling hardness associated with plant defense against pathogens. Beside this marker we also found other markers; 1073587|F|0-24, and 5,357,783 under both heat and HD environments have the same candicey gene TraesCS6D02G105100, which encodes for NB-ARC protein, indicating that this gene also affects hardness under stress condition. However, we detected a low expression level for this gene in all tissues.
We detected a significant association for both kernel weight and hardness on the long arm of chromosome 5D under a stress environment. Previously, we detected a significant peak for kernel size under HD environment on the long arm of chromosome 5D [19]. We also found MTAs that appeared under stress and that were associated with kernel shape traits on the long arm on chromosome 5D [19]. These results indicated the occurrence of alleles that contribute to hardness, kernel weight, and kernel shape traits on chromosome 5D under stress environments and this could explain the relationship between hardness and shape traits under these conditions. However, we could not identify common markers between hardness and kernel weight despite of the phenotypic correlation between them and with the other shape traits. This may be due to the complexity of the MSD accessions, which have huge diversity resulting from the diverse D genome sources. In addition, the candidate gene TraesCS5D02G036300 associated with marker 1,025,407 encodes a zinc finger RING-type domain-containing protein. A recent study in wheat indicated that overexpression of this gene increases tolerance to drought and salinity [39].
On chromosome 1A, the marker 3,947,128 associated with the candidate gene TraesCS1 A02G008000 encodes a protein associated with triggering symptoms of celiac disease in humans. Celiac disease is a multisystem immune-based disease triggered by the ingestion of gluten in genetically susceptible individuals [40]. We noticed a high expression level for this marker in the kernel, endosperm, and aleurone, same as the puroindoline expression level. Ribeiro et al. [40] found a negative but non-significant correlation between hardness and amount of toxic epitopes potentially associated with celiac disease.
We also identified MTAs that were only associated with the hardness indexes HI and HDI, which indicate an association with the change of hardness in response to stress condition. We found that all the A, B, and D genomes contributed to the hardness change under stress environments. Among them, we found that marker 1,099,828 on chromosome 3B contributed to hardness, with PVE of 12 and 10% for HI and HDI, respectively. The candidate gene TraesCS3B02G533100 encodes an NAM gene. Uauy et al. [41] reported that NAM-B1 underlies a high-protein-content locus (GPC-B1) that originated from wild emmer (durum) wheat, Triticum turgidum. This gene accelerates senescence, leading to increased protein content and nutrient content in the kernel (through a sink-source relationship). This gene shows expression on the seed coat, aleurone, and spike.
Puroindoline genes are major genes that control kernel hardness under optimum condition. However, it remains unclear whether puroindoline genes also control kernel hardness under stress condition or whether other genetic factors are responsible. Interestingly, we identified a significant peak on chromosome 4D under heat and HD stress environments as well as hardness indexes HI and HDI. This indicates the occurrence of genetic factors on chromosome 4D that contribute to hardness maintenance under stress environment. We also identified significant peaks on chromosomes 6D and 7D under stress, indicating that the associated loci contribute to kernel hardness under stress condition. In contrast with Wilkinson et al. [42], who found a puroindoline-like gene sequence on the long arm of a chromosome homoeologous to chromosomes 7A, 7B, and 7D, we did not identify such a gene sequence on chromosome 7D, indicating that the occurrence of other genetic factors on chromosome 7D contribute to kernel hardness.
We detected several markers on chromosome 4D. Marker 998809|F|0-7 was associated with stress conditions and had a PVE of 14.4 and 11.1% under heat and HD environments, respectively. This marker encodes an ATP/ADP transporter gene. Wang et al. [43] explained that ATP/ADP transporter genes increase the starch content in transgenic Arabidopsis. This could explain why the kernels of some genotypes became softer under stress conditions. Additionally, we found marker 1062681|F|0-26 contributes for hardness and HI as well as heat and HD environments, with PVE 20.5, 17.3, and 17.4%, respectively. The associated candidate gene TraesCS4D02G047400 encodes glutamine synthetase. This enzyme plays an important role in nitrogen-use efficiency and in the uptake and assimilation of nitrogen [44]. These findings indicate that these MTAs are associated with nitrogen-use efficiency genes under stress environment, which is an important stay-green stress-tolerance mechanism. The stay-green trait results from a balance between N demand by the kernels and the N supply during kernel filling [45]. This suggests that hardness resulted from a stay-green mechanism associated with glutamine synthetase. Pinto et al. [46] showed an association of the stay-green trait in wheat with stress tolerance and identified markers on chromosomes homoeologous to 4A and 4B. Here we found the association on chromosome 4D, demonstrating that the stay-green influence hardness. We observed the same expression as puroindoline in the seed in addition to the spike. Muhitch [47] found glutamine synthetase in the developing kernel of maize plants, and it was abundant in the pedicel, pericarp, and kernel glumes.
We found marker 1001438|F|0-46 to be associated with heat conditions, and its contribution to hardness was high (PVE of 18.9%). The associated gene was close to the Rht-D1 region (2.4 Mb), indicating that dwarfing genes could be one of the factors that affects kernel hardness under heat conditions. Wang et al. [14] studied a recombinant inbred line population derived from crosses between accessions with soft and extra soft kernels and found a QTL on chromosome 4D that was close to the semi-dwarf gene Rht-D1.
Based on what is known about the abovementioned candidate genes, we hypothesize that the expression of those genes would affect kernel hardness. For instance, the candidate sucrose transporter gene, which plays a role in transporting sucrose, promotes starch accumulation in the kernels, and the candidate ATP/ADP transporter gene also increases starch content, which leads to decreased hardness under stress. The candidate NAM gene, which plays a role in source-sink nutrient and protein contents, leads to increased protein content, which increases kernel hardness. Moreover, the candidate glutamate synthase gene plays a role in nitrogen-use efficiency, increasing the uptake and assimilation of nitrogen resulting in high protein content, which in turn produces harder kernels.
We hypothesize that the stress-tolerant plants have a mechanism that utilizes the abovementioned genes efficiently, thereby maintaining the ratio between the starch and protein contents and leading to stabilization of kernel hardness and yield. However, further study will be needed to validate this hypothesis.

Conclusions
Heat and combined heat-drought (HD) environments increase hardness. We observed that the MSD line (MSD187) with tolerance potential to heat and HD environments had more stable hardness than the sensitive line (MSD259). We identified MTAs on chromosome 4D that associate with hardness under stress environment and hardness indexes (HI and HDI). This indicates that occurrence of MTAs contributes to the hardness changes under stress conditions. Along with dissecting candidate genes that can contributes for hardness change specially under stress environment, our findings could play an important role in understanding the factors that control the changes of the hardness under stress environments and their relationships to the stress tolerance. This will help breeders to develop stress-tolerant cultivars that maintain high yield and stable hardness, capable of resisting the adverse effects of global warming, thereby improving food security.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/agronomy11061061/s1, Figure S1. Manhattan plot and Q-Q plot for kernel hardness indices (a) Hardness heat indices HI, (b) Hardness drought indices HDI, (c,d) Q-Q plot for hardness indices. Figure S2. Relationship between all markers on chromosome 5D associated with hardness and kernel weight and shape related traits. Figure S3. Expression of the candidate genes and Pina, Pinb, GSP-1 in different tissues of wheat based on RNA-sequencing data collected from different experiments. Table S1. Correlation between hardness, kernel weight and shape related traits under optimum, heat and HD. Table S2. Marker trait associations of hardness of multiple synthetic derivatives for heat index (HI) and heat-drought index (HDI).