Genome-Wide Association Study (GWAS) Identifies Key Candidate Genes Associated with Leaf Size in Alfalfa (Medicago sativa L.)

Abstract

Leaf size significantly impacts photosynthetic capacity and forage yield in alfalfa, a major legume forage crop. Therefore, elucidating the genetic factors governing leaf development is critical for breeding improved alfalfa varieties. In this study, a genome-wide association analysis (GWAS) was performed to dissect the genetic architecture of leaf length (LL) and leaf width (LW) using 220 alfalfa accessions phenotyped over three years. Substantial variation for both traits was observed across environments, with coefficients of variation ranging from 10.09–16.53%. GWAS identified 26 significant SNPs associated with leaf morphology spread across seven chromosomes. Each SNP accounts for 9.7–15.6% of the phenotypic variance. Haplotype analyses confirmed positive correlations between the number of superior alleles and both LL and LW. BLAST searches revealed six candidate genes involved in leaf development within 20 kb flanking regions of significant SNPs. Our results provide novel marker-trait associations and candidate loci to facilitate molecular breeding efforts to optimize leaf size and improve productivity in alfalfa. This study establishes a foundation for integrating favorable alleles into future alfalfa varieties.

1. Introduction

Photosynthesis is a fundamental process that underpins plant growth, development, and organic matter accumulation [1]. As the primary photosynthetic organ, leaves play a pivotal role in carbon fixation [2]. Leaf morphology and area impact key physiological processes, including photosynthetic rate, transpiration, and carbon sequestration, which collectively determine biomass yield [3,4,5]. Modifying leaf size presents a promising route to improve productivity, particularly for forage crops such as alfalfa, where both stems and leaves are harvested. Elucidating the genetic factors governing leaf size will provide molecular markers to guide targeted breeding efforts to optimize leaf traits and improve crop performance.

Leaf morphological development is a complex process involving many functional genes and transcriptional regulators [6]. Overexpression of genes such as Auxin Regulated Gene involved in Organ Size (ARGOS) and Small Auxin-up RNA (SAUR) results in larger-than-normal leaves [7,8]. In addition, miRNAs play key regulatory roles in leaf development and morphogenesis across plant species [9]. For instance, miR319 regulates the transcription factor TEOSINTE BRANCHED1/CYCLOIDEA/PROLIFERATING CELL FACTORS (TCP), which is crucial for cell proliferation and leaf growth. Overexpression of miR319 or loss of TCP function induces enlarged and curled leaves [10]. In Arabidopsis, the REVOLUTA transcription factor involved in leaf development is regulated by miR165 [11]. Quantitative trait locus (QTL) mapping and gene cloning have also uncovered loci and candidate genes associated with leaf size variation in major grass crops, including maize [12], rice [13], wheat [14], and barley [15].

Alfalfa (Medicago sativa L.), known as the “Queen of Forages”, is a globally cultivated forage crop valued for its high biomass yield and nutritional quality [16,17]. The change in its growth and development determines its yield and economic benefit [18]. Increasing biomass is a key breeding objective for alfalfa. However, alfalfa poses challenges for genetic analysis and plant improvement as an autotetraploid species with a complex and self-incompatible genome [19].

In the post-genomic era, the ability to identify gene function remains a major challenge in molecular biology [20]. This is crucial for advancing the process of crop breeding. At present, the establishment of multi-omics databases has accelerated the research process of molecular breeding. Sun et al. [21] developed Milletdb, which is a millets multi-omics database containing a large amount of data. This database can provide effective services for functional genomics and population genetics research in millets. In addition, QTL mapping and genome-wide association studies (GWAS) are two primary strategies used to dissect the genetic architecture of complex quantitative traits, such as biomass yield [22]. QTL mapping employs biparental populations to identify genomic regions harboring loci that influence trait variation. Subsequently, significant QTL can be validated for marker-assisted selection or integrated with GWAS to refine candidate regions [23]. Compared with QTL mapping, GWAS offers a higher mapping resolution by exploiting historical recombination events across diverse germplasms. It also bypasses time-consuming population development. GWAS-derived markers linked to traits of interest can be directly implemented in genomic selection [23]. Recently, GWAS has been extensively used to understand the genetic underpinnings of leaf morphology, which is crucial for alfalfa productivity. For example, Chiteri et al. [24] used GWAS to identify four candidate genes associated with multiple leaf traits in mung beans, which are suitable candidate genes for further study of their roles in leaf development, growth, and function. GWAS has also revealed leaf-size loci in other crops, including rice [25], maize [26], poplar [27] and wheat [28]. The application of GWAS to alfalfa will likely provide new insights into leaf trait genetics to guide breeding for optimized leaf morphology and improved biomass yield.

