Identification of genetic loci associated with five agronomic traits in alfalfa using multi-environment trials

The use of multi-environment trials to test yield-related traits in a diverse alfalfa panel allowed to find multiple molecular markers associated with complex agronomic traits. Yield is one of the most important target traits in alfalfa breeding; however, yield is a complex trait affected by genetic and environmental factors. In this study, we used multi-environment trials to test yield-related traits in a diverse panel composed of 200 alfalfa accessions and varieties. Phenotypic data of maturity stage measured as mean stage by count (MSC), dry matter content, plant height (PH), biomass yield (Yi), and fall dormancy (FD) were collected in three locations in Idaho, Oregon, and Washington from 2018 to 2020. Single-trial and stagewise analyses were used to obtain estimated trait means of entries by environment. The plants were genotyped using a genotyping by sequencing approach and obtained a genotypic matrix with 97,345 single nucleotide polymorphisms. Genome-wide association studies identified a total of 84 markers associated with the traits analyzed. Of those, 29 markers were in noncoding regions and 55 markers were in coding regions. Ten significant SNPs at the same locus were associated with FD and they were linked to a gene annotated as a nuclear fusion defective 4-like (NFD4). Additional SNPs associated with MSC, PH, and Yi were annotated as transcription factors such as Cysteine3Histidine (C3H), Hap3/NF-YB family, and serine/threonine-protein phosphatase 7 proteins, respectively. Our results provide insight into the genetic factors that influence alfalfa maturity, yield, and dormancy, which is helpful to speed up the genetic gain toward alfalfa yield improvement.


Introduction
Alfalfa (Medicago sativa L.) is an autotetraploid (2n = 4x = 32) forage crop with a non-redundant haploid genome size of 816 Mb (Shen et al. 2020). Alfalfa is grown for hay, silage, and pasture due to its high nutritive value beneficial to livestock performance. Breeding of alfalfa is complicated and slowed by its tetraploid genome, outcrossing pollination (exhibiting severe inbreeding depression), and high heterozygosity, which prevents the creation of inbred lines. Alfalfa breeding programs have been focused on disease resistance and winter hardness. Genetic improvement for yield has been limited due to its complexity and polygenic inheritance and interaction with environmental factors (Putnam et al. 2007). To date, the development of alfalfa new varieties relies on phenotypic recurrent selection in field environments. Breeding process by phenotypic 1 3 121 Page 2 of 23 recurrent selection takes long periods to successfully develop new varieties (Bertrand et al. 2018).
Recurrent selection methods would be most effective when integrated with marker-assisted selection (MAS). MAS is a promising breeding strategy for improving complex traits (Xu and Crouch 2008). However, identification of marker loci tightly linked to target traits is the first step toward MAS. The high-density single nucleotide polymorphism (SNP) markers at the whole genomic level provide new platforms toward genome-wide association studies (GWAS) . Previous studies on GWAS in alfalfa have identified genetic markers associated with yield (Sakiroglu and Brummer 2017), plant height (Robins et al. 2007), plant branching (Wang et al. 2020), drought tolerance (Yu 2017), salt tolerance , and forage quality (Lin et al. 2021).
Multi-environment trials (MET) are crucial steps in plant breeding programs to estimate environmental effect. Phenotypic variation often occurs in the same genotype across environments, which is known as genotype-by-environment interaction (GEI). AS perennial forage crop, alfalfa yield is determined from multi-harvests per year for multiple years. To obtain the adjusted means of entries from MET, mixed linear model (MLM) is often used for estimation of GEI (Smith et al. 2005). Modeling for MET by MLM can use two approaches: single-stage and stagewise analyses. The singlestage analysis is the gold standard for MET data where all sources of variation are accounted simultaneously in a single model to estimate variance components and to calculate the best linear unbiased estimates (BLUEs) or best linear unbiased predictions (BLUPs) (Smith et al. 2001). Previous studies by comparing single-stage and stagewise analyses have found that both approaches are suitable for MET data analysis (Piepho et al. 2012;Damesa et al. 2019;Buntaran et al. 2020). However, single-stage analyses have several practical disadvantages including higher computational effort and complex variance-covariance structure for GEI. Additionally, single-stage model can be highly complex as it includes all specific details of spatial design in heterogeneous MET. Stagewise analysis can overcome these issues by improving computational efficiency (Möhring and Piepho 2009;Damesa et al. 2017). Stagewise analysis splits the analysis into first-stage or single-trial analysis and secondstage or joint analysis. In the first-stage analysis, each trial is analyzed separately to calculate BLUE values by accession focusing on the modeling of spatial variance. The BLUE values from the first-stage analysis are the input for the second stage of the stagewise analysis where accession means from all trials are subjected to a joint analysis to compute marginal means (Möhring and Piepho 2009).
The objective of this study was to identify genetic markers that influence five agronomic traits: biomass yield, plant height, fall dormancy, dry matter content, and maturity stage in alfalfa. To this end, we used MET and collected phenotypic data for five agronomic traits in a diverse panel of 200 alfalfa varieties/accessions planted in three locations from 2018 to 2020. Phenotypic data were calculated by MLM and used for GWAS in combination with a marker matrix from genotyping by sequencing (GBS). Marker-trait association identified a total of 84 markers associated with five traits. The annotation of associated markers reveals several putative candidate genes. The information obtained from this study provides new insight into genetic factors that influence biomass yield, maturity stage and fall dormancy of alfalfa.

Plant materials and field experiments
Plant materials were described previously (Lin et al. 2021). In brief, a diverse panel of 200 varieties and accessions from industry, and government were used. Collectively, S&W Seed Co., Alforex Seed Co., Legacy Seed Co, and Blue River Hybrid Seed Co. contributed with 52 varieties and 148 accessions were selected form the USDA-ARS National Plant Germplasm System. The diverse panel of germplasm was planted in three locations: Kimberly, Idaho (ID) (42°33′0.16″N 114°20′17.33″W), Union, Oregon (OR) (45°12′21.4″N 117°52′35.6″W) and Prosser, Washington (WA) (46°17′31.45″N 119°44′10.26″W). During the early stages of testing, augmented randomized complete block design (ARCBD) is preferred because selection decisions are based on data from multiple field trials with relatively large numbers of entries with limited replication (Lin and Poushinsky 1983). In ARCBD, a set of control cultivars are replicated in each block. The replicated controls were used to estimate block effects and error to adjust the values of each new genotype being tested. This design is especially useful in initially steps of breeding programs with many genotypes and limited seed for replications (Federer 1961). In this study, since seeds provided by USDA germplasm center are limited, we decided to perform ARCBD with 3 replicates in three locations (one replicate/location). The ARCBD had 11 blocks with different spatial distributions. Field trials in ID and WA consisted of four columns and 55 rows, while OR consisted of two columns and 110 rows (Table 1 and Table S1). Each block contained 20 plots with a plot size of 1 m by 4.6 m. Twenty entries were planted by block including two controls (Vernal and Hi-Gest360) and 18 new entries. The controls were replicated in each block.

