Deciphering Mineral Homeostasis in Barley Seed Transfer Cells at Transcriptional Level.

In addition to the micronutrient inadequacy of staple crops for optimal human nutrition, a global downtrend in crop-quality has emerged from intensive breeding for yield. This trend will be aggravated by elevated levels of the greenhouse gas carbon dioxide. Therefore, crop biofortification is inevitable to ensure a sustainable supply of minerals to the large part of human population who is dietary dependent on staple crops. This requires a thorough understanding of plant-mineral interactions due to the complexity of mineral homeostasis. Employing RNA sequencing, we here communicate transfer cell specific effects of excess iron and zinc during grain filling in our model crop plant barley. Responding to alterations in mineral contents, we found a long range of different genes and transcripts. Among them, it is worth to highlight the auxin and ethylene signaling factors Arfs, Abcbs, Cand1, Hps4, Hac1, Ecr1, and Ctr1, diurnal fluctuation components Sdg2, Imb1, Lip1, and PhyC, retroelements, sulfur homeostasis components Amp1, Hmt3, Eil3, and Vip1, mineral trafficking components Med16, Cnnm4, Aha2, Clpc1, and Pcbps, and vacuole organization factors Ymr155W, RabG3F, Vps4, and Cbl3. Our analysis introduces new interactors and signifies a broad spectrum of regulatory levels from chromatin remodeling to intracellular protein sorting mechanisms active in the plant mineral homeostasis. The results highlight the importance of storage proteins in metal ion toxicity-resistance and chelation. Interestingly, the protein sorting and recycling factors Exoc7, Cdc1, Sec23A, and Rab11A contributed to the response as well as the polar distributors of metal-transporters ensuring the directional flow of minerals. Alternative isoform switching was found important for plant adaptation and occurred among transcripts coding for identical proteins as well as transcripts coding for protein isoforms. We also identified differences in the alternative-isoform preference between the treatments, indicating metal-affinity shifts among isoforms of metal transporters. Most important, we found the zinc treatment to impair both photosynthesis and respiration. A wide range of transcriptional changes including stress-related genes and negative feedback loops emphasize the importance to withhold mineral contents below certain cellular levels which otherwise might lead to agronomical impeding side-effects. By illustrating new mechanisms, genes, and transcripts, this report provides a solid platform towards understanding the complex network of plant mineral homeostasis.

Introduction

Recognized for long, iron and zinc are essential cellular elements of importance to public health [1, 2]. However, the human diet-dependency on staple crops which are poor in minerals provokes inadequate intake of iron and zinc. Human micronutrient deficiency is a global health challenge and is categorized among the most important disease precursors. Iron and zinc deficiencies in plants are also agricultural challenges [3, 4]. Opposing to the higher levels of plant micronutrients is not only the artificial selection for higher yield but also the daily increasing levels of atmospheric carbon dioxide [5]. To combat micronutrient deficiency, we have to boost the micronutrient levels in crop varieties. The strategy although have some severe negative feedback loops resulted by metal ion-enhanced oxidative stress [2]. Therefore, it is very decisive to fully unravel the components of mineral homeostasis. Microarray transcript profiling has been used in the search for important components of mineral homeostasis. This introduced the heavy metal ATPase Hma3, cation diffusion facilitator 1, and nicotianamine synthase Nas3 as the responsible genes for zinc hyperaccumulation in Arabidopsis halleri [6]. Another cross-species comparison introduced the zinc transporter ZIP4, metal tolerance proteins MTP1 and 8, HMA3 and 4, iron regulated transporter IRT3, and natural resistance-associated macrophage protein NRAMP3 as well as lignin biosynthesis to be involved in metal hyperaccumulation by Thlaspi caerulescens [7]. Kobayashi et al. [8] highlighted the induction of methionine cycle and phytosiderophore biosynthesis genes under iron deficiency in rice. Mineral loading components of rice seeds have also been mapped to the members of Zip, Yellow stipe-like (Ysl), and Nramp gene families [9]. A comprehensive microarray study indicated the importance of iron for plastids and Fe-S cluster biosynthesis as well as detoxification mechanisms of metals like zinc under iron deficiency [10]. Moreover, global repression in photosynthesis related genes has been emphasized in iron depleted cyanobacteria [11]. Recently, RNA-Sequencing (RNA-Seq) has been applied to dissect plant mineral homeostasis [12, 13] and most interesting, recalled the repression of photosynthesis and tetrapyrrole biosynthesis related genes in iron deficient Arabidopsis [14].

In spite of our current knowledge, plant mineral homeostasis has yet to be fully explored. Utilizing microarray, we previously mapped the expression of diverse family members including Zip, Nramp, Ysl, Hma, Mtp, and Ca2+/H+-exchanger (Cax) genes to the transfer cells, aleurone layer, endosperm, and embryo of the barley grain [15]. However, microarray studies are limited to our knowledge of previously identified genes, they do not provide transcript splice-specific information, and suffer from narrow detection dynamic-range of expressional changes. To further extend, field-grown barley plants at the height of grain filling were subjected to foliar application of iron or zinc and a transcriptome-wide expression analysis was carried out for the transfer cells. Transfer cells bridge the maternal flux of materials towards the endosperm. As the edible part of grain, endosperm has been the target for micronutrient biofortification in crops. One advantage of the current study was the system-wide function and interaction analysis of the differentially expressed genes in order to identify highly potent candidates and strategies for biofortification. By using RNA-Seq, we not only provide experimental evidence to our recent homeostasis models [2] but also introduce a whole range of new aspects and components cooperating in plant mineral homeostasis. This includes different chromatin remodeling factors, auxin transporters, auxin and ethylene signaling components, clock genes, storage proteins, retroelements, and targeted sorting of protein. The results imply the involvement of different regulatory and non-regulatory mechanisms in fast response to the excess of micronutrients and recommend the bypassing of negative feedback loops and/or the use of upstream regulatory factors for biofortification crops.

Results and Discussion

Data analysis and statistics

