XinJiangDaYe genome assembly files (https://figshare.com/projects/whole_genome_sequencing_and_assembly_of_Medicago_sativa/66380) were used to derive the genome sequence information of alfalfa. The hidden Markov model (HMM) data of the ABC1K domain (PF03109) in the PFAM database (http://pfam.xfam.org/) were acquired to identify MsABC1K gene family members. The HMMER 3.0 software was utilized to retrieve MsABC1K protein sequences. The identified MsABC1Ks were uploaded to the National Center for Biotechnology Information (NCBI) Conserved Domain Database website (https://www.ncbi.nlm.nih.gov/cdd) to confirm the presence of conserved structural domains.

Based on genome annotation files, the TBtools software was utilized to visualize the chromosomal localization of the MsABC1K gene (Chen et al., 2020). pI values, molecular weight (MW) values, cDNA sequence length, the grand average of hydropathicity (GRAVY), and protein length of MsABC1Ks were predicted using ExPaSy online tools (https://prosite.expasy.org/) (Artimo et al., 2012).

The neighbor-joining method along with 1,000 bootstrap replicates was utilized to generate a phylogenetic tree in the MEGA software (Kumar et al., 2018). Sequences of 17 AtABC1K proteins were collected from Arabidopsis genome assembly TAIR10 (http://www.arabidopsis.org/). Subsequently, rice genome assembly v6 (http://rice.plantbiology.msu.edu/) was used to acquire the protein sequences of 16 OsABC1Ks. Next, maize genome assembly 5b.60 (http://www.maizesequence.org/) was utilized to collect the protein sequences of 19 ZmABC1Ks.

Multiple Expectation Maximization for Motif Elicitation (MEME) suite (http://meme-suite.org/tools/meme) was employed to assess conserved motifs in predicted alfalfa MsABC1Ks (Bailey et al., 2009); the obtained motifs were visualized using the TBtools software. The maximum pattern number determined was set as 10, with default values for the other parameters. Using the XinJiangDaYe genome annotation file, the MsABC1K gene structure was visualized using the TBtools software (v2.083). To identify cis-elements, the 2,000-bp promoter region located upstream of the ATG start codon of the alfalfa MsABC1K family members was uploaded to the PlantCARE website (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

BUSCA (http://busca.biocomp.unibo.it/) was utilized to predict the subcellular localization of the 22 MsABC1K genes. The secondary structure of the MsABC1K protein was predicted using the online tool SOPMA (http://npsa-pbil.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa_sopma.html). The online software ExPaSy SWISS-MODEL (https://swissmodel.expasy.org/interactive) was utilized to predict the 3D structure of the 22 MsABC1K proteins.

“Zhongmu No.1” alfalfa cultivar seeds were cultivated at the Institute of Grassland Research, Chinese Academy of Agricultural Sciences. The plants were cultured in a greenhouse at 25°C (day)/22°C (night) under a light/dark photoperiod (16/8 h). Samples of roots, leaves, stems, and flowers from a 2-month-old alfalfa plant were utilized to estimate tissue-specific expression. One-month-old alfalfa plants were treated with 300 mmol/L NaCl and 20% PEG-6000 to simulate salt stress and drought stress, respectively, and leaf samples were obtained at 0 h (CK) and 3, 6, and 12 h after treatment. Three biological replicates were used for each sample. The samples were placed in liquid nitrogen after sampling, frozen at −80°C, and then stored for quantitative reverse transcription–polymerase chain reaction (qRT−PCR) experiments.

The HiPure Polysaccharides & Polyphenolics Plant RNA Kit (Genepioneer, Beijing, China) was employed to extract total RNA. Sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA). The mRNA obtained from total RNA was purified with poly-T oligo-attached magnetic beads. cDNA was acquired by reverse transcription of the enriched mRNA. The purified double-stranded cDNA fragments were subjected to end repair, modified with A base addition, ligated to Illumina sequencing adapters, and sequenced using the Illumina NovaSeq X Plus system at the Nanjing Genepioneer Technology Company (Nanjing, China). The gene expression level was calculated based on the Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) value. The TBtools software was utilized to visualize data.

The TransZol Up Plus RNA Kit (Transgene ER501, Beijing, China) was utilized to extract total RNA from different treatments. Reverse transcription was achieved using the HiScript II Q RT SuperMix for qPCR (+gDNA wiper) Kit (Vazyme R223, Nanjing, China). The Primer Premier 6 software was used to design specific primers (Singh et al., 1998). qRT-PCR amplification was conducted using the QuantStudio 5 PCR system (Thermo, Waltham, MA, USA). The MsActin gene was set as an internal control. qRT-PCR primers are shown in Supplementary Table S1 . Experiments were performed with three replicates, and the 2^−ΔΔCt method was utilized to estimate the relative gene expression levels (Zhao et al., 2021; Wu et al., 2023). SPSS version 22.0 was used to perform statistical analyses. One-way ANOVA and Duncan’s test were used to determine differences in the parameters. Significant differences in mean values were assessed using Student’s t-test. * and ** indicate significant differences at p < 0.05 and 0.01, respectively.

Twenty-two MsABC1K genes in the alfalfa genome were recognized according to the HMM profile of the ABC1 domain (PF03109) sequence alignment and conserved structural domain analysis. The validated ABC1K genes were termed MsABC1K1 to MsABC1K15 based on their phylogenetic relationship with AtABC1Ks. Allele genes located on the homologous chromosome were differentiated using lowercase letters a, b, and c. Supplementary Table S2 shows the names, IDs, and amino acid sequences of the genes.

The annotation data of the alfalfa genome showed that the 22 MsABC1K genes were distributed on 17 alfalfa chromosomes ( Figure 1 ). Chr1.1, chr1.4, chr2.3, chr3.1, chr3.2, chr3.3, chr3.4, chr5.1, chr6.4, chr7.2, chr8.1, chr8.3, and mitochondrial chromosome 10913 contained only one MsABC1K gene. Three MsABC1K genes were present on chr4.4, and two MsABC1K genes were located on chr1.2, chr5.4, and chr6.1.

Physicochemical feature analysis showed that the coding sequence (CDS) of the alfalfa ABC1K genes was 963–2,877 bp in length and encoded 320–958 amino acids. The pI values and MW of the MsABC1K members were 5.29–9.58 and 36.26–106.13 kDa, respectively. The grand average of hydropathicity (GRAVY) values of the 22 MsABC1K proteins were negative, thus suggesting that the proteins were hydrophilic ( Supplementary Table S3 ).

Among the 22 MsABC1Ks, the subcellular localization prediction showed that MsABC1K1, MsABC1K2a, MsABC1K5b, MsABC1K9, and MsABC1K13 were located in chloroplasts and that MsABC1K2b/3a/5a/6/7/8b/10b/11/12/14/15 were located in the organelle membrane and endomembrane system. MsABC1K3b/5c/8a/10a/10c and MsABC1K4 were localized in the mitochondria and nucleus, respectively ( Supplementary Table S3 ).

The 3D structures of all 22 MsABC1K proteins were predicted using ExPaSy SWISS-MODEL ( Supplementary Figure S1 ). The MsABC1K proteins were mainly constructed based on α-helices and random coils, and most of the MsABC1K proteins had more than 50% distribution of α-helices ( Supplementary Table S4 ).

To determine the evolutionary relationship of MsABC1K proteins with ABC1K proteins of other species, the complete protein sequences of 22 MsABC1Ks, 17 AtABC1Ks, 16 OsABC1Ks, and 19 ZmABC1Ks were utilized to construct an unrooted phylogenetic tree ( Supplementary Table S5 ). The 22 MsABC1K members were classified into four groups ( Figure 2 ). Among these groups, subfamily I had the most members, including MsABC1K1/2a/2b/3a/3b/4/5a/5b/5c/6/7/8a/8b. Subfamily II had four members, including MsABC1K9/10a/10b/10c. Subfamily III had three members, namely, MsABC1K11/12/13. Subfamily IV had two members, MsABC1K14 and MsABC1K15.

To discover the conserved motifs of MsABC1Ks, the MEME online tool was utilized to analyze the motif distribution of all 22 MsABC1Ks, and 10 conserved motifs were detected ( Figure 3 ). Motifs 1, 3, 4, 5, 7, and 9 exhibited high conservation across all MsABC1Ks. Motif 2 was not detected in MsABC1K12 or MsABC1K14. Motif 6 was not found in the subfamily I member MsABC1K5a; the subfamily III members MsABC1K11, MsABC1K12, and MsABC1K13; or the subfamily IV member MsABC1K15. Motif 8 was not found in the subfamily I member MsABC1K5a. Motif 10 was not detected in the subfamily I members MsABC1K5a, MsABC1K7, MsABC1K8a, or MsABC1K8b.

To determine the MsABC1K gene structural composition, we analyzed the intron–exon structure of MsABC1Ks. MsABC1K gene structure showed differences as well as similarities among the four subfamilies ( Figure 4 ). Most of the genes (18 of 22, 81.82%) had seven exons or more, and only three genes, namely, MsABC1K10a, MsABC1K10b, and MsABC1K10c, had four exons. In subfamily I, the number of exons was 7–22, which widely fluctuated as compared to that in the other subfamilies. In subfamily II, MsABC1K9 had 10 exons, which varied from that noted in MsABC1K10a/MsABC1K10b/MsABC1K10c. Subfamily III members had 9–18 exons, while subfamily IV members had 11 or 16 exons.

By employing a 2.0-kb promoter sequence located upstream of the MsABC1K gene start codon, we discovered the cis-elements are potentially involved in the transcriptional regulation of MsABC1K. We speculated that most cis-regulatory elements showed an association with light responsiveness in the alfalfa MsABC1K gene promoter region, some of which were specifically related to hormone and drought responsiveness ( Figure 5 ). Thus, MsABC1K gene expression is potentially controlled by diverse light-responsive regulatory elements, phytohormones, and abiotic stresses. Different MsABC1K genes had similar cis-acting elements, which reflected the evolutionary conservation of the MsABC1K gene family.

To clarify the MsABC1K expression pattern in different alfalfa tissues, transcriptome sequencing was performed to analyze MsABC1K gene expression in alfalfa stems, roots, leaves, and flowers ( Figure 6 ). Twenty-one MsABC1K genes were expressed in all four alfalfa tissues; however, MsABC1K2b was expressed only in flowers. MsABC1K gene expression abundance significantly varied among the different tissues. The MsABC1K genes MsABC1K1/3a/3b/4/5a/6/7/8a/8b/9 exhibited preferential expression in leaves and flowers, while the MsABC1K genes MsABC1K5b/5c/12/15 exhibited preferential expression in stems and leaves. Additionally, the MsABC1K genes MsABC1K2a/10a/10b/10c/11/13/14 showed similar expression patterns in different tissues. These results demonstrate that MsABC1K genes perform distinct functions during normal development and growth.

To examine whether MsABC1K genes participate in salt stress response in alfalfa, we conducted qRT-PCR to detect the leaf expression profiles of 12 MsABC1Ks belonging to four different subfamilies following treatment with 300 mM NaCl for various time periods: 0 h, 3 h, 6 h, and 12 h ( Figure 7 ). Six MsABC1K genes (MsABC1K2a/10a/11/13/14/15) showed a significant increase in their expression levels; among these genes, the expression level of MsABC1K14 was 16-fold at 12 h after NaCl treatment. Four MsABC1K genes (MsABC1K4, MsABC1K5b, MsABC1K7, and MsABC1K9) exhibited significantly increased expression in the early stage, which finally decreased; in contrast, MsABC1K1 and MsABC1K12 gene expression was remarkably suppressed after salt stress exposure. Moreover, following treatment with 300 mM NaCl, 12 MsABC1K genes displayed relatively complicated expression patterns at 3 h, 6 h, and 12 h, possibly indicating a different regulatory function of each MsABC1K gene in alfalfa after exposure to salt stress.

To better comprehend the response of MsABC1K genes under drought stress, we conducted qRT-PCR analysis of alfalfa following treatment with 20% (w/v) PEG-6000. As shown in Figure 8 , 12 MsABC1K genes responded to drought stress. These genes were clustered into two groups according to expression patterns: group I included two genes (MsABC1K1 and MsABC1K12) with significant downregulation in drought stress, and group II included 10 MsABC1K genes with upregulated expression in drought stress. The 10 MsABC1K members in the second group showed peak expression levels at varying time points following PEG-6000 treatment. In this case, the expression abundance of the MsABC1K genes MsABC1K2a/4/5b/7/9/10a/11/13/15 peaked at 3 h following drought stress exposure, while MsABC1K14 expression abundance peaked at 6 h following drought stress treatment. These findings confirmed that MsABC1K genes are generally involved and have critical functions in alfalfa’s response to abiotic stress.