Phenotypic data collection
Five traits including maturity stage measured as mean stage by count (MSC), dry matter content (DM), plant height (PH), biomass yield (Yi), and fall dormancy (FD) were measured during 2018-2020 with 3-5 harvests/year in three locations. (Table S2). One week before the harvest, the average PH in each plot was measured. A total 19 data sets of PH were collected. The yield (Yi) was measured by fresh weight of aboveground biomass a total of 31 times. Additionally, total yield was analyzed as the sum of all harvests by year by location. In each year, fresh biomass of the first cutting was dried in an oven at 60 °C for 5 days to obtain the DM as a proportion of dry weight (DW) over fresh weight (FW): Alfalfa maturity stage was measured by mean stage by count (MSC). MSC was obtained by sampling one handful of shoots by accession cut at the soil surface and the shoots were separated and counted belonging to each developmental stage from 0 to 9 (Fick and Mueller 1989). MSC was calculated as the average of the individual stage categories present in the sample, weighted for the number of the shoots at each stage: where S is the stage number (0-9); N is the number of shoots in stage S ; C is the total number of shoots in the sample.
One week before the harvest, PH was measured in cm from the basal part of the plant to the tip of the tallest plant for six times per plot without stretching the stems according to Van Riper and Owen (1964). Fall dormancy scores were obtained by measuring plant height in cm 25 days after the fall cutting. PH in the fall was transformed to fall dormancy rates using the intercept and slope of a linear regression model using standard FD control cultivars included in the panel: Maverick (FD 1), Vernal (FD 2), Hi-Gest 360 (FD 3), Legend (FD 4), Archer (FD 5), and ABI 700 (FD 6) ( Table S3). The MET dataset was analyzed by mixed linear models to obtain adjusted mean values (Table S4).

Single-trial analysis
Phenotypic traits were statistically analyzed using a single-trial and stagewise approaches. The responsive variable was fitted using a mixed linear model to obtain the best linear unbiased estimate (BLUE) values of all entries in a single-trial analysis using control cultivars as fixed effects, block position as random effects, and column and row position in the field as residuals: where y ij is a vector ( n i × 1 ) of response variable observed in ith accession at the jth block, is the overall mean, i is a vector of accessions as fixed effects, (c) ij is the fixed effect of the ith control accessions nested in the plot jth, b j is the random effect of the jth block, and ij ∼ N 0, 2 is the vector ( n i × 1 ) of residuals for the trial. Block and control plot effects are used to adjust response values. Errors derived from control plots were used to obtain valid differences among new entries. To improve the estimation of residuals field position coordinates of plots were used in a spatial modeling with an autoregressive model of order one (AR1 × AR1) as described by Stefanova et al. (2009). The residual vector was modeled as follows: where Σr and Σc are the r × r and c × c correlation matrices for the row and column dimensions, and ⊗ is the Kronecker product of correlation matrices to obtain a single autocorrelation parameter. Spatially correlated data assumes that the correlation between measurements increased as the separation between them decrease: (3)  Test   ID  22  198  2  198  11  55  4  10  OR  22  198  2  198  11  110  2  10  WA  22  198  2  198  11  55  4  10 where 1 , 2 and 3 are Pearson's correlation values between adjacent rows and columns and Pearson's values are greater in close positions ( 1 > 2 > 3 ) . Single-trial analysis allows to estimate harvest-specific means (BLUEs) and variances.

Stagewise analysis
Stagewise analysis used BLUE by accession obtained from multi-environment datasets to obtain best linear unbiased prediction (BLUP) values. The environment was defined as the unique interaction of location by phenotyping time (year and cut), presenting environmental variations such as temperature, daylength, and rainfall (Smith et al. 2001). BLUPs were obtained by comparing multiple variance-covariance structures for multi-environment datasets as random effects using ASReml-R software (Butler et al. 2018). The correlation structure for the effects of a common accessions across multiple environments was specified in a nested model: where y ij is the response variable observed in ith accession at the jth environment, is the overall mean; E j is the fixed effect of the environment jth, G(E) ij is a vector of random genetic effects of accession ith nested in the environment jth; assuming a distribution of G(E) ij identical and independently where G E is a s × s variance-covariance matrix across s environments and m is the identity matrix for m accessions. Various forms of G E may be used, for example, the simplest G E structure is DIAG_G E , which assumes different genetic variance for each environment and assumes the genetic effects are independent: where the diagonal elements ( 2 GE ) are given by the genetic variances for each environment explained by the GEI and the off-diagonal elements are the covariance between pairs of environments. A more realistic approach in MET is the unstructured ( US_G E ) variance-covariance matrix with no assumption of equal variances and/or covariances: where the diagonal elements were described above, and the covariances are different in each pair of environments. CORGH is equivalent to the US matrix, but it is ⎤ ⎥ ⎥ ⎥ ⎦ parametrized in terms of correlations instead of covariances between environments. Factor analytic (FA-k) variance models provide a good approximation to the unstructured variance model reducing the number of parameters required (Thompson et al. 2003). FA-k models are based on the number of the k factors used to model accurately the observed variance-covariance relationships among and within environments. FA-k models can be specified in ASReml-R software (Butler et al. 2018) where the variance was given as the direct product of an FA covariance matrix for environments ( FA_G E ) and an accession effect correlation matrix parameterized as: where is the matrix of environment loadings with dimension s × k ; is an off-diagonal matrix with site-specific genetic variances with dimension s × s ; ⊗ is the Kronecker product and m is the identity matrix for m accessions. Factor loadings are part of the outcome from FA, which serves as a data reduction method designed to explain the correlation between observed variables using a smaller number of factors. For example, in a FA1 model ( k = 1 ) with four environments ( s = 4 ) the and matrices are: where 1s is the factor loading for the environment s and the diagonal elements of are the site-specific genetic variances. The resulting FA_G E variance-covariance matrix becomes the direct product of an FA1 matrix for environmental effects: where the diagonal elements of FA_G E are the squared loadings plus site-specific variance (e.g., 2 11 + Ψ 1 ) and the genetic covariance between sites 1 and 2 is the product of their loadings ( 11 12 ). The FA1 model assumes a single latent variable explaining the covariance, FA2 assumes two latent variables, and so on. In comparison with US variance-covariance formulation, FA-k models use fewer parameters in a high number of environments: where p FA and p US are the number of parameters for FA-k and US models, respectively, s is the number of environments, and k is the number of factors (Thompson et al. 2003). FA-k models are convenient when the number of environments increases, avoiding model convergence and inaccurate or uncertain predictions.

