The alfalfa genome used in this study was obtained from the Alfalfa Genome Project (https://figshare.com/projects/whole_genome_sequencing_and_assembly_of_Medicago_sativa/66380) (Chen et al., 2020b). This assembly represents a phased genome comprising four haplotypes across 32 chromosomes, which facilitates the resolution of duplicated and homeologous loci in autotetraploid alfalfa. The genome sequence of Arabidopsis thaliana was retrieved from the TAIR database (https://www.arabidopsis.org/), while the genomes of other plant species, including Medicago truncatula, Glycine max, and Oryza sativa, were downloaded from Ensembl Plants (https://plants.ensembl.org/index.html). The DCL, AGO, and RDR gene family members in A. thaliana were identified based on previously published studies (Supplementary Table S1) (Podder et al., 2023). To identify putative DCL, AGO, and RDR genes in the alfalfa genome, we performed BLASTP searches using known A. thaliana protein sequences as queries, with an E-value threshold of ≤ 1×10⁻¹⁰. The resulting candidate genes were further validated using the Conserved Domain Database (CDD, https://www.ncbi.nlm.nih.gov/cdd/) to confirm the presence of characteristic domains specific to each gene family: DEAD, Helicase-C, Dicer-dimer, PAZ, and DSRM domains for DCL proteins; ArgoN, PAZ, MID, and PIWI domains for AGO proteins; and RdRP domain for RDR proteins. Any sequences lacking these essential conserved domains were excluded from further analysis.

Using TBtools v2.156 (Chen et al., 2020a), we extracted information on protein sequence length, molecular weight, and theoretical isoelectric point (pI) for the identified MsDCL, MsAGO, and MsRDR proteins from the GFF annotation file of the alfalfa reference genome. The subcellular localization of these proteins was predicted using the online tool WoLF PSORT (https://wolfpsort.hgc.jp/).

Multiple sequence alignments of DCL, AGO, and RDR protein sequences from Medicago sativa, A. thaliana, Glycine max, and Oryza sativa were performed using MEGA11. The gene IDs for DCL, AGO, and RDR family members in G. max and O. sativa were obtained from previously published studies (Kapoor et al., 2008; Liu et al., 2014). A phylogenetic tree was constructed using the neighbor-joining (NJ) method with 1,000 bootstrap replicates to assess tree reliability. The resulting tree was visualized using the online tool Evolgenius (https://evolgenius.info//evolview-v2/#login). To identify conserved motifs in MsDCL, MsAGO, and MsRDR proteins, the MEME Suite (https://meme-suite.org/meme/index.html) was used with default parameters, except that the maximum number of motifs was set to 10. Gene structure information for MsDCL, MsAGO, and MsRDR was extracted from the GFF annotation file of the alfalfa genome. Visualization of the results was performed using TBtools.

To investigate the collinearity relationships of DCL, AGO, and RDR genes both within alfalfa and across different species, the MCScanX tool was employed (Wang et al., 2012). Chromosomal localization of MsDCL, MsAGO, and MsRDR family members was visualized using TBtools. Tandem duplication events were defined as two or more genes with sequence similarity exceeding 70% located within a 200 kb region on the same chromosome, a commonly used operational criterion in plant gene-family studies for capturing tandem arrays that may span larger intergenic regions. Pairwise sequence similarity between genes was assessed using MEGA11. To identify tandemly duplicated MsDCL, MsAGO, and MsRDR genes, their chromosomal positions were analyzed and compared accordingly.

TBtools was used to extract the 2,000 bp upstream promoter regions of the MsDCL, MsAGO, and MsRDR genes. Cis-acting regulatory elements within these sequences were identified using the PlantCARE online database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). After identification, the distribution and functional classification of these elements were visualized using TBtools.

RNA-seq datasets from six alfalfa tissues (roots, elongated stems, pre-elongated stems, leaves, flowers, and nodules; project ID: SRP055547) and from plants subjected to abiotic stress (salt and drought treatments; run accessions: SRR7160313–SRR7160357) were obtained from the NCBI Sequence Read Archive (SRA) (O’Rourke et al., 2015; Dong et al., 2021). Cleaned reads were aligned to the Xinjiang Daye reference genome using TopHat2. Gene expression levels were quantified using FPKM (Fragments Per Kilobase of transcript per Million mapped reads). Differential expression analysis under stress conditions was performed with DESeq2, applying thresholds of padj < 0.05 and |log₂FC| ≥ 1 to identify significant changes. The expression profiles of MsDCL, MsAGO, and MsRDR family members were subsequently analyzed and visualized using TBtools.

Due to limited material availability, the alfalfa cultivar Zhongmu No. 4 (Medicago sativa L.) was used for RT-qPCR assays in this study, and seeds of this cultivar were provided by the Institute of Animal Science, Chinese Academy of Agricultural Sciences. Prior to cultivation, the seeds were stratified at 4 °C for three days. Seedlings were then grown in a greenhouse under controlled conditions for two weeks, with a 16 h light/8 h dark photoperiod, 70–80% relative humidity, and day/night temperatures of 24 °C/20 °C. Stress treatment was conducted using a hydroponic method. Two-week-old seedlings were transferred to containers containing a 1/2 strength Hoagland nutrient solution for a 24-hour acclimation period. Subsequently, the nutrient solution was completely replaced with a 1/2 strength Hoagland solution containing 15% PEG 6000 (for drought stress) or 200 mM NaCl (for salt stress) to initiate the stress treatment. The control group continued to be maintained in a standard 1/2 strength Hoagland solution. The volume of nutrient solution in each container was kept consistent across all containers and was not replenished throughout the entire 24-hour treatment period to maintain a constant stress intensity. For both treatments, leaf samples were collected at 0, 6, 12, and 24 hours post-treatment, with the 0-hour samples serving as controls. Each treatment included three biological replicates, with five seedlings pooled per replicate. Untreated control plants were grown under the same greenhouse conditions. Total RNA was extracted from all samples using TRIzol reagent (Invitrogen, USA), following the manufacturer’s instructions. First-strand cDNA was synthesized using the EasyScript First-Strand cDNA Synthesis Kit (TransGen Biotech, China). Gene-specific primers were designed using Primer 5.0 (Supplementary Table S2). RT-qPCR was performed using SYBR Premix Ex Taq (Takara, Japan) on a 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Each sample was analyzed in technical triplicate, and gene expression levels were normalized against the alfalfa ACTIN gene. The expression stability of the ACTIN gene under salt and drought stress conditions was verified in our preliminary experiments, a finding that aligns with previous studies in alfalfa where ACTIN has been confirmed as a reliable internal reference gene under comparable abiotic stress treatments (Li et al., 2014). Therefore, it was selected as the internal control for normalization of RT−qPCR data in this study. Relative expression was calculated using the 2^−ΔΔCT method.

Based on conserved domain structures and previously characterized RNA silencing pathway genes in Arabidopsis thaliana, a BLAST search was performed against the genome of Medicago sativa cv. Xinjiang Daye. A total of 31 MsDCL, 82 MsAGO, and 52 MsRDR genes were identified in the alfalfa genome. These genes were designated as MsDCL1 to MsDCL31, MsAGO1 to MsAGO82, and MsRDR1 to MsRDR52, respectively, according to their chromosomal locations. The detailed features of these genes and their encoded proteins are listed in Supplementary Tables S3 and S4.

Within the MsDCL gene family, MsDCL6 encodes the shortest protein (1,140 amino acids, 129.52 kDa), whereas MsDCL20 encodes the longest (1,889 amino acids, 212.23 kDa). The theoretical isoelectric point (pI) of MsDCL proteins ranges from 5.87 (MsDCL19) to 7.99 (MsDCL24), and their instability indices range from 40.00 (MsDCL6) to 45.48 (MsDCL11). For the MsAGO gene family, MsAGO73 encodes the smallest protein (636 amino acids, 71.92 kDa), while MsAGO66 encodes the largest (1,129 amino acids, 124.64 kDa). The pI values of MsAGO proteins range from 8.43 (MsAGO82) to 9.77 (MsAGO59), and the instability index ranges from 37.70 (MsAGO31) to 51.83 (MsAGO54). Within the MsRDR family, MsRDR52 encodes the shortest protein (794 amino acids, 90.52 kDa), while MsRDR14 encodes the longest (1,396 amino acids, 159.48 kDa). The pI values of MsRDR proteins range from 6.21 (MsRDR29) to 9.11 (MsRDR12), and their instability indices range from 36.05 (MsRDR44) to 47.98 (MsRDR47). Comparative analysis revealed that MsDCL proteins are generally longer and heavier than MsAGO and MsRDR proteins, whereas MsAGO proteins tend to exhibit higher isoelectric points. Subcellular localization prediction (Supplementary Table S4) showed that most MsDCL proteins are targeted to the chloroplast (23) or nucleus (8). Among the MsAGO proteins, 66 are predicted to localize in the nucleus, with others distributed across the chloroplast (8), mitochondrion (5), cytosol (2), and peroxisome (1). Most MsRDR proteins are also predicted to be nuclear-localized (46), with a few targeted to the cytosol (4), chloroplast (1), or mitochondrion (1).

All identified RNA silencing pathway genes are distributed across the 32 assembled chromosomes of the alfalfa genome, with three genes located on unanchored scaffolds (21603, 32402, and 43330). Specifically, 30 MsDCL genes (MsDCL1–MsDCL30) are distributed across 20 chromosomes, with none detected on chromosomes chr4.1–4.4, chr5.1–5.4, or chr6.1–6.4. MsDCL31 is located on unanchored scaffold 32402. Similarly, 81 MsAGO genes (MsAGO1–MsAGO81) are located on 27 chromosomes, while MsAGO82 is present on scaffold 21603, with no MsAGO genes found on chr7.1–7.4 or chr8.3. The MsRDR family consists of 51 genes (MsRDR1–MsRDR51) spread across 14 chromosomes, and MsRDR52 is located on scaffold 43330. Notably, MsRDR genes are absent from chr2.1–2.4, chr4.1–4.4, chr5.1–5.5, chr6.2, chr7.1–7.4, and chr8.4 (Figure 1).

Phylogenetic trees of the DCL, AGO, and RDR gene families were constructed based on protein sequences from M. sativa, A. thaliana, G. max, and O. sativa (Figure 2). The DCL genes were grouped into four distinct clades (I–IV) according to their evolutionary relationships. Each of the four A. thaliana DCL genes was assigned to a separate clade. Among the M. sativa genes, clades I–IV contained 4, 19, 4, and 4 MsDCL genes, respectively. Similarly, the AGO genes were divided into three major clades (I–III). Consistent with previous classifications (Vaucheret, 2008; Zhang et al., 2015), AtAGO1/5/10 were grouped into clade I, AtAGO2/3/7 into clade II, and AtAGO4/6/8/9 into clade III. The M. sativa genes were distributed as follows: 33 MsAGO genes in clade I, 11 in clade II, and 38 in clade III. The RDR gene family was classified into four clades (I–IV). Clade III included the three A. thaliana genes AtRDR3, AtRDR4, and AtRDR5, while clades I, II, and IV each contained a single A. thaliana RDR gene. In M. sativa, clades I–IV contained 40, 4, 4, and 4 MsRDR genes, respectively.

To investigate the structural characteristics of RNA silencing pathway genes in alfalfa, conserved motifs were identified using the MEME online tool. A total of 10 conserved motifs were detected for each gene family, as shown in Supplementary Figure S1. The types and numbers of motifs varied across individual gene members. For the MsDCL proteins (Supplementary Figure S2), most members contained all 10 motifs. Notably, MsDCL6 lacked motif 2; MsDCL15–23, MsDCL25, and MsDCL28 lacked motif 10; MsDCL31 lacked motifs 3, 6, 7, 8, and 10. All other MsDCL proteins contained the full set of motifs. Regarding MsAGO proteins (Supplementary Figure S3), 76 proteins contained all 10 motifs. However, MsAGO48 and MsAGO82 lacked motif 2; MsAGO73 lacked motif 5; MsAGO21 lacked motifs 1 and 2; MsAGO59 lacked motifs 1, 2, 4, and 8; and MsAGO19 retained only motifs 5, 6, 7, and 10. Among the MsRDR proteins (Supplementary Figure S2), 40 members possessed all 10 motifs. MsRDR1, MsRDR4, MsRDR23, MsRDR33, MsRDR48, and MsRDR51 lacked motif 8. MsRDR46, MsRDR47, and MsRDR52 were missing motifs 3, 4, 6, 7, and 9, while MsRDR45 lacked motifs 2, 3, 4, 6, 7, and 9.

In addition to conserved motifs, protein domain architectures were analyzed. All MsDCL proteins possessed five core domains: Ribonuclease_3, Helicase_C, Dicer_dimer, PAZ, and DSRM. However, the DEAD domain was absent in MsDCL19–22, which instead contained the RESIII domain, a feature unique to these proteins (Supplementary Figure S2). All MsAGO proteins contained the Piwi domain (Piwi-like in MsAGO19 and MsAGO59), PAZ, ArgoL1, and ArgoL2 domains. The ArgoN domain was missing in MsAGO34 and MsAGO73. The Rib_recp_KP_reg domain was exclusively found in MsAGO33 and MsAGO37, while the Gly-rich_Ago1 domain was unique to MsAGO66, MsAGO67, MsAGO69, MsAGO70, MsAGO72, MsAGO75, and MsAGO76. The Herpes_TAF50 domain was present only in MsAGO31 (Supplementary Figure S3). All MsRDR proteins contained the canonical RdRP domain. Notably, MsRDR14 harbored two unique domains, HLH and ZapB, while MsRDR48–51 were characterized by the exclusive presence of the RRM_SF domain (Supplementary Figure S2). Therefore, the presence of family-specific or member-specific domains suggests that MsDCL, MsAGO, and MsRDR genes may have functional specialization or diversification within the RNA silencing pathways.

Gene duplication plays a pivotal role in the evolution of gene families by facilitating the emergence of novel genes and functions. Both segmental and tandem duplications are key mechanisms driving gene family expansion. In this study, a total of 20 tandem duplication events were identified among the MsDCL, MsAGO, and MsRDR gene families. For example, tandem duplications were observed between MsRDR2–MsRDR11 on chr1.1, with sequence similarities ranging from 80.6% to 98.5%, and between MsAGO4 and MsAGO5 on chr1.3, with 100% sequence similarity. Among these tandem duplication events, 3 involved MsDCL genes, 9 involved MsAGO genes, and 8 involved MsRDR genes. Notably, seven tandem duplication events of MsRDR genes were detected across chr1.1, chr1.2, chr1.3, and chr1.4, with events 2, 4, 6, and 8 each encompassing more than seven MsRDR genes (Figure 1; Supplementary Table S5). Further details on these tandem duplication events are provided in Supplementary Table S5.

In addition to tandem duplications, a total of 164 segmental duplication events were identified across the alfalfa genome (Supplementary Table S6). Among these, 28 were associated with MsDCL genes (Supplementary Figure S4A), 99 with MsAGO genes (Figure 3), and 37 with MsRDR genes (Supplementary Figure S4B). All segmental duplications involving MsDCL and MsRDR genes occurred between homologous chromosomes. For instance, segmental duplication was detected between MsDCL2 on chr1.2 and MsDCL5 on chr1.4, as well as between MsRDR1 on chr1.1 and MsRDR12 on chr1.2 (Supplementary Figure S4). In contrast, MsAGO genes exhibited segmental duplications on both homologous and non-homologous chromosomes. For example, a duplication event was identified between MsAGO10 on chr2.1 and MsAGO28 on chr4.1 (Figure 3).

To further elucidate the evolutionary mechanisms of the DCL, AGO, and RDR gene families, we performed a comparative collinearity analysis between M. sativa and three representative plant species: A. thaliana, M. truncatula, and G. max (Figure 4). Between M. sativa and A. thaliana, a total of 31 segmentally duplicated gene pairs were identified, including 5 involving MsDCL genes, 22 involving MsAGO genes, and 4 involving MsRDR genes. In the comparison with M. truncatula, 85 segmental duplication events were detected, comprising 19 MsDCL-related pairs, 53 MsAGO-related pairs, and 13 MsRDR-related pairs. Likewise, synteny analysis with G. max revealed 132 segmentally duplicated gene pairs, including 30 related to MsDCL genes, 79 to MsAGO genes, and 23 to MsRDR genes.

Cis-acting regulatory elements in the promoter regions of MsDCL, MsAGO, and MsRDR genes were identified using the PlantCARE database. These elements were classified into three major categories: plant growth and development, hormonal and stress responses, and light responsiveness (Supplementary Table S7). Sixteen hormone- and stress-responsive cis-elements were selected for further analysis (Figure 5). Among these, the MsAGO family contained all 16 types, while the MsDCL and MsRDR families collectively harbored 14 of the 16 types (absent: AuxRE and WUN-motif). In the MsDCL family, ABRE was the most frequently identified cis-element, suggesting its potential involvement in the abscisic acid signaling pathway. In contrast, ARE was the most abundant in the promoters of MsAGO genes, whereas the promoters of MsRDR genes were enriched for the CGTCA and TGACG motifs, implying their potential responsiveness to methyl jasmonate signaling. These results, based on bioinformatics predictions, indicate that MsDCL, MsAGO, and MsRDR genes may play distinct roles in hormone-mediated stress responses, though their specific functions require further experimental validation.

To investigate the tissue-specific expression patterns of RNA silencing pathway genes in alfalfa, RNA-seq data from six distinct tissues—roots, leaves, flowers, elongated stems, pre-elongated stems, and nodules—were analyzed. Among the 165 identified RNA silencing pathway genes, 133 genes exhibited detectable expression in at least one of these tissues. The remaining 32 genes, including six MsDCL genes (MsDCL5, MsDCL6, MsDCL7, MsDCL18, MsDCL23, MsDCL31), 20 MsAGO genes (MsAGO4, MsAGO10, MsAGO13, MsAGO16, MsAGO19, MsAGO21, MsAGO23, MsAGO26, MsAGO34, MsAGO39, MsAGO44, MsAGO52, MsAGO53, MsAGO59, MsAGO65, MsAGO68, MsAGO71, MsAGO74, MsAGO78, MsAGO81), and six MsRDR genes (MsRDR7, MsRDR12, MsRDR27, MsRDR28, MsRDR44, MsRDR51), showed no detectable expression, implying that these genes may function in other tissues or in response to specific biotic or abiotic stresses. Among the 133 expressed genes, several displayed strong tissue specificity. Notably, MsAGO46, MsAGO50, and MsAGO61 were exclusively expressed in flowers; MsRDR24 was detected only in leaves; MsDCL12 and MsAGO24 were specifically expressed in nodules; and MsRDR26 was confined to pre-elongated stems. Moreover, six genes—MsDCL4, MsDCL17, MsAGO7, MsAGO8, MsRDR10, and MsRDR32—were uniquely expressed in roots (Figure 6; Supplementary Table S8). While many genes exhibited expression across multiple tissues, their expression levels varied considerably, indicating potential functional divergence and tissue-specific regulatory mechanisms.

To confirm the consistency of expression patterns with RNA-seq data for genes in the RNA silencing pathway under abiotic stress, six genes (MsDCL11, MsAGO54, MsAGO60, MsAGO73, MsAGO76, and MsRDR37) were selected based on their expression changes under salt and drought conditions and verified by RT-qPCR under corresponding treatments (Supplementary Figure S5; Figure 7). MsDCL11, MsAGO60, and MsAGO76 were significantly upregulated at 6 h and 12 h, followed by downregulation at 24 h. MsAGO73 and MsRDR37 showed a similar trend but responded more rapidly, with their expression beginning to decline as early as 12 h. In contrast, the expression of MsAGO54 decreased continuously and progressively throughout the 24 h salt treatment. Under drought stress, all six genes displayed a similar expression pattern, characterized by an initial significant upregulation followed by a gradual decline. Among them, MsAGO54 and MsAGO73 responded more rapidly, with their expression starting to decrease by 12 h, while the other genes maintained an upward trend within the first 12 h. Notably, the expression pattern of MsAGO76 and MsDCL11 was consistent under both salt and drought conditions.