*Article* **Genome Wide Association Study Uncovers the QTLome for Osmotic Adjustment and Related Drought Adaptive Traits in Durum Wheat**

**Giuseppe Emanuele Condorelli <sup>1</sup> , Maria Newcomb <sup>2</sup> , Eder Licieri Groli <sup>1</sup> , Marco Maccaferri <sup>1</sup> , Cristian Forestan <sup>1</sup> , Ebrahim Babaeian <sup>3</sup> , Markus Tuller <sup>3</sup> , Jeffrey Westcott White <sup>4</sup> , Rick Ward <sup>2</sup> , Todd Mockler <sup>5</sup> , Nadia Shakoor <sup>5</sup> and Roberto Tuberosa 1,\***


**Abstract:** Osmotic adjustment (OA) is a major component of drought resistance in crops. The genetic basis of OA in wheat and other crops remains largely unknown. In this study, 248 field-grown durum wheat elite accessions grown under well-watered conditions, underwent a progressively severe drought treatment started at heading. Leaf samples were collected at heading and 17 days later. The following traits were considered: flowering time (FT), leaf relative water content (RWC), osmotic potential (ψs), OA, chlorophyll content (SPAD), and leaf rolling (LR). The high variability (3.89-fold) in OA among drought-stressed accessions resulted in high repeatability of the trait (h <sup>2</sup> = 72.3%). Notably, a high positive correlation (r = 0.78) between OA and RWC was found under severe drought conditions. A genome-wide association study (GWAS) revealed 15 significant QTLs (Quantitative Trait Loci) for OA (global R<sup>2</sup> = 63.6%), as well as eight major QTL hotspots/clusters on chromosome arms 1BL, 2BL, 4AL, 5AL, 6AL, 6BL, and 7BS, where a higher OA capacity was positively associated with RWC and/or SPAD, and negatively with LR, indicating a beneficial effect of OA on the water status of the plant. The comparative analysis with the results of 15 previous field trials conducted under varying water regimes showed concurrent effects of five OA QTL cluster hotspots on normalized difference vegetation index (NDVI), thousand-kernel weight (TKW), and/or grain yield (GY). Gene content analysis of the cluster regions revealed the presence of several candidate genes, including bidirectional sugar transporter SWEET, rhomboid-like protein, and S-adenosyl-L-methionine-dependent methyltransferases superfamily protein, as well as DREB1. Our results support OA as a valuable proxy for marker-assisted selection (MAS) aimed at enhancing drought resistance in wheat.

**Keywords:** drought; durum wheat; osmotic adjustment; QTL

#### **1. Introduction**

Drought is one of the most devastating abiotic stressors limiting crop yield, adaptability, and quality [1,2]. Recent global climate models predict a consistent rainfall reduction in temperate drylands [3–5], hence destabilizing food systems and global food security [6]. The plant reaction to drought is mediated by complex molecular systems linked to the transcriptome [7–10], as well as hormone signaling and metabolism [11–13]. In particular, drought is the major abiotic stress curtailing yield and lowering quality [14,15] in

**Citation:** Condorelli, G.E.; Newcomb, M.; Groli, E.L.; Maccaferri, M.; Forestan, C.; Babaeian, E.; Tuller, M.; White, J.W.; Ward, R.; Mockler, T.; et al. Genome Wide Association Study Uncovers the QTLome for Osmotic Adjustment and Related Drought Adaptive Traits in Durum Wheat. *Genes* **2022**, *13*, 293. https:// doi.org/10.3390/genes13020293

Academic Editor: Patrizia Galeffi

Received: 11 December 2021 Accepted: 29 January 2022 Published: 2 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

durum wheat (*Triticum turgidum* ssp. *durum*; 2n = 28, AABB), the most cultivated wheat in the Mediterranean regions [16], whose genome sequence was recently assembled *de novo* [17]. Among the strategies adopted by plants to withstand water scarcity [18,19], osmotic adjustment (OA) plays a major role in enhancing drought resistance through an active accumulation of solutes in response to a water potential reduction, thereby preserving cellular turgor [18,20–27]. Active OA maintains relative water content at low leaf water potential in order to sustain plant growth without impairing normal cellular functions [28]. Plants accumulate low-molecular weight organic solutes, such as soluble sugars [24,29] and proline [27,30–32], both of which increase under water stress, hence enhancing OA and contributing to maintain photosynthesis, as well as stomatal conductance, at lower water potentials. To date, the dissection of the genetic basis of OA has received limited attention, mainly due to the difficulty in measuring this trait in more than a limited number of accessions [33], an essential prerequisite to properly map and characterize the effects of the QTLs [34–36] underscoring OA variability. In cereals, the dissection of the OA QTLome has been attempted in rice [37,38] and barley [39,40] based on the evaluation of biparental recombinant inbred lines (RIL) populations, hence surveying only a limited amount of genetic variability as compared to that surveyed in GWAS studies. Herein, we report the results of the first large-scale genetic dissection of the OA QTLome in wheat via GWAS based on the field evaluation of 248 durum wheat elite accessions grown under conditions of progressively increasing drought and previously tested for grain yield in 15 field trials carried out under a broad range of water regimes in Mediterranean countries [41]. Three major QTL clusters were identified, where OA was unrelated to flowering time while being positively associated with the water status of the plant and grain yield as reported in [41], supporting the beneficial role of OA in enhancing drought resistance, most likely through an avoidance strategy. A comparative analysis with the sequence information available for these regions in durum [17] and bread wheat [42] revealed a number of putative candidate genes.

#### **2. Materials and Methods**

#### *2.1. Plant Material and Field Management*