DNA isolation and GBS library construction
DNA isolation and library construction were completed as previously described (Lin et al. 2021

Sequence analysis and genotyping
Raw sequence fastq files were aligned to the Medicago sativa (Zhongmu No.1) reference genome (Shen et al. 2020) using Bowtie2 (Langmead and Salzberg 2012). Next-Generation Sequencing Experience Platform (NGSEP) 4.0 was used for variant calling and filtering with the following standards: ploidy was set as 4; the maximum base quality score was set at 30; minimum allele frequency (MAF) greater than 0.05; a missing value of the samples genotyped lower than 20%; minimum genotyping quality was set at 30 (Duitama et al. 2014). After filtering, the variant call format (VCF) file was imputed using the hidden Markov model in NGSEP and transformed into GWASpoly format using the function VCFConverter of NGSEP software.

Population structure analysis
The genotypic data were used for principal component analysis (PCA). First, the vcf file containing 97,345 SNPs was transformed to bed format using the plink software (Purcell et al. 2007). Then, eigenvector file containing the data of PC1, PC2, and PC3 was generated using GCTA (Yang et al. 2011). The PC values were included as covariates in the analysis of GWAS.

Genome-wide association study and functional annotation
GWAS were performed using the R package GWASpoly with Q + K MLM. The Q + K MLM incorporates correlations among individuals including kinship matrix (K) and population structure (here, principal components) as covariates (Q) as follows (Rosyara et al. 2016): where y is a vector of observed phenotypes ( n × 1 ); X is an incidence matrix of fixed effects ( n × p ); is a vector of fixed effects ( p × 1 ); Z is an incidence matrix mapping observed phenotypes to genotypes (accessions) ( n × t ); S is a structure incidence matrix ( t × d ); is a vector of SNPs effects ( d × 1 ); Q is an incidence matrix for population size ( n × q ); v is a vector of fixed effects for the principal components ( q × 1 ); u is a vector of random effects of the mixed model ( t × 1 ), with a covariance proportional to K matrix, Var[u] = 2 g K ; K is the ( t × t ) kinship matrix inferred from genotypes and is a residuals vector (Rosyara et al. 2016). K matrix in MLM accounts for phenotypic covariance due to genetic relatedness to reduce the number of false positive associations. GWAS were performed with eight different models: general, diplo-general, additive, 1-dom-ref, 1-domalt, 2-dom-ref, 2-dom-alt, and diplo-additive (Rosyara et al. 2016). The − log 10 p-values were used for identifying SNPs significantly associated with the traits with the Bonferroni correction at the cutoff value of 5%. Significantly associated markers were annotated against the Uniprot100 database (Bateman 2019) using the reference transcriptome dataset for M. sativa ).

Linkage disequilibrium
LD between pairs of SNPs was estimated using squared allele-frequency correlations ( r 2 ) calculated with the R package ldsep (Gerard 2021). The SNPs with minor allele frequency lower than 5% were excluded in the estimates of LD. Decay of LD with distances on base pairs (bp) based on the M. sativa reference genome was evaluated by nonlinear model according to (Remington et al. 2001). The expected r 2 value under drift-recombination equilibrium was E r 2 = 1∕(1 + C) , where C = 4ad , where a is an estimated regression coefficient and d is the physical distance in bp. Assuming a low level of mutation and finite sample size n , the expectation becomes (Hill and Weir 1988): To evaluate consistency of LD, global LD decay and estimated LD values for each chromosome were defined as the distance at r 2 = 0.1.

Genotyping and SNP distribution and LD decay
Of the raw sequencing reads, 90.1% were successfully aligned to the Medicago sativa reference genome (Shen et al. 2020). After filtering, 97,345 SNPs were obtained with an average depth of 30×. Over 38,000 (39.04%) biallelic SNPs were localized in coding regions and 60.96% were in noncoding regions. The SNPs were not evenly distributed across chromosomes, as SNP density was higher close to the ends of the chromosomes, and lower close to the possible centromeric regions. A large (3 Mb) gap was found in the centromeric region of Chromosome 1. All chromosomes except chromosome 6 contained over 100 SNPs per megabase (Mb) on average, with the highest density observed on chromosomes 4 and 7 with 133.17 and 133.13 SNPs per Mb, respectively, while chromosome 6 contained 96 SNPs per Mb (Figure S1a). A high proportion of heterozygous markers (0.94) is shown in Figure S1b. The percentages of five allele types were AAAA = 2.73%, AAAB = 21.85%, AABB = 50.54%, ABBB = 21.75%, and BBBB = 2.71% ( Figure S1b).

Phenotypic variations
Five agronomic traits were collected in three locations over three years, multiple cutting times per year (Table S2). The results of statistics of MSC, DM, PH, FD, and Yi are presented in Table S4. Low values of MSC were found in WA and ID with the lowest value in WA_2020_5 (1.94). For DM, the highest mean value was found in WA_2018_1 (21.56%) and the lowest value in WA_2019_1 (17.01%). The highest mean value of PH was obtained in OR_2019_1 (104.14 cm) and the lowest mean values in WA_2020_5 (55.88 cm). The average PH of all accessions before the first cut was higher than those collected in the later cuts, which showed a similar trend as the yield data. Biomass yield was analyzed as single harvest and annual forage yield was the sum of all harvests by year by location. In 2018, the highest annual forage yield was found in OR (11.89 Mg ha −1 year −1 ) and lowest in ID (8.44 Mg ha −1 year −1 ). In 2019, the highest mean value of annual forage yield was found in ID (21.91 Mg ha −1 year −1 ) and lowest mean value in WA (16.07 Mg ha −1 year −1 ). In 2020 annual forage yield was collected in OR and WA with missing data from ID, the highest mean value was found in WA (17.08 Mg ha −1 year −1 ) and lowest in OR (15.02 Mg ha −1 year −1 ). Data analysis of biomass yield by single harvest showed the highest mean value for the first cut in ID in 2019 ID_2019_1 (7.37 Mg ha −1 ) followed by OR_2019_1 (6.59 Mg ha −1 ), while the lowest mean value was obtained in WA_2020_5 (1.04 Mg ha −1 ) and ID_2018_3 (2.07 Mg ha −1 ). Generally, yields of the first cuts were highest in the same year at the same location, then decreased gradually in the later cuts. FD is classified according to plant regrowth in early autumn into fall dormant (1-3), semi-dormant (4-6), and low-dormant (7-9) accessions. FD was different according to location, but percentages of FD grade were similar in ID and WA. 11% and 13.1% of accessions were classified to class 3, 78.9% and 80.5% to class 4, and 10% and 6.3% to class 5 in ID and WA, respectively. Accessions in OR have wider FD classification: 1.5%, 8.4%, 30.5%, 46.8%, 12.1% and 0.5% were classified into classes 1, 2, 3, 4, 5, and 6, respectively. In general, OR presented the highest values of FD compared with WA and ID.

