Skip to main content

Identification of genetic loci associated with forage quality in response to water deficit in autotetraploid alfalfa (Medicago sativa L.)



Alfalfa has been cultivated in many regions around the world as an important forage crop due to its nutritive value to livestock and ability to adapt to various environments. However, the genetic basis by which plasticity of quality-relevant traits influence alfalfa adaption to different water conditions remain largely unknown.


In the present study, 198 accessions of alfalfa of the core collection for drought tolerance were evaluated for 26 forage quality traits in a field trial under an imposed deficit irrigation gradient. Regression analysis between quality traits and water stress revealed that values of fiber-related traits were negatively correlated with values of energy-related traits as water deficit increased. More than one hundred significant markers associated with forage quality under different water treatments were identified using genome-wide association studies with genotyping by sequencing. Among them, 131 markers associated with multiple traits in all the water deficit treatments. Most of the associated markers were dependent to the levels of water deficit, suggesting genetic controls for forage quality traits were dependent to the stress treatment. Twenty-four loci associated with forage quality were annotated to functional genes that may play roles in cell development or in response to water stress.


This study addressed the genetic base of phenotypic variation of forage quality traits under water deficit. The SNP markers identified in this study will be useful in marker-assisted selection for the genetic improvement of alfalfa with enhanced drought tolerance while maintaining forage quality.


Alfalfa, “Queen of the Forages”, is the most productive and highest quality forage crop. Alfalfa quality has been determined by many factors, including protein, fiber and lignin contents, relative feed value, total digestible nutrients, and other physical and chemical factors. Alfalfa quality is directly influenced by its feeding value from animal performance. Fiber contents such as acid detergent fiber (ADF) and neutral detergent fiber (NDF) are important factors affecting the forage quality. Alfalfa forage contains 35–55% NDF, which contributes ~ 20–30% of the digestible energy value of alfalfa, the rest coming from non-fiber components. Relative feed value (RFV) is a tool that indexes alfalfa quality primarily based on its NDF content. The RFV index estimates digestible dry matter (DDM) of the alfalfa from ADF, and calculates the dry matter intake (DMI) potential from NDF. However, RFV has a significant shortcoming because it does not take into account how variations in NDF digestibility affect the energy content or intake potential of alfalfa. Even when harvested at an immature stage, the digestibility of alfalfa fiber can be very different [1]. In 2004, scientists at the University of Wisconsin designed another index, relative forage quality (RFQ) for estimating forage quality. The RFQ uses fiber digestibility and the total digestible nutrients of the forage to estimate intake [2]. The RFQ is an improvement over RFV as it better reflects the performance on animal fed. The RFQ emphasizes fiber digestibility while RFV uses digestible dry matter intake. The RFV continues to be widely used as an index to assess quality, compare forage varieties, and price forages. However, differences in the digestibility of the fiber fraction can result in a difference in animal performance when forages with a similar RFV index are fed.

Alfalfa forage quality and yield is affected by environmental factors, such as soil salinity and water supply. Alfalfa yield significantly declines when irrigation was not adequate [3]. Irrigation in an arid climate can affect nutritive value of alfalfa hay. Drought and high salinity are major factors that affect plant growth in the arid and semi-arid regions. Plants cope with these challenges by stress-avoidance or stress-tolerance. Stress-tolerant plants have evolved certain adaptive mechanisms within phenotypic plasticity to achieve different degrees of tolerance, whereas stress avoidance is the ability of plants to minimize the adverse effect.

The extent of phenotypic plasticity is primarily determined by genetic changes. Plants evolved mechanisms that facilitate adaptation to environmental changes under selective pressure [4]. It is unclear whether plant phenotypic plasticity is controlled by specific genes or a result of epistatic interaction during the selection of individual traits. However, genetic diversity and heterozygosity enhance adaptability to variable environments. It is important to identify plant functional traits and plasticity which will help plants adapt to global climate change [5].

Identification of molecular markers is the first step in marker-assisted breeding for genetic improvement. Single nucleotide polymorphism (SNP) is a type of markers that widely exist throughout the genome. One of the high-throughput and highly efficient approaches to discover SNPs is genotyping by sequencing (GBS), which was primarily used for phylogenesis and genome-wide association studies (GWAS) in maize [6]. It based on the high-throughput next generation sequencing [7]. Compared to the commercial SNP arrays, GBS has more advantages as its low cost, time saving and easy automation [8]. GWAS requires a large number of markers for mapping complex traits at the whole genome level. Over 15,000 SNPs were used to identify significant markers associated with cell wall biosynthesis and biomass yield in M. sativa [9]. In alfalfa, Li et al. [10] used GBS markers to construct a high-density linkage map representing high synteny between linkages of M. sativa and its wild relative, M. truncatula.

Unlike traditional genetic mapping, GWAS use thousands of SNPs throughout the genome to identify quantitative trait loci (QTL) associated with traits of interest using linkage disequilibrium [11]. Several factors can potentially influence GWAS power to identify significant associations, such as phenotypic variation, individual number, allele frequency, and population structure [12]. Mixed models have been used to correct the population structure and reduce the false positives of marker-trait association [13]. GWAS has been successfully used for mapping agronomic traits in some major crops, such as rice, maize, wheat, sorghum and soybean [14].

Alfalfa is an autotetraploid species (2n = 4X = 32) and alfalfa plants are highly heterozygous. It is a considerable challenge to develop markers with allele dosage in such a complex genome [15]. Since alfalfa cultivars are genetically broad-based synthetic populations, they provide an ideal system in which GBS, GWAS, and genomic selection (GS) can be applied. GWAS have been used for mapping quantitative trait loci associated with biomass yield, biotic and abiotic stresses in alfalfa [16]. Genetic markers associated with forage quality were identified by GWAS [17]. In addition, nearly 10,000 SNP markers were used for GS modeling, and further showed that GS increased genetic gain of biomass yield in alfalfa [18].

