High Resolution Melting and Insertion Site-Based Polymorphism Markers for Wheat Variability Analysis and Candidate Genes Selection at Drought and Heat MQTL Loci

The practical use of molecular markers is facilitated by cost-effective detection techniques. In this work, wheat insertion site-based polymorphisms (ISBP) markers were set up for genotyping using high-resolution melting analysis (HRM). Polymorphic HRM-ISBP assays were developed for wheat chromosomes 4A and 3B and used for wheat variability assessment. The marker sequences were mapped against the wheat genome reference sequence, targeting interesting genes. Those genes were located within or in proximity to previously described quantitative trait loci (QTL) or meta-quantitative trait loci (MQTL) for drought and heat stress tolerance, and also yield and yield related traits. Eighteen of the markers used tagged drought related genes and, interestingly, eight of the genes were differentially expressed under different abiotic stress conditions. These results confirmed HRM as a cost-effective and efficient tool for wheat breeding programs.


Introduction
Wheat is among the most important and widely grown crops worldwide [1] and one of the most important grain food crops in the human diet (https://www.fao.org). Wheat development and yield can be affected by abiotic stresses as drought [2][3][4][5] and heat [6,7], whose frequencies would be increased by the strong effects of the predicted climate change and global warming [8][9][10]. In fact, drought is considered one of the most limiting environmental factors [3,11,12], strongly affecting the growth [13,14] and production of crops, with significant reductions in the final yield of cereals, including wheat [13]. Heat stress usually affects crops during the post-anthesis period, with negative effects on final production [15] and end-use quality products [16]. These abiotic stresses are important 4A [33,61] harbors several QTLs related to biotic and abiotic stresses tolerance, agronomic traits as grain yield and quality, and regulation of physiological traits as plant height, maturity, or dormancy [62][63][64][65][66][67][68][69][70]. This chromosome represents an important target in plant breeding, marker design for variability analysis and candidate genes assessment. In this study, ISBP markers from wheat chromosome 3B [30] and 4A were used to develop and validate HRM assays, and to assess the genetic variability in a wheat collection. These markers were also used to target meta-quantitative trait loci (MQTL) [71,72] related to drought and heat stresses, as well as yield and yield related traits. Candidate gene analyses were performed for the ISBP markers and the genes were validated by gene expression analyses carried out among different drought and heat stress conditions.
Genomic DNA was isolated from young leaf tissue according to the cetyl trimethyl ammonium bromide (CTAB) method of Murray and Thompson [73], as optimized by Hernández et al. [74]. The quality and concentration of samples were assessed by electrophoresis in a 0.8% agarose gel.

Insertion Site-Based Polymorphism Markers Development
ISBP markers were initially developed based on the wheat chromosome 4A survey sequencing [61] and confirmed in the bread wheat reference genome sequence RefSeq v1 [33].
The assemblies corresponding to the 4A wheat chromosome survey sequencing were generated using the "Newbler v2.7" software package (Roche Diagnostics Corporation, Basel, Switzerland) using default parameters. IsbpFinder [30] was run on the assemblies obtained, and ISBP markers were located on the 4AS and 4AL chromosome arms assemblies. The corresponding 45 ISBP primers were designed using Primer3 (http://primer3.sourceforge.net) and mapped to the bread wheat reference genome RefSeq v1 [33]. Marker set up was carried out using 6 durum and bread wheat lines representative of variability (Supplementary Materials Table S1).
ISBP amplicons obtained by standard PCR (55 • C annealing, [30]) using 5 T. aestivum varieties (Supplementary Materials Table S1) were purified by Exonuclease l (Exo I, New England Biolabs, Inc., Ipswich, MA, USA) and SAP treatment (5 µL DNA, 1U Exol, 1xSAP buffer, 1U SAP in 9 µL at 37 • C for 1 h). The purified fragments were then sequenced on an ABI PRISM ® 3730XL (Applied Biosystems, Foster City, CA, USA) genetic analyzer using the forward and reverse ISBP primers [30] and using the ABIPRISM BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). HRM analyses were carried out for the same 5 T. aestivum varieties (Supplementary Materials Table S1) using 6 ISBP primer pairs previously developed for wheat chromosome 3B [30]: HRM3B_273339424, HRM3B_609364064, HRM3B_124761338, HRM3B_203288704, HRM3B_465802537, and HRM3B_331497483. This analysis was carried out in a Rotorgene TM 6000, model 5-Plex real time PCR (Corbett Research, Mortlake, NSW, Australia). The PCR reaction volume was of 15 µL and the mixture composition was: Master Mix of Type-it HRM PCR Kit (Qiagen, CA, USA); 0.7 µM of each primer and 2 µL of genomic DNA (30 ng). The PCR protocol consisted on an initial denaturation step at 95 • C for 10 min; 45 amplification cycles of denaturation at 95 • C for 10 s, annealing at 55 • C for 15 s and a final extension at 72 • C for 20 s; the high-resolution melting was set out by ramping from 65 to 95 • C, with fluorescence data acquisition at 0.1 • C increments (waiting for 2 s every acquisition). HRM results were compared with the amplicon sequencing for the same five samples.
The high-resolution melting analysis was later performed to assign HRM pattern types to the ISBP markers, using 19 durum and bread wheat lines (Supplementary Materials Table S1). The PCR protocol consisted of an initial denaturation step at 95 • C for 5 min; 50 cycles of denaturation at 95 • C for 20 s, annealing at 60 • C for 20 s, and a final extension at 72 • C for 20 s. HRM analysis was undertaken once amplification was completed by ramping from 65 to 95 • C, with fluorescence data acquisition at 0.1 • C increments (waiting for 2 s every acquisition). Results were analyzed by using the Rotorgene software version 1.7 (Qiagen, the Netherlands), and HRM curves were normalized according to the manufacturer's instructions.

Candidate Genes and Gene Expression Analyses
ISBP markers were mapped in wheat chromosomes 4A and 3B, and then were compared to the wheat MQTLs described in Acuña-Galindo et al. [71]. To obtain the MQTLs positions in wheat chromosomes, marker sequences [71] were extracted from "Graingenes" (graingenes.org) and "Wheat SSR DB markers" (wheatssr.nig.ac.jp), and then located by mapping flanking markers (using BLAST) against RefSeq v1 [33]. Only the markers with a corresponding amplicon shorter than 500 bp and a perfect BLAST match (no gap, no mismatch) were considered.
Markers sequences in chromosomes 4A and 3B were blasted against the RefSeq1 v1 [33] using the parameters "-task", "blastn-short" and "-ungapped". The resulting hits were then processed to pair forward and reverse sequence hits with an amplicon <1000 base pairs (bp). For subsequent analysis, paired sequences were ordered by number of mismatches, so markers position was inferred from the position of pairs with lower number of mismatches (0 in most cases). To identify candidate genes associated to each marker, the results were filtered and the hits with best e-value were selected for each molecular marker. The candidate genes were manually selected within a window of +/−20 kb of the marker's hit in the pseudomolecule [33] gene model annotation. Due to the reduced gene density found for wheat chromosome 3B ISBP markers (only three genes were found), that window was extended to +/−300 kb for this chromosome. Genes described as "uncharacterized protein" were then manually annotated. Sequences were obtained in Ensembl Plants (T. aestivum RefSeq v1.1) (https://plants.ensembl.org/), and then searched in UniProt (https://www.uniprot.org/). The annotated hits with e-value 0.0 and a score >2000 were selected, except for the gene TraesCS4A01G410700, which possesses a short length (207aa).
To overview the results from gene expression analyses, heatmaps were drawn using the data retrieved from Wheat Expression (www.wheat-expression.com/) and the 'NMF 0.21.0 R package [75]. The information used was generated by Liu et al. [22], Ma et al. [26], and Galvez et al. [27]. Liu et al. [22] experimental seedling samples grown in controlled conditions were associated to NCBI SRA ID SRP045409 (control (IS), heat and drought (PEG induced drought) stress for 1 and 6 h (PEG1 and PEG6), respectively). Ma et al. [26] experimental samples grown in a shelter corresponded to NCBI SRA ID SRP102636 (anther stage irrigated leaf phenotype (AD_C), anther stage drought-stressed leaf phenotype (AD_S), tetrad stage irrigated developing spike phenotype (T_C), and tetrad stage drought-stressed developing spike phenotype (T_S)). Galvez et al. [27] flag leaf samples from field grown plants used have NCBI SRA ID SRP119300 (irrigated (IF), mild stress (MS), and severe stress (SS) flag leaf samples). Transcripts Per Kilobase Millions (TPMs) of genes under every condition were calculated as mean value of TPMs of its constitutive experiments. A differential gene expression (DGE) analysis was performed using RevSeqv1 [33] gene models through two bioinformatic pipelines: Kallisto (version 0.43.0) with the R library "sleuth"(version 0.28.1), and STAR with the R library DESeq2 (version 1.14.1). A consensus threshold for the two pipelines (|lg2FC, β| > 1 and p-adjust, Q-value < 0.05 [27]) was used.

Wheat Variability Assessment
The PCR amplification protocol used was the previously described for HRM pattern type assessment. Samples genotyping was performed using Melt and HRM analysis options of the Rotor-Gene TM 6000 software. PCR was repeated three times to ensure the amplifications. Results were then corroborated with Rotor-Gene™ ScreenClust HRM™ Software. In those cases where the number of genotypes assigned was unclear, ScreenClust HRM™ Software was also used for the final decision. A binary matrix for genotyping results was created and then analyzed using two different software, PhylTools v.1.32 (Wageningen Agricultural University, The Netherlands) and PowerMarker v.3.25 [76]. The first one was used to calculate the genetic distances for haploid data with "individuals" as hierarchy level and the Nei index [77]. PowerMarker was used to obtain the statistics mean allele number, mean gene diversity, and Polymorphism Index Content (PIC). The unweighted pair group method with arithmetic mean (UPGMA) clustering was carried out using the NEIGHBOR module in the Phylip 3.695 package [78] with default parameters. The selected tree method was UPGMA. The final dendrogram was drawn using MEGA v.6.0 [79] software with the results from the genotyping.
The goodness of fit of the UPGMA tree was calculated by the Cophenetic Correlation Coefficient (CCC) using a visual basic program for Microsoft Excel 2000 [80]. The CCC calculated from the linear regression between the corresponding values of the original distance matrix and the cophenetic matrix derived from the calculation of the UPGMA tree.

Markers Sequence Validation and HRM Pattern Types Assignment
To validate markers sequences, difference plots and HRM normalized curves from the high-resolution melting analyses were obtained for each of the six ISBP markers for wheat chromosome 3B and compared with their corresponding amplicon sequences. The HRM profiles were successfully validated. Sequence polymorphism from one to five nucleotides were detected by HRM analyses.   The HRM pattern type assignment was based on the pattern of normalized high-resolution melting curves obtained, and their potential to genotype a high number of varieties. The ISBP HRM patterns were classified into four different types (Table 1, Figure 2): (i) pattern type A, excellent markers to genotype a high number of wheat varieties simultaneously. The HRM curves are very different and easily distinguishable into classes ( Figure 2a); (ii) pattern type B, good markers to differentiate several groups of varieties (genotypes) in the same run. HRM curves can be differentiated in an easy way ( Figure 2b); (iii) pattern type C, good HRM markers, but not recommended for genotyping a broad variety of samples. The differentiation between HRM curves is not clear in all cases ( Figure 2c); and (iv) pattern type D, assigned to markers which are not recommended for HRM genotyping, due to a low amplification efficiency or to a gradient of HRM curves that hinders precise classification into classes ( Figure 2d). Twelve ISBP markers were classified as HRM pattern type A; 4 markers showed HRM pattern type B; 10 markers had HRM pattern type C; and 13 were considered with pattern type D ( Table 1). Four of the 45 wheat chromosome 4A ISBP primer pairs designed (Table 2), were discarded for the rest of analyses due to non-reliable PCR amplifications. The HRM pattern type assignment was based on the pattern of normalized high-resolution melting curves obtained, and their potential to genotype a high number of varieties. The ISBP HRM patterns were classified into four different types (Table 1, Figure 2): (i) pattern type A, excellent markers to genotype a high number of wheat varieties simultaneously. The HRM curves are very different and easily distinguishable into classes ( Figure 2a); (ii) pattern type B, good markers to differentiate several groups of varieties (genotypes) in the same run. HRM curves can be differentiated in an easy way ( Figure 2b); (iii) pattern type C, good HRM markers, but not recommended for genotyping a broad variety of samples. The differentiation between HRM curves is not clear in all cases ( Figure 2c); and (iv) pattern type D, assigned to markers which are not recommended for HRM genotyping, due to a low amplification efficiency or to a gradient of HRM curves that hinders precise classification into classes (Figure 2d). Twelve ISBP markers were classified as HRM pattern type A; 4 markers showed HRM pattern type B; 10 markers had HRM pattern type C; and 13 were considered with pattern type D ( Table 1). Four of the 45 wheat chromosome 4A ISBP primer pairs designed (Table 2), were discarded for the rest of analyses due to non-reliable PCR amplifications.
Chr: chromosome location; Pattern: high resolution melting pattern type (A: excellent marker to genotype different wheat varieties simultaneously; B: good marker to genotype several wheat varieties simultaneously; C: good marker but not recommended for a broad variety of wheat samples; D: marker not recommended for HRM genotyping; and n/s for markers non-suitable for HRM); Amp. size: amplicon size in base pairs (bp); *: ISBP markers used in the durum wheat collection variability assessment.

Candidate Gene Analysis
After mapping ISBP markers and comparing them with the MQTLs positions previously described in Acuña-Galindo et al. [71], some markers for the wheat chromosome 4A were found in the proximity of interesting QTLs or within MQTLs related to drought and heat stress tolerance, as well as QTLs for yield components (Figure 3a). Two ISBP markers, HRM4A_317085557 and HRM4A_460238681, were located within MQTL30 [71], related to the physiological drought trait and root vigor. Markers HRM4A_617938526 and HRM4A_618105078 were placed in MQTL31 [71], in proximity to QTLs related to drought and heat stresses. The marker HRM4A_660524139 was placed close to two of the QTLs located within MQTL31, associated to traits for yield component and coleoptile vigor, and also heat stress. Marker HRM4A_583704598 was placed between MQTL30 and MQTL31, in proximity to 2 QTL related to drought and heat stresses. Finally, markers HRM4A_681664894 and HRM4A_683608822 were placed in MQTL32 [71], close to QTLs related to drought; and markers HRM4A_702156718, HRM4A_714743756, and HRM4A_716986193 were found close to a QTL located in MQTL32, related to drought tolerance (Figure 3a).

Candidate Gene Analysis
After mapping ISBP markers and comparing them with the MQTLs positions previously described in Acuña-Galindo et al. [71], some markers for the wheat chromosome 4A were found in the proximity of interesting QTLs or within MQTLs related to drought and heat stress tolerance, as well as QTLs for yield components (Figure 3a). Two ISBP markers, HRM4A_317085557 and HRM4A_460238681, were located within MQTL30 [71], related to the physiological drought trait and root vigor. Markers HRM4A_617938526 and HRM4A_618105078 were placed in MQTL31 [71], in proximity to QTLs related to drought and heat stresses. The marker HRM4A_660524139 was placed close to two of the QTLs located within MQTL31, associated to traits for yield component and coleoptile vigor, and also heat stress. Marker HRM4A_583704598 was placed between MQTL30 and MQTL31, in proximity to 2 QTL related to drought and heat stresses. Finally, markers HRM4A_681664894 and HRM4A_683608822 were placed in MQTL32 [71], close to QTLs related to drought; and markers HRM4A_702156718, HRM4A_714743756, and HRM4A_716986193 were found close to a QTL located in MQTL32, related to drought tolerance (Figure 3a).  Two of the wheat chromosome 3B ISBP markers (HRM3B_465802537 and HRM3B_609364064) were mapped within MQTL26 [71], in proximity to QTLs controlling yield component and biomass traits, and also related to heat stress (Figure 3b).
Additionally, as result of the candidate gene analysis, the developed ISBP markers for wheat chromosome 4A mapped in the proximity of 61 genes (Supplementary Materials Table S3 and Figure  S1). The chromosome position for the ISBP and their closest genes are shown in Figure 4a. Two of the wheat chromosome 3B ISBP markers (HRM3B_465802537 and HRM3B_609364064) were mapped within MQTL26 [71], in proximity to QTLs controlling yield component and biomass traits, and also related to heat stress (Figure 3b).
Additionally, as result of the candidate gene analysis, the developed ISBP markers for wheat chromosome 4A mapped in the proximity of 61 genes (Supplementary Materials Table S3 and Figure S1). The chromosome position for the ISBP and their closest genes are shown in Figure 4a. Based on gene expression differences under different drought stress conditions, we filtered 23 genes (37.7% of the total genes) with a TPM value above 2.5 (Table 2 and Figure 5). An expression heatmap, using all available public studies RNASeq data in wheat drought responses [22,26,27] is shown in Figure 5. Based on gene expression differences under different drought stress conditions, we filtered 23 genes (37.7% of the total genes) with a TPM value above 2.5 (Table 2 and Figure 5). An expression heatmap, using all available public studies RNASeq data in wheat drought responses [22,26,27] is shown in Figure 5. severe stress field condition [27]; IS: seedling control; PEG1: seedling 1 h PEG stress; PEG6: seedling 6 h PEG stress [22]; AD_C: anther stage irrigated shelter phenotype; AD_S: anther stage drought stressed shelter phenotype; T_C: tetrad stage irrigated shelter phenotype; and T_S: tetrad stage drought shelter phenotype [26].
The wheat chromosome 3B ISBP markers mapped in the proximity of 49 genes (Supplementary  Materials Table S4 and Figure S2). The closest genes are shown next to the corresponding marker in Figure 4b. Seventeen of these genes (34.69% of the total) showed TPM values above 2.5 (Table 3 and Figure 6). The drought responsive genes mapped by ISBP markers located in proximity or within QTLs or MQTLs are shown in Table 4. severe stress field condition [27]; IS: seedling control; PEG1: seedling 1 h PEG stress; PEG6: seedling 6 h PEG stress [22]; AD_C: anther stage irrigated shelter phenotype; AD_S: anther stage drought stressed shelter phenotype; T_C: tetrad stage irrigated shelter phenotype; and T_S: tetrad stage drought shelter phenotype [26].
The wheat chromosome 3B ISBP markers mapped in the proximity of 49 genes (Supplementary  Materials Table S4 and Figure S2). The closest genes are shown next to the corresponding marker in Figure 4b. Seventeen of these genes (34.69% of the total) showed TPM values above 2.5 (Table 3 and Figure 6). The drought responsive genes mapped by ISBP markers located in proximity or within QTLs or MQTLs are shown in Table 4.  [27]; IS: seedling PEG shock control; PEG1: seedling 1 h PEG stress; PEG6: seedling 6 h PEG stress [22]; AD_C: anther stage irrigated shelter phenotype; AD_S: anther stage drought stressed shelter phenotype; T_C: tetrad stage irrigated shelter phenotype; and T_S: tetrad stage drought shelter phenotype [26]. severe stress field conditions [27]; IS: seedling PEG shock control; PEG1: seedling 1 h PEG stress; PEG6: seedling 6 h PEG stress [22]; AD_C: anther stage irrigated shelter phenotype; AD_S: anther stage drought stressed shelter phenotype; T_C: tetrad stage irrigated shelter phenotype; and T_S: tetrad stage drought shelter phenotype [26].   After differential gene expression analysis we obtained 5 DE genes in chromosome 4A and 3 in chromosome 3B, which were up and down-regulated by the PEG drought treatment [22]. The 4A chromosome DE genes TraesCS4A01G003500 and TraesCS4A01G043500, were up and down regulated under PEG6 drought treatment, respectively. Genes TraesCS4A01G047000, TraesCS4A01G410700 were up-regulated, while the gene TraesCS4A01G069200 was down-regulated under PEG6 treatment. Chromosome 3B gene TraesCS3B01G221100 was down-regulated under PEG6 treatment, while genes TraesCS3B01G290200 and TraesCS3b01G290300 were up-regulated (Supplementary Materials Table S5).

Wheat Variability Assessment by High Resolution Melting Analysis
To assess the polymorphism levels in a wheat collection, the 45 ISBP markers developed for the wheat chromosome 4A (Table 1), were evaluated. Thirteen of them were selected based on their reproducibility and polymorphism (Table 1) and were used in HRM analyses to study the genetic diversity among durum and bread wheat lines in panel 1 (Supplementary Materials Table S1).
The number of alleles detected for these markers ranged from 2 to 7 (mean = 3.38), and the polymorphism index content varied between 0.24 and 0.68 (mean = 0.52) ( Table 5). The cluster analysis shows three clearly differentiated clusters (Figure 7). The first one contains 30 T. aestivum lines, while the remaining 8 bread wheat lines (TaesIN-13, TaesLI-06, TaesLI-07, TaesIF-06, TaesIF-08, TaesIF-07, TaesLI-05, and TspeBO-01) were placed within the second cluster. This cluster also contains Triticum durum and Triticum dicoccoides accessions, while Triticum monococcum and Triticum urartu were placed in the third cluster (Figure 7). The cophenetic correlation coefficient obtained for the UPGMA tree was 0.86.  Seven of these ISBP markers (Table 1) selected by their efficiency and polymorphism for durum wheat, were used for the variability assessment of durum wheat lines in panel 2 (Supplementary  Materials Table S2). The UPGMA tree for the wheat panel 2 is shown in Supplementary Materials Figure S3a,b. This cluster analysis resulted in 9 differentiated clusters. Five lines (BGE002866, BGE013055, BGE020464, BGE013652, and BGE013722) were not placed in any of these clusters. The observed distribution of durum wheat lines across the clusters could be associated in some cases to the geographical location and agroclimatic areas (Supplementary Materials Figure S3c). There are two clear clusters (clusters 4 and 6) where species from southern Spain are presented in a larger proportion. Furthermore, wheat lines placed in cluster 1 are mainly located in norther temperate zones without dry season and temperate summer; while cluster 4 contains landraces mainly located in southern template areas with dry and cool summer. The CCC obtained for the UPGMA tree was 0.67. The number of alleles detected in this analysis ranged from 3 to 6 (mean allele number = 3.86), and the PIC mean value was 0.486 (Supplementary Materials Table S6).

Discussion
Insertion site-based polymorphism markers have been described as a useful tool for wheat genomic studies [35] and attractive alternative to markers as SSR or SNP, due to the high repetitive content of some cereal genomes [81,82]. Due to their straightforward design and high polymorphism, there are previous studies which developed and used specific ISBP markers for several wheat chromosomes and different purposes (i.e., Barabaschi et al. [83] for the bread wheat chromosome 5A, focused on polymorphism assessment; Lucas et al. [84], who used markers derived from the wheat chromosome 1A to map this chromosome, and also with marker assisted breeding purposes; Li et al. [36] applied wheat chromosome 3B ISBP linked to mildew resistance genes in durum wheat; or Sehgal et al. [41], who used chromosome 3A ISBP markers for gene discovery and physical mapping). In this work, new ISBP markers were developed for the wheat chromosome 4A, which contains interesting genes related to biotic and abiotic stresses, as drought tolerance [65,68,69]. It also harbors loci related to essential agronomic traits as yield and grain quality [62,64,66]. The ISBP markers resulted highly polymorphic ( Table 5). Thirteen of the developed ISBP markers were used for wheat variability assessment and showed melting curve polymorphisms, seven of them with a PIC value higher than 0.50. Thus, they resulted highly resolutive tools for wheat variability assessment using a cost-effective technique as high resolution melting analysis (HRM). This technique is described as an optimized methodology for melting curve assessment, which allows the determination of melting temperature and profile of an amplicon [85]. Some studies have highlighted some advantages for the HRM technique, as its reduced cost per sample in comparison with other techniques used for SNP detection [86][87][88]. It is worth noting that the required system consists of standard and affordable RT-PCR equipment, which is suitable for in-house genotyping and adequate for small/medium breeders. Other advantages for HRM are the excellent results for the detection of homozygous and heterozygous variants [46,[89][90][91]; its use for gene mapping, SNPs and mutations [46,90,92]; and its efficiency for the identification of species and closely related varieties [88,[93][94][95][96][97]. In this regard, our results from ISBP markers sequence validation, as well as HRM genotyping, support the efficiency of HRM analysis in wheat varieties differentiation. Our results are in agreement with Dong et al. [53], who highlighted that HRM does not require any digestion or gel electrophoresis, so it provides a worthwhile approach for SNP/indel genotyping of different varieties without prior sequence knowledge, as required by other methods. Results from the wheat genetic diversity assessment confirmed that HRM is a convenient way for a first screening to determine variability groups, prior to resequencing only representative varieties as the basis to develop other SNP platforms. Nevertheless, some HRM limitations have been pointed out. Distefano et al. [88] highlighted that sometimes, the HRM profiles could be similar preventing the differentiation of some of the genotypes. Regardless of this, Wu et al. [98] proposed that this issue can be solved using mixed strategies. Accordingly, our results show that a combined use of different ISBP markers can differentiate all the wheat lines studied.
Twenty ISBP markers for wheat chromosomes 4A and six for chromosome 3B mapped to interesting candidate genes, mainly related to drought and heat stresses and yield components (Tables 2 and 3). They were validated using data from available RNASeq public studies in wheat drought responses [22,26,27], as they showed differences on their expression under different stress conditions.
In chromosome 4A, the ISBP marker HRM4A_109848074 mapped next (97 bp) to the gene TraesCS4A01G098300, which encodes a xyloxyltransferase 1, and participates in carbohydrates metabolism in the development of cellular walls [99]. This process is markedly affected by water stress [100,101] and this can be observed in Figure 5, where this gene decreases its expression under severe stress field conditions. This result is also in agreement with results found by Abebe et al. [102], who analyzed spikes of barley grown under controlled drought conditions, and also found that this gene was down-regulated. The marker HRM4A_2791416 mapped 814 bp upstream to the gene TraesCS4A01G003600, which encodes an alpha/beta-hydrolases superfamily protein with functional adaptability in plants [103]. This gene reduces its expression as drought stress increases under field conditions ( Figure 5). Marker HMR4A_67413676, mapped 975 bp to the gene TraesCS4A01G069200, which encodes an armadillo repeat-only protein. These kind of repeat proteins participate in the coordination of protein interactions during stress and hormonal signalling in plants [104]. In agreement with this, the gene was downregulated under PEG6 drought treatment (Supplementary Materials  Table S5). The marker HMR4A_683608822, which was located within the drought stress tolerance MQTL32 [71] (Table 5 and Figure 3a), mapped 2685 bp upstream to the gene TraesCS4A01G410700, which encodes a ras-related protein RABC2a. The function of this gene has been related to ABA induced stress tolerance in barley [105]. Accordingly to a drought stress response role, this gene was upregulated under PEG6 drought treatment (Supplementary Materials Table S5). The marker HMR4A_36371442 mapped 2736 bp upstream to the gene TraesCS4A01G043500 and was downregulated under PEG6 drought treatment (Supplementary Materials Table S5). Contrary to this, this gene increases its expression under severe stress field conditions ( Figure 5). This gene encodes a STAS domain containing-protein, which plays a role to membrane attachment of many anion transporters in transport activity and regulation in plants [106]. In fact, it has been demonstrated its crucial role in the activity of sulfate transporter in Arabidopsis thaliana [107], providing key amino acids in the sulfate transport activity [108]. The importance of this activity should be noted, since sulfate is an element that has been described as an essential component in the structure of plant enzymes and reserve proteins in grain [109]. Further, marker HRM4A_716986193, which was located close to a QTL within MQTL32 [71] ( Table 5 and Figure 3a), mapped in proximity (3576 bp upstream) to the gene TraesCS4A01G671200LC. This gene encodes a peptidase M20/M25/M40 family protein, which is involved in drought stress responses [110]. In fact, proteolysis under drought conditions allows a reorganization in the plant's metabolism, and also increases plants drought tolerance [20,111,112]. According to this, the results showed increased expression of this gene under different drought stress treatments ( Figure 6). This is also in agreement with the results showed by Simova-Stoilova et al. [113], who assessed wheat leaves under severe soil drought and found an increase in peptidase activity. Thus, this drought responsive genes represent an interesting candidate for a known drought and heat stress tolerance MQTL.
Additionally, within the window of +/−20 kb used for the candidate gene analysis in chromosome 4A, there were two interesting genes: TraesCS4A01G003500 (5141 bp upstream from marker HRM4A_2791416) and TraesCS4A01G047000 (14,885 bp upstream from marker HRM4A_38654555). The first gene encodes a thionin like-protein gene, which plays an important role in the growth and development of the plant and its defense against pathogens [114]. It was found differentially expressed under PEG drought treatment, being upregulated under PEG6 drought treatment (Supplementary  Materials Table S5). The gene TraesCS4A01G047000 encodes a formin-like protein, which plays a primary role in the organization of plant's structure [115]. In agreement with our results where this gene was found upregulated under PEG6 drought treatment (Supplementary Materials Table S5), formin-like proteins showed variations in their expression under drought conditions in wheat [115]. Some HRM/ISBP markers mapped to previously described MQTL loci [71] (Table 5 and Figure 3a), which were mainly associated to drought and heat stresses tolerance. Within these markers, HRM4A_618105078, which tags the MQTL31 [71] and mapped close (2442 bp upstream) to the gene TraesCS4A01G497800LC can be highlighted. This gene encodes a receptor-like protein kinase, which is involved in abiotic stress responses [116], matching the description assigned to the MQTL.
Additionally, within the window of +/−300 kb, four wheat chromosome 3B ISBP markers can be highlighted. The marker HRM3B_124761338 mapped 6356 bp downstream to the gene TraesCS3B01G138700, which encodes a ribonucleoside-diphosphate reductase, an essential enzyme for DNA synthesis [117,118]. This gene increases its expression under severe field stress conditions, and it decreases in response to a PEG drought treatment ( Figure 6). Marker HRM3B_609364064 mapped 31,936 bp downstream to the gene TraesCS3B01G575000LC, encoding myosin-1. Plant myosins have a functional role in organelle movement in response to biotic and abiotic stresses [119]. This response is shown in Figure 6, where this gene significantly increments its expression under PEG1 and PEG6 drought treatments. Marker HRM3B_465802537 mapped to two interesting genes, the gene TraesCS3B01G290200 (83,072 bp downstream) and the gene TraesCS3B01G290300 (125,805 bp downstream), which were both found upregulated under PEG6 drought treatment (Supplementary  Materials Table S5), and contrary to this, decreased their expression under severe stress field conditions ( Figure 6). The gene TraesCS3B01G290200 encodes a glycosyltransferase, an enzyme which possesses a main role in plant's stress tolerance and defense [120], and in agreement with our results, it has been previously found upregulated in wheat leaf under drought conditions [121]; and the gene TraesCS3B01G290300 encodes an ABC transporter B family protein, which is significantly involved in organs growth, plant nutrition and development and plant responses to abiotic stresses [122]. In fact, as Rampino et al. [123] highlighted, and in agreement with our results (Supplementary Materials  Table S6), the up-regulation of this gene in wheat under heat and drought conditions confirms the implication of this gene family in drought responses. Therefore, this family protein has been related to grain formation in wheat [124], which is consistent with the location of marker HRM3B_465802537 within MQTL26 [71], mainly related to yield components. Thus, this marker can be useful in wheat breeding, for the marker-assisted selection of this interesting gene and MQTL. Finally, the marker HRM4A_273339424 mapped 89,674 bp upstream to the gene TraesCS3B01G221100, which encodes a protein kinase superfamily protein. This protein's family is involved in plants' responses to abiotic stresses and plants' development [125]. According to this, this gene was found downregulated under PEG6 drought treatment (Supplementary Materials Table S5), and it also shows differences in its expression across different stress conditions ( Figure 6). Therefore, these results agree with Wei et al. [126], who highlighted that kinase proteins are involved in various responses with exposure time of drought.
According to our results, the developed HRM-ISBP markers can be used in wheat breeding programs to genotype interesting genomic regions in a cost-effective manner. These markers can useful resources for marker-assisted selection (MAS) focused on abiotic stress responses and yield components, to tag interesting known QTLs and MQTLs related to drought and heat stresses tolerance, and also yield-related traits.

Conclusions
In this work, highly polymorphic ISBP markers for wheat chromosome 4A were successfully developed and applied in a genetic variability assessment of a collection of durum and bread wheats, using the high-resolution melting analysis technique. These HRM-ISBP markers represent cost-effective and efficient tools for wheat breeding programs focused on variability assessments. The obtained results provide an interesting framework for wheat genetic studies and varieties selection. These HRM-ISBP markers have also been shown useful for tagging interesting genes associated to drought and heat stresses tolerance, some of which showed differential expression patterns under stress conditions. In addition, some of these markers can be applied in breeding through marker-assisted selection of QTL and MQTL related to abiotic stresses as drought and heat, and also yield and yield related traits. The resources and results presented here can also facilitate the understanding of important traits in other species with large genomes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/9/1294/s1; Table S1. List of durum and bread wheat lines used for variability analysis by high resolution melting. 1 : marker set up analysis; 2 : marker's sequence validation by HRM analysis; 3 : HRM pattern types obtained; IFAPA: Instituto Andaluz de Investigación y Formación Agraria, Pesquera, Alimentaria y de la Producción ecológica; INRA: Institut National de la Recherche Agronomique; and WGRC: Wheat Genetics Resource Center; Table S2. Durum wheat lines used for variability analysis by high resolution melting (HRM); Table S3. Candidate genes tagged by the developed HRM-ISBP markers for wheat chromosome 4A. Genes with expression differences are shown in bold and DE genes are indicated with "*". The positive and negative values in the "Dist" column indicate if the corresponding gene is downstream or upstream of the respective marker. Chr: chromosome location. Dist (bp): distance from the gene to the marker in base pairs; Table S4. Candidate genes tagged by the HRM-ISBP markers for wheat chromosome 3B. Genes with expression differences are shown in bold and DE genes are indicated with "*". The positive and negative values in the "Dist" column indicate if the corresponding gene is downstream or upstream of the respective marker. Chr: chromosome location. Dist (bp): distance from the gene to the marker in base pairs; Table S5. Differential expression (DE) analysis significance parameters for DE genes. Significant values (|lg2FC, β|> 1 and p-adjust, Q-value < 0.05) are shown in bold; Table S6. Genetic parameters for wheat chromosome 4A ISBP markers used in durum wheat variability assessment. HRM type: high resolution melting pattern type; No. alleles: number of alleles found with the marker; PIC: polymorphism index content; and *: one of the genotypes was described as "null genotype"; Figure S1. Gene expression analysis under different water stress conditions for all candidate genes located within a +/−20 kb window to the wheat chromosome 4A ISBP markers. Genes with differences on their expression are shown in bold and DE genes are indicated with "*". IF: irrigated field conditions; MS: mild stress field conditions; SS: severe stress field conditions; IS: seedling PEG shock control; PEG1: seedling 1 h PEG stress; PEG6: seedling 6 h PEG stress; AD_C: anther stage irrigated shelter phenotype; AD_S: anther stage drought stressed shelter phenotype; T_C: tetrad stage irrigated shelter phenotype; and T_S: tetrad stage drought shelter phenotype; Figure S2. Gene expression analysis under different water stress conditions for all candidate genes located within a +/−300 kb window to the wheat chromosome 3B ISBP markers. Genes with differences on their expression are shown in bold and DE genes are indicated with "*". IF: irrigated field conditions; MS: mild stress field conditions; SS: severe stress field conditions; IS: seedling PEG shock control; PEG1: seedling 1 h PEG stress; PEG6: seedling 6 h PEG stress; AD_C: anther stage irrigated shelter phenotype; AD_S: anther stage drought stressed shelter phenotype; T_C: tetrad stage irrigated shelter phenotype; and T_S: tetrad stage drought shelter phenotype; Figure S3. Phylogenetic UPGMA tree showing the relationships between 76 durum wheat lines genotyped with 7 wheat chromosome 4A ISBP markers. a) wheat lines are colored based on their geographic zone (Supplementary Materials Table S2) (Center: green; North: blue; North-east: light blue; North-west: dark blue; South: red; South-east: orange; South-west: maroon; and East: purple); b) wheat lines are colored for species (T. turgidum subsp durum appears in green, T. turgidum subsp turgidum in orange and T. turgidum subsp. dicoccon in purple). The colored and vertical lines indicate the differentiated clusters (cluster 1-yellow; cluster 2-blue; cluster 3-dark blue; cluster 4-pink; cluster 5-purple; cluster 6-red; cluster 7-green; cluster 8-grey; cluster 9-light blue; and black lines show the wheat lines that have not been included within any cluster); and c) each dot represents a wheat line, the colors are the same in "a)". Funding: This research was funded by project P18-RT-992 from Junta de Andalucía (Andalusian Regional Government), Spain (Co-funded by FEDER), and by the Spanish Ministry of Science and Innovation project PID2019-109089RB-C32.