For this study, 248 durum wheat elite accessions (Durum Panel) were chosen at the University of Bologna (Table S1). Most (189) of these accessions were originally assembled by [43] to represent a large portion of the genetic diversity (Table S2) present in the major improved durum wheat gene pools adapted to Mediterranean environments. The field trial was conducted at the University of Arizona Maricopa Agricultural Center (33.070 ◦N, 111.974 ◦W, elevation 360 m) on a Casa Grande Soil (fine-loamy, mixed, superactive, hyperthermic Typic Natrargids) (Figure S1). The Durum Panel was planted on 28 November 2017 according to a row-column experimental design with two replicates. Each accession was evaluated in two-row plots (3.5-m long, 0.76-m row spacing) with an average plant density of 22 plants/m<sup>2</sup> . Orita and Tiburon, both representing the Arizona's "Desert-Durum"® wheat, were chosen as border plots. Before planting, granular nitrogen at 112 kg ha−<sup>1</sup> was incorporated into the soil. Sprinkler irrigation was used to germinate seeds and establish the crop, followed by subsurface drip irrigation matching evapotranspiration for optimal plant growth, once or twice a week as needed. The pressurized subsurface drip irrigation system was installed before planting when one dripline with emitters spaced every 0.30 m was buried at ~0.10 m depth along each seed row. The final irrigation event was on 11 March 11 2018 (i.e., 103 days after sowing, DAS), when ~50% of the accessions had flag leaf sheaths opened (i.e., at Zadoks growth stage 47) [44]. From here on, the whole experiment was subjected to a progressive water deficit until 2–3 April 2018, when plants at the anthesis halfway stage (Zadoks growth stage 65, on average) were harvested to measure total above-ground biomass.

#### *2.2. Meteorological Data and Soil Moisture Monitoring*

Daily and hourly meteorological reports for the growing season were obtained from the Arizona Meteorological Network [45]. In addition, high temporal resolution meteorological data, particularly air temperature, relative humidity, and photosynthetic photon flux density (PPFD), for the experimental site were recorded at 5-s intervals with an automated weather station (Clima Sensor US, Adolf Thies GmbH & Co. KG, Göttingen, Germany) and a quantum sensor (SQ-214, Apogee Instruments, Inc., Logan, UT, USA). These data were made available by the TERRA Phenotyping Reference Platform [46]. Vapor pressure deficit (VPD) was calculated as the difference between the saturation and actual vapor pressure [47]. The soil volumetric water content (VWC) was monitored in and between seed rows with time-domain reflectometry (TDR) sensors (True TDR-315, Acclima, Inc., Meridian, ID, USA) installed at three locations within the experiment at 1, 10, and 50 cm depths at each location. Additional soil sensors were installed between rows at 15 cm depth to measure the soil matric potential (Tensiomark, ecoTech Umwelt-Meßsysteme GmbH, Bonn, Germany). All soil sensors recorded data at 15-min intervals throughout the entire growing season. Based on the characterization of the soil hydraulic and physical properties of the experimental site under the TERRA-REF project, the volumetric water contents corresponding to the permanent wilting point (θPWP) and the field capacity (θFC) at 10–15 cm depth were approximately 0.110 and 0.282 m3/m<sup>3</sup> , respectively. The VWC dynamics at the three measurement locations, for the entire growing season, are depicted in Figure S2. The Durum Panel accessions were monitored regularly for above-ground diseases and pests, which remained below threshold levels, hence not requiring control treatments, while growing degree days (GDD) were monitored until harvest (3-4 April 2018) at 125 days after sowing (DAS) (Table S3). Growth stages of each accession were defined based on the basis of the Zadoks scale [44] at 92, 93, 98, 101, 111, and 118 DAS, and flowering time (FT) was recorded when more than 50% of ears in the plot had flowered (anthesis half-way). Plants were harvested 125 days (DAS) to allow for planting the next phenotyping experiment; therefore, biomass data indicate the status at a point in time rather than direct estimates of final yields.

#### *2.3. Evaluation of RWC, ψs, OA, LR, Leaf Chlorophyll Content (SPAD), and Biomass*

The entire Durum Panel was evaluated for leaf relative water content (RWC) and osmotic potential (ψs) in well-watered (12 March 2018, 104 DAS) and severe drought (27 March 2018, 119 DAS) conditions. At the first sampling (fully-irrigated conditions) awns were visible on approximately 50% of accessions, while the second sampling was carried out under severe drought conditions when most accessions were at early grainfilling (Figure S3). Fully expanded flag leaves of eight different plants were sampled in each plot (experimental unit) at dawn from 6:00 to 7:00 a.m. Leaves were immediately placed in sealed plastic bags, stored in portable coolers (4 ◦C) to minimize water loss due to evaporation, and transported to the lab where leaves were removed from the bags. After cutting the leaf tips (5 cm), the remaining leaf portion (average length 15 cm) was cut in the middle to obtain two homogeneous pieces of similar weight, then mixed and stored in Falcon (50 ml) conical centrifuge tubes. One batch was used to measure OA following the "Rehydration method" described in Reference [23]. Leaves were rehydrated for 4 h in distilled water to reach full turgor, then dried, and stored in a freezer (−20 ◦C). After thawing, the cell sap was collected using a garlic press, and 10 µl were transferred onto a paper sample disc covering the sampling cuvette of a vapor pressure osmometer (Wescor 5520, Logan, UT, US), previously calibrated using the 290, 1000, and 100 mmol kg−<sup>1</sup> standards. After each measurement, the osmometer cuvette was rinsed using deionized water. Finally, the resulting osmolality (mosmol kg−<sup>1</sup> ) was converted to osmolarity (MPa) using the following formula: <sup>ψ</sup>s (MPa) = <sup>−</sup>c (mosmol kg−<sup>1</sup> ) <sup>×</sup> 2.58 <sup>×</sup> <sup>10</sup>−<sup>3</sup> [48], and OA was measured as the difference between the ψs at full turgor in control and in water-stressed conditions: ψs (control)—ψs (water stress). The other batch was used to measure RWC. Fresh leaves were weighed (FW) then submerged in distilled water in the Falcon tubes

and stored at 4 ◦C for rehydration (10 h). Rehydrated leaves were wiped thoroughly with blotting paper and weighed (turgid weight: TW). Then, leaves were oven-dried at 65 ◦C for three days prior to measuring the dry weight (DW). In the end, RWC values were computed as follows: [(FW-DW)/(TW-DW)] × 100 [49] (Figure 1). Leaf rolling (LR) was visually estimated at midday (112 DAS) with a 0 (no leaf rolling) to 9 (all leaves severely rolled) score when the majority of the accessions showed a LR > 5. Finally, the chlorophyll content was assessed (114 DAP) based on Soil-Plant Analysis Development (SPAD) estimates obtained with a non-destructive chlorophyll meter SPAD-502Plus (Konica Minolta Sensing, Inc., Osaka, Japan) as an indicator of leaf chlorophyll content and nitrogen (N) status. At the end of the field trial, plants within the entire two-row plots were cut on 3–4 April 2018 with mechanical harvester (Carter Mfg. Co. equipment, Donalsonville, GA, US), while subsamples of 2–3 plants were collected to evaluate moisture content in order to estimate dry biomass.