Despite its agronomic significance, the genetic underpinnings and molecular mechanisms of leaf morphology in alfalfa are less explored than those in other crops. GWAS offers a powerful approach to elucidate the genetic basis of complex quantitative traits in alfalfa, including leaf morphology. The ubiquity, abundance, and precise genomic localization of SNP markers make them highly amenable to GWAS and subsequent applications in marker-assisted breeding [29]. In this study, we performed an SNP-based GWAS on 220 alfalfa accessions to identify novel loci that control leaf morphology. Our objectives were to identify candidate genes influencing leaf traits and provide informative markers to facilitate molecular breeding efforts aimed at optimizing leaf characteristics and plant productivity in alfalfa.

2. Materials and Methods

2.1. Plant Materials and Growth Conditions

The plant materials used in this study comprised 220 alfalfa accessions collected worldwide [30]. The phenotypes of different individuals of the same germplasm resources used by us were more consistent, and there were great differences among germplasm resources. Field experiments were conducted at the International Agricultural High-tech Industrial Park of the Chinese Academy of Agricultural Sciences in Langfang City, Hebei Province, China (39.59° N, 116.59° E). This region is characterized by a warm, temperate continental monsoon climate with an average annual temperature of 11.9 °C and precipitation of 554.9 mm. The soil was of the medium loam type with a pH of 7.37, containing 1.69% organic matter. In 2017, 220 materials were planted in a greenhouse, and in April 2018, a single plant with the same growth in each material was selected and transplanted into the field. Because alfalfa is a perennial crop, it remains in the field for many years after planting. The experimental design was a completely randomized block with three replicates per accession and five plants per plot. No supplemental fertilization or irrigation was applied. Weeding was performed manually as needed. Prior to phenotypic data collection, all plants were uniformly trimmed to homogenize growth.

2.2. Phenotypic Data Collection and Analysis

Leaf length (LL) and leaf width (LW) were measured in 2018, 2019, and 2020 at the stage when the first flower appeared for each genotype. For consistent sampling across years, the middle leaflet of the third or fourth fully expanded trifoliate leaf, counted from the shoot tip, was collected from each plant. Three leaves were measured per plant. LL was quantified at the longest part and LW at the widest part of each leaf using a ruler. The average of the three leaves represents the LL and LW for an individual plant. Both single-year and three-year average data were analyzed for LL and LW; single-year phenotype dates were denoted as 18LL, 18LW, 19LL, 19LW, 20LL, and 20LW, and three-year average phenotype dates were denoted as LL-mean and LW-mean. Statistical analysis of the phenotypic data was performed using Excel 2016, SPSS 19.0, and R 4.1.3 software. Origin 2022 software was used to perform Pearson correlation analysis and visualization of leaf size data. We used the R package “lme4” to calculate the generalized heritability (H²) and related contents. H² was calculated as follows:

$$H^2=\frac{V_g}{V_g+\frac{V_e}{L}}\times 100\%$$

where V_g denotes the genetic variance, V_e denotes the residual variance, and L denotes the number of environments.

2.3. Sequencing and SNP Calling

