Transcriptomic Analysis of Soil Grown T. aestivum cv. Root to Reveal the Changes in Expression of Genes in Response to Multiple Nutrients Deficiency

Deficiency of necessary macronutrients, i.e., Potassium (K), Magnesium (Mg), Nitrogen (N), Phosphorus (P), and Sulfate (S) in the soil leads to a reduction in plant growth and yield, which is a result of changes in expression level of various genes. This study was performed to identify the differentially expressed genes and its associated metabolic pathways occurred in soil grown wheat root samples excavated from the control and treated fields. To identify the difference in gene expression levels due to deficiency of the said nutrients, a transcriptomic, meta-analysis was performed on array expression profile data. A set of 435 statistically significant probes encoding 398 Nutrient Deficiency Response Genes (NRGs) responding at-least one nutrients deficiency (ND) were identified. Out of them 55 NRGs were found to response to minimum two ND. Singular Enrichment Analysis (SEA) predicts ontological based classifications and functional analysis of NRGs in different cellular/molecular pathways involved in root development and growth. Functional annotation and reaction mechanism of differentially expressed genes, proteins/enzymes in the different metabolic pathway through MapMan analysis were explored. Further the meta-analysis was performed to revels the active involvement each NRGs in distinct tissues and their comparative potential expression analysis in different stress conditions. The study results in exploring the role of major acting candidate genes such as Non-specific serine/threonine protein kinase, Xyloglucan endotransglucosylase/hydrolase, Peroxides, Glycerophosphoryl diester phosphodiesterase, S-adenosylmethionine decarboxylase proenzyme, Dehydrin family proteins, Transcription factors, Membrane Proteins, Metal binding proteins, Photosystem proteins, Transporter and Transferase associated in different metabolic pathways. Finally, the differences of transcriptional responses in the soil-grown root of T. aestivum cv. and in-vitro grown model plants under nutrients deficiency were summarized.


