Comparative transcriptome analysis reveals gene network regulating cadmium uptake and translocation in peanut roots under iron deficiency

Iron (Fe) is an essential element for plant growth and development, whereas cadmium (Cd) is non-essential and highly toxic. Previous studies showed that Fe deficiency enhanced Cd uptake and accumulation in peanuts. However, the molecular mechanism underlying the increased Cd accumulation in Fe-deficient peanut plants is poorly understood. We employed a comparative transcriptome analysis approach to identify differentially expressed genes (DEGs) in peanut roots exposed to Fe-sufficient without Cd, Fe-deficient without Cd, Fe-sufficient with Cd and Fe-deficient with Cd. Compared with the control, Fe deficiency induced 465 up-regulated and 211 down-regulated DEGs, whereas the up- and down-regulated DEGs in Cd exposed plants were 329 and 189, respectively. Under Fe-deficient conditions, Cd exposure resulted in 907 up-regulated DEGs and 953 down-regulated DEGs. In the presence of Cd, Fe deficiency induced 1042 up-regulated and 847 down-regulated genes, respectively. Based on our array data, we found that metal transporter genes such as CAX4, COPT1, IRT1, NRAMP5, OPT3, YSL3, VIT3 and VIT4 might be involved in iron homeostasis. Moreover, combined with quantitative real-time PCR, IRT1, NRAMP3, NRAMP5, OPT3, YSL3, ABCC3, ZIP1, and ZIP5 were verified to be responsible for Cd uptake and translocation in peanut plants under iron deficiency. Additionally, a larger amount of ABC transporter genes was induced or suppressed by iron deficiency under Cd exposure, indicating that this family may play important roles in Fe/Cd uptake and transport. The up-regulated expression of NRAMP5 and IRT1 genes induced by iron deficiency may enhance Cd uptake in peanut roots. The decrease of Cd translocation from roots to shoots may be resulted from the down-regulation of ZIP1, ZIP5 and YSL3 under iron deficiency.


Background
Cadmium (Cd) is a non-essential and highly toxic heavy metal that is commonly released into the arable soil due to anthropogenic activities. The concentration of Cd in leaves greater than 5-10 μg g − 1 dry mass (DM) is toxic to non-tolerant crop plants [1]. For humans and animals, Cd may damage the mitochondrial and induce cell death by apoptosis and/or necrosis, leading to tissue inflammation and fibrosis [2]. Cd exposure has been shown to be associated with cancers of the prostate, lungs and testes in humans [3]. Due to its highly bioavailability, Cd is easily absorbed and accumulated in plants, and subsequently transferred to humans/animals via food chains. Thus, the presence of Cd in the arable soil can cause serious risks to human health. The minimization of Cd in edible parts of crops is an important demand especially in Cd-contaminated soil.
Iron (Fe) is an essential element that has multiple physiological functions including chlorophyll biosynthesis, photosynthesis, respiration, and redox reactions in plants. Despite its abundance in the earth's crust, Fe often precipitates into insoluble Fe(III) oxides under aerobic conditions, especially in high-pH and calcareous soils [4]. Consequently, dissolved Fe in the soil solution is usually far lower than that required to sustain plant growth [5]. Iron deficiency causes chlorosis, growth retardation, and reduced crop productivity. It has become an important yield-limiting factor for crops growing in calcareous and alkaline soils.
Peanut (Arachis hypogaea L.) is one of the world's fourth largest oilseed crop that is sensitive to Fe deficiency [6]. Peanut was also demonstrated to have a high capacity for accumulating Cd in both the seed and vegetative tissues depending on cultivars [7][8][9]. In previous studies, we have found that the uptake and accumulation of Cd in peanut plants is considerably enhanced by iron deficiency [10][11][12][13]. However, the mechanism underlying iron deficiency-induced increase of Cd accumulation in peanuts has not been well understood.
Generally, the accumulation of Cd in the shoot of plants is controlled by several biological processes, including (i) symplastic uptake by root epidermal cells, (ii) radial transport to the vascular cylinder, (iii) xylem loading, and (iv) root-to-shoot translocation [11]. Most of these processes are regulated by the transporters of several essential metals such as Fe, Mn, and Zn. During the past few years, several Fe transporters belonging to natural resistance associated macrophage proteins (Nramp), Zrt/Irt-like proteins (ZIP) and P1B-ATPases, have been identified to be responsible for the transport of Cd in plants [14]. Thomine et al. [15] found that iron deficiency could induce the expression of AtNramp3 and AtNramp4 in Arabidopsis. The increase of Cd uptake and accumulation under iron deficiency has been confirmed to be mediated by Fe 2+ transporters such as OsIRT1 and OsIRT2 in rice [16]. Similarly, our results showed that the expression of AhIRT1 and AhNramp1 genes in peanut roots can be induced by iron deficiency, which is associated with Cd uptake and accumulation [13]. He et al. [17] demonstrated that Fe supply prevents Cd uptake in Arabidopsis by inhibiting AtIRT1 expression and enhancing antagonism between Fe and Cd uptake.
Taken together, we hypothesize that iron deficiency-induced expression of Fe 2+ transporter genes may be responsible for the increased uptake and accumulation of Cd in iron-deficient peanut plants. To test this hypothesis, a comparative transcriptome analysis was carried out on iron-sufficient and -deficient peanut plants under Cd exposure. The aims were: (i) to obtain the gene expression patterns in the roots under iron deficiency; (ii) to identify the key genes related to Fe/Cd uptake and translocation; and (iii) to elucidate the gene regulatory network that are responsible for Cd uptake and translocation in peanut plants under iron deficiency. The results presented here will be useful for better understanding the mechanism of Cd accumulation induced by Fe deficiency and facilitate gene function studies in peanut.

Results
The accumulation and translocation of cd and Fe in plants Figure 1 shows plant growth as well as the accumulation and translocation of Cd and Fe in peanut plants. Compared with normal Fe supply, Fe-deficiency inhibited plant growth and resulted in leaf chlorosis ( Fig. 1 a and  b), while Cd did not change plant growth in both Fe treatments (Fig. 1b). It was also observed that Fe-deficiency significantly increased root Cd concentrations but decreased Cd concentration in xylem sap and the percentage of Cd in shoots ( Fig. 1 b and c), whereas Cd concentrations in the shoot were not affected (Fig.  1c). These results showed good agreement with previous findings [10][11][12][13], suggesting that Fe deficiency increased uptake but decreased root-to-shoot translocation of Cd in peanut plants. In comparison with Fe-deficient plants, Fe-sufficient plants exhibited lower root Fe concentration and higher proportion of Fe in shoots (Fig. 1d). The reverse trend in the percentage of Cd and Fe in shoots indicates that the two ions might compete with each other during the translocation from roots to shoots.

RNA sequencing analysis of eight cDNA libraries
To assess the global transcriptome profile of peanut roots in response to Fe-deficiency and/or Cd exposure, RNA-Seq analysis was performed on peanut roots exposed to Fe-sufficient without Cd (Fe 50 ), Fe-deficient without Cd (Fe 0 ), Fe-sufficient with Cd (Fe 50 Cd) and Fe-deficient with Cd (Fe 0 Cd). Two biological replicates were performed per treatment, with a total of eight cDNA libraries constructed. A total of 58,774,869, 54,456,644, 55,033,082, and 58,221,919 raw reads were produced from two biological replicate libraries of Fe 50 , Fe 0 , Fe 50 Cd, and Fe 0 Cd respectively. After removing low quality reads and reads containing adapter sequences, a total of 57,331,325, 53,124,279, 54,059,155, and 56,589,717 clean reads remained for Fe 50 , Fe 0 , Fe 50 Cd, and Fe 0 Cd, respectively. The percentage of clean reads in all eight libraries was more than 97.14%, and the percentage of Q20 exceeded 97.44% (Table 1). Pearson's correlation coefficients for all biological replicates were 0.949-0.959 (Additional file 1: Table S1), suggesting the data were highly reproducible.
The high-quality clean reads were mapped to the A. duranensis genome (reference genome) using HISAT (Hierarchical indexing for spliced alignment of transcripts). Ultimately, more than 79% of the clean reads were successfully mapped for all cDNA libraries, and over 77% were observed to be unique mapped reads ( Table 2), suggesting that the samples were comparable.