Based on field phenotype observations, DNA was extracted from 100 mg of fresh young leaf tissue from a representative plant (a single plant in each material that was relatively consistent with the other individuals) selected from each accession using the CWBIO Plant Genome DNA Kit (CoWin Biosciences, Beijing, China). Library construction and whole-genome resequencing were performed at Beijing Berry and Kang Biotechnology Co., Ltd. (Beijing, China) to generate ~10 Gb of raw sequencing data per sample, with ≥85% of reads exceeding the Q30 quality threshold. The raw sequencing data were deposited at the National Genomics Data Center (NGDC, https://bigd.big.ac.cn/, accessed on 8 November 2023, BioProject: PRJCA004024).

The SNP calling pipeline was described previously [30]. Briefly, resequencing reads were aligned to the Zhongmu-4 reference genome, and SNPs were identified using the criteria of minor allele frequency (MAF) ≥0.05, missing rate <10%, and minimum sequencing depth >5. This yielded a final set of 875,023 high-quality SNPs for downstream GWAS.

2.4. Genome-Wide Association Study (GWAS) and Haplotype Analysis

GWAS utilized 875,023 previously identified high-quality SNPs. GWAS was implemented in TASSEL 5 software using the general linear model (GLM) approach [31]. The R² in the GWAS results directly derived from Tassel can represent the phenotypic variance explained (PVE) by a single SNP. We then used the R package “CMplot” to visualize the GWAS results. Significance thresholds were set at a logarithm of the odds (LOD) score of ≥6. Manhattan and quantile-quantile (QQ) plots showed GWAS results across the genome and assessed expected and observed p-value distributions.

GWAS results were analyzed in combination with haplotypes. First, we exported the haplotype of the selected SNPs corresponding to each material from the genotype HapMap file and then combined it with the phenotypic data of each material to identify the haplotypes with better phenotypic characteristics as superior alleles. The significance between groups was tested using the t-test. Here, we used LL-mean and LW-mean as phenotypic data for haplotype analysis. For further analysis of the superior alleles, we selected 13 Chinese cultivars from 220 accessions to determine the proportion of superior alleles in these 13 materials. These 13 materials are CF032020, CF002722, CF020901, CF000715, CF020976, CF030056, PI502646, PI499544, PI502647, PI491400, PI491401, PI430638 and PI430636 (Table S1).

2.5. Candidate Gene Analysis

Genes within the significantly correlated loci were identified using the reference haploid genome of Zhongmu-4. All genes within 40 kb (20 kb upstream and downstream) of significant SNPs were identified according to the linkage disequilibrium (LD) of the association panel. These genes were searched using BLASTP at the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/, accessed on 8 November 2023). Genes with known roles related to plant development and leaf morphology were prioritized as candidates influencing associated leaf traits.

3. Results

3.1. Phenotypic Data Analysis

LL and LW exhibited extensive variation across the 220 alfalfa accessions in all environments (

Table 1. Analysis of single- and three-year mean phenotypic data for leaf length (LL) and leaf width (LW).

Category	Average	SD	Median	Minimum	Maximum	Skewness	Kurtosis	SE	CV (%)	p-Value
18LL	1.76	0.22	1.8	1	2.3	−0.16	0.36	0.01	12.56	<0.001
19LL	2.53	0.31	2.6	1.6	3.3	−0.49	0.5	0.02	12.1	<0.001
20LL	2.3	0.27	2.3	1.5	3	−0.14	0.22	0.02	11.6	<0.01
LL-mean	2.17	0.22	2.19	1.09	2.59	−0.29	−0.43	0.01	10.09	<0.001
18LW	0.83	0.1	0.8	0.5	1.1	−0.14	0.02	0.01	12.08	<0.001
19LW	1.44	0.24	1.5	0.6	2	−0.4	0.46	0.02	16.53	<0.001
20LW	1.09	0.14	1.1	0.7	1.6	0.32	1.08	0.01	13.22	<0.001
LW-mean	1.1	0.15	1.12	0.49	1.46	−0.71	0.8	0.01	13.25	<0.001

). The coefficient of variation ranged from 10.09% (LL-mean) to 16.53% (19LW), indicating substantial phenotypic diversity. In addition to 20LW, the skewness and kurtosis values fell between −1 and 1 for all traits (Table 1). The normal distribution plots further confirmed that LL and LW followed typical quantitative genetic patterns (Figure 1). Pearson correlation analysis revealed strong positive correlations between LL and LW within and across years (Figure 2). The correlation between the LL-mean and LW-mean was 0.78, implying that a shared genetic basis may underlie the two traits. Broad-sense heritability estimates were moderate for LL (65.09%) and LW (56.62%), suggesting that although genetic factors play an important role, environmental factors also substantially influence the observed phenotypic variation.

3.2. Genome-Wide Association Studies

