Integrative analysis of hexaploid wheat roots identifies signature components during iron starvation.

Iron (Fe) is an essential micronutrient for all organisms. In crop plants, Fe deficiency can decrease crop yield significantly; however, our current understanding of how major crops respond to Fe deficiency remains limited. Herein, the effect of Fe deprivation at both the transcriptomic and metabolic level in hexaploid wheat was investigated. Genome-wide gene expression reprogramming was observed in wheat roots subjected to Fe starvation, with a total of 5854 genes differentially expressed. Homoeologue and subgenome-specific analysis unveiled the induction-biased contribution from the A and B genomes. In general, the predominance of genes coding for nicotianamine synthase, yellow stripe-like transporters, metal transporters, ABC transporters, and zinc-induced facilitator-like protein was noted. Expression of genes related to the Strategy II mode of Fe uptake was also predominant. Our transcriptomic data were in agreement with the GC-MS analysis that showed the enhanced accumulation of various metabolites such as fumarate, malonate, succinate, and xylofuranose, which could be contributing to Fe mobilization. Interestingly, Fe starvation leads to a significant temporal increase of glutathione S-transferase at both the transcriptional level and enzymatic activity level, which indicates the involvement of glutathione in response to Fe stress in wheat roots. Taken together, our result provides new insight into the wheat response to Fe starvation at the molecular level and lays the foundation to design new strategies for the improvement of Fe nutrition in crops.

Introduction

Iron (Fe) is among the essential micronutrients in plants that participate as catalytic cofactors in several key processes including photosynthesis, respiration, chlorophyll biosynthesis, and nitrogen fixation (Kim and Rees, 1992; Morrissey and Guerinot, 2009; Li ). The bioavailability of Fe in the soil is strongly dependent on its solubility, aerobic and calcareous soil condition, pH, and the presence of natural ligands secreted by plant roots (Marschner, 1995; Morrissey and Guerinot, 2009; Thomine ). To circumvent the above challenges for Fe uptake by the roots, plants have adapted two modes of uptake strategies. Strategy I, mostly predominant in eudicot species, primarily relies on the enrichment of rhizospheric regions with protons (H+) and other reducing agents (Hell and Stephan, 2003; Santi ; Santi and Schmidt, 2009; Kobayashi and Nishizawa, 2012). In contrast, graminaceous species follow the Strategy II mode of uptake, which involves the transport of Fe in the complex-chelated form (Mori ; Kobayashi and Nishizawa, 2012; Connorton ). The primary components involved in chelation are the phytosiderophores (PSs) secreted by plant cells in the rhizospheric region mainly by efflux transporters (Morrissey and Guerinot, 2009; Kobayashi and Nishizawa, 2012). One of the main components of these secreted siderophores involved in Fe chelation are referred to as mugineic acids (MAs). The transporters involved in the secretion of the MAs are identified as transporter of mugineic (TOM) acid proteins (Nozoye ). The complex chelated form of Fe3+–PS is taken up by the specific root transporter referred to as yellow stripe-like transporter (YSL) proteins (Curie ; Gross ; Yordem ). Subsequently, Fe is transported in the plant organelles by multiple partners including specialized long-distance, tissue-specific transporters as reported in monocots (e.g. rice and maize) and eudicots (e.g. Arabidopsis) (Kim ; Waters and Sankaran, 2011). In other monocots such as maize and rice, the presence of genetic components for both Strategy I and Strategy II was reported (Ishimaru ; Inoue ; Lee ; Zanin ). At the metabolome level, analytical approaches utilizing GC-MC and LC-MS were used to study the components of Fe starvation in plants (Kabir ; Palmer ). However, the molecular and metabolic activity predominant for Fe uptake by the roots of hexaploid wheat still remains to be elucidated.

Microarrays have been successfully used to investigate the global transcriptional changes in Arabidopsis plants grown in Fe-starved conditions (Thimm ; Buckhout ; Yang ). Significant changes in the expression pattern of the genes primarily involved in the Strategy II mode of uptake were observed in roots of maize and rice subjected to Fe deficiency stress (Quinet ; Bashir ; Kobayashi ; Li ). The transcriptome analysis of Fe-starved plants is characterized by an important change in gene expression of several transcriptional regulators such as transcription factors (TFs) (Colangelo ; Kobayashi ; Long ; Li ; Connorton ) or key genes related to phytohormone homeostasis (Schmidt ; Hindt and Guerinot, 2012). Fe starvation also leads to changes in the abundance of transcripts related to plant metabolism and genes involved in signalling pathways that modulate nutrient uptake. Additionally, genes involved in ethylene/auxin signalling or linked with certain other macronutrients such as nitrogen, sulfur, and phosphorus are also significantly expressed (Romera and Alcantara, 1994, Zheng ; Romera ; Borlotti ; Zuchi ; Lin ; Zanin ; Garnica ). Therefore, transcriptome analysis is a powerful approach to help understand the network involved in plant response to Fe.

Hexaploid wheat is an important crop that is also a good, affordable source of nutrition. Fe deficiency strongly affects the growth of cereal crops and productivity (Yousfi ; Bocchini ). The development of transcriptome technology such as RNA sequencing (RNA-seq) and the availability of the genome sequence for hexaploid wheat (genomes A, B, and D) (International Wheat Genome Sequencing Consortium, 2014) combined with metabolomic approaches will help improve our knowledge of wheat response to nutritional stress (e.g. –Fe) (Borrill ). Our current knowledge on how the wheat genome responds to Fe-limiting conditions remains obscure and unanswered. In addition, how the different genomes could impact the homoeologue-based expression of the transcripts, categorically under different stress conditions, remains to be explored (Ramirez-Gonzalez ). In the current study, an RNA-seq-based approach was utilized to address the molecular components involved in wheat roots during Fe-starved conditions. Further, a genome-based expression study was performed that provided a preliminary clue on coordinated expression of the homoeologues. Metabolic profiling results support the data from the transcriptome study and pinpoint the role of the glutathione-mediated response. Collectively our results provide the first insight into molecular and biochemical responses of hexaploid wheat under Fe starvation.

Materials and methods

Plant materials, starvation conditions, and plant sampling

Hexaploid bread wheat variety ‘C-306’ was used for experiments. Wheat grains were subjected to stratification in the dark at 4 °C overnight and were allowed to germinate for 5 d on Petri plates containing 2–3 layers of wet Whatman filter paper. The endosperms were removed from the developing seedlings once they started browning. Subsequently, the seedlings were transferred in a PhytaBox™ and grown for 7 d in Hoagland’s nutrient solution. A total of 48 seedlings were divided equally in four PhytaBoxes, used for each treatment. After 7 d, nutrient solutions were replaced on the basis of different treatments. For Fe starvation (–Fe), 2 µM Fe (III) EDTA was used as the Fe source. For control plants (+Fe), concentrations of nutrients were unchanged in the above-mentioned Hoagland’s solution containing 20 µM Fe(III) EDTA. Treated seedlings were grown in the described medium for 20 d in a growth chamber set at 20±1 °C, 50–70% relative humidity, and a photon rate of 300 μmol quanta m−2 s−1 with a 16 h day/8 h night cycle. For sampling, roots and shoots were collected at different time points after starvation (5, 10, 15, and 20 d). On the basis of a distinct phenotype, samples collected at 20 days after starvation (DAS) were used for transcriptome analysis. A total of four biological replicates (each containing 10–12 seedlings) were used. Subsequently, RNA extractions were performed independently for each of the pools. Prior to RNA sequencing, quality of RNA was checked and extractions derived from two replicates were pooled together, thereby generating two experimental samples for each of the respective conditions. The remaining samples were snap-frozen in liquid nitrogen and stored at –80 °C. To clearly observe the primary root and the first-order lateral root, individual plants were moved onto a 150 mm wide Petri plate filled with distilled water and characteristics were manually examined. To ascertain root characteristics, 6–8 wheat seedlings were used for each of the treatments with two experimental replicates.

RNA extraction and Illumina sequencing

Total RNA was extracted from the treated root samples along with the control by Trizol (Invitrogen, ThermoFisher, USA) as per the manufacturer’s instruction. The quality and quantity were checked on a 1% denaturing RNA agarose gel and Nanodrop, respectively. Subsequently, all the RNAs used for library preparations were checked to ensure that their RNA integrity number was ≥8.5 using Bioanalyzer (Agilent, USA). The RNA samples that passed quality control were then processed for library preparation. The paired-ended libraries were prepared from the total RNA using the Illumina TruSeq stranded mRNA library prep kit as per the instructions (Illumina Inc., USA). The generated mean of the fragment size distribution of the libraries was 559, 584, 54, and 604 bp for the samples. The generated libraries were sequenced on the NextSeq 500 using 2×150 bp chemistry. The raw reads were processed further before the sorted high-quality reads were mapped to the reference genome.