Estimated performance in multi-environment trials
All traits collected were modeled by single-trial and stagewise approaches using the ASReml-R software. In the single-trial analysis, BLUEs were obtained using accessions as fixed effects, block position as a random effect, and column and row position as residuals. Three models were tested with different assumptions of residual variance: model one used autoregressive information of order one in rows and columns (AR1), model two used autoregressive information of order one just in rows (AR1-R), and model three used symmetric autoregressive information in rows and columns (SAR). Goodness-of-fit for single-trial models was compared by Akaike's information criteria (AIC) and BLUEs values (ST0) were obtained from the model with the lowest AIC. Model three (SAR) had the best performance, with the lowest AIC values in 38 datasets, followed by model two (AR1-R) and model one (AR1) with 24 and 23 datasets with the lowest AIC, respectively. Among those, model three was the best for most datasets in MSC, PH, Yi, and FD (Table S5).
The broad-sense heritability (H 2 ) was calculated in single-trial models according to the method of Lush (1940) ( Table S6). FD has medium to high H 2 in all environments (0.45-0.92), indicating that most of the phenotypic variations of FD were genetically controlled. DM, MSC, PH, and Yi had some environments with low H 2 (> 0.1). DM has low H 2 in WA_2018_1; MSC has low H 2 in OR_2018_1, WA_2019_1 and ID_2018_1; PH has low H 2 in ID_2019_1, OR_2018_1, OR_2019_1 and WA_2020_1; Yi has low H 2 in ID_2018_1 and OR_2018_2. Higher heritability was found in DM of WA_2020_1 (0.53), MSC of WA_2020_2 (0.53), PH of WA_2019_2 (0.85), and Yi of OR_2019_2 (0.70). Total annual yield has high heritability in OR_2019 (0.73) and lower heritability in ID_2018 (0.12).
In the stagewise approach, estimated performance was obtained as BLUP values using adjusted means from single-trial analysis and weights, calculated as reciprocal of the variances of the means. Five stagewise models including diagonal (DIAG), unstructured with correlations (CORGH), and factor analytic with k-number of factors modeled (FA1, FA2, FA3) were compared based on the number of estimated variance parameters and AIC (Table 2). Generally, FA-k models were computationally more efficient than models with a high number of variance parameters like CORGH. The highest AIC values were obtained with DIAG while FA2 and FA3 were the best models in all traits. FA3 showed the lowest AIC for the traits of MSC, PH, and Yi, while traits of DM and FD were better modeled with FA2. BLUP values were obtained from the model with the lowest AIC score. For each trait, BLUPs were predicted by accession by location by harvest (ST1), by accession by location (ST2), or by accession (ST3).
Pearson correlation analysis of ST1, ST2 and ST3 BLUP values showed that FD, PH, and Yi traits were highly correlated while DM while MSC presented medium to low Pearson correlations (Table S7 and Fig. 1). Overall, predicted means of PH and FD had the highest Pearson correlation (0.879), PH and Yi had a Pearson correlation of 0.773, and FD and Yi had a Pearson correlation of 0.737 (Fig. 1). GEI analysis showed a similar distribution for FD in ID and WA with a mean value of 3.92 (± 0.40) and 3.83 (± 0.37) while FD in OR has a mean value of 3.62 (± 0.81). PH and Yi were significantly different by location; PH had the highest mean value in ID (28.65 ± 2.67), followed by OR (25.39 ± 3.02), and WA (19.01 ± 2.30). Similarly, Yi had the highest mean value in ID (4.37 ± 0.45), followed by OR (4.05 ± 0.43), and WA (3.34 ± 0.65) (Fig. 1).

Genome-wide association studies
BLUE and BLUP values from single-trial and stagewise analyses were used to identify significant SNPs associated with agronomic traits in different environments. One hundred and forty-two significantly associated markers (Bonferroni correction > 0.05) were identified using all GWASpoly models and all traits. General and Diplo-general GWASpoly models had the highest number of significant markers with 33 and 29 markers, respectively. Twenty-one markers were identified by two or more GWASpoly models on the same trait and the same environment, these markers were considered duplicated markers (Table S8). We removed duplicate markers maintained 84 unique markers for further analysis. Of those, 55 markers were located in coding regions and 29 markers were in intergenic regions. Percentage of explained variance was calculated as R 2 for each marker and ten markers had R 2 values > 0.1 (Table S8). Markers with the most significant values were Chr2_50327189 (MSC), Chr3_92058964 (MSC), Chr5_16888188 (FD) with a − log 10 p-values of 8.96, 9.43 and 8.39, respectively. By comparing GWAS results between single-trial and stagewise, we identified most significant markers by ST1 (72). ST0 dataset identified 44 significant markers, while stagewise analysis ST2 and ST3 identified 25 and one significant markers, respectively. The use of ST1 BLUP values from MET data significantly increased the number of significant markers identified by GWAS (Table S9).

Markers associated with yield and plant height
In total, 13 yield-associated markers were identified in all datasets (Table S9, Figure S2 and Figure S3); these included four markers located on chromosome 3, one on chromosome 4, four on chromosome 5, three on chromosome 6 and one on chromosome 8. Markers Chr3_5642102, Chr5_81433854, and Chr6_93711284 were consistently identified two to seven times in different yield datasets (Table S9). Ten markers were in the coding regions, and they were annotated as functional proteins with different roles Table 2 Goodness-of-fit for the range of genetic variance models tested in multienvironmental trials Model tested were diagonal (DIAG), US-correlations (CORGH), and factor analytic with k-number of factors modeled (FA1, FA2, FA3). DF is the number of estimated variance parameters in the genetic variance matrix G by trait and AIC is Akaike's information criteria. Traits tested were maturity stage measured as mean stage by count (MSC), dry matter content (DM), plant height (PH), fall dormancy (FD), yield (Yi) and total yearly yield. Gray cells correspond to the best model chosen by the lowest AIC scores  (Table 3). Markers Chr6_37188399 and Chr8_63386124 were annotated as THOC4 and Tubby-like F-box with functions as a transcriptional coactivator and transcription factor, respectively. Markers Chr3_93233760, Chr5_19865639, and Chr6_79943661 were annotated as BABAM1, TRN2, and KEA5 involved in plant development. Marker Chr3_55269659 was annotated as TMV resistance gene associated with biotic resistance. Markers Chr3_5642102, Chr3_45718711, and Chr5_81433854 were annotated as FHY3, PP7 and Calcineurin-like phosphoesterase, respectively. These genes are activated in response to abiotic stress (Table 3).
In total, 20 PH-associated markers were identified in multiple datasets (Table S9 and Figure S4). Markers Chr3_67223769, Chr5_18432131, Chr7_3866977, Chr7_3918330, Chr7_3918343, Chr7_3918348, and Chr8_54366238 were identified in two harvests. Markers Chr4_8551406 and Chr6_15133674 were identified in three harvests (Table S9). Thirteen out of 20 markers associated with PH were in coding regions and annotated as proteins with different roles. Markers Chr4_79237342 and Chr4_79237343 were located at the same locus annotated as transcriptional regulator GRWD1. Chr4_50606506 and Chr8_63386102 were in loci annotated as NFYB/HAP3 and Tubby-like F-box with functions as transcription factors. Markers Chr2_49623557 and Chr6_15133674 were at loci annotated as O-FucT and F-box/LRR with roles in cell adhesion and plant development, respectively. Markers Chr1_11444898 and Chr8_25420777 were at loci annotated as GMI1 and SURPG related to DNA and RNA biogenesis. These markers had relatively higher values