Transcriptome sequencing was performed to analyze the cellular response to foliar applications of iron and zinc in transfer cells of immature barley-grains. On average, 44.5% of the reads were mapped onto the barley genome (S1 Table). More than 50% of the 17.2 million exonic single-hit mapped reads evidenced exon-exon junctions per sample-replicate (Fig 1 and S1 Table). The reproducibility measure of slope of inter-replicate regression lines was improved from 0.86 to 0.95 after data correction. Data reproducibility was further checked using a sensitive measure of deviation from perfect reproducibility [16]. The averages of ratios' deviations improved from 55% in uncorrected data to 44% in corrected data. Reference gene-based correction showed no advantage compared to the total exon read count-based correction (Fig 2), revealed by bias structure analysis as introduced in [17]. This is expected from analyzing single cell-type samples.

Fig 1

Mapping statistics.

On average, 17 million exonic-reads were mapped per replicate. Fe and Zn represent iron and zinc treatments. Samples collected at 6 h and 24 h after the treatments are labeled by 6 and 24. UT is for untreated sample. R1, R2, and R3 represent biological replicates.

Fig 2

Bias structure analysis.

The reference gene based-correction showed no advantage over the genome-wide correction. This can be explained by the single cell-type origin of the samples. We selected 18 reference genes with coefficient variations below 10% after analyzing 43 publicly-available barley microarray datasets including 891 samples and representing 22840 probes (S1 File). The inter-replicate variation and inter-treatment bias were calculated using the selected genes as described previously [17]. We considered all possible inter-treatment comparisons in our RNA-Seq data. There was no difference in inter-treatment bias among the single reference gene based-correction, exon read count based-correction measured as RPM, and uncorrected data. Sequential application of the corrections [17] was also not efficient. However, the exon read count based-correction was used for data analysis due to the large proportion of expelled inter-replicate variations. RPM: Read count Per Million Mapped reads, Std.: Standard deviation.

On average, 17 million exonic-reads were mapped per replicate. Fe and pan class="Chemical">Zn represent iron and zinc treatments. Samples collected at 6 h and 24 h after the treatments are labeled by 6 and 24. UT is for untreated sample. R1, R2, and R3 represent biological replicates.

Statistical comparisons were carried out at both gene and transcript levels. On average, more than 90% of the differentially expressed genes and transcripts showed no significant changes at their transcript and gene levels, respectively. While this can be explained by switching, cumulative effects, and independent regulation of alternative transcript isoforms, it also reveals a complementary advantage of performing comparisons at both levels. Furthermore, we examined the expression changes of 30 differential expression events including genes and transcripts using Real-Time PCR. All changes, except one, were confirmed qualitatively (Fig 3).

Fig 3

Real-Time PCR confirmed the expression changes.

Gene- and transcript-specific primers were used to examine the expression changes of 30 selected events found by RNA-Seq. Accession numbers are shown below the gene names. The relative expression fold changes are shown individually. Except for the Cca1 transcript TCONS_00166188, the expression changes were all in agreement with the RNA-Seq results (see S1 Fig). As the most stable gene (S1 File), the vacuolar ATP synthase was used to normalize the Real-Time PCR data. We also found identical expression patterns when normalizing the data to barley Gadph (S1 Fig). Three biological replicates (1–3) of samples including untreated sample (UT), 6 h after treatments (6Fe and 6Zn), and 24 h after treatments (24Fe and 24Zn) are shown in different colors. See S2 File for the functions of genes.

Responding to the mineral supply, transcript isoform switching was a major regulapan class="Gene">tory mechanism and covered 19% of the differentially expressed transcripts (Table 1). Isoform switching in 18 genes (S3 File) including sucrose-phosphate synthase and indole-3-acetic acid-amido synthetase occurred among alternative transcripts coding for identical proteins. This mechanism can facilitate the utilization of one specific protein product through different antagonistic signaling pathways.

Table 1

Statistics of the comparisons.

Gene-level comparisons	Gene/Transcript	DE	Rep	Ind	IS	IS1	IS2
24Fe/UT	16250	302	166	136	2	2	0
24Zn/UT	16171	106	85	21	0	0	0
6Zn/6Fe	16476	78	26	52	1	1	0
24Zn/24Fe	16758	105	95	10	2	1	1
Isoform-level comparisons	Gene/Transcript	DE	Rep	Ind	IS	IS1	IS2
24Fe/UT	18367	249	139	110	72	63	13
24Zn/UT	18696	484	239	245	113	94	19
6Zn/6Fe	18047	410	215	195	58	48	10
24Zn/24Fe	19389	585	296	289	86	71	11

Number of significant changes (DE), repressed events (Rep), induced events (Ind), events with switched isoforms (IS), isoform switching among transcripts coding different protein isoforms (IS1), and isoform switching among transcripts only different in untranslated regions (IS2). UT: untreated samples, 6 & 24: samples of 6 h and 24 h after treatments, Fe & Zn: iron and zinc treatments.

Number of significant changes (DE), repressed events (Rep), induced events (pan class="Chemical">Ind), events with switched isoforms (IS), isoform switching among transcripts coding different protein isoforms (IS1), and isoform switching among transcripts only different in untranslated regions (IS2). UT: untreated samples, 6 & 24: samples of 6 h and 24 h after treatments, Fe & Zn: iron and zinc treatments.

The effects of iron and zinc on upstream signaling pathways