RNA-seq analysis

Adaptor clipping and quality trimming of the raw reads were performed using Trimmomatic-0.35. The sequenced raw reads were processed to obtain high-quality clean reads. Ambiguous reads (reads with unknown nucleotides ‘N’ >5%) and low-quality sequences (>10% quality threshold QV <20 phred score) were removed. A minimum length of 100 nt after trimming was applied. Finally, high-quality (QV >20), paired-end reads were used for reference-based read mapping. The genome of Triticum aestivum L. was taken as reference for analysis. The Gene Feature Format files were downloaded from Ensembl Plants (TGACv1.37). The reads were mapped to the reference genome using TopHat v2.1.1 with default parameters. The Cufflinks v2.2.1 program assembled the transcriptome data from RNA-seq data and quantified their expression. Mapped reads were subjected to Cufflinks, followed by Cuffmerge and Cuffdiff. Log2 fold change (FC) values >1 were considered up-regulated whereas those with an FC <1 were considered as down-regulated. These genes were further categorized on the basis of statistical significance (P<0.05) and the false discovery rate (FDR 0.05) after Benjamin–Hochberg corrections for multiple testing for their significant expression.

Comprehensive gene annotation of wheat sequences was done using the KOBAS 3.0 (Xie ) annotate module by alignment with rice sequences (BLASTP, cut-off 1e-5). MapMan was used to visualize the pathways involving wheat differentially expressed genes (DEGs). Pathway enrichment analysis was performed using the KOBAS standalone tool. MeV was employed to construct heatmaps for selected DEGs using the normalized expression values of genes. The data generated from this study have been deposited in the NCBI Sequence Read Archive (SRA) database and are accessible with the submission ID-SUB5206887 and BioProjectID-PRJNA529036.

Gene annotation filtering and functional enrichment analysis

