Genome-wide identification, characterization and expression analysis of HAK genes and decoding their role in responding to potassium deficiency and abiotic stress in Medicago truncatula

Identification and sequence analysis of MtHAKs

MtHAKs sequences were obtained from the Medicago truncatula genome databases (HAPMAP, https://medicagohapmap2.org). The amino acid (aa) sequences of Arabidopsis (TAIR, http://www.arabidopsis.org/) and rice (TIGR, http://rice.plantbiology.msu.edu/) HAKs were used as the reference sequences for searching predicted homolog sequences in M. truncatula using the HMMER3.0 software (http://hmmer.org/). Subsequently, the genes were screened using a threshold of <1e−100 E-value (full sequence and best one domain). Candidate protein members were verified using the SMART database (http://smart.embl-heidelberg.de/) and NCBI-Conserved Domain Database (CDD, https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) (Zhao et al., 2021), with proteins with shorter aa length (<400 aa) and those containing incomplete K⁺ transporter domains being discarded. The longest gene was chosen for further analysis only if it had alternative splicing variants. Subcellular localization of MtHAK proteins was predicted using the WOLF PSORT software (https://www.genscript.com/wolf-psort.html) and the TMHMM Server 2.0 online tool (https://services.healthtech.dtu.dk/service.php?TMHMM-2.0) was used for predicting the protein transmembrane helices.

Construction of MtHAKs phylogenetic tree

HAK protein sequences of Arabidopsis and rice were retrieved from the NCBI database (https://www.ncbi.nlm.nih.gov) (Table S1), while multiple sequence alignment was conducted using the ClustalW program (Version 2.1; http://www.clustal.org/). MEGA7.0 was used to construct the phylogenetic tree using the neighbor-joining method along with the bootstrap replicates being up to 1,000 (Liu et al., 2019; Liu et al., 2020).

Gene structure and conserved motif analysis

Gene structure and conserved motifs were visualized using the TBtools (Toolkit for Biologists integrating various biological data-handling tools) software (Chen et al., 2020). The conserved and identified motifs of protein sequences were predicted via the MEME (Multiple Expectation Maximization for motif Elicitation) program (Version 5.1.1), with the maximum protein motif number being set as 10, and the other parameters set as default (http://meme-suite.org/tools/meme) (Bailey et al., 2009).

Chromosomal location and synteny analysis

The MtHAK chromosomal location was illustrated by the circos diagram by annotating genes to their specific chromosomal location in their genome sequences by using the TBtools software. These syntenic analyses were carried out by using the MCScanX with gene duplication parameters (Wang et al., 2012).

Analysis of cis-acting regulatory elements in MtHAKs promoter regions

Putative cis-acting regulatory elements were analyzed using the PlantCARE online software (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). The 2.0 kb promoter sequences located upstream of the transcription starting site in each MtHAK gene were extracted from the M. truncatula genome database.

Analysis of microarray expression profile

The microarray data of the expression profiles of MtHAKs in the roots, vegetative bud, stem, petiole, leaf, flower, pods, and seeds and their responses to abiotic stress were obtained from the MtGEA (Benedito et al., 2008). When a gene corresponded to multiple probes, the maximum value of the probe was selected for the subsequent analysis. The normalized microarray data was used to create the heatmap through the TBtools software, based on the mean value of each gene expression in all the analyzed organs. The expression patterns of MtHAKs in response to salt, drought, and cold stresses were obtained from the NCBI under GEO accession number GSE136739 (Song et al., 2017). The expression abundance of each MtHAK gene was represented by fragments per kilobase million (FPKM). The relative stress-induced expression levels were calculated by comparing with the control samples. The clustered heatmap was generated using the TBtools software and based on their relative expression.

Stress treatment and qRT-PCR

For K⁺ deficiency stress treatment, two-week-old seedlings were grown in 1/2 Hoagland nutrient medium without K⁺ for 0 (control), 1, 6, 12, 24, and 48 h, respectively. For the salt stress treatment, two-week-old seedlings were grown in 1/2 Hoagland nutrient medium containing 300 mM NaCl for 0, 1, 6, 12, 24, and 48 h, respectively. For drought stress treatment, two-week-old seedlings were grown in 1/2 Hoagland nutrient medium containing 18% PEG6000 for 0, 1, 6, 12, 24, and 48 h, respectively. The root samples were subsequently cut, then snap frozen in liquid nitrogen, and finally stored at −80 °C until further use. The qRT-PCR analysis was performed in triplicates for each of the biological replicates. Their relative expression levels were calculated using the 2^−ΔΔCt analysis method (Liu et al., 2019; Zhao et al., 2022). The expression levels of the control samples were normalized to one, with the MtActin gene being used as the internal control. Standard deviations and the significant differences were indicated by error bars and an asterisk (*) (p < 0:05), respectively.

Identification of HAK members in M. truncatula

To identify M. truncatula HAK genes, we conducted a genome-wide search using the HMMER3.0 Software (http://hmmer.org/) based on the M. truncatula genome sequences along with the Arabidopsis and rice HAK genes as subjected queries. Then, we identified 20 nucleotide sequences with a typical canonical K⁺ transporter domain (Pfam accession no. PF02705) using the Pfam and SMART databases, and they were subsequently designated as MtHAK1 to MtHAK20 depending on their chromosomal positions (Table 1). Detailed information on the 20 HAK genes is listed in Table S1. The number of protein transmembrane segments (TMS) ranged between 10 and 13, with the most common being 12–13 (70%). All the examined HAK proteins were predicted to be mainly localized in the plasma membrane using a PSORT analysis (http://www.psort.org). The protein length of the 20 identified HAK proteins ranged from 619 aa (MtHAK3) to 856 aa (MtHAK2) with an average length of 778 aa. Their relative molecular weights (MW) varied from 69.03 kDa (MtHAK3) to 95.67 kDa (MtHAK2). The isoelectric points (pI) ranged from 5.44 (MtHAK8) to 9.39 (MtHAK19).

HAKs phylogenetic relationship among M. truncatula, Arabidopsis and rice

To analyze the evolutionary relationships of the MtHAK proteins, we conducted phylogenetic analyses of 60 HAK amino acid sequences (20 – M. truncatula, 13 – Arabidopsis, and 27 – rice) to construct a phylogenetic tree using the neighbor-joining method. According to the evolutionary tree, we classified all HAK members into four major groups: Groups I–IV. Furthermore, we classified the MtHAK proteins into three clusters (from I to III): Cluster I (MtHAK6, 14, 15, 17, and 18), Cluster II (MtHAK4, 7, 9, 10, 11, 12, 13, 19, and 20), and Cluster III (MtHAK1, 2, 3, 17, 5, and 8) (Fig. 1). All members in group IV belong to rice. The most members existed in Cluster II in M. truncatula, thus comprising 45% of all MtHAKs. The phylogenetic tree showed that MtHAKs were most closely related to Arabidopsis KUPs than those of rice HAKs, thereby indicating that MtHAKs might share evolutionary functional similarities with Arabidopsis KUPs. All MtHAKs in cluster I were distributed together with the already-identified AtHAK5, which suggested that they may be crucial for K⁺ uptake from a low-K⁺ level soil (Lara et al., 2020). Among cluster II members, MtHAK4 and MtHAK19 shared high sequence identity with AtKUP2 (Elumalai, Nagpal & Reed, 2002), and AtKUP4 (Rigas et al., 2001; Vicente-Agullo et al., 2004), respectively, thus implying they are likely to be involved in plant development processes. Additionally, among cluster III, MtHAK1 and MtHAK8 clustered together with AtKUP7 (Han et al., 2016), thereby suggesting their role in K⁺ acquisition and translocation under low K⁺ concentration.

Gene structure and motif composition of MtHAK genes

MtHAK proteins were listed in order based on the phylogenetic analysis (Fig. 2A), which was consistent with the results in Fig. 1. Closely related members shared similar exon/intron structures, which were related to their biological functions. Gene structures of the MtHAKs were abtained based on the arrangement of the untranslated region, exon, and intron sequences generated using the TBTools software. As shown in Fig. 2B, the exon number of MtHAK genes varied from 8 to 10, and the longest exon existed in the end of a gene except for MtHAK2, which is consistent with previous reported data (He et al., 2012; Hyun et al., 2014). Additionally, most MtHAKs in the same cluster shared high exon-intron structure similarity (Fig. 2B).

To study the structural features, we analyzed conserved protein motifs of MtHAKs using the MEME program. We identified the conserved protein motifs varying from 29 to 50 aa in lengths and designated them as motifs 1–10. Conserved protein motif information is shown in Table S2. The highly conserved K⁺ transport domain (GVVYGDLGTSPLY), included in motif 9, existed in all MtHAK proteins besides MtHAK10 (Fig. 2, Table S2). Motifs 1, 2, 3, 4, 5, 6, 7, 8, and 10 were almost evenly distributed along with a feature domain of K⁺ transporter (Fig. 2C, Table S2) in all the MtHAK proteins. Therefore, the motifs of conserved K⁺ transporter and similarities of gene structure in the same cluster together implied the closing function among these HAK members.

Chromosomal distribution and synteny analysis of MtHAK genes

All identified MtHAK genes were mapped onto chromosomes from the M.truncatula genome database to identify and locate their chromosomal distribution. Results showed that MtHAKs were distributed on seven of the eight chromosomes, with chromosome 8 containing the highest number of six MtHAK genes (Fig. 3). Five MtHAK genes were located on chromosome 5, three on chromosome 2, two on chromosome 4 and 6, one on chromosome 3 and 7, and no gene was allocated on chromosome 1 (Fig. 3). These results indicated that MtHAKs were scattered randomly onto different chromosome locations.

We further performed synteny analysis between M. truncatula and Arabidopsis to verify the evolutionary relationships and history of the MtHAKs. Subsequently, we found seven collinear gene pairs between M. truncatula and Arabidopsis in the dataset (Fig. 3 and Table S3). This indicated that these identified genes might already have existed before protein structure divergence, thereby further implying a strong phylogenetic relationship. Furthermore, only one gene pair (MtHAK2/MtHAK5) existed as paralogs in M. truncatula.

Analysis of cis-acting elements in the promoter region of MtHAK genes

To further investigate the gene function and regulatory mechanism of MtHAKs, we analyzed the 2 kb regions upstream of the translation start site of the 20 MtHAK genes using the PlantCARE database. We identified 73 putative cis-elements in the MtHAK promoters based on functional annotation, and the major types of cis-elements are shown in Fig. 4 and Table S4. Post analysis, we identified cis-elements corresponding to different plant hormones like auxin (TGA-element and AuxRE-core), gibberellin (GARE-motif and P-box), MeJA (TGACG-motif and CGTCA-motif), ethylene (ERE-box), ABA (ABRE), and salicylic acid (TCA-element), in the promoter regions of all MtHAKs genes except MtHAK20, thereby suggesting that MtHAKs expression may be regulated by different phytohormones. Furthermore, we also found abiotic stress-responsive elements, including STRE, ARE, WRE3, WUN-motif, MBS, LTR, DRE-core, DRE1, and TC-rich repeats, in all the MtHAKs promoter regions except MtHAK20. Additionally, zein metabolism regulation element (O2-site), endosperm expression element (GCN4-motif and AACA-motif), palisade mesophyll cells element (HD-Zip 1), meristem expression element (CAT-box and CCGTCC-motif), and seed regulation element RY-element were also abundant in the promoter of MtHAKs except MtHAK20. However, the MtHAK20 promoter region was abundant in light-responsive elements (Table S4).

Spatial expression profiles of MtHAK genes

To gain further insights into the potential biological function of MtHAK genes, we used the publicly available microarray data of the Medicago truncatula Gene Expression Atlas (MtGEA, https://mtgea.noble.org/v3/) to investigate the temporal and spatial expression pattern of the MtHAKs. MtHAK4 showed relatively high expression in all tissues, while that of MtHAK18 was low in all tissue (Table S5). Notably, MtHAK6 and MtHAK16 were expressed preferentially in the roots, thereby implicating their role in K⁺ uptake from the soil (Fig. 5, Table S5).

Both cluster III genes, MtHAK2 and MtHAK3 exhibited similar expression patterns and relatively high expression in leaves. MtHAK13 was exclusively and highly expressed in floral organs, whereas MtHAK8 showed the same in immature seeds (Fig. 5, Table S5). Interestingly, MtHAK5 and MtHAK12 exhibited high and gradually increased expression patterns during the reproductive stages and finally peaked at 24 days after pollination (DAP) Contrastingly, MtHAK15 was specifically highly expressed in immature seeds (10 DAP) with the expression pattern gradually decreasing along with seed maturation. Therefore, the spatial and temporal expression profiles indicated the functional diversity of MtHAK genes in Medicago trunculata development.

Expression patterns of MtHAK genes under K⁺ deficiency

Due to the major function of the HAK family being K⁺ transport, we investigated the expression profiles of MtHAK genes in the roots under K⁺ deficient conditions using qRT-PCR. As shown in Fig. 6, among the 20 MtHAK genes, we obtained eight genes that showed upregulated expression patterns post K⁺ deficiency treatment. MtHAK6, MtHAK7, and MtHAK17 expression slightly increased and finally peaked at 48 h post treatment. MtHAK15 and MtHAK18 showed nearly the same expression pattern at the five different time points. MtHAK9, MtHAK10, and MtHAK11 transcripts were strongly upregulated at 6 h, then peaking at 12 h and 24 h, and finally went down at 48 h. Therefore, these results suggested that these MtHAK genes were K⁺ deficiency-responsive. Furthermore, it was noteworthy that MtHAK6 was highly and specifically expressed in Medicago trunculata roots and also significantly upregulated in response to K⁺ deficiency.

Expression patterns of MtHAK genes under salt and drought stresses

Several HAK genes have been reported to participate in abiotic stresses (Elumalai, Nagpal & Reed, 2002; Vicente-Agullo et al., 2004; Chen et al., 2015; Shen et al., 2015). To verify this hypothesis, we evaluated the expression profiles of eight K⁺ deficiency responsive genes via qRT-PCR under salt and drought stress treatments. The results revealed that all eight genes were induced by salt and drought stresses to different extents (Figs. 7 and 8).

We determined the expression profile of MtHAK genes in Medicago trunculata roots at different times (0, 1, 6, 12, 24, and 48 h) under salt treatment (300 mM NaCl in nutrient solution). The results showed that expression of MtHAK7, MtHAK9, MtHAK15, MtHAK17, and MtHAK18 exhibited significant upregulation. Interestingly, MtHAK7 and MtHAK18 were quickly and continuously upregulated from 1 h and subsequently increased at 48 h (Fig. 7).

Additionally, we analyzed the expression profiles of MtHAK genes in Medicago trunculata roots under drought treatment simulated by 18% PEG6000 at different times (0, 1, 6, 12, 24, and 48 h). Under drought treatment, all selected genes besides MtHAK9 were upregulated, albeit to different levels at different times (Fig. 8). In particular, MtHAK10, MtHAK15, and MtHAK18 rapidly responded to dehydration at 1 h. Contrastingly, MtHAK17 was moderately upregulated from 6 to 48 h. Both MtHAK6 and MtHAK7 exhibited highly induced expression at 24 h.

Interestingly, we found that MtHAK15, MtHAK17, and MtHAK18 were strongly upregulated by both salt and drought stresses. The expression level of MtHAK18 increased rapidly at 1 h as compared to the control, under both salt and drought treatments (Figs. 7 and 8).