A panel of germplasm composed of 200 alfalfa accessions selected from the USDA-ARS National Plant Germplasm System (NPGS) alfalfa collection was used for evaluation and mapping loci associated with drought tolerance. A majority of germplasm were collected from Northwest regions of US and Canada including British Columbia, Saskatchewan, Manitoba, Idaho, Montana, Nebraska, Washington, and North and South Dakota. The remaining accessions were collected from different countries including Afghanistan, Algeria, Bulgaria, China, Germany, India, Lebanon, Oman, Russia, Spain, Turkey, and Yemen (Table S1).

Alfalfa accessions were planted on the Roza farm at the Irrigated Agriculture Research and Extension Center of Washington State University, Prosser, WA, in 2015 and 2016. Field experiments used a randomized complete block design with three replications. Each block contains nine plants with 15 cm between rows and 5 cm between plants. They were watered regularly until the plant size was uniform (day 35–40). Subsequently, the control plants were watered regularly while drought treatment was applied to the stressed plants by withholding water during the dry season (June–September). Periodic irrigation was applied as necessary. The exact watering interval was determined during the trial through visual evaluation of wilting of the plants and measurement of soil moisture. The soil moisture was measured using ProCheck with the GS3 sensor (Decagon Devices Inc., Pullman, WA) to ensure that water deficit conditions (i.e., that the average soil water potential was below −0.03 Mpa) were maintained until harvesting.

Early in the flowering stage, aboveground biomass was harvested, with three cuttings per year, and biomass yield was measured. A statistical analysis was conducted for each harvest. The equality of variance and means were analyzed by ANOVA. Meanwhile, a normality of distribution was tested for each harvest. Pearson correlation was calculated between harvests using trait means. All analyses were conducted using JMP13 (https://www.jmp.com/en_us/home.html).

High molecular weight DNA was extracted from fresh young leaves of alfalfa plants using the Qiagen DNeasy 96 Plant kit (Qiagen, CA), following the manufacturer's protocol. DNA was quantified using a Nanodrop 1000 spectrophotometer (Thermo Scientific, http://www.thermoscientific.com) using the absorbance at 260 nm and quality checked using the ratio of 260/280 nm. High molecular weight DNA was subsequently used for library preparation. Two GBS libraries were prepared by digestion with the EcoT22I restriction enzyme and the digested DNA was ligated to unique nucleotide adapters (barcodes) followed by PCR amplification as described by Elshire et al. (2011). Libraries were placed in two lanes of an Illumina Hi-Seq2000 instrument at Cornell University Sequencing facility (Ithaca, NY) to be sequenced using 100-base single-end sequencing.

Sequence reads with barcodes were collapsed into a set of unique sequence tags using FASTQ. The master tags were then aligned to the reference genome of M. truncatula (Mt4.0 v1) according to Glaubitz et al. (2014). A tags-on-physical-map (TOPM) file containing the genomic position of each tag with a unique alignment was generated based on the physical map of the reference genome. The number of times each tag appeared in the sequencing data was scored using the barcode information. The tags by taxa (TBT) file were recorded according to the master tags listed for each sample. The information recorded in the TOPM and TBT files was then used to discover SNPs at each tag locus. The proportion of taxa covered by the tag locus, minor allele frequency, and inbreeding coefficient (FIT) were used for filtering the SNPs. After filtering, a total of 10,327 SNPs, with a mean individual depth of 27 X were obtained. The Row data of GBS were submitted to the NCBI Sequence Read Archive with bioproject ID: PRJNA287263 and biosample accession numbers: AMN03779142-SAMN03779330.

GBS markers were further filtered with a cutoff value of 0.05 for minor allele frequency (MAF) and a missing value of 0.25. The remaining high quality SNPs were used for marker-trait association analysis with TASSEL (Bradbury et al., 2007). Marker-trait association was analyzed by a mixed linear model (MLM). A kinship (K) matrix was used for controlling population structure during association mapping. A false discovery rate (FDR) of 0.05 was used as a multiple-test correction for significantly-associated markers (Benjamini and Hochberg, 1995). To assign the associated markers to known genes, a BLAST search was carried out in the DNA database of Phytozome against the Medicago truncatula genome sequence, Mt4.0 v1 (http://phytozome.jgi.doe.gov/jbrowse/index.html?data=genomes%2FMtruncatula&loc) and the National Centre for Biotechnology Information (NCBI, http://www.ncbi.nlm.nih.gov/), using the flanking sequence of each significant marker as a query. Known genes linked to the significant loci were then assigned as putative candidates based on the annotation of gene functions.

Phenotypic data were analyzed using ANOVA for each cutting. The result is presented in Table 1. Biomass yields for all three cuttings from water deficit and control treatments were significantly different (P < 0.0001) among accessions. General speaking, biomass yield was significant reduced in water deficit than those of control. The parameter estimates showed variations between cuttings with highest value of 3.5 lbs in the 1st cutting of control plants and lowest of 0.48 lbs in the 3rd cutting of those of water deficit. The normality of distributions of three cuttings from control and water deficit treatments was analyzed and the result is presented in Figure 1. The mean values were higher in the 1st cuttings, with a lower mean value for the 3rd cuttings of both treatments. Correlations between cuttings were also analyzed and the result is presented in Table 2. Higher correlations (r = 0.47–0.82) was observed between cuttings of control plants, while lower correlations were found between treatments. However, higher correlation (r = 0.59) was also found between 1st cuttings of control and stressed plants (Table 2). This probably due to a small treatment effect on biomass was obtained during the 1st cuttings of the water stress treatment.

Biomass yield data of three harvests were used for marker-trait association with the same genotyping data set. To compare the genotypic variation of biomass under drought and control conditions in different seasons, each harvest was analyzed independently. Under the control condition, only four markers (S1_6017132, S1_276968267, S1_21394479, and S1_144642771) were significantly associated with biomass yield (Figure 2, Table 3). Whereas, a total of 28 markers were found to be significantly associated with biomass yield under water stress (Figure 3, Table 3). Among them, 18 markers were identified in all three harvests (Table 3, highlighted with yellow) and six were identified in two harvests (Table 3, highlighted with blue). Only four markers were identified in one harvest (Table 2, unhighlighted). Since a high level of synteny was observed between the alfalfa linkage maps and the M. truncatula physical map (Li et al., 2014), we aligned GBS sequence tags to the reference genome of M. truncatula (Mt 4.0 v1). Based on the M. truncatula physical map, the associated markers were located at 22 loci throughout all chromosomes (Table 3, Figure 3). Since several makers were located at the same region within 60 bases of one another, we treated them as single locus, although they are listed as different markers as they have different SNP sites (Table 3). Nine markers were not mapped on the M. truncatula physical map. We performed BLAST searches against the DNA databases of NCBI and Phytozome using their flanking sequences as queries. Seven markers were reassigned to specific chromosomes based on the BLAST searches (Table 3, “^*”). The remaining two markers were still unable to be aligned to the reference genome (Table 3, “U”).

The markers' p-values ranged from E-4 to E-12 with marker S1_92116984 being the lowest (2.33E-12) and marker S1_41503160 being the highest (5.02E-4). Most significant markers were located on chromosomes 1, 3, and 4. Similar trend was found for the markers' r² values which ranged from 0.10 to 0.32 (Table 2).

To identify possible candidate genes underlying the resistance loci, The BLASTN algorithm was performed using the flanking sequences of the significant markers as queries against the Phytozome (http://phytozome.jgi.doe.gov/pz/portal.html) M. truncatula genome sequence (Mt4.0 v1) and the NCBI (www.ncbi.nlm.nih.gov/) nucleotide databases. Four markers (S1_20593175 to 182) at the same locus on chromosome 1 were linked to a leucine-rich repeat receptor-like kinase (LRR; Table 2, Figure 3). Another marker S1_368369179 on the same chromosome was linked to exportin 1-like protein. On chromosome 2, markers S1_50692984 and S1_324999215 were linked to a translation initiation factor (IF2) and a nucleotidyltransferase (NTF), respectively. Markers S1_140803955 and S1_122448799 on chromosome 4 were linked to prolyl 4-hydroxylase subunit and phospholipase-like protein, respectively. Two markers, S1_326753690 and S1_526753692, located at the same locus on chromosome 6 were linked to a nodule cysteine-rich (NCR) secreted peptide. On chromosome 7, markers S1_212244210 and S1_352766048 were linked to poly(A) RNA polymerase (PAPD7) and ribosome biogenesis factor (LTV1). On chromosome 8, Markers S1_255810503, S1_251014656 and 328564258 were linked to late nodulin, NBS-LRR type disease resistance protein and serine/threonine-protein kinase (KIN4), respectively. Two markers, S1_276968267 and S1_276968408 at the same locus were linked to B3 DNA-binding domain protein.

To understand genetic factors by which biomass yield is affected by water stress, phenotypic data of biomass yield under well-watered and water stress treatments were applied separately for GWAS in combination with the same genotypic data set. Different results were obtained between the two treatments. A total of 28 significant markers were associated with biomass yield in plants under water deficit, while only four markers were significantly associated with the same trait under well-watered condition. Among them, all markers identified were different expect one (S1_276968267) on chromosome 3. This marker was significant in all three harvests under water deficit and in the third harvest of well-watered condition. Different sets of markers were identified by different treatments, suggesting that different mechanisms may exist to control the yield trait under water deficit and well-watered conditions. Although the candidate gene, B3 DNA-binding domain proteins were identified in both treatments, the genetic positions of the markers linked to this gene were differently located (one on chromosome 1 and another on chromosome 3). B3 DNA-binding domain protein is a transcription factor that has been reported to play an essential role in plant growth and development under cold stress (Yamasaki et al., 2004). However, a superfamily of B3 DNA-binding domain proteins exists in plants. The ones identified in the present study by different treatments may present different members of the family with different functions.

It has been reported that a number of genes and pathways involved in the yield and yield components under abiotic stress in plants (Mickelbart et al., 2015, for review). Our result of identification of multiple loci genome-wide supports that biomass yield under water stress in alfalfa is a complex trait with polygenic inheritance. Multiple loci associated with yield trait were consistent at different harvesting periods under water deficit. For instance, among 28 markers identified within the stress treatment, 18 and 6 were identified by three and two harvests, respectively. Only 4 were identified by single harvest. It may also suggest that the experimental data and statistical procedures used in the present study are reliable.

Our results indicated that yield under drought is controlled by many loci over all chromosomes and this has been found in previous studies as well (Zhang et al., 2015). In the previous report, we also used GWAS and identified a group of loci associated with drought resistance index, a similar trait to biomass under drought, and the loci identified in that study were also located throughout the 8 chromosomes of alfalfa (Zhang et al., 2015). Markers identified between the two studies are not identical, however, likely due to the differences between field conditions used in the present experiment and the greenhouse conditions used previously. Whether markers identified by greenhouse or field experiments are useful for MAS or not will need further validation, which we are going to perform using different platform such haplotyping and other technologies.

The polygenic control of biomass under drought in alfalfa has been reported by different group. Recently, Ray et al. (2015) identified 25 QTLs associated with biomass under drought in two alfalfa mapping populations. They were located on all chromosomes except chromosome 7. In the present study, we identified two loci on chromosome 7. This inconsistency may due to different mapping approaches having been applied in our and Ray et al.'s studies. Ray et al. used a traditional QTL mapping approach with 600 single-dose markers, whereas we used a whole genome approach with more than 4.6 thousand single nucleotide polymorphic (SNP) markers. The SNP markers have higher coverage because of their high abundance and widespread distribution throughout the genome. Moreover, drought tolerance is a complex trait and is affected by genetic and environmental interaction. It is not surprising to have different results from different mapping populations and in different environments.

Several marker loci identified in the present study were linked to genes playing roles in stress responses. Among them, the S1_20593175 locus on chromosomes 1 was linked to leucine-rich repeat receptor-like kinase. Leucine-rich repeat receptor-like protein kinases are plasma membrane proteins that are important in the transduction of various plant environmental and developmental signals (Tichtinsky et al., 2003). It has been reported that leucine-rich repeat receptor-like kinase 1 plays a role in the abscisic acid (ABA) signaling pathway in Arabidopsis (Osakabe et al., 2005). ABA is an important phytohormone that regulates plants' adaptive responses to abiotic stresses, such as drought and high salinity. It has been shown that ABA synthesis and accumulation are promoted by abiotic stress in plant cells, and ABA triggers stomatal closure in guard cells and regulates the expression of many genes with function in dehydration tolerance of plants (Finkelstein et al., 2002). Another drought responsive gene, a phospholipase-like protein, was linked to the S1_122448799 on chromosome 4. Phospholipases are involved in phospholipid-based signaling representing a component of drought-responsive signal transduction pathways (Hirayama et al., 1995; Sang et al., 2001). A serine/threonine-protein kinase was linked to S1_328564258 on chromosome 8. Plant-specific serine/threonine kinases are involved in plant response to abiotic stresses and abscisic acid signaling. It has been reported that the serine/threonine-protein kinase SOS2 is involved in the calcium-mediated regulation in plant salt tolerance in Arabidopsis (Guo et al., 2001; Zhu, 2016). Interestingly, Marker S1_251014656 on chromosome 8 was linked to a NBS-LRR type disease resistance gene. The plant NBS_LRR gene family contains a large class of disease resistance genes (DeYoung and Innes, 2006). The identification of a marker linked to the NBS-LRR genes in the present study suggests that it also involve in drought response in alfalfa. A similar finding of drought-related role for a NBS-LRR has also been reported in Arabidopsis where overexpression of the NBS–LRR gene ADR1 enhanced drought tolerance (Chini et al., 2004). It has been suggested that a signaling network exists between disease resistance and drought tolerance, and ADR1 may be involved in signal transduction in this network (Chini et al., 2004). The remaining genes had no direct function in drought tolerance, some may involve in the regulation of downstream drought tolerance genes. For instance, B3 DNA-binding domain protein linked to S1_276968267 and S1_276968408 on chromosome 3 may play a role in the transcriptional regulation for plant growth and development in response to abiotic stress (Yamasaki et al., 2004). Another example is the association of translation initiation factor IF2 with drought tolerance in the present analysis. It has been reported that phosphorylation of the α-subunit of translation initiation factor 2 is involved in the regulation of protein synthesis in response to a variety of environmental stresses (Jiménez-Díaz et al., 2013).

In conclusion, we used GBS for genome-wide genotyping and association and identified 22 marker loci associated with biomass yield under water deficit in the field in a panel of germplasm representing a diverse selection from the alfalfa gene pool. These marker loci were not presented in the well-watered control except one. The comparison of associated markers between control and stress treatments and among three sets of data obtained by three harvests showed that most of them were identical, indicating the consistency of the statistical analyses and stress treatments in the present study. The comparable results of the present study and the previous reports consistently demonstrated that loci associated with yield under drought were located throughout all 8 chromosomes, reflecting the polygenic nature of this complex trait. The BLAST search of the flanking sequences of the associated loci against DNA databases identified several stress-responsive genes linked to the drought associated loci, including serine/threonine-protein kinase, leucine-rich repeat receptor-like kinase, B3 DNA-binding domain protein, translation initiation factor IF2 and phospholipase-like protein. With further validation, those markers closely linked to drought resistance can be used for MAS to accelerate the development of new alfalfa cultivars with improved resistance to drought and other abiotic stresses.

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.