In the present study, we evaluated 26 forage quality traits in a panel of 198 alfalfa accessions of the core collection for drought tolerance obtained from the USDA-ARS Western Regional Plant Introduction Station. Majority of the accessions were collected in 1980s, in Canada and Northern US and dryland regions of other countries. The plants were evaluated in the field under three irrigation regimes: well-watered, mild and severe water deficits. To investigate the genetic base of the forage quality and its interaction with water stress treatments, we applied high-throughput genome sequencing GBS followed by GWAS, an integrated framework merged a QTL mapping approach to investigate genomic architecture of phenotypic plasticity of alfalfa quality traits under a gradient of water deficits. The ultimate goal is to identify genetic markers associated with forage quality traits in alfalfa under a deficit irrigation gradient, and use the closely linked markers for marker-assisted selection in breeding for high quality alfalfa varieties based on genetic potential and to reduce the confounding of environmental conditions with traditional breeding methods.


Phenotypic variations of forage quality traits

The analysis of variance (ANOVA) for 26 quality traits was carried out among the panel of germplasm and the result is presented in Table 1. The sum of squares varied from 0.29 in net energy of lactation (NEL) to 130,770.31 in RFV. The differences of most of the traits are statistically significant with the probability < 0.0001. Regression analyses of phenotypic variations in different water conditions showed significant effects of water deficit on forage quality (Fig. 1). The values of fiber-related traits, including ADF, aNDF, dNDF30, dNDF48 decreased as water deficit applied (Fig. 1a 1–4). Water deficit also decreased the contents of fat, RUP, IVDDM30 and NDFD48 (Fig. 1b 1–4), and slightly decreased TDNL, protein, IVDDM48 and lignin (Fig. 1c 1–4). Whereas the values for energy-related traits including: ME, NEL, TDN, NEM, NEG, RFV, ENE, DDM, and NFC increased by water deficit (Fig. 1d, e, and f). Water deficit slightly increased RFQ and ash (Fig. 1g 1–2).

Table 1 Analysis of variances of forage quality traits in the panel of 198 alfalfa accessions
Fig. 1
figure 1

Regression plots for 26 forage quality traits showing phenotypic variations under well-watered, mild and severe water deficits (X-axis, from left to right). Phenotypic values of each trait are presented by Y-axis

Phenotypic plasticity was estimated by calculating the plasticity index (PI) for each trait in the given accession under well-watered and water deficit conditions as described in the section of Methods. Overall, higher PIs were found in water deficit conditions compared to well-watered control except lignin, fat and protein (Fig. 2). Within the water stress treatments, higher PIs were found in the mild stress than those in the severe stress for most of the quality traits. Among them, the highest PIs were found in RFQ with 0.55, 0.77 and 0.75 for the control, mild and severe drought conditions, respectively. The lowest PIs appeared in DMM with 0.12, 0.18 and 0.16 for control, mild and severe drought, respectively. The rest of the traits showed higher PIs in severe drought compared to the mild treatment. The PI values were very similar between DMI and aNDF, ENE and NEL, and TDN and ME.

Fig. 2
figure 2

Phenotypic plasticity index of 26 alfalfa quality traits in response to well-watered (CTL), mild water deficit (Mild) and severe water deficit (Severe) conditions

The correlation coefficients (r) between fiber-related traits and energy-related traits were decreased as drought increased. Among them ADF, aNDF, RUP, lignin dNDF30 and dNDF48 were negatively correlated with other quality traits and the correlations were increased as drought increased (Fig. 3 a, b and c, blue panels). Positive correlations were found between energy-related traits (Fig. 3b and c, red panels on the bottom right). The correlations were increased as drought increased. Highest r values were found between these traits when the plants were under severe drought (Fig. 3c). However, no significant change was found between fat and any other traits by drought treatments (Fig. 3b and c).

Fig. 3
figure 3

Correlation coefficient among 26 quality traits under well-watered control (CTL), mild water deficit (Mild) and severe water deficit (Severe) in alfalfa

Cluster analysis of germplasm using forage quality traits

The mean values of 26 forage quality traits were used for cluster analysis. Two large clusters and 14 subclusters were classified as showing in Fig. 4. The first large cluster contained 8 subclusters. Most of germplasm in this cluster were collected from cultivars from US and Canada and their quality traits such as crude protein and RFV were relatively higher (Table S1), so we named it as the higher forage quality cluster (Fig. 4, top cluster). Two checks, Rambler and Saranac, susceptible to salt/drought are in this cluster (Fig. 4, subclusters 5 and 8). The bottom cluster was furtherly classified into 6 subclusters containing germplasm collected worldwide, including old cultivars and landraces with relatively lower forage quality (Fig. 4, bottom cluster). Three salt/drought resistance checks, Malone, Mesa Sirsa and Wilson are in this cluster (Fig. 4 subcluster 9). There was a trend that alfalfa germplasm with resistance to salt/drought had lower forage quality, while higher quality was found in the susceptible alfalfa germplasm.

Fig. 4
figure 4

A hierarchical cluster obtained using farthest neighbor method with phenotypic values of all quality traits evaluated in the present study. Accessions (PIs) were clustered into 2 clusters (High and low quality clusters) and 14 subclusters. The high quality cluster contains 8 subclusters with relatively higher quality. The low quality cluster contains 6 subclusters with relatively lower quality. Two checks, “Saranac” and “Rambler”, susceptible to salt/drought stress were clustered into the high quality cluster (Subclusters 5 and 8, respectively), and three drought/salt resistance checks, “Mesa Sirsa”, “Wilson” and “Malone” were clustered into the low quality cluster (Subcluster 9)

Genome-wide association for forage quality