Our data addresses a key role for chromatin modifications in mineral homeostasis. The SKB1- mediated histone dimethylation interferes with iron uptake in Arabidopsis [18]. Accordingly, barley Skb1 showed higher expression in iron-treated plants compared to zinc-treated plants. Further, we found metal ion-triggered expression changes in 38 chromatin remodeling genes, providing genome-wide insights into the importance of fast chromatin-level responses to mineral availability. Comparing iron to zinc treatments, 27 of these genes were differentially expressed (S2 Table). Both regulatory and non-regulatory mineral homeostasis genes including Ysl, Fro2, Nas, Hma, bHlH39, Irt1, Ferritin, and many others follow light-dark cycles [2]. Most interesting, the histone methyltransferase Sdg2 which regulates the core clock components [19] was induced by iron and zinc (see S4 File for accession numbers). Disclosed by expression changes in the oscillator master-regulator genes (S2 Fig), SDG2 introduces a signaling pathway by which the metals can govern the plant circadian rhythm. As an ethylene-signaling negative regulator, the histone acetyltransferase Hac1 [20] was also induced by iron. Auxin and ethylene are key players in plant iron homeostasis mainly by working through the promotion of uptake mechanisms [2]. Our results revealed the suppression and enhancement of ethylene signaling by iron and zinc, respectively. After iron treatment, ethylene signaling was restricted by induction of known repressors; two hypersensitive to phosphate starvation genes Hps4 [21] and one serine/threonine-protein kinase gene Ctr1 [22] (Fig 4A). In contrast, the positive regulators Win1 and 5'-3' exoribonuclease Xrn4 [23] had higher levels of expression after zinc treatment (Fig 4A). We also found the EIN3-interacting transcription factor EER4 with higher expression in zinc-treated plants compared to iron-treated plants. The most likely scenario is the activation of iron deficiency response by the zinc repletion [see 2], underlined by the decreased iron content of seeds (Table 2). Therefore, iron-treated plants did not need rapid recovery from ethylene response which was consistent with the observed repression of EIN3-binding F-box gene Ebf1 that confers a fast feedback against ethylene [see 24].

Fig 4

Hormonal signaling pathways.

(A) Ethylene signaling. In contrast to iron-treated plants where the ethylene response was suppressed, zinc treatment triggered the signaling pathway. (B) Auxin (IAA) signaling. The treatments altered the cellular compartmentalization of auxin as well as its signaling pathways. The data reveals a very early-stage induction in auxin signaling after the treatments followed by negative feedback loops enabling auxin sequestration within organelles and cellular auxin efflux. This occurred likely due to iron surplus in the iron-treated plants as well as zinc repletion and/or continuous iron deficiency after foliar application of zinc. Accession numbers of the genes are available in S4 File. Different transcript isoforms of genes are shown by capital letters, immediately after the name of genes. Fe and Zn represent iron and zinc treatments. The samples of 6 h and 24 h after treatments are labeled by 6 and 24. UT stands for untreated plants. Comparisons of 24Fe/Untreated sample and 24Zn/Untreated sample are shown as 24Fe and 24Zn, respectively. Zinc treatment was compared with iron treatment which is shown as 6Zn/6Fe or 24Zn/24Fe. See S2 File or the text for detailed functions of genes. AA: amino acid, IBA: indole-3-butyric acid.

Table 2

Seed mineral content.

Samples	Nutrients
	Cu	Fe	Zn	Mg	Mn	P	S	K
6Zn	3.2	20	23.4	794	9	2682	759	5684
24Zn	3.7	19.3	84.2	759	10.4	2574	839	5776
6Fe	3.5	67.8	21.7	594	5	2542	706	5854
24Fe	3.7	67.7	23	790	8.2	2524	759	5032
UT	3.9	22.1	25.3	852	9.3	2762	812	6026

Total seed mineral contents (μg/g) of the samples were determined by ICP-MS. The quantities represent the average of four technical replicates with a coefficient variation of 1–2%. 6Fe: six hours after iron treatment, 6Zn: six hours after zinc treatment, 24Fe: 24 hours after iron treatment, 24Zn: 24 hours after zinc treatment, UT: untreated sample.

Total seed minerpan class="Chemical">al contents (μg/g) of the samples were determined by ICP-MS. The quantities represent the average of four technical replicates with a coefficient variation of 1–2%. 6Fe: six hours after iron treatment, 6Zn: six hours after zinc treatment, 24Fe: 24 hours after iron treatment, 24Zn: 24 hours after zinc treatment, UT: untreated sample.

Recalling the importance of auxin in mineral homeostasis [2], its signaling cascade has yet to be established. Here, we introduce different components of auxin signaling involved in mineral homeostasis (Fig 4B). In general, the elevated iron and zinc levels were in accordance with suppression of auxin signaling (Table 2 and Fig 4B). ABCB/MDR transporters facilitate the sequestration within cellular compartments or efflux across the plasma membrane (PM) of auxin [25-27]. As auxin transporters, we found a putative Golgi ABCB2 as well as the PM ABCB19 induced by the treatments (Fig 4B). Resembling zinc-mediated iron deficiency signaling, the auxin precursor exporter Abcg37 (Pdr9) [26] had higher expression level in iron-treated plants compared to zinc after 6 h. It seems, the cells activated auxin signaling soon after the treatments followed by auxin sequestration within cellular compartments or efflux across PM. In parallel with the increased auxin efflux from cytosol, Cand1 which is involved in AUX degradation [28] was repressed by iron. AUX proteins are short-lived repressors of early auxin-response [29]. The auxin signaling repressor ECR1 [30] was also found to work in parallel after zinc treatment (Fig 4B). Accordingly, the auxin-response factors Arf12-like, Arf4, Arf25, and Arf9 were all repressed. The mutant arf12/25 plants show lower cellular iron, zinc, and auxin contents [31].

In addition to the auxin efflux and distribution, the reversible catpan class="Chemical">alysis of less active auxin-amino acid (IAA-AA) conjugates [32] participates in mineral homeostasis. The influx of metals like zinc into the ER has been proposed to inhibit the IAA-leucine resistance 1 (ILR1)-mediated IAA-AA hydrolysis which can be activated by IAR1-facilitated metal efflux [32-35]. The expressional adjustments of a putative ER ZIP protein coding gene Iar1 likely released metals copper and zinc for vital biological processes in iron-treated cells and suppressed IAA-AA hydrolysis after zinc treatment (Fig 4B). A putative cytoplasmic IAA-amido synthetase GH3.5 (Jar1) switched transcript isoforms coding for an identical amino acid sequence when comparing iron and zinc treatments. Also, two probable chloroplast- and ER-localized IAA-amido synthetase IaasGH3.2 coding genes [36] were repressed which underlines decreased auxin levels in iron-treated plants.

Iron deficiency promotes the formation of branched root hairs [37, 38] and laterpan class="Chemical">al root elongation is triggered after iron supply [39]. In this context, we noticed root growth and development factors Nac1 [40], the fusion gene nodulin/glutamine synthase [41], the auxin-induced in root cultures Air12 [42], and the root hair defective Rhd3 [43]. Their expression changes imply the enhancement of iron uptake/deficiency signaling and root elongation mechanism at early and later stages after iron treatment, respectively. Most likely, these genes are upstream factors responding to mineral fluctuations in different tissues.

Excess zinc and iron trigger stress-response factors and affect photosynthesis and respiration

More than 20% of the responding genes were stress related, including DNA damage/repair, programmed cell death and autophagy, retroelements, and members of storage proteins. Difpan class="Chemical">ferent upstream stress response genes were employed to decrease ROS production and ROS-mediated cellular damage (S3 Fig and S4 Table, see S1 Text for functional details). The induction of 2Fe-2S binding domain-free coding transcripts in two xanthine dehydrogenase genes may be an adaptive strategy against iron deficiency after zinc treatment. These isoforms should also be resistant to conversion from the dehydrogenase to the oxidase form which produces superoxide and hydrogen peroxide [see 44].

Surprisingly, we found seed storage protein encoding genes induced by the treatments (Fig 5). Compared to the average amino acid composition of plant proteins, we noticed 750%, 350%, and 50% higher glutamine, proline, and cysteine frequencies in the upregulated storage proteins analyzed by the CLC software. Hydrolysis of hordeins and partial deamidation of glutamine residues leads to small peptides carrying ionized carboxyl groups with radical scavenging and metal chelation activities [45]. Therefore, we hypothesize a function for the induced storage proteins in both metal chelation and metal toxicity resistance. Considering the function of retroelements in plant adaptation to environmental cues like metal depletion [46, 47], we also noticed a fast activation of retroelements in response to metal surplus (Fig 6).

Fig 5

Differential expression of seed storage proteins.

Genes encoding seed storage proteins including B, C, D, and G hordeins (H), were upregulated by the treatments. See S2 File for the accession numbers and gene names. 6Fe: 6 h after iron treatment, 6Zn: 6 h after zinc treatment, 24Fe: 24 h after iron treatment, 24Zn: 24 h after zinc treatment, UT: untreated plants. For example, 24Fe/UT represents the comparison of 24Fe with UT.

Fig 6

Retroelements were activated by foliar applications of iron and zinc.

See S3 Table for the accession numbers and other details of events shown by numbers of 1–34. 6Fe: 6 h after iron treatment, 6Zn: 6 h after zinc treatment, 24Fe: 24 h after iron treatment, 24Zn: 24 h after zinc treatment, UT: untreated plants. For example, 24Fe/UT represents the comparison of 24Fe with UT.

Genes encoding seed storage proteins including B, C, D, and G hordeins (H), were upregulated by the treatments. See S2 File for the accession numbers and gene names. 6pan class="Chemical">Fe: 6 h after iron treatment, 6Zn: 6 h after zinc treatment, 24Fe: 24 h after iron treatment, 24Zn: 24 h after zinc treatment, UT: untreated plants. For example, 24Fe/UT represents the comparison of 24Fe with UT.

See S3 Table for the accession numbers and other details of events shown by numbers of 1–34. 6Fe: 6 h after pan class="Chemical">iron treatment, 6Zn: 6 h after zinc treatment, 24Fe: 24 h after iron treatment, 24Zn: 24 h after zinc treatment, UT: untreated plants. For example, 24Fe/UT represents the comparison of 24Fe with UT.

Our results also emphasize that photosynthesis and respipan class="Species">ration were affected by zinc-mediated iron deficiency (S2 Text and S5 Fig). The data provides additional evidence to the interplay between iron and oxygen during the evolution as discussed in [2] by indicating a limited aerobic metabolism, i.e., suppression of Krebs cycle and electron transfer chains (S2 Text), to ensure less oxidative stress under excess iron levels. Taken together, the described negative side-effects of excess iron and zinc have implications for crop biofortification. Not only is it costly for the cells and potentially imposes a yield penalty, but the elevated cellular metal contents also trigger negative feedbacks, which work against the desired mineral levels and prioritize a minimum buildup of oxidative stress in biofortified crops.

The S-adenosylmethionine cycle and mineral homeostasis

Nicotianamine is a key player in pan class="Chemical">metal transport and chelation [48]. Most upstream in the nicotianamine biosynthetic pathway, glutamate carboxypeptidase 2 (AMP1) releases folate monoglutamate from folate polyglutamate [49]. The amp1 mutants are insensitive to nitric oxide [50]. Therefore, repression of Amp1 can be explained as a negative feedback against nitric oxide-triggered iron uptake under continuous iron deficiency after 24 h of zinc treatment (Table 2 and Fig 7). In contrast, iron treatment enhanced the transcription of Amp1 to promote S-adenosylmethionine (SAM) biosynthesis. The repression of chloroplast folate-biopterin transporter Fbt1 [51] and methionine synthase Ms1 [52] and the induction of homocysteine S-methyltransferase Hmt3 [53, 54] could also trigger the production of nicothianamine precursor “SAM” (Fig 7). The gene coding for SAM decarboxylase (AdoMetDC), which utilizes SAM for biosynthesis of polyamines [55], was also repressed by iron. Therefore, SAM was preferably channeled to the nicotianamine biosynthesis after iron treatment. In contrast, zinc treatment repressed the nicothianamine and polyamine biosynthesis pathways in favor of glutathione prouction (Fig 7).

Fig 7

SAM plays an important role in mineral homeostasis.

