Comparative Transcriptional Profiling of Two Wheat Genotypes, with Contrasting Levels of Minerals in Grains, Shows Expression Differences during Grain Filling

Wheat is one of the most important cereal crops in the world. To identify the candidate genes for mineral accumulation, it is important to examine differential transcriptome between wheat genotypes, with contrasting levels of minerals in grains. A transcriptional comparison of developing grains was carried out between two wheat genotypes- Triticum aestivum Cv. WL711 (low grain mineral), and T. aestivum L. IITR26 (high grain mineral), using Affymetrix GeneChip Wheat Genome Array. The study identified a total of 580 probe sets as differentially expressed (with log2 fold change of ≥2 at p≤0.01) between the two genotypes, during grain filling. Transcripts with significant differences in induction or repression between the two genotypes included genes related to metal homeostasis, metal tolerance, lignin and flavonoid biosynthesis, amino acid and protein transport, vacuolar-sorting receptor, aquaporins, and stress responses. Meta-analysis revealed spatial and temporal signatures of a majority of the differentially regulated transcripts.


Introduction
Micronutrients play an important role in metabolism, such as in the production and functioning of enzymes, hormones and other substances. Therefore, an adequate intake of trace elements is necessary for proper growth and development [1]. Prevalence of mineral deficiency, especially of iron and zinc, has been estimated in more than two billion people worldwide (http://www.unicef. org/). In developing countries, every second pregnant woman and about 40% of preschool children are estimated to be affected by the mineral micronutrient deficiency (http://www.who.int/ nutrition/topics/ida/en/). The extent of reliance on staple crop based diet is one of the major reasons for mineral micronutrient malnutrition in human beings, mainly in resource-poor countries [2]. Development of mineral enriched cereal grains through breeding or biotechnological interventions is important for addressing micronutrient insufficiency [3].
Bread wheat (Triticum aestivum L.) is a globally important cereal crop, accounting for 20% of the world's daily food supply [4]. Most high yielding wheat cultivars have low content of grain mineral nutrients [5,6]. The primitive wheat genotypes constitute a pool of significant variability for grain micronutrients, which can be utilized for breeding and/or genetically engineering wheat cultivars [7]. Landraces, the primitive cultivated genotypes, often exhibit elevated grain mineral concentration and biotic and abiotic stress tolerance [8,9]. For example, a primitive cultivar, T. aestivum L. IITR26 exhibits higher grain mineral concentration than the post-green revolution modern cultivar, T. aestivum Cv. WL711 [10,11,12]. We have earlier reported a distinct level of mineral distribution in maternal and filial grain tissues of the two wheat genotypes, IITR26 and WL711 [10,11]. IITR26 is more efficient in accumulating micronutrients (e.g. Fe, Zn and Mn) in grain tissues than WL711. Furthermore, comparatively higher occurrence of nutritionally important minerals was observed in the endosperm of IITR26 grains, than WL711 [10,11]. Differential expression profiling of genotypes with contrasting characters is an efficient tool for understanding the molecular basis of such phenotypic differences [13]. Therefore, comparative transcriptional profiling in the two wheat genotypes, IITR26 and WL711, during grain development was done in this study.
The three major grain developmental phases in wheat are 6-10 days after anthesis (DAA), 12-21 DAA and 28-42 DAA [14] (http://www.wheatbp.net). During grain filling, mineral nutrients are transported by the maternal tissues to filial grain tissues via endosperm cavity and transfer cells [15,16]. Grain growth rate has been recorded as the highest during 7 th to 28 th DAA and lowest during 28 DAA till maturity [17]. The grain developmental stages at, 14 and 28 DAA are major transition phases in wheat grain development with distinctive pattern of transcript abundance [14,18]. After 28 DAA, the grain starts to desiccate and gradually attains physiological maturity at about 42 DAA [14].
Microarray is a useful technique for accurate and high throughput gene expression analysis, though it is relatively less sensitive in detecting rare transcripts [19,20]. Microarray technology has been extensively utilized to examine transcriptome related to grain development in wheat [14,21,22,23,24,25,26,27]. We used Affymetrix GeneChip Wheat Genome Arrays to examine the transcriptome changes in developing grains of two genotypes, IITR26 and WL711, at 14 and 28 DAA. The study identifies several transcripts showing significant differential expression during grain filling, which could be candidate genes for mineral accumulation in grains at a higher level. The data set can be further examined for understanding the regulation of mineral accumulation and sequence polymorphism in the candidate genes.

Plant material
The two wheat genotypes, a wheat cultivar (T. aestivum Cv.WL711) and a landrace (T. aestivum L. IITR26), with contrasting levels of minerals in grain tissues [11], were grown in the experimental field of National Agri-Food Biotechnology Institute (NABI), SAS Nagar (Mohali) Punjab, India (310 m above sea level; Latitude 30u 47' North; Longitude 76u 41' East). The ears of each genotype were tagged at anthesis, in three replicates. The spikes were harvested at 14 and 28 DAA, frozen in liquid nitrogen immediately after harvesting and stored at 280uC for RNA extraction.

Mineral concentration analysis
About 25 seeds of 14 DAA and 28 DAA were lyophilized (VirTis, sentry 2.0, USA) and crushed in a clean mortar -pestle. Approximately 100 mg powder was used for acid digestion in the Microwave Reaction System (Mars 6, CEM Corporation, USA). Mineral concentration was estimated in the digested samples using inductively coupled plasma mass spectrometry (ICP-MS; 77006AgilentTechnologies, Santa Clara, CA), following standard protocol.

RNA preparation
For RNA extraction, developing grains were extracted from the first and second floret of spikelet from the middle portion of the ear. Three grains were taken for RNA extraction in each sample. Total RNA was extracted using TRIzol (Invitrozen) and then purified with RNeasy Plant Mini Kit (Qiagen) following the manufacturer's protocol. On-column DNase (Qiagen) digestion was performed as instructed in the manual. RNA integrity was determined using the RNA 6000 Nano assay bioanalyzer (Agilent).

Microarray hybridization and data analysis
Double stranded cDNA synthesis, in vitro transcription to synthesize biotin labeled aRNA, purification and fragmentation of aRNA, and hybridization of arrays were performed according to Affymetrix technical manual. The Affymetrix GeneChip Wheat Genome Array contains 61,290 probe sets, representing about 25 K unigenes. The hybridized chips were washed, stained and scanned using the GeneChip scanner to generate the CEL files.
The CEL files were imported into GeneSpring GX v12 (Agilent Technologies). Signal intensities were recorded for all the 61,290 probe sets. The data has been deposited at NCBI (http://www. ncbi.nlm.nih.gov), with accession number GSE56032. The signal intensities were normalized by using Robust Multi-array Average (RMA) algorithm [28]. The Principal Component Analysis (PCA) in GeneSpring GX v12 established that the three biological replicates were located close to one another. The high correlation coefficient was observed among the three replicated samples, indicating less genetic background noise (Figure 1a). To correct the variability in the normalized expression values, the probe sets with coefficient of variation ,50% were retained, and the rest were discarded. One-way analysis of variance (ANOVA, p#0.01) was done with multiple testing correction algorithms using the Benjamini Hochberg approach with false discovery cut off, q # 0.01. Volcano plot with t-test unpaired computation was performed (p#0.01), and the transcripts showing statistically significant expression difference of two fold (log2) or above were selected for further analysis (Figure 1b). A log2 ratio = 1 indicates a change in the expression of two fold.

Functional annotation and gene ontology analysis
Functional annotation of the differentially expressed probe sets (with log2 fold change of $2 at 14 and/or 28 DAA) was done by employing HarvEST Assembly: XW Affymetrix Wheat 1 version 1.59 (www.harvest.ucr.edu). Blast2GO annotation of the differentially expressed transcripts was performed for each probe set query sequence to enrich the annotation and to generate the combined graph of biological process. Hierarchical clustered heat map was produced for the differentially expressed transcripts, using MEV (version 4.6.2) [29]. Enrichment of Gene Ontology (GO) terms in significantly differentially expressed genes ($2 log2 fold change) was evaluated using AgriGO analysis tool (http://bioinfo.cau.edu. cn/agriGO) with Fisher tests and Bonferroni multiple testing correction (P,0.05) [30].

In silico analysis of wheat transcripts
The differentially regulated transcripts were categorized by using MapMan software version 3.6.0RC1 (http://mapman. gabipd.org/web/guest/mapman). Spatial and temporal distribution pattern of the differential transcripts was examined at Genevestigator (https://www.genevestigator.com/gv/plant.jsp) using expression data of 1532 samples in T. aestivum microarray database of Affymetrix.

Quantitative real time PCR
Total RNA was extracted from the developing grains (14 DAA) of the two genotypes (IITR26 and WL711), in three biological replicates. First strand cDNA was synthesized by using Super-Script III First Strand Synthesis Kit (Life Technologies). Primers were designed using Primer Express software. Quantitative real time PCR was performed in three biological replicates using SYBR Green (Qiagen, USA) fluorescence dye and analyzed by 7500 Fast Real-Time PCR System (Applied Biosystems). The qRT-PCR was performed using gene specific primers (Table S1 in File S1), and analyzed as described previously [31]. The wheat 18S gene was used as the internal control to normalize the expression data.

Mineral concentration
The two genotypes, IITR26 and WL711, were analyzed for difference in mineral concentration in developing grains. The concentration of nutritionally important micronutrients, such as iron, zinc and manganese, was found significantly higher in the developing grains of IITR26, as compared to WL711 ( Figure 2). This is in agreement with our previous reports which show a

Comparative transcriptional profiling and data analysis
Transcriptome analysis of developing grain was compared between two wheat genotypes (IITR26 and WL711). A total of 433 and 315 probe sets (with log2 fold change of $2 at p#0.01) were identified as differentially expressed between IITR26 vs WL711 on 14 and 28 DAA, respectively ( Figure 3a). Putative gene function was assigned to each probe set (Table S2 in File S1) using public databases as mentioned in materials and methods. A total of 149 differentially regulated probe sets could not be assigned to any function, out of them 56 were common between 14 and 28 DAA. Among these, the GO enrichment patterns show disproportionately represented transcripts of genes involved in the biological process of response to stress, polysaccharide metabolism, multiorganism processes, and the biological function of nutrient reservoir and hydrolase activities in IITR26 developing grains (Table 1). In the developing grains of WL711, the GO biological processes-lipid transport, DNA replication, DNA metabolic process, starch metabolic process, and the functions-lipid binding, enzyme regulation, hydrolase and peptidase inhibitor activities were over-represented ( Table 2). The functionally annotated transcripts were categorized into 15 bins by using MapMan software ( Figure 3b) to understand global changes in developing grains of the two contrasting genotypes. Hierarchically clustered heat maps were produced for the 580 differentially expressed transcripts ( Figure S1 in File S1). Meta-analysis indicated spatial and temporal distribution of the differentially regulated probe sets, assigned to putative locations of 22 anatomical plant parts based

Candidate genes involved in the transport and accumulation of minerals and other solutes
Transcripts related to metal homeostasis showed differential expression levels in the developing grains of the two genotypes, IITR26 and WL711 (Table 3). Metallothionein (MT) was upregulated in IITR26 developing grains as compared to WL711. MT is a cystein rich protein having the capacity to bind with metals through its thiol group. It plays an essential role in metal ion homeostasis, transport and storage in the cell [32,33,34]. In the seed, its abundance has been reported in aleurone and embryo [35], whereas endosperm shows comparatively less expression. This coincides with aleurone and embryo being hot spots for mineral localization [10]. Elevated transcription of Manganese ion and Transition metal ion binding proteins was observed in IITR26, mainly during 14 DAA (Table 3). Transition metal binding protein, an orthologue of Farnesylated protein 3, is involved in executing metal ion homeostasis and, biotic and abiotic stress tolerance [34]. Its expression has been recorded mostly in seed tissues ( Figure 4).
A transcript of the NAM-1 gene was found up-regulated in developing grains of IITR26 during grain filling. RNAi mediated suppression of NAM in wheat resulted in lower grain mineral and protein concentrations [36,37]. Further, expression of the probe set was found limited to endosperm in meta-analysis ( Figure 4). Seed specific NAM-1 could be one of the candidate genes which play role in mineral accumulation in wheat, during grain filling.
Late embryogenesis abundant (LEA) proteins play crucial role in conferring tolerance to metal stress [38,39]. LEA proteins are suggested to be involved also in sequestering divalent metals, and micronutrient trafficking in phloem tissues [40,41]. Enhanced expression of LEA-12 could facilitate higher level of mineral accumulation in the grains of IITR26.
Elevated expression of lignin biosynthesis related genes have been observed, mainly during metal stress, in metal hyperaccumulator plants [42]. The up-regulated Cinnamoyl CoA reductase gene, an enzyme of the lignin branch biosynthetic pathway, presumably provides metal tolerance in developing grains of IITR26. Up-regulation of Flavonoid biosynthetic pathways genes has been notified under metal and biotic stresses [43]. Flavonoids, such as quercetin, have capacity of chelation of transitional metals [44]. Flavanone 3-hydroxylase has sites for putative iron binding [45]. Genevestigator analysis revealed the expression of Flavanone 3-hydroxylase and Dihydroflavonol-4reductase genes largely in seed tissue ( Figure 4). Their up-  regulation, during early grain filling stage in IITR26, could be helpful in favorably accumulating more transition metals [46]. In WL711 developing grains, transcripts for Protein phosphatase 2c-PP2c were significantly up-regulated (Table 3). PP2c is involved in the regulation of cation transporters and membrane polarization, influencing ion homeostasis in plant cells during stress and during phloem transport [47,48]. The Calreticulin-1 gene was up-regulated in WL711 developing grains, and its presence has been found mostly in seed tissues ( Figure 4). It codes for multifunctional protein that binds Ca ions (as a second messenger) with low affinity, but high capacity, and its role has been demonstrated in drought and metal stress tolerance [49]. A probe set representing Glycolate oxidase gene, an important player in oxalate biosynthesis, was significantly up-regulated in the developing seeds of WL711 (Table 3). Oxylate has been notified to play important role in calcium regulation and metal detoxification, and its access reduces calcium bioavailability [50]. Nicotianamine (NA) is a metal chelator. NA is synthesized by trimerization of three molecules of S-adenosyl-Lmethionine by nicotianamine synthase (NAS). The role of NAS-3 has been reported in mineral homeostasis in developing organs in maize [51]. Induced expression of NAS3 has been recorded during excessive mineral, while expression is down regulated under mineral deficient condition [52]. The up-regulation of NAS-3, at later grain filling stage (Table 3), could be indicative of mineral sufficiency within the grain of WL711, as compared to IITR26.
In wheat, aleurone is a single cell layer storage tissue which accumulates minerals, phytic acid and proteins in globoids or protein storage vacuoles (PSV), within the cell [53]. High resolution X-ray micrograph of aleurone cell revealed the well defined edges of globoids, surrounded by a membrane, tonoplast [54]. Aquaporins, localized in the tonoplast of PSVs [55,56],      facilitate transmembrane transport of water and solutes, and provide cellular adaptation to metal stress by water homeostasis [57,58]. About three fold increased expression of aquaporin-Tonoplast intrinsic protein (TIP-3) was recorded in IITR26 developing grains (Table 3). Meta-analysis showed its abundance in seed tissues as compared to other organs ( Figure 4). Some other aquaporins, Plasma membrane intrinsic protein (PIP-1) and Nodulin 26-like intrinsic protein (NIP-1), play an important role in osmo-adaptation and homeostosis of cells to facilitate transport of water and phloem derived nutrients to grain tissues [59,60]. However, NIP-1 has been notified as less efficient in transporting solute and ions [61]. Enhanced expression of TIP-3 and PIP-1 could promote higher level of mineral sequestration in IITR26 grains ( Table 3). The aleurone globoids or PSVs store large number of organelle-specific proteins, synthesized by endosplasmic reticulum associated ribosomes and transported through vacuolar sorting [62]. Vacuolar sorting receptor-1 (VSR-1) is involved in sorting of proteins in PSV [63]. In meta-analysis, the predominance of VSR-1 was observed in endosperm ( Figure 4). The enhanced expression of VSR-1 in IITR26 (Table 3) could be suggestive of higher storage activity [64] in IITR26 grains. Several probe sets related to lipid transport and binding (LTP family protein) were significantly up-regulated in the developing grains of WL711 (Table 3). The expression of many of them has been recorded predominant in grain, as compared to other anatomical part (Figure 4). The enhanced expression of GDSL esterases/lipases, that catalyze hydrolysis of ester bonds-primarily in triglycerides, indicates relatively more active lipid metabolism [65] in WL711 developing grains. LTPs regulate lipid transport, homeostasis and lipid-mediated cellular processes [66]. Within grain, several LTPs have been observed up-regulated in aleurone cells around the embryo, scutellum and vascular bundle of the embryo, resulting maximum accumulation of lipids in embryo in grains [67,68]. Their role in providing protection to embryo from pathogen attack has been speculated [65,67,68]. Unlike WL711, transporters for amino acid and protein (SecE) [69] were expressed at higher level in IITR26 during grain filling (Table 3). This is in agreement with previous reports showing significant correlation between grain protein and micronutrient content [15,70,71]. The results also coincide with our observation of proportional association between elemental images of sulphur (possibly sulphur-rich peptides) and iron in cross-sections of the two contrasting genotypes [10].

Candidate genes representing stress responses and metal tolerance
Blast2GO results of combined graphs distributed the differentially expressed transcripts in different levels of biological processes, mapping the transcripts related with metal/cation transport, localization and responses etc. The results revealed transcripts representing stress response to metal ion activity in developing grains of IITR26, such as peroxins (Ta.9561.3.-S1_a_at), VSR-1 (Ta.9898.1.A1_at), PIP-1 (Ta.2895.1.S1_x_at), Ankyrin repeat protein (Ta.23465.1.S1_at) and Flavonoid (Ta.8881.1.S1_at) ( Figure 5). Peroxisome biogenesis proteins (peroxins) are involved in docking cargo-receptor complexes at the peroxisomal membrane [72]. Several studies revealed essential role of peroxins in defense mechanisms conferring resistance against pathogen attack [73,74]. VSR-1 is a candidate gene controlling the defense network, and its up-regulation has been reported in metal tolerant and disease resistant plants [75,76]. Involvement of plasma membrane aquaporins (PIP-1) has been discussed in metal homeostasis and defense response in plants [77,78,79]. Flavonoids enhance protective abilities of the plant against biotic stresses [80,81]. Ankyrin repeat proteins have been reported to mediate metal tolerance and enhanced disease resistance in plants [82,83].
Several genes associated with abiotic and biotic stress were found differentially regulated in developing grains of both the genotypes (Table 4), presumably having coordination among the growth, maintenance and stress during seed development [84]. For example, the transcripts for b-purothionin, cysteine-rich proteins, and Defensin, up-regulated in WL711, are known to play significant role in protecting wheat grains from microbial pathogens [68,85]. In IITR26 developing grains, transcripts of Puroindolin-A (PINa), a cystine-and tryptophan-rich protein, was significantly predominant. However, the probe set for PINb did not pass the quality criteria. A relationship has been incurred between grain hardness and puroindoline content [86]. However, both the genotypes, IITR26 and WL711, exhibit similar grain hardness indices (,90) [12,87,88]. This indicates contribution of PINa in grain related traits, other than hardness. Role of PINa  protein has been established in plant defense as antimicrobial molecule [85]. Transcripts of many other antimicrobial molecules, Chitinase, Glycin rich proteins, Wheatwin, VAMP, Bowman-Birk and Kunitz-type protease inhibitor, were significantly up-regulated in the developing grains of IITR26 (Table 4). A positive correlation has been inferred between degree of biotic stress resistance and glucan synthase activity, on the contrary glucanase activity is negatively correlated [89]. The enhanced glucan synthase activity and reduced glucanase activity in IITR26 developing grains further corroborates enhanced resistance to pathogen attack in IITR26 (Table 4). Several fold reduction in superoxide dismutase activity in the developing grain of IITR26 (Table 4), could be indicative of less reactive oxygen species (ROS) formation [90]. Mampman analysis of cellular response overview reflected up-regulation of redox related transcripts in WL711 during grain filling ( Figure 6). This is in agreement with the previous report of induced redox state in susceptible wheat as compared to the biotic stress tolerant line [91]. Meta-analysis determined dominant expression potential of several stress related transcripts in grain tissues (Figure 4), as compared to other plant organs, suggesting their role in seed development. The similarity search in Genevestigator, using the differentially expressed transcripts, revealed perturbations in which spikelet samples of wheat genotypes, with stress resistant (CS-7EL) and susceptible (CS) backgrounds (GEO accession GSE21386), have been compared ( Figure S3 in File S1). The results coincide with the previous reports examining interactions between metal ion excess and biotic stress factors [92,93], and the responses to excess metal ions resemble to that of biotic and abiotic stresses [94,95]. Our results support the elemental defense hypothesis which suggests a crosstalk between metal accumulation and biotic and abiotic stress responses [96,97,98,99]. However, it needs to be investigated and validated in other wheat genotypes, with contrasting levels of minerals in grains, if the crosstalk substantiates the enhanced level of mineral sequestration in grains.
Pectin methylesterase, which catalyses the de-esterification of pectin, is regulated by pectin methylesterase inhibitor. Its role has been discussed in wounding, osmotic stress, senescence and seed development [100]. Significant up-regulation of several pectin methylesterase inhibitor related probe sets has been observed in WL711 (Table 4). Down regulation of pectin methylesterase inhibitor activity might inhibit de-esterification of pectin, which could result into reduced metal tolerance in the tissue [101].
The developing grains of IITR26 exhibited up-regulation of transcripts for F-box and NAC genes (Table 4). F-box proteins, apart from its role in seed development [102], provide tolerance during metal stress in plant tissues [103]. NAC proteins play a key role during stress signaling pathways and provide tolerance against various abiotic and biotic stresses [104]. Furthermore, the abscisic acid (ABA)/stress-induced proteins [105], such as HVA22 were expressed at higher level in IITR26 developing grains ( Table 4).
The results are suggestive of IITR26 being comparatively more tolerant to stress during grain filling, which eventually could be a positive factor for metal accumulation.

Quantitative RT-PCR and similarity search meta-analysis
The relative fold change in expression of a few genes related to the transport and accumulation of metal was examined by real time PCR analysis (Figure 7). The expression of Metallothionein, NAM-1 and LEA-12, and protein, SecE, were found high in the developing grains of IITR26. The expression profile of these genes was in agreement with the microarray results. Hence, q-RT PCR results validated the microarray data, as reported in several published reports.
Similarity search meta-analysis was performed against 1532 experiments listed in Genevestigator, by using the metal related transcripts differentially expressed between IITR26 and WL711 ( Table 3). The perturbation showing maximum similarity with our data comprised the transcriptional comparison between the developing grains of two wheat genotypes: LOK-1 and WH291 (Figure 8a). The mineral concentration analysis revealed relatively higher concentrations of micronutrients in mature grains of LOK-1 than WH291 (Figure 8b), which coincides with up-regulation of the metal related transcripts in early-stage developing grains (Figure 8a), supporting our data.

Conclusion
The microarray based comparative transcriptome profiling of developing grains of two genotypes contrasting for grain mineral concentration identified several differentially expressed genes. Prevalence of transcripts related to metal transport and accumulation, metal tolerance, vacuolar sorting receptor-1, aquaporins (TIP-3 and PIP-1), amino acid and protein transporters, lignin and flavonoid biosynthesis, and stress responses, was noticed in IITR26 during grain filling. For some of the differentially expressed probe sets, gene functions are yet to be determined (Table S2 in File S1). Thus, the differential expression analysis of microarray data revealed several candidate genes which may facilitate the elevated levels of minerals in the grains. However, the possibility of involvement of the differential transcripts in establishing other traits, for which the two genotypes differ, may not be denied. The differential expression data provided in this study would be useful for designing programs for functional genomics, exploration of sequence polymorphism in the candidate genes, and understanding other metabolic pathways common with mineral accumulation in wheat grains.