The combinations of the filtered 10,327 GBS markers and phenotypic data of 26 quality traits were analyzed by GWAS using a mixed linear model. The profile of marker-trait association for well-watered control (A), mild (B) and severe (C) water deficits were illustrated in quantile-quantile plots (QQ) (Fig. 5). As illustrated in Fig. 5, a consistent difference between expected (X-axis) and observed (Y-axis) p-value across the whole genome was implied by deviation from the X = Y. The significances of marker-trait association were presented in negative log P-values on the Y-axis. Only a small number of true associations were shown among majority of unassociated SNPs. Overall, lower significance was obtained for all traits in the control except ash (Fig. 5a). Most significant marker-trait association were found in quality traits under the mild stress (Fig. 5b). The level of significances was lower by the severe stress compared to that of the mild stress (Fig. 5c).

Fig. 5
figure 5

Quantile-quantile plots of marker-trait association from GWAS for forage quality traits under well-watered (a), mild (b) and severe (c) water deficits in the alfalfa association panel. The expected (solid lines) against observed (dot lines) -log10 p-values are presented on X-axis and Y-axis, respectively. Each color curve represents a quality trait as showing at the bottom of the figures

Of 26 traits analyzed, most significant marker-trait associations were found in ash, NDFD48, dNDF30, dNDF48, NFC and TDNL. (Fig. 6 b, e, h, k, n and q). Similar association profiles were found in the severe stress, but the marker’s significances were lower than those under the mild stress (Fig. 6 c, f, i, l, o and r). Whereas, no or less significances were shown for the same markers under control condition (Fig. 6 a, d, g, j, m and p). Among markers identified, the highest significant markers were associated with ash and they were located on chromosomes 2, 6, 7 and unknown chromosome (U) (Fig. 6b). Significant markers associated with NDFD48 were also identified on same chromosomes (Fig. 6e) under mild stress but not in control and severe stress (Fig. 6 d and f). Significant markers were associated with dNDF30 and dNDF48 and they were located on chromosomes 1, 2, 3 and 8 (Fig. 6 h and k). Markers associated with NFC and TDNL were found in mild drought and they were located on chromosomes 2, 6 and 7 (Fig. 6, n, q).

Fig. 6
figure 6

Manhattan plots of marker-trait association of six most significant quality traits under well-watered control (CTL), mild water deficit (Mild) and severe water deficit (Severe). The X-axis presents chromosome positions of loci based on the reference genome of M. truncatula (Mt4.0, v1). The Y-axis shows negative log (P-values) of marker-trait association. Chromosome numbers were assigned and illustrated at the bottom of the figures

Common markers identified among multiple quality traits by different drought treatments

Despite different loci identified among quality traits, common markers were found among multiple traits (Table 2). Marker S1_110050725 on chromosome 4 identified in CTL for ADF was also significantly associated with other 10 traits including DDM, ENE, IVDMD30, IVDMD48, ME, NEG, NEL NEM, Protein and TDN (Table 2, top row). Similarly, markers S1_305729816 for DMI1 in CTL was also associated with 6 other traits: IVDMD30, IVDMD48, NDFD48, RFQ, TNDL and crude protein contents. Marker S1_231443201 identified in ADF shared its association with DDM, ENE, ME, NEG, NEL, NEM, TDN and TDNL. Four markers (S1_197238737–90) at the same locus on chromosome 6 and unknown marker S1_292679040 identified for ash in the mild stress were also associated with 6 other traits. Markers S1_351118210 and S1_276968305 identified in CTL for IVDMD48 and ash, respectively, were significantly associated with 5 other traits. Marker S1_174013573 identified for DMI1 were also associated with DMI, fat, RFQ and RFV in the severe stress. Eight markers associated with ash also associated with NDFD48, NFC and TDNL (Table 2). Additionally, nine, eighteen and fifty-five markers identified in four, three and two traits, respectively. The remaining markers were associated with one trait (Table 2, bottom part). Interestingly, most of the high significant markers with lower p value and higher R2 are among of common markers, suggesting these markers may have major effects on forage quality under drought. Among those, ten markers were associated with three traits (NDFD48, NFC and TDN). The p values of these markers ranged from 4.01E-08 (S1_21394491) to 5.79E-16 (S1_197238737) and the marker’s R2 ranged from 0.22 to 0.38, respectively (Table 2). To oversee the genetic architecture of the population under different treatments, we compared markers significantly associated with CTL, mild and severe drought treatments with those identified markers using mean values of all treatments. Among 68 markers identified in the control, 17 were also identified in the mean (Fig. 7a). Of 70 significant markers identified in the mild stress, only 10 were also identified in the mean (Fig. 7a). Among 67 markers identified in severe drought, 20 were also found in the mean (Fig. 7a). We have also compared the common markers identified among the three treatments directly. There were 3 common markers between each pair of treatments (Fig. 7b). Only 2 markers were found in all three treatments (Fig. 7b).

Table 2 Significant markers associated with forage quality traits under well-watered control (CTL), mild water deficit (Mild) and severe water deficit (Severe) in the panel of 198 accessions. Chromosome numbers are based on M. truncatula
Fig. 7
figure 7

A Vann chart of significant loci associated with forage quality resulting from GWAS for quality traits in alfalfa under well-watered (CTL), mild and severe water deficits compared with mean (a) and without mean (b). The numbers of significant loci identified under each treatment were compared with those of mean values of all treatments, showing the numbers of common (overlapped) and specific loci for different treatments

Identification of functional genes closer to the significant marker loci