Markers associated with dry matter content and maturity stage
Five markers were significantly associated with DM (Table S9 and Figure S5). Markers Chr2_34932941 and Chr2_34932942 were adjacently located in a noncoding region. Similarly, Chr1_91305390 was in a noncoding region as well. Marker Chr4_62801591 explained 10% of phenotypic variation of DM in ID 2019 and was located at the locus annotated as vesicle transport protein GOT1B-like. Marker Chr4_42520012 was in a coding region at the locus annotated as RGA2 (Table S8 and Table 3). Eighteen markers were significantly associated with MSC (Table S9 and Figure S6). Four markers, Chr2_56825927, Chr4_5209571, Chr2_63316866, and Chr1_99859764, explained more than 10% phenotypic variations (Table S8). Eleven out of 18 markers associated with MSC were in coding regions and annotated as proteins with different roles. Markers Chr1_53838421, Chr7_89237636, and Chr8_37920842 were at loci annotated as C3H1, DDX42, and PAT1, respectively. These genes are involved in DNA and RNA metabolism. Chr5_73969365 was located at the locus annotated as GIP1 that may act as a coactivator to transcription factors. Marker Chr2_63316866 was located at the locus annotated as PEPC (Table 3).

Markers associated with fall dormancy
Fall dormancy datasets were analyzed using two different approaches: 1) raw plant height measured 25 days after the fall cutting was analyzed and 2) PH data were transformed to FD classes using a linear regression with standard controls with known FDs (Table S3). Datasets from both approaches used the mixed linear models to obtain adjusted means and were subsequently analyzed by GWAS. More markers were identified using the first approach compared to the FD regression. Only three unique markers were associated with FD in WA_2018_3 (Chr6_989737, Chr6_37344769, Chr7_93097726), with higher − log 10 p-values in the second FD dataset, while 28 markers are significant in the first approach (Table S10). Therefore, we used markers from the first approach as markers associated with FD to be annotated (Table 3 and Figure S7).
Twenty out of 28 markers were in coding regions. Two markers Chr1_24328719 and Chr7_93097726 were located at two loci annotated as BAHD. Marker Chr4_11411532 was located at the locus annotated GST. Markers Chr6_989737, Chr6_106457549, and Chr8_34732089 were located at loci annotated as RST1, PPIase, and PHB, respectively. These genes are involved in plant development and stress response. Ten SNPs significantly associated with FD were located at the same locus annotated as NFD4 gene (Fig. 2a, b). The NFD4 gene codes for a protein named Nuclear Fusion Defective 4-like protein containing 573 amino acids with a nodulinlike domain in the N-terminal region (Emanuelsson et al. 2000;Portereiko et al. 2006). According to our alfalfa pan-transcriptome report , the NFD4 gene has 23 isoforms highly conserved from the position 16,888,888,195, henceforth named by their three last positions: 158 to 195, respectively (Fig. 3a). In our dataset, five SNPs were located within the coding region of NFD4 with five missense substitutions: (158) S135A, (163) L136F, (172) A139S, (188) N145D, and (195) V147A (Fig. 3b). Additionally, different alleles of  were associated with different FD scores. For example, Chr5_16888158, genotypes with allele type 1 (AAAB) had a FD mean of 5.07 (± 0.14) while genotypes with allele types 2 or 3 (AABB or ABBB) have a FD mean of 3.79 (± 0.34) and 3.81 (± 0.34), respectively (Fig. 2c).

Linkage disequilibrium analysis of significant markers
The functional relationship between LD (estimated as the square of the correlation coefficient ( r 2 ) between pairs of SNPs) and physical distances on the M. sativa genome was determined by fitting a nonlinear model. The extent of LD rapidly decayed to 23.3 kb at r 2 = 0.1 across all the alfalfa accessions but varied among chromosomes (Fig. 4). The average r 2 ranged between 108.9 and 0.22 kb for Chr4 and Chr8, respectively (Fig. 4).
LD analysis among markers associated with DM, fall dormancy, maturity stage, plant height, and biomass yield were performed using adjacent SNPs in a 10 kb window by Haploview v4.2 (Barrett et al. 2005). Twenty-two blocks were identified to harbor multiple markers at the same locus. Chromosomes 2 and 3 had one LD block, chromosomes 6 and 7 had two LD blocks, chromosomes 1 and 4 had three LD blocks, and chromosomes 8 and 5 had four and six LD blocks, respectively ( Figure S8). Marker Chr3_95636007 associated with PH was in a linkage group of 12 SNPs in a region of 57 nucleotides. Markers Chr5_10731696 and Chr5_10731717 associated with Yi were in a linkage group of 7 SNPs in a region of 984 nucleotides. Marker Chr5_40950448 associated with MSC was in a linkage group of 8 SNPs in a region of 4502 nucleotides. Marker Chr8_37920842 associated with MSC was in a linkage group of 14 SNPs in a region of 205 nucleotides. Markers Chr8_63386102, Chr8_63386124, and Chr8_63400012 associated with PH and Yi were in a linkage group of 5 SNPs in a region of 13,910 nucleotides.

