MYB transcription factors in alfalfa (Medicago sativa): genome-wide identification and expression analysis under abiotic stresses

Identification of the MYB gene family in alfalfa

To identify more comprehensive MYB TF family genes, the sequences of 168 A. thaliana, 430 Glycine max, 185 M. truncatula, and 130 O. sativa MYB proteins were downloaded from the Plant Transcription Factor Database (http://planttfdb.cbi.pku.edu.cn/) (Jin et al., 2017). These sequences were used as queries in a BLAST search against the M. sativa proteome sequences provided by the Noble Research Institute (https://www.alfalfatoolbox.org/) (Altschul et al., 1990), with an E-value cut-off of 0.00001. The obtained MYB sequences were confirmed based on the presence of intact MYB domains using the Pfam Program (http://pfam.xfam.org/), and the expectation cut-off (E value) 1.0 was set as the threshold for significance (El-Gebali et al., 2019). Moreover, the remaining sequences were analyzed by the cluster database at high identity with tolerance (CD-HIT) web server (http://www.bioinformatics.org/cd-hit/) (Li & Godzik, 2006), using default parameters to remove redundant data. Subsequently, non-redundant sequences were renamed as predicted MsMYB genes. The grand average of hydropathicity (GRAVY) index values and theoretical isoelectric point (pI) of these predicted MsMYB proteins were determined by the ProtParam Tool (https://web.expasy.org/protparam/) (Gasteiger et al., 2003). Additionally, subcellular localization of MsMYB genes was predicted using the Target P 1.1 server (Emanuelsson et al., 2000), and validated in WoLF-PSORT (https://www.genscript.com/wolf-psort.html) (Horton et al., 2007).

In silico functional analysis of MsMYB genes

The MsMYB protein function was predicted by (gene ontology) GO annotation, using the web-accessible Blast2GO v4.1 annotation system (https://www.blast2go.com/) (Gotz et al., 2008). Briefly, the MsMYB protein sequence was used to search for similar sequences against the NCBI non-redundant (Nr) database using the Blast tool in the Blast2GO software, with an expectation value of 10⁻³. Next, mapping and annotation were performed on Blast2GO using default parameters. Moreover, GO functional classification was performed by WEGO 2.0 (Ye et al., 2018). In addition, GO enrichment analysis for MsMYB genes was conducted on agriGO v2.0 (http://systemsbiology.cau.edu.cn/agriGOv2/) with default parameters (Tian et al., 2017), and the A. thaliana gene model (TAIR9) was selected as the reference.

Conserved motif and phylogenetic analysis of the MsMYB genes

In order to analyze the sequence features of MYB repeats in R2R3-MYB proteins, the amino acid sequences of R2 and R3 repeats of all R2R3-MYB proteins in M. sativa were extracted, and multiple sequence alignments of these identified R2R3-MYB proteins were performed using MUSCLE with default parameters. The sequence logos for R2 and R3 repeats were generated using the WebLogo (http://weblogo.berkeley.edu/logo.cgi) with default settings (Crooks et al., 2004). In addition, the motifs conserved among MsMYB members were identified using the multiple expectation maximization for motif elicitation (MEME) v4.11.1 online program (http://meme-suite.org/index.html) (Bailey et al., 2009). The maximum number of motifs was set to 45, according to a previous report (Zheng et al., 2017), and other parameters used the default settings. Moreover, to clarify the evolutionary relationship of these predicted MsMYB TF family genes, multiple sequence alignments of MYB protein sequences were performed using Clustal W (http://www.clustal.org/clustal2/) with default parameters (Larkin et al., 2007). The subsequent phylogenetic analysis relied on the neighbor-joining (NJ) method, as implemented in the MEGA v6.0 software (https://www.megasoftware.net/) (Tamura et al., 2013), and a bootstrap analysis was applied based on 1,000 replicates. Additionally, nonsynonymous (Ka) and synonymous (Ks) substitution rates were calculated to explore the mechanism of gene divergence after duplication, and Ka and Ks were computed using DnaSP 5 software (Librado & Rozas, 2009).

In silico expression analysis of the MsMYB genes during plant development

Genome-wide transcriptome data from M. sativa in different tissues during development were downloaded from the CADL-Gene Expression Atlas (https://www.alfalfatoolbox.org/atlasCADL/) provided by the Noble Research Institute (O’Rourke et al., 2015). The transcriptome data were derived from six tissues, including leaf, flower, pre-elonged stem, elonged stem, root, and nodule. Subsequently, these expression data were analyzed and clustered using the hierarchical cluster program MEV 4.9.0 (http://www.mybiosoftware.com/mev-4-6-2-multiple-experiment-viewer.html) to draw a heatmap of MsMYB genes in different tissues during development.

In silico transcriptome analysis of the response of the MsMYB genes to abiotic stress

In previous studies, a total of four transcriptome sequencing projects were performed to obtain additional genetic information for alfalfa relating to the response to abiotic stress, including cold (SRR7091780–SRR7091794, Zhou et al., 2018), and ABA, drought, and salt treatments (SRR7160313–SRR7160357, Luo et al., 2019a, 2019b). In this study, the expression level of some MsMYB genes (with the absolute value of fold change ≥ 2) was obtained after the local nucleotide blast (BLASTN) against these four transcriptome data was performed (Camacho et al., 2009). The clustering analysis and the heat map generated were performed by the hierarchical cluster program MEV 4.9.0.

In silico analysis of cis-regulatory element

Sequences of 2,000 bp from promoters of these changed MsMYB genes during abiotic stresses were analyzed for potential cis-regulatory elements, by querying them through the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). A total of six cis-regulatory elements were recorded, including abscisic acid responsive (ABRE), CGTCA-motif (methyl jasmonate responsive), G-box (light inducible), low-temperature responsive (LTR), MBS binding site (drought responsive), and TC-rich repeats (defense and stress responsive).

Plant materials and stress treatments

The plant materials used were alfalfa variety Zhongmu No. 1, which was cultivated and provided by the Qingchuan Yang Laboratory of the Beijing Institute of Animal Sciences, Chinese Academy of Agricultural Sciences. The experimental samples in this study were obtained by a hydroponic experiment. Before germination, vernalization of seeds was conducted at 4 °C to maintain the consistency of germination. After 3 days of seed germination, the seedlings were transferred into a nutrient solution (1/2 MS, pH = 5.8) and grown under a 16 h light/8 h dark cycle at 22 °C. Different stress treatments of seedlings were performed when the third leaf of the alfalfa appeared (approximately 10 days after germination). For cold treatment, the seedlings were placed in an artificial climate incubator and frozen under a 16 h light/8 h dark cycle at 4 °C. Four cold-treatment time points were used (2, 6, 24, and 48 h), along with one control. In the drought or salt treatments, the seedlings were transferred into a nutrient solution containing 400 mM Mannitol or 250 mM NaCl, and grown under a 16 h light/8 h dark cycle at 22 °C. There were eight treatments in these two experiments, which included seven treatment time points (1, 3, 6, 12, 24 h, and stress removal 1 and 12 h) and one control. For ABA treatment, the seedlings were transferred into a nutrient solution containing 10 μM ABA, and a total of three ABA treatments were performed, including 1, 3 and 12 h treatment time points. In order to be consistent with the experimental materials of transcriptome sequencing in our laboratory and obtain tissues closely related to abiotic stress, the whole seedling was harvested for cold treatment, and the root tip was harvested for the ABA, drought, and salt treatments. A total of six seedlings were collected and pooled into a frozen tube pipe, one for each treatment at the corresponding time point, and these samples were flash-frozen in liquid nitrogen and stored at −80 °C.

Quantitative real-time PCR analysis

Ribonucleic acid was extracted from the whole seedlings (cold treatment) or root tips (ABA, drought, and salt treatments) of control and treated seedlings using the Trizol method (Sangon Biotech, Shanghai, China), according to the manufacturer’s instructions. The concentration of each sample was determined using a NanoDrop ND1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). A one μg aliquot of DNase-digested total RNA was used to synthesize the single-strand cDNA using a FastKing RT Kit (Tiangen Biotech, Beijing, China), following the manufacturer’s protocol. The subsequent qRT-PCR was performed using a TB Green™ Premix Ex Taq™ Kit (TaKaRa, Dalian, China) on a CFX96 Real-Time PCR Detection System (Bio-Rad, Los Angeles, CA, USA). Each 20 μL reaction contained 10 μL TB Green Premix, 0.2 μM of each primer, and two μL cDNA. The reactions were initially denatured (95 °C/30 s), then subjected to 40 cycles of 95 °C/5 s, 60 °C/30 s. According to the transcriptome sequence of Zhongmu No. 1, gene-specific primers were designed using Primer Premier 6 software (Premier Biosoft International, Palo Alto, CA, USA), and are shown in Table S1. Each of the three biological replicates were supported by three technical replicates, and the relative expression levels were normalized to the expression of the Medicago actin gene (AA660796) (Liu et al., 2016, 2017a) and calculated using the 2^−ΔΔCt method.

Protein–protein and protein-DNA interactions predictions

The predictions of protein–protein and protein-DNA interactions for MsMYB genes were carried out by the online server Arabidopsis Interactions Viewer (http://bar.utoronto.ca/interactions/cgi-bin/arabidopsis_interactions_viewer.cgi), which queries a database of 70,944 predicted and 39,505 confirmed A. thaliana-interacting proteins. The expression patterns of genes that interact with MsMYB genes during abiotic stresses were analyzed using the hierarchical cluster program MEV 4.9.0, and the correlation coefficient between the expression levels of the MsMYB gene and its interacting genes was calculated. In addition, the MYB-core motifs (C/TNGTTA/G) were searched in the upstream sequence (2,000 bp) of the homologous gene of interacting genes in M. truncatula, and these interacting genes were obtained by the prediction of protein-DNA interactions. To verify the reliability of these predicted protein–protein interactions, a total of five predicted protein–protein interactions were selected (Table S2). The recombinant pLexA and pB42AD plasmids were co-transformed into yeast strain EGY48 and plated on medium lacking Ura, His, and Trp (-U-H-T) at 30 °C. The colonies were transferred to induction medium lacking Ura, His, Trp, or Leu (-U-H-T+X-gal and -U-H-T-L+X-gal) for interaction screening.

Identification and classification of the alfalfa MYB gene family

To identify MYB genes in the M. sativa genome, a BLASTP search was performed using 168 A. thaliana, 430 Glycine max, 185 M. truncatula, and 130 O. sativa MYB protein sequences as queries. Then, a Pfam search was used to ensure that they contained the MYB DNA-binding domain, for verification of the identity of these MYB sequences. Additionally, these putative MYB genes were inspected using CD-HIT to ensure they possessed complete open reading frames and maps to unique genomic loci. As a result, a total of 265 non-redundant MYB sequences were obtained and named MsMYB001 to MsMYB265, based on the order of their serial number in the database of the M. sativa proteome sequences. These MYB TFs were classified into four distinct groups, namely, “R1-MYB,” “R2R3-MYB,” “R1R2R3-MYB,” and “atypical-MYB” genes, based on the presence of one, two, three, or four MYB repeats, respectively. Our analysis revealed that the R2R3-MYB subfamily consisted of the highest number of MYB genes, and the number of these four types of MYB genes was 50 (18.87%), 186 (70.19%), 26 (9.81%), and 3 (1.13%), respectively (Table 1). The classification information of each gene is shown in Table S3. In a previous study, 17 MsMYBs were predicted and analyzed (Dong et al., 2018). In order to verify the accuracy of predicted MsMYB genes in our study, a BLAST search for 17 reported MsMYBs in the 265 predicted MsMYBs was performed. As a result, 15 of the 17 reported MsMYBs were found in these 265 MsMYBs, and the details are shown in Table S4. However, the other two reported MsMYBs were not found, which may be due to the inconsistencies in the sequencing of samples between these two sequencing projects provided by Dong et al. (2018) and the Noble Research Institute.

In addition, the physiochemical properties of these MsMYB proteins were also analyzed, including protein length, molecular weight, pI, GRAVY, and subcellular localization (Table S3). The 265 predicted MsMYB proteins ranged from 79 (MsMYB259) to 1,879 AA (MsMYB173) in length, with an average length of 360.76 AA (Table 1). The molecular weight of these MsMYB proteins ranged from 9,130.13 D in MsMYB047 to 208,302.91 D in MsMYB173. Moreover, the mean of pI and GRAVY was 7.19 and −0.71, respectively. We also predicted the subcellular localization of these MsMYB proteins using several localization predictor software programs. Our analysis revealed that 232 (87.55%) MsMYB proteins were found to be nuclear-localized (Table S3). The remaining 33 MsMYB proteins were predicted to be localized in the chloroplast (15), cytoplasm (11), plasma membrane (2), vascular bundle (2), endoplasmic reticulum (1), golgi apparatus (1), and mitochondria (1).

In silico functional classification of MYB transcription factors

For the functional annotation of all MsMYB genes, GO functional analysis was performed by the Blast2GO program. In this study, a total of 23 GO categories were assigned to the 265 MsMYB genes (Fig. S1; Table S5). “Cell” gene (62, 23.4%) was the dominant category in the cellular component category, followed by “cell part” and “organelle” (61 for both, 23.0%). In the molecular function category, a total of 235 (88.7%) genes were assigned to “binding,” including “organic cyclic compound binding” (234, 88.3%), “heterocyclic compound binding” (234, 88.3%), and “ion binding” (4, 1.5%). Regarding the biological process category, “cellular process” (50, 18.9%) and “metabolic process” (48, 18.1%) were the most dominant groups. Additionally, we performed a GO enrichment analysis of the functional significance, using the agriGO website to identify the significantly enriched GO terms among these 265 MsMYB genes, with a p score cut-off of 0.05. As a result, a total of eight GO terms were considered to be significantly enriched among these genes (Table 2), and four and three GO terms belonged to “molecular function” (F) and “biological process” (P), respectively. However, only one GO term was significantly enriched in “cellular component” (C), namely “nucleus” (GO:0005634).

Phylogenetic analysis of MsMYBs

In this study, we performed phylogenetic analysis of M. sativa and M. truncatula MYB proteins using the NJ method to evaluate the evolutionary significance of the MYB genes. The circular phylogenetic tree was constructed based on the alignment of the amino acid sequences of 310 MYB genes, including 265 from M. sativa and 45 from A. thaliana, L. japonicas, and M. truncatula (15 genes for each species) (Table S6). Based on the topology and robustness of the NJ phylogenetic tree results in M. truncatula (Wang & Li, 2017), we resolved these 265 MsMYB genes into 12 subgroups (A–L), which ranged in size from three (H) to 52 (L) (Fig. 1).

In silico analysis of the conserved MYB domains in MsMYB genes

In order to investigate the sequence features of R2R3-MYB proteins in M. sativa, multiple sequence alignments were performed with Clustal W, and a highly conserved 104 AA region among all R2R3-MYB proteins was identified. Subsequently, the WebLogo program was used to study the variation within the conserved motifs of R2R3-MYB genes in alfalfa. As a result, the R2 and R3 MYB repeats of the MsR2R3-MYBs contain five highly conserved Trp (W) residues that play a key role in sequence-specific binding of DNA (Fig. 2). Of these five Trp residues, the third conserved tryptophan residue (W₄₆) in the R2 repeat was not completely conserved in all of the MsR2R3-MYBs. Moreover, we also found that some amino acids showed high conservativeness, such as Glu (E), Asp (D), Leu (L), Arg (R), Lys (K), Ser (S), Cys (C), Gly (G), and Asn (N). Furthermore, the protein motif organization among all 265 MsMYBs was investigated using MEME. The MEME results showed that the width of 45 identified motifs range from 14 to 50 AA, and the maximum E-value of these motifs was 1.0e-041 (Table S7). In order to validate the clade definition based on phylogenetic analysis, the protein motif distribution for each MsMYB was analyzed. As a result, one or more motifs outside of the MYB domain were shared by most members within the same clade, which provides further support for the results of the phylogenetic analysis. The distribution of the conserved motifs within each clade is shown in Fig. S2. In order to explore the evolutionary patterns, divergence and selection pressure of the MsMYBs, a total of 44 homologous pairs were identified using phylogeny-based and bidirectional best-hit methods. The Ka/Ks ratio is widely applied to measure genetic evolution and selection pressure. Scatter plot statistics showed that all gene pairs had evolved mainly under the influence of purifying selection except for three homologous pairs (MsMYB082/084, MsMYB235/253, and MsMYB256/257) (Fig. 3).

Expression profile analysis of the MYB genes in M. sativa tissues

Microarray datasets from the alfalfa B47 genotype were downloaded from the CADL-Gene Expression Atlas and used to assess the transcript abundance profiles of MsMYB-encoding genes in six major tissues: leaf, flower, pre-elonged stem, elonged stem, root, and nodule. We only examined the transcript abundance of 211 MsMYB genes because the remaining 54 genes were not represented in the dataset. In order to visualize the datasets, an expression heat map was obtained by the program MEV 4.9.0. According to the expression patterns of these genes over six tissues of the alfalfa B47 genotype, 211 MsMYB genes were divided into eight subgroups, namely A to H (Fig. 4). Subgroup A contains eight genes, and they showed the highest transcript accumulation level in nodule tissue. Subgroup B includes 40 MsMYB genes, which showed the lowest transcript accumulation in the root and nodule. In contrast, some of these genes showed the highest transcript accumulation level in the leaf, flower, or pre-elonged stem. Moreover, a total of 37 MsMYB genes were steadily expressed in all of the six tissues tested, and these MsMYB genes belong to subgroup C. Approximately 33 of these 211 MsMYB genes (subgroup D) showed the highest transcript accumulation level in flower tissue, 10 (subgroup E) showed the highest transcript accumulation in leaf tissue, 35 (subgroup F) showed the highest transcript accumulation in root or nodule tissue, and 30 (subgroup G) showed the highest transcript accumulation in the stem, including pre-elonged and elonged stems. The remaining 18 MsMYB genes (subgroup H) showed the lowest transcript accumulation in the leaf and flower.

Expression analysis of MsMYB genes in response to abiotic stresses

In order to investigate the expression levels of MsMYB genes under abiotic stresses, local nucleotide blast (BLASTN) against four transcriptome datasets reported by our laboratory was conducted. As a result, we determined that the expression of 51 and 52 MsMYB genes changed significantly during cold and three other stresses, respectively, and these 51 genes were present in all four treatments. Because tissue samples used for cold sequencing differed from the other three stresses, the expression patterns of these genes during cold treatment were analyzed separately from the other three treatments. As shown in Fig. 5, a total of 51 MsMYB genes were divided into four subgroups, namely A–D. Compared to the control samples, 11 MsMYB genes (subgroup B) were inhibited, and the remaining genes were induced, during cold stress. Interestingly, a total of 21 genes belonging to subgroup A and C were induced earlier by cold stress than those in subgroup D (19 genes). Moreover, the expression patterns of 52 MsMYB genes during ABA, drought, and salt treatments were also analyzed. Our analysis revealed that these MsMYB genes were also divided into four subgroups: A, B, C, and D. As shown in Fig. 6, subgroups A, B, and D were highly expressed in all three treatments. Additionally, subgroup C was highly expressed in response to drought and salt stresses, but most genes were not expressed in response to ABA treatment.

Cis-regulatory element in MsMYB gene promoters

Cis-regulatory elements control expression patterns of stress-responsive genes, and these elements are located upstream of gene-coding sequences and provide binding sites for TFs. Thus, we investigated the distribution of six cis-regulatory elements in 52 stress-responsive MsMYB gene promoters. As a result, a total of 107 ABRE cis-elements were identified, which was more than the number of the other five cis-elements (Fig. 7). Meanwhile, we found at least one cis-element in each MsMYB gene promoter, and only the MsMYB040 promoter contained all six cis-elements.

Gene expression analysis qRT-PCR validation

To further confirm the RNA-seq responses of these 52 MsMYB genes to abiotic stresses, qRT-PCR was performed for 14 MsMYBs that were significantly induced genes under various stress treatments (Figs. 5 and 6; Table S1), including ABA, cold, drought, and salt. The expression patterns of most of the MsMYB genes in the qRT-PCR analysis were consistent with the RNA-Seq analysis, but the magnitude of the fold changes varied between RNA-seq and qRT-PCR experiments (Fig. 8). Moreover, these results showed that all 14 genes were induced to different degrees by abiotic stress.

Protein–protein and protein-DNA interactions predictions

As a TF, the MYB gene can regulate plant development and responses to various environmental changes by controlling the expression of downstream genes, and is also regulated by the upstream genes. Therefore, the prediction of protein–protein and protein-DNA interactions was performed in this study to identify some candidate genes that putatively interact with the MsMYB genes, and a total of 178 and 957 predicted interactions were recognized, respectively. Subsequently, these predicted interaction genes were screened with the absolute value of fold change ≥ 2 as the threshold. As a result, a total of 170 protein–protein and 914 protein-DNA interactions were predicted (Tables S8 and S9, respectively), and the correlations between the MsMYB gene and its putative interactor’s genes are shown in Tables S10 and S11, respectively. Of the 914 predicted protein-DNA interactions, only one homologous gene promoter of interacting genes in M. truncatula does not contain any MYB-core motifs (C/TNGTTA/G) (Table S9). Additionally, there were 27 and 105 more reliably predicted protein–protein (fold change ≥ 2, correlation coefficient ≥ 0.8) and protein-DNA (fold change ≥ 10, correlation coefficient ≥ 0.8) interactions, and the expression patterns of these predicted interaction genes were affected to varying degrees by abiotic stresses (Figs. 9 and 10; Figs. S3 and S4). We found that most predicted protein–protein and protein-DNA interactions showed significant positive correlations during ABA or cold treatments, but negative correlations under drought or salt stresses (Fig. 11; Fig. S5). Furthermore, the yeast co-transformants of all five predicted protein–protein interactions grew on SD/-U-H-T medium (Fig. S6), and three yeast co-transformants turn blue on SD/-U-H-T+X-gal and SD/-U-H-T-L+X-gal induction mediums (Fig. 12), indicating that MsMYB043 and MsMYB253 strongly interacted with MSAD320162, and MsMYB253 also interacted with MSAD308489. However, the remaining two predicted protein–protein interactions didn’t show a significant relationship.