Using the flanking sequences of the significant markers, we performed a BLAST search against the reference sequences to identify potential candidate genes underlying the significant marker loci. Of markers identified, 23 were found to be in close vicinity to known genes in the M. truncatula genome (Table 2). Among those identified under well-watered condition, marker S1_305729816 associated with DMI1 was in close vicinity to a programmed cell death. Marker S1_371197359 associated with ash was overlapped with reticulon-like protein B2. Marker S1_153379823 associated with NFC was adjacent to cadmium/zinc transporting ATPase. A number of markers were identified under mild water deficit, among them, marker S1_292679040 was adjacent to cytochrome P450 family protein; marker S1_300540244 was adjacent to UDP-glucosyltransferase; S1_311409163 was adjacent to exostosin; S1_62932631 was adjacent to helix loop helix DNA-binding domain protein; S1_309228882 was adjacent to carotenoid cleavage dioxygenase; S1_130145842 was adjacent to RING zinc finger protein; S1_294827153 was adjacent to auxin-binding protein ABP19a; S1_274178328 was adjacent to E3 ubiquitin protein ligase XBOS32; S1_319319782 was adjacent to IQ calmodulin-binding motif protein; and S1_380790483 was adjacent to Feronia receptor-like kinase. Three markers, S1_20593175, S1_368369179 and S1_92116984 associated with dNDF30 and dNDF48 were in close vicinity to LRR-receptor-like kinase, HASTY1 and OPT family oligopeptide transporter, respectively. Of those identified under severe water deficit, markers S1_108148173, S1_259198355 and S1_261323300 were adjacent to TIR-NBS-LRR resistance protein, interferon-induced guanylate-binding protein and peroxidase family protein, respectively. Another marker S1_310998070 associated with fat was adjacent to RNA-binding (RRM/RBD/RNP motif) protein.


Mild drought intends to decrease fiber content and improve digestibility in alfalfa

Production of alfalfa for high quality requires an understanding of how environmental and forage management practices influence crop growth and development. Causations exists between the environment, plant response, and nutritive value. In general, yield and forage quality are inversely related. Under regular management, alfalfa yield increases most rapidly in early spring and summer, whereas the quality decreases during the hot summer. Cutting frequently for high-quality forage always induce low yield [19]. In previous study, we identified 22 markers associated with alfalfa yield under different water deficit conditions [20]. Of those, 15 markers were identical with the quality-related markers identified in this study, indicating that these loci might participate in controlling both the yield and quality traits in alfalfa. Interestingly, 14 of them are associated with dNDF30, implying the correlation between digestible neutral detergent fiber and yield under drought stress. It was reported that digestible NDF of forage can decrease over 40% in the maturity phase [21]. Lignification process is becoming active during maturity, lignin, cellulose and other complex carbohydrates are enriched and bound together to form vascular and xylem tissue, support plant growth and nutrient transportation system [22]. However, Lignin is an essentially indigestible compound. Lignin in the cell wall reduces the digestibility of cellulose and hemicellulose by rumen microbes [23]. Given that alfalfa quality is negatively correlated with yield during maturity, it is likely to suggest that the loci associated with NDF may also affect alfalfa yield under water deficit.

Any factor that retards plant development tends to promote the maintenance of forage quality. If a plant is stressed during growth, a shorter, finer-stemmed, leafier alfalfa is often produced. On the other hand, high temperature accelerates growth, tends to have a negative impact on forage quality. Alfalfa is relatively drought tolerant because its deeper root systems allow alfalfa to absorb deep soil water and quickly recover from drought conditions. However, when transpiration exceeds water absorption, a stress is imposed on the plant influencing metabolism, development, growth, and ultimately yield. Water deficit promotes a reduction in vegetative growth and promote early maturity. It has been suggested that mild drought stress may be beneficial for forage quality as drought-stressed alfalfa will accelerate its shift to reproductive growth [24]. Furthermore, the greater proportion of leaves in short-term of drought stress improve forage feed quality and digestibility [24]. However, if drought stress has been too severe, and for an extended period, plant stress is permanent and may not be recovered.

Alfalfa fiber is consisted of three components: cellulose, hemicellulose and lignin. Increasing fiber content of a forage generally decreases its energy content. Of the fiber fractions, cellulose is the major compound digested by the animal while lignin is virtually indigestible in both the rumen and lower intestines. In our study, drought decreased significantly both ADF, aNDF and lignin, which in turn increased energy-related traits such as TDN, ENE, DDM, NFC, RFV and RFQ. Cell wall remodeling is a common response of plants to abiotic stresses. Cellulose content in cell wall was significantly reduced as biomass composition drastically altered under drought stress. Drought stress increased cellulose conversion rates by enzymatic saccharification, affecting cell wall structural rigidity. Under drought stress, both cell wall composition and the extent of cell wall plasticity significantly variated among genotypes. However, only weak correlations were found between different levels of drought resistance, suggesting their independent genetic control.

Genetic architecture of forage quality under well-watered and water deficit conditions

Among markers associated with forage quality under different irrigation episodes, a small number of the markers were in common between well-watered and water deficit conditions, while most of them responded dependently to the treatments (Fig. 6), suggesting their dependent genetic control. However, when phenotypic mean was used for GWAS, similar association patterns were found amongst energy-related traits, including DDM, TDN, ENE, ME, NEM, NEG and NEL, and traits of DMI, DMI1, RFV and RFQ. The genetic responses to mean values of these traits may suggest common genetic bases among them. This is logical as all these were energy-related traits.