Identification of differentially expressed genes (DEGs)
A total of 63,191 genes including 34,553 known genes and 18,030 novel genes, 8017 lncRNA, 35 misc. RNA, 580 tRNA and 1976 pseudogenes were identified in eight cDNA libraries. Pairwise comparison analysis for each gene were performed between Fe-sufficient and Fe-deficient treatments (Fe 0 /Fe 50 and Fe 0 Cd/Fe 50 Cd), or between Cd-absent and Cd-present treatments (Fe 50 Cd/ Fe 50 and Fe 0 Cd/Fe 0 ). DEGs were identified by the threshold of P adj -value < 0.05. As a result, a total of 3024 genes were differentially regulated in the four comparisons, of which 676, 1889, 518, and 1860 DEGs were identified in Fe 0 vs Fe 50 , Fe 0 Cd vs Fe 50 Cd, Fe 50 Cd vs Fe 50 and Fe 0 Cd vs Fe 0 , respectively (Fig. 2a). Regardless of the absence or presence of Cd, 189 genes were differentially expressed between Fe-sufficient and Fe-deficient treatments, suggesting that these genes were specifically induced or supressed by Fe deficiency. It was also found that 190 genes were specifically regulated by Cd exposure under both Fe treatments (Fig. 2a).

Gene ontology (GO) analysis of DEGs
GO assignments were used to classify the functions of DEGs, and the results of significantly enriched GO terms (P adj -value < 0.05) were shown in Fig. 3. A total of 99 Fe deficiency-responsive DEGs (Fe 0 vs Fe 50 ) were assigned into 7 enriched GO terms consisting of 3 biological process (response to chemical, protein ubiquitination, protein modification by small protein conjugation) and 4 molecular function (nucleic acid binding transcription factor activity, transcription factor activity, ubiquitin-protein transferase activity, ubiquitin-like protein transferase activity) (Fig. 3a). Meanwhile, 66 Cd-responsive DEGs were assigned into 9 enriched GO terms, including 4 biological process (cell wall organization or biogenesis, external encapsulating structure organization, cell wall organization, cell wall modification), 2 cellular components (cell wall, external encapsulating structure) and 3 molecular function (carbohydrate binding, copper ion binding, pectinesterase activity) (Fig. 3b). A total of 857 DEGs between Fe 0 Cd and Fe 0 were assigned into 46 enriched GO terms consisting of 28 biological process (response to stress, response to oxidative stress, metal ion transport, response to chemical, etc.) and 18 molecular function (nucleic acid binding transcription factor activity, transcription factor activity, heme binding, tetrapyrrole binding, etc.) subcategories (Fig. 3c), and 832 DEGs of Fe 0 Cd vs Fe 50 Cd were assigned into 36 enriched GO terms including 17 biological process (response to stress, response to oxidative stress, response to chemical, etc.), 2 cellular components (cell wall, external encapsulating    Table S2). The four comparisons differed from each other in metabolic pathways of DEGs (Fig. 4). The top five pathways of Fe deficiency-responsive DEGs were phenylpropanoid biosynthesis, plant-pathogen interaction, diterpenoid biosynthesis, carotenoid biosynthesis, and amino sugar and nucleotide sugar metabolism (Fig. 4a). Among these pathways, phenylpropanoid biosynthesis was identified as significantly enriched (P adj -value < 0.05). For Cd-responsive DEGs, pentose and glucuronate interconversions, phenylpropanoid biosynthesis, glutathione metabolism, galactose metabolism, and diterpenoid biosynthesis were the top five categories, and no enriched pathway was identified (Fig. 4b). A larger number of enriched pathways were found in DEGs of Fe 0 Cd vs Fe 50 Cd (pentose and glucuronate interconversions, glutathione metabolism, flavonoid biosynthesis, and phenylpropanoid biosynthesis) (Fig. 4c) and Fe 0 Cd vs Fe 0 (phenylpropanoid biosynthesis, cysteine and methionine metabolism, diterpenoid biosynthesis, MAPK signaling pathway, pentose and glucuronate interconversions, phenylalanine metabolism, and ABC transporters) (Fig. 4d).

