Aluminum (Al) is a light metal and the third most abundant element in the earth's crust (Ma, 2005). When the pH in the soil is lower than 5.0, Al is dissolved from the harmless form (oxide or aluminosilicate form) into the soil solution, mostly in the form of Al³⁺, and results in phytotoxic to most herbaceous plants even at low concentrations. In the tropical and subtropical regions, Al toxicity has been considered as the major factor limiting crop production in acidic soils, which account for 40% of the world's arable land (Kochian et al., 2005). Traditionally, the application of large quantities of lime always used to alleviate the soil Al toxicity and then sustain the crop production. However, this practice is expensive and being unsustainable and not environmentally friendly. Thus, understanding the nature of Al tolerance mechanisms in plants and then developing cultivars with improved tolerance to acidic soil stress is an appealing approach to addressing this issue.

Up to data, they are two main Al tolerance mechanisms in plants, namely, the exclusion mechanism and tolerance mechanism, have been proposed. The secretion of certain organic acids, such as citrate, oxalate and malate, are induced by Al is one of the best proved characterization in exclusion mechanism in plants. In Al stress conditions, organic acids able to form strong complexes with Al³⁺ and then preventing the binding of Al to cellular components to alleviate Al toxicity (Ma et al., 2001; Ma, 2005). Many agriculturally important Al-tolerant plant species, such as wheat (Triticum aestivum) (Delhaize et al., 1993), snapbean (Phaseolus vulgaris) (Miyasaka et al., 1991), maize (Zea mays) (Pellet et al., 1995), Cassia tora (Ma et al., 1997a), soybean (Glycine max) (Yang et al., 2000), buckwheat (Fagopyrum esculentum) (Ma et al., 1997b), taro (Colocasia esculenta) (Ma, 1998), and rye (Secale cereal) (Li et al., 2000) release one or two of these three organic acids in response to Al stress. Moreover, the importance of high level of ascorbic acid to Al tolerance also has been indicated in transgenic tobacco (Nicotiana tabacum) (Yin et al., 2010). Through genetic and molecular analyses, various functional genes have been identified and confirmed as important components in Al tolerance. ALMT1 is an Al-activated malate transporter gene identified from wheat. Overexpression of this gene in barley confers an Al-activated efflux of malate and results in the Al tolerance both in hydroponic culture and acid soil (Delhaize et al., 2004). When a Pseudomonas aeruginosa derived citrate synthase (CS) gene transformed into tobacco genome, higher citrate synthase activity, citrate efflux and greater Al resistance are observed in transgenic lines (de la Fuente et al., 1997). AtALS3 (aluminum-sensitive 3) encodes a phloem-localized ABC transporter-like protein, which is required for Al resistance/tolerance in Arabidopsis by redistributing accumulated Al³⁺ away from sensitive tissues, such as root, and thus reducing the toxic effects of Al (Larsen et al., 2005). In rice (Oryza sativa), the disruption of OsASR5 gene resulted in hypersensitivity to Al toxicity, and which may function as a transcription factor to protect rice cells from Al toxicity by regulating the expression of various genes (Arenhart et al., 2013). Recently, Yang et al. (2014) has shown auxin is responsible for the Al-induced inhibition of root growth and acts as the downstream of ethylene-regulated TAA1 expression in the root-apex transition zone.

Considering its complexity, it is essential to interpret the functional elements and molecular constituents involved in Al tolerance mechanisms on a whole-genome level in plants. Using microarray technology, a large number of Al-responsive genes in many plant species including Arabidopsis (Kumari et al., 2008), soybean (Duressa et al., 2010, 2011; You et al., 2011), wheat (Houde and Diallo, 2008), and Medicago truncatula (Chandran et al., 2008) have been identified. In addition, the recently developed high-throughput sequencing (RNA-Seq) has clear advantages over microarray methods and has been considered as the ideal option to discover new genes and estimate transcript abundance at genome-wide scale, especially useful for species without genome sequence (Trapnell et al., 2012; Zeng et al., 2015). Based on RNA-Seq platforms, genome-scale transcriptome analyses have been used to identify Al-stress-responsive genes in rice (Arenhart et al., 2014), buckwheat (Fagopyrum tataricum) (Yokosho et al., 2014; Zhu et al., 2015), maize (Mattiello et al., 2014), Anthoxanthum odoratum (Gould et al., 2015), and Medicago truncatula (Chen et al., 2012). These Al-stress-responsive genes identified by RNA-Seq are involved in many physiological and metabolic processes, such as protection against cell wall toxicity and oxidative stress, organic acid metabolism, and exudation, Al transportation, and hormone signal transduction.