In the GWAS, we only found nine associated markers that have consistent effects across water deficit treatments (Fig. 6b). The rests were differentially associated with respective treatments. Interestingly, about 2 folds of markers were associated with mild water deficit compared to those identified by severe water deficit (Fig. 6b), suggesting that mild water deficit affect genetic responses more favorably for forage quality traits than those in the severe stress and the control. In previous study, we identified 22 markers associated with biomass yield under water deficit in the same panel of accessions [20]. These markers were not found in the well-watered control. Similarly, in the present study, markers associated with forage quality under water deficit were not found in the control condition. In another study of genome-wide association and genetic selection on alfalfa quality traits, 10 SNP markers associated with acid detergent lignin (ADL), NDFD and CP were identified on chromosomes 1, 2, 4, 5, 7 and 8 in alfalfa [17]. Protein content is another important factor influencing alfalfa quality. In this study, we identified 14 markers associated with CP and RUP in all water conditions, indicating these markers may play common roles in protein-synthesis regardless water conditions. However, based on genetic positions, these markers were not overlapped with those reported previously [17]. This is probably due to that the genetic background of alfalfa germplasm and stress treatments used in two studies were different. Since drought tolerance is a complex trait and affected by genetic and environmental interaction (G x E), the allelic effect of associated causal variants may be influenced by the treatment of the stresses. Therefore, we cannot directly address whether conditionally neutral alleles accumulate genetic variation at a faster rate than constitutively expressed genetic variation. For example, the number of significant markers were significantly reduced when severe water deficit occurred compared to mild stress and well-watered control. This may indicate that the plants shut down some metabolic pathways to save energy to accomplish drought avoidance under severe drought stress. It was also possible that power of QTL detection was lower in severe drought conditions because of the lower variations.

Putative candidate genes associated with forage quality

Among 23 annotated genes associated with forage quality traits, three genes were identified under well-watered condition (Table 2). The programmed cell death (PCD) protein was associated with DMI1, CP, RFQ, TDNL, NDFD48, IVTDMD-30 and IVTDMD-48. PCD in plants is a crucial component of development and defense mechanisms. It is an important process during the secondary cell wall formation in plants [25]. Its associations with multiple traits in the present study suggested that PCD involved in regulating forage quality in multiple ways. Reticulon-like protein B2 (RTNLB2) was associated with ash and TDNL. It has been reported that the RTNLB2 is located in endoplasmic reticulum and plays a role in regulating receptor transport to plasma membrane in Arabidopsis [26]. Another putative candidate, cadmium/zinc transporting ATPase (cadA) was associated with NFC. The cadA is located in vacuole and involved in cadmium and zinc or cobalt transport and may have a detoxification function through a vacuolar sequestration in Arabidopsis [27]. Fourteen genes were identified under mild water deficit (Table 2). Of those, cytochrome P450 was associated with ash, dNDF48, NDFD48, NCF and TDNL. P450 family protein is a large enzymatic protein family in plants and play a role in plant development and biotic and abiotic stresses responses [28]. Many P450 monooxygenases such as cinnamate 4-hydroxylase (C4H) and ferulate 5-hydroxylase (F5H) are key enzymes in lignin biosynthesis [29]. Sakiroglu et al. (2012) used F5H as one of the candidate genes in the lignin synthesis pathway and analyzed the haplotypes in exons 1 and 2 of F5H in three subspecies of diploid alfalfa and found that exon 2 had greatest number of SNPs. However, weak association was found between exon 2 of F5H and cell wall biosynthesis [30]. UDP-glucosyltransferase (UGT) was associated with ash, NDFD48, NFC and TDNL. UGT plays a role in abscisic acid (ABA) homeostasis which regulates the plant response to environmental stresses such as drought, cold and salinity [31]. In previous study, another gene encoding sugar transferase, identified on chromosome 8, is also related to several quality traits, such as ADF, NDF and xylose [9]. A RING zinc finger protein (RZFP) was associated with ash, NFC and TDNL under mild water deficit. It has been reported that overexpression of AtRZFP enhanced salt and osmotic tolerance through enhancing ROSs scavenging, maintaining Na+ and K+ homeostasis, adjusting the stomatal aperture to reduce water loss, and accumulating soluble sugars and proline to adjust the osmotic potential [32]. Proline and soluble sugar contents were increased when overexpression AtRZFP in Arabidopsis [32]. An E3 ubiquitin protein ligase XBOS32 was associated with ash and NFC. The role of E3 Ub-ligase in controlling protein turnover has been suggested by modifying UPS-related proteins and contributes to nuclear proteome plasticity during plant responses to environmental stress signals [33]. An IQ calmodulin-binding motif protein was associated with ash. Gao et al. [34] reported that a gene of osa-miR369c encoding IQ calmodulin-binding motif protein affected the regulation of plant growth under several abiotic stresses such as temperature, drought and salinity in rice. The identification of calmodulin-binding proteins in the present study supports the assumption that this regulator is important player in response to abiotic stress through the calcium-signaling pathway [35]. Five genes were identified under severe water deficit (Table 2). Among them, the TIR-NBS-LRR protein was associated with ash and NFC. The plant TIR-NBS-LRR gene family contains a large class of disease resistance genes [36]. The identification of the TIR-NBS-LRR gene associated with drought in present investigation suggested that this gene may play a role in drought response. There is evidence to suggest that overexpression of the NBS–LRR gene ADR1 enhanced drought tolerance in Arabidopsis and the ADR1 may play a role in signal transduction in a cross-talk in signaling network between disease resistance and drought tolerance [37]. An interferon-induced guanylate-binding protein (IIGBP) was associated with ash and NFC. The IIGBP is a GTPase induced by interferon and plays a role in directing inflammasome subtype-specific responses and their consequences for cell-autonomous immunity against a wide variety of microbial pathogens [38]. A peroxidase family protein was associated with ash and NFC. The peroxidase responses are directly involved in the protection of plant cells against adverse environmental conditions. Several roles have been attributed to plant peroxidases in response to biotic and abiotic stresses. A type III peroxidase RCI3 participated in the induction of HAK5, which is a high-affinity uptake transporter of potassium [39]. Moreover, peroxidases may have a cell wall cross-linking activity during plant defense mechanisms [40]. An RNA-binding (RRM/RBD/RNP motif) protein was associated with fat under severe water deficit. RNA-binding proteins (RBP) play important roles in post-transcriptional gene regulation. Recent investigation of plant RBPs demonstrated that, in addition to their role in diverse developmental processes, they are also involved in adaptation of plants to various environmental conditions [41]. Although the remaining genes identified under water deficit do not have direct roles in stress responses, they involve in diverse processes in cell developments. For instance, Auxin-binding protein (ABP) was associated with ash, NFC and TDNL under mild water deficit. It has been suggested that ABP1 in Arabidopsis is involved in a broad range of cellular responses to auxin, acting either as the main regulator of the response, such as interface for entry into cell division or, as a fine-tuning device as for the regulation of expression of early auxin response genes [42].