Verification of the DEG results
To further verify the transcriptome data, ten DEGs involved in metal transport, including IRT1, NRAMP3, NRAMP5, OPT3, YSL3, CAX4, HMA5, ABCC3, ZIP1, and ZIP5, were selected for RT-qPCR analysis. The RT-qPCR results, presented in Fig. 5, showed a good agreement with the RNA-Seq data (

Discussion
Although iron deficiency has been demonstrated to considerably enhance the uptake and accumulation of Cd in peanut plants [10][11][12][13], limited information is available about the physiological and molecular mechanisms underlying iron deficiency-induced increase of Cd accumulation in peanut. RNA-seq analysis has been used for revealing the molecular mechanisms of Cd uptake, translocation and accumulation in many plant species [18][19][20]. The complete genome sequencing of A. duranensis and A. ipaensis, the diploid ancestors of cultivated peanut [21], has shed light to genomic studies on cultivated peanut. However, sequence information of peanut in response to iron deficiency and Cd exposure is scarce.
In this study, we obtained 57,331,325, 53,124,279, 54,059,155, and 56,589,717 clean reads from the peanut roots treated with Fe 50 , Fe 0 , Fe 50 Cd, and Fe 0 Cd respectively. More than 79% of clean reads for all cDNA libraries were functionally annotated in the A. duranensis genome [21], and more than 77% were unique mapped reads (  (Fig. 2a). Of them, 54 DEGs were identified to highly similar with transporters that may be involved in the uptake and translocation of Fe/Cd in plants.
These results provide clues to mechanisms underpinning Cd uptake and accumulation in Fe-deficient peanut plants.
Genes involved in Fe uptake and translocation were greatly induced by Fe deficiency in peanut roots. As an Fe-efficient plant that develops strategy I mechanism in response to Fe deficiency, peanut takes up Fe 2+ by increasing Fe 2+ transporter coupled with an increase of ferric reductase activity and rhizosphere acidification by releasing protons from the roots under Fe-deficient conditions [6]. We found that both the Fe 2+ transport gene IRT1 and ferric reductase FRO1 were up-regulated by Fe deficiency in peanut roots (Table 3, Fig. 5). The result indicates that IRT1/FRO1 system constitutes the major pathway for Fe entry into root epidermal cells, and induction of these genes might improve Fe nutrition under Fe-deficient stress. A transporter gene of NRAMP family, NRAMP5 (NRAMP1) [13], was also highly induced in roots under Fe-deficient conditions ( Table 3, Fig. 5). The NRAMP5 (NRAMP1) have been shown to function in Fe transport in several plants such as Arabidopsis [22], Noccaea caerulescens [23], Malus baccata [24], peanut [25] and rice [26]. The induction of NRAMP5 (NRAMP1) under Fe-deficiency indicates that this gene is involved in the uptake of Fe by peanut roots. Interestingly, both the IRT1 and NRAMP5 are shown to be involved in Cd transport into root cells across membrane [13,16,27]. The current results confirmed that iron deficiency-induced expression of IRT1 and NRAMP5 may be responsible for the increase of Cd uptake and accumulation in Fe-deficient peanut plants [13].
We also found that large number of ABC transporter genes were significantly induced or suppressed by Fe deficiency and/or Cd exposure (Table 3). Although the detailed functions of ABC transporter genes is poorly understood, some members have been verified to play roles in the uptake and translocation of Fe/Cd in plants [28]. ABCB19, an auxin transporter that mediates long-distance polar auxin transport in stems and roots [29,30], was induced by Fe deficiency and/or Cd exposure (Table 3). Fe deficiency can increase the levels of auxin in the roots, which may promote the formation of root hair [31]. ABCC3 (MRP3), like most ABCC transporters such as such as ABCC2 and ABCC3, is vacuolar membrane-localized protein involved in the vacuolar transport of PC-Cd complexes [32]. The induction of ABCC3 in peanut roots by Fe deficiency with Cd exposure (Table 3, Fig. 5) may contribute to vacuolar Cd sequestration, enhancing Cd detoxification and reducing root-to-shoot translocation of Table 3 DEGs possibly involved in metal transport in Fe-sufficient and Fe-deficient peanut plants exposed to 0 or 2 μM CdCl 2 for seven days Cd. ABCC8 homologous to MRP6 in Arabidopsis, which was shown to be part of a cluster (AtMRP6, AtMRP3 and AtMRP7, as well as SAT3) involved in metal tolerance [33]. Besides ABC transporters, another peptide transporter OPT3 was significantly induced under Fe-deficient conditions (Table 3, Fig. 5). In Arabidopsis, OPT3 was shown to load iron into the phloem, facilitate iron recirculation from the xylem to the phloem, and regulate both shoot-to-root iron signaling and iron redistribution from mature to developing tissues [34,35]. Reduced expression of OPT3 induces an over accumulation of Fe in roots and leaves, partially due to an elevated expression of the IRT1 [36]. Therefore, the strong induction of OPT3 in the root of Fe-deficient peanut seedlings suggests that the gene play a role in the redistribution of Fe between vegetative tissues. Moreover, OPT3 was also found to have an impact on the uptake and translocation of Cd in Arabidopsis [34,35]. Mendoza-Cózatl et al. [35] observed that an OPT3 mutant of Arabidopsis, opt3-2, over-accumulates Cd in seeds and roots, but under-accumulates Cd in leaves. However, Zhai et al. [34] demonstrated that OPT3 is not involved in the phloem-based remobilization of Cd from source to sink tissues, despite it can mediate Cd transport in vitro. They speculated that increased Cd in seeds of the mutant is due to the impact of reduced expression of OPT3 on other transporters such as YSL1 or YSL3 [34]. YSL3 Table 3 DEGs possibly involved in metal transport in Fe-sufficient and Fe-deficient peanut plants exposed to 0 or 2 μM CdCl 2 for seven days (Continued) is a plasma-localized transporter delivering a broad range of nicotianamine-metal complexes. SnYSL3 gene is shown to be involved in the translocation and detoxification of Cd in Solanum nigrum, and its overexpression in Arabidopsis increased the translocation of Cd and Fe from roots to shoots [37]. In the present study, YSL3 was down-regulated in peanut roots under Fe-deficient stress ( Table 3, Fig. 5), which may be partially responsible for the reduced root-to-shoot translocation of Cd in Fe-deficient plants. Fig. 5 The qRT-PCR analysis of metal transport-related genes in peanut roots exposed to 0 or 2 μM Cd under Fe-sufficient (Fe50) and Fe-deficient (Fe0) conditions. The relative expression of each gene was calculated as the 2 −ΔΔCT value and normalized by geometric mean of three stably expressed reference genes. Data are presented as means with SE (n = 3). Different letters above error bars indicate that values are significantly different (P < 0.05) according to Duncan's multiple range test after one-way analysis of variance Several other transporter genes such as CAX4-like, COPT1, VIT3 and VIT4, were highly induced or suppressed by Fe deficiency (Table 3, Fig. 5). CAX4 is a cation/H + antiporter that plays a key role in mediating cations, such as Ca 2+ and Cd 2+ , influx into the vacuole [38]. Thus, the induction of CAX4-like by Fe deficiency in peanut roots may enhance vacuole sequestration of Cd and consequently, reducing root-to-shoot Cd translocation. COPT1 participates in copper acquisition and accumulation and regulates root elongation and pollen development [39]. In the current study, we found that COPT1 was up-regulated by Fe deficiency, but down-regulated in Fe-deficient seedlings exposed to Cd. We inferred that COPT1 indirectly effect Fe and/or Cd acquisition by regulating root elongation. Regardless of Cd exposure, the vacuolar iron transporters, VIT3 and VIT4, were suppressed under Fe deficiency condition. These alterations could reduce Fe sequestration in the vacuolar of root cells, and therefore enhance Fe transport from roots to shoots. NRAMP3 localizes at the vacuolar membrane, and is able to release Fe and Cd from the vacuolar under limited Fe conditions [40,41]. The induction of NRAMP3 may improve the remobilization of Fe from vacuoles of root cells in Fe-deficient peanut plants particularly under Cd stress.
Apart from abovementioned genes, three genes belonging to ZIP family, ZIP1, ZIP5, and MTP11, differentially expressed between the control and Cd exposure with Fe deficiency (Table 3, Fig. 5). ZIP1 shows a high degree of homology to the AtZIP2 from Arabidopsis [42] and OsZIP1 from rice [43], while ZIP5 is homologous to AtZIP1 [42]. Both the AtZIP2 and OsZIP1 show a broad substrate specificity for divalent cations including Cd [42,43]. In Arabidopsis, AtZIP2 is a plasma membrane localized transporter that may contribute to Mn/Zn movement in the stele to the xylem parenchyma, for subsequent xylem loading and transport to the shoot [44]. Similarly, OsZIP1 transcripts were localized to the epidermis and stele of roots of zinc-deprived plants, suggesting the involvement of this transporter in zinc absorption or transfer from the vascular tissue [43]. AtZIP1 is a vacuolar transporter that may play a role in remobilizing Mn from the vacuole to the cytoplasm in root stellar cells, and may contribute to radial movement to the xylem parenchyma [44]. Thus, the reduced expression of ZIP1 and ZIP5 in the roots of Fe-deficient peanut plants under Cd stress may inhibit root-to-shoot Cd translocation. In the case of MTP11, it was shown to function in Mn transport and tolerance by sequestering Mn into the pre-vacuolar compartments in Arabidopsis [45]. Whether MTP11 is associated with Fe and/or Cd transport need to be further investigated.

Conclusions
In conclusion, the current comparative study revealed that CAX4, COPT1, IRT1, NRAMP3, NRAMP5, OPT3, YSL3, VIT3 and VIT4 might be involved in iron homeostasis in Fe-deficient peanut plants. More importantly, some genes such as IRT1, NRAMP3, NRAMP5, OPT3, YSL3, ABCC3, ZIP1, and ZIP5, were identified to be responsible for the uptake, distribution, and translocation of Cd in peanut plants under iron deficiency. Based on our array data, we proposed a model to explain why iron deficiency induced an increase of Cd uptake but a decrease of Cd translocation from roots to shoots in peanut plants. The up-regulated expression of NRAMP5 and IRT1 genes induced by iron deficiency may enhance the uptake of cadmium by peanut roots. Iron deficiency-induced down-regulation of ZIP1, ZIP5 and YSL3 might result in a decrease of Cd xylem (or phloem) loading in root stele, and consequently, reducing root-to-shoot Cd translocation in peanut plants. Additionally, although detailed information is still unclear, larger amounts of ABC transporter genes were induced or suppressed by iron deficiency under Cd exposure, indicating that further study of this family would be helpful to understand the mechanism of Fe/Cd uptake and transport.

Plant growth and treatment
Peanut plants (Arachis hypogaea cv. Fenghua 1) were cultured as the conditions previously reported by Chen et al. [13] in a growth chamber of the Huaibei Normal University. Seeds (obtained commercially from the Peanut Institute of Shandong Province, Qingdao) were sterilized with 1% sodium hypochlorite for 10 min, and then they were rinsed with tap water for 24 h and germinated in well-washed sand. After 5 days of emergence, the uniform sized seedlings were selected and transferred to polyethylene pot (six plants per pot) containing 3.5 L of nutrient solution (pH 5.8) [11]. Seven-day-old seedlings were pretreated with (Fe 50 ) or without (Fe 0 ) 50 μM FeEDTA for 7 d in basal nutrient solution. Thereafter, 0 or 2 μM CdCl 2 were added in nutrient solution for each Fe treatment, and continuously cultured for one week. The Cd level (2 μM) was arranged according to the Cd concentration in the soil solution of a Cd-contaminated farmland in China (soil moisture at 60% of field capacity) [7]. The experiment was arranged as a completely random design with nine replications (pots). During the growing period, the nutrient solution was renewed twice a week, and pots were randomly arranged and moved daily to minimize position effects.
Root samples for RNA-seq (two biological replicates) and qRT-PCR (three biological replicates) analysis were collected separately from Fe-sufficient and Fe-deficient plants after Cd treatment. Each biological replicate contains a pool of six different plants growing in each pot. All samples were immediately frozen in liquid nitrogen and stored at − 80°C.

Determination of cd and Fe in plants
Three pots of plants for each Fe treatment under Cd exposure were used for collecting xylem sap excluded from the cut surface as the method described by Uraguchi et al. [46]. After weighing, the collected sap was used for determining the Cd concentration by graphite furnace atomic absorption spectrometry (GF-AAS). Roots and shoots were separated and rinsed three times with deionized water. Thereafter, plant samples were oven-dried at 70°C for constant weight. The dried tissues were weighed and ground into powder. The concentration of Cd and Fe in plant samples was measured by flame AAS after digested with HNO 3 -HClO 4 (3:1, v/v).
The translocation of Cd and Fe from roots to shoots was indicated as the percentage of metal in shoots, which was calculated as follows [7]: Total RNA was extracted by using Trizol® Reagent (Invitrogen) and purified using the RNeasy Plant Mini kit (Qiagen) according to the manufacturer's instructions. The purity and integrity of RNA were analyzed using NanoPhot-ometer® spectrophotometer (IMPLEN, CA, USA) and Agilent 2100 Bioanalyzer (Agilent, USA), respectively. Eight cDNA libraries named Fe 50 _1, Fe 50 _2, Fe 0 _1, Fe 0 _2, Fe 50 Cd _1, Fe 50 Cd _2, Fe 0 Cd_1 and Fe 0 Cd_2 were constructed as the method previously described by Hu et al. [47]. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq™ 2500 platform and 125 bp/150 bp paired-end reads were generated. All library construction and sequencing were done at the Novogene Bioinformatics Institute (Beijing, China).

Bioinformatics analysis Data filtering and mapping
To obtain high-quality clean reads, the raw reads from RNA-seq were filtered by removing the adaptor sequences, ambiguous 'N' nucleotides and low-quality sequences. Clean reads were aligned to the peanut reference genome (NCBI) using HISAT2 (v2.2.4, http:// www.ccb.jhu.edu/software/hisat) with default parameters. The mapped reads of each sample were assembled by StringTie (v1.3.3b) [48] in a reference-based approach.

Identification of DEGs
The numbers of reads mapped to the reference were counted using featureCounts v1.5.0-p3 [49]. The gene expression levels were represented by the expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM), which was calculated on the basis of the length of the gene and reads count mapped to this gene. Differential expression analysis of two groups was performed using the DESeq2 R package (1.16.1) according to the method described by Love et al. [50]. The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P adj -value < 0.05 found by DESeq2 were assigned as differentially expressed.

Gene ontology (GO) and pathway analysis of DEGs
Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented by the cluster-Profiler R package, in which gene length bias was corrected.
GO terms with corrected P value less than 0.05 were considered significantly enriched by differential expressed genes. KEGG pathways were retrieved from KEGG web server (http://www.genome.jp/kegg/). The clusterProfiler R package was used to test the statistical enrichment of differential expression genes in KEGG pathways.

qRT-PCR validation
Ten DEGs were randomly selected for qRT-PCR validation. Primer sequences of these genes as well as reference genes are listed in Additional file 3: Table S3. Total RNA (0.5 μg) from each root sample was reverse transcribed into cDNA using Prime Script® RT reagent Kit (Takara, Dalian, China) and random primers following manufacturer's instructions. Quantitative PCR reactions were performed in 20 μl reaction volumes using a SYBR Premix EX Taq Kit (Takara) according to the manufacturer's instructions. Reactions were carried out on an ABI7300 (Applied Biosystems, CA, USA). Each biological replicate was technically replicated three times. The relative expression levels of the selected genes were calculated using the 2 -ΔΔCT method and normalized by geometric mean of three stably expressed housekeeping genes (AhADH3, Ah60S and Ahactin) [51,52].