SAM biosynthesis was enhanced by iron. HMT catalyzes SAM biosynthesis by methyl transfer to homocysteine [53, 54]. The suppression of PA biosynthesis could results in a large pool of substrate “SAM” in favor of NA biosynthesis after iron treatment. Not only PA biosynthesis but also NA biosynthesis was suppressed by zinc. Concurrently, glutathione synthetase (Gsh-s) was induced. Accession numbers of the genes are available in S4 File. Fe and Zn represent iron and zinc treatments. The samples of 6 h and 24 h after treatments are represented by 6 and 24. UT stands for untreated plants. Comparisons of 24Fe/Untreated sample and 24Zn/Untreated sample are shown as 24Fe and 24Zn, respectively. Zinc treatment was compared with iron treatment which is shown as 6Zn/6Fe or 24Zn/24Fe. HMT: homocysteine S-methyltransferase, F-MGlu: folate monoglutamate, F-PolyGlu: folate polyglutamate, NA: nicotianamine, PA: polyamine, SAH: S-adenosylhomocysteine, SAM: S-adenosylmethionine.

SAM biosynthesis was enhanced by pan class="Chemical">iron. HMT catalyzes SAM biosynthesis by methyl transfer to homocysteine [53, 54]. The suppression of PA biosynthesis could results in a large pool of substrate “SAM” in favor of NA biosynthesis after iron treatment. Not only PA biosynthesis but also NA biosynthesis was suppressed by zinc. Concurrently, glutathione synthetase (Gsh-s) was induced. Accession numbers of the genes are available in S4 File. Fe and Zn represent iron and zinc treatments. The samples of 6 h and 24 h after treatments are represented by 6 and 24. UT stands for untreated plants. Comparisons of 24Fe/Untreated sample and 24Zn/Untreated sample are shown as 24Fe and 24Zn, respectively. Zinc treatment was compared with iron treatment which is shown as 6Zn/6Fe or 24Zn/24Fe. HMT: homocysteine S-methyltransferase, F-MGlu: folate monoglutamate, F-PolyGlu: folate polyglutamate, NA: nicotianamine, PA: polyamine, SAH: S-adenosylhomocysteine, SAM: S-adenosylmethionine.

Ethylene-insensitive3 like 3 (pan class="Gene">EIL3) may co-regulate cellular metal and SAM reserves. EIL3 is involved in high-affinity sulfate uptake [54] and can be induced by iron deficiency [56]. Eil3 accordingly showed higher expression in iron-treated plants compared to zinc-treated plants after 6 h, perhaps to meet the lower seed sulfur content of iron-treated plants (Table 2). This was followed by the enhanced SAM production in iron-treated plants after 24 h. We also found the transcription factor Vip1 repressed in the zinc-treated plants. Arabidopsis plants defective in VIP1 have higher contents of sulfur, glutathione and cysteine [57]. Accordingly, the seed sulfur content was elevated by the zinc treatment (Table 2).

Alterations in metal-trafficking components

Members belonging to a diverse set of metal transporter families were in action against the excess minerpan class="Chemical">al supply (Fig 8). The treatments compelled the cells to restrain their uptake routes of iron or zinc. Vacuolar sequestration was also activated. Concurrently, uptake and internal-source utilization of other metals were enhanced. Indicating the plasticity of plant mineral homeostasis in terms of regulation and negative feedbacks to the metal repletion, this sheds light on the potential usage of upstream regulatory components for crop biofortification strategies.

Fig 8

Expression changes in metal transporters.

The treatments compelled the cells to restrain the iron or zinc uptake routes across the plasma membrane and chloroplast envelope and to sequester the minerals within vacuoles. In parallel, uptake and internal-source utilization of other metals were enhanced. This indicates a strong dependency of cellular levels of minerals on each other and sheds light on the negative feedback loops to avoid excess mineral contents. Accession numbers of the genes are available in S4 File. Different transcript isoforms of genes are shown by capital letters, immediately after the name of genes. Fe and Zn represent iron and zinc treatments. The samples of 6 h and 24 h after treatments are represented by 6 and 24. UT stands for untreated plants. Comparisons of 24Fe/Untreated sample and 24Zn/Untreated sample are shown as 24Fe and 24Zn, respectively. Zinc treatment was compared with iron treatment which is shown as 6Zn/6Fe or 24Zn/24Fe. DMA: deoxymugineic acid, GS: glutathione, NA: nicotianamine.

The Arabidopsis proton pumps, Aha2 is working in iron uptake from the rhizosphere [58]. Here, we report the barley Aha2 gene repressed in the transfer cells by both zinc and iron treatments (Fig 8). This introduces the apoplastic alkalinization as a widespread mechanism in plants which serves as the first barrier to lessen the uptake under metal repletion. As negative feedback loops, furthermore, a coordinated action of metal uptake and storage mechanisms were opposed to the enhanced cellular metal contents. HMA transporters are elaborated zinc and copper transporters [59-63]. The PM Hma2 was repressed by zinc and iron which supports the rice HMA2 influx activity reported in [64]. It needs further evidence on whether HMA2 can transport iron or the iron-mediated suppression was an inevitable regulatory response. We also found the vacuolar zinc transporter of Cation Diffusion Facilitator Mtp1 [65] with higher expression in plants treated by zinc compared to iron after 6 h (Fig 8).

We noticed different members of the ABC transporters influenced by the treatments. In pan class="Species">Arabidopsis, Mrp3 can be induced by cadmium, nickel, arsenic, cobalt and lead but not by zinc and iron [66]. Most interestingly, our transcript-level comparisons revealed the contribution of Abcc3 (Mrp3) to the iron homeostasis (Fig 8). The higher expression of Abcc8-like in zinc-treated plants compared to iron-treated plants was likely to sequester zinc within vacuoles. Its homologue in Arabidopsis, MRP6, is involved in cadmium detoxification [67]. Recalling the glutathione S-Mn/Cu conjugate transport activity of Arabidopsis ABCC1 and ABCC2 [68, 69], the expression changes of Abcc2 (Mrp2, Fig 8) explain the additional supply of manganese and copper to recover the distressed photosynthesis (consult S2 Text and Table 2). Furthermore, the protease CLPC1 which is a key factor to supply chloroplasts with iron [70] had higher expression after zinc treatment compared to iron (Fig 8). Finally, the repression of Abcb27 as the homologue of aluminum sensitive ALS1 in rice [71], was conceivably to hinder zinc-overloading of chloroplasts (Fig 8).