In the present study, we evaluated 26 forage quality traits in a panel of 198 alfalfa accessions in the field trial under deficit irrigation gradience. Our results showed that water deficit decreased fiber contents and enhanced energy-related traits. The highest correlation coefficient was obtained between RFQ and the quality mean, supporting that the RFQ is more accurate in estimating overall forage quality compared to the RFV. Only a small number of markers were commonly associated with all treatments. Most of the associated markers were dependent on water deficit treatments, suggesting diverse genetic controls for forage quality traits in different levels of drought stress. Although GWAS on forage quality have been reported, we are the first to address the genetic base of phenotypic plasticity of forage quality traits under water deficits. The information gained from the present study will be useful for the genetic improvement of alfalfa by enhancing drought/salt tolerance while maintaining forage quality.


Plant materials

A panel of germplasm composited of 198 alfalfa accessions with potential drought tolerance were obtained from the USDA-ARS Western Regional Plant Introduction Station. Majority of the germplasms were collected in 1980s, in Canada and Northern US including British Columbia, Saskatchewan, Manitoba, Idaho, Montana, Nebraska, Washington, and North and South Dakota. The objective of the initial collecting project was to sample alfalfa stands that had survived 25 or more years in drought stressed environments. The remaining accessions were from different countries, including twelve collected from Afghanistan, two from China and Russia, and one from each of the following countries, Algeria, Bulgaria, India, Lebanon, Germany, Spain, Turkey, Oman and Yemen (Table S1). All plant accessions used in this study were provided by the USDA-ARS Western Regional Plant Introduction Station.

Field experiments

Field experiments were conducted as previously described (Yu, 2017). Briefly, single representative plants of each of the 198 alfalfa accessions were clonally propagated. The cloned plants were transplanted to the field of the Roza Farm at the Irrigated Agriculture Research and Extension Center, Washington State University, Prosser WA, in 2016. Well-watered control, mild and severe drought stresses were applied to the plots as described (Yu, 2017).

Alfalfa sampling and forage quality measurements