Alfalfa (Medicago sativa L.) is the most extensively cultivated forage legume and plays vital ecological and economic roles in agricultural systems around the world (Liu et al., 2016). However, alfalfa is very sensitive to soil acidity, which greatly limits its productivity and persistence performance (Khu et al., 2013). Thus, a better understanding of the molecular mechanisms involved in alfalfa responses to Al stress would be critical for Al-tolerant alfalfa breeding programs. Previous studies have shown that overexpression of endogenous malate dehydrogenase or bacterial CS result in enhanced organic acid synthesis, Al secretion and Al resistance (Tesfaye et al., 2001; Barone et al., 2008). In addition to these transgenic studies, a proteomic analysis of alfalfa after Al treatment has revealed that leaf proteins responding to the Al stress are mainly involved in energy metabolism and antioxidant/reactive oxygen species (ROS) detoxification pathways (Rahman et al., 2014). Given the relatively low-throughput characteristics of proteomics-based approaches and the complex nature of the Al stress and resistance mechanisms, it is necessary to identify and functionally characterize Al-responsive genes in alfalfa on a genome-wide scale. Currently, to our knowledge, the genome-wide transcriptomic analysis of the Al-responsive genes in alfalfa has not been reported, especially within the root tips, where the primary site for Al detoxification and accumulation. Thus, in the present study, we carried out the first global transcriptome analysis of the alfalfa root tips during Al treatments using the Illumina RNA-Seq method. The results obtained in this study will extend the knowledge of the genetic basis of alfalfa response to Al stress at the transcriptional level.

Plants frequently encounter Al stress in acid soils and have thus evolved a series of responses and adaptive mechanisms to cope with Al stress. Understanding the molecular mechanisms underlying Al stress responses is important for Al-tolerant crop breeding program. Using microarray and newly developed deep sequencing technologies, the transcriptomic responses of many plants to Al stress have been comprehensively documented; however, most of these efforts have focused on model plants, such as Arabidopsis (Kumari et al., 2008), soybean (You et al., 2011), maize (Maron et al., 2008; Mattiello et al., 2014), and rice (Arenhart et al., 2014). In the present study, we analyzed the transcript profiles of alfalfa roots in response to Al stress by using the Illumina deep sequencing system and identified a total of 75,903 unigenes in the four sample libraries, which was more than previously reported for alfalfa root transcriptome analyses (Postnikova et al., 2013). Of these unigenes, more than 85% had significant similarity (BLAST, E ≤ 10⁻⁵) to genes in the public databases (Table 2). Using the more stringent criteria of both FDR ≤ 0.0001 and an expression difference greater than four-fold, our results detected 2464 Al stress-related DEGs after treatment with Al for 24 h compared with the control, thus indicating that these genes responded to Al stress in alfalfa. We further classified these Al stress-related DEGs into three groups on the basis of their expression patterns. As shown in Figure 4, most DEGs were up-regulated, a result consistent with those from previous reports in other plants under Al stress (You et al., 2011; Yokosho et al., 2014; Chen et al., 2015). The number of genes down-regulated by Al remained nearly constant over the course of the treatment, whereas the number of up-regulated genes had a dramatically different trend in which many responsive genes were observed at the early (4 h) and/or late (24 h) Al time points compared with those observed at the middle time point (8 h) (Figures 2A,B). These results contrast with those from previous microarray and RNA-Seq analyses of plant response to Al and other types of abiotic stress (Kumari et al., 2008; Xu et al., 2014), suggesting a diverse and complex mechanism by which alfalfa responded to Al stress. Furthermore, the qRT-PCR results showed a significantly positive correlation (R² = 0.9092) between the fold-changes of the gene expression ratios obtained by RNA-Seq and those obtained by qRT-PCR (Figure 7), indicating that our RNA-Seq experimental results are valid, and that RNA-Seq data is an accurate method to identify transcripts that are differentially regulated in response to Al. These results will greatly aid in our understanding of the molecular processes associated with Al stress responses and provide further insight for the Al-tolerant alfalfa breeding programs.

Previous studies have shown that some plant species can form sufficiently strong complexes with Al³⁺ to protect the roots from Al stress by releasing organic acids, such as citrate, oxalate, and malate (Kochian et al., 2015). In the present study, the KEGG pathway enrichment analysis of the DEGs indicated that some key genes related to citrate biosynthesis in the TCA cycle were also significantly enriched after Al treatment (Figure 6). These key genes included citrate synthase (CS), phosphoenolpyruvate carboxylase (PEPC), and malate dehydrogenase (MDH). Alterations in the activities of these enzymes may lead to accumulation of citrate in the cytoplasm (Ma et al., 2001). Overexpression of CS genes increases citrate efflux in cultured carrot cells (Koyama et al., 1999), Arabidopsis (Koyama et al., 2000), canola (Anoop et al., 2003), and tobacco (de la Fuente et al., 1997). When soybean roots are exposed to Al, the activities of mitochondrial MDH and CS have been found to increase in an Al-dose-dependent manner (Xu et al., 2010). In the present study, all of the key genes encoding enzymes involved in citrate biosynthesis mentioned above were identified, and the expression levels of the 4 PEPC, 3 MDH, and 3 AC genes were significantly increased during the Al treatment process (Figure S6A). These results suggested that organic acid production in response to Al occurs in the alfalfa roots.

Currently, two main types of Al resistance mechanisms have been documented that allow plants to cope with Al toxicity: one is Al exclusion mechanism, which prevents Al from entering the root apex (both apoplasm and symplasm) and the Al internal detoxification mechanism, in which Al enters the plant and is detoxified and sequestered (Kochian et al., 2015). Both mechanisms use organic acid anions. In the present transcriptome analysis, although the genes encoding TCA cycle enzymes were up-regulated, the related organic acid transporter genes, which are induced in other plants under Al stress, such as aluminum sensitive malate transporter (ALMT) and NRAMP aluminum transporter 1 (NRAT1), were not found. One gene belonging to the multidrug and toxin extrusion (MATE) family, which has been widely reported to function as citrate transporters in the induction of Al tolerance of many plants (Delhaize et al., 2012), was identified, but its transcript abundance was down-regulated (Figure S6B). These results further support the idea that there is no apparent correlation between the Al-induced expression of organic acid biosynthetic enzymes and increased exudation of organic acids (Chandran et al., 2008) and that the biosynthesis rather than exudation of organic acids is more critical for Al response in alfalfa roots. However, one gene and 8 genes that showing similarity to the aluminum-sensitive 1 (ALS1) and major facilitator superfamily (MFS) protein genes, respectively, were identified (Figure S6B). Previous studies have shown that these genes are involved in internal Al detoxification mechanisms by sequestering the Al-organic acid anion complexes inside the vacuoles of root cells (Huang et al., 2012; Yokosho et al., 2014). The up-regulation of these two types of genes in our RNA-Seq data suggested their involvement in the internal Al detoxification mechanism, although further functional analysis is required. In addition to the organic acid transporters, transporters for other small molecules/ions were also found to be up-regulated in alfalfa roots (Figure S6B), such as sugar transporters, sulfate transporters, vacuolar iron transporters, zinc transporters, and nitrate transporters, thus indicating that the uptake and translocation of other nutrients is affected by Al stress. Similar results have also been observed in buckwheat (Yokosho et al., 2014) and hydrangea (Chen et al., 2015), and additional investigations are necessary for further understanding this mechanism of Al responses in alfalfa.

Ribosomes are essential ribonucleoprotein complexes that are engaged in protein synthesis and thus are indispensable for metabolism, cell division, and growth (Wang J. Y. et al., 2013). In addition to their housekeeping functions, increasing evidence has suggested that ribosomal proteins also play more regulatory roles in leaf development, auxin responsiveness, wounding, and biotic/ abiotic stress responses (Liu et al., 2016). It has been reported that ribosomal protein genes can be differentially regulated by various environmental conditions. For example, up-regulation of several ribosomal proteins genes has been observed in maize exposed to UV-B light (Casati and Walbot, 2003), whereas the DEGs that respond to Al stress in the maize root are down-regulated (Maron et al., 2008). In Anthoxanthum, both up-regulated and down-regulated DEGs were identified in tolerant and sensitive genotypes when exposed to Al treatment (Gould et al., 2015). Microarray analyses have revealed that the transcript abundance for 5 ribosomal genes is increased after 6 h Al treatment, whereas after 48 h exposure, the transcripts for 5 ribosomal genes increase and the transcripts for 3 ribosomal genes decrease in abundance, thus indicating that there may be increased demand for specific ribosomal components during Al exposure (Kumari et al., 2008). In the present study, the KEGG pathway with the largest number of significantly enriched DEGs was “ribosome.” In contrast to previous studies in which fewer Al-response DEGs were detected, many more (157) DEGs were identified in our study, and most are components of large or small ribosomal subunits (Figure S6C). Among these DEGs, the transcript abundances for 93 DEGs were up-regulated during the 24 h treatment, and the remaining DEGs were down-regulated, which suggested a high biological importance for ribosomal genes in response to Al stress in alfalfa.

Heat Shock Proteins (HSPs) are important in stabilizing, folding, and degrading damaged proteins. In addition, numerous studies have determined that HSPs also function as molecular chaperones and protect plants from damage under stress conditions, such as Al- (Kumari et al., 2008; Duressa et al., 2011), NaCl- (Jiang and Deyholos, 2006), and heat temperature-induced stress (Sung et al., 2001). In the present study, there were 28 transcripts encoding HSPs (14 HSP70s, 7 HSP20s, 5 HSPs, and 2 HSP80s) that were enriched in the “protein processing in endoplasm reticulum” KEGG pathway (Figure S6D), which was the second most significantly enriched term. HSP70s is one of the most important members of the HSP protein family. The major role of HSP70 is to prevent protein aggregation and assist in protein refolding, import, translocation, signal transduction, and transcriptional activation (Wang Y. et al., 2013). In this transcriptome analysis, more than 53.8% of the HSPs were HSP70s, and 8 of these genes were up-regulated, whereas 6 were down-regulated, thus indicating that HSP70 may play important roles in Al stress responses in alfalfa.

Plant hormones are signaling molecules that regulate a wide range of metabolic and development processes at extremely low concentrations. Ethylene mediates root growth inhibition, whereas increased auxin biosynthesis in roots further facilitates ethylene-dependent growth inhibition (Swarup et al., 2007). Al has also been shown to affect root growth by modifying the levels of auxin (Ponce et al., 2005) and ethylene (Sun et al., 2007). It is well known that auxin response factors (ARFs) are a family of transcription factors that specifically bind to auxin-response elements of primary/early auxin response genes, whereas auxin-responsive proteins (AUX/IAA) repress the activity of ARFs, and thus this interaction plays a pivotal and concerted role in regulating the auxin response pathway (Dharmasiri and Estelle, 2004). Small auxin up RNA (SAUR) genes can be readily induced by exogenous auxins and function as negative regulators of auxin synthesis and transport (Kant et al., 2009). Ethylene responsive factors (ERFs) are plant-specific transcription factors belonging to the AP2/ERF family and have been reported to bind to the GCC-box found in the promoter of ethylene-inducible genes, thus acting as transcriptional activators of the ethylene response cascade and protecting cells from damage caused by metal stresses in plants (Makhloufi et al., 2014). In this transcriptome analysis, genes encoding enzymes involved in auxin homeostasis and response pathways, such as SAUR (CL3536.Contig2_All) and AUX/IAA (Unigene9889_All) were down-regulated, whereas the ARF (Unigene30886_All) was up-regulated through the Al treatment. In addition, the transcript abundances of most ERF genes were also increased (Figure S6E). These results suggested that both the inhibition of auxin production and the downstream signal transduction of these two hormones may be adopted in alfalfa in response to the root growth inhibition caused by Al stress.

ZL and WL conceived and designed the study. WL, LY, ZZ, LM, and YL performed the experiments. WL, CX and ZL analyzed the data. WL, CX and YW wrote the paper. All of the authors read and approved the final manuscript.

The raw sequence data supporting the results of this article are available in the Short Read Archive (SRA) (accession number SRS654614), http://www.ncbi.nlm.nih.gov/sra.