As a subgroup of the oligopeptide transporter family (OPT), YSL proteins are metal transporters [72-76]. We found a putative ER/vapan class="Chemical">cuolar YSL15 coding gene induced by iron. Its homologue in rice, OsYSL15, transports Fe (III)-deoxymugineic acid and displays inward transport in yeast experiments [74, 75] making vacuolar sequestration of iron unlikely. It seems HvYSL15 releases metals out of the ER in parallel with Iar1 (Figs 4B and 8). The Chloroplast/PM YSL16 coding gene was repressed by both zinc and iron. Rice YSL16 imports Cu-nicotinamine complex across PM [77]. Therefore, this response is likely to preserve the chloroplast copper. The grain copper content was decreased by the treatments (Table 2).

The Arabidopsis pan class="Gene">NRAMP3 and 4 export iron and manganese out of the vacuoles [78, 79]. We found a putative PM/tonoplast NRAMP3/4-like isoform repressed by iron treatment. Another transcript isoform of this gene had higher expression level in iron-treated plants compared to zinc. While Thlaspi caerulescens NRAMP4 transports iron and zinc, TcNRAMP3 only transports iron [80]. In contrast, AtNRAMP2 does not transport iron [81]. Accordingly, a putative tonoplast NRAMP2 displayed lower expression by zinc compared to iron. This gene likely releases the vacuolar zinc. As we noticed after zinc application, iron deficiency downregulates the vacuolar iron importer nodulin-like 2 in Arabidopsis [82]. The nodulin-like 1 had accordingly lower expression in zinc-treated plants compared to iron (Fig 8). Furthermore, a putative PM zinc transporter 2 (Zip2) was repressed by iron (Fig 8). Recalling the zinc and manganese transport activities of AtZIP2 [83, 84], the response is supported by the reduced manganese levels after iron treatment (Table 2). It seems HvZIP2 is not very specific and translocates other metals like iron as well. Additionally, the barley PM ZIP7 [85] had higher expression in iron-treated plants compared to zinc after 24 h. As for Nramp2, we noticed alternative transcripts of Zip7, encoding for small-size proteins. Whether these microproteins contribute to the transport activity or act in negative regulation of NRAMP2 and ZIP7 need further experiments.

Among the toxic metal transporters of cyclic nucleotide-gated channels [86-88], the Cngc1 and Cngc14 were differentially expressed. Furthermore, two dentin sialophosphoprotein-like genes had higher expression levels in iron-treated plants compared to zinc after 6 h. In human, the distribution of dentin sialophosphoprotein regulates mineral deposition and when mutated, it resembles the abnormal tooth biomineralization by defective cyclin M4 (cnnm4) [89]. The PM ancient conserved domain protein CNNM4 has magnesium efflux activity in mice [90]. Accordingly, iron induced putative tonoplast and Golgi/ER CNNM4 proteins with CBS domain (Fig 8). The CBS domain is known from bacterial CorC and CorB magnesium/cobalt efflux transporters [91, 92]. The suppression of vacuolar magnesium storage in zinc-treated plants after 24 h was in agreement with the repression of Cbl3. The tonoplast CBL3-CIPK23 calcium signaling promotes vacuolar magnesium sequestration in Arabidopsis [93]. This emphasizes that the vacuolar-function is central in the metal and calcium cross homeostasis proposed in [2]. Shown in Fig 8, vacuoles are playing crucial role in mineral homeostasis by storing and releasing of micronutrients in repleted and depleted cells, respectively. Accordingly, the membrane protein YMR155W-like coding gene was downregulated by iron. When mutated, it affects vacuolar membrane fragmentation and shows metal resistance phenotype in yeast [94, 95]. Therefore, the repression of Ymr155W likely enhances metal sequestration through vacuolar organization. We also found the vacuolar protein sorting 4 with higher expression in iron-treated plants compared to zinc. Arabidopsis VPS4 is involved in the protein trafficking and maintenance of large central vacuole [96]. The same expression pattern was found for the gene RabG3F that together with monensin sensitivity 1(Mon1) is involved in prevacuolar compartments-to-vacuole trafficking and vacuole biogenesis [97]. MON1 is also involved in iron homeostasis of animal cells [98, 99].

Human RNA poly(rC)-bpan class="Chemical">inding proteins (PCBPs), as iron chaperones, are involved in delivering iron to non-heme iron containing proteins [100, 101]. Here, we report expression changes in PCBPs and thereby their potential role in plant mineral homeostasis not only by iron chaperoning but also through post-transcriptional regulation [see 102]. Recently, the mediator subunit 16 has been identified as a very upstream positive regulator of iron uptake in Arabidopsis [103]. Here, we report a sharp repression for Med16 after iron treatment implying a common regulator of iron homeostasis in graminaceous and non-graminaceous plants. Additionally, an iron-induced hemerythrin-like gene possibly worked as the repressor of iron uptake factors. The hemerythrin-like FBXL5 senses iron and oxygen availability and thereby negatively regulates iron absorption in mammalian cells [104-106]. Finally, we noticed the serine/threonine-protein kinase TOR and replication protein RPA which recently have been introduced as iron homeostasis components [107]. Rpa2, a homologue of GmRpa2A, was repressed after both of the treatments. Two other Rpa2 genes were induced by zinc. Taken together, the adaptability of cellular mineral homeostasis to avoid excess metal levels highlights the advantage of regulatory factors over the downstream transporters in crop biofortification.

Intracellular protein-sorting regulates metal homeostasis