Shoot samples were collected from the field plots with each water stress treatment and they were subjected to quality analyses. Plant samples were dried in oven at 60 °C. They were then ground in Wiley Mill (Thomas Scientific, US) prior to the final grinding in Cyclotec 1093 sampling mill (Foss, Hillerød, Denmark) through a 1 mm screen. Sample powders were loaded and measured by Near Infrared Reflectance Spectroscopy (NIRS). Spectra were collected by a scanning monochromator (FOSS NIR Systems 6500, Silver Spring, MD, USA) in the spectral range from 400 to 2500 nm. A published NIRS Consortium equation 13AH50.2-eqa (NIRS Forage and Feed Testing Consortium, was used to predict quality factors.

Statistical analysis

Phenotypic data were subjected to an analysis of variance with the random effect of genotype and the fixed effect of drought treatment. The regression plots (Fig. 1) and correlation coefficient (Fig. 3) were obtained using JMP13 (SAS Institute, NJ). Cluster analysis was performed using a nonlinear mapping method to investigate the relationships among 198 accessions using the combination of quality traits in the field experiments. Correlation analysis was done between the traits evaluated using the JMP Genomics (SAS Institute, NJ).

To estimate phenotypic plasticity, a plasticity index was calculated according to Valladares et al., [43] as follow:

$$ PI=\left({M}_{max}-{M}_{min}\right)/{M}_{max} $$

Where PI is plasticity index, Mmax is the highest value of the treatment average and Mmin is the lowest value of treatment average for a specific trait in the population.

Genotyping by sequencing

Leaf samples were collected from individuals and used for DNA extraction using the Qiagen DNeasy 96 Plant kit, according to the manufacture’s protocol (Qiagen, CA). A methylation sensitive restriction enzyme, EcoT221 was used for DNA digestion, followed by library construction. Genomic libraries was sequenced using Illumina HiSeq2000. FastQC (v0.11.2) was used for initial quality check of the sequence reads ( Process_Radtags built in Stacks was used for deconvoluting and cleaning sequencing reads [44]. The resulting reads with high quality were then aligned to the M. truncatula reference genome (Mt4.0 v1) ( truncatula) using the Burrow Wheelers Alignment tool (Version 0.5.9) with default alignment parameters [45]. Loci with missing > 50%, MAF < 5% were removed. After filtering, 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.

Genome-wide association analysis

The filtered marker data were used for GWAS using TASSEL according to Bradbury et al. [46]. A mixed linear model was used for GWAS as previously described [20]. The Benjamini false discovery rate (FDR) of 0.05 was used as a threshold for identifying significant marker-trait association [47].

Availability of data and materials

All data generated or analyzed in this study are included in this article and the supplemental files. The raw GBS data were submitted to the NCBI database with the bioproject ID: PRJNA287263.



Abscisic acid


Acid Detergent Fiber


Neutral Detergent Fiber analyzed with amylase


analysis of variance


Crude Protein




Dry Matter


Digestible Dry Matter


Dry Matter Intake using NDF


Dry Matter Intake using NDF and NDFD


Estimated Net Energy


False discovery rate


Genotyping by sequencing


Genomic selection


Genome-wide association studies


30-h In Vitro Digestible Dry Matter


48-h In Vitro Digestible Dry Matter


Metabolizable Energy


Net Energy for Maintenance


Net Energy for Gain


30-h Digestible NDF


48-h Digestible NDF


48-h NDFD


Net Energy for Lactation


Nonfibrous Carbohydrates


Near-infrared spectroscopy


Programed cell death


Plasticity index


Quantitative trait loci


Relative Forage Quality Index


Relative Feed Value Index


Rumen Undegradable Protein


Single nucleotide polymorphism


Trait Analysis by aSSociation, Evolution and Linkage


Total Digestible Nutrients


Total Digestible Nutrients for legume


  1. Goeser JP, Combs DK. An alternative method to assess 24-h ruminal in vitro neutral detergent fiber digestibility. J Dairy Sci. 2009;92(8):3833–41.

    Article  CAS  PubMed  Google Scholar 

  2. Undersander D, Moore JE. Relative forage quality. Focus on Forage. 2002. Accessed 2002.

  3. Guitjens JC. Alfalfa irrigation during drought. J Irrig Drain Eng. 1993;119(6):1092–8.

    Article  Google Scholar 

  4. Schlichting CD. The evolution of phenotypic plasticity in plants. Annu Rev Ecol Evol Syst. 1986;17:667–93.

    Article  Google Scholar 

  5. Gratani L. Plant phenotypic plasticity in response to environmental factors. Adv Botany. 2014.

  6. Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, Mitchell SE. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One. 2011;6(5):e19379.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Nadeem MA, Nawaz MA, Shahid MQ, Doğan Y, Comertpay G, Yıldız M, Hatipoğlu R, Ahmad F, Alsaleh A, Labhane N, Özkan H. DNA molecular markers in plant breeding: current status and recent advancements in genomic selection and genome editing. Biotechnol Biotechnol Equip. 2018;32(2):261–85.

    Article  CAS  Google Scholar 

  8. Poland JA, Rife TW. Genotyping-by-sequencing for plant breeding and genetics. The Plant Genome. 2012;5(3):92–102.

    CAS  Google Scholar 

  9. Sakiroglu M, Brummer EC. Identification of loci controlling forage yield and nutritive value in diploid alfalfa using GBS-GWAS. Theo Appl Genet. 2017;130(2):261–8.

    Article  CAS  Google Scholar 

  10. Li X, Wei Y, Acharya A, Jiang Q, Kang J, Brummer EC. A saturated genetic linkage map of autotetraploid alfalfa (Medicago sativa L.) developed using genotyping-by-sequencing is highly syntenous with the Medicago truncatula genome. G3-Genes Genom Genet. 2014;4(10):1971–9.

    CAS  Google Scholar 

  11. Scheben A, Batley J, Edwards D. Genotyping-by-sequencing approaches to characterize crop genomes: choosing the right tool for the right application. Plant Biotechnol J. 2017;15(2):149–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Korte A, Farlow A. The advantages and limitations of trait analysis with GWAS: a review. Plant Methods. 2013;9(1):29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Zhang Z, Ersoz E, Lai CQ, Todhunter RJ, Tiwari HK, Gore MA, Bradbury PJ, Yu J, Arnett DK, Ordovas JM, Buckler ES. Mixed linear model approach adapted for genome-wide association studies. Nat Genet. 2010;42(4):355.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Alqudah AM, Sallam A, Baenziger PS, Börner A. GWAS: fast-forwarding gene identification and characterization in temperate cereals: lessons from barley–a review. J Adv Res. 2020;22:119.

    Article  PubMed  Google Scholar 

  15. Harper AL, Trick M, Higgins J, Fraser F, Clissold L, Wells R, Hattori C, Werner P, Bancroft I. Associative transcriptomics of traits in the polyploid crop species Brassica napus. Nat Biotechnol. 2012;30(8):798.

    Article  CAS  PubMed  Google Scholar 

  16. Hawkins C, Yu LX. Recent progress in alfalfa (Medicago sativa L.) genomics and genomic selection. Crop J. 2018;6(6):565–75.

    Article  Google Scholar 

  17. Biazzi E, Nazzicari N, Pecetti L, Brummer EC, Palmonari A, Tava A, Annicchiarico P. Genome-wide association mapping and genomic selection for alfalfa (Medicago sativa) forage quality traits. PLoS One. 2017;12(1):e0169234.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  18. Li X, Wei Y, Acharya A, Hansen JL, Crawford JL, Viands DR, Michaud R, Claessens A, Brummer EC. Genomic prediction of biomass yield in two selection cycles of a tetraploid alfalfa breeding population. Plant Genome. 2015;8(2).

  19. Orloff S, Putnam D. Cutting schedule strategies to maximize returns. 2006. Accessed 11 Dec 2006.

  20. Yu LX. Identification of single-nucleotide polymorphic loci associated with biomass yield under water deficit in alfalfa (Medicago sativa L.) using genome-wide sequencing and association mapping. Front Plant Sci. 2017;8:1152.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Hoffman PC, Lundberg KM, Bauman LM, Shaver RD. The effect of maturity on NDF digestibility. Focus on Forage. 2003;5(15):1–3.

    Google Scholar 

  22. Barros J, Serk H, Granlund I, Pesquet E. The cell biology of lignification in higher plants. Ann Bot. 2015;115(7):1053–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Jung HJ, Samac DA, Sarath G. Modifying crops to increase cell wall digestibility. Plant Sci. 2012;185:65–77.

    Article  PubMed  CAS  Google Scholar 

  24. Cassida K. Managing alfalfa during drought. Michigan State University Extension 2012. Accessed 20 Jul 2012.

  25. Serk H, Gorzsás A, Tuominen H, Pesquet E. Cooperative lignification of xylem tracheary elements. Plant Signal Behav. 2015;10(4):e1003753.

    PubMed  PubMed Central  Google Scholar 

  26. Lee HY, Bowen CH, Popescu GV, Kang HG, Kato N, Ma SS, Dinesh-Kumar S, Snyder M, Popescu SC. Arabidopsis RTNLB1 and RTNLB2 reticulon-like proteins regulate intracellular trafficking and activity of the fls2 immune receptor. Pant Cell. 2011;23:3374–91.

    Article  CAS  Google Scholar 

  27. Mishra S, Mishra A, Küpper H. Protein biochemistry and expression regulation of cadmium/zinc pumping ATPases in the hyperaccumulator plants Arabidopsis halleri and Noccaea caerulescens. Front Plant Sci. 2017;8:835.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Xu J, Wang XY, Guo WZ. The cytochrome P450 superfamily: key players in plant development and defense. J Integr Agri. 2015;14:1673–86.

    Article  CAS  Google Scholar 

  29. Gou M, Ran X, Martin DW, Liu CJ. The scaffold proteins of lignin biosynthetic cytochrome P450 enzymes. Nat plants. 2018;4(5):299–310.

    Article  CAS  PubMed  Google Scholar 

  30. Sakiroglu M, Sherman-Broyles S, Story A, Moore KJ, Doyle JJ, Brummer EC. Patterns of linkage disequilibrium and association mapping in diploid alfalfa (M. sativa L.). Theo Appl Genet. 2012;125:577–90.

    Article  CAS  Google Scholar 

  31. Dong T, Xu XY, Park Y, Kim DH, Lee Y, Hwang I. Abscisic acid uridine diphosphate glucosyltransferases play a crucial role in abscisic acid homeostasis in Arabidopsis. Plant Physi. 2014;165:227–89.

    Article  CAS  Google Scholar 

  32. Zang D, Li H, Xu H, Zhang W, Zhang Y, Shi X, Wang Y. An Arabidopsis zinc finger protein increases abiotic stress tolerance by regulating sodium and potassium homeostasis, reactive oxygen species scavenging and osmotic potential. Front Plant Sci. 2016;7:1272.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Serrano I, Campos L, Rivas S. Roles of E3 ubiquitin-ligases in nuclear protein homeostasis during plant stress responses. Front Plant Sci. 2018;9:139.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Gao P, Bai X, Yang L, Lv D, Li Y, Cai H, Ji W, Guo D, Zhu Y. Over-expression of Osa-MIR396c decreases salt and alkali stress tolerance. Planta. 2010;231:991–1001.

    Article  CAS  PubMed  Google Scholar 

  35. Ranty B, Aldon D, Galaud JP. Plant calmodulins and calmodulin-related proteins: multifaceted relays to decode calcium signals. Plant Signal Behav. 2006;1:96–104.

    Article  PubMed  PubMed Central  Google Scholar 

  36. DeYoung BJ, Innes RW. Plant NBS-LRR proteins in pathogen sensing and host defense. Nat Immunol. 2006;7:1243–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Chini A, Grant JJ, Seki M, Shinozaki K, Loake GJ. Drought tolerance established by enhanced expression of the CC–NBS–LRR gene, ADR1, requires salicylic acid, EDS1 and ABI1. Plant J. 2004;38:810–22.

    Article  CAS  PubMed  Google Scholar 

  38. Kim BH, Chee JD, Bradfield CJ, Park ES, Kumar P, MacMicking JD. IFN-induced guanylate binding proteins in inflammasome activation and host defense. Nat Immunol. 2016;17:481–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Kim MJ, Ciani S, Schachtman DP. A peroxidase contributes to ROS production during Arabidopsis root response to potassium deficiency. Mol Plant. 2010;3(2):420–7.

    Article  CAS  PubMed  Google Scholar 

  40. Chen EL, Chen YA, Chen LM, Liu ZH. Effect of copper on peroxidase activity and lignin content in Raphanus sativus. Plant Physi Bioch. 2002;40:439–44.

    Article  CAS  Google Scholar 

  41. Lorkovic ZJ. Role of plant RNA-binding proteins in development, stress response and genome organization. Trends Plant Sci. 2009;14:229–36.

    Article  CAS  PubMed  Google Scholar 

  42. Tromas A, Perrot-Rechenmann C. Recent progress in auxin biology. Comptes Rendus Biologies. 2010;333(4):297–306.

    Article  CAS  PubMed  Google Scholar 

  43. Valladares F, Martınez-Ferri E, Balaguer L, Perez-Corona E, Manrique E. Low leaf-level response to light and nutrients in Mediterranean evergreen oaks: a conservative resource-use strategy? New Phytol. 2000;148:79–91.

    Article  CAS  PubMed  Google Scholar 

  44. Catchen JM, Amores A, Hohenlohe P, Cresko W, Postlethwait JH. Stacks: building and genotyping loci de novo from short-read sequences. G3: genes. Genomes, Genetics. 2011;1(3):171–82.

    CAS  Google Scholar 

  45. Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, Buckler ES. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633–5.

    Article  CAS  PubMed  Google Scholar 

  47. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc. 1995;57:289–300.

    Google Scholar 

Download references


We thank Mrs. Martha Rivera for her technical assistance and data collection.


This work was supported partially by The United State Department of Agriculture NIFA_AFRP Grant Number 2015–70005-24071 and the Agriculture Research Service base fund.

Author information

Authors and Affiliations



LXY planned and designed the research, did genotyping and analyzed the data. BB and SF performed field work and forage quality tests. LXY, SL, JH and SN wrote the manuscript. LXY, SL prepared the figures. LXY, SL and CCM revised the manuscript. The authors read and approved the final manuscript

Corresponding author

Correspondence to Long-Xi Yu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

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.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1

. Alfalfa accessions used in GWAS for salt tolerance.

Additional file 2.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lin, S., Medina, C.A., Boge, B. et al. Identification of genetic loci associated with forage quality in response to water deficit in autotetraploid alfalfa (Medicago sativa L.). BMC Plant Biol 20, 303 (2020).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: