Iron deficiency triggers rapid changes in transcription and splicing patterns
To gain insights into transcriptional and post-transcriptional alterations in response to Fe deficiency, Arabidopsis plants were subjected to 0.5, 6, and 12 hours of growth in Fe-free media and subsequent transcriptomic profiling against plants that were transferred to fresh Fe-replete nutrient solution using the RNA-seq methodology. The Illumina HiSeq 4000 sequencing system was adopted to acquire approximately 70-80 million reads for each library with a read length of 100 base pairs (Supplementary Table S1). On average, a total of circa 27,000 genes was expressed in both roots and shoots, of which subsets of 1,161 (roots) and 1,027 (shoots) were defined as differentially expressed between Fe-deficient and Fe-sufficient control plants at one or more of the three sampling time points with relevant expression levels (RPKM > the square root of the mean expression value of the whole dataset), P < 0.05, and a fold-change > 2 after normalization using the TMM (Trimmed Mean of M-values) method (Fig. 1a). In both roots and shoots, differentially expressed genes (DEGs) showed a relatively small overlap among the experimental time points, indicative of a highly dynamic and temporally distinctive response to Fe deficiency (Fig. 1c).
Changes in the pattern of alternative splicing (AS) upon Fe deficiency were analysed with the aid of the Read Analysis and Comparison Kit in Java (RACK J) software toolbox39. Applying the same thresholds used for the classification of DE genes and considering only genes in which the AS feature was observed in all three replicates, a subset of 7,517 genes was classified as differentially alternatively spliced (DAS), exceeding the number of differentially expressed genes (DEGs) by more than fourfold (Fig. 1f). The largest fraction of DAS genes (82%) produced transcripts that harboured differential intron retention (DIR) features, a smaller subset (15%) carried differential donor/acceptor (DDA) sites, and 2.8% of the DAS transcripts exhibited differential exon skipping (DES) after exposure to Fe deficiency (Fig. 1e). A considerable subset of DAS genes carried two or more different AS features. Genes with DAS features showed distinct overlaps among the different experimental timepoints which depended on the type of AS with DIR being more conserved than DEG and the overlap of DDA genes comparable to that observed for DEGs (Fig. 1c, d). With the exception of DES, in which enhanced and repressed features varied over time, DAS showed a neutral outcome of splicing efficiency upon exposure to Fe deficiency, with about similar proportions of genes carrying enhanced or reduced splicing features (Supplementary Fig. S1).
Fe deficiency induces a series of concatenated responses in roots and shoots
Transfer of plants to Fe-deplete media triggered the consecutive induction of a suite of distinctly timed responses (Fig. 2). Induction of these processes was almost simultaneously observed in roots and shoots, suggesting that minor changes in Fe supply suffice to sense Fe deficiency in all plant parts. After six hours of exposure to Fe-deficient conditions, induction of the basic Fe uptake machinery comprising the Fe chelate reductase FRO2 and the Fe2+ transporter IRT1, vacuolar sequestration of excessive cytosolic levels of secondary-substrate metal cations such as Mn2+ and Zn2+, and downregulation of Fe import into the vacuole were the most prominent transcriptionally regulated responses in roots (Fig. 2a). Induction of these processes was accompanied by increased expression of a suite of genes encoding regulators such as the IRONMAN peptides IMA1 and IMA2 and the transcription factor POPEYE (PYE) (Fig. 2a). A module that was induced slightly later comprised genes mediating the mobilization of Fe in the rhizosphere via secretion of protons and Fe-mobilizing coumarins (Fig. 2b). Induction of Fe mobilisation was accompanied by increased expression of the transcription factors MYB72 and MYB10, which were shown to be crucial for plant survival in soils with severely restricted Fe availability (Fig. 2 b)40. Notably, sampling at very early stages of Fe deficiency (0.5 h) revealed expression changes in a direction which was antagonistic to that observed at later stages, suggesting that rearrangements of the transcriptional machinery cause transient perturbations in gene expression that later result in robust transcriptional regulation of these genes (Supplementary Tables S2-5).
In shoots, exposure to Fe-deplete media for 6 hours induced Fe transport across the plasma membrane via OPT3 and IRT3 and expression of a variety of regulators such as PYE, the E3 ligase BRUTUS (BTS), subgroup Ib bHLH proteins, and IMA1-IMA3 (Fig. 2c). Similar to root cells, sequestration of transition metals into the vacuole was induced at this time point. In addition, altered transcription of genes involved in ROS homeostasis (CGLD27, ENH1) was observed in shoots (Supplementary Table S2). After 12 hours, increased trans-plasma membrane transport of Fe-nicotianamine (NA) by YSL1 and Cu2+ via COPT2 can be inferred from the induction of the cognate genes at this time-point. Also, chlorophyll metabolism (HEMA1, NYC1) and ROS homeostasis (NEET) were affected by growth on Fe-deplete media (Fig. 2d; Supplementary Table S2).
In addition to genes directly involved in Fe homeostasis, in roots, and to a lesser extent in shoots, short-term exposure to Fe deficiency caused a dramatic increase in the transcription of genes associated with jasmonic acid (JA) biosynthesis and signalling (Fig. 3; Supplementary Table S2 and S3), a response which was strictly restricted to the 6-hour time point. In particular, the first steps of JA biosynthesis and the expression of a suite of transcriptional regulators, referred to as jasmonate ZIM-domain (JAZ) proteins, were strongly induced in response to Fe deficiency (Fig. 3).
DAS orchestrates rerouting of the central carbon metabolism
Glycolysis is a key process in energy production and provides metabolic intermediates for biosynthetic processes, storage, and anaplerotic processes. A subset of 75 genes associated with glycolysis was regulated by the Fe regime, almost exclusively by DAS (Supplementary Table S4). Notable exceptions from this pattern were the first and the last step in glycolysis, catalysed by phosphofructokinase (PFK) and pyruvate kinase (PK), which are robustly upregulated at later stages of iron deficiency in transcriptional surveys9,11. Analysis of PFK1 expression by RT-qPCR showed increasing induction over the first three days of Fe deficiency in roots and, to a much lesser extent, in shoots (Fig. 4a, c). Expression of the Fe status marker bHLH38 was increased over the first three days of Fe deficiency with a higher transcript level in shoots, indicating that the lower induction of PFK1 in shoots was not associated with a healthier Fe status of leaf cells (Fig. 4c). Since PFK is catalysing the rate-limiting step in glycolysis, it can be assumed that glycolytic flux is particularly increased in root cells. The final step in glycolysis, the PK-mediated conversion of phosphoenolpyruvate (PEP) to pyruvate, was repressed in roots via enhanced IR of the cytosolic PK isoform At5g08570. In shoots, transcripts of this isoform carried reduced IR upon exposure to Fe deficiency. RT-qPCR analysis revealed decreased transcript abundance of this isoform in both roots and shoots one day after exposure to Fe-deplete media (Fig. 4d). Downregulation of another PK-encoding gene in roots (At3g49160) was reported in several transcriptomic studies7,11, further supporting the supposition that PK activity is decreased upon Fe deficiency.
For several steps of root glycolysis, DAS appear to precede changes in enzyme activity or abundance. Fructose-bisphosphate aldolase (FBA) catalyses the reversible aldol cleavage of fructose-1,6-bisphosphate, yielding dihydroxyacetone phosphate. In roots, the cytosolic isoform FBA8 exhibited both repressed IR and DA features (Fig. 4a). Also, the cytosolic phosphoglycerate mutase (PGM) isoforms iPGAM1 and iPGAM2 showed reduced IR upon Fe deficiency (Fig. 4a). In a previous study, we observed increased phosphorylation of FBA8 and iPGAM1 protein and accumulation of the (cytosolic) FBA4 protein in roots upon prolonged Fe deficiency23,41, indicating Fe deficiency-induced increase in enzyme activity at later stages of Fe deficiency. Moreover, FBA and PGM accumulated in Fe-deficient roots of Beta vulgaris and Cucumis sativa19,42, suggesting that increased abundance of these enzymes in response to Fe deficiency is conserved across strategy-I species.
The anti-directional regulation of PFK and PK in roots of Fe-deficient plants implies imbalances of glycolysis intermediates under conditions of Fe starvation. To further investigate this matter, we determined the concentrations of glucose-6-phosphate (G6P), phosphoenolpyruvate (PEP), and pyruvate over the first three days of Fe deficiency via UHPLC-MS analysis. In roots, but not in shoots, the level of both PEP and pyruvate was increased three days after the onset of Fe deficiency, indicating severe perturbances of pyruvate metabolism at later stages of Fe deficiency (Fig. 4b). Other glycolysis intermediates such as G6P or fructose-6-phosphate (F6P) did neither show pronounced changes in roots nor in shoots upon short-term or extended exposure to Fe deficiency (Fig. 4b), suggesting that the concentrations of theses intermediates are in steady-state during Fe deficiency.
The TCA cycle malfunctions in roots of Fe-deficient plants
A subset of 34 of TCA-related genes was responsive to the Fe regime (Supplementary Table S5). Similar to what has been observed for glycolysis, genes encoding enzymes of the TCA cycle showed limited transcriptional control and were predominantly regulated by DAS (Supplementary Table S5). In roots, all TCA cycle metabolites under study accumulated three days after the onset of Fe deficiency, with citrate, malate, and succinate being most abundant (Supplementary Fig. S3). In particular, citrate levels were strongly increased; no such accumulation was observed in shoots. Citrate is synthesized in the first committed and pace-making step of the cycle by condensation of oxaloacetate and acetyl CoA, catalysed by citrate synthase (CSY). However, the mitochondrial isoforms CSY4 and CSY5 were neither regulated by DAS nor by DE, suggesting other causes for the increased citrate levels. Aconitase (ACO), catalysing the subsequent conversion of citrate to isocitrate, harbours an active [Fe4S4]2+ cluster, which may compromise ACO activity in Fe-deficient plants and contribute to the accumulation of citrate. In line with this assumption, ACO2 and ACO3 were transcriptionally downregulated in roots three days after transfer to Fe-deficient media7,9,43, indicating restricted conversion of citrate into isocitrate. Also other enzymes of the TCA cycle such as citrate synthase and isocitric dehydrogenase were shown to affect by the iron regime13, further supporting the notion that the TCA cycle is compromised or truncated under conditions of Fe deficiency. This scenario is supported by a lack of consistent and mostly repressive DAS-regulation of all steps of the TCA cycle (Supplementary Fig. S3).
The routes of pyruvate metabolization differ between roots and shoots
Pyruvate derived from glycolysis can be either converted to acetyl-CoA by pyruvate dehydrogenase (PDH) in the mitochondrial matrix, or decarboxylated to acetaldehyde through pyruvate decarboxylase (PDC) in the cytosol. PDH is a multienzyme complex composed of three enzymes, pyruvate dehydrogenase (E1), dihydrolipoyl transacetylase (E2), and dihydrolipoyl dehydrogenase (E3). In roots, isoforms of all three enzymes showed generally reduced IR upon Fe deficiency (Fig. 5a). In shoots, a more complex, mostly repressive regulation was observed (Fig. 5b). Alcohol dehydrogenase (ADH) reduces pyruvate-derived acetaldehyde to ethanol and NAD+. In roots, short-term exposure to Fe deficiency resulted in repressed IR of both PDC1 and ADH1 transcripts, indicative of increased ethanolic fermentation (Fig. 5a). Also, the class III type alcohol dehydrogenase ADH2 (HOT5) and the putative ADH At4g22110 showed DAS features that were mostly repressed upon Fe starvation. By contrast, in shoots ADH1 carried chiefly enhanced DIR features (Fig. 5b). Determination of ADH1 transcript levels in roots by RT-qPCR revealed a steep increase after two days of Fe deficiency (Supplementary Fig. S4), matching transcriptomic studies carried out at later stages on roots of Fe-deficient plants7,44. These data suggest that the ADH-mediated fermentation route is supported in roots, but not in shoots of Fe-deficient plants.
An alternative fate of acetealdehyde is the conversion to acetic acid, catalysed by a group of aldehydedehydrogenases (ALDHs) in a reaction yielding NADH. In shoots, three NAD-dependent ALDHs (the mitochondrial family 2 proteins ALDH2B7 and ALDH2B4, and the cytoplasmatic family 3 protein ALDH3H1) showed reduced IR upon Fe starvation (Fig. 5b). In roots, ALDH2B7 transcripts carried enhanced IR at two time points after the onset of Fe-deficient conditions (Fig. 5a), indicating repressed acetic acid formation. It thus appears that different fermentation routes are engaged in roots and shoots, possibly driven by differences in redox regulation and energy status of leaf and root cells.
Besides reduced PK activity, the build-up of toxic pyruvate concentrations in roots of Fe-deficient plants is circumvented by rapid metabolization of its precursor, PEP. PEP is a potent inhibitor of PFK activity and, if present at high levels, decreases glycolysis rates45-47. Increased β-carboxylation of PEP via PEP carboxylase (PPC) in roots is a hallmark of Fe-deficient roots of strategy I plants17. In roots, PPC1 transcripts carried reduced IR features (Fig. 5a), transcription of the gene was induced at later stages of Fe deficiency7. The reverse reaction, the conversion of oxaloacetate to PEP by PEP carboxykinase (PCK), showed enhanced IR (Fig. 5a,b). In concert with these results, PCK1 was found to be downregulated in roots (but not in shoots) of plants subjected to Fe deficiency for three days7,48.
DE and DAS concertedly regulate the ferrome
Strikingly, the vast majority of DAS-regulated genes did not change significantly in expression over the experimental period. In total, less than 1% of the DAS genes were also differentially expressed (Fig. S1e). Plotting fold-changes of the different DAS types versus DE revealed a moderate correlation between DIR and DE (Fig. 6a). By contrast, regulation of gene activity by DDA and DE was almost mutually exclusive (Fig. 6b).
The principle of mutual exclusivity does not seem to apply when only genes that were shown to robustly respond to Fe deficiency at the transcriptional level in previous studies are considered, a subset that has been referred to as the ‘ferrome’49. To revise and update this definition, we surveyed recent public transcriptomic datasets of Fe-deficient plants derived from RNA-seq profiling studies. For roots, seven datasets were analysed, and genes that were found to be differentially expressed in at least four of these studies (three out of five in shoots) were included in the Arabidopsis ferrome. This procedure yielded subsets of 108 and 100 genes in roots and shoots, respectively; a suite of 39 genes was robustly differentially expressed in both roots and shoots (Supplementary Table S6). From 114 DE-regulated ferrome genes, a subsection of 42 (37%) was additionally regulated by DAS (Supplementary Table S2), a percentage that is substantially higher than the average of 3-5% observed in roots at the different sampling points. Most of the genes involved in Fe uptake were transcriptionally regulated; however, some regulators (e.g., GRF11, BSTL2), genes encoding enzymes such as the β-glucosidase BGLU42 which is critical for the deglycosylation of coumarins prior to secretion50, and the Zn/Fe transporter IRT3, were exclusively regulated by DAS. All of these genes were shown to be transcriptionally induced at later stages of Fe deficiency (Supplementary Table S2).
Unexpectedly, in both roots and shoots the percentage of DE-regulated ferrome genes increased over time. While after 0.5 hours the majority of Fe-responsive genes was DAS-regulated, the proportions of DE and DAS were about equal at 6 hours and massively shifted towards transcriptional regulation at 12 hours (Fig. 7a, b). A comparison with a previous survey of genes responsive to a 3-day-period of Fe starvation revealed a steep decrease in the percentage of DAS-regulated genes at later stages of Fe deficiency (Fig. 7b). Comparing different subgroups within the ferrome showed that transcripts encoding enzymes tended to have a higher percentage of DAS than transporters and regulators. After 3 days, transcriptional regulation of metabolic genes was still lower than in the other categories, but the participation of DAS in the overall gene regulation was dramatically reduced to about relative to what was observed in the short term (Fig. 7b). Such decreased participation of DAS in the regulation of Fe homeostasis was also observed when other gene ontologies were considered (Fig. 7c). This analysis also revealed that at later stages of Fe deficiency transcriptional regulation of ferrome genes is prioritized over other processes, indicating that a more severely disturbed Fe metabolism prioritises expression of a specific subset of genes to secure efficient Fe acquisition. It may also be assumed that at this stage, early DAS events have been established as stable changes in protein abundance, which is not monitored in the two transcriptomic studies considered here.
Gene architecture predefines the type of gene regulation
From the survey of Fe-responsive genes, it appears that the various Fe-responsive processes are preferentially controlled by either DE (e.g., JA-biosynthesis and signalling), DAS (glycolysis and TCA cycle), or both DE and DAS (Fe uptake), regulatory patterns that are possibly associated with the amplitude of the response (Supplementary Tables S2-5). To investigate whether genomic traits influence or govern the type of gene regulation, we first analysed the number and length of introns of Fe-responsive genes. Intron length was only moderately correlated with the probability of IR. Introns with a length of more than 200 bp had a slightly higher chance of being retained upon Fe deficiency than shorter introns, reaching a peak at about 700 bp (Supplementary Fig. S2a). The number of introns, on the other hand, appeared to be a strong predictor for AS. When normalized to the average intron number across the genome, the number of genes producing transcripts with DIR, DDA, or DES features showed a steep incline which saturated at an intron number between 10 and 20 (Supplementary Fig. S2b, c). Notably, intron-rich genes harbouring >50 introns were rarely found to be subject to DAS (Supplementary Fig. S2b, c).
As further parameters that presumptively affect the type of gene regulation, we investigated the influence of 5’ and 3’ splice site strength, promoter length, and the number of transcription factor binding sites (TFBSs) on DE and DAS. Plotting the minimum average splice site strength [ASS; (min 5’ + min 3’)/2] of Fe-responsive genes against the intron number showed that weak splice sites were generally associated with lower intron numbers (Fig. 8a-c), a finding that may be causal rather than merely correlative. Considering DE- and DAS-regulated genes separately revealed that the average intron number was much lower in DE than in DAS genes, while ASS appears to be biased towards higher values (i.e., stronger splice sites) in this group relative to DAS-regulated genes. With few exceptions, genes with intron numbers >20 were regulated by DAS (Fig. 8b). However, highly negative ASS values do not appear to compromise transcriptional regulation, indicating that DE can be promiscuously employed to regulate genes hardwired for being DAS-regulated. Factors that may support transcriptional regulation such as promoter length and the number of transcription factor binding sites (TFBSs) did not pose major effects on gene expression. While promotor length was correlated with the number of TFBSs, no bias towards higher values was observed for DE-regulated genes, suggesting that these parameters are not decisive for the type of regulation (Fig. 8c, d). It thus appears that DE is not dictated by genomic features and can be called into play when more robust regulation is required.
The type of gene regulation is genomically hardwired across biological process
To further investigate the mechanisms underlying the regulation of Fe-responsive genes, we compared the genetic architecture of genes involved in various processes that are affected by the Fe regime. Homing in on individual genes of the ferrome, a comparison of the DE- and DAS-regulated subsets supported the trend observed for all Fe-responsive genes, i.e., a marked shift towards stronger splice sites for DE-regulated genes (Fig. 9). However, regulation by DE was also observed for genes with very weak (i.e., highly negative) splice sites and high intron numbers (Fig. 9a). The difference in ASS between DE- and DAS- regulated genes was obvious at all experimental time points, and was independent of the plant part under study (i.e., roots vs. shoots) and the direction of regulation (Supplementary Fig. S3).
In addition to the ferrome, Fe-responsive genes involved in glycolysis (75 out of a total of 133), the TCA cycle (34/63), JA-signalling (25/36), transcriptional regulation (432/1717), and pre-mRNA splicing (163/396) were analysed. This comparison revealed pronounced differences in splice site strength and intron number among the various groups that appear to be typical of a particular process. When compared with all Fe-responsive entities (n = 9,190), genes from the ferrome group, JA-related genes, and genes encoding transcription factors showed a lower-than-average number of introns and a higher-than-average ASS, while genes related to glycolysis, the TCA cycle, and splicing-related genes exhibited an opposite pattern (Fig. 10).
Separating exclusively DE- or DAS-regulated genes revealed pronounced differences between the two groups, with DE-regulated genes exhibiting much stronger splice sites and lower intron numbers (Fig. 8b). Strikingly, DAS-regulated transcription factors harbour features of DE-controlled genes, possibly reflecting an evolutionary trend of this group towards this type of regulation. It is further interesting to note that – while containing transcription factors and metabolic genes to almost equal portions – ferrome genes behaved similar to genes encoding transcription factors, suggesting that the massive changes in abundance in response to environmental cues necessitates this type of regulation in most genes in these two groups to allow for an adequate and efficient response to environmental cues.
Although promoter length and the number of TFBSs do not seem to affect the mode of gene regulation, these parameters can be supportive to or even critical for gene expression. For example, in yeast, promoters of stress-responsive genes were found to be longer than those of housekeeping genes51, suggesting that such responsiveness requires a more elaborate interaction between trans-acting factors and cis-regulatory elements. In the present study, the average promoter length (1,595 bp) and the average number of TFBSs (33.5) of all 9,190 Fe-responsive genes were not significantly different from those determined for genes of the categories ferrome, glycolysis, and TCA cycle, but slightly higher (38.1) for genes encoding transcription factors (Fig. 9). Crucially, splicing-related genes had on average significantly shorter promotors (1,290 bp) and less TFBSs (27.5), indicating a trend against regulation by DE of these genes. Another deviation from the overall average was observed for JA-related genes, which had longer promoters and more TFBSs (Fig. 11). The promoter architecture of JA-related genes may reflect the extent of recruitment by different signalling pathways. For example, genes encoding enzymes mediating the first steps of JA synthesis (i.e., LOX2, LOX3, LOX4, AOS, and AOC3) had an average promotor length of 2,764 bp and 63.6 TFBSs, almost twice of the average of Fe-responsive genes. JA synthesis is important for a plenitude of processes, which is possibly mirrored by a large number of regulatory cis-consensus motifs on the promotors of these genes. By contrast, clade Ib bHLH proteins, which are highly responsive to the Fe regime but have not been associated with other responses, have short promotors (on average 918 bp) and relatively few TFBSs (24.5). Collectively, the data suggest that a particular process is prone to a certain type of gene regulation, which is, by an appreciable extent, hardwired by the structure of the genes involved in this process.