**Figure 1.** Major "Rehydration method" steps used at the Maricopa Agricultural Center (MAC) to measure osmotic adjustment (OA) and leaf relative water content (RWC): (**A**) Sampling of eight fully expanded homogeneous flag leaves before dawn. (**B**) Stacking the eight leaves and cutting off the tips. The remaining leaf parts (ca. 15 cm long) were cut in the middle to obtain two homogeneous pieces of similar weight, then mixed and inserted in Falcon 50 ml Conical Centrifuge Tubes. (**C**) Weighing of the leaf samples for RWC. (**D**) Collection of leaf cell sap for OA analysis using a garlic press. (**E**) Calibration of the osmometer (Wescor 5520) with sodium chloride solution of increasing concentration. (**F**) Pipetting ca. 10 µl of leaf cell sap onto a paper disc placed on the sampling cuvette of the osmometer (Wescor 5520).

#### *2.4. Statistical Analysis*

The *lme4* package (*r-project*) and custom R scripts were used to conduct a spatial adjustment analysis considering row and column as random effects, as well as a moving mean of specific size. *R-project* was used to calculate repeatability values (*h 2* ) and Pearson's correlation *r* coefficients among traits. *Minitab 18* software [50] was used to calculate the global percentage of phenotypic variation (global QTL model, *R <sup>2</sup>*%) explained by all QTLs for each trait.

#### *2.5. SNP Genotyping, Population Structure, and GWAS Model*

Durum panel genomic DNA was extracted using the NucleoSpin®8/96 Plant II Core Kit from Macherey Nagel and sent for SNP genotyping to [51]. The Illumina iSelect 90K wheat SNP assay [52] was used, and genotype calls were acquired as reported in [53]. Markers were assigned on the basis of the tetraploid wheat consensus map reported in [54]. Haploview 4.2 software [55] was used to calculate Linkage Disequilibrium (LD) decay among markers for the A and B genomes, and only Single Nucleotide Polymorphisms (SNPs) with minor allele frequency (MAF) > 0.05 were considered. LD decay pattern based on the consensus genetic distances was inspected considering squared allele frequency correlation (r<sup>2</sup> ) estimates from all pairwise comparisons among intra-chromosomal SNPs in TASSEL (Trait Analysis by aSSociation, Evolution, and Linkage) *5.2.37*. The Hill and Weir