Significant sets of DEGs under Fe starvation were further mapped using Gene Ontology (GO) and MapMan. The GO annotation was downloaded from ensembl plants (https://plants.ensembl.org/biomart/martview). Mercator (Lohse ) was used to build a MapMan mapping file for TGACv1 sequences, and DEGs were visualized with the MapMan v3.1.1 tool (Thimm ). For functional categorization of DEGs that were positively and negatively correlated with Fe starvation, BINGO version 3.0.3 (Maere ) was used to perform GO enrichment analysis with a hypergeometric test, and significant terms with an FDR value <0.05 were considered. For GO mapping, the GO_full.obo ontology file was downloaded from the GO consortium. MapMan classification was used to categorize DEGs into TFs. The results were visualized as a network using EnrichmentMap version 2.2.1 (Merico ), and gene expression overviews in various pathways were visualized with the MapMan tool.

Homoeologue-specific expression analysis for genome bias

To identify wheat homoeologous triplets; ensembl biomart TGACv1 was used to extract all possible homoeologous relationships. This led to 86 830 pairwise homoeologous relationships. An in-house script was used to select only the triplets where the contribution from each genome was 1:1:1 (A:B:D), thus generating 16 850 triplets. Further, triplets resulting from potential translocation events were not considered; that is, only homoeologue triplets from the same chromosome (e.g. the 1A, 1B, 1D triplet is accepted, whereas, 2A, 3B, 3D is rejected) were taken for analysis. Finally, 15 604 triplets (15 604×3=46 812 genes) were used for studying genome induction bias in response to Fe stress. Paired-end reads were aligned to the reference (selected scaffolds from the genome that harbour 15 604 triplets) using TopHat v2.1.1 with a specific argument (--b2-very-sensitive) (Powell ), which leads to more stringent alignments as required for homoeologues. The Cufflinks pipeline was used to obtain FPKM (fragments per kilobase of transcript per million mapped read) values and DEGs.

Further, relative abundance levels and expression bias for homoeologues were studied by considering the homoeologous gene triplets with expression of FPKM ≥1 each in both the control and Fe-starved conditions. For this, the normalized relative expression for each homoeologue within a triad was calculated. For example, the relative expression from A will be represented as:

Thus, the relative expression levels of A, B, and D homoeologues within each triad were calculated similarly. Seven homoeologue expression bias categories were defined as described earlier by Ramirez-Gonzalez . In total, seven categories defined as one balanced category and six unbalanced homoeologue-suppressed/homoeologue-dominant (from either of the genomes) were listed for the ideal relative expression of A, B, and D. Eucledian distance (using the cdist function from the rdist package, R3.3.2) was calculated between normalized relative expression for each triad and the seven ideal categories. The shortest distance was used as a deciding factor to group the triads into the seven respective categories.

Quantitative real-time PCR (qRT-PCR)

For validation of the gene expression, qRT-PCR analysis was performed. Total RNA (2 µg) isolated from the above experiments was used to validate the expression by using qRT-PCR. Genomic DNA present in trace amounts was removed by DNase I treatment using a Turbo DNA-free kit (Invitrogen, ThermoFisher, USA). Further, cDNA was synthesized from 2 µg of DNA-free RNA using Superscript III first strand (Invitrogen, ThermoFisher, USA) with random hexamer primers following the manufacturer’s guidelines. For qPCRs, gene-specific primers of each gene were designed from conserved regions of all three homeologue sequences (Supplementary Table S1 at JXB online). qRT-PCR was performed using QuantiTect SYBR Green RT-PCR mastermix (Qiagen, USA) with programs recommended by the manufacturer in the ABI 7700 sequence detector (Applied Biosystems, USA). ADP-ribosylation factor (ARF) and actin were used as internal controls. Two independent experimental replicates with four technical replicates were performed for each sample. The relative amount of gene expression was calculated using the 2−ΔΔCT method.

Metabolite extraction and GC-MS profiling

Extraction of total metabolites was performed similarly to as previously described (Wang ). Roots of plants grown under –Fe and +Fe were sampled in triplicate and dried for 1 week. A 50 mg aliquot of each crushed sample was extracted in a mixture of solvents (pre-cooled; 300 µl of methanol, 30 µl of 2 mg ml−1 nonadecanoic acid methylester, and 30 µl of 0.2 mg ml−1 sorbitol) for 15 min (70 °C, 1000 rpm). Further, at room temperature, 200 µl of chloroform was added and shaken for 5 min (37 °C, 1000 rpm). To obtain phase separation, 400 µl of H2O was added to each sample, vortexed, and centrifuged (10 min, 13 000 rpm). From the upper polar phase, ~200 µl was finally aliquoted for complete drying.

For GC-MS analysis, metabolites were subjected to methoxyamination and trimethylsilylation. Dried polar phases were shaken for 1.5 h at 30 °C in 40 µl of MeOX (40 mg ml−1 methoxyaminhydrochloride in pyridine) followed by 30 min shaking at 37 °C in 80 µl of BSTFA mixture [70 µl of N,O-bis(trimethylsilyl)trifluoroacetamide+10 µl of alkane mix]. The derivatized metabolites were subjected to GC-MS analysis (Agilent Technologies 7890, USA) coupled with MS. Measurement from an injection volume of 1 µl was taken in split-less mode in a DB-5 column (30 m×0.25 mm, 0.25 µm film thickness, Agilent) using helium as carrier gas. Metabolites were separated as described by Wagner . Qualitative analysis of chromatograms was performed on the MassHunter Qualitative analysis Sp1 workstation (Agilent, USA). Identification and annotation of each compound was supervised manually using AMDIS software and the NIST08 database (https://www.nist.gov/srd/nist-standard-reference-database-1a-v17). Data were normalized to sample weight and internal control (sorbitol). Statistical analysis was performed as described previously (Quanbeck ). The log2 ratio of metabolite abundances in –Fe was plotted against +Fe. Delta method approximation was used to calculate the SEs of the log ratio, SE log ratio=1/ln 2√[(SET/T)2+(SEC /C)2], where SET and SEC are standard errors of average –Fe and +Fe metabolite abundances.

Measurement of glutathione S-transferase (GST) activity and Fe mobilization assay for PS estimation

Activity measurement of GST was performed in the wheat roots subjected to –Fe for 10, 15, and 20 d of treatment along with the control plants (no stress) as mentioned in ‘Plant materials, starvation conditions, and plant sampling’. The estimation was done using the glutathione S-transferases assay kit (Sigma, USA). Briefly, 100 mg of tissue was used for total protein extraction. An equal amount of total protein (25 µg) was used as a source of enzyme, and 1-chloro-2,4-dinitrobenzene (DNB) was used as a substrate. The resulting GST–DNB conjugate was measured at 340 nm wavelength during the time course of the reaction. The direct increase in absorption was measured and GST activity was calculated as described in the manufacturer’s instruction kit.

For PS release, Fe mobilization assay was performed as described earlier (Takagi ). Briefly, 10 seedlings (each for –Fe and +Fe for 10 d treatment) were used for the experiment. After treatment, the seedlings were subjected to PS release 2 h after onset of the light period in aerobic conditions for 3 h in 20 ml of deionized water with 10 mg l−1 Micropur (Katadyn, Switzerland). Next, 2 ml of freshly precipitated solution of Fe(OH)3 and 0.5 ml of 0.5 M Na acetate buffer (pH 5.6) was added to 8 ml of the collection solution. The solution was shaken for 2 h and subsequently was filtered (Whatman #1) into 0.2 ml of 6 N HCl. Ferric iron was reduced by addition of 0.5 ml of 8% hydroxylamine-hydrochloride and heating to 60 °C for 20 min. The total concentration of ferrous iron was calculated by measuring absorbance at 562 nm after adding 0.2 ml of 0.25% ferrozine and 1 ml of 2 M Na-acetate buffer (pH 4.7).

Elemental analysis using inductively coupled plasma-mass spectroscopy (ICP-MS) and nitrate estimation

Elemental analysis in roots and shoots was performed using ICP-MS. Metal analysis was performed as described earlier (Bhati ; Aggarwal, 2018). Briefly, the samples were ground to a fine powder and subsequently subjected to microwave digestion with HNO3 as described earlier (SuraPure™, Merck). The respective metal standards were also prepared for analysis. Three independent replicates from the experiments were used for metal analysis.

Nitrate content in wheat roots was measured according to a method described previously (Cataldo ). Briefly, 1 g of fresh tissue was homogenized in 6 ml of deionized water and centrifuged at 30 000 g for 15 min. Then 100 µl of supernatant was added to 400 µl of salicylic acid (w/v dissolved in concentrated H2SO4). After mixing, the reaction was kept at room temperature for 20 min. NaOH at 2 N (9.5 ml) was then added slowly to raise the pH above 12. The samples were allowed to cool and readings were taken at 410 nm in a spectrophotometer.

Results

Fe starvation affects wheat growth capacity and nutrient uptake

Fe starvation is known to affect plant growth capacity. In order to determine the effect of Fe starvation, 1-week-old wheat seedlings, grown on complete medium (in the presence of Fe), were transferred to Fe starvation medium for further growth. At 20 DAS, the wheat seedlings started showing strong phenotypic symptoms including visible chlorosis, and therefore a detailed study was performed for this specific time point. In response to Fe starvation, plants showed a decrease in the shoot biomass with an enhanced chlorosis phenotype and shortening of the root system compared with control wheat seedlings (Fig. 1A, B). Roots of Fe-starved wheat seedlings showed a decrease in the number of lateral roots and a significant reduction in the primary root length in comparison with control plants (Fig. 1C, D). Earlier studies have suggested that changes in the root system and Fe supply not only affect the Fe accumulation capacity but also impact the uptake of other nutrients such as Zn, Cd, etc. (Sperotto ; Shukla ). Therefore, the effect of Fe starvation stress on the uptake of Zn, Mn, Cu, and Mg in wheat undergoing Fe stress was studied (Table 1). Out data indicated increased uptake of nutrient elements such as Zn, Mn, Cu, and Mg in roots, but accumulation in shoots was either unaltered or decreased. This further supports the importance of studying the root response to Fe starvation to gain insight into the molecular mechanism evolved by wheat to cope with nutritional stress.

Fig. 1.

Effect of Fe starvation (–Fe) on the growth parameters of wheat seedlings post 20 DAS. (A) Phenotype of wheat seedlings exposed to Fe starvation. (B) Total biomass of roots and shoots of wheat seedlings after 20 DAS. A total of 12–15 seedlings were collected for calculating the fresh tissue weight (in g). (C) Number of first-order lateral roots in roots subjected to –Fe conditions and control plants (+Fe). (D) Primary root length of wheat roots. A total of 10–12 seedlings were used for measuring the total primary root length of wheat seedlings under –Fe and +Fe condition. #indicates a significant difference at P<0.05. (This figure is available in colour at JXB online.)

Table 1.

Metal concentration (mg g–1 DW of tissue) in roots and shoots of wheat seedlings subjected to –Fe stress

	Roots					Shoots
Treatments	Fe	Zn	Mn	Mg	Cu	Fe	Zn	Mn	Mg	Cu
+Fe	106.73±14.2	14.6±3.2	15.8±2.0	542±70	7.4±1.1	72.9±1.7	26.8±.09	17.4±.05	1031±2.3	10.4±0.6
–Fe	66.8±17	40.4±12	38.2±8.9	912±213	17.9±3.7	28.6±1.6	22.3±1.05	22.1±1.9	1013±48	5.1±0.3

Differential expression analysis, homoeologue (A, B, and D) induction, and expression bias during –Fe response

The effect of Fe starvation on the wheat root transcriptome has not been investigated to date. To perform this study, RNA-seq technology was used to identify the changes in the transcripts of wheat roots, where plants were grown in the presence or absence of optimal Fe (20 DAS). Our analysis resulted in 87 million quality filtered reads, with an average of nearly 22 million reads from each sample (>87% of reads had a quality score >30). Filtered reads from the four libraries had a mapping rate ranging from 81.7% to 85.4% when mapped against release-37 of the wheat genome using TopHat (Supplementary Table S2). As a quality check, a strong correlation within the two biological replicates from each condition was observed, while a clear variation was seen between the two conditions (Fig. 2A). We thereafter analysed the expression values as FPKMs, calculated by using Cufflinks software. DEGs were then identified by calculating logFC and performing statistical tests between FPKM values from control and stressed samples using CuffDiff. A total of 50 610 genes with an FPKM of ≥1 in at least one of the two conditions were considered to be ‘expressed transcripts’. In all, 7221 genes had logFC >0 and 8010 had logFC <0 (Fig. 2B). On setting up a criterion of logFC >1 for up-regulated genes, and of < –1 for down-regulated genes, and an FDR of <0.05, a total of 3478 genes were highly expressed, whereas 2376 were down-regulated under –Fe conditions in wheat roots (Fig. 2C). Interestingly, 45 genes were also induced exclusively during starvation conditions when compared with the control samples (Fig. 2C; Supplementary Table S3).

Fig. 2.

Analysis of RNA-seq data from the wheat roots during Fe starvation. (A) Principal component analysis of samples from control (Control_0, Control_1) and –Fe (Single1_0, Single1_1) conditions. (B) Volcano plot of DEGs; the x-axis shows the log2 fold change difference in the expression of genes in iron-starved conditions with respect to control, and the y-axis indicates the negative log of the P-value (pval) for the differences in expression. Genes without significant differences are indicated by grey dots. Significant genes with logFC>0 are represented by red dots, and those with logFC<0 are represented by green dots in the scatter plot. (C) Number of DEGs under –Fe. Up, up-regulated under –Fe (logFC>1); down, down-regulated under –Fe (logFC< –1); Exc–Fe, exclusively expressed in –Fe; ExcCtrl, exclusively expressed in control conditions.

Our data allowed us to analyse the chromosomal distribution of the DEGs under the –Fe condition. While all chromosomes contributed to the DEGs, the highest number of genes was mapped on chromosome 2 of the wheat genome (Fig. 3A). Equal representation of transcripts was observed for chromosomes 7 and 5. The remaining chromosomes, 1, 3, 4, and 6, showed 13, 15, 12, and 10% distribution of DEGs, respectively. In polyploidy crops such as wheat, homoeologue induction bias could impact plant response to various stresses (Liu ; Powell ). To determine the extent of induction bias (from A, B, and D subgenomes) in wheat during Fe starvation, homoeologue-specific expression analysis was performed. Starting from a list of 8473 gene triplets present/expressed selected as ‘accepted triplets’ (Supplementary Table S4), most of the homoeologue triplets (8349) showed no significant bias in expression (A=B=D). Of these, 8321 triplets appeared to be unaffected by –Fe stress, while 22 and six homoeologous triplets were up- and down-regulated, respectively. Homoeologue expression bias was observed in 124 homoeologue triplets. Eighty-seven triplets had only one of the homoeologues differentially expressed (up- or down-regulated). Subgenome-specific contribution towards this bias within triplets is depicted in the left panel of Fig. 3B. These include 47 in the category ‘1UP’ with A>B=D, B>A=D, and D>A=B, and40 in the category named ‘1DOWN’ having AFe starvation. A few of the prominent transcripts exclusively induced by these two genomes include transcripts related to MYB TFs, metal transporters, zinc transporters, RINGH-H2 type proteins, genes belonging to major facilitator superfamily proteins, etc. Additionally, 37 triplets had two of the homoeologues differentially expressed while the third showed normal expression even under Fe stress (Fig. 3B; right panel). The ‘2UP’ category comprises AB>D, AD>B, and BD>A, while ‘2DOWN’ comprises ABFe starvation, the additive homoeologue contribution from either the A or B subgenome was the highest.

Fig. 3.

Genomic distribution and homoeologue bias studies during Fe starvation. (A) Chromosomal distribution of DEGs. (B) Pie charts showing (left panel) subgenomic distribution of genome induction bias in triads where one of the homoeologues was up-/down-regulated. A, B, and D depict the subgenome to which the differentially expressed homoeologue belongs; right panel: distribution of triads for which two of the homoeologues were differentially expressed and the third one had normal expression. AB refers to the triads for which up-/down-regulation was observed in the homoeologues belonging to A and B subgenomes, while the D subgenome homoeologue behaved as normal with respect to control. (C) Sankey diagram depicting the homoeologue expression bias in control conditions and Fe starvation. Homoeologue triads were classified into the seven defined categories based on relative normalized expression within each triad. Nodes flowing from control to –Fe (Fe starvation) represent the triads with the same (flow to same category) as well as changed (flow to a different category) expression patterns across both conditions. Distinct colours represent the flow of triads belonging to the seven categories from the control condition into the same category under –Fe or transition into a different category under –Fe.

Table 2.

List of genes showing genome induction bias contribution from either the A or B subgenome

Gene	LogFC	RAP-DB description
A genome
TRIAE_CS42_3AL_TGACv1_197291_AA0665720	3.01	Myb transcription factor domain containing protein.
TRIAE_CS42_2AL_TGACv1_094922_AA0304900	2.8	Similar to Prolyl endopeptidase (EC 3.4.21.26) (Post-proline cleaving enzyme) (PE).
TRIAE_CS42_5AL_TGACv1_379416_AA1256390	2.44	Similar to Solute carrier family 35, member F1.
TRIAE_CS42_2AL_TGACv1_094606_AA0300340	2.18	RmlC-like jelly roll fold domain containing protein.
TRIAE_CS42_3AS_TGACv1_211026_AA0683370	3.36	Protein of unknown function DUF1399 family protein.
TRIAE_CS42_7AL_TGACv1_557470_AA1781720	3.47	Heavy metal-transporting P1B-ATPase, Root-to-shoot cadmium (Cd) translocation
TRIAE_CS42_6AL_TGACv1_471682_AA1512850	2.26	Similar to CRT/DRE binding factor 1.
TRIAE_CS42_6AS_TGACv1_485332_AA1543320	1.85	Zinc/iron permease family protein.
TRIAE_CS42_2AL_TGACv1_092944_AA0267670	2.43	Heavy metal transport/detoxification protein domain containing protein.
TRIAE_CS42_5AL_TGACv1_375378_AA1220660	2.21	Similar to SUSIBA2 (WRKY protein).
TRIAE_CS42_1AL_TGACv1_000708_AA0017440	2.57	Chlorophyll a-b binding protein 2, chloroplast precursor (LHCII type I CAB-2) (LHCP).
TRIAE_CS42_2AL_TGACv1_093154_AA0273510	1.86	C2 domain containing protein.
TRIAE_CS42_5AL_TGACv1_374025_AA1188170	2.66	Nodulin-like domain containing protein.
TRIAE_CS42_3AL_TGACv1_197522_AA0666570	2.22	Lipase, GDSL domain containing protein.
TRIAE_CS42_7AL_TGACv1_556100_AA1755640	2.34	Bifunctional inhibitor/plant lipid transfer protein/seed storage domain containing protein.
TRIAE_CS42_4AS_TGACv1_306527_AA1009640	2.2	Multi antimicrobial extrusion protein MatE family protein.
TRIAE_CS42_7AS_TGACv1_570336_AA1834060	1.98	Major facilitator superfamily protein.
TRIAE_CS42_2AL_TGACv1_093456_AA0280470	2.15	Glycolipid transfer protein domain domain containing protein.
TRIAE_CS42_1AL_TGACv1_002205_AA0039730	–2.09	Similar to IAA8 (Fragment).
TRIAE_CS42_4AL_TGACv1_290815_AA0989640	–2	Lipase, class 3 family protein.
TRIAE_CS42_5AL_TGACv1_374561_AA1203290	–2.24	Cinnamyl alcohol dehydrogenase (EC 1.1.1.195).
TRIAE_CS42_5AL_TGACv1_374359_AA1197930	–2.36	Similar to Lipoxygenase L-2 (EC 1.13.11.12).
TRIAE_CS42_6AL_TGACv1_472625_AA1524260	–1.89	RAG1-activating protein 1 homologue domain containing protein.
TRIAE_CS42_7AL_TGACv1_558250_AA1791610	–2.32	Similar to Pleiotropic drug resistance protein 3.
TRIAE_CS42_2AL_TGACv1_097246_AA0323500	–1.72	Similar to Peroxidase (EC 1.11.1.7).
TRIAE_CS42_7AL_TGACv1_556210_AA1758330	–2.21	Similar to Kaurene synthase A (Fragment).
TRIAE_CS42_1AL_TGACv1_000555_AA0014640	–1.84	No apical meristem (NAM) protein domain containing protein.
TRIAE_CS42_6AL_TGACv1_471077_AA1502240	–1.79	Similar to OSIGBa0145M07.8 protein.
TRIAE_CS42_7AL_TGACv1_559906_AA1801190	–1.68	Similar to H0801D08.12 protein.
TRIAE_CS42_7AL_TGACv1_558101_AA1790150	–1.76	Similar to Acyl-ACP thioesterase (Fragment).
TRIAE_CS42_6AL_TGACv1_472321_AA1520860	–2.06	Similar to Subtilisin-like protease (Fragment).
B genome
TRIAE_CS42_5BS_TGACv1_423346_AA1374840	2.19	Similar to Calmodulin NtCaM13.
TRIAE_CS42_3B_TGACv1_224030_AA0791180	2.2	Similar to IN2-2 protein.
TRIAE_CS42_2BL_TGACv1_129296_AA0377070	2.46	Similar to OSIGBa0127A14.7 protein.
TRIAE_CS42_6BL_TGACv1_499646_AA1588130	2.32	TGF-beta receptor, type I/II extracellular region family protein.
TRIAE_CS42_7BS_TGACv1_592587_AA1940840	2.25	Similar to RING-H2 finger protein ATL1R (RING-H2 finger protein ATL8).
TRIAE_CS42_2BS_TGACv1_146290_AA0461590	2.64	NA
TRIAE_CS42_5BL_TGACv1_404610_AA1306090	1.84	Similar to Senescence-associated protein SAG102.
TRIAE_CS42_5BL_TGACv1_405319_AA1324840	2.77	Similar to Transporter associated with antigen processing-like protein.
TRIAE_CS42_7BL_TGACv1_577614_AA1879540	2.69	Peptidase A1 domain containing protein.
TRIAE_CS42_4BS_TGACv1_329166_AA1098520	2.32	Similar to Alcohol dehydrogenase.
TRIAE_CS42_4BL_TGACv1_321683_AA1064050	2.56	Protein of unknown function DUF1262 family protein.
TRIAE_CS42_7BL_TGACv1_591489_AA1920550	NA	Similar to zinc transporter 4.
TRIAE_CS42_7BL_TGACv1_577301_AA1871590	2.41	Delta-tonoplast intrinsic protein.
TRIAE_CS42_4BS_TGACv1_329309_AA1100040	NA	Similar to Major facilitator superfamily antiporter.
TRIAE_CS42_2BL_TGACv1_129348_AA0379680	2	Hypothetical conserved gene.
TRIAE_CS42_2BS_TGACv1_148847_AA0495340	4.34	Helix-loop-helix DNA-binding domain containing protein.
TRIAE_CS42_2BL_TGACv1_130820_AA0418390	1.88	NA
TRIAE_CS42_1BL_TGACv1_031794_AA0120680	–2.79	Divalent ion symporter domain containing protein.
TRIAE_CS42_5BL_TGACv1_406838_AA1350360	–2.29	Similar to anther-specific proline-rich protein APG.
TRIAE_CS42_2BL_TGACv1_130182_AA0405730	–2.09	Protein of unknown function DUF3741 domain containing protein.
TRIAE_CS42_4BS_TGACv1_330685_AA1108490	–2.83	Serine/threonine protein kinase-related domain containing protein.
TRIAE_CS42_2BS_TGACv1_147909_AA0488910	–2.28	Delayed-early response protein/equilibrative nucleoside transporter family protein.
TRIAE_CS42_3B_TGACv1_224656_AA0799480	–1.83	Glyoxalase/bleomycin resistance protein/dioxygenase domain containing protein.
TRIAE_CS42_7BL_TGACv1_577476_AA1876230	–2.93	Protein kinase, catalytic domain domain containing protein.
TRIAE_CS42_6BL_TGACv1_499809_AA1592210	–1.88	Cupredoxin domain containing protein.
TRIAE_CS42_7BL_TGACv1_576755_AA1853830	–2.04	Cellulase (EC 3.2.1.4).
TRIAE_CS42_2BL_TGACv1_130686_AA0416110	–2.3	Similar to OSIGBa0096P03.3 protein.
TRIAE_CS42_2BS_TGACv1_146035_AA0453750	NA	Conserved hypothetical protein.
TRIAE_CS42_1BL_TGACv1_030699_AA0098220	–1.91	Lipase, class 3 family protein.
TRIAE_CS42_2BL_TGACv1_129634_AA0391150	–1.78	2OG-Fe(II) oxygenase domain containing protein.

Expression changes within triads with respect to iron stress were observed to identify the presence of differential regulation among homoeologues from each triad, thus contributing towards genome induction bias. The table gives a descriptive list of homoeologues that were up-/down-regulated under iron starvation only in the A subgenome, followed by those only in the B subgenome

For polyploid genomes such as wheat, the interaction of its subgenomes is known to affect the final phenotype or contribute towards trait development (Borrill ). To check the homoeologue/subgenome expression bias under –Fe conditions, the effect on the expression of transcripts from the genomes was determined by comparing the relative normalized expression for each homoeologue within a triad. This resulted in seven combinations including one balanced category and six homoeologue-specific dominance or suppression categories. Our analysis reflected that most of the triads fell under the balanced category (Table 3; Fig. 3C). This category was represented by 77% and 77.89% of the total triads for control and Fe starvation conditions, respectively. Triads with unbalanced expression varied in the range of 0.94 to 7.09 for all six subcategories across the two conditions (Table 3; Fig. 3C). Our analysis revealed that maximum genome-specific expression bias was observed for A and B in both the conditions. Interestingly, the D genome was least suppressed, with a representation of 5.15% of the total triads taken into consideration in control when compared with 4.93% in Fe starvation. A total of 89% triads showed a conserved balanced/unbalanced contribution across both conditions. Overall, a higher relative abundance of the D genome (control, 33.94%; –Fe, 33.95%) was noted as compared with the A (control, 33.06%; –Fe, 33.00%) and B (control, 33.04%; –Fe, 33.95%) genomes.

Table 3.

Percentage of homoeologue triads categorized into ideal genome expression bias categories in control and Fe-starved conditions

	Control	–Fe
Balanced	77.07%	77.89%
A-suppressed	6.90%	6.90%
B-suppressed	7.09%	6.82%
D-suppressed	5.15%	4.93%
A-dominant	1.12%	0.94%
B-dominant	1.24%	1.18%
D-dominant	1.43%	1.32%

Identification of DEGs pinpoints the prolific expression of genes involved in the Strategy II mode of Fe uptake

To identify the transcripts differentially expressed in roots in response to Fe starvation, the top 50 up-regulated and the top 50 down-regulated genes were shortlisted (Fig. 4). Expression of most of the highly up-regulated transcripts ranged from 12 to 4.8 log2FC (Supplementary Table S3), indicating their higher fold accumulation under Fe starvation compared with the control conditions. As observed in the top 50 up-regulated genes, the predominance of pathway genes encoding the components for Strategy II mode Fe uptake was observed across all the 3478 up-regulated genes. Categorically, these highly induced genes are known to be involved during the Strategy II mode of Fe uptake, including the subfamily of nicotinamine synthase (NAS) and deoxymugineic acid (DMA) biosynthesis genes (Fig. 4, left panel). Other important genes include genes belonging to the major facilitator superfamily including an ABC transporter, zinc-induced facilitator-like (ZIFL) transporters, and sulfate transporters. Therefore, under Fe starvation, the DMA biosynthesis genes were highly induced. Similarly, YSL genes were significantly induced under Fe starvation conditions. All the homoeologues of TaYSL9 and TaYSL1A showed especially high expression under starvation conditions (Supplementary Table S5). Along similar lines, genes coding for NRAMP also showed high transcript abundance in roots subjected to Fe starvation (Supplementary Table S5). The validation of the expression of a few Strategy II uptake genes was also done by qRT-PCR. Our results for qRT-PCR validate our inference from the RNA-seq analysis. Almost all of the Strategy II genes tested for their expression showed a very high fold expression in Fe-starved roots as compared with the control (Supplementary Fig. S1). The highest expression was obtained for ZIFL4, DMAS1, NAAT1, and NAS1 which are the prime components for Strategy II-mediated uptake. Other transcripts coding for thaumatin-like proteins, etc. were also induced under Fe-limiting conditions.

Fig. 4.

Top up-regulated and down-regulated genes in –Fe conditions, annotated via KOBAS using rice RAP-DB/RefSeq annotations as reference. The heat map shows the top 50 genes that are highly up-regulated (red left panel) and down-regulated (green right panel) identified in wheat roots under –Fe conditions with respect to control. For expression analysis, FPKM values were obtained using Cufflinks, and CuffDiff was used to identify DEGs by calculating significant changes in transcript expression between the stressed and normal samples (FDR ≤0.05).

Transcripts of genes involved in the Strategy I mode of Fe uptake were also analysed (Supplementary Table S6). An important component of the Strategy I pathway, the H+-ATPase (AHA) subfamily genes, was not differentially expressed under Fe-limiting condition (Supplementary Table S6). This was in agreement with the qRT-PCR analysis where either no change or down-regulation of AHA genes was observed (Supplementary Fig. S1) in wheat roots. Metallo-reductases are important components of Strategy I represented by ferric-chelate reductase (FRO). Most of the transcripts coding for wheat FROs do not show significant changes in starved roots as compared with the control except for one transcript. Interestingly, iron-regulated transporters (IRTs; TRIAE_CS42_7DS_TGACv1_622068_AA2032200; TRIAE_CS42_4AL_TGACv1_289466_AA0971640) were significantly expressed in Fe-starved wheat roots, suggesting their involvement in wheat under metal stress. The high expression of TaIRT1 was also confirmed by qRT-PCR (Supplementary Fig. S1). This suggests that IRT function might be conserved among the plant species.

In addition to this, 45 genes showed exclusive transcript abundance under Fe starvation (Supplementary Table S3). No transcript for any of these genes was detected in the control root samples, suggesting their high specificity for Fe response. Some of the transcripts responding to Fe starvation code for metallothionein, metal transporters, vacuolar iron transporters (VIT1), and NAS2. This suggests that in wheat few components involved in Fe uptake could actually respond exclusively to Fe starvation. The fold expression levels of genes which were highly down-regulated ranged from –7.1 to –2.85-fold at the level of the log2 scale (Supplementary Table S3). Interestingly, multiple genes encoding nitrate transporters were highly down-regulated under Fe-limiting conditions (Fig. 4, right panel). Two genes coding for cytochrome P450 also showed down-regulation. Among the others, genes coding for histone deacetylase, peptidases A1 domain-containing protein, cinnamyl alcohol dehydrogenase, and dirigent protein show significant down-regulation.

Functional enrichment network of Fe starvation-related genes

GO annotations and classification of DEGs were performed to obtain an overview of processes that are representative of cellular, molecular, and biological functions. The analysis was further extended to cluster analysis using Cytoscape plugins, BINGO, and Enrichment Map (Maere ; Merico ). A total of 5854 DEGs were found to be differentially expressed in response to Fe starvation; significant GO categories were assigned to all DEGs. The DEGs annotated for GO terms were visualized using the WEGO tool (Fig. 5A; Supplementary Table S7). The most enriched GO terms in the ‘cellular component’ category were membrane and intracellular organelle. Other over-represented terms included catalytic activity and ion as well as organic compound binding for the ‘molecular function’ category, while metabolic processes related to nitrogen as well as other cellular processes were observed for the ‘biological process’ category (Fig. 5A). Overall, catalytic activity and binding activities were the most significantly enriched GO terms in the –Fe condition (Fig. 5B). Further mapping of up-regulated genes to databases such as the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (Xie ) and MapMan (Thimm ) revealed enrichment related to phenyl-propanoid biosynthesis, amino acid biosynthesis, carbon metabolism, and glutathione metabolism (Fig. 5C; Supplementary Fig. S2). The role of glutathione in response to Fe starvation is intriguing and deserves further investigations.

Fig. 5.

Gene Ontology (GO) categorization of the differentially expressed genes and its analysis. (A) WEGO plot describing GO annotation and classification of DEGs, with the left y-axis showing the percentage of genes belonging to the respective GO terms (red bars for down-regulated genes and grey bars for up-regulated genes) and the right y-axis depicting the number of both up- and down-regulated genes. Percentage and number of genes were calculated for the three broad main categories listed on the x-axis. (B) Enriched GO terms in DEGs under the –Fe conditions; the y-axis depicts the significance of GO term enrichment. (C) The top 20 enriched KEGG pathways represented by the up-regulated genes under iron starvation condition. The x-axis shows names of the pathways; the y-axis represents the number of genes enriched in respective pathways.

The GO category ‘integral to membrane’ was found to be associated with pathways clustered into ‘metal ion and trans-membrane transport’- and ‘photosynthesis’- related GO categories (Fig. 6A). Clustered pathways involved in ‘response to nutrient levels’ and ‘sulfur amino acid metabolism’ were also enriched (Fig. 6A). A few other up-regulated clustered pathway genes involved fall into the category ‘response to phosphate starvation’. In the enriched network of down-regulated DEGs, genes involved in Fe and other metal homeostasis were associated with two different clusters, mainly including genes related to lipid, ketone, and carboxylic metabolism, in addition to those involved in nitric acid and the salicylic acid response (Fig. 6B). A second cluster contains transcripts related to purine and adenine nucleotide binding and pyrophosphatase activity. Thus, the whole set of significant up- and down-regulated DEGs were clustered in distinct cellular processes, suggesting a differential transcriptional response under Fe starvation.

Fig. 6.

Co-expression/hub genes and function enrichment network for identified DEGs. Function enrichment network for DEGs associated under iron starvation in wheat roots with high significance (FDR ≤0.05) for (A) up-regulated and (B) down-regulated genes. Enriched GO functional categories are clustered with correlated DEGs and represented by node circles.

Identification of transcriptional regulatory genes during Fe starvation

To address the important question of how the –Fe signal is connected to the transcriptional machinery in wheat, genes encoding TFs were identified. Previously, the TFs involved in the response to Fe starvation were identified in the model plants Arabidopsis and rice, including POPEYE (PYE), basic helix–loop–helix (bHLH), ethylene insensitive-3 (EIN3), and EIN3-like1 (EIL1) (Long ; Bauer ; Ogo ; Ivanov ; Bashir , Mai ). Consistently, our transcriptome analysis showed that all homoeologues for PYE encoding a bHLH transcription factor and BRUTUS-like (BTS) known to encode a putative hemerythrin E3 ligase protein in Arabidopsis and rice were highly expressed in wheat roots under Fe starvation (Supplementary Table S8). Analysis of these categories of genes led to identification of 41 significantly up-regulated TF family members (log2FC >1) (Supplementary Table S8). The TF family members belong to categories such as APETALA2 ethylene-responsive element-binding proteins (AP2/EREBP), WRKY, C2H2, zinc finger proteins (including C3HC4), NAM, bHLH, hemerythrin, and U-box (Fig. 7A). Out of these, transcripts coding for AP2/EREBP followed by WRKY and bHLH were the most abundant in wheat. Genes encoding AP2-EREBP, WRKY, and the C2H2 type of TFs also showed high predominance (Supplementary Table S8; Supplementary Fig. S3). However, genes coding for TFs such as MADS, PHD, and homeobox (HB) domain-containing proteins were highly represented in the list of down-regulated transcripts. Finally, to ascertain regulatory functions under Fe starvation, enrichment mapping was performed for the significantly expressed DEGs which were classified into TFs. Interestingly, genes involved in regulation of gene expression were clustered with primary and nucleotide metabolism-related genes in both positive and negative correlations. This functional cluster was positively correlated with TFs involved in GO categories such as ‘vacuolar transport’ and ‘DNA topoisomerase III activity’. However, two separate clusters ‘cellular response to auxin stimulus and signalling’ and ‘cell-fate specification’ were negatively correlated with ‘regulation of gene expression’ and ‘primary and nucleotide metabolism’ (Fig. 7B). This suggested an important regulatory role for TFs involved in auxin signalling and cell-fate specification in Fe starvation response control, thereby modulating the network of Fe homeostasis.

Fig. 7.

Transcriptional factors (TFs) significantly associated with Fe starvation (FDR ≤0.05) in wheat roots. (A) List of TFs differentially expressed in response to Fe starvation stress. Blue bars represent up-regulated and orange bars represent down-regulated TFs. (B) Co-expression/hub genes and network analysis using Fe-responsive TFs (FDR ≤0.05) in wheat roots. Circles indicate the processes associated with it (green for down-regulated genes; blue for up-regulated genes). Enriched GO functional categories are clustered with correlated TFs and represented by node circles.

GSTs are involved in the response to Fe starvation in wheat

Glutathione-mediated conjugation of multiple metabolites plays an important role during metal stress (Zhang ). Our transcriptome data revealed the enhanced expression of multiple GSTs in Fe-starved wheat root as compared with control (Supplementary Table S9). To correlate the expression response with its enzymatic activity, the temporal response of GST was measured in wheat roots under –Fe conditions. Using the GST functional assay, the activity was determined in wheat roots of plants grown for 10, 15, and 20 d under Fe starvation. Our activity assays showed a significant increase in the GST activity under –Fe conditions, which peaked at 10 d and 20 d after the beginning of Fe starvation compared with control plants (Fig. 8A). Therefore, our experiments validate the increased GST transcript abundance with enhanced glutathione activity. These results indicate an important role for glutathione in response to Fe starvation in wheat roots.

Fig. 8.

Measurement of GST activity and metabolite profiling of wheat roots subjected to Fe starvation. (A) Glutathione S-transferase activity of wheat roots under Fe starvation (–Fe) and control (+Fe) conditions. (B) Metabolite profiling of amino acids, sugars, polyols, organic acids, and related compounds. Change in abundance of significant (P<0.05) metabolites identified by GC-MS in Fe-starved roots. Variation in abundance of each metabolite is represented by log2 fold values of the response ratio (–Fe/+Fe) of metabolite concentrations. Values are means of three biological replicates, with the bar representing the log ratio of the SE; * indicates a significant difference at P<0.01; #indicates a significant difference at P<0.05.

Fe starvation causes an accumulation of organic acids and polyhydroxy acids in wheat

To obtain a comparative insight into the metabolite profile of wheat roots of plants grown in the absence or presence or Fe, GC-MS profiling analysis was performed. Metabolites were extracted from the roots of three replicate pools of plants each containing eight seedlings. Qualitative processing of each chromatogram for peak area and identification was performed with MassHunter version B.05.00 software coupled with the NIST11 compound library. The compound annotation was determined by comparing individual resolved peaks with library searches based on mass spectra and compound chromatographic retention indices. Interestingly, analysis resulted in the identification of 54 metabolites, and a further 39 annotated metabolites were analysed for their response ratio (Supplementary Table S10). To compare the change under Fe starvation, metabolite abundances (–Fe roots/control roots) were calculated and expressed as log2FC values. Fe starvation significantly affected accumulation of 22 metabolites that include organic acids, polyhydroxy acids, amino acids, and some of the sugars, fatty acids, and phosphates (Fig. 8B). Amongst treatment-specific changes, a few organic acids such as fumaric acid, acetic acid, and malonic acid showed a higher level of accumulation in Fe-starved roots as compared with the control. In contrast, citric acid, malic acid, valeric acid, and aconitic acid were significantly lower in Fe-starved roots when compared with control samples. In Fe-starved roots, the accumulation of amino acids, mainly l-valine, showed a significant increase, while hydroxylbutyric acid and pyroglutamic acid were decreased in comparison with control roots. Polyhydroxy acids such as gluconic acid and glyceric acid were also significantly high in Fe-starved roots, whereas the level of hexonic acid and arabinoic acid was found to be low. Taken together, our results showed that during Fe starvation wheat roots undergo reprogramming for metabolic changes to maintain the Fe flux.

Discussion

How plants maintain nutrient homeostasis is a fascinating question in plant biology. On this theme, Arabidopsis thaliana, with its fully sequenced small genome, has provided some basic information. Developing our knowledge on nutrient homeostasis in crops, especially those having complex genomes, including hexaploid wheat, is more challenging. Recently, the release of the genome sequence, generation of transcriptome data, and the advancement of different functional tools have provided much needed impetus in this direction (Borrill ). The current study was undertaken to gain insight into the response of hexaploid wheat exposed to Fe starvation which is known to severely affect crop production. In this study, the transcriptome of wheat (cultivar C306) in response to Fe starvation was generated. Our analysis revealed that wheat: (i) utilizes primarily the Strategy II mode of Fe uptake; (ii) accumulates transcripts coding for the methionine salvage pathway coupled with enhanced GST activity; and (iii) accumulates specific metabolites including fumarate, acetate, malonate, and xylofuranose, to efficiently mobilize soil Fe. Interestingly, our systematic analysis of expression data revealed that transcripts show an induction bias for A and B genomes of wheat in response to the –Fe condition. Overall, this work provides the first comprehensive insight at the molecular level for wheat roots during Fe deficiency.

Wheat plants subjected to Fe starvation show physiological defects such as a decrease in root growth. This phenotype was consistent with a previous report showing a negative impact on the root growth of wheat seedlings under Fe stress (Garnica ). Our analysis also suggested down-regulation of a few nitrate transporters for the time point studied. Homology searches indicated the identified nitrate transporters in wheat as TaNRT2.2 and TaNRT2.3 (Buchner and Hawkesford, 2014). These observations suggest that Fe starvation represses nitrate transporters, which in turn leads to lower accumulation of nitrate levels in the roots (Supplementary Fig. S4). Although the impact of Fe starvation on nitrate metabolism has not been studied in detail yet, a similar decrease in nitrate levels was also observed in cucumber shoots subjected to –Fe conditions (Borlotti ). This observation further reinforces the existence of an interaction between macro- and micronutrients in plants, which has recently gained attention (Rouached and Rhee, 2017; Bouain ).

Due to the low availability of micronutrients in soil, plants are equipped to recruit components that could participate in the two major strategies for Fe uptake (Kobayashi and Nishizawa, 2012). Cereals such as maize and rice predominantly utilize the Strategy II mode of uptake, unlike Arabidopsis, that employs the Strategy I mode of Fe uptake (Hell and Stephan, 2003; Li ; Zanin ). Nevertheless, a comprehensive study to identify the molecular players involved during Fe-limiting conditions in wheat is lacking. Our RNA-seq-based analysis strongly suggested an increase in the transcript abundance of genes for the Strategy II mode of Fe uptake. With the exception of conserved genes such as IRT1 (iron regulated transporter), prime genes for Strategy I uptake mechanism such as FRO (ferric chelate reductase) and proton H-ATPase (AHA-like) were also present during starved condition but at very low abundance, and no FRO activity has been reported for cereals. During this study, wheat roots subjected to Fe starvation also do not show any FRO activity (data not shown). Therefore, our results reveal that wheat utilizes Strategy II-mediated uptake of Fe and induces IRT genes that might be conserved for their function across the plant kingdom.

The series of events during the Strategy II uptake mechanism used by Gramineae involves secretion of PSs that facilitate conjugation of Fe to form Fe complexes. These Fe–PS complexes are imported in roots for remobilization to the different developing organs. As represented in Fig. 9, the transport of these Fe complexes occurs via membrane transporters via ZIFL genes coding for TOM proteins (Nozoye , 2015). Multiple wheat ZIFL genes show high transcript accumulation and are the closest homologues of OsTOM1 from rice, thereby hinting at its role in Fe acquisition. Utilizing gain- and loss-of-function approaches, OsTOM1 has been demonstrated to be a DMA effluxer for its role in enhancing mobilization and thereby improving Fe uptake (Nozoye ). Altogether, based upon the high similarity and response under Fe starvation, it is tempting to propose wheat ZIFL (ZIFL4) as a functional transporter of DMA. However, other ZIFL genes were also significantly expressed, but a functional wheat transporter for DMA needs to be deciphered. PS secretion was also observed during our analysis of the roots of C306 exposed to Fe starvation. In this study, a high PS release of ~40–45 nmol g−1 root biomass (3 h)−1 was observed, reinforcing our RNA-seq data showing that DMA biosynthesis genes and its transporters were highly up-regulated. Control plants show accumulation of ~1–2 nmol g−1, suggesting that wheat releases a basal level of PSs to maintain the constant flux of Fe in roots. Upon careful analysis of the RNA-seq data, it was observed that multiple such efflux transporters are represented in the list that are highly abundant under Fe starvation. During our study, we identified TaDMAS1 coding for DMA synthase that was highly responsive to Fe starvation (Supplementary Table S5). Previously, TaDMAS1 was reported as a gene that was broadly expressed across the tissue and was regulated during Fe-starved conditions. These observations support the notion that TaDMAS1 has potential to enhance the seed Fe storage capacity (Beasley ). During our study, GO term enrichment was also observed for the DMA biosynthesis pathway, transmembrane transporters, and cellular response that reinforce the importance of these functional categories. Surprisingly, genes pertaining to photosynthesis were also observed. This anomaly could be explained by the presence of an active photosystem in the basal roots of our collected wheat tissue, similar to that observed in Arabidopsis due to the changes in the auxin/cytokinin ratio (Kobayashi ). Secondly, during this study, the wheat roots under Fe starvation show significant accumulation of two Golden-2 like (GLK) TFs (~logFC 1.54 and 1.55; Supplementary Table S3). GLKs are known to improve the phototrophic performance of roots and thereby enhance root photosynthesis (Kobayashi ). These reasons may account for the enhanced accumulation of transcripts related to photosynthesis along with other metabolic-related genes.

Fig. 9.

Schematic representation describing the core components involved in Fe starvation. The red font indicates the genes/metabolites which were highly up-regulated/showed high accumulation during our study, whereas the green font indicates down-regulated/low accumulation. Methionine salvage pathway bins indicate the gene expression levels for the transcript indicated next to it. Red and blue bins near pathways and steps represent the up- and down-regulation of related genes, respectively. FBP, fructose-1,6-bisphosphatase I; PFP, diphosphate-dependent phosphofructokinase; ALDO, fructose-bisphosphate aldolase; GAPDH, glyceraldehyde 3-phosphate dehydrogenase; PGAM, 2,3-bisphosphoglycerate-dependent phosphoglycerate mutase; PDC, pyruvate decarboxylase; GLYK, d-glycerate 3-kinase; GlyA, glycine hydroxymethyltransferase; ltaE, threonine aldolase; BMT2, homocysteine S-methyltransferase; metE, homocysteine methyltransferase; metK, S-adenosylmethionine synthetase; SamDC, S-adenosylmethionine decarboxylase; SRM, spermidine synthase; MTN, 5'-methylthioadenosine nucleosidase; mtnK, 5-methylthioribose kinase; mtnA, methylthioribose-1-phosphate isomerase; DEP1, enolase-phosphatase E1; mtnC, enolase-phosphatase E1; mtnD, 1,2-dihydroxy-3-keto-5-methylthiopentene dioxygenase; TyrAT, tyrosine aminotransferase; NAS, nicotianamine synthase; NAAT, nicotianamine aminotransferase; DMAS, 3''-deamino-3''-oxonicotianamine reductase; YSL, Yellow Stripe Like; maiA, maleylacetoacetate isomerase; AAT, aspartate aminotransferase; TAT, tyrosine aminotransferase; HPD, 4-hydroxyphenylpyruvate dioxygenase; GGT, gamma-glutamyltranspeptidase; AOS3, hydroperoxide dehydratase; OPR, 12-oxophytodienoic acid reductase; ACX, acyl-CoA oxidase; MFP2, enoyl-CoA hydratase; ACAT, acetyl-CoA acyltransferase 1.

Wheat is a hexaploid crop with three genomes and therefore studying the homoeologue induction and expression bias could provide an insight into the regulation of the transcript under a given biotic and abiotic stress. Therefore, we subjected our RNA-seq data to these analyses. In-depth analysis confirms the minimal suppression of the expression bias for the D-genome-derived transcripts as compared with A or B. These observations support the previous analysis wherein suppression of the D genome was significantly less frequent in multiple tissue (Ramirez-Gonzalez ). Previously, homoeologue expression was also studied for wheat under infection with Fusarium pseudograminearum, wherein the expression bias was observed for the B and D subgenomes (Powell ). Utilizing RNA-seq data, changes in the pattern of homoeologous gene expression were also reported in bread wheat (Leach ). Therefore, such studies are important to pin-point actively expressed/induced homoeologues that could be targets for genetic improvement. However, the genome-dependent PS release and its influence by the genotypic variability in the complex wheat genome need to be addressed. Previous studies have shown that the PS release mechanism in wheat is effective when the three genomes A, B, and D contribute synergistically to induce its biosynthesis and release. In general, genome induction bias for the A and B genome was observed for the hexaploid wheat under Fe stress (Fig. 3B). Interestingly, these genomes were also suggested to be important during starvation of other micronutrients such as Zn (Tolay ). These observations provide a clue to the importance of the A and B genomes during micronutrient starvation in hexaploid wheat when compared with their progenitor. The PS release intensity is also dependent on the complementary action of the genotypic ploidy with the following pattern: T. aestivum (BBAADD)>T. dicoccum (BBAA)>T. monococcum (AA)>Ae. tauschii (DD) (Ma ; Tolay ).

Multiple genes for the salvage pathway were highly up-regulated, including NAS1, NAS2, and NAAT1 (Supplementary Tables S3, S5). The contribution of these transcripts arises from different genomes, suggesting that all the genomes of wheat are capable of responding to the starvation condition. Based on our RNA-seq analysis, it is likely that the gene coding for TaYS1A and TaYSL9 could be the putative transporter for the Fe–siderophore complex. Previous studies have indicated the high response of TaYS1A under Fe starvation in roots and shoots (Kumar ). Of the closest orthologues of wheat, YS1A, HvYS1A, and ZmYS1A, HvYS1A has been shown to be involved in specific transport of Fe, whereas ZmYS1 was reported to have broad substrate specificity (Schaaf ). Our RNA-seq data also revealed a conserved function of genes from monocots and dicots during Fe regulation, and known important Fe regulators such as PYE and BRUTUS also showed differential expression in wheat roots (Ivanov ). TFs spanning multiple gene families such as MYB (Rubio ), bHLH (Colangelo ; Jakoby ; Ogo ), C2H2, NAC, AP2-ERB (Kim ), and WRKY (Devaiah ) were characterized in Arabidopsis and rice for their involvement in nutrient uptake and were found to be highly represented in root of Fe-starved wheat. Similarly, hemerythrin motif-containing proteins were also overexpressed in our study. As in the case of other graminaceous plants such as rice, hemerythrin domain-containing proteins including the OsHRZ homologue for Arabidopsis BRUTUS was also differentially perturbed in –Fe (Kobayashi ). Overall, based on the information obtained from our studies, future research focus should be able to gain better insights into the molecular functions of these genes using a TILLING population in wheat and/or a heterologous system.

Changes in metabolite content in response to Fe deficiency can lead the plant to undergo adaptive processes to maintain Fe homeostasis (Schmidt ). A schematic representation of data provides a comprehensive summary of the metabolic and transcriptional changes that occur during Fe starvation in wheat roots (Fig. 9). The mobilization and transport of Fe inside the plant tissue depend on its speciation and complexation form. Depending on the tissue, Fe could form a complex with multiple metabolites such as nicotianamine, malate, or citrate that facilitate the transport through the phloem and xylem (Palmer ; Grillet ). Our analysis of the metabolome of roots suggested the predominance of additional organic compounds such as fumarate, acetate, and glycerate, thereby suggesting that they could also be involved in Fe complexation to assist in long-distance transport and remobilization to different tissues (Fig. 9). These metabolites are reported for their ability to complex with the Fe(III) state during transport in xylem vessels as detected during X-ray absorption near-edge structure spectroscopy-based studies (Terzano ). Similarly, the accumulation of glycerate could be responsible for added thrust towards the glycerate–serine pathway during Fe starvation (Fig. 9). The glycerate–serine pathway is important to generate methionine and the precursor (S-adenosylmethionine) of DMA biosynthesis. Consistent with the metabolome data, the transcripts encoding these proteins were highly up-regulated in our study as presented by the bins (Fig. 9). Members of ABC-B transporters were also reported to participate in oxidative stress that was linked with metal stress including Fe starvation (Kispal ; Schaedler ). During our analysis, we have identified a TAP-type, ABCB subfamily transporter gene referred to as TaABCB25 that was highly up-regulated during Fe starvation (Fig. 4). The closest orthologues for wheat ABCB25 transporters were identified in rice (OsABCB25) and Arabidopsis (AtABCB27) (Supplementary Fig. S5A). The temporal induction of TaABCC25 was also re-confirmed in roots subjected to Fe starvation (Supplementary Fig. S5B). Arabidopsis thaliana ABCB27 was shown to have the ability to transport glutathione conjugate metabolite precursors for Fe–S cluster assembly (Schaedler ). Based on the presence of domain structure, these ABC-B transporters from rice, Arabidopsis, and wheat qualify as half-size transporters (Supplementary Fig. S5C). In graminaceous species, glutathione-related activity was observed to be highly up-regulated under Fe-deficient conditions, supporting our observation in the current study (Bashir ; Nishizawa, 2007). Enhanced GST activity in wheat under Fe starvation and the abundance of TaABCB25 in roots warrants further study of this gene for its role in mobilizing micronutrient uptake. Glutathione and GST play an essential role during Fe starvation-related responses in Arabidopsis and also quench reactive molecules to protect the cell from oxidative damage and prevent chlorophyll loss (Ramirez ; Shanmugam ; Kumar and Trivedi, 2018). Our analysis also suggests that high transcript and enzyme activity of GST could be linked to primary metabolism since it accumulates fumarate (Fig. 9). In Arabidopsis, it has been shown that GST catalyses glutathione-dependent isomerization of maleyl-acetoacetate into fumaryl-acetoacetate, eventually leading to accumulation of fumarate and acetoacetate (Dixon ). Additionally, GST is also linked to the conjugation of the reduced glutathione (GSH) forms to various substrates to make them water soluble. The sulfur-derived metabolites and glutathione are well known to be involved in biosynthesis of methionine, a well-known precursor for DMA biosynthesis (Forieri ). Altogether, this supports the notion that in wheat, the role of glutathione is important during Fe starvation, and ABC transporters might play a significant role in mobilization of metabolites/organic acids or specific siderophores.

In conclusion, in this study, the core components for the Fe starvation response in hexaploid wheat have been identified so as to provide a better understanding of the molecular events that participate during Fe starvation in wheat roots. In particular, we provide lines of evidence for the role of GST in the response to –Fe in roots. Complementary approaches, analytical and transcriptome analysis, reinforce the importance of primary metabolism for reprogramming organic acids and amino acids, thereby leading to Fe homeostasis in wheat. The information here not only helps to design strategies to improve plant response to Fe starvation in wheat, but will also foster Fe uptake and accumulation studies, which are required to boost productivity and grain nutritional quality through Fe biofortification programmes.

Supplementary data

Supplementary data are available at JXB online.

Fig. S1. qRT-PCR validation of selected genes from the DEGs in Fe-deficient roots after 20 d of starvation.

Fig. S2. Overview of genes modulated by Fe starvation in wheat roots.

Fig. S3. MapMan visualization depicting the differentially expressed transcription factor families.

Fig. S4. Estimation of nitrate levels under Fe starvation using the salicylic acid method.

Fig. S5. Characterization of the wheat ABCB25 transporter.

Table S1. List of primers used in the current study.

Table S2. Summary of filtered and mapped reads for each sample.

Table S3. DEGs in response to Fe starvation in wheat roots.

Table S4. List of 8473 homoeologue triads that were used for homoeologue induction and expression bias analysis.

Table S5. Expression profiles of genes/gene families involved in the Strategy II mode of Fe uptake.

Table S6. Expression profiles of genes/gene families involved in the Strategy I mode of Fe uptake.

Table S7. Gene Ontology analysis of up- and down-regulated genes in response to Fe starvation.

Table S8. Expression profiles of different transcription factor genes that are differentially up- and down-regulated.

Table S9. Expression profiles of genes involved during the glutathione-mediated detoxification process.

Table S10. GC-MS analysis of wheat roots subjected to 20 d of Fe starvation.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.