Our analysis also showed the importance of endocytosis based regulation in mineral homeostasis. Auxilin-related protein 2 which disturbs endocytosis and aluminum loading through uncoating of clathrin-coated vesicles in Arabidopsis [108], showed higher expression by iron compared to zinc, and likely halts iron uptake. The myosin-Vb gene had similar expression pattern. Together with RAB11A, animal myosin Vb retards recycling of the transferrin receptor to PM and affects iron import [109]. Responsible for endocytosis of animal transferrin receptor [110], the epidermal growth factor receptor substrate 15-like was properly repressed by iron. Conceivably, a similar mechanism functions in planta against and in favor of iron uptake after 24 h of iron and zinc treatments, respectively. In addition to the degradation-by-internalization of iron transporter IRT1 [111], the internalization of rice manganese transporter NRAMP3 has been reported under high manganese concentrations [112]. Accordingly, two protein transport Sec23A-like genes were induced by the treatments. The zinc-induced Sec23A had lower expression in zinc-treated samples compared to iron after 6 h. Yeast SEC23A is involved in internalization of PM proteins [113]. Therefore, plant SEC23A proteins are good candidates for disassembling of PM metal transporters. Here, we also highlight the metallophosphoesterase gene MPP_Cdc1_like which positively regulates manganese uptake in yeast [114]. It was upregulated after 24 h of zinc treatment. This is in accordance with the increased manganese content seen after 24 h of zinc application (Table 2). Compared to the zinc-treated plants, two other transcripts of MPP_Cdc1_like had higher expression levels in iron-treated plants after 6 h. Considering both the yeast cdc1 mutant complementation by mutations defective in vacuolar sorting of ubiquitinated proteins [114] and the CDC1 facilitated ER-to-PM/cell wall protein transport in yeast and human [115, 116], we propose that the gene MPP_Cdc1_like triggers the docking/integration of PM metal transporters like NRAMP3. We also highlight the homologue of AtExo70E1 (Exoc7) which was induced by iron. It may facilitate polar accumulation of transporters in the PM [117]. In Arabidopsis, overexpression of the FYVE1 accumulates PM IRT1 in apolar mode impairing the rate of metal uptake [118]. Phosphatidylinositol-3 and 4-phosphate are also involved in intracellular polarized and apolar trafficking of metabolites and proteins [119]. Not only the phosphatidylinositol-3-phosphate binding WD repeat and FYVE domain-containing protein 3 was induced in zinc-treated plants but also the induction of genes involved in the phosphatidylinositol-3-phosphate biosynthesis and recycling was observed (S5 Table). Moreover, the genes responsible for its usage were repressed. Therefore, higher phosphatidylinositol-3-phosphate seems to be a key trafficking factor under zinc repletion.

Conclusions

Cellular shortage and surplus of metal are both hindering cellular metabolism due to the diverse structural/catalytic roles and the enhancement of oxidative stress, respectively. We observed a high degree of complexity in iron and zinc homeostasis which seems to be the aftermath of the engagement of metals at a vast numbers of metabolic pathways. Therefore, genetic engineering for sustainably improved crop fitness against micronutrient depletion and repletion are laborious goals. The cellular impact of excess micronutrients including changes in hormonal and circadian clock signaling, affected photosynthesis and respiration, negative feedback loops and triggered oxidative stress are not trivial in terms of affecting the plant fitness. Tissue specific genetic manipulation can leave out some of the complexities. However, this can suffer from the weak tissue-specific gene promoters which we personally have realized by independent transgenic experiments conducted in our laboratory. Moreover, all the regulatory mechanisms acting on the genes, transcripts, and proteins of interest should be taken into account in the efforts to bypass the negative feedback loops. Here, we highlight the docking and recycling mechanisms and uneven distribution of mineral transporters in plasma membrane. Overexpression of influx and efflux metal-transporters in a polar localized fashion can be very efficient to enhance the directional flow of minerals. Specific targeting of the upstream regulatory genes can be another strategy. However, all the scenarios demand a high resolution map of metal homeostasis. Here, our attempt furnishes the current knowledge with the vast reaches of transcriptomic maps in iron and zinc homeostasis of barley transfer cells. This report not only covers the known grass and non-graminaceous mineral homeostasis components, but also introduces new interactors into the plant mineral homeostasis. This includes different adaptation mechanisms plant cells use against alterations in mineral availability.

Materials and Methods

Experimental design

Field-grown barley (Quench cv.) plants at the growth stage 18±2 days after anthesis were treated by foliar application of 15 mM pan class="Chemical">FeSO4.7H2O or 0.5% of ZnSO4.7H2O. Representing the control (UT) and 24 h after treatments (24Fe and 24Zn), immature seeds were collected from untreated and treated plants simultaneously at 8:30 a.m. Samples of 6 h after treatments (6Fe and 6Zn) were collected the day before at 2:30 p.m. We did not compare 6 h after treatments with either of untreated plants or samples collected 24 h after the treatments due to the possibility of interaction between mineral and diurnal-based changes. All four preferred comparisons (24Fe/UT, 24Zn/UT, 24Zn/24Fe, and 6Zn/6UT) were performed on transcript data as well as gene quantities. Three spikes from individual plants were collected representing three biological replicates. Samples were stored at -80°C before RNA extraction. Pooled RNAs extracted from three seeds representing top, middle, and bottom part of each spike was applied as one biological replicate. Finally, a single Illumina sequencing run including three channels each containing one of the three replicates from all the five samples was applied.

Total RNA extraction and sequencing

Laser capture microdissection [15] was used to collect transfer cells from 70 μm grain cryosections. Twenty sections from each seed and 60 in total were laser-bombarded for each of the biological replicates. Arcturus PicoPure RNA isolation kit (cat no. 12204–01) was used to isolate total RNA. RNA quantity and quality was determined by Agilent 2100 Bioanalyzer. For Real-Time PCR, 4 μl of the RNA samples were kept and the rest (≈ 28 μl) were used in sequencing. Due to the low quantity of samples ranging from 3.6 to 18 ng (S7 Table), linear amplification of RNAs was performed using NuGEN Ovation® RNA-Seq kit (Fig 9). In parallel with the amplification, the semi-random primers of the kit also enriched the RNA population of samples against ribosomal RNA. Paired end 2 × 101 bp sequencing was performed at AROS Applied Biotechnology A/S (Aarhus, Denmark) using HiSeq 2000 Illumina sequencing platform and Truseq technology. CASAVA 1.8.2 was used to convert sequencing results into fastq format. The data are available under the GenBank accession number of SRA297575. The ABI Applied Biosystems 7900 HT Thermal Cycler used for Real-Time PCR (see S5 File for primers). The random-primed first-strand cDNAs were synthesized using Superscript II (Invitrogen Life Technologies, Carlsbad, CA, USA) according to the manufacturer’s instructions.