Stagewise analysis increases the number of significant markers
Multi-environment trials are crucial steps in plant breeding programs to estimate genotypic values across a wide range of environments. The analysis of MET data requires the combination of several approaches including data manipulation and modeling to obtain the adjusted means. Mixed linear model has been more frequently used to analyze MET Fig. 2 Significant SNPs associated with fall dormancy located in the coding region of nuclear fusion defective 4-like (NFD4) gene. a Manhattan plot of chromosome 5 in response to FD in WA. b Chromosome location and gene structure of NFD4. Blue boxes represent exons. c Boxplots of FD scaled scores corresponding allele types at the NFD4 locus in Idaho (ID) and Washington (WA). Allele types are coded as 0, AAAA; 1, AAAB; 2, AABB and 3, ABBB data. The adjusted means can also be obtained by general linear models (GLM). However, GLM focuses on modeling the main response, but is unable to identify genotype stability to environmental changes. Some examples of GLM for MET data include regression on the mean model (Finlay and Wilkinson 1963) or additive main effects and multiplicative interactions (AMMI) model (Crossa et al. 1990). MLM has become popular for analyzing MET data because it can incorporate random factors and use incomplete data from different designs such as augment and incomplete block design. MLM can model GEI in terms of heterogeneity of variances and covariances (Isik et al. 2017). In the present study, the MET data were analyzed in single-trial and stagewise using MLM. Eight mixed linear models were used to estimate the adjusted means. The use of adjusted means such as BLUEs and BLUPs allows inclusion of spatial and environmental variations in multi-environment trials, in which, environment was defined as the interaction of location by phenotyping time (Smith et al. 2001).
In this work, each phenotypic dataset was first analyzed using a single-trial approach (ST0). ST0 corresponds to the first analyses of individual trials to obtain the adjusted means values of accessions in each environment. During modeling, accessions were assumed as fixed effects, blocks as random effects and row and column positions as residuals connecting spatial analysis with field design. For single-trial analysis, autoregressive (AR1) and symmetric autoregressive (SAR) models provided the best performance in most datasets tested based on AIC values. Those models are known as spatial correlation models, and they are specifically designed to model spatially autocorrelated data based on neighborhood relationships (Wolfinger 1996;Gilmour et al. 1997). Spatial models are especially useful in unbalanced or nonreplicated experiments such as ARCBD. Models AR1 and SAR follow the assumption that phenotypic data from different entries are uncorrelated, but covariance structure holds for spatial distribution. Previous reports in data modeling of spring barley and winter wheat showed that the use of spatial models in augmented designs outperformed the basic model without spatial information, and the model AR1 showed the best fit in most of the cases (Müller et al. 2010).
Stagewise analysis used BLUE values of the single-trial analysis, weights calculated as reciprocal of the variances of the means, and variance-covariance matrix from the individual environment (Piepho et al. 2012). Stagewise adjusted means (BLUPs) were estimated as accession by location by harvest (ST1), accession by location (ST2), or accession (ST3). Five stagewise models were compared (DIAG, CORGH, FA1, FA2, and FA3) and FA2 and FA3 produced the lowest AIC values. FA-k models incorporated spatial variation for individual trials and perform extremely well in terms of priding good fit to the data in comparison with other MET models (Kelly et al. 2007). Our finding are consistent with previous results where FA-k models improve selection accuracy up to 6% (Smith et al. 2001(Smith et al. , 2015. There is no single model and data analysis fitting all experiments even with the same experimental design. For that reason, is fundamental to identify appropriate statistical models to obtain correct analyses of data. Our approach was to improve data modeling to obtain BLUPs to be used in GWAS and increase the accuracy to identify associated markers. By comparing GWAS results across environments, we identified 15 common markers in three or more environments. Marker Chr6_93711284 associated with yield was identified in seven environments in WA. In a similar way, seven markers (Chr5_16888158, 160, 161, 163, 169, 170, and 188) associated with FD were identified in three environments in WA. We expected that most markers were associated with a trait across most environments as previous studies in alfalfa (Yu 2017;Lin et al. 2021); however, it was not observed in our experiments. Possible explanations include the presence of multiple loci controlling yield-related traits throughout the year, with allele-by-cut interactions (Robins et al. 2007), or the effect of environmental variation during each harvest affecting the phenotypic response and therefore affecting the GWAS results. Adjusted means of multiple harvest by location (ST2) in five agronomic traits and adjusted means of the overall response across locations (ST3) in PH identified 25 and one significant marker, respectively. Those 26 markers are important because they were associated with adjusted means values from different environments, suggesting robust markers.

Nuclear fusion defective 4-like protein is associated with fall dormancy in alfalfa
Fall dormancy is a measure of fall growth of alfalfa. This trait is fundamental for the survival of perennial plants via winter hardiness. There is a positive relationship between winter injury scores and shoot height in autumn (Weishaar et al. 2005); therefore, fall-dormant genotypes typically pay for cold-hardiness with a penalty of a low spring-summer yield (Stout and Hall 1989). It has been reported that environmental factors, such as photoperiod, and soil moisture, affect FD (Shih et al. 1967). The mechanism of genes regulation of FD via photoperiod and temperature in alfalfa remains unknown. Several transcriptomic analyses of alfalfa identified differentially expressed genes under cold acclimation related with transduction of light and photoperiod signals, receptor kinases, accumulation of raffinose-family oligo-saccharides, and biosynthesis of amino acids as proline (Zhang et al. 2015;Du et al. 2017;Liu et al. 2019) However, these studies only show the transcription changes in two varieties: Maverik (FD 1) and CUF101 (FD 9) limiting the molecular mechanisms related with FD. Our diverse panel contains 60 varieties with reported FD which allowed to identify new markers associated with FD (Table S3).
Marker Chr6_106457549 located in a locus annotated as peptidyl-prolyl cis-trans isomerase with roles in plant growth and development, hormone signaling, and stress response (Singh et al. 2020); marker Chr8_34732089 at a locus annotated as prohibitin with roles in the production of reactive oxygen species and induced premature leaf senescence (Ahn et al. 2006); and marker Chr4_11411532 at a locus annotated as glutathione s-transferase-related protein induced by abiotic stress that protects the cell from oxidative damage (Kumar and Trivedi 2018). Additionally, we identified ten markers (Chr5_16888158, 160, 161, 163, 169,170, 172, 188, 190 and 195) located at the same locus in exon 2 of  (Figs. 2c, 3b). NFD proteins have been reported to have roles in karyogamy specifically in the fusion of the polar nuclei during female gametophyte development and response to salt stress (Portereiko et al. 2006); NFD4 has also been predicted to be a nodulin protein located in mitochondria (Emanuelsson et al. 2000;Portereiko et al. 2006). Markers associated with FD in the NFD4 gene were consistently identified multiple times in ID and WA (Table S9). Additionally, allele type had great repercussion in phenotypic responses. Allele type 1 and 3 increased the FD scores at SNPs Chr5_16888158, Chr5_16888160, Chr5_16888161 and Chr5_16888169. While allele type 2 had no affect on FD. Allele type 3 increased the FD scores at SNPs Chr5_16888163 and Chr5_16888170 (Fig. 2c). Despite clear association of SNPs in this gene, the specific role of NFD4 on fall dormancy remains unknown.

Comment and specific markers associated with yield and plant height
Despite different requirements for the varieties and multiple traits tested during alfalfa breeding, there are two traits that are the most important for alfalfa growers: high biomass yield and good forage quality. However, the improvement in alfalfa biomass yield has been hampered since 1990 (USDA 2022) for multiple reasons, including longer selection cycles, low breeding investment, outbreeding mating system with inbreeding depression, tetrasomic inheritance, or the fruitless to measure yield as harvest index since the entire aboveground portion of the alfalfa plant is harvested (Annicchiarico et al. 2015). Alfalfa breeders have been more successful in improving pest resistance, winter tolerance, and forage nutritive value, but be less successful in yield (Li and Brummer 2012). The lack of knowledge of the genetic basis of yield-related traits hinders the genetic gain of alfalfa production. Agronomic traits such as biomass yield, plant height, maturity and fall dormancy are important traits for yield. Plant height and biomass yield are quantitative traits that exhibit continuous variation because they are controlled by many genes and interact with the environment (Annicchiarico et al. 2015). In the present study, we observed considerable variation in yield across locations and harvests, implying that environmental factors play important roles in yield responses. The strong positive correlations between PH and Yi in all environments suggested that the two traits may be controlled by common genetic factors.
In this analysis, markers Chr8_63386102 and Chr8_63386124 associated with PH and Yi were located at the same locus annotated as Tubby-like F-box protein, and marker Chr8_63400012 associated with PH was located at 13 kb in the adjacent gene annotated as PLP3B coding for a thioredoxin domain-containing protein. Tubby-like F-box proteins are bipartite transcription regulators involved in multiple processes (Lai et al. 2004). Tubby-like proteins have been reported as potential candidates in the dry matter content of alfalfa (Sakiroglu and Brummer 2017).
PLP3B is involved in activity regulation of target proteins through changes in the redox state of thiol groups. In Brassica rapa L. ssp. chinensis, PLP3B was downregulated in response to oxidative stress by cadmium exposure (Yu et al. 2017b).
Additional markers were associated with only one trait. Marker Chr4_79237342 associated with PH was located at the locus annotated as glutamate-rich WD repeat-containing protein 1 (GRWD1). GRWD1 is a transcription regulator controlling ribosome biogenesis and histone methylation. Marker Chr3_5642102 was significantly associated with yield in OR during the 2020 first cut. The marker is at the locus annotated as FAR-RED ELONGATED HYPOCO-TYL 3 (FHY3) protein, which is a positive regulator of phytochrome A. The fhy3 mutant displays an elongated hypocotyl elongation and anthocyanin accumulation (Wang and Deng 2002). FHY3 contains three domain structures with an N-terminal C2H2 zinc finger domain, a central putative core transposase domain, and a C-terminal SWIM motif. The N-terminal domain is essential for direct DNA binding, and the central and C-terminal domain is essential for the transcriptional activation of FHY3 (Lin et al. 2008). Marker Chr5_19865639 associated with annual forage yield in WA 2019 is at a locus annotated as protein TORNADO 2 (TRN2). TRN2 is a member of the tetraspanin family with role in leaf patterning: lamina venation, symmetry, and lateral growth (Cnops et al. 2006) and mutations in TRN2 affects mainly the activity of the shoot apical meristem which is fundamental for development and growth of the aerial plant body. They may be genetically controlled by common loci as well as specific loci with different mechanisms.

Markers associated with dry matter content and maturity stage
Dry matter content is everything in the sample other than water including protein, fiber, fat, minerals, etc. High DM is desirable because moisture affects the weight of the forage but does not provide nutrient value to the animal. Additionally, alfalfa dry matter content reveals the actual amounts of various nutrients available in the forage (Putnam et al. 2007). Thus, dry matter content has become the point of reference to optimize the forage and animal performance. Previous works identified three SNPs associated with DM in alfalfa. SNPs were in genes annotated as Tubby-like proteins, E3 SUMO-protein ligase SIZ1 and RNA-dependent RNA polymerase family protein (Sakiroglu and Brummer 2017). In this work, we identified five significant markers associated with DM in alfalfa, but only two markers were in coding regions. Marker Chr4_62801591 was in the coding region of the glutelin precursor accumulation4 (GPA4) gene that encodes the protein Golgi transport 1B (GOT1B). GOT1B is a protein known to be involved in vesicle transport from the endoplasmic reticulum (ER) to the Golgi apparatus, and its misfunction affected the anterograde transport of secretory proteins. Mutants gpa4 produce a retention of proguletins in the ER lumen in endosperm cells in mature seeds in rice forming proteins bodies within the ER lumen (Wang et al. 2016). In rice, the accumulation of proguletins changed the seed phenotype producing floury endosperm, dry seeds and less dry weight of grains (Ren et al. 2014). In alfalfa, glutelin proteins also are predicted to be localized in the cell wall of most lignified tissues (Sergeant et al. 2019) and SNPs in GPA4 could affect the lignification, callose deposition and finally DM. Marker Chr4_42520012 was located in the coding region of the RGA2 gene, which encodes a CC-NBS-LRR avirulence protein RGA2, with a reported role in resistance to leaf rust caused by Puccinia triticina in wheat (Loutre et al. 2009) or late blight Phytophthora infestans in potato (Song et al. 2003). However, NBS-LRR genes have additional roles in abiotic stress (Chini et al. 2004). For example, methylation of RGA2 gene causes downregulation in root of waterlogging-sensitive wheat preventing the active transport of substances across a membrane (Pan et al. 2020). In alfalfa, NBS-LRR genes have been associated with biomass yield under drought stress (Yu 2017) and digestible neutral detergent fiber at 30 and 48 h (Lin et al. 2020). This new report suggests that NBS-LRR proteins could be involved in multiple responses in alfalfa.
Alfalfa maturity stage was measured as mean stage by count according to Fick and Mueller (Fick and Mueller 1989). Alfalfa maturity is a key factor that controls alfalfa quality and developmental stage, influencing the time of harvest. There is positive correlation between MSC and harvest interval, amylase-treated neutral detergent fiber, yield and acid detergent lignin and inverse correlated with crude protein and neutral detergent fiber digestibility . Therefore, high values of MSC decrease alfalfa digestibility. GWAS of alfalfa maturity stage identified 18 significant markers linked to genes with roles in photosynthesis, DNA binding, mRNA binding, and cell wall organization. Marker Chr2_63316866 was in the gene phosphoenolpyruvate carboxylase (PEPC), which plays a key role during C4 photosynthesis, pH regulation, and stomatal opening (Cousins et al. 2007). In alfalfa PEPC plays a key role in CO 2 fixation and nitrogen assimilation in root nodules (Pathirana et al. 1992) through malate formation which involves in metabolization of fatty acids and amino acids as asparagine, the major form of N allocated in shoots (Fischinger and Schulze 2010). In plants, C3H1 activity of PEPC is associated with tolerance to different abiotic stresses such as high temperature (Bandyopadhyay et al. 2007), drought (Gu et al. 2013), or salt stress (Azimi et al. 2012). Marker Chr1_53838421 was in a region annotated as C3H1-type domain-containing protein. C3H1 is a RNA binding protein involved in mRNA decay (Adachi et al. 2014). In plants C3H1 is a seed specific gene controlling seed development and plant growth. Overexpression of C3H1 decrease the root development and the seed yield in tobacco (Ji et al. 2022). Marker Chr5_73969365 was in a gene coding for RNA binding protein (RBP) with GIP1_N and DUF1296 domains. RBP acts as a transcription factor involved in plant development. GIP1 may act as a coactivator that regulates transcription factors involved in lateral organ development of plants, such as bZIP transcription factors and LBD18 (Lee et al. 2014). Alfalfa breeders are interested in low values of MSC because they are related with other quality traits such as low lignin levels, lower amylase-treated neutral detergent fiber or high crude protein (Sulc et al. 2019). To the best of our knowledge, there are no other reports of molecular markers associated with MSC, and the associated TFs reported in this work may play an important role during alfalfa development, affecting alfalfa maturity and alfalfa quality during the early stages of development.

Linkage disequilibrium and QTL effects
In this work we analyzed LD decay in the diversity panel of germplasm and found that the average LD decay was 23.3 Kbp at r 2 = 0.1 . This result was slightly lower than that reported by Li et al. (2014), where the LD decay was 26 Kbp at r 2 = 0.2 . The rapid LD decay might be caused by the heterozygosity of autotetraploid alfalfa in the diversity panel of germplasm used in this experiment. The high percentage of heterozygous markers (94.5%) with more than half as diallelic duplex SNPs (AABB) was observed ( Figure S1b). LD blocks < 14 Kbp were found in markers associated with five agronomic traits. The pairwise analysis between markers in the same regions on chromosomes 6 and 8 showed r 2 values close to 1, indicating tight linkage between them. However, LD was not observed between markers in different chromosome regions. Larger LD block was found at chromosome 8 between the position 63,386,102 to 63,400,012 where three markers associated with PH and Yi. LD block was not detected for the rest 44 markers using a window of 10 kb, suggesting that these loci are likely to contribute to phenotypic effects independently. The effect values of the detected QTLs were calculated using the partial R 2 for each significant marker using the function of Fit. QTL of the R package GWASpoly (Table S8). Fit.QTL uses backward elimination to remove markers that are not predictive for the target variable or have the least predictive power (Rosyara et al. 2016). This approach was used to select a subset of markers fitting a multiple linear regression model with all the independent variables (markers). The marker with the highest p-value is removed from the model and fitted with a new model. The process is repeated marker by marker in the model until reached a threshold. In this study, relative high values (> 12%) of the phenotypic variance explained were obtained for markers associated with MSC (Chr2_56825927, Chr2_63316866, and Chr4_5209571), and PH (Chr6_15133674, Chr8_63400012, Chr1_11444898). One explanation is that the Beavis effect might cause overestimating QTL effects in our population of 200 individuals (Xu 2003).
GWAS uses an MLM to estimate correlations between phenotypes and genotypes while taking into account the relatedness between individuals. However, most GWAS are focused on obtaining the best DNA markers by different quality control filters, forgetting that half of the success in GWAS relies on the phenotypic data modeling. In the present study, we analyzed five agronomic traits: maturity measured as mean stage by count, fall dormancy, biomass yield, dry matter content, and plant height in a diverse panel of 200 individuals. Phenotypic values were modeled using single-trial and stagewise approaches to obtain estimated mean values with normal distributions. We conclude that stagewise is superior to single-trial analysis, and FA-k models were superior to diagonal, unstructured, and US-correlation models. FA-k models with more than one factor had the lowest AIC value when the number of environments was higher. BLUE/BLUP values used in GWAS identified 114 non-redundant SNP markers significantly associated with five agronomic traits: yield, plant height, fall dormancy, dry matter content, and maturity stage. We propose a clear workflow to identify significantly associated SNP markers by GWAS in augmented experiment design in multi-environmental trials using a stagewise approach (Fig. 5). Ten SNPs markers were significantly associated with fall dormancy at the same locus annotated as the NFD4 gene. The next step is the experimental validation of the most promising markers and to implement them in a MAS program. For example, our group has successfully discovered and validated Verticillium wilt (VW) associated markers in alfalfa. GWAS and QTL mapping allowed to identify of two NBS-LRR genes associated with resistance to VW resistance (Yu et al. 2017a(Yu et al. , 2020, and the effect of resistance genes MsVR38 and MsVR39 was validated in silico, in vitro and in planta (Lin et al. 2022). The allele types at the NFD4 locus could be useful in classifying FD scores. In this case we propose amplify, clone, and sequence the exon 2 of NFD4 to validate the haplotypes in accessions with high and low FD. Additionally, the haplotypes can be used to develop DNA markers for marker-assisted selection in alfalfa breeding programs. The associated markers and linked genes provide new insight toward molecular breeding in improving the target agronomic traits of alfalfa.
Author Contribution statement LXY and SN conceived this study. SL and CM performed the genetic and association analysis. CM, SL, and LXY wrote the manuscript. SN, DC, GW, GS, SF, and DL conducted field trials and collected the phenotypic data from each location.

Fig. 5
Flowchart of strategies to identify significantly associated markers in multi-environment trials (MET). Single-trial and stagewise approaches were applied using ASReml-R software. Single-trial model used entries as fixed effects, blocks as random effects and field position coordinates of plots as residuals to obtain the specific means as best linear unbiased estimation (BLUE) values (ST0). Control plots (Checks) and spatial distribution information were incorporated as an autoregressive (AR) model. The Stagewise model used ST0 BLUE values, weights, and variance-covariance (VCOV) matrices to obtain adjusted means as best linear unbiased predictions (BLUPs). Stagewise models include several structures for modeling genotypeby-environment interactions (GEI) including factor analytic (FA-k) form. Adjusted means were estimated by accession by location by harvest (ST1), accession by location (ST2), or by accession (ST3). ST0, ST1, ST2, and ST3 were analyzed by GWASpoly using a Q + K model. Significantly associated markers were annotated using the reference transcriptome dataset (RTD) of Medicago sativa