Genome‐wide association analysis of salinity responsive traits in Medicago truncatula

Abstract Salinity stress is an important cause of crop yield loss in many parts of the world. Here, we performed genome‐wide association studies of salinity‐stress responsive traits in 132 HapMap genotypes of the model legume Medicago truncatula. Plants grown in soil were subjected to a step‐wise increase in NaCl concentration, from 0 through 0.5% and 1.0% to 1.5%, and the following traits were measured: vigor, shoot biomass, shoot water content, leaf chlorophyll content, leaf size, and leaf and root concentrations of proline and major ions (Na+, Cl−, K+, Ca2+, etc.). Genome‐wide association studies were carried out using 2.5 million single nucleotide polymorphisms, and 12 genomic regions associated with at least four traits each were identified. Transcript‐level analysis of the top eight candidate genes in five extreme genotypes revealed association between salinity tolerance and transcript‐level changes for seven of the genes, encoding a vacuolar H+‐ATPase, two transcription factors, two proteins involved in vesicle trafficking, one peroxidase, and a protein of unknown function. Earlier functional studies on putative orthologues of two of the top eight genes (a vacuolar H+‐ATPase and a peroxidase) demonstrated their involvement in plant salinity tolerance.

intracellular compartments are high. Proline also acts as a reactive oxygen species (ROS) scavenger and molecular chaperone to stabilize proteins and bio-membranes under stress (Ashraf & Foolad, 2007;Matysik, Alia, & Mohanty, 2002). Proline biosynthesis serves as a redox buffer by consuming two NADPH per proline molecule, which helps to utilize excess electrons generated in the chloroplast under stress (Hare & Cress, 1997). Although generally regarded as a beneficial osmoticum, the link between proline and stress tolerance is somewhat unclear. In some studies, proline was found to accumulate more in tolerant plant genotypes than in sensitive genotypes under salt stress, consistent with an active role in stress tolerance (Jain, Nainawatee, Jain, & Chowdhury, 1991;Misra & Gupta, 2005;Ranganayakulu, Veeranagamallaiah, & Sudhakar, 2013). In other studies, however, proline levels are not positively correlated with salinity tolerance, but instead appear to increase as a consequence of cellular damage (Aziz, Martin-Tanguy, & Larher, 1998;Lacerda, Cambraia, Oliva, & Ruiz, 2003;Moftah & Michel, 1987). Nonetheless, at least one genome-wide association study (GWAS) has been performed to identify genes controlling proline level under dehydration stress in Arabidopsis (Verslues, Lasky, Juenger, Liu, & Kumar, 2014).
In leaves, high Na + concentrations cause stomatal closure due to the osmotic effect of the solute, which consequently causes a reduction in the rates of photosynthesis and growth (Brugnoli & Lauteri, 1991;Munns & Tester, 2008). Another consequence of decreased stomatal aperture is that intercellular CO 2 limitation causes NADP + pool depletion, photorespiration enhancement, and ROS accumulation (Hossain & Dietz, 2016).
To improve plant salinity tolerance, multiple approaches have been pursued. Buoyed by discovery of genes involved in plant salinity responses, genetic engineering has been used to generate a variety of transgenic plants in attempts to increase stress tolerance. Although positive outcomes have repeatedly been reported for transgenic plants under controlledgrowth conditions, less success has been found under field conditions (Fita, Rodríguez-Burruezo, Boscaiu, Prohens, & Vicente, 2015;Flowers, 2004;Hanin, Ebel, Ngom, Laplaze, & Masmoudi, 2016). On the other hand, conventional breeding for salinity tolerance, through selection and introgression, has met with limited success in rice and wheat, possibly because of the complex nature of this trait (Ashraf & Foolad, 2013;Munns, James, & Läuchli, 2006;Varshney, Bansal, Aggarwal, Datta, & Craufurd, 2011).
In recent years, with advances in next generation sequencing technology, great progress has been made in identifying quantitative trait loci (QTL), associated with salt tolerance (Ashraf & Foolad, 2013;Hamwieh et al., 2011;Thomson et al., 2010). For example, a Na + /H + antiporter (GmCHX1) was identified as a major salt-tolerance gene in the soybean wild accessions by combining whole-genome sequencing with QTL mapping (Qi et al., 2014). Later, GWAS indicated that this gene is a major contributor to the variance in salinity sensitivity among multiple soybean ecotypes; a conclusion supported by marker-assisted selection that resulted in salt-tolerant lines of soybean Patil et al., 2016;Zeng et al., 2017). Because cultivated soybean genotypes that contain nonfunctional GmCHX1 are generally very sensitive to salinity stress, it was postulated that loss of function of this gene in soybean may improve growth and seed production in nonsaline environments (Qi et al., 2014). Similar studies in Arabidopsis identified the sodium transporter, AtHTK1, as a major contributor to variance in shoot sodium accumulation under salinity stress in wild populations (Baxter et al., 2010;Busoms et al., 2015;Rus et al., 2006). Another GWAS study in Arabidopsis (Julkowska et al., 2016) indicated that a kinase, LRR-KISS (Leucine-Rich Repeat Kinase family protein Induced by Salt Stress), underlies variance in plant growth under salinity stress.
Apart from the progress noted above for soybean and Arabidopsis, little is known about the genes shaping natural variation in salinity tolerance in other species. A few traits affected by salinity stress have been investigated via GWAS with high density markers, including traits related to seed germination, growth rate, transpiration rate, and tissue Na + /K + contents in rice (Al-Tamimi et al., 2016;Patishtan, Hartley, Fonseca de Carvalho, & Maathuis, 2017;Shi et al., 2017), and root growth in Arabidopsis (Kobayashi et al., 2016). However, no overlapping genes or common mechanisms have been identified among these studies, and none of the identified genes have been validated with respect to salinity tolerance. Apart from soybean, no GWAS analysis of salinity traits with high-density markers has been performed in other legumes.
Medicago truncatula is a model legume species for genetic and genomic research (Barker et al., 1990;Kang, Li, Sinharoy, & Verdier, 2016;Young & Udvardi, 2009 (Young et al., 2011). In the current study, we characterized a spectrum of salinity stress-related traits in a diverse subset of the M. truncatula HapMap panel consisting of 132 M. truncatula genotypes.
Plants were grown in soil and were subjected to a step-wise increase in NaCl concentration, from 0 through 0.5% and 1.0% to 1.5%. Traits that were measured in salt-stressed (and control) plants included a qualitative, visual salt-tolerance score, relative shoot biomass reduction (%), leaf chlorophyll content reduction (%), leaf size reduction (%), shoot water content (%), and leaf and root concentration of proline and major ions (Na + , Cl − , K + , Ca 2+ , Mg 2+ , NH 4 + , PO 4 3+ , SO 4 2− , NO 3 − , and malate). GWAS was carried out using 2.5 million high-quality SNP markers. By performing GWAS on complex, emergent traits like growth/biomass as well as on potentiallyunderlying, fundamental traits or "phenes" such as ion and metabolite concentrations, which reflect cellular ion homeostasis and metabolism, we hoped to break-down salinity tolerance into its component parts while at the same time identifying genetic loci for tolerance in a more refined and robust manner.

| Plant material
Seeds of all the lines used in this study were obtained from the M. truncatula HapMap germplasm resource center (http://www. medicagohapmap.org/hapmap/germplasm). One hundred thirty-two lines (Table S1) were selected based on population structure for maximum variance from the 220 lines that were used previously for the GWAS of drought-related traits (Kang et al., 2015).
M. truncatula seeds were scarified on sand paper (p800) and germinated on wet filter paper for 48 hr at 4°C by overnight storage in the dark at room temperature, followed by planting. Plants were grown individually in 2″ × 7″ plastic cones (Stuewe & Sons Inc.) containing a 4:1 (v/v) mixture of Metromix 350 and sand, in a growth chamber with 16 hr day/8 hr night, 22°C, and 40% relative humidity.
For each experiment, 12 seedlings were planted for each line, six for control conditions and six for salinity-stress treatment. The plants were placed in the growth chamber in a randomized complete block design. Light density (a combination of fluorescent and incandescent) at plant level was approximately 200 μmol m −2 s −1 . Plants were watered with one fourth B&D nutrient solution (Broughton & Dilworth, 1971) containing 8 mM of N (2 mM KNO 3 , 3 mM NH 4 NO 3 ).
Three replicates were performed for the entire experiment.

| Salinity stress treatment
Ten days after seedling transfer to pots (day 10), plants were watered with 0.5% NaCl (86 mM) in one half B&D nutrient solution, which was applied to both the tray holding the cones and the top of the soil with a squeeze bottle to avoid leaf damage. Excess solution in the trays was removed 2 hr after the salinity solution treatment. Sodium chloride concentration in the nutrient solution was increased to 1.0% (172 mM) at day 15 then 1.5% (257 mM) at day 20, in the same way. No extra watering between treatments was necessary because of reduced transpiration. Roots and shoots were harvested 5 days following the 1.5% NaCl treatment.

| In-vivo leaf chlorophyll content measurement
Four days following the 1.5% NaCl treatment, in-vivo leaf chlorophyll content was measured on salt-stressed and control plants, using a Chlorophyll Meter SPAD-502plus (Spectrum Technologies Inc., http://www.specmeters.com/). The terminal leaflet of the uppermost fully-expanded leaf was used for the measurement. Each leaflet was measured twice at different positions, avoiding the midrib and edges, and the average of the two readings was recorded.

| Leaf size measurement and plant harvest
Five days following the 1.5% NaCl treatment, the two youngest fullyexpanded leaves, from both control and NaCl-treated plants, were harvested and immediately used for individual leaf size measurement with the Li-3000A (LI-COR) portable area meter. These two leaves were then flash-frozen in liquid nitrogen, and stored at −80°C until being dried at −40°C in a lyophilizer. After the two youngest fullyexpanded leaves were harvested, the rest of the shoot was harvested separately. Shoots were dried in a 55°C oven and weighed. Total shoot dry weight was the sum of the two youngest fully-expanded leaves and the rest of the shoot. Roots were harvested at the same time as the shoots, rinsed well, flash-frozen in liquid nitrogen, and then stored at −80°C until being dried at −40°C in a lyophilizer. Dried tissues were ground to powder in a bead beater (BioSpec, https:// www.biospec.com/).

| Tissue ion (including malate) measurement
Approximately 5 mg of dry ground leaf or root tissue (six plants pooled together) was dissolved in 5-14 ml of MQ water, vortexed, and then shaken at 200 rpm for 1 hr. The solution was then filtered through a 0.2 μm nylon filter (F13-2020, VWR) and the filtrate was used for measurement of total ions with ion chromatography. Chromatographic separation was achieved on a Thermo Scientific ICS-5000 IC system (Thermo Fisher Scientific, USA) using a Dionex CS12A Ion Pac analytical column with a AG12A guard column for cations, or a Dionex AS11HC analytical column with a AG11HC guard column for anions.

| Tissue proline quantification
Tissue proline levels were analyzed with a biochemical assay modified from Bates et al (Bates, Waldren, & Teare, 1973) and Hamid et al (Hamid et al., 2003). Proline concentration was determined using a standard curve using L-proline.

| Filtering of SNPs
SNPs were filtered with TASSEL 5.2.7 (Bradbury et al., 2007), with a minimum allele frequency of 5% and minimum counts of 112 (85% of 132). Missing SNPs were not imputed because of the high density of the existing SNPs. After filtering, a total of 2,528,531 SNPs remained and were used in the association analysis.

| Cladogram
The cladogram tree was generated in TASSEL 5.2.7 (Bradbury et al., 2007) with the neighbor-joining method using 40,000 randomly selected SNPs (5,000 SNPs each chromosome).
2.9 | Q matrix deduction Q matrix was deduced with the STRUCTURE program (Pritchard, Stephens, & Donnelly, 2000) using 40,000 SNPs (5,000 random SNPs from each chromosome). We evaluated K = 1-9 to infer the optimal value of K (i.e., the number of clusters) from the simulation summary using the methods of Pritchard et al. (2000) and Evanno, Regnaut, and Goudet (2005).
2.10 | Genome-wide genotype-phenotype association analysis The least square means of the phenotypic data collected in three replicates were calculated in R with library "lsmeans" and used in GWAS.
The proline and ion measurements were carried out only on the last replicate, and the standard-calibrated values obtained were used directly in GWAS. The mixed linear model (MLM) and general linear model (GLM) embedded in TASSEL (Bradbury et al., 2007) were used to test for association between SNPs and phenotypic traits. For the MLM analyses, a Kinship matrix was calculated with centered-IBS in TASSEL (Bradbury et al., 2007), and no compression was applied.
Quantile-quantile (QQ) plots were generated in R (for GLM and MLM-Q) or using 10% of all SNPs randomly-selected from the eight chromosomes in TASSEL (for MLM + Q). Linkage disequilibrium among SNPs was also performed in TASSEL 5.2.7 (Bradbury et al., 2007).

| Total RNA extraction and qRT-PCR
Total RNA from shoot and root was extracted using RNeasy Plant Mini Kit (Qiagen) followed by genomic DNA removal and column purification. Reverse transcription was performed with oligo dT 20 and Super Script III Reverse Transcriptase as described previously (Kakar et al., 2008). Amplification of templates followed standard PCR protocols with SYBR Green-based detection system. Transcript levels were normalized using the geometric average of three housekeeping genes, Ubiquitin-Conjugating enzyme E2 (TC106312), Polypyrimidine Tract-Binding protein (TC111751), and Ubiquitin (TC102473). The primers used for all the target and housekeeping genes are listed in Table S6.

| Statistical analysis
Principal component analysis (PCA) was performed in R using the "PCA" function in the package "FactoMineR" (Lê, Josse, & Husson, 2008). Association analysis among traits was performed in SAS 9.4 (SAS Institute, Cary, NC) using the CORR procedure. False discovery rate (q value) was calculated using the R package "qvalue" (Storey & Tibshirani, 2003). Significant tests were performed in Excel 2013 using student's t-test, two-tails, assuming equal variance.

| Phenotypic data collection and correlations among traits
Preliminary experiments with 30 M. truncatula HapMap lines were carried out to establish a salinity stress regime that affected plant growth without overwhelming plants completely. As a result, a protocol with sequential increases in NaCl concentration in the nutrient solution was chosen: 0.5% NaCl (86 mM) for 5 days, 1% (171 mM) for 5 days, and 1.5% (257 mM) for 5 days. These treatment conditions resulted in significant leaf chlorosis, leaf size reduction, and growth retardation compared with non-treated control plants ( Figure 1). On the basis primarily of plant vigor and the degree of leaf chlorophyll loss, we scored stressed plants from 1 to 5, with 5 being the most tolerant ( Figure 1a).
In addition, shoot biomass, leaf chlorophyll content, and leaf size were (c) Leaves. In b and c, the plant/leaf on the left was non-stressed. All photos were taken 4 days after the 1.5% NaCl treatment, 1 day before plant harvest distribution was observed; leaf K + /Na + ratio was an exception ( Figure 2). Notably, concentrations of sodium, chloride, calcium, and proline were much higher in leaves than roots of NaCl-treated plants, whereas potassium concentration exhibited the opposite response to salinity ( Figure 2g-i,m-o; Figure S1). As a consequence, the K + /Na + ratio was approximately 2-fold higher in the roots than in the leaves ( Figure 2j,p).
To reveal relationships between traits, correlation analyses were performed. In general, shoot and leaf traits were positively correlated but, with the exception of proline concentration, shoot and root traits were less correlated (Table 1). Relatively strong and significant correlations were found between all of the following traits: shoot tolerance scores; reduction in leaf size, shoot dry weight, and leaf chlorophyll; leaf concentrations of sodium, chloride, and proline; and water content of the shoot (Table 1). Compared with chloride, sodium, calcium, and the K/Na ratio, potassium concentration in leaves were less tightly correlated to tolerance scores (r = −0.35). Sodium and chloride concentrations were extremely highly correlated to each other, both in leaves and roots, with correlation coefficients of 0.99 and 0.89, respectively. Interestingly, proline levels in leaves were negatively correlated to tolerance scores and shoot water contents, but were positively correlated to leaf sodium concentrations, whereas root proline levels were positively correlated with salinity tolerance ( Figure 3).
Thus, salinity-tolerant M. truncatula genotypes tended to accumulate less proline in the leaf but more proline in the root under salinity stress, compared with the sensitive genotypes.
Several other ions were analyzed, and leaf magnesium, leaf sulfate, leaf ammonium, leaf nitrate, root malate, and root phosphate were found to be significantly correlated with salinity tolerance, either positively or negatively (Table S2).
3.2 | Genome-wide association analysis and GO enrichment of "top-suspect" genes Genome-wide association analysis using mixed linear model or general linear model was performed in TASSEL (Bradbury et al., 2007;Zhang et al., 2010). In the cladogram tree, the 132 lines used in the current study were clearly separated into two major groups ( Figure   S2). This is consistent with the cluster number determined by the program structure (Pritchard et al., 2000), which was also two ( Figure   S3). Because of the existence of population structure in the GWAS population, QQ plots generated with GLM or MLM without Q-value correction were inflated for majority of the traits ( Figure S4). To control for false positives, the population structure (Q value) was included in the MLM; GWAS Manhattan plots and QQ plots were generated and the lowest P values ranged from 10 −6 to 10 −13 (Figures 4, Figure S4). From the Manhattan plots, chromosome 2 appeared to be a "hot spot" for SNPs with lowest P values in multiple traits. For each trait, we collected the 200 SNPs with the lowest P values, the genes in each SNP's vicinity, as well as the closest Arabidopsis orthologues and M. truncatula microarray probe-sets (Table S3).
To gain an overview of genes that may contribute to salinity tolerance, we performed gene ontology (GO) enrichment analysis of all genes in the vicinity of the 200 SNPs with the lowest P values for 14 traits (Table S3). We used the GO of Arabidopsis orthologues for this analysis because they are better curated than for Medicago.

| Identification of potential causative SNPs and genes
In analyzing the top 100 SNPs associated with each of the 14 traits, we found a large portion of them to be linked to multiple traits (Table   S3), forming distinct "hotspots" on chromosomes especially on chromosome 2 ( Figure 5). Considering the tight association among different traits (Table 1), we selected the top SNPs that were tightly associated with multiple traits, as well as the genomic regions that contained these SNPs. In doing so, we identified 12 genomic regions (QTLs) harboring SNPs that are in tight association with at least four traits (rank 200 or higher; Figure 6). In parallel, we performed PCA analysis of the first 12 traits (ToleranceScore, ST_DWred, LF_ChlRed, LF_SizeRed, ST_Water%, LF_Proline, LF_Sodium, LF_Chloride, LF_Potassium, LF_KNaRatio, LF_Calcium, and RT_Proline) that are tightly correlated in Table 1 and performed GWAS analysis using PC1 (that explained 55% of the variance) as a new trait, which identi-  Table S3). The central SNP is 2:14451314, which has low association P-values in seven traits. A total of six genes are located in this region ( Figure 6b; Table S3). In addition to this region on chromosome 2, three other regions were identified (Figure 6a,c,d).  Table S3). Other less significant regions are highlighted in Table S3.

| Evaluation of potential salinity-tolerance genes
We investigated further potential salinity-tolerance conferred by "sus- In an effort to further narrow down our list of suspect genes, we analyzed how these genes are regulated under salinity stress in M. truncatula, using published data from a study that employed hydroponics and relatively short-term (up to 48 hr) salinity stress (Li, Su, data, and all but three of these genes were regulated under one or more of these three stress conditions (Table 2).
In addition to examining published gene expression data, we chose eight of the top suspect genes ( Table 2 in bold) and performed In analyzing SNP variance in these eight genes, we extracted all the SNPs that reside within the genes and 1000 bp upstream and downstream of the start and stop codons, and selected those that varied consistently between tolerant and sensitive lines (Table S4). A total of 46 SNPs matching these criteria were identified, with Medtr7g081010 (H + -ATPase) containing the most, whereas there were none in Medtr2g436940 (unknown protein). One missense mutation was found in Medtr4g046713 (peroxidase family protein), which caused an amino acid switch from Tyrosine to Serine. In addition, 10 SNPs occurred in introns and two in UTRs.
Finally, it is interesting to note that the profile of these 46 SNPs in the M truncatula A17 reference genome is highly similar to the two sensitive lines (HM081 and HM152), with only five being different (highlighted in Table S4). A17 was not used in the GWAS analyses but was characterized for salinity sensitivity as a check line; it has an average salinity sensitivity score of 2.1 and was, therefore, relatively sensitive to salinity stress ( Figure S7).

| Correlation among salinity responsive traits
Salinity is a complex abiotic stress, with ionic and osmotic stress components that trigger a variety of responses in plants. In the current study, we applied long-term, gradual salinity stress to soil grown M. truncatula plants and characterized a wide-range of physiological and biochemical traits in conjunction with GWAS. In the M. truncatula HapMap panel, we observed substantial variation in responses to salinity stress (Figure 1), and almost all of the traits characterized demonstrated a normal distribution ( Figure 2). With multiple traits characterized, we were able to analyze the relationships among different traits as rarely done before. Sodium and chloride levels in leaves (r = −0.55, P < 0.001 for both) but not in roots (P > 0.2) were highly (negatively) correlated with salinity tolerance scores, with tolerant genotypes containing lower concentrations of sodium and chloride (Table 1). Therefore, it appears that maintaining low sodium and chloride concentrations in the shoot is an important strategy used by M. truncatula plants to survive and grow under salinity stress. This is consistent with earlier reports that photosynthetic organs are more sensitive than roots to high salt (Munns et al., 2006).  reporter when the opposite trend was found (Aziz et al., 1998;Heuer, 2010;Lacerda et al., 2003;Moftah & Michel, 1987). However, all previous studies used a small number of genotypes and only focused on the shoot. Here, we demonstrated that leaf proline levels are negatively correlated with tolerance scores (r = −0.59, P < 0.0001) and shoot water content (r = −0.65, P < 0.0001), whereas root proline levels are positively correlated with tolerance (r = 0.33, P = 0.0002; Figure 3; Table S2). In other words, salinity-tolerant M. truncatula genotypes tended to accumulate more proline in the roots but less proline in the leaves than sensitive genotypes. This phenomenon has not been reported before and we postulate that it may reflect different roles of proline in roots and leaves. When plants are under salinity stress, the primary stress encountered by photosynthetically-active leaves is redox stress (ROS) associated with insufficient ability to channel high-energy electrons into anabolic pathways due to stomatal closure and reduced carbon fixation (Chaves, Flexas, & Pinheiro, 2009;Hossain & Dietz, 2016). Proline biosynthesis consumes electrons (NADPH), helping plants to cope with increased ROS production when stomata close (Hare & Cress, 1997). The situation in roots appears to be quite different, where proline levels correlate positively with salinity tolerance. Rather than being subject to increased ROS associated with light-energy harvesting without carbon fixation in leaves when stomata close, roots are prone to dehydration stress in saline soils. Proline biosynthesis in roots presumably contributes to overall osmolyte production in saline soils, which helps root cells take up the water necessary for transpiration, stomatal opening, photosynthesis, and so forth. By producing more proline in roots, salinity-tolerant plants may be better able than sensitive plants to support shoot functions and carbon supply back to the root for further proline/osmolyte synthesis. In this sense, proline biosynthesis in the root may be a primary salinity tolerance mechanism in Medicago. Further studies are required to test this hypothesis, for instance by down-regulating proline biosynthesis via RNAinterference in transgenic "hairy roots" but not in shoots of Medicago.
It would also be interesting to determine whether proline accumulates differentially in roots and/or shoots of plants with differential tolerance to other kinds of stress.

| Identification and validation of suspect drought tolerance genes
To reduce false positives in our genome-wide association studies, we used MLM and included population structure as co-variants to counter bias due to any population structure. With this stringent MLM + Q + K analysis, we identified SNPs with P values as low as 10 −13 (leaf K/Na ratio), which is far below the P-value threshold after stringent Bonferroni correction (1.98 × 10 −8 with 95% confidence). Taking advantage of the many salinity-related traits that were measured, we focused on SNPs that were repeatedly linked to multiple traits, rather than simply selecting SNPs based on applying thresholds to individual traits. With this approach, 12 genomic regions with clear borders were identified ( Figure 6, Table S3). A total of 74 genes and 214 SNPs were present in these 12 regions, with 94 SNPs within genes. Of these genes, 39 had corresponding transcript information in microarray data sets, and the majority were found to be either regulated by short-term salinity stress or drought/desiccation stress in earlier studies (Li et al., 2009;Verdier et al., 2012;Zhang et al., 2014; Table 2). The tight LD among the top overlapping 214 SNPs ( Figure 8) across multiple chromosomes indicates possible co-evolution of the salinity tolerance genes. This phenomenon has been previously reported in human and mouse (Kulminski, 2011;Petkov et al., 2005). Note, however, that although multiple traits may map to the same genomic region, they are not necessarily associated with the same SNPs ( Figure 6, Table   S3). From GO enrichment analysis, we found that genes in the categories "response to stress" and "response to stimulus" are highly enriched with FDR < 0.001. Taken together, it is evident that the current GWAS studies identified a set of genes that are enriched in salinity/dehydration stress responses.
In examining the expression patterns of the eight top suspect genes by qRT-PCR analysis of the five extremely tolerant and sensitive lines, we found significant regulation of these genes in response to gradual salinity stress ( Figure 9). Association of salinity-stress sensitivity with gene transcript level changes was evident in seven of these eight genes including one H + -ATPase (Medtr7g081010), two transcription factors (Medtr2g436900, Medtr2g437020), two genes likely to be involved in vesicle trafficking (Medtr2g028750, Medtr2g028790), one gene involved in cellular redox homeostasis (Medtr4g046713), and one with unknown function (Medtr2g436940).
The majority of these genes tended to have stable or increased expression under salinity stress in the tolerant lines but sharply decreased expression in the sensitive lines, especially in the root ( Figure 9). The overall gene regulation patterns in the root bear similarity to those observed previously under short-term salinity stress (Li et al., 2009) (https://mtgea.noble.org/v3/index.php). SNP variances correlating with salinity stress sensitivity were identified in these genes, including one that causes a missense mutation, and 12 causing modifications in the introns and the UTR regions (Table S4).
Vesicular trafficking has been demonstrated to play important roles in plant adaptation to salinity stress (Baral et al., 2015;Garcia de la Garma et al., 2015;Hamaji et al., 2009). Vesicle trafficking is likely to be involved in deployment of specific ion, water, or metabolite transporters to the cell or organellar membranes for ion or water homeostasis (Baral et al., 2015). Vesicle trafficking may also be  (Li et al., 2009); b Zhang et al. (Zhang et al., 2014); c Verdier et al. (Verdier et al., 2012). *, **, and ***indicate significant differences at p < 0.05, 0.01 and 0.001, respectively.

FIGURE 9
Transcript levels of potential causative genes in tolerant and sensitive genotypes. Relative transcript levels were determined by qRT-PCR. The tolerance scores for HM091, HM010, HM198, HM081, and HM152 were 4.2, 3.8, 3.5, 1.7, and 1.5, respectively. Transcript levels are expressed relative to the mean of three housekeeping genes (Ubiquitin-Conjugating enzyme E2, Polypyrimidine Tract-Binding protein and Ubiquitin). n = 3. Error bars represent standard errors. Significance tests were between control and salinity stress treatments only in each genotype. Significat differences at *P < 0.05, **P < 0.01, and ***P < 0.001. Ctl: control. Salt: salinity-treated involved in ROS signaling and facilitating plasma membrane area reduction during plasmolysis caused by osmotic stress (Baral et al., 2015;Leshem, Seri, & Levine, 2007). Here, we identified two genes that are involved in vesicle trafficking: one synaptobrevin-like protein (Medtr2g028790), which is an R-SNARE (Soluble N-ethylmaleimidesensitive factor Attachment Protein) that mediates vesicle fusion, and one nonclathrin coat protein zeta1-COP (Medtr2g028750), which is a component of the COPI coatomer that coats vesicles transporting proteins between the Golgi complex and the endoplasmic reticulum.
Both genes had clear transcription reduction in the root under salinity stress in the sensitive genotypes but not in the tolerant genotypes in M. truncatula (Figure 9). Earlier studies demonstrated that modification of various SNAREs could either enhance or reduce salinity tolerance in Arabidopsis depending on the experimental set-up (Hamaji et al., 2009;Leshem et al., 2007;Tarte et al., 2015). However, the involvement of non-clathrin coat proteins in salinity stress responses has never been studied. Our GWAS results suggest a potential role of non-clathrin coat protein zeta1-COP in this process, which, thus, deserves further investigation.
Vacuolar H + -ATPase (V-ATPase) generates and maintains a proton gradient across the tonoplast that provides the driving force for secondary transport processes, including storage of excess ions in the vacuole under salinity stress (Silva & Gerós, 2009). Yeast two-hybrid experiments showed that SOS2 (salt overly sensitive 2) interacts directly with V-ATPase subunit B1 and B2 (AT4G38510, putative orthologue of Medtr7g081010), and mutant plants with reduced V-ATPase activity were extremely salt sensitive (Batelli et al., 2007). In alfalfa, over-expression of V-ATPase subunit B from the halophyte, Suaeda corniculata (a putative orthologue of Medtr7g081010) resulted in plants more tolerant to salt and saline-alkali stresses (Wang et al., 2016). Similarly, overexpression of the putative wheat orthologue (GenBank accession number: EF105343) of the M. truncatula V-ATPase subunit B (Medtr7g081010) significantly improves germination rate under salinity and overall salt tolerance in Arabidopsis (Wang, He, Zhao, Shen, & Huang, 2011). In the current study, we found that the expression of V-ATPase subunit B2 (Medtr7g081010) remained stable in roots of tolerant lines, but decreased significantly in sensitive lines. In addition, we identified multiple SNPs in the UTR region, the intron, and 1 kb vicinity of this gene that differentiate the tolerant from the sensitive genotypes. Work is ongoing to locate the causative SNP(s) underlying the differential regulation of this gene under salinity stress.
The importance of ROS scavenging and the two transcription factors (mTERF protein and C2 domain protein) in plant salinity stress responses have been reported in several previous studies (Abogadallah, 2010;Quesada, 2016;Robles, Micol, & Quesada, 2015;Xu, Leister, & Kleine, 2017;Yokotani et al., 2009). In the current study, we identified a peroxidase (Medtr4g046713) as a gene associated with salinity tolerance as it is the only predicted gene in a genomic region that harbors 26 SNPs linked to six distinct traits (Figure 6i).
The Arabidopsis orthologue, rare cold inducible gene 3 (RCI3, AT1G05260), has been shown to increase salinity tolerance when overexpressed and to decrease salinity tolerance when mutated Apart from the genes discussed above that are linked to multiple salinity-related traits, we are interested in genes associated with single traits that have homologs in other species that have been reported to play important roles in salinity tolerance (Table S3), including potassium channels/transporters, AAA family ATPases, cation/H + exchanger, cation-chloride cotransporter, LEA proteins, and a delta-1-pyrroline-5carboxylate synthetase (P5CS), that is rate-limiting for proline biosynthesis. As one example, the critical role of maintaining potassium homeostasis in plants during salinity stress has been intensively studied (Shabala & Pottosin, 2010). It has been shown that overexpression of a K + /H + antiporter in tomato enhanced its salinity tolerance (Huertas et al., 2013).
Recently, GWAS of salinity stress-related traits with relatively high precision have been performed in Arabidopsis, rice, and rapeseed (Al-Tamimi et al., 2016;Kobayashi et al., 2016;Shi et al., 2017;Wan et al., 2017). We compared our GWAS results with these studies to identify potentially-conserved salinity tolerance genes (orthologues).
Among the top genes identified in various GWAS of salinity traits, only one, AT2G44480 (Table S5 in Kobayashi et al., 2016), which encodes a beta glucosidase, is homologous to a top gene identified in the current GWAS. Intriguingly, M. truncatula has three tandem repeats of this gene (Medtr4g066280, Medtr4g066330, and Medtr4g066340, Figure 6J), suggesting that gene expression levels may impact salt tolerance.
In summary, genome wide association studies (GWAS) of multiple physiological and biochemical traits related to salinity-stress adaptation in M. truncatula were carried out using 2.5 million SNPs. Correlation analysis among the traits showed that salinity tolerance in the population is negatively correlated to sodium, chloride, and proline concentration in the leaf, while root proline levels were positively associated with salinity-stress tolerance. GWAS analyses were performed with a mixed linear model, and 12 genomic regions that are associated with multiple traits each were identified; these regions harbor a total of 214 SNPs and 74 genes. These genomic regions were confirmed by a GWAS analysis using PC1 as a virtual trait, generated from principle component analysis of the 12 traits that are tightly correlated. Significant linkage disequilibrium was observed among these SNPs, possibly indicating co-evolution of salinity tolerance alleles.
Examination of gene expression of eight top suspect genes revealed associations between salinity tolerance and transcript level changes under salinity stress. Earlier functional studies on orthologues of two of the top eight genes (a vacuolar H + −ATPase and a peroxidase) confirmed their involvement in plant salinity tolerance in other species.
Functional annotation of potential salinity tolerance genes points to the importance of transcriptional regulation, vesicle trafficking, and ROS scavenging under saline conditions. medicagohapmap.org/hapmap/germplasm). The authors would like to thank Shulan Zhang, Elis Acosta-Torres and Alicia Smith for assistance during data collection, and Wenchao Zhang for assistance in data analysis. This work was supported by The Samuel Roberts Noble Foundation.