Genome-Wide Association of Stem Carbohydrate Accumulation and Remobilization during Grain Growth in Bread Wheat (Triticum aestivum L.) in Mediterranean Environments

Water deficit represents an important challenge for wheat production in many regions of the world. Accumulation and remobilization of water-soluble carbohydrates (WSCs) in stems are part of the physiological responses regulated by plants to cope with water stress and, in turn, determine grain yield (GY). The genetic mechanisms underlying the variation in WSC are only partially understood. In this study, we aimed to identify Single Nucleotide Polymorphism (SNP) markers that account for variation in a suite of WSC and GY, evaluated in 225 cultivars and advanced lines of spring wheat. These genotypes were established in two sites in the Mediterranean region of Central Chile, under water-limited and full irrigation conditions, and assessed in two growing seasons, namely anthesis and maturity growth periods. A genome-wide association study (GWAS) was performed by using 3243 SNP markers. Genetic variance accounted for 5 to 52% of phenotypic variation of the assessed traits. A rapid linkage disequilibrium decay was observed across chromosomes (r2 ≤ 0.2 at 2.52 kbp). Marker-trait association tests identified 96 SNPs related to stem weight (SW), WSCs, and GY, among other traits, at the different sites, growing seasons, and growth periods. The percentage of SNPs that were part of the gene-coding regions was 34%. Most of these genes are involved in the defensive response to drought and biotic stress. A complimentary analysis detected significant effects of different haplotypes on WSC and SW, in anthesis and maturity. Our results evidence both genetic and environmental influence on WSC dynamics in spring wheat. At the same time, they provide a series of markers suitable for supporting assisted selection approaches and functional characterization of genes.


Introduction
Bread wheat (Triticum aestivum L.) is one of the world's most important cereal crops in terms of cultivated area and economic activity [1]. It is expected that global demand for bread wheat will increase by 70%, by 2050, in line with the projected growth in world population [2]. Drought and high temperatures are the most important abiotic factors affecting nearly 200 million ha of global wheat production [3]. Wheat is very sensitive to climatic and environmental conditions [4], which is reflected in different yields when the same genotype is grown under water-limited and full irrigation conditions [5].
In Mediterranean-climate regions, temperate cereals like wheat are usually exposed to a severe water deficit during the grain filling period, the so-called terminal droughtstress [6][7][8][9]. This water deficit leads to reductions in canopy photosynthesis, as well as lower levels of assimilate transfer to the grain, leading to smaller kernel weights and GY [7,10].
Water-soluble carbohydrates (WSCs) are accumulated in stems during pre-anthesis (and up to 20% of them after anthesis in wheat) as reserves, and later these are remobilized to the grain [11][12][13], reducing the negative effects of water deficit. Indeed, the contribution of stem WSCs to grain growth could be very significant under water-deficit conditions [14]. A wide genetic variability in stem WSCs has been reported for wheat. Considering that variability, high and low WSC genotypes with differential yield, stem biomass, or yield tolerance index have been identified, utilizing recombinant inbred lines [12], cultivars and advanced semi-dwarf lines [13,15,16], lines deriving from synthetic hexaploids [14], and double-haploids populations [17]. The variation in total WSC has been attributed to genotypic variation in components of the fructan biosynthesis [16,18].
Physiological traits associated with tolerance to water deficit, such as accumulation and remobilization of stem WSCs, are included in some breeding programs, to increase yield in rainfed environments [19,20]. These traits are quantitative, determined by multigenic control, and strongly influenced by the environment [21]. The "dissection" of complex traits has been enhanced by the development of analysis technologies focused on identifying quantitative trait loci (QTL) underlying phenotypic variation [22][23][24][25]. Analysis of candidate genes and genome-wide association studies (GWASs) have been performed in bread wheat to identify SNP markers controlling traits such as dehydration tolerance [26,27], leaf rust resistance [28], grain yield [29], carbon isotope discrimination [30], kernel and root morphology [31,32], resistance to low nitrogen [33], and leaf traits [34], among others. For instance, a GWAS analysis of stem WSCs at their maximum concentration (~14 days post-anthesis) was performed in 166 Chinese bread wheat cultivars (grown in four environments), using 18,207 SNP markers, and this identified 52 significant marker-trait associations on all wheat chromosomes except for 2A, 2D, 4D, 5B, 6A, and 6D [35]. Another study conducted in a large collection (n = 990) of bread wheat genotypes in Australia revealed significant associations between stem WSCs at post-anthesis and molecular markers (SNPs and diversity arrays technology (DArT)) located on chromosomes 1A, 1B, 1D, 2D, and 4A [36]. Recently, Fu et al. [37] conducted a GWAS on 166 wheat accessions, from China and other countries, utilizing SNP arrays. They identified 19-36 significant loci associated with WSC (measured by near-infrared spectroscopy) at 10 to 30 days post-anthesis. These loci were distributed across the seven groups of chromosomes.
In the present study, we performed a GWAS and a haplotype analysis, to identify SNP markers accounting for variation in a suite of stem traits, including weight, the WSC content at anthesis and maturity, and the apparent remobilization of stem WSCs. We analyzed 225 advanced spring wheat lines, which were established at two sites in the Mediterranean-climate region of Central Chile, under water-limited and full irrigation conditions, respectively. These wheat genotypes were obtained from different breeding programs developed by research institutions in Chile, Uruguay, and Mexico. We hypothesized that, in a diverse set of accessions, a high-density genotyping with SNP markers and the linkage disequilibrium present in specific regions of the genome would allow identification of new genome-phenotype associations for stem carbohydrate traits related to tolerance to water deficit. This study complements our previous research [15,16,30] on the genetic and genomic characterization of physiological and agronomical traits of bread wheat under water deficit. At the same, considering the importance of the genotype by environment interaction on wheat yield and adaption, this study aimed to broaden our current knowledge about the marker-trait associations, particularly under specific conditions representative of the South American Mediterranean climate. Our study analyzed that sort of associations, integrating two analytical approaches, single-marker and haplotype association tests. This topic has been scarcely described for other studies dealing with WSC and yield traits in wheat. The SNP-markers identified by these approaches contribute to defining candidate genes suitable for improving, by traditional or biotechnological breeding methods, both adaptive and agronomic traits in wheat. This will be required to keep its productivity under more challenging environmental conditions.

Grain Yield, Stem Carbohydrate Accumulation, and Remobilization Were Influenced by Genotype and Site
We compared the 225 wheat cultivars and advanced lines grown at two sites during 2012, to analyze the effect of contrasting water regimes on genotype performance. These sites, Cauquenes and Santa Rosa, represented water-limited (rainfed) and full irrigated conditions, respectively. Large differences in GY were observed between both sites, with an average of 3140 ± 550 and 9820 ± 970 kg ha −1 at Cauquenes and Santa Rosa, respectively ( Figure 1a). The opposite pattern was observed for SWa, WSCa, and ARWSC (Figure 1b-f), indicating an increase in these traits in response to a lower water availability (data distributions for the other traits are in Supplementary Materials Figure S1). Mean values for the different traits in each site conditions are given in Supplementary Materials Table S1. Analysis of variance revealed a significant genotype by site effect for GY, SWm, WSCm, WSCCm and WSCCa/m 2 ( Table 1). For the other traits, the genotype (G) and site (S) main factors were both significant. The relative importance of the genotype by site (GxS) interaction was reflected in the percentage of phenotypic variance that this factor accounted for, ranging from 12.4 to 18.3% for SWm and WSCCa/m 2 , respectively. The genetic variance, estimated by σ 2 G , ranged from 5.0 to 51.6% of the phenotypic variation for WSCm and SWa, respectively.

Significant SNPs Were Detected in Association with Stem Carbohydrate Accumulation and Remobilization
The single-SNP tests identified ninety-six significant associations across all sites ( Table 2). The total number of significant SNPs varied according to the trait analyzed, ranging from 2 (for SWa and SWm) to 23 (for WSCa), but it also varied depending on the site, from 5 (Cauquenes 2011) to 36 (Santa Rosa FI 2011). The identified SNPs were specific for each site and year, as is depicted in Supplementary Materials Figure S2. Significant SNPs were distributed across almost all the wheat genome, with a higher representation of 6B and 6D chromosomes, which had 10 and 14 SNPs, respectively. Details of the significant SNPs per chromosome are provided in Supplementary Materials Table S2. The highest proportion (64%) of significantly associated SNPs was located within intergenic regions, whereas nearly a third of them (34%) were part of gene coding regions ( Figure 2a); for this last category, 37% encoded proteins without a characterized function and 63% were proteins participating in energy/metabolism (18%), protein syn-thesis/modification (12%), transcription (12%), and signal transduction, among others ( Figure 2b). A list of SNP-markers significantly associated (p-value < 0.001 and q-value < 0.15) with the assessed traits is presented in Table 3 (a complete list is included in Supplementary Materials Table S3). For WSCa in Cauquenes, in 2012, the detected SNPs belonged to genes encoding a serine/threonine-protein kinase, a guanosine nucleotide diphosphate dissociation inhibitor, a strictosidine synthase, an ABC transporter C family member 10, and a hypersensitive induced response protein 3. These proteins are part of functional categories, such as protein synthesis or modification, energy or metabolism, cell transport, and pathogen defense. Other SNPs, located in genes encoding uncharacterized proteins, were also identified for WSCa in Cauquenes, 2012, and ARWSC in Santa Rosa, 2011. Finally, for ARWSC and WSCCa/m 2 in Santa Rosa, 2011, some significant SNPs were also located in intergenic regions. In terms of the specific effect of each SNP on trait variation (R 2 ), it ranged from 8.4% (0.0837) to 14.9% (0.1491).

A Rapid Linkage Disequilibrium (LD) Decay Was Observed across Chromosomes
The extent of linkage disequilibrium (LD) was analyzed across each chromosome. On average, LD over physical distance decayed to r 2 ≤ 0.2 at 2.52 kbp. The variation in the decay depended on the specific chromosomes, with the most rapid decay observed on chromosome 2D (r 2 ≤ 0.2 at 1.09 kbp) and the slowest decay on chromosome 6D (r 2 ≤ 0.2 at 8.31 kbp). Graphs showing the LD decay per chromosome are included in Supplementary Materials Figure S3. The LD observed for two representative genomic regions containing significant SNPs is depicted in the heatmaps in Figure 3. The first corresponds to an intergenic region (1870 bp) located on chromosome 6D. This region contained the markers iniaGBS52443 (SNP1382) and iniaGBS52444 (SNP1383), which were significantly associated with WSCCa and WSCCa/m 2 ( Figure 3a, Table 3, and Supplementary Materials Table S3). The second region (768 bp) is on chromosome 3A, and it encodes a serine/threonine-protein kinase. A marker within this gene, iniaGBS56801 (SNP1507) was associated with WSCa ( Figure 3b). These SNPs, in an intergenic and gene region, were part of moderate-to-high LD and low LD regions, respectively.

Single-SNP Markers Are Part of Significant Haplotypes in Association with Traits
As an insight to determine the effect of haplotypes on the studied traits, we carried out an SNP-haplotype analysis for three gene regions harboring significant single-SNP markers. For the serine/threonine-protein kinase gene (chromosome 3A), we identified four significant haplotypes that were associated with SWa and another three associated with SWm, WSCCm, and WSCCa in Cauquenes, 2012 ( Table 4). Analysis of the intergenic region on chromosome 6D identified three significant haplotypes for SWa, three for SWm, and one for WSCCm/m 2 , in Santa Rosa, 2011 (full irrigation). Finally, for the gene encoding an endoglucanase, present on chromosome 7B, we identified two significant haplotypes in association with WSCa and WSCCa/m 2 , in Santa Rosa, 2011 (mild water stress). From these results, we further considered SWa as a study case, to analyze the specific effect of haplotypes on phenotype. In Santa Rosa, 2011 (full irrigation), haplotypes H1 (ACACCAGGGCCT) and H3 (CCGTTAGGGTTC), belonging to the intergenic region, differed significantly in terms of their mean values for SWa (Figure 4a). Similarly, for Cauquenes 2012, the H1 (GTGTTCAT) and H4 (GTAGCTGT) haplotypes in the serine/threonine-protein kinase gene had a significantly lower mean SWa than H2 (GTAGCTAC) and H3 (TTAGCCAT) (Figure 4b).  Figure 4).   Figure 4).

Discussion
The GWAS analysis of GY and nine stem WSC traits was performed in a panel of 225 advanced wheat lines and cultivars grown in environments subjected to terminaldrought-stress and full-irrigation conditions during two growing seasons. Although the selected genotypes had a relatively narrow range of heading dates (from 79 to 84 days), a significant effect of genotype, site, or GxS interaction was observed on all the analyzed traits (Table 1). These effects confirmed the high phenotypic variability for physiological and agronomic traits estimated previously, when a larger set of genotypes (378 lines, including those included in this study) was analyzed [15]. Furthermore, in an elite wheat-breeding population assessed in Australia [38], genotype and site, either as principal or combined (GxS) factors, were also shown to be significant for WSCC and WSCC/m 2 (at approximately 180 degree-days post-anthesis) across environments differing in water availability. In the current work, the differential influence of genetics relative to other factors (including environmental ones) was represented by a low-to-moderate σ 2 G /σ 2 P ratio, which ranged from 0.05 to 0.52 for WSCm and SWa, respectively (Table 1). These estimates are lower than others mentioned by Dong et al. [35] for WSC, indicating that this trait is largely genetically determined (broad-sense heritability of 0.78-0.90). However, distinct genotypes, experimental conditions, genetic analyses, and methods for measuring WSC would explain the observed difference between studies.
Genotyping and analysis of LD detected a relatively rapid decay in this parameter over physical distance (r 2 ≤ 0.2 at 2.52 kbp). A similar pattern was reported by Sela et al. [39] for emmer wheat (Triticum dicoccoides), which decayed from r 2 = 0.5 to r 2 = 0.1 within a 2 kbp span. These authors suggested a significant effect of recombination operating as a diversifying force. The breakdown in LD is primarily driven by recombination, so the extent of LD is expected to vary as an inverse relation to the local recombination rate [40]. In this study, we considered a density of 0.22 SNPs per Mbp (Supplementary Materials  Table S4), after SNP filtering. Considering an average of 14.6 SNPs per Mbp estimated in a recently published wheat genome sequence [41], further studies should be able to confirm our estimation for LD decay.
The GWAS analysis using 3243 SNP-markers identified 96 significant associations for the assessed traits (Table 2). Significant single-SNPs were specific to each site, water regime, and growing season (Supplementary Materials Figure S2), showing an experiment-specific pattern of results, and the presence of an overall GxS and genotype by time (GxT) effect (be-tween growing seasons). A similar significant involvement of GxS in phenotype variability in wheat was identified for WSC and WSCC in experiments conducted in China [35] and Australia [36], respectively. In our study, the number of significant associations depended on the trait (Table 2) and the chromosome (Supplementary Materials Table S2). The contribution of particular chromosomes within the wheat genome towards yield (e.g., GY) and WSC traits is variable within the literature describing GWAS results from trials established on sites with water limitations [30,35,37,42,43]. Indeed, most of the wheat chromosomes were represented across these studies, and the effect of individual significant SNPs on phenotypic variation, expressed in terms of R 2 , was spread across a wide range of values (up to 16%).
Our association tests showed that most of the significant SNPs (64%) were within intergenic regions (Figure 2). This is consistent with the overall distribution of SNPs determined from information contained in public databases (e.g., National Center for Biotechnology  Information, NCBI), where about 60% of SNPs are in genomic regions outside coding regions [44]. Such polymorphisms could be part of cis-regulatory elements and modulate gene expression at a transcriptional level. In contrast, 34% of significantly associated SNPs were located within genes which, for example, encode proteins such as serine/threonine protein kinase, strictosidine synthase, ABC transporter, and hypersensitive-induced response protein, among others. These types of proteins have been implicated in responses to stress, membrane transport, signal transduction, and cell death. Other GWAS analyses of wheat WSC traits have identified candidate genes from two main categories: carbohydrate metabolism and stress responses [35,37]. Differences among studies are expected due to the diversity of experimental factors and the genotyping approaches implemented. Previously, we studied the stem carbohydrate dynamics and expression of genes involved in fructan accumulation and remobilization during grain growth [16]. In specific, we detected significant differences (at the transcript level) between tolerant and susceptible genotypes, under water stress, in the regulation of genes encoding fructan biosynthesis enzymes, such as fructan 1-fructosyltransferase (1-FFT), fructan:fructan 6G-fructosyltransferase (6G-FFT), sucrose:fructan 6-fructosyltransferase (6-SFT), fructan 1-exohydrolase w2 (1-FEH w2), and fructan 1-exohydrolase (1-FEH). These genes have been mapped in different chromosomes. The 1-FFT-A1 gene was located on chromosome 4A [45]; the 1-FEH w1, w2, and w3 genes were mapped on chromosomes 6A, 6B, and 6D [46]; and 1-SST genes were isolated and located on chromosomes 4A, 7A, and 7D [47] of Chinese bread wheat cultivars. However, after the quality control performed on genotyping data, no SNPs from these genes were represented in the datasets utilized in our association tests. The use of genotyping 90 K arrays has recently confirmed the presence of significantly associated SNPs, belonging to those genes, validating the importance of genes regulating fructan metabolism in the response of wheat to water stress [37].
From single-SNP associations (Table 3), a set of genes was identified as candidates for further functional characterization. The first of these is a gene that encodes a serine/threonine protein kinase. This enzyme has important functions in abiotic and biotic stress-signaling pathways [48]. In wheat, this type of protein has been identified in response to drought [49]. Members of the wheat SnRK gene family (e.g., TaSnRK2.4) have been implicated in the regulation of enhanced osmotic potential, strengthening tolerance to drought, salt, and freezing stresses [50]. The second gene encodes a guanosine nucleotide diphosphate dissociation inhibitor (GDI), a protein interacting with Rab GTPases, which mediate the targeting and fusion of vesicles in the exocytic and endocytic pathways [51]. In plants, GDI has been related to regulation of the growth of pollen tubes and root hairs, as well as tolerance to oxygen deprivation [52]. The third gene encodes a strictosidine synthase (STR), an enzyme that plays an important role in the biosynthesis of terpenoid indole alkaloids, and in rice it is important for anther development and pollen-wall formation [53]. In Catharanthus roseus, STR is involved in the response to salinity, low temperature stress [54], and PEG-induced drought stress [55]. The fourth gene encodes an ABC transporter C family member 10, a protein that uses ATP to transport substances across the cell membrane. In wheat, this class of transporter has been localized in the vacuolar membrane and is involved in grain formation and mycotoxin tolerance [56]. The fifth gene encodes the hypersensitive induced response protein 3, which is part of the hypersensitive reaction in plants. It is related to the development of spontaneous lesions in leaves in defense against pathogen attack, and specifically, this protein triggers cell death in maize, barley, and rice [57]. Finally, the last of these candidate genes encodes an endoglucanase belonging to pathogenesis-related (PR) defense proteins, which are accumulated in response to viruses, bacteria, and fungi, as well as to abiotic factors. In wheat, in particular, endoglucanase is produced in response to drought [58].
We performed an SNP-haplotype analysis for three gene regions harboring significant single-SNP markers. This approach allowed for the detection of wider genomic regions associated with the different phenotypes analyzed (Figures 3 and 4 and Table 4). Haplotype-based analysis has been proposed as one of the most important developments in the genetic analysis of crop genomes [59]. Among other features, haplotype-analysis could capture epistatic interactions between SNPs at a locus, reduce the number of tests, and provide more power than a single marker when an allelic series exists at a locus [60]. In our study, significant haplotypes were identified in the genome in both coding and selected intergenic regions. At the same time, the haplotypes included SNPs that were identified by single-marker associations, evidencing the consistency between both approaches. Information derived from both of the detection techniques allowed for the definition of regions with high LD that were significantly associated with WSC traits, and this indicated significant differences in the performance of genotypes containing specific haplotypes (Figure 4). A haplotype-based approach also can be an alternative for detecting genometrait associations when detection by single-SNP analysis produces a limited number of significant associations.

Association Population and Growth Conditions
A collection of 225 different genotypes of spring wheat, generated by breeding programs in Uruguay (Instituto Nacional de Investigación Agropecuaria, INIA), Chile (Instituto de Investigaciones Agropecuarias, INIA), and Mexico (International Wheat and Maize Improvement Centre, CIMMYT), was utilized as an association population. This set is part of the 384 lines described and analyzed previously by del Pozo [15], in terms of their performance under water-limited conditions in a Mediterranean climate. The 225 genotypes were selected by excluding lines with extreme phenological variation, considering an average for heading date of between 79 and 84 days.
Wheat genotypes were assessed in two sites located in the Mediterranean climate region of central Chile: Cauquenes (35 • 58 S, 72 • 17 W, 177 m.a.s.l.) and Santa Rosa (36 • 32 S, 71 • 55 W, 220 m.a.s.l). Cauquenes is a rainfed area and represented the water-limited conditions; annual precipitation was 580 and 600 mm in 2011 and 2012, respectively. Santa Rosa represented full irrigation (FI) and moderate water stress (MWS) conditions, which were achieved through support irrigation. Trials were assayed during two consecutive growing seasons, 2011 and 2012, except for the MWS trial, which was only set up during 2011. Cauquenes corresponds to the Mediterranean drought-prone area of Chile. According to del Pozo et al. [15], the average annual temperature is 14.7 • C, the minimum average is 4.7 • C (July), and the maximum is 27 • C (January). The average evapotranspiration is 1200 mm, and the annual precipitation was 410 and 600 mm in 2011 and 2012, respectively. No irrigation was applied in this site during 2011 and 2012. Santa Rosa corresponds to a high-yielding area. The average annual temperature in this area is 13.0 • C, the minimum average is 3.0 • C (July), and the maximum is 28.6 • C. The annual precipitation was 736 and 806 mm, in 2011 and 2012, respectively. Details about the trial establishment, experimental design, fertilization, control of weeds, irrigation, and soil water monitoring were previously detailed by del Pozo et al. [15]. Briefly, plants were grown in 2 × 1 m plots in rows separated by 0.2 m, with seeds sown at 20 g/m 2 . The experiment was established as an alpha-lattice design with 20 incomplete blocks, and 20 genotypes per block. A replicate comprised 400 plots and 384 genotypes. Two replicates were included in both trials, except for the MWS treatment in 2011, at both Cauquenes and Santa Rosa, where only single replicates were established.