formula [56] was used in R-project to define the confidence interval (CI) for QTLs in accordance with the curve fit and the distance at which LD decays below r<sup>2</sup> = 0.3 [57]. Haploview 4.2 tagger function set to r<sup>2</sup> < 1.0 was used to calculate a kinship matrix (K) of genetic relationships among individual accessions of the Durum Panel with all non-redundant 7,723 SNPs. Kinship based on Identity-by-State (IBS) among accessions was calculated in TASSEL 5.2.37. In addition, a subset of non-redundant 2,382 SNP markers (*r <sup>2</sup>* < 0.5) was used to evaluate the population structure (Q) in STRUCTURE 2.3.4. software [58] using the corresponding tagger function in Haploview 4.2 software [55]. Numbers of hypothetical subpopulations ranging from k = 2 to 10 were assessed using 50,000 burn-in iterations, followed by 100,000 recorded Markov-Chain iterations, in five independent runs for each k in order to estimate the sampling variance (robustness) of population structure inference. Then, the rate of change in the logarithm of the probability of likelihood (LnP(D)) value between successive k values was considered (∆k statistics) [58], together with the rate of variation (decline) in number of accessions clearly attributed to subpopulations (accessions with Q membership's coefficient ≥ 0.5). Finally, the level of differentiation among subpopulations was measured using the Fixation Index (Fst) among all possible population pairwise combinations [59]. Subsequently, 17,721 SNPs with MAF > 0.05, imputed with LinkImpute (LDkNNi) [60,61], were used for GWAS-Mixed Linear Model [MLM; [62,63] in TASSEL. MLM was specified as follows: y = Xβ + Zu + e [64], where y is the phenotype value, β is the fixed effect due to the marker, and u is a vector of random effects not accounted for by the markers; *X* and *Z* are incidence matrices that related y to β and u, while e is the unobserved vector of random residual. In this study, both Kinship matrix (K) and Structure Population (Q) were included as random effects in the model (MLM-Q+K), while flowering time was included as a covariate considering GWAS QQ-plot results (Figures S4 and S5). Then, GWAS *p-*values and r<sup>2</sup> effects were analyzed, and QTL significance was determined as follows: "highly significant" refers to *p-*value < 0.0001 and "significant" refers to *p-*value < 0.001. The QTL confidence interval (CI, in cM) was measured on the basis of the average genetic distance at which LD decayed below r<sup>2</sup> of 0.3 [56], a threshold frequently adopted in GWAS [54,57,65]. Considering a LD of r<sup>2</sup> = 0.3, the corresponding inter-marker genetic distance was 3.0 cM [57], and the CI of ± 3.0 cM based on map positions of QTL tag-SNPs was chosen. Finally, Minitab 18 software [50] was used to calculate the proportion of variance for phenotypic traits explained by selected SNPs.

#### *2.6. Identification of Candidate Genes*

The physical position of each QTL was determined by the position of the flanking SNP markers after their alignment on the *Triticum turgidum* ssp. *durum* reference genome of (cv. Svevo) [66]. Genes within the confidence intervals associated with the eight main QTL hotspots were retrieved from the EnsemblPlants database [67], together with their functional annotation and the amino acid sequences of putative proteins. Gene Ontology (GO) term enrichment was determined by comparing the genes included in each QTL to the number of genes annotated in each GO term with g:Profiler web software [68]. Statistical significance of terms for genes in the physical intervals was assessed using the hypergeometric statistic for every term and the g:SCS correction method for multiple testing. Durum wheat GO annotation was retrieved from the Ensembl plant genome database. To identify the most important metabolic pathways associated to eight QTLs, genes within cluster intervals were aligned to KEGG (Kyoto Encyclopedia of Genes and Genomes) database using Reference [69]. Genes annotated within the intervals were compared with their orthologs from *Triticum aestivum* (cv. Chinese Spring; IWGSC RefSeq v1.0) [70]. Identification of candidate genes was further supported by a knowledge network (proteins, biological pathways, phenotypes, and publications) created using the KnetMiner program, using the bread wheat orthologs [71], and by the analysis of temporal and spatial gene expression at the Wheat Expression Browser and ePLANT databases [72], as of September 2021.

#### **3. Results**

#### *3.1. Population Structure*

The Durum Panel showed a clear population genetic structure with an optimal number of eight (k = 8) subpopulations on the basis of pairwise comparisons among and within subgroups with 155 accessions (62.5%) clearly grouped into one of the eight main gene pools at a Q membership coefficient ≥ 0.5, while the remaining 93 were considered as admixed. The Fixation Index (Fst) and Neighbor Joining tree [73] highlighted a high genetic diversity between the old Italian accessions (S1) and the modern French, North American, Canadian and Australian cultivars (S8), while a considerable admixture among subgroups characterized the ICARDA, CIMMYT, and Italian groups. Subgroups details are reported in Table S2.

#### *3.2. Phenotypic Analysis*

Osmotic potential in well-watered (control) conditions (ψs-c; h<sup>2</sup> = 0.57) ranged from −1.44 to −0.74 MPa, with an average of −1.13 MPa, while, in water-stressed conditions, (ψs-s; h<sup>2</sup> = 0.58) ranged from <sup>−</sup>2.63 to <sup>−</sup>1.56 MPa, with an average of <sup>−</sup>2.00 MPa. The difference between osmotic potential measured at full turgor in well-watered (control: ψs-c) and in water stressed (ψs-s) conditions was considered to compute OA (h<sup>2</sup> = 0.72), which showed a normal distribution and ranged from 0.38 to 1.48 MPa, with an average of 0.95 MPa. Figure S3 reports flowering time distribution, while Figures S6 and S7 report box plots and histogram distributions, for OA, ψs-c and ψs-s, RWC-c and RWC-s, LR, and SPAD. RWC-c (h<sup>2</sup> = 0.29) ranged from 89.9 to 101.3%, with an average of 95.7%, while RWC-s (h<sup>2</sup> = 0.78) ranged from 45.2 to 76.9%, with an average of 62.2% (Figures S6 and S7). Leaf rolling (LR; h<sup>2</sup> = 0.84) at 112 DAP ranged from 2.86 to 9.60, with an average of 6.13 (Table 1), while leaf chlorophyll content (SPAD; h<sup>2</sup> = 0.76) at 114 DAP ranged from 31.9 to 48.8, with an average of 42.0. The Pearson's correlation coefficient was positive between OA and RWC-s (r = 0.78), while a negative association was found between OA and LR (r = −0.25), RWC-s and ψs-s (r = −0.49), and RWC-s and LR (r = −0.30) (Table 2, Figure 2).

**Figure 2.** Scatter plot showing Pearson correlation between (**A**) osmotic adjustment (OA) and relative water content under drought (RWC-s), (**B**) OA and osmotic potential under drought (ψs-s), and (**C**) ψs-s and RWC-s. R-project [74].

**Table 1.** Summary statistics and heritability (%) for osmotic adjustment (OA), osmotic potential under drought stress (ψs-s), osmotic potential under well-watered conditions (ψs-c), relative water content under drought stress (RWC-s), relative water content under well-watered conditions (RWC-c), leaf rolling (LR), and chlorophyll content (SPAD) in a panel of 248 durum wheat elite advanced lines and cultivars from worldwide.


<sup>a</sup> megapascal (MPa); <sup>b</sup> % RWC.

**Table 2.** Pearson's correlation plot among osmotic adjustment (OA), osmotic potential (ψ<sup>s</sup> ) under full (--c) and deficit irrigation (--s), relative water content (RWC) under full (--c) and deficit irrigation (--s), leaf rolling (LR), and chlorophyll content (SPAD).


\*\*\* *p-*value < 0.0001, \*\* 0.0001 < *p-*value < 0.001, \* 0.001 < *p-*value < 0.01. *R-project* [74].

#### *3.3. Genetic Analysis*

The rate of linkage disequilibrium (LD) decay of the 248 durum wheat elite accessions of the Durum Panel is reported in Figure 3. The average QTL confidence interval (CI) was determined on the basis of the average genetic distance at which LD decayed below r<sup>2</sup> of 0.3 multiplied by 2, corresponding to 2.12 cM (CI = ± 1.06 cM from the QTL tagSNP). Fifteen flowering time QTLs were identified and are reported in Table S4. Major QTLs for flowering time included those on chromosome arms 2AS (*QFT.ubo-2A.1* and *QFT.ubo-2A.2*), on 4AS (*QFT.ubo-2A.1*) and 6BL (*QFT.ubo-6B.1*). Among others, *Ppd-A1* was clearly identified by *QFT.ubo-2A.2* = IWA2526. Using FT as covariate, GWAS analysis (MLM-Q+K) identified 70 significant QTLs (log *p-*value > 3.00) for ψs-c, ψs-s, OA, RWC-s, LR, and/or SPAD, organized into QTL clusters. A larger portion of ψs QTLs were detected under drought (62.5%) as compared to well-watered conditions (37.5%). In particular, two major ψs-s QTLs were observed on chromosomes 1B (*Qψsc.ubo-1B.2*) and 6A (*Qψsc.ubo-6A.1*), with a log *p-*value of 4.68 and 6.04, and with R<sup>2</sup> of 5.84 and 7.88%, respectively. A total of 15 OA QTLs were mapped on chromosome arms 1AL, 1BL, 2AS, 2AL, 2BL, 4AL, 4BS, 5AL, 6AL, 6BL, and 7BS, with the three major ones being those on chromosomes 2B (*QOA.ubo-2B.2*) and 6A (*QOA.ubo-6A.1*and *QOA.ubo-6A.2*) with a log *p-*value of 4.13, 4.01, and 4.45, and R<sup>2</sup> of 4.37, 4.23, and 4.78%, respectively (Table 3). Adopting flowering time as a covariate effectively removed the effects associated to FT on OA, except for *QFT.ubo-2A.2* = *QOA.ubo-2A.1* and *QFT.ubo-6B.1* = *QOA.ubo-6B.1*. Nine RWC-s loci were mapped on five chromosome arms (1BS, 2AS, 4AL, 6AL, and 6BL), and two major QTLs were observed on 4AL (*QRWCs.ubo-4A.2* and *QRWCs.ubo-4A.3*), with a log *p-*value of 4.83 and 4.27, and R<sup>2</sup> of 3.95 and 3.84, respectively (Table 3).

**Table 3.** Significant GWAS-QTLs for osmotic adjustment (OA) and RWC-s (*p-*value < 0.001). QTL intervals were defined based on a confidence interval of ± 3.0 cM from the map position of the QTL tagging-SNPs. The rows with grey background indicate the QTLs affecting both OA and RWC-s. Position and peak marker of each QTL region are based on the tetraploid wheat consensus map [54].


*QRWCs.ubo-6B.2* IWB71722 6B 644758469 114.3 3.00 2.44 A/G −2.46

**Figure 3.** The rate of linkage disequilibrium (LD) decay of the 248 durum wheat elite accessions (Durum Panel). The Hill and Weir formula [56] was used to describe the LD decay of *r 2* . The LD among SNPs in the Durum Panel was estimated using Haploview 4.2 [54]. The blue curve represents the model fit to LD decay (non-linear regression of *r <sup>2</sup>* on distance). A confidence interval of 2.12 cM for the QTLs is shown when LD (*r 2* ) is 0.3 (red line). R-project [74].

Nine LR loci were mapped on seven chromosome arms (1BL, 2AL, 3AS, 3AL, 3BL, 6AS, and 6BL), and one major QTL was observed on chromosome arm 2AL (*QLR.ubo-2A.1*), with a log *p-*value of 4.08, and R<sup>2</sup> of 4.92. As to SPAD, 21 QTLs were mapped on 12 chromosome

arms (1AL, 2AL, 2BS, 3BS, 3BL, 4AL, 4BS, 5AL, 6BS, 6BL, 7AL, and 7BS), and three major QTLs were observed on chromosome arms 1AL (*QSPAD.ubo-1A.1*), 4BS (*QSPAD.ubo-4B.1*), and 5AL (*QSPAD.ubo-5A.1*), with a log *p-*value of 6.08, 6.61, and 6.81, and R<sup>2</sup> of 7.87, 8.69, and 8.99, respectively. The global R<sup>2</sup> of the multiple QTL model was 58.0% for ψs-c, 56.5% for ψs-s, 63.6% for OA, 25.7% for RWC-s, 44.1% for LR, and 50.2% for SPAD. On the basis of their concurrent allelic effects on OA and other related traits, eight major QTL clusters were identified: (i) *DR\_QTL\_cluster\_1#* (RWC-s and ψs-s) on 1B, (ii) *DR\_QTL\_cluster\_2#* (OA and ψs-c) on 2B, (iii) DR\_QTL\_cluster\_3# (OA and RWC-s) on 4A, (iv) *DR\_QTL\_cluster\_4#* (OA and SPAD) on 5A, (v) *DR\_QTL\_cluster\_5#* (OA and RWC-s) on 6A, (vi) *DR\_QTL\_cluster\_6#* (OA and RWC-s) on 6B, (vii) *DR\_QTL\_cluster\_7#* (OA, RWC-s, and SPAD) on 6B, and (viii) *DR\_QTL\_cluster\_8#* (OA and SPAD) on 7B (Table 4 and Figure 4).

**Table 4.** List of GWAS-QTL clusters identified in the Durum Panel and significantly associated with osmotic adjustment (OA), RWC under drought stress (RWC-s), osmotic potential under well-watered conditions (ψs-c) and drought stress (ψs-s), and leaf rolling (LR). The co-localization with previously known normalized difference vegetation index (NDVI), chlorophyll content (SPAD), root growth angle (RGA), thousand kernel weight (TKW), and grain yield (GY) QTLs is reported.


<sup>a</sup> Reference [60], <sup>b</sup> Reference [75], <sup>c</sup> Reference [76], <sup>d</sup> Reference [41], <sup>e</sup> Reference [77] and <sup>f</sup> Reference [78].

**Figure 4.** Chromosome position on the durum consensus map [54] and R<sup>2</sup> of the QTL hotspots for osmotic adjustment (OA) and/or relative water content (RWC) on chromosome arms 2BL (*DR\_QTL\_cluster\_2#*), 4AL (*DR\_QTL\_cluster\_3#*), and 6AL (*DR\_QTL\_cluster\_5#*), as well as their overlaps with QTLs previously reported in literature: <sup>a</sup> Reference [60], <sup>b</sup> Reference [75], <sup>c</sup> Reference [76], <sup>d</sup> Reference [41], <sup>e</sup> Reference [77], and <sup>f</sup> Reference [78]. The QTLs for flowering time (FT) are shown in green.

These regions were selected for a more detailed analysis and comparative analysis with previously reported results on grain yield in both durum and bread wheat, as discussed hereafter. Durum wheat genes within the confidence intervals of the eight selected QTL hotspots were retrieved from EnsemblPlants database, together with their functional annotation (Table S5). Gene Ontology (GO) enrichment analysis (Figure 5) and KEGG pathways reconstruction (Table 5) were used to further functionally characterize the genes included in the eight QTL clusters. In parallel, the bread wheat orthologous genes were identified, as well (see Table S5), for the comparison of genes annotated in the syntenic regions of the two wheat species. KnetMiner [71] knowledge networks, constructed using bread wheat orthologs, were integrated to identify putative candidate gene(s) within the confidence interval of each QTL. The confidence interval of *DR\_QTL\_cluster\_1#* on chromosome 1B corresponds to a physical interval of approximately 7.0 Mb with 46 highconfidence (HC) genes in the Svevo genome (Table S5). Among the genes included in the interval, no GO terms were significantly enriched, while KEGG mapping annotated 21 genes to nine functional categories (Table 5). The two most notable candidates in the confidence interval are *TRITD1Bv1G127690*, which encodes a transmembrane protein with transporter activity homologous of the *Arabidopsis* Major facilitator superfamily MEE15, and *TRITD1Bv1G126800*, which encodes for a seven transmembrane MLO-like protein. The confidence interval of the *DR\_QTL\_cluster\_2#* on chromosome 2B corresponds to a 3.2 Mb interval, which contains 63 high-confidence (HC) genes. GO terms associated to "stress response" and "antioxidant activity" were enriched among these genes (Figure 5), due to the presence of 10 peroxidase encoding genes in the QTL interval. KEGG mapping confirmed their annotation in secondary metabolism pathway (Table 5), acting in the phenylpropanoid biosynthesis. In addition to these peroxidase encoding genes that could act in drought-stress response and adaptation, the two most notable candidates in the interval are *TRITD2Bv1G263980*, encoding for a protein kinase and *TRITD2Bv1G264060*, which encodes a DDB1- and CUL4-associated factor-like protein 1. Both genes are located at the confidence interval boundaries and could be functionally related to OA-related aspects (Table S5). Notably, the comparison of the syntenic physical region in *T. aestivum* Chinese Spring evidenced several gaps in the corresponding Svevo region: for 30 HC genes annotated in the Chinese Spring syntenic region, their Svevo orthologs are indeed included in unmapped scaffolds and could, thus, be included among the list of putative candidates. Even if none of these 30 missing genes are drought-stress responsive or OArelated functional annotation, the presence of these gaps in *T. durum* genome assembly and gene annotation could clearly impair candidate gene discovery. The confidence interval of *DR\_QTL\_cluster\_3#* on chromosome 4A corresponds to approximately 3.5 Mb with 33 HC genes (Table S5), mapped to seven KEGG functional categories (Table 5) and enriched in manganese transport-related GO terms (Figure 5). *TRITD4Av1G256080* and *TRITD4Av1G256120* indeed encode for membrane protein of ER body-like proteins, likely working as metal transporters. Additional genes annotated in the QTL interval include *TRITD4Av1G255460*, *TRITD4Av1G255480*, and *TRITD4Av1G255490*, encoding for three Glutathione S-transferases, *TRITD4Av1G255990*, encoding for an RNA-binding family protein with RRM/RBD/RNP motifs, and *TRITD4Av1G256200*, which encodes for a 5' methylthioadenosine/S-adenosylhomocysteine nucleosidase. In addition, for this QTL hotspot, a gap between *T. durum* and *T. aestivum* chromosome assemblies was found, hence impairing a more proper identification of candidate genes. The confidence interval of *DR\_QTL\_cluster\_4#* on chromosome 5A spans a physical interval of approximately 4.0 Mb, which contains 39 HC genes in the Svevo genome (Table S5). Lyase activity was the unique GO term enriched among the genes within the QTL (Figure 5), while KEGG mapping assigned 17 of them to eight functional groups (Table 5). Interestingly, *TRITD5Av1G246840*, encoding for a putative Phospholipase D potentially implicated in multiple plant stress responses, and *TRITD5Av1G247330*, a putative Lipoxygenase required for jasmonic acid accumulation, were mapped to the "environmental information processing" category. In addition, *TRITD5Av1G247220*, encoding for an UDP-N-acetylglucosamine (UAA) trans-

porter family protein, could be involved in the osmosensory signaling pathway and cell wall organization. The confidence interval of the *DR\_QTL\_cluster\_5#* on chromosome 6A spans approximately 2.0 Mb, with 32 HC genes (Table S5). While the "manganese binding" function was the only GO-enriched term (Figure 5), six of the genes included in the region encode for putative Cinnamoyl CoA reductases and were mapped to the "secondary metabolism"/"phenylpropanoid biosynthesis" pathways (Table 5). The two most notable candidates in the interval are *TRITD6Av1G217800* and *TRITD6Av1G218080*, which encode for two F-box protein PP2, and *TRITD6Av1G217670*, encoding for DREB1, a CRT-binding factor. The confidence interval of *DR\_QTL\_cluster\_6#* on chromosome 6B spans a 5.9 Mb interval, which contains 40 HC genes in the Svevo genome (Table S5). No GO term was found significantly enriched, while KEGG mapping annotated 20 genes to 12 different functional categories (Table 5). In particular, three genes encoding for MYB transcription factors potentially regulating different aspects of stress response were assigned to the "genetic information processing" category. Moreover, *TRITD6Bv1G133070*, orthologs of the CCT motif-containing response regulator protein coding gene of *Arabidopsis*, appears as an even more interesting candidate gene at this locus. The confidence interval of the *DR\_QTL\_cluster\_7#,* located on chromosome 6B, spans approximately 1.2 Mb, with 14 HC genes, including *TRITD6Bv1G207930*, encoding for a protein kinase family protein/WD-40 repeat family protein 3. Both GO-enrichment and KEGG mapping did not identify other genes in this QTL interval (Table 5 and Figure 5). Finally, the confidence interval of *DR\_QTL\_cluster\_8#* on chromosome 7B corresponds to a physical interval of approximately 1.5 Mb, which contains 20 HC genes in the Svevo genome, lacking GO enrichment or predominant KEGG annotations (Table 5 and Figure 5), but including *TRITD7Bv1G002000*, a gene encoding for a photosynthetic NDH subcomplex B3 (Table S5).

**Figure 5.** Gene Ontology (GO) enrichment analysis. Dot plot shows GO terms of biological processes (BP), molecular functions (MF) and cellular compartment (CC) identified using g:Profiler [68] to be enriched (adjusted *p-*value < 0.05) among the genes included in each QTL interval. The size of the dots is based on gene count enriched in the pathway, and the color of the dots represents the adjusted *p-*values.


**Table 5.** Summary of KEGG functional pathways mapped for genes included in each QTL interval, grouped based on metabolic activities.

#### **4. Discussion**

A number of authors have proposed OA as an important adaptive mechanism to support higher crop yield under stressful environmental conditions, as reviewed in [24,79]. Notably, grain yield differences have been shown to be positively correlated to OA in cereals [21,30,80–83], hence indicating OA as a valuable proxy to predict grain production [24]. This notwithstanding, the QTLome regulating OA in wheat and other crops remains basically unknown, the main reason being the difficulty to adequately screen the large number (>200) of (i) RILs of the mapping populations and/or (ii) accessions of GWAS mapping panels required for a meaningful QTL discovery. In field conditions, the collection of leaves and their processing must be completed rapidly to minimize the bias introduced by the time of sample collection in an adequately large number of genotypes, an essential prerequisite for identifying and accurately mapping QTLs [19,34]. The QTLome dissection of OA in cereals was first attempted in rice [37,38] and barley [39,40]. In bread wheat, Reference [84] mapped an osmoregulation gene locus [85] located on the short arm on chromosome 7A. However, OA and osmoregulation differ and have different functional meanings. While OA refers to a lowering of osmotic potential (ψs) due to an accumulation of osmolytes in response to a water deficit, osmoregulation refers to the ψs regulation by the addition/removal of osmolytes until the intracellular potential is approximately equal to that of the medium surrounding the cell [20]. The gene described by Morgan regulates turgor pressure and water content by osmotic adjustment [84,85], hence the term osmoregulation. In this study, OA was measured according to the "rehydration method" [23,86]. Although this method was criticized by [85], other authors consider it optimal for screening large populations [20,26,86–89] in view of its merits in terms of labor and cost-effectiveness as compared to the other methods [23]. In our experience, the rehydration of the leaf samples greatly facilitated (i) the cell sap extraction especially in samples collected in water-stressed plants and (ii) the OA screening of the 248 diverse accessions of the Durum Panel. Collectively, this resulted in high OA repeatability (h<sup>2</sup> ) and a positive and negative correlation with RWC-s and LR, respectively. The positive correlation between OA

and RWC clearly indicates an active physiological role of the former to maintain a more favorable water status of the plant, playing a key role for avoiding and mitigating the negative effects of water loss under drought. Overall, our results validate the effectiveness of the "Rehydration method" as an ideal option for handling the large number of samples required for the genetic dissection of the OA QTLome.

#### *4.1. GWAS Mapping and Comparative Analysis with Previous QTL Studies in Durum Wheat*

Overall, eight major QTL hotspots were detected on the long arm of chromosomes, 1BL, 2BL, 4AL, 5AL, 6AL, 6BL, and 7BS. The use of flowering time (FT) as covariate for the GWAS analysis reduced the bias caused by the photoperiod-response (*Ppd*) locus and other loci that affect FT, hence allowing a more accurate estimate of QTL effects on a per se basis rather than due to effects related to variability in phenology. Notably, none of the eight major QTL hotspots evidenced by GWAS analysis overlapped with the osmoregulation gene locus described by Reference [85] in bread wheat. *DR\_QTL\_cluster\_1#*, *DR\_QTL\_cluster\_2#*, and *DR\_QTL\_cluster\_5#* overlapped with Normalized Difference Vegetation Index (NDVI) loci identified in 2017 on the same Durum Panel under similar drought conditions using Unmanned Aerial Vehicles (UAV-Sequoia and UAV-Red-Edge), as well as ground-based platforms [60]. Additionally, *DR\_QTL\_cluster\_3#* and *DR\_QTL\_cluster\_5#* overlapped with chlorophyll content (SPAD) loci under drought described in Reference [60]. Both NDVI and SPAD have long been recognized for their ability to estimate crop biomass and predict grain yield [90–94]. *DR\_QTL\_cluster\_2#* overlapped with grain yield, thousand-kernel weight, and NDVI loci previously reported in a durum wheat elite population tested in contrasting thermo-pluviometric conditions [76]. *DR\_QTL\_cluster\_5#* co-mapped with *QRga.ubo-6A.2*, one of the most important loci for root growth angle in durum wheat [75], with thousand-grain weight, particularly under low water availability environments, as well as with grain yield, in the 183 elite accessions of the Durum Panel that were previously evaluated in 15 field trails under a wide range of water regimes [41]. *DR\_QTL\_cluster\_3#* and *DR\_QTL\_cluster\_7#* co-mapped with a major grain yield QTL reported by Reference [78] in an RIL population developed from the hexaploid wheat cross between Chinese Spring × SQ1 evaluated across a broad combination of 24 site × treatment × year combinations. The concurrent effects on grain-yield related traits reported herein fully support the conclusions of Reference [24] on OA being a valuable proxy with a positive effect on crop yield under water-limited conditions and not merely for survival under severe drought. These QTL hotspots will further enhance drought tolerance in durum wheat.

#### *4.2. Candidate Genes*

By combining the physical confidence interval position of the QTL hotspots, functional prediction of annotated genes, and biological data mining, we investigated candidate genes from selected QTLs involved in OA and/or drought resistance. Putative drought candidate genes encoding for proteins involved in drought stress responses, as well as grain development, were mapped within the eight major selected QTL hotspots, even if some highlighted gaps in the Svevo genome assembly could hinder candidate gene identification at two QTL hotspots. Among the identified candidates, the seven transmembrane MLO-like protein (*TRITD1Bv1G126800*; *DR\_QTL\_cluster\_1#*) was shown to act in drought and salt stress responses through signaling of the phytohormone abscisic acid (ABA) [95], with biological knowledge networks analysis strongly supporting its role in oxidative stresses, salt, and drought tolerance (Figure 6A). Similarly, the phospholipase D (PLD) coding gene (*DR\_QTL\_cluster\_4#*) is also involved in ABA responses [96]. Interestingly, genes associated with manganese transport and binding were identified in both *DR\_QTL\_cluster\_3#* and *DR\_QTL\_cluster\_5#*. Moreover, exogenous application of Mn was recently shown to reduce the negative effects caused by drought, harsh temperature, and salinity, increasing ROS detoxification and secondary metabolism [97,98]. Despite the clear enrichment for stress response-related GO terms among genes at *DR\_QTL\_cluster\_2#*, the identification of 10 tandemly duplicated peroxidase encoding genes acting in the phenylpropanoid biosynthetic make the construction of a knowledge network and putative candidate gene more complex and less reliable. Conversely, of considerable interest is, instead, the dehydrationresponsive element-binding protein DREB (*DR\_QTL\_cluster\_5#*; Figure 6B), that belongs to a family of plant-specific transcription factors that can specifically bind to DRE/CRT elements in the response to abiotic stresses, such as drought, salt, and low temperature [99,100], reviewed in [101]. In addition, *TRITD6Bv1G133070* (*DR\_QTL\_cluster\_6#*), orthologs of the *Arabidopsis* CCT motif-containing response regulator protein, was shown to be involved in both the generation of circadian rhythms and long-term drought adaptation [102]. Finally, *TRITD7Bv1G002000* (*DR\_QTL\_cluster\_8#*, affecting both OA and SPAD) encodes for chloroplast NAD(P)H dehydrogenase complex, involved in cyclic electron flow around photosystem I to produce ATP [103].

**Figure 6.** KnetMiner network views displaying knowledge networks of selected candidate genes Table S5. *TRITD1Bv1G126800* (**A**), encoding for a seven transmembrane MLO-like protein (*DR\_QTL\_cluster\_1#*), and *TRITD6Av1G217670* (**B**), which encodes for a dehydration-responsive element-binding protein DREB (*DR\_QTL\_cluster\_5#*). Networks were constructed using the *Triticum aestivum* orthologous genes *TraesCS1B02G216300* and *TraesCS6A02G381200*, respectively, and can be accessed using the following links: https://knetminer.com/beta/knetspace/network/528cbd3a-52 d5-40b5-91be-f59323db55a3 (accessed on 1 April 2021) and https://knetminer.com/beta/knetspace/ network/44b31582-bfaa-495f-9272-87a4f06d40a6 (accessed on 1 April 2021).

#### **5. Conclusions**

This study is the first to report QTLs for OA via GWAS mapping in wheat. From a methodological standpoint, our results support the validity of the "Rehydration method" as the fastest and most effective protocol for large-scale screening of OA under wellwatered and drought conditions. The genetic variants within the Durum Panel evaluated herein allowed for the detection of significant loci for OA, ψs, RWC, LR, and SPAD, with eight multiple concurrent QTL hotspots, all unrelated to phenology, hence being more valuable from a breeding standpoint. Importantly, five of these clusters (*DR\_QTL\_cluster\_1#*,*DR\_QTL\_cluster\_2#*, *DR\_QTL\_cluster\_3#*, *DR\_QTL\_cluster\_5#*, and *DR\_QTL \_cluster \_7#*) co-located with QTLs for grain yield and/or grain yield-related traits described in previous multi-environmental studies and, in one case (*DR\_QTL\_cluster\_5#*), co-located with a major QTL controlling root growth angle which has been demonstrated to play a relevant role in maintaining the water access in deep soils during wheat terminal drought [104]. The candidate genes identified within the confidence intervals of selected drought response-specific QTL hotspots provide useful markers for future breeding schemes. Our results support the role of OA as an important drought-stress adaptive mechanism with beneficial effects on the plant water status in durum wheat.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/genes13020293/s1, Figure S1: field trial experiment, Figure S2: volumetric water content analysis, Figure S3: Zadoks growth stage intervals & flowering time, Figure S4: Q-Q (quantilequantile) of the GWAS analysis for osmotic adjustment (OA), Figure S5: Manhattan plot of the GWAS analysis for osmotic adjustment (OA), Figure S6: box plots for (A) relative water content (RWC) and (B) osmotic potential (ψs) under full and deficit irrigation, Figure S7: histograms, Table S1: list of the 248 durum wheat accessions, Table S2: genetic population structure of the Durum Panel through model-based clustering in STRUCTURE 2.3.4 software, Table S3: growing degree days, Table S4: list of QTLs for flowering time, Table S5: candidate genes.

**Author Contributions:** E.B., J.W.W., M.M., M.N., M.T., N.S., R.T., R.W. and T.M. designed the experiment. M.N., E.B., J.W.W. and R.W. managed the field experiment, G.E.C., E.L.G. and M.N. collected the leaf samples in the field-grown durum panel under well-irrigated and water-stressed conditions. G.E.C. and E.L.G. conducted the laboratory measurements. G.E.C., E.L.G. and M.M. produced the genotypic data. C.F. and M.M. conducted the candidate gene analysis. G.E.C., C.F., M.M. and R.T. analyzed the data, interpreted the results, and developed the first draft. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study is part of the TERRA REF experiment, funded by the Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy, under Award Number DE-AR0000594. This study was supported in part by USDA-ARS agreement no. 58-2020-6-019 from the USDA Agricultural Research Service.

**Institutional Review Board Statement:** Only Plants, No animals or humans were tested, so it can be eliminated.

**Informed Consent Statement:** Only Plants, No animals or humans were tested, so it can be eliminated.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This study is the result of the collaborative project among the Department of Agricultural and Food Sciences (DISTAL) of the University of Bologna (Italy), Maricopa Agricultural Center (MAC) of the University of Arizona (USA), and the U.S. Arid-Land Agricultural Research Center of USDA ARS (USA). The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. The authors are thankful for all the support and technical assistance provided by John T. Heun, Pedro Andrade-Sanchez, and Zbigniew S. Kolber. The authors are grateful to the late Abraham Blum for his encouragement to pursue this study and for his valuable suggestions.

**Conflicts of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

#### **Abbreviations**

CI, confidence interval; DAS, days after sowing; DW, dry weight; Fst, fixation index; FT, flowering time; FW, fresh weight; GDD, growing degree days; GO, Gene Ontology; GWAS, genome-wide association study; GY, grain yield; K, kinship matrix; LD, linkage disequilibrium; LR, leaf rolling; MAF, minor allele frequency; MAS, marker-assisted selection; MLM, mixed linear model; NDVI, normalized difference vegetation index; OA, osmotic adjustment; PPFD, photosynthetic photon flux density; Q, structure population; QTL, quantitative trait locus, RGA, root growth angle; RIL, recombinant inbred lines; RWC, relative water content; SNP, single nucleotide polymorphism; SPAD, soil plant analysis development; TDR, time-domain re-flectometry; TKW, thousand kernel weight; TW, turgid weight; VPD, vapor pressure deficit; VWC, volumetric water content; ψs, osmotic potential; θPWP, permanent wilting point; θFC, field capacity.

#### **References**