GWAS was performed using 875,023 high-quality SNPs distributed across the eight chromosomes of the alfalfa genome. Manhattan plots were generated, depicting SNP p-values (negative log-transformed) against their genomic positions for each trait (Figure 3). QQ plots were used to compare the observed and expected p-value distributions (Figure S1). Using a significance threshold of LOD score ≥ 6, a total of 26 significant SNPs were identified on seven chromosomes associated with LL and LW (

Table 2. Significant SNP markers related to leaf size.

Trait	Marker	Variant	Superior Allele	Chromosome	p-Value	LOD	R2 (%)
18LL	chr4__21547816	T/C	T/T	4	7.88 × 10−8	7.103	10.9
	chr4__25161830	A/G	A/A	4	8.01 × 10−7	6.096	10.7
	chr1__32276127	G/A	G/G	1	9.97 × 10−7	6.001	10.5
18LW	chr1__62417155	C/T	C/T	1	1.34 × 10−7	6.873	11.6
	chr8__32102377	A/G	A/G	8	3.68 × 10−7	6.435	12.4
	chr4__45272900	C/T	C/T, C/C	4	6.31 × 10−7	6.2	12
	chr4__2117557	A/G	A/G	4	6.69 × 10−7	6.175	10.4
	chr6__104023158	G/A	A/G	6	7.77 × 10−7	6.11	10.3
19LL	chr4__10421186	A/T	A/T, A/A	4	1.27 × 10−7	6.897	13.5
	chr5__24099977	T/G	T/T	5	2.09 × 10−7	6.68	11.4
	chr8__28919244	T/G	T/T, G/T	8	8.55 × 10−7	6.068	12.2
19LW	chr4__10421186	A/T	A/A, A/T	4	1.69 × 10−9	8.773	14.7
	chr5__24099977	T/G	T/T	5	4.12 × 10−8	7.385	11.1
	chr3__54671421	A/G	A/G, A/A	3	1.60 × 10−7	6.795	11.5
	chr5__33941257	C/G	C/G	5	8.43 × 10−7	6.074	11
	chr8__16727552	T/G	G/T, T/T	8	9.89 × 10−7	6.005	10.4
20LL	chr5__56222287	C/A	C/C	5	1.44 × 10−7	6.841	11.2
20LW	chr7__25704991	T/C	T/T	7	9.88 × 10−7	6.005	9.8
LL-mean	chr8__34096759	T/C	C/T	8	7.76 × 10−7	6.11	10.1
	chr3__49868041	C/T	C/C	3	9.14 × 10−7	6.039	11.5
LW-mean	chr6__22371428	G/A	A/G	6	2.48 × 10−10	9.605	15.6
	chr5__46799336	C/T	C/C	5	2.49 × 10−8	7.603	12.4
	chr1__12214599	G/T	G/T, G/G	1	1.64 × 10−7	6.785	12.5
	chr4__23137976	T/C	T/T	4	4.26 × 10−7	6.371	11.9
	chr3__18740637	G/A	A/G, A/A	3	6.65 × 10−7	6.177	11.4
	chr6__23303002	A/C	A/C	6	9.44 × 10−7	6.025	9.7

). Specifically, the significant SNPs detected across the chromosomes were as follows: three on chromosome 1, three on chromosome 3, seven on chromosome 4, five on chromosome 5, three on chromosome 6, one on chromosome 7, and four on chromosome 8. Notably, no associations were observed on chromosome 2. The proportion of PVE by the 26 significant SNPs ranged from 9.7% to 15.6%. For LL, chr4__10421186 exhibited the highest PVE (R² = 13.5%). For LW, chr6__22371428 exhibited the highest PVE (R² = 15.6%).

The GWAS results were further analyzed, and two colocalization SNPs were obtained. chr4__10421186 was colocalized in 19LL and 19LW, and chr5__24099977 was colocalized in 19LL and 19LW (Figure 3c,d). Among them, the PVE of chr4__10421186 in 19LL and 19LW was second only to chr6__22371428 in the LW-mean among all SNPs. This suggests that these markers may have a significant effect on alfalfa leaf size. To better identify SNP loci present across various phenotypes, we compiled the statistics of loci with LOD values of 4 or higher for all observed phenotypes (Figure S2). In total, 157 SNPs were detected in two or more phenotypes. The phenotypes with the most colocalized SNPs were 18LL and 18LW, with 33 colocalized SNPs. One SNP site (chr4__2185950) was co-located across the five phenotypes. This SNP site had an LOD of 4.22–5.96 across the five phenotypes.

To validate the GWAS results, haplotype analyses were performed for all 26 significant SNPs. Significant phenotypic differences were detected at each locus among the haplotypes (Figures S2 and S3). Superior alleles were defined as those that were associated with increased LL and LW values. The relationship between superior allele dosage and leaf morphology was examined by counting the number of superior alleles within each accession. Positive correlations were evident between superior allele number and both LL and LW based on the three-year phenotype averages (Figure 4a,b). Accessions with superior alleles tended to exhibit longer and broader leaves. Significant differences in LL and LW were also apparent among the groups harboring different superior allele numbers (Figure 4). Collectively, these validations support the accuracy of the GWAS outputs for further genetic dissection of leaf morphology.