Fig 9

Library preparation for sequencing.

Double-stranded cDNA profile after in-vitro linear amplification, fragmentation, and adaptor ligation is shown. Fragmentation enriched samples for ≈ 220 base-long fragments. Ligation of adaptors at both ends of the fragments changed the average length to ≈ 330 base pairs.

RNA-Seq data analysis

The quality of the sequencing reads were checked by FastQC version 0.10.1. Based on the overall high quality of the sequencing reads (S6A Fig), we did no trimming except the TopHat automatic filtering. Unmapped reads were then trimmed using Trimmomatic-0.30 [120]. The 5ʹ and 3ʹ end of reads were hotspots of Kmers and imbalanced nucleotide content (S6B and S6C Fig). Therefore, we trimmed for the first 30 and the last 11 nucleotides of the reads. In addition, a quality filtering was applied for both ends of the reads with nucleotides showing lower than 23 quality score in Phred-33 scoring system. Additionally, we performed a sliding trimming with a window size of four and average quality score of 29. The minimum acceptable read-length was 30 bases. The trimmed reads were mapped and unmapped reads were further trimmed for the 31–35th and 80–90th nucleotide positions of original reads. Second-round trimmed reads were mapped and the remaining unmapped reads were excluded. Bowtie2-2.1.0 [121] was used to index the barley genome release 18 (EnsemblPlants: http://plants.ensembl.org/info/data/ftp/index.html). Tophat-2.0.8b performed a gapped alignment to predict exon-exon junctions and detect alternative spliced isoforms [122]. To map the original reads, the inner distance between pairs was defined to 40 bases with standard deviation of 141 bases. These could cover from a complete match to 181 bases distance between pairs. The mean length of the sequenced fragments was ≈ 300–330 with inner distance of 100–120 bases between pairs (Fig 9). We allowed three and two mismatches per original and trimmed reads and defined the minimum and maximum intron length of 35 and 80000 bases, respectively. Single-hit-report option was also activated. To enrich the gene model annotations of barley genome release 18, the detected splice-sites after each mapping step were used for the next mapping. After the three rounds of mapping, the cuffmerge utility of Cufflinks-2.1.1 was used to include all the known and identified novel transcripts in the analysis by building a new genome annotation file [see 123]. The Tophat outputs for the three rounds of mapping were merged for each sample-replicate using Samtools-0.1.19 [124]. Finally, the merged new annotation file and merged accepted hits were applied to count the expression quantities using SAMMate-2.7.4 [125].

Statistical analysis

The inter-treatment comparisons were carried out using the edgeR [126]. Global dispersion was applied and the significant expression changes were selected with minimum two-fold changes using the corrected p-vpan class="Chemical">alue cut-off of 0.01. Multiple testing correction was applied to correct the p-values [127]. Considering the outlier values in weak expression events, a low quantity filtering was performed before statistical analysis. Events were not excluded when the average quantity in each comparison (two samples) and the averages of three possible two-replicate combinations, at least in one of samples, were ≥ 5 reads per million.

Functional annotation

We performed a manual functionpan class="Chemical">al annotation of novel transcripts and checked the previously annotated events. This was limited to the events which statistically showed significant changes in the comparisons. We extracted the functional information of genes from Ensemblplants (at http://plants.ensembl.org/index.html) and then, we carried out an extensive homology search for the novel and also annotated transcripts at both nucleotide and protein levels. We utilized different data bases including NCBI (http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastHome), TAIR (http://www.arabidopsis.org/), Rice database (http://rice.plantbiology.msu.edu/), and UniProtKB/Swiss-Prot (http://www.uniprot.org/uniprot/). To assign a function, we followed the minimum similarities of 90% at the protein level and 80% at the nucleotide level, and an E value below 10−30. With the ambition to use the most recent functional studies, we did an extensive literature search to get a more complete picture of the function of genes. Subcellular localization of proteins was predicted by the PlantLoc server [128]. All the data are available in S2 File.

Real-Time PCR results.

(PDF)

Click here for additional data file.

Differentially expressed components of circadian rhythm.

(PDF)

Click here for additional data file.

Differentially expressed upstream regulatory genes involved in cellular stress response.

(PDF)

Click here for additional data file.

Differentially expressed genes and transcripts involved in protein folding and degradation.

(PDF)

Click here for additional data file.

Chloroplastic tetrapyrrole, heme, and chlorophyll biosynthetic pathways were influenced after the treatments.

(PDF)

Click here for additional data file.

Quality check of sequencing reads.

(PDF)

Click here for additional data file.

Selected reference genes for bias structure analysis.

We analyzed 43 publicly-available microarray datasets including 891 difpan class="Chemical">ferent barley samples to find the most stable reference genes.

(PDF)

Click here for additional data file.

Differentially expressed events found in different comparisons.

(XLSX)

Click here for additional data file.

Genes showing alternative isoform switching among transcripts which code for an identical amino acid sequence.

(XLSX)

Click here for additional data file.

Accession numbers of the genes discussed in the main text.

(PDF)

Click here for additional data file.

Real-time PCR primers.

(PDF)

Click here for additional data file.

Read mapping statistics.

(PDF)

Click here for additional data file.

Differentially expressed genes involved in chromatin remodeling.

(PDF)

Click here for additional data file.

Differentially activated retroelements.

(XLSX)

Click here for additional data file.

Differentially expressed genes involved in DNA-damage response.

(PDF)

Click here for additional data file.

Differentially expressed genes involved in phosphatidylinositol-3-phosphate biosynthesis, recycling, and usage.

(XLSX)

Click here for additional data file.

Differentially expressed genes positively involved in biotic stress response.

(PDF)

Click here for additional data file.

Quality check of RNA samples.

(PDF)

Click here for additional data file.

Iron- and zinc-triggered stress responses.

It describes some of the stress-response genes found differentially expressed after the treatments.

(PDF)

Click here for additional data file.

Excess zinc and iron affect photosynthesis and respiration.

(PDF)

Click here for additional data file.