INTRODUCTION
Potassium (K), Magnesium (Mg), Nitrogen (N), Phosphorus (P), and Sulfate (S) are necessary macronutrients present in the soil for plant growth and development. Each nutrient acts in a specific manner in a variety of physiological and morphological aspects of plant growth. The plant growth and development is a result of different cellular and metabolic processes (Marschner and Rimmington, 1988;Maathuis, 2009). The K is an essential micronutrient which takes up about 2-10% of total dry weight in plants (Leigh and Wyn Jones, 1984). Earlier study reported that the K plays an important role in enzyme activation, protein synthesis, photosynthesis, turgor pressure, osmoregulation, and electrical neutralization (Römheld and Kirkby, 2010). Hence, the deficiency of the K brings many negative impacts on plant growth due to inappropriate osmotic pressure, nutrition imbalance and changes in photosynthesis, protein synthesis and other molecular effects. The deficiency of the Mg leads to improper functioning of chlorophyll, required for capturing the sun's energy during photosynthesis process (Cakmak and Yazici, 2010). The cellular molecules such as adenosine triphosphate (ATP), amino acids, nucleic acids, chlorophyll, phospholipids and other plant hormones are made up of N, P, and S (Takehisa et al., 2013). Among these nutrients, the larger quantity of the N is required for the production of said bio-molecules therefore it leaches faster from the soil. It is essential in the development of chlorophyll and plays an important role in the stimulus of photosynthesis. The P is crucial macronutrient at the early stage of plant growth and development as it is required for maintaining the various processes within the plant. Similarly, the S also plays an important role in plant growth because of its involvement in the production of various metabolites such as amino acids, proteins, sulpho-lipids and production of chlorophyll (Sorin et al., 2015).
Wheat (T. aestivum) is a cereal grain crop cultivated worldwide and became second largest cereal crop in terms of yield after maize in 2016 (Crosson and Frederick, 2016). Recently, yield losses were observed in this crop due to different abiotic and biotic stress. Under abiotic stress, the ND affect majorly in the quality and quantity of the crop. Identification and understanding of the cellular mechanism of signaling, sensing and adaptation in response to ND is crucial to enhance productivity and quality of the wheat (Gill et al., 2016;Gupta et al., 2016bGupta et al., , 2017. Considering the complexity of above said mechanisms and processes, it is crucial to understand the functional elements and molecular constituents involved during these NDs at transcriptome level. The transcriptomic analysis provides genome-wide measurement of messenger RNA (mRNA) expression levels based on microarray technology.
Microarray technology is rigorously used to uncover the regulatory, signaling, sensing and adaptation mechanism of transcripts in response to nutrients deficiency (mainly N, P, and K) in model plants (Wasaki et al., 2003(Wasaki et al., , 2006Wu et al., 2003;Armengaud et al., 2004;Misson et al., 2005;Lian et al., 2006;Li et al., 2010a,b;Krapp et al., 2011;Cai et al., 2012;Ma et al., 2012;Park et al., 2012). These in-vitro grown customized experiments elucidate that several transcriptional responses occur when plants were subjected to various nutrients deficiency. Gruen et al. (2015) performed gene expression profiles of the soil-grown wheat root to analyze the effect of combined nutrient deficiencies (https:// www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE61679). They have generated transcriptomic data of wheat root samples grown in limited nutrients field known as soil grown, which is different from in-vitro grown samples. The major difference in the said soil grown experiment with other in-vitro grown samples are (i) difference in types of soil, (ii) availability of inadequate nutrients in soils, and (iii) variation in environmental conditions. Hence the soil grown experiment conducted by Gruen et al. is found to be more close to the actual cultivation field conditions. This inspires us to carry a differential expression study on said soil grown data and to compare and verify its results with other invitro grown samples. In the present study, nutrient deficiency responsive genes (NRGs) in response to distinct NDs conditions were explored. Along with GO term, MapMan based analysis of each NRGs helps in the identification of function and their role in metabolic pathways. The meta-analysis was performed to uncover the NRGs involvement in distinct tissues of wheat and expression potential in different ND and stress conditions using 109 wheat root specific samples (Hughey and Butte, 2015). Used wheat root specific samples were retrieved from NCBI-GEO and EBI-ArrayExpress databases. This meta-analysis gives insights into the features of NRGs, which play crucial role in biological systems of wheat plant under abiotic stress. Finally, this study fills the gaps of earlier transcriptomic studies performed on model plants and soil-grown wheat, with an aim to identify the genes responsible in nutrient uptake mechanism.

Microarray Data and Experiment
Microarray data, (Accession no. GSE61679) was downloaded in order to identify the transcriptional responses in the root of wheat, grown in soil with limited nutrients availability. To generate transcriptomic data, the root material was excavated in triplicate at booting stage (At this stage, development of head within the sheath of the flag leaf becomes visible) from the control fields with all nutrients supply known as control sample and treated samples with limited nutrient supply in the fields. These results of 6 treatments (one control and five treated samples) for each replicate, altogether 18 samples of root material were collected for all the 3 replicates. The root growing environmental conditions fluctuates as per day/night weather and the average temperature reported was 18 • C (day)/10 • C (night). Further, the total RNA was extracted using modified protocols of Verwoerd et al. (1989) with an additional phenol-chloroform-isoamylalcohol extraction of the aqueous phase after the first centrifugation and subsequent treatment with RNAse-free DNAse [Promega, Madison, USA]. In addition, silica-membrane purification kit was used to purify the total extracted RNA [RNeasy Plant mini kit, Qiagen, Maryland, USA]. Further, the biotin labeling of cDNAs was performed and the biotin-modified cRNAs were obtained from 500 ng total RNA using the Affymetrix 3 ′ IVT Express kit [Affymetrix User Manual P/N 702646 rev.7]. The hybridization performed in GeneChip R Hybridization, Wash, and Stain Kit was used for cRNA hybridization according to the user manual P/N 702731 of GeneChip R Expression wash. After hybridization of the gene chips were washed in the Affymetrix GeneChip fluidics station 450 [Fluidics script EukGE-WS2v5]. Gene chip scanning was performed using the Affymetrix R GeneChip R Scanner 3000, with default settings of Affymetrix R GeneChip R Command Console R Software(AGCC). The whole experiment was executed in Affymetrix R GeneChip R Wheat Genome array platform (Accession no. GPL3802). This array contains 61,127 probe sets representing 55,052 transcripts for all 42 chromosomes of the wheat genome (http://www.affymetrix.com).

Microarray Data Analysis
The images of probe set signal values were examined through Affymetrix R GCOS 1.0 (MAS 5.0) software to generate summarized raw signal values in .CEL files. The raw data contains three biological replicates for each nutrient deficiency, including a control sample. A candidate gene was considered as expressed, when it appears in all these three replicates after RMA normalization, which is an inbuilt algorithm in GeneSpring GX 13.0 [Agilent Technologies, Inc]. For Affymetrix gene expression data, GeneSpring tool has two basic workflows i.e., (i) Analysis: Biological significance and (ii) Data import wizard. In our analysis, we have used 6 basic steps i.e., (1) Summary Report, (2) Experiment Grouping, (3) Quality check (QC) check on samples, (4) Filter Probe Sets, (5) Statistical Analysis, and (6) Fold Change of Analysis: Biological Significance workflow. Initially the summary report and grouping (Control vs. Nutrients deficient) for all sample replicates were generated. The QC of data was carried out using the Boxplot, Hierarchical clustering and Principal Component Analysis (PCA) (Raychaudhuri et al., 2000). The PCA calculates the principal components and contribution of an observation to component which is also known as contribution ratio, can be obtained by the ratio of the squared factor score of the observation by the eigenvalue associated with that component. Formally, the contribution of observation i to component j is, known as Ctr i, j and obtained by: where δ j eigenvalue of jth component and f 2 i, j is factor score for the principal components (Abdi and Williams, 2010). The boxplot shows the shape of the data distribution along with its central value (median) and variability (inter-quartile range). The hierarchical clustering is used to look the patterns of gene expression values in array data and PCA computes all possible principal components and shows it in a 3D scatter plot. The quantitative analysis (QA) was carried out through an initial filtering using percentile shift method with default parameter values. Thereafter in subsequent steps, the statistical and fold change analysis was carried out on filtered probe sets. A probe/gene was considered as differentially expressed OR denoted as NRGs, when it satisfied (i) Fold change (FC) value FC value = +2 for up-regulated or FC value = −2 for down-regulated genes (calculated for each gene of ND samples based on control sample) and (ii) corrected P-value = 0.05 (P-value computation was carried out through t-test utilizing Benjamini-Hochberg FDR as a multiple testing correction methods). Identification significant probe sets encoded genes IDs, name and their description were manually extracted from Genome array file (Accession no: GPL3802) and also revalidated using PLEXdb (http://www.plexdb.org/) and ARTRAdb (http:// artra.kazusa.or.jp/artra/ARI3_101/index.html) databases.

Evaluation of Gene Ontological Relation of NRGs
To understand various biological processes, molecular mechanism resulted due to the different ND in the wheat root, SEA of NRGs were performed through AgriGO analysis toolkit (http://bioinfo.cau.edu.cn/agriGO/). SEA analysis was performed using Fisher Test with 0.05 level of significance and Benjamini-Hochberg-Yekutieli Multi-test adjustment method using T. aestivum genome as a reference genome (Du et al., 2010;Gupta et al., 2016a).

MapMan Based Analysis
MapMan analysis used to find the role of specific gene in a metabolic pathway using prior MapMan annotations scheme (http://mapman.gabipd.org/web/guest/mapmanstore) and inhouse PERL programs as in Yadav et al. (2016) and Thimm et al. (2004). The macronutrients controlling genes were classified into a functional category within the hierarchy of MapMan pathway scheme (also known as BINCODE and file was downloaded from MapMan recourse). Further, the pathway enrichment analysis of functional category was performed by calculating the cumulative hypergeometric P-value. It is a probability of a gene group overrepresented within a functional BINCODE at a rate higher than chance expectation. The calculation of the hypergeometric Pvalue for 398 differentially expressed genes over 55052 sample genes in 9 different situations (i.e., Up-regulated genes under K, N, P, Mg, and S deficiency and Down-regulated genes under K, N, P, and Mg deficiency) was carried out in which 4182 genes are found to be associated with different metabolic pathways. Several tests were executed for all sets of functional gene categories at different hierarchical levels. The significant hypergeometric Pvalue of a gene/a genes group was used for generation of heatmap in the R.

Investigation and Meta-Analysis of NRGs
Gene investigation for all NRGs through Java based tool Genevestigator was carried out. This tool contains 2101 samples of different array experiments performed for wild type and non-wild type wheat in response to different stress conditions in its database [Genevestigator, Inc]. These wheat specific microarray experiments samples were fetched from NCBI-GEO and EBI-Array express. Since our analysis belongs to wheat root, exposed under different nutrients deficiency and generated through Affymetrix R GeneChip R wheat genome Array platform. Therefore, only the root specific 109 samples were found to be selected for this analysis. Genevestigator tool platform provides two types of analysis (i.e., Single Experiment Analysis and Compendium Wide Analysis). The Single Experiment Analysis is used to visualize the expression of NRGs across individual experiments and Compendium Wide Analysis is used to easily visualize the expression across thousands of experiments in a single analysis . In both analysis, NRGs were used as input genes to perform meta-analysis while comparing with selected samples. Functional annotation of characterized/uncharacterized proteins was performed using standard protocol used by Gupta et al. (2016b).

Identification of Nutrient Deficient Response Genes
The triplicate samples of roots of wheat were excavated at booting stage from the control soil as well as K, Mg, N, P, and S nutrients limited plots and total RNA was extracted by Guren et al. (2015). Further, Affymetrix R GeneChip R wheat genome Array platform generates total 18 raw signal values file (.CEL files) and used in GeneSpring tool for biological significance analysis (Bolstad et al., 2003). Normalization and other basic steps of biological significance analysis as discussed in methodology section were executed. Boxplot indicates shape of data distributions across all the samples after normalization ( Figure S1). The dendrogram resulted from hierarchical clustering indicates the similar expression patterns of microarray data ( Figure S2). The PCA was performed with an aim to simplify the description and analyze the structure of the samples of the data sets. In this analysis, the contribution ratio was also calculated for all treated samples with respect to control sample component. The projection of each sample on the principal component is represented through a 3D scatter plot ( Figure 1A). It indicates one point for each array and color them accordingly by grouping as per experiment factors through which one can visualize the separations between groups of replicates. In this analysis, each sample replicate falls within a group or clustered together, while getting separated from other sample groups (Raychaudhuri et al., 2000). The Mg deficiency group replicates cluster is found to have similar pattern as in the control replicates cluster ( Figure 1A). The QA of probes using percentile shift method resulted in 52,083 probes qualifying 0.75 quintile threshold value. Normalized probe sets were statistically tested using t-test (FDR = 0.05) and fold change analysis (FC ≥ +2.0 or FC ≤ −2.0) to check whether the candidate is differentially expressed in at least 1 out of 5 ND conditions against the control sample. This analysis identified NRGs, which qualify the defined screening threshold ( Figure 1B) and their corresponding FC values are listed in Table S1.  (Figure 2A) in root were expressed while comparing with other ND. Higher number gene expression lead to important changes in metabolic process such as alteration in lipid metabolism, compound rearrangement in cell wall, carbon supply acceleration for organic acid synthesis as earlier reported in rice root (Wasaki et al., 2003;Li et al., 2010b). The second highest numbers of genes were expressed under ND_N ( Figure 2B). These genes are associated with diverse cellular function i.e., organ development, protein biosynthesis etc. (Gruber et al., 2013). Venn diagram illustrating the number of common genes and unique genes expressed among N P K Mg, and S deficiency and their sub sets are shown in Figure 2C. Out of 435, total 58 probes representing 55 transcripts were found to be responsible for at least two NDs, which indicates their importance in comparison to others genes ( Table S2). The expression patterns of 55 transcripts across all NDs depicted in terms of heatmap (Figure 3). The heatmap diagram indicates that 7 genes (i.e., BQ170823, CA746776, CA745229, AJ698954, CA745977, CA747262, and CA744247) are up-regulated in four NDs, which infers that these genes are crucial for different biological processes and metabolic pathways associated with abiotic stresses.

Singular Enrichment Analysis of NRGs
The SEA of 435 probes was executed to annotate their Gene Ontology (GO) terms. Only 291 probes were satisfying the significance criteria (given in material and methods section) and GO term were assigned for each probe (Table S1), while for remaining 144 probes GO terms were not identified (Du et al., 2010). Each probe name was retrieved from Affymetrix R GeneChip R wheat genome Array platform file and revalidated through comparative analysis with Rice and Arabidopsis ontology using PLEXdb and ARTRAdb databases. Analysis illustrates that number of NRGs participation in biological processes, molecular function and cellular components decrease according to their sequence. In the biological process, most of the NRGs were classified into cellular process and metabolic process Viz. cellular biogenic amine biosynthesis, response to stress, response to stimulus, inorganic anion transport and cellular glucan metabolism ( Figure 4A). Within the molecular function category, heme binding, iron ion binding, peroxidase activity, nicotianamine synthase activity, antioxidant activity, oxidoreductase activity, adenosylmethionine decarboxylase activity are the subcategories of GO terms in which genes are playing role ( Figure 4B). In cellular component category, a significant number NRGs were classified into vesicle, cytoplasmic membrane-bounded vesicle, membrane and integral to membrane.

MapMan Based Pathway Analysis
To investigate the effect and changes occurs in wheat root due to ND through different metabolic pathways, a MapMan based pathway analysis was performed (Thimm et al., 2004). In this analysis, MapMan wheat genome array file containing BINCODE, name, and description of pathway (http://mapman. gabipd.org/web/guest/mapmanstore) was used to annotate the   pathways influenced during different ND conditions using inhouse PERL scripts. The mappings of transcriptome data define the functional categories and identify the significant role of different pathways. The plant root expansion is directly related to growth of shoots and formation of booting stages. During these stages, essential macronutrients such as K, P, S, Mg, N etc. are gradually accumulated through root. Table 1 shows the name and description of genes participating in different pathways/process with their corresponding hypergeometric Pvalue in Tables S3, S4. Total 110 BICODEs were associated with 398 NRGs and falls within superset of 17 BINCODE and 19 BINCODE or pathway (highly enriched) for up-regulated and down-regulated genes respectively, illustrated in Figure 5. This figure indicates most of the down-regulated genes under the P deficiency, participates in metabolic pathways such as carbohydrate (CHO), nitrogen, polyamine, co-factor, vitamin metabolism and synthesis pathways of cell wall, cell wall cellulose, DNA/chromatin, RNA, glycerol, lipid, peroxides, fatty acid and their elongation, although a few of these show active involvement in abiotic stress (i.e., cold, drought and salt), cytochrome P450 pathway and major intrinsic protein transport. Further, the pathway analysis of up-regulated genes under the P deficiency inferred their involvement in the transport of phosphate, synthesis of ribosomal protein, polyamine and lipid metabolism, RNA-mediated regulation of transcription and MYB-related transcription factor activity ( Figure 5A). The pathways analysis of down-regulated genes under the N deficiency indicates their role in Photosystem-I, II, and III light reactions, lipid degradation, glycolipid synthesis and transport of potassium ( Figure 5B). Up-regulated genes (CA683846, CK163751) under the N deficiency are involved in the reaction mechanism of amino acid metabolism and its degradation, while down-regulated genes (BJ219294, BJ247599) participated in the activation/deactivations of CHO metabolism. The pathways in which up-regulated differentially expressed genes play an important role under the K deficiency is CHO metabolism, hormone metabolism, signaling G-proteins, lipid metabolism/degradation and metal handling process ( Figure 4A). Gene CA650490 is involved in the signaling pathways of plant defense hormone Jasmonate acid biosynthesis and its degradation (Turner et al., 2002). The ethylene synthesis is involved in the inhibition of cell division, DNA synthesis, growth of the meristems and axillary buds (Burg, 1973). Up-regulated genes (BT00925, BE446498, and CA650490) under Mg deficiency participated in amino acid and hormone metabolism while only 3 down-regulated genes (BG908438, BJ307743, and CA725003) are associated with lipid metabolism. Only 6 up-regulated genes under the S deficiency are found to be associated with the metal binding process, polyamine metabolism and cellular developments (Table S3). Overall MapMan analysis predicts the metabolic pathways and associated genes plays

Meta-Analysis of NRGs
Meta-analysis provides a significant analysis of the genes of interest based on a large-scale systematic combination of normalized and quality-controlled expression data with experimental context variables using different ontologies i.e., anatomy development, perturbation, or genetic background (Hughey and Butte, 2015). Meta-analysis of NRGs was performed using the Genenvestigator tool, which has wild type and non-wild type 2101 samples of Affymetrix R GeneChip R Wheat Genome Array in its knowledgebase. These samples were generated under 164 different conditions in seedling, inflorescence, shoot and rhizome across the world. The root specific 109 samples for single experiment analysis and compendium wide analysis was selected (Table S5) for the said meta-analysis. Single experiment analysis annotates the percentage expression of each probe across the root specific wheat samples. The percentage of gene expression potential of each NRGs across the selected 109 root specific samples were generated through Single experiment analysis ( Figure S3). Figure S3 indicates that large numbers of NRGs were differentially expressed while wheat root expose to NDs, biotic and abiotic stress and dehydration. Further, the compendium analysis facilitates gene expression intensity on a particular tissue with their favorable conditions of gene regulation. A hierarchical clustering method using Euclidean distance function was used to annotate the tissue-specific expression with respect to all samples . Clustering of all probes using Genevestigator reveals the level of expression patterns of NRGs in 22 different tissues of the wheat (Figure S4), which is very crucial to identify the function of particular genes. The higher level of expression of a candidate gene in particular tissue will be helpful in inferring its role. Similarly, tissue based expression potential of 58 probes playing a crucial role in at least two ND conditions were also annotated and depicted in Figure 6. Moreover, in compendium analysis, the signature tool was used to perform comparative analysis in order to identify other stress conditions, which causes the similar gene expression signature using selected samples. The signature tool uses selected standard datasets, to find out whether the NRGs expression trend in selected samples (at a given condition) is similar or contradictory to the expression trend of standard datasets ( Figure S5). Similar analysis was also performed for 58 NRGs and same condition is observed as in figure Figure S5. The comparative expression potential of 55 genes in stress responsive conditions and is shown in Figure 7. Further, it has been observed out of 398 NGRs, only 137 genes are found to encode the protein sequences, which were annotated using standard sequence analysis protocols. Out of 137 proteins, 42 characterized and 95 uncharacterized proteins were identified ( Table 2). In addition, the expression signature of these proteins in specific macronutrient deficiency conditions were termed as "Low, " "High, " and "Medium" ( Table 2). These terms were calculated according to respective color values obtained through Genevestigator tool, as shown in Figure 7. Taking into the account of NDs conditions, the scale "low" is denoted =3 out of 5 experiments were represented with black color. Similarly for "Medium" =2 and for "High" =1 out of 5 black color representations were considered in ND. Likewise in stress conditions =8 out of 14 were termed as "Low, " whereas =5 was considered as Medium and =2 as High respectively (the black color box indicates that gene is not expressed in respective condition). Similar analysis was also performed for other stress conditions and the results are summarized in Table 2.

Expression and Functional Analysis of Transcript Encoding Key Characterized Enzymes Involved in ND and Stress
ND conditions generate adverse environmental growth conditions for plant, which results in inducing the degradative processes. During such conditions/treatments, plant cells acquire ability to respond to these diverse stresses by the means of various flexible and balanced stress response networks. These networks involves various metabolites such as reductionoxidation (redox), stress hormones, signaling pathways, reactive oxygen species (ROS), starch biosynthesis and growth regulators, as well as calcium and protein kinase cascades etc. (Gill and Tuteja, 2010;Foyer et al., 2016;Gill et al., 2016;Batra et al., 2017). In our study, 16 transcripts were found to encode peroxidase family enzymes (EC 1.11.1.7), which is an oxidoreductase that uses several inorganic and organic substances as hydrogen donors in the presence of H 2 O 2 (Dawson, 1988;Siegel, 1993).
Generally, peroxidase enzymes are located in the cell wall of plant cells and act as a crucial enzymatic system involved ROS generation (Bela et al., 2015;Camejo et al., 2016). High expression of peroxidase family enzymes indicates that it used soil inorganic and organic substances to provide the required energy to plant during low nutrient treatment of wheat root (Francoz et al., 2015). Further analysis suggested that 3 transcripts encode plant cell wall enzyme xyloglucan endotransglucosylase/hydrolases (XTH), which is capable of joining xyloglucan chains to oligosaccharides. XTH (EC 2.4.1.207) is involved in the alteration of the load-bearing cellwall components and these properties make it very important in the regulation of growth and development of shoot and root ( Table 2; Maris et al., 2011;Hara et al., 2013). Expression of XTH transcripts during ND condition affect the growth of wheat root as the earlier study by Maris et al. in Arabidopsis root also corroborate with our finding (Maris et al., 2009). The up-regulated expression of the transcript, code Plant Potassium Transporter (HKT1) indicates the transportation of the K in plant cells is required for a wide variety of functions like maintenance of electrical potential gradients across cell membranes, generation of turgor and activation of numerous enzymes. Most of these functions directly depend on the transportation of potassium concentration across membranebound K transport proteins (Britto and Kronzucker, 2008). Interestingly, 2 transcripts expression encode cold responsive protein that indicates that Wlt10 gene of the wheat root is activated during ND condition, which enables root to survive at freezing temperatures (Thomashow, 1998;Gupta et al., 2015a). The expression of S-adenosylmethionine decarboxylase proenzyme (SAMDC, EC 4.1.1.50), as a key rate-limiting enzyme in the polyamine (PA) biosynthesis is required for plant growth and development (Hu et al., 2005). Moreover, PA works as a metabolite and play protective as other defensive roles in response to stress (Liu et al., 2015). Three transcripts FIGURE 6 | Tissue specific potential gene expression of 55 NRGs against selected root samples. The expression potential of the genes in a tissue of wheat is measured in term of percentages value, which is shown in color key.
encode nonspecific serine/threonine protein kinase enzyme (EC 2.7.11.1) which belongs to the family of transferases (Prasad et al., 2013;Gupta et al., 2015b). Generally, these enzymes are involved in plant response to abiotic stresses and abscisic acid (ABA)-dependent plant development through phosphorylation process (Kulik et al., 2011). Expressions of Late embryogenesis abundant protein (LEA) have been linked to the survival of the wheat root under the intervals of ND and stress through defensive enzymatic function (Gupta et al., 2016a). In spite of decades of effort, the molecular-level mechanisms defining this protective function remain unknown for this protein so far (Amara et al., 2014). Moreover, some other transcripts encode characterized proteins/enzymes such as pollen allergen, cytochrome P450, vacuolar invertase, glutamine-dependent asparagine synthetase, MYB-related protein, leucine-rich repeat protein, tonoplast intrinsic protein, glutamine synthetase, and histone H2A protein play various important roles in ND and stress conditions.

Expression and Functional Annotation of Transcript Encoding Key Uncharacterized Enzymes Involved in ND and Stress
Functional analysis of uncharacterized proteins was performed on the basis of similarity search against protein database, i.e., PDB, NCBI, and UniProt as well as conserved domainbased protein family annotations (Gupta et al., 2016b). Thus, we have annotated and identified the functions for 106 uncharacterized proteins expressed during different ND. The uncharacterized enzymes with similar sequence to characterize proteins tend to have a similar function, which is already explained characterized enzyme section. The functional analysis of remaining uncharacterized proteins/enzymes, expressed moderately or high during ND and stress conditions are described (    Result of Compendium Wide analysis provides the expression level of a probe across selected 109 samples of root in response different stress. Level of gene expression measured in specified experimental sets termed as "Low," "Medium," and "High" as shown in Figure 7. Taking into the account of ND conditions, the scale "Low" is denoted ≥3 out of 5 experiments were represented with black color. Similarly for "Medium" ≥2 and for "High" ≤1 out of 5 black color representations were considered in ND. Likewise in stress conditions ≥8 out of 14 were termed as "Low," whereas ≥5 was considered as Medium and ≤2 as High respectively (the black color box indicates that gene is not expressed in respective condition).
this enzyme has been characterized through microarray analysis and various studies showed enhanced expression of its transcripts (Uhde-Stone et al., 2003;Misson et al., 2005;Morcuende et al., 2007). Similarly, in our study, up and down-regulated expression of 3 transcripts was observed, which encodes GPX-PDE under the P deficiency while in the N deficiency down-regulated 2 transcripts encodes GPX-PDE. It plays important role in hair growth and development and in the P deficiency induced phospholipid degradation pathway in wheat roots and also found in Lupinus albus (Cheng et al., 2011). AWPM-19-like membrane family protein encoded by 2 transcripts in our study plays important roles in plant development and various stress , respectively as shown in the color bar. Similarity search was done using subsets of individual stress and nutrient deficiency conditions such as heat, cold, drought/dehydration, salt, submergence, and shift from aerobic to anaerobic germination cold and drought. The expression data were obtained using Genevestigator (Zimmermann et al., 2008).
responses . It has a similar sequence of LEA protein, which is induced by ABA during freezing tolerance (Koike et al., 1997).

DISCUSSION
Transcriptomic and Meta-Analysis of Array Data of T. aestivum cv. root Since past one and half decade several transcriptome analysis were performed in different plant species in response to diverse micronutrient deficiency to mine the important genes along with their role in different metabolic pathways (Wasaki et al., 2003(Wasaki et al., , 2006Li et al., 2010a,b;Cai et al., 2012;Ma et al., 2012;Park et al., 2012;Takehisa et al., 2013Takehisa et al., , 2015Yu et al., 2016;Forieri et al., 2017;Gupta et al., 2017). In continuation, the current study reports 398 differentially expressed genes through a comparative transcriptomic study of micronutrient deficient samples grown in soil having environment growth conditions. The study identified 58 probes (55 genes) are commonly expressed under at least two ND conditions (Figure 3). The differentially expressed genes were identified through normalization, QC, QA of array data (Raychaudhuri et al., 2000;Bolstad et al., 2003). Further, the GO-based analysis of differentially expressed NRGs mimics their association with distinct GO terms. Through the GO term it can be interpreted that; one NRG can participate in different biological functions and vice-versa (Gupta et al., 2016a). Accordingly, the study identified 291 NRGs associated with different GO terms indicating their roles in regulation of different reaction mechanisms in response to macronutrients stress (Table S1). Taking an example: the peroxidase precursor gene, which generally participate in peroxidase activity (GO: 0004601) also plays an important role in peroxidase reaction (GO: 0006804) in cytoplasmic membrane-bounded vesicle (GO: 0016023) of plant cells. Figure 4 represents different significantly enriched GO terms (P-value = 0.05) obtained through functional characterization of NRGs involved in biological processes and molecular functions. The major affected biological processes are cellular biogenic amine biosynthesis, inorganic anion transport, response to stress, response to stimulus and cellular glucan metabolism. 17 NRGs involved in the biosynthesis of cellular amines and poly amines, are found to play a crucial role in different signaling networks and metabolic process of plants (Hussain et al., 2011;Berkowitz et al., 2016). Five NRGs which are involved in the movement of inorganic anions i.e., in or out of the cell or within a cell, or between plant cells are being observed as transporters which indicates their important functions during P, N, and K deficiency (Sharma, 2016;Dong et al., 2017;Forieri et al., 2017). Forty NRGs are participating in response to stress and stimulus indicates their involvement in different metabolic and signaling pathways. Eleven NRGs are found to be involved in peroxidase reaction and oxidative stress. The biosynthesis of glucan/carbohydrate and its metabolism are the key processes mediating plant responses to abiotic stresses. Under such stress, plants generally remobilize glucan/carbohydrate to supply energy and carbon at times when photosynthesis process may become limited (Thalmann and Santelia, 2017). Nicotianamine synthase activity (13 NRGs) and metal ion binding activity (16 NRGs) are the main sub category under molecular functions. Nicotianamine is an intermediate for the biosynthesis of mugineic acidfamily phytosiderophores in the crop plant and work as key substance for iron metabolism. Nicotianamine synthase catalyzes the formation of nicotianamine from S-adenosylmethionine (Higuchi et al., 1996). Heme and tetrapyrrole binding may occur during the transportation of the macronutrients for plant organs. Similar GO based analysis was also reported in earlier studies (Ma et al., 2012;Takehisa et al., 2013Takehisa et al., , 2015Yu et al., 2016;Shi et al., 2017;Zhu et al., 2017), using the ontological terms to describe the function of differentially expressed genes. The NRGs were mapped through MapMan scheme ( Figure 5) and results indicate the number of significant genes, which are associated with different metabolic pathways in macronutrients deficiency. The study was performed with an assumption that the cellular processes and metabolic pathways involved in supply of water, nutrients factors are in general common for availability and deficiency of K, N, P, Mg, and S. Interestingly, the metabolic pathway and cellular process are also involved in amino acid metabolism, protein synthesis, metal binding, hormone synthesis, lipid synthesis and its degradation regulated by similar genes under most of the ND cases. Different cellular metabolic processes, i.e., DNA synthesis/repair, RNA regulation and transcription of transcription factors and related proteins, nitrogen/ammonia metabolism, abiotic stress and cell wall synthesis have been found to be affected under the P deficiency. Similarly Photosystem process and related light reaction have been found to be affected in the N deficiency (Table S4). Further, a meta-analysis of each transcript was performed to identify the differential expression in response to particular ND and stress condition. This analysis summarized the results of single experiment analysis and compendium wide analysis by taking consideration of all 109 wheat root specific microarray data samples using Genevestigator tool (Thimm et al., 2004). Metaanalysis checks the expression potential of each NRG under Mg, N, P, K, and S deficiency and other stress such as dehydration, drought and salt etc. (Figure 7). Moreover, the hierarchical clustering based tissue-specific percentage expression potential of each transcript explains the tissue-specific role of candidate genes ( Figure S4). Further functional role of different NRGs encoding uncharacterized proteins/enzymes ( Table 2) were identified and explained. Thereafter, the functional roles of 35 transcripts out of 55 transcripts were annotated through standard protocols (Gupta et al., 2016b). Remaining 20 transcripts are yet to be functionally characterized in wheat varieties, through in-vitro validation studies.

Comparative Analysis with Transcriptome Studies Performed In-vitro in Model Plants
Recently, Ruan et al. (2015) performed in-vitro transcriptome studies in seedlings of two wheat genotypes i.e., low-K tolerant "Tongzhou916" and low-K susceptible "Shiluan02-1" using Affymetrix R GeneChip R analyzer in response to low K deficiency. They found 2713 and 2485 differentially expressed genes in Tongzhou916 and Shiluan02-1, respectively. Among 2713 genes of Tongzhou916, 1321 genes are up-regulated and 1392 genes are down-regulated; whereas out of 2485 genes of Shiluan02-1, 1177 genes are up-regulated and 1308 genes are down-regulated (Ruan et al., 2015). While comparing with the present study involving differentially expressed genes under K deficiency, only 42 and 3 genes are found to be upregulated and down-regulated respectively. Comparative analysis showed that probes i.e., Ta.21127.1.S1_at, Ta.1207.1.S1_x_at, or Ta.1207.1.S1_s_at encoded for Nitrate/Chlorate transporter (Gene ID: CA607301), and OPR1 12-oxophytodienoate reductase (Gene ID: CA650490), which were commonly expressed in both the study. It indicates their role as the K deficiency affected genes, in different stages of wheat growth. Despite this fact, we observed commonly affected metabolic pathways i.e., Jasmonic acid and related, ethylene-related, transporter of the peptide, oligo-peptides, and potassium are affected in both soil-grown and in vitro grown wheat roots and seedlings under the K deficiency. Moreover, two earlier studies by Armengaud et al. and Ma et al. in model plant Arabidopsis and Rice under the K deficiency respectively, also corroborates the role of jasmonic acid signaling (Armengaud et al., 2004;Ma et al., 2012). In-vitro transcriptomic studies performed (Li et al., 2009(Li et al., , 2010b in the root of rice under low phosphorus stress shows that the differentially expressed genes also participate in different functional category i.e., metabolism, secondary metabolism signaling, stress DNA synthesis, transcription, transport, lipids metabolism, the P transport. Similarly, in our study, a number of genes are found to be participating in above said functions, which infers that similar reaction mechanism might be occurring under the P deficiency whether they are soil-grown or in vitro. It can also be inferred that the less number of gene expressions occurring in soil-grown sample in comparison to in vitro grown sample, may be due to the relatively uncontrolled growth conditions in soil grown study in comparison to the in vitro grown study.

CONCLUSIONS
The study reports transcriptome analysis using Affymetrix R GeneChip R Wheat Genome Array to illustrate the changes in gene expression profiles in soil grown nutrients deficient roots of wheat. Microarray data analysis facilitated the discovery of a large number genes expressed in wheat roots grown with limited K, N, P, Mg, and S nutrient controlled plots. The Hierarchical clustering and PCA were carried out to ensure the overall quality of the normalized samples ( Figure 1A). Further, the differentially expressed 398 NRGs were identified through statistical analysis. The NRGs responsible for more than one ND were identified and the functional mechanism for each of such candidate was also explored. The SEA provides a GO based classification of NRGs and hence classified them accordingly into respective biological processes, cellular components, and molecular functions (Figure 4). MapMan based analysis illustrates the specific role of each gene in different metabolic pathways. Different metabolic pathways (Table 1) occurring in wheat root got affected due to multiple NDs which deciphers the role of an individual gene in multiple pathways. In addition, a meta-analysis of NRGs was carried out to provide a concrete understanding of gene expression potential with respect to selected root specific 109 samples (Table S5), generated under biotic/abiotic stress. Further, the study describes the major role of characterized/uncharacterized protein/enzyme under nutrients deficiency and abiotic stress condition. Finally, the variation in genes expression level in vitro and soil-grown environment was discussed for the entire set of candidate genes. Overall, the study illustrates set of known genes along with a set of novel candidate genes expressed under nutrients deficiency, which opens a new arena for plant biotechnologist to study and identify their specific functions.

ETHICS STATEMENT
This research does not carry out any experiment on human and animals. All data used in this in silico work were collected from the open sources. Hence, authors declare that there is no compliance with ethical standards.

AUTHOR CONTRIBUTIONS
Study was conceived and designed by: SG and PV; Work was performed and analyzed by SG, BY, and PV; Manuscript prepared by: SG, UR, SF, and PV.

ACKNOWLEDGMENTS
The authors are grateful to Gruen, A. and his group members at Plant Biology and Crop Science Group, Rothamsted Research, Harpenden, UK for providing their microarray data to carry out this work. Authors are also thankful to Indian Institute of Information Technology, Allahabad for providing the required infrastructure and computational facilities to complete this work.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017. 01025/full#supplementary-material Figure S1 | Boxplots of gene expression levels in the wheat in response to different nutrients deficiency. The left y-axis shows the normalized intensity values and x-axis represents the three replicates (R1, R2, and R3) for Control, Potassium (K), Magnesium (Mg), Nitrogen (N), Phosphorus (P), and Sulfur (S). Figure S2 | Hierarchical clustering of all array samples data. The array data were examined by the hierarchical clustering method with two views of the clusters: group of treatment view and gene tree view. The group of tree consists of 6 experimental groups i.e. Control, Potassium deficiency (ND_K), Magnesium deficiency (ND_Mg), Phosphorus deficiency (ND_P), Nitrogen deficiency (ND_N), and Sulfur deficiency (ND_S). Figure S3 | Percentage expression potential of each differentially expressed probes in 109 roots specific samples of wheat available in different expression databases. Figure S4 | Tissue specific potential gene expression of all probe sets against selected root samples. The expression potential of the genes in a tissue sample of wheat is measured in term of percentages value, which is shown in color key. Figure S5 | Heat map diagram for all probe sets while comparing with root specific samples. The green and red colors specify down-regulation (log2 [2.5]) and up-regulation (log2 [-2.5]) respectively as shown in the color bar. Similarity search was done using subsets of individual stress and nutrient deficiency conditions such as heat, cold, drought/dehydration, salt, submergence, and shift from aerobic to anaerobic germination cold and drought. The expression data were obtained using Genevestigator (Zimmermann et al., 2008).
Table S1 | List of NRGs expreesed in wheat root under different nutrients deficiency conditions and their fold change value and other annotated explanation. Table S2 | Probes with annoteted Gene ID and description, playing role in atleast two nutrients deficiency conditions. Blank box represent that particular gene dose not expressed in respective deficiency.  Table S4 | Identified hypergeometric P-values of genes associated with different metabolic pathways under all Nutriant dificiencies. The P value is calculated to identify the probability of perticular gene/gene group is over-represented within a functional bin at a rate higher than chance expectation. The gene/gene group has P value <0.05 was considered for futher analysis.