To assess the potential application of the identified superior alleles for alfalfa improvement, their frequencies were examined in 13 Chinese cultivars of the 220 accessions (Figure 5). For LL, the superior allele frequencies at the six significant SNPs ranged from 50 to 92%, while the frequencies were below 40% at three SNPs (chr8__34096759, chr3__4986804, and chr4__25161830). Regarding LW, superior allele frequencies exceeded 60% at 11 SNPs, with fixation (100%) at four of these (chr4__45272900, chr3__54671421, chr8__16727552, and chr1__12214599). In contrast, the frequencies were below 40% at six SNPs, including one absent superior allele (chr5__33941257). Overall, 11 and 12 of the 13 cultivars carried superior alleles at over 50% of the LL- and LW-associated loci, respectively. However, PI430638 contained less than 50% superior alleles for both traits. These results highlight specific superior alleles at low frequencies in the current germplasm that could be targeted to optimize leaf morphology in future alfalfa breeding.

3.3. Candidate Gene Analysis for Leaf Development in Alfalfa

Genes within 20 kb regions flanking significant SNPs were identified using the Zhongmu-4 reference genome. In total, 17 significant markers were associated with 32 annotated genes in the syntenic Medicago truncatula genome (Table S2). Based on functional annotations, six candidate genes were implicated in leaf development (

Table 3. Prediction and functional annotation of candidate genes related to leaf development.

Trait	Marker	Gene Model	Position				BLAST-P
			Chromosome	Start-Pos	End-Pos	Stand	Annotation	E-Value	%ID
18LL	chr4__21547816	Msa0540020	4	21527816	21567816	-	TATA-binding protein-associated factor BTAF1 isoform X1	0	95.68%
18LW	chr1__62417155	Msa0034440	1	62397155	62437155	-	putative histone acetyltransferase	4.00 × 10−76	78.53%
	chr8__32102377	Msa1189740	8	32082377	32122377	-	GDSL esterase/lipase At5g22810	3.00 × 10−87	70.44%
	chr4__45272900	Msa0549940	4	45252900	45292900	-	ATP-dependent Clp protease proteolytic subunit 6, chloroplastic	0	95.99%
19LW	chr4__10421186	Msa0534800	4	10401186	10441186	+	protein FAR1-RELATED SEQUENCE 5-like	0	86.65%
	chr3__54671421	Msa0362890	3	54651421	54691421	-	E3 ubiquitin-protein ligase KEG	0	96.87%

). The associated genes included putative histone acetyltransferase (Chr1__62417155), E3 ubiquitin-protein ligase KEG (Chr3__54671421), TATA-binding protein-associated factor BTAF1 isoform X1 (Chr4__21547816), ATP-dependent Clp protease proteolytic subunit 6 (Chr4__45272900), protein FAR1-RELATED SEQUENCE 5-like (Chr4__10421186), and GDSL esterase/lipase At5g22810 (Chr8__32102377).

4. Discussion

Leaf size is a key trait that influences light capture and photosynthetic efficiency, making it an important target for crop yield improvement [32]. Positive correlations have been reported between increased leaf size and yield in rice, wheat, and other species [33,34]. As a vital forage crop, alfalfa is widely cultivated because of its high protein and fiber content in livestock feed [35]. Since leaves are harvested with stalks, leaf size directly impacts yield and forage quality. The genetic basis of leaf morphology in alfalfa remains poorly understood, posing challenges owing to its autopolyploid genome and high heterozygosity [19]. Advancements in next-generation sequencing have enabled high-throughput SNP discovery at a reduced cost. SNPs, insertion/deletions (indels), and structural variations (SV) as three ways of sequence variation are the basis of phenotypic variation in species [36]. It is of great importance to understand the distribution and impacts of these variants in germplasm collections in order to utilize them in breeding programs [37].

Moreover, the recent completion of the alfalfa genome sequence provides an invaluable genomic resource to facilitate GWAS and QTL mapping [19,38]. GWAS and QTL mapping are powerful approaches to elucidate the genetic architecture of complex quantitative traits in plants. For example, Chen et al. [39] performed a GWAS on tea plants and identified six candidate genes associated with leaf size and other morphological traits. The expression of the two-leaf size-related candidates was validated using RT-qPCR, which confirmed their functional roles. In maize, Miculan et al. [20] uncovered 25 candidate genes linked to leaf development through GWAS, including those involved in vacuolar function, cell wall processes, and vesicle trafficking. The GWAS outputs can provide critical markers for accelerating molecular breeding via marker-assisted selection.