Agronomic and Physiological Traits
Stem weights at anthesis (SWa) and maturity (SWm) were evaluated in five main stems per plot per genotype. The concentration of WSC in stems was determined at anthesis (WSCa) and maturity (WSCm), using the anthrone reactive method [61]. The WSC content (WSCC) was determined per unit of stem at anthesis (WSCCa = WSCa × SWa) and maturity (WSCCm = WSCm × SWm), and on a growth-area basis (m 2 ) for the same phenological stages (WSCCa/m 2 and WSCCm/m 2 ). The apparent WSC remobilization (ARWSC) was calculated as the difference in WSC content between anthesis and maturity. Grain yield (GY) was determined by harvesting the whole plot.

Statistical and Genetic Analysis
An approach based on a mixed linear model was applied to assess the genetic variation underlying the different traits. Initially the normal distribution was evaluated. A transformation utilizing natural logarithms was used to normalize WSC, WSCC, and WSCC/m 2 , for both the anthesis and maturity measurements. An ANOVA was carried out for each trait, using the following simplified model: where y ijkl is the response measured in the l-th plot of the i-th genotype within the j-th site and the k-th replicate, µ is the overall mean, G i is the random effect of genotype i, S j is the fixed effect of site j, R(S) k(j) is the fixed effect of replicate k within site j, GxS ij is the random interaction effect of genotype i with site k, and e ijkl is the random error in replicate k at site j within genotype i. ANOVA was executed by using the GLM procedure in SAS 9.2 [62]. Additionally, variance components were estimated by Restricted Maximum Likelihood, using the MIXED procedure. An estimation of the relative importance of genetic variation (σ 2 G ) on phenotypic variation (σ 2 p ) was estimated by the σ 2 G /σ 2 p ratio, where σ 2 p = σ 2 G + σ 2 GxS + σ 2 e , and σ 2 GxS and σ 2 e are the variance components estimated for the GxS interaction and residual, respectively. These analyses were done only with data from Cauquenes and Santa Rosa, for the 2012 season, considering that both datasets included two full replicates.

Genotyping
Genotyping by sequencing (GBS) was performed, using the Illumina HiSeq 2000 sequencing platform and the bioinformatic approach for SNP identification described previously by Lado et al. [63]. SNP-markers with a minor allele frequency < 5% and departure from Hardy-Weinberg equilibrium, as well as those with more than 5% missing data, were removed for the subsequent analysis stages, utilizing JMP Genomics 7.1 (SAS Institute, 2013a). The depurated database included 3243 SNP-markers selected from an initial dataset comprising 8746 SNPs.

Linkage Disequilibrium
Linkage disequilibrium (LD) was estimated among pairwise combinations of SNPs per chromosome. It was expressed in terms of the squared correlation of allele frequencies r 2 . The r 2 value between pairs of SNP markers, within each chromosome, was estimated, using TASSEL 5.2 [64], with the option of sliding window approach (120 SNPs per window). To assess the extent of LD, the decay of LD within physical distance (base pairs) between SNPs, within each chromosome, was evaluated by nonlinear regression analysis of r 2 values [65]. Analysis was performed by applying the NLIN procedure in SAS 9.4 [62].

Association Analysis
SNP-marker trait association tests were performed using a mixed linear model (MLM), adjusted for population structure and genetic relatedness, utilizing using TASSEL 5.2 [64]. The population structure matrix (Q) was obtained by principal component analysis (PCA), whereas the genetic relatedness (kinship) was incorporated as a scaled identity-by-descent matrix (K). Accession means were adjusted by Best Linear Unbiased Predictor, using Proc MIXED in SAS v9.2 [62]. The significance associated with the results of the SNP-marker trait association tests were subjected to multiple testing correction by false discovery rate (FDR) implemented through the program QVALUE in R [66,67].

Haplotype Analysis
Association tests between traits and selected genomic regions containing significantly associated single-SNPs were carried out, using the Haplo.stats package in R [68]. Haplotypes were estimated for each selected genomic region, and their frequencies were determined by using the modified expectation maximization method of haplotype inference [69]. Haplotypes and associations with traits were determined by utilizing genotyping, trait, and PCA matrices for the genotype population. Singleton alleles were discarded, as well as haplotypes with a frequency < 0.05. A global score statistic and haplotype-specific scores were derived from generalized linear models. A correction for multiple testing was also done by FDR, similar to single marker association tests. Significant differences among least-square means of haplotype classes were determined by ANOVA and the Tukey HSD test, using the GLM procedure in SAS 9.4 [62].

Gene Models and Annotations
Gene models and annotations for genes containing significant SNPs were obtained from the Triticum aestivum genome sequence available at the Ensembl platform [70]. Functional classes were defined using Quick GO [71].

Conclusions
Wheat is a key crop for world food sustainability, so genetic improvement in yield for sites with limited water availability must be a priority under the climate-change scenarios projected for this century. In that framework, identification and characterization of the genetic mechanisms underlying adaptive and productive traits are fundamental. In this study, we identified a series of SNPs controlling the variation in WSC traits in genotypes grown under different water regimes, in the Mediterranean-climate region of Central Chile. Most of these significant SNPs are part of gene-encoding proteins involved in the defensive response to drought and biotic stress. These SNPs represent a starting point to define candidate genes for functional characterization of their biological role under water stress. Our results also showed a complex and polygenic inheritance of the assessed traits. In that sense, the relative influence of genetic or environment (sites and irrigation conditions) factors was variable, depending on traits. Results from this study provided a set of SNPs that could be utilized to support marker-assisted selection suitable for future genetic improvement of spring wheat lines under water-limited areas.
Supplementary Materials: The following are available online, at https://www.mdpi.com/2223 -7747/10/3/539/s1. Table S1: General statistics for the assessed 225 wheat genotypes. Table S2: Number of significantly associated SNPs by chromosome and trait. Table S3: List of significantly associated SNP-markers. Table S4: Distribution of analyzed SNP-markers across wheat chromosomes. Figure S1: Comparative distribution of data for Cauquenes and Santa Rosa site, in growing season 2012. Figure S2: Venn diagrams for significantly associated SNPs according to sites and growing seasons. Figure S3: Decay of linkage disequilibrium (LD) by chromosome.
Author Contributions: A.d.P., I.M., and A.Y. planned and designed the research; F.P.G. and A.Y. analyzed data; F.P.G., A.d.P., and A.Y. interpreted data; A.Y., A.d.P., and I.M. conducted fieldwork and data and sample collection; A.Y. processed and analyzed samples; F.P.G., A.d.P., and A.Y. wrote the manuscript; A.d.P. and I.M. applied for funding. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: The data are available from the corresponding author upon reasonable request.

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