We previously combined QTL mapping and RNA-seq to identify potential leaf trait candidate genes [40,41]. Here, we conducted the first GWAS for leaf morphology in alfalfa using a diverse panel evaluated over three years in Hebei Province, China. The GLM selected in this study is a model with a wide detection range that can detect many SNPs associated with target traits [39]. We identified 26 significant SNPs and six candidates influencing LL and LW. Unlike traditional bi-parental linkage mapping, GWAS exploits historical recombination within diverse germplasms without developing specialized crosses. It can also directly tag causative polymorphisms for downstream functional validation [42]. Notably, our GWAS candidates did not overlap with those of previous QTL studies, which likely reflects the population- and environment-specific associations. Therefore, the newly identified loci provide complementary insights into leaf genetic architecture to guide alfalfa improvement.

LL and LW exhibited year-to-year variation, likely attributable to environmental fluctuations, including climate and field effects [41]. In particular, the LL and LW in 2018 were significantly smaller than those in the other two years, possibly because 2018 was the first year of alfalfa transplantation. The coefficients of variation for LL and LW ranged from 10.09% to 16.53% across the environments, demonstrating substantial phenotypic diversity in these traits among the alfalfa accessions. This variation provides useful genetic potential for optimizing leaf morphology through breeding. Furthermore, we observed a high positive correlation between LL and LW, which is consistent with previous reports on alfalfa [43,44].

Haplotype analysis, an essential molecular marker, plays a crucial role in molecular breeding applications and is widely used for crops such as wheat [45], maize [46], and rice [47]. In alfalfa, a larger leaf size is considered a favorable phenotype [40]. Our haplotype analysis of the 26 significant GWAS SNPs revealed that T/T was the most frequent superior allele, occurring in 20.59% of the SNPs. This was followed by A/G and A/A, accounting for 17.65% and 14.71%, respectively. This information can guide haplotype-based selection in alfalfa breeding programs. We found that the Chinese cultivars of the panel carried superior alleles at only 22.22–77.78% of LL loci and 41.18–88.24% of LW loci, indicating potential targets for introgression via marker-assisted selection. Furthermore, superior allele dosage was positively associated with leaf size, supporting the feasibility of pyramiding favorable alleles to increase leaf area.

This study identified 26 SNPs that were significantly associated with alfalfa leaf development. Six candidate genes were identified and annotated. Chr4__45272900, located on chromosome 4, is linked to ATP-dependent Clp protease proteolytic subunit 6, which is chloroplastic, indicating its association with or localization in chloroplasts. This gene encodes an essential housekeeping enzyme in plant chloroplasts that is involved in chloroplast development and function [48,49]. Chr3__54671421, located on chromosome 3, is linked to the E3 ubiquitin-protein ligase KEG. This gene plays a key role in plant development and has various effects on leaf development [50], such as regulation of chloroplast development [51], leaf senescence [52], and stomatal development [53]. Chr1__62417155, located on chromosome 1, is linked to a putative histone acetyltransferase. Jenna et al. [54] found that the histone acetyltransferase GCN5 has specific roles in the leaf tissue of Arabidopsis, influencing cell growth and division in rosette leaves in complex and sometimes opposite ways. Chr4__10421186, located on chromosome 4, is linked to the FAR1-RELATED SEQUENCE 5-like protein. FAR1 has been reported to be a positive regulator of chlorophyll biosynthesis in Arabidopsis thaliana [55] and is related to chlorophyll synthesis [56]. Furthermore, chr8__32102377 on chromosome 8 and chr4__21547816 on chromosome 4 were linked to genes annotated as GDSL esterase/lipase At5g22810 and TATA-binding protein-associated factor BTAF1 isoform X1, respectively. Although these genes have not been directly linked to leaf development, they are known to be involved in developing plant meristems [57,58].

In summary, we conducted a comprehensive GWAS of 220 alfalfa variants grown in Langfang, Hebei Province, between 2018 and 2020. Our analysis identified 26 significant SNPs with six candidate genes related to leaf development. This study offers valuable resources for alfalfa leaf development and breeding.