Introduction

Nitrogen (N) is an essential nutrient for plant growth and it is an integral component of several plant compounds. N fertilizers are used excessively to increase crop yield1. In potato, application of high dose of nitrogen is very common in fields to achieve high tuber yield. For example, in India, on an average 180–250 kg N/ha is applied in potato fields to get 30–40 tonnes per hectare tuber yield2. Of total N application, nearly half amount of N is utilized by potato crop2, whereas excessive N increases environmental pollution, deteriorates soil health and water quality, and also leads to higher cost of production3. Hence, understanding about genes involved in N metabolism in potato, especially low N stress, is essential to modulate plants to optimize the applied N. Several agronomic and soil-based approaches have been demonstrated in crops on N fertilization in plants including potato4,5,6. Unlike the model plant, Arabidopsis thaliana, and model crops rice and maize where rich source of genomics-based information is available on N metabolism research3,7, knowledge about potential genes involved in N deficiency and sufficiency in potato is very limited. Hence, identifying genes underlying N deficiency is essential to strengthen future research on N metabolism in potato.

Advancement has been made in N metabolism research in plants1. The N metabolic pathways are controlled by various genes, which express differently under varied N regimes. The key enzymes/genes involved in N metabolism (uptake, translocation, assimilation/utilization and remobilization) pathways are nitrate transporter (NRT), ammonium transporter, nitrate reductase, nitrite reductase, glutamine synthetase, glutamate synthase, glutamate dehydrogenase, asparagine synthetase, aspartate aminotransferase, alanine aminotransferase; and many other genes/regulatory molecules alter N pathways and determine gene expression as well as phenotypes3. Previous studies have investigated roles of nitrate transporters, signal transduction molecules and transcription factors in N metabolism in plants3. Moreover, RNA-sequencing technology has been proven to discover numerous genes/factors involved in N gene networks in several crops for multiple traits such as role of N starvation in rice8,9; identifying nitrogen deficiency genes in cucumber10; effect of N limitation condition in maize7, N use efficiency genes in Brassica juncea11, N responsive regulatory elements in potato12, N utilization genes in tea13 and so on.

Recently, we have proposed an integrated approach of breeding, genomics and physiology to improve NUE in potato by translating knowledge from other plants like Arabidopsis thaliana, rice, maize and wheat14. Furthermore, we have also demonstrated precision phenotyping of potato in aeroponic culture with accurate macro- and micro-nutrients supplies including minimum N dose while maintaining tuber yield15. With the availability of potato genome sequences16 and increasing sequencing facilities at low cost, it is now possible to unveil genes underlying N metabolism in potato. To our knowledge, in potato Gálvez et al.12 have identified N responsive genes and regulatory motifs in field-grown plants using RNA sequencing technology.

Hence, the aim of this study was to discover genes/regulatory elements associated with N deficiency (low N) versus sufficient N (high N, control) in potato plants grown in aeroponic culture under controlled conditions by RNA-sequence-based transcriptomes analysis. Differentially expressed genes (DEGs), heat map, Venn diagram, scatter plot, volcano plots, GO characterization, KEGG pathways analysis, motifs identification and potential genes to be involved in N metabolism pathways in potato were analyzed. Further, selected DEGs were validated by real-time qPCR analysis. Our results usher an insight on enrichment of genes in potato under N deficiency, which could be utilized for N metabolism related research in potato.

Results

Transcriptome data generation and processing

A popular Indian potato variety Kufri Jyoti was grown in three replicates in aeroponic culture under controlled conditions with varied N supplies (low N: 0.2 milli molar N; and high N: 4 milli molar N) (Fig. 1). To obtain an overview of genes and regulatory elements involved in N metabolism in potato, especially low N stress versus high N (control), plant tissues (shoots, roots and stolons) were analyzed by RNA-sequencing. High quality data (QV > 20) ranging between 3.25 to 4.93 Gb per sample were generated using Illumina technology. The quality reads of all six tissues showed 83.60-86.50% mapping to the reference potato genome (Table S1 and Fig. S1).

Figure 1
figure 1

Phenotypes of potato plants grown in aeroponic culture with low N (0.2 milli Molar) and high N (4 milli Molar) supply showing shoots (a,c) and roots (b,d) biomass under different N regimes. (a,b: low N; c,d: high N).

Full size image

Differentially expressed genes (DEGs) analysis

Individual transcriptomes were assembled and then DEGs were identified using cufflinks and cuffdiff, respectively. DEGs were analyzed separately for shoot, root and stolon tissues between low N stress and high N (control) (Table S2). Complete list of DEGs is provided in supplementary datasets (Table S3: shoots; Table S4: roots; and Table S5: stolons). Significant DEGs were identified based on the statistical significance (p ≤ 0.05) and ≥ 2 log2 fold change (FC) for up-regulated, and ≤ −2 log2 FC values for down-regulated genes. In shoots, a total of 19730 genes were differentially expressed, of which significant DEGs were 761 (up-regulated) and 280 (down-regulated); whereas, of 20736 DEGs in roots, significant DEGs were 572 (up-regulated) and 292 (down-regulated). In stolons, of 21494 DEG, 688 and 230 DEGs were significantly up-regulated and down-regulated, respectively. In addition, many DEGs were exclusively expressed in the tissues. Heat maps of the top 50 DEGs are shown in Fig. 2 (shoots), Fig. 3 (roots), and Fig. 4 (stolons). Where, red and green colours represent up-regulated and down-regulated DEGs, respectively in low N compared to high N. Scatter plot (Fig. S2) and volcano plot (Fig. S3) show significant up-regulated and down-regulated DEGs in shoot, root and stolon tissues.

Figure 2
figure 2

Heat maps of top 50 differentially expressed genes in shoots of potato plants grown in aeroponic culture with low N (0.2 mM) and high N (4 mM, control) supply. In the heat maps, each horizontal line refers to a gene. Relatively up-regulated genes are shown in red colour, whereas down-regulated genes are shown in green colour under low N stress compared to control (high N).

Full size image
Figure 3
figure 3

Heat maps of top 50 differentially expressed genes in roots of potato plants grown in aeroponic culture with low N (0.2 mM) and high N (4 mM, control) supply. In the heat maps, each horizontal line refers to a gene. Relatively up-regulated genes are shown in red colour, whereas down-regulated genes are shown in green colour under low N stress compared to control (high N).

Full size image
Figure 4
figure 4

Heat maps of top 50 differentially expressed genes in stolons of potato plants grown in aeroponic culture with low N (0.2 mM) and high N (4 mM, control) supply. In the heat maps, each horizontal line refers to a gene. Relatively up-regulated genes are shown in red colour, whereas down-regulated genes are shown in green colour under low N stress compared to control (high N).

Full size image

Venn diagram analysis showed that only eight up-regulated and four down-regulated DEGs were common in shoot, root and stolon tissues (Fig. 5). In general, gene categories such as glutaredoxins, transcription factors (Myb, WRKY, LOB domain and heat shock transcription factors), nitrate transporters, aquaporin and ABC transporter, sulfate/bicarbonate/oxalate exchanger and transporter sat-1, urea active transporter, sodium/proline symporter, ferric-chelate reductase, miraculin, and methylketone synthase Ib were differentially expressed (up-/down-regulation) in all three tissues. Particularly, genes with known function namely extension (PGSC0003DMG400000783), ferric-chelate reductase (PGSC0003DMG400000184), GTP binding protein (PGSC0003DMG402012350), UPF0497 membrane protein (PGSC0003DMG400027047) were up-regulated, methylketone synthase Ib (PGSC0003DMG400025825), CLC-Nt2 protein (PGSC0003DMG400020215), stem-specific protein TSJT1 (PGSC0003DMG400012822) genes were down regulated in all three tissues. Besides, many genes showed tissue-specific response under N stress in potato (supplementary dataset: Table S15).

Figure 5
figure 5

Venn diagrams of up-regulated and down-regulated DEGs showing distribution of tissues-specific and tissues-independent genes in potato plants (shoot, root and stolon).

Full size image

Identification of potential genes under N stress

Selected top 20 significantly up-regulated ( ≥ 2.0 log2 FC; p < 0.05) and down-regulated ( ≤ −2.0 log2 FC; p < 0.05) DEGs under N stress are summarized in Table 1 (shoots), Table 2 (roots) and Table 3 (stolons). These DEGs recorded high levels of gene expression in all the tissues such as shoots (up-regulated between 4.61 to 2.72 log2 FC, and down-regulated between −4.31 to −2.35 log2 FC); roots (up-regulated: 6.96 to 3.12 log2 FC, and down-regulated: −5.67 to −2.83 log2 FC); stolons (up-regulated: 10.07 to 5.70 log2 FC, and down-regulated: −8.93 to −3.22 log2 FC).

Table 1 Selected top 20 differentially expressed genes (DEG) (p < 0.05) in potato shoots under low N stress versus high N (control).
Full size table
Table 2 Selected top 20 differentially expressed genes (DEG) (p < 0.05) in potato roots under low N stress versus high N (control).
Full size table
Table 3 Selected top 20 differentially expressed genes (DEG) (p < 0.05) in potato stolons under low N stress versus high N (control).
Full size table

Analysis of the DEGs across the tissues showed that glutaredoxin/glutaredoxin family genes were over expressed ( > 3 log2 FC) in both shoots and roots, and played an important role in N stress tolerance in potato. High-affinity nitrate transporter genes were also highly expressed in roots under N stress. Similarly, in stolon tissues ion transporters (e.g. urea active transporter, sodium/proline symporter) were highly up-regulated (>5 log2 FC), nitrate transporter and ion homeostasis genes like ferric-chelate reductase were down-regulated (< −3 log2 FC), and metallocarboxypeptidase inhibitor was over expressed (>5 log2 FC). In roots, ferric chelate reductase was up-regulated; however, sulfate/bicarbonate/ oxalate exchanger and transporter sat-1, nitrate transporter were under expressed (< −3 log2 FC). Other important genes like phospho-lipid related genes (e.g. tartrate-resistant acid phosphatase type 5, glycerophosphodiester phosphodiesterase) were under expressed in shoots, whereas GDSL esterase/lipase was down regulated (< −3 log2 FC) in roots; and protein phosphatase 2 C) were up-regulated in roots under N stress.

Transcription factors (TFs) play key roles in plant adaptation to stress metabolism. TFs (e.g. TAGL11) was up-regulated (>3 log2 FC) in roots; whereas, WRKY transcription factor 16 was under expressed (< −3 log2 FC) in shoots under N stress. TF like heat stress transcription factor A-6b was up-regulated (>5 log2 FC) in stolons. Besides, DNA binding proteins (e.g. BTB/POZ domain-, LOB domain-containing protein, F-box) were down-regulated (< −3 log2 FC) in stolons. Amino acid is an essential component of protein and its role was observed more in roots. The amino acid metabolism genes (e.g. malate synthase) were up-regulated (>3 log2 FC) in roots; whereas, gamma-glutamyl-gamma-aminobutyrate hydrolase was down-regulated (< −3 log2 FC) in roots. However, gene 2-oxoglutarate-dependent dioxygenase showed both up-regulation and down-regulation in roots under N stress.

Starch is the major food source of potato. A few classes of genes were only observed in stolon tissues. Starch/sugar metabolism related genes like glucose-6-phosphate/phosphate translocator 2 were highly up-regulated (>5 log2 FC) in stolons. Inhibitor proteins like cysteine protease inhibitor 1, tuber storage protein patatin-04/09, ethylene biosynthesis related genes ACC oxidase were highly over expressed (>5 log2 FC) in stolons. Terpenoids biosynthesis genes such as sesquiterpene synthase were under expressed (< −3 log2 FC) in stolons. Stress responsive genes like dehydration-responsive protein RD22, hydroxyproline-rich glycoprotein were under expressed (< −3 log2 FC) in stolons.

Gene ontology (GO) characterization

All DEGs were functionally assigned with the GO terms, in which the molecular function GO domain was found to be highest in roots, shoots and stolons. Complete list and details are provided in supplementary datasets (shoots: Table S6; roots: Table S7 and stolons: Table S8) and result is summarized in Table S9. In shoots, 25142 molecular function followed by 18805 biological process and 17921 cellular component GO terms were observed. Whereas, in roots 26132 molecular function followed by 19402 biological process and 18210 cellular component GO terms were found. In stolon also, molecular function (27286) was predominant followed by biological process (19753) and cellular component (18578). Overall a few GO terms such as cell, cell part, membrane, membrane part, catalytic activity, binding, metabolic process, cellular process were highly enriched in both up-regulated and down-regulated DEGs under low N and high N in the tissues. WEGO plot showing the GO terms of DEGs in different tissues are depicted in supplementary figures (Fig. S4: shoots, Fig. S5: roots, and Fig. S6: stolons). A few GO terms were observed exclusively in certain tissues. Analysis across the tissues showed that GO terms such as molecular carrier activity (molecular function) and rhythmic process (biological process) were observed only in down-regulated shoots, whereas GO term locomotion (biological process) was observed only in down-regulated roots. None of the gene term was found exclusively only in up-regulated genes. This indicates role of various genes in N stress metabolism in potato.

KEGG pathways analysis

The DEGs were processed in KAAS and classified them into 24 KEGG functional pathways categories, which include annotation of 5385, 5572 and 5594 DEGs in shoots, roots and stolons, respectively (supplementary datasets-Table S10, Table S11, Table S12; and supplementary file: Table S13). Overall, maximum annotated genes were found to be associated with KEGG pathways like signal transduction; translation; carbohydrate metabolism; transport and catabolism; folding, sorting and degradation; amino acid metabolism; energy metabolism; lipid metabolism; and environmental adaptation etc. The KEGG pathways such as lipid metabolism, metabolism of other amino acids, glycan biosynthesis and metabolism, metabolism of terpenoids and polyketides, biosynthesis of other secondary metabolites were highest in stolons than other tissues. Interestingly, carbohydrate metabolism, amino acid metabolism, signal transduction, membrane transport, and environmental adaptation KEGG pathways were highest in roots; whereas, energy metabolism and metabolism of cofactors and vitamins pathways were highest in shoots. KEGG pathways of N metabolism related genes are shown in Fig. S10. This signifies importance of various gene networks of N metabolism in potato.

Conserved motifs analysis

Conserved motifs were analyzed in highly up-regulated (≥3 log2 FC) and down-regulated (≤−3 log2 FC) DEGs with known functions in the tissues using MEME (version 5.1.0) software. In shoots, six motifs were found with significant E values and all six were over-represented in glutaredoxin genes and some motifs in other genes (e.g. motif 1: TCCCATTTGGAWGTGTATCGAGCTCGTAAAYTATAGGGTTTGCT; motif 2: TGGTGAAAAT CACCACCGCRCTTGDTGCTCCCAACTTCATCACCATATCC; and motif 3: CTCATTAGCAC CWCCAACTAACTCTTTCCCTATAAATATTGCTGG). In roots, only motif 1 (AAGAGARRNA RWGRAAGAARAARVARA) and motif 2 (TTGGTGTTGSDVHTGGKGGA) were observed with significant E values and both motifs were over-represented in the genes studied. In stolon tissues, all six motifs were found with significant E value (e.g. motif 1: ATCYYTCYMWCWMY TCTCTWYTWTWTTCCWBATCTMCMTSWWMRCYTCTC; motif 2: TGRTKGCTYYWY SASMCATGCCKMMMCAAGAARCAKTRSWRARAWTKRAG; and motif 3: YTCYCCWT YNMCWMCYCHWYCTCCMCTYM) and motif 3 was over-represented in the genes studied, whereas metallocarboxypeptidase genes contained all six motifs. Motif details are shown in supplementary files (Fig. S7: shoot; Fig. S8: root; and Fig. S9: stolon).

Validation by RT-qPCR analysis

Selected 12 genes (2 genes from each tissue of root, shoot, and stolon tissues of both low N and high N treatments) were validated by RT-qPCR analysis. Results were found in agreement with the RNA-seq-based gene expression pattern with minor variations in the log2 FC values (Table S14).

Plant biomass, total chlorophyll and total N content analysis

Phenotypic traits on per plant basis were estimated and showed higher biomass with high N supply. Traits such as plant height, leaf length, leaf width, tuber yield, tube number, root dry weight, shoot dry weight, total chlorophyll content and total N content were found significantly higher in high N than low N supplied plants (Fig. 6).

Figure 6
figure 6

Plant biomass, tuber traits, total chlorophyll and total N content in potato plants of cv. Kufri Jyoti grown in aeroponic culture with contrasting N regimes till 90 days crop duration. DW: Dry weight, FW: Fresh weight, LN: Low N (0.2 mM N), HN: High N (4 mM N), Values on bar are mean ± standard error; * and ** indicate statistical significance at P ≤ 0.05 and P ≤ 0.01, respectively.

Full size image

Discussion

The present study provides an overview enrichment of genes associated with N metabolism in potato to provide adaptive strategies against N stress management through genomics intervention. Genes were analyzed in root, shoot and stolon tissues of potato plants grown in aeroponic culture for N metabolism such as uptake and utilization/assimilation and remobilization. This study revealed significant increase in plant biomass (plant height, leaf length, leaf width, shoot dry weight and root dry weight), fresh tuber yield, tuber number, total chlorophyll content and total N content with increasing N dose. Our results also concur with earlier findings on trends of variation in phenotypic traits under low N stress in potato15, for example, Gálvez et al.12 reported increase in plant biomass, tuber yield, leaf chlorophyll content on increasing N supply in potato.

RNA-seq approach has been applied in many crops for multiple traits to identify gene networks. In this study, we identified maximum number of up-regulated DEGs in shoots (761) followed by stolons (688) and roots (572), whereas maximum number of down-regulated DEGs were observed in roots (292), shoots (280) and stolons (230). In terms of fold change gene expression values, maximum up-regulation (10.073 log2 FC) and down-regulated (−8.933 log2 FC) were observed in stolon tissues followed by roots (6.969 and −5.674 log2 FC) and shoots (4.619 and −4.317 log2 FC). Besides, Venn diagram analysis unveiled tissue specific and common genes. The GO characterization shows predominance of molecular function in the tissues and some exclusive GO terms were also noticed. The GO terms like cell, cell part, membrane, membrane part, catalytic activity, binding, metabolic process, and cellular process were highly enriched in all the tissues. KEGG pathways analysis shows that the highest numbers of annotated genes belong to signal transduction, translation, and carbohydrate metabolism, and some genes in N metabolism. In motifs analysis, we identified more conserved motifs in shoots and stolons than roots, as reported earlier patatin gene by Gálvez et al.12 Further, RT-qPCR analysis of some selected genes validated the RNA-seq results. In potato, one study identified motifs associated with N genes in field-grown plants under with- and without-N supplies through RNA-seq approach12. Our study provides deeper insights into N stress responsive genes in potato, where plants were grown in aeroponic culture under controlled conditions with more precision on nutrients supply and plant growth as well as tuberization.

We identified a number of DEGs under low N stress in different potato tissues (shoot, root and stolon). Of which a few potential DEGs having higher gene expression (>2 log2 FC up-regulation or <−2 log2 FC down-regulation, as given in Tables 1, 2 and 3) are discussed here. Interestingly, glutaredoxins were the most highly up-regulated DEGs in shoots (also up-regulated in roots) during low N stress in potato compared to high N, and these were observed in all three tissues but more numbers of genes in shoots. A large number of glutaredoxin family proteins exist in plant, which are a group of oxidoreductases that control reactive oxygen species. Glutaredoxins families perform diverse functions in plants such as development, defense response, abiotic stress response (salt, drought, heavy metal, low nitrogen and temperature), redox signalling, flower development, hormonal regulation, ion homeostasis and adaptation to environments17. Another study shows differential expression of glutaredoxins in response to biotic and abiotic stresses during suberization in potato18, and overexpression of the CC-type glutaredoxin OsGRX6 affects hormone signalling and nitrogen status in rice19. Our study suggests that changes in expression of glutaredoxins in shoots play key roles in adaptation to N deficiency in potato.

In this study, transcription factors (Myb, WRKY, TAGL11, F-box, RING-H2 finger protein ATL18, heat shock TFs, heat stress transcription factor A-6b, BTB/POZ domain- and LOB domain-containing proteins) played crucial roles in regulation of genes during N deficiency in potato. Among them Myb, WRKY, LOB domain and heat shock TFs were common in all three tissues, whereas TFs like BTB/POZ domain and RING-H2 finger protein ATL18 were found in roots and stolons, and F-box TF was observed in shoot and stolon, and TAGL11 (root) and heat stress transcription factor A-6b (stolon) were tissue specific. The Myb and WRKY TFs were down-regulated in shoots under N stress; whereas, TAGL11 TF was up-regulated in roots under N stress. TFs are one of the key regulators in plant metabolism and constitute an important part of gene networks and signalling molecules for biotic and abiotic stress responses20. A study in maize shows that transcription factors GLK5, MADS64 and bZIP108, asparagine synthetase, protein kinase and a protein phosphatise are involved in N limitation conditions21. Besides, WRKY is one of the largest TFs gene families, which regulates signalling pathways and plays critical roles in plants in response to biotic and abiotic stresses22. The F-box proteins are encoded by a large number of gene families and function in various cellular processes including cell cycle, transcription and signalling molecules23. Furthermore, involvement of zinc finger protein 3 has been proven in salt stress and other osmotic response in Arabidopsis thaliana24. The roles of heat shock transcription factors (Hsfs) have been shown in response to various stresses, particularly adaptation to abiotic stresses. Transcriptomes and genomes have been investigated to study Hsfs in potato25. Besides, BTB/POZ domain-containing protein and LOB domain-containing protein were also highly down-regulated under low N stress in potato stolon tissues. Our study highlights that variation in the gene expression pattern of TFs and other binding proteins indicate their implication in N stress mechanism in potato.

Another gene FLOWERING LOCUS T was down-regulated under N deficiency exclusively in shoots only and implies to early plant maturity under N stress condition in potato. The FLOWERING LOCUS T protein is the main component of signalling molecules for flowering in plants, as studied in Arabidpsis thaliana, rice and tomato. Navarro et al.26 demonstrate that potato flowering and tuber development are controlled by FT-like paralogues namely StSP3D and StSP6A, which respond to different environmental cues. A study shows that CONSTANS molecules regulate StSP6A gene-driven tuberization in potato26. Another plant peptide hormones like CLAVATA3(CLV3)/EMBRYO SURROUNDING REGION-related (CLE7) family was the most down-regulated under N stress in roots only. The CLE molecules have been showed to play vital role in signalling for meristematic differentiations of roots and shoots, and play diverse roles in developmental process particularly symbiosis, parasitism and response to abiotic stress27. The varied gene expression of these genes indicates that these genes are important in potato under N deficiency.

Among the N metabolism genes, high-affinity nitrate transporters were the most down-regulated in roots under low N and suggest their crucial roles in N uptake mechanism for adaptation to N stress in potato. Overall, nitrate transporters were present in all three tissues, showing importance of N uptake in potato. Given that nitrate transporters perform very essential roles in N uptake in plants and they include mainly four gene families NRT1 (NPF), NRT2, CLC, and SLAC1/SLAH28. Nitrate transporters function for nitrogen uptake and root system architecture, protein storage, source to sink relationship, ionic balance, response to biotic and abiotic stresses, and carbon-nitrogen balance. Several reviews discussed the roles of nitrate transporters in regulation of plant growth and development and a potential target for marker or genomics-assisted breeding, genetic manipulation via genome editing or transgenic technologies29,30.

Down-regulation of aquaporin TIP1;3 gene was observed in stolons that showed its role in N stress in potato. In general, aquaporin and ABC transporter genes were found in all three tissues. The aquaporins are commonly known as water channel proteins, which regulate movement of water in plants. Aquaporins play important role in movement of solutes, small molecules and metal ions in response to biotic and abiotic stresses in plants31. Besides the urea transporters of Arabidopsis thaliana, several other transporters are known in plants like aquaporins and ABC transporters32. Furthermore in this study, genes like sulfate/bicarbonate/oxalate exchanger and transporter sat-1, urea active transporter, sodium/proline symporter, ferric-chelate reductase were also up-regulated or down-regulated in all the tissues. The genes such as sulfate/bicarbonate/oxalate exchanger and transporter sat-1 were down-regulated in roots, whereas urea active transporter, sodium/proline symporter were up-regulated in stolons under N stress. Previous study also found changes in sulphate-related genes in potato under varied N supplies12. Besides, ferric-chelate reductase and iron-regulated transporter 1 were down-regulated under N stress in potato roots. Role of ferric reductase oxidase enzyme is to reduce ferric Fe (III) into ferrous Fe (II) and thus play important roles in Fe homeostasis in plants33. Nevertheless, plant root architecture plays an important role in nutrient uptake in crops especially in potato34, where the underground plant part tuber is economically important. Our results indicate potential roles of these N metabolism associated genes especially transporters under N stress tolerance in potato.

Starch is the main component of potato tubers and has driven attentions for food as well as non-food/commercial uses. Despite numerous researches in potato biotechnology, detailed information on gene network in starch metabolism is limited35. In this study, glucose-6-phosphate/phosphate translocator 2 (GPT2) was the most up-regulated gene under N stress in stolons and it substantiates roles of carbohydrate metabolism related genes in potato tuberization. Moreover, this gene was found (up-regulation) in both root and stolon tissues, which could be associated with tuber development. Signalling role of GPT2 has been observed in seeds, seedlings and mature leaves of Arabidopsis thaliana in response to environments36. Glucose-6-phosphate is translocated via GPT during starch biosynthesis, in which inorganic phosphate is released or used as substrate during pentose phosphate pathways. To illustrate, a plastidic GPT has been purified from maize endosperm and corresponding cDNA was isolated from pea roots and potato tubers37. Another gene like glucose-1-phosphate adenylyltransferase was highly up-regulated exclusively in stolon tissues only under N stress, which is involved in starch and sucrose metabolism in potato35. Moreover, genes like malate synthase, branched-chain-amino-acid aminotransferase and 2-oxoglutarate-dependent dioxygenase were observed in only in root and stolon tissues but had varied fold change. In particular, malate synthase and branched-chain-amino-acid aminotransferase genes were highly up-regulated in roots, while 2-oxoglutarate-dependent dioxygenase gene was both up-regulated and down-regulated in roots under N stress. Role of 2-oxoglutarate-dependent dioxygenase gene SlF3HL has been confirmed recently for chilling stress tolerance in tomato38. Malate synthase plays an important role in starch synthesis in potato39, while branched-chain amino acid aminotransferase enzyme catalyzes conversion of branched-chain amino acids and α-ketoglutarate into branched chain α-keto acids and glutamate40. Thus, our study implicates that starch metabolism associated genes play vital role in potato tuberization process.

In addition, cytochrome P450 hydroxylase, gamma-glutamyl-gamma-aminobutyrate hydrolase, and flavonol synthase/flavanone 3-hydroxylase were either up- or down-regulated only in root and stolon tissues, but they were highly down-regulated in roots under N stress. The cytochrome P450 hydroxylases function in a variety of metabolic pathways in plants and involved in the jasmonic acid and ethylene signaling pathways, enhances plant resistance to biotic and abiotic stresses in soybean41 and under N stress in cucumber10. The gamma-glutamyl-gamma-aminobutyrate hydrolase catalyzes amino acid metabolism (gamma aminobutyric acid, GABA). The study shows that GABA concentration increases under various biotic and abiotic stress conditions in plants such as temperature, salinity, dehydration, low oxygen, mechanical damage etc; and play key roles in plant development during C:N balance42. Another study shows role of flavanone 3-hydroxylase in biosynthesis of phenolics (taxifolin and catechin) in spruce to confer defense against bark beetle and fungus associates43. Role of flavonol-specific genes like flavonol synthase/flavanone 3-hydroxylase has also been investigated in phenylpropanoid biosynthesis pathways in cucumber, which shows higher anthocyanin content under N deficiency10. The ACC oxidase was over-expressed in stolons under N stress and its role in ethylene biosysnthesis is well known under N limitation in cucumber10. The roles of these genes strengthen their importance in potato under N stress.

Some more genes like miraculin was found in all tissues, whereas patatin-05 and class III peroxidise were observed in root and stolon, and class I chitinase gene was differentially expressed in stolon only. All these genes were highly up-regulated in stolons under N stress. Patatin is the major protein of potato tubers. Reduced activities of defence-related enzymes like chitinase, chitosanase and peroxidase were reported in Arabidopsis thaliana at low N44. In stolons, dehydration-responsive protein RD22, hydroxyproline-rich glycoprotein and sesquiterpene synthase were highly down-regulated under N stress in potato. Recently, dehydration-responsive protein RD22 has been demonstrated for role under salt stress in soybean seedlings45. Hydroxyproline-rich glycoproteins function in plant cell wall in response to pathogens attack, as its role has been investigated in resistance to downy mildew in pearl millet46. Sesquiterpene synthase has been reported in many plants to play key roles in terpeniod metabolites synthesis and provide adaptation to adverse conditions under biotic and abiotic stresses in Santalum album47. Thus, these genes play vital roles in adaptation to N deficiency in potato.

Several other genes like protein kinases and phospholipids showed differential expression under N deficiency in potato. These genes were either over- or under-expressed in the tissues specific response such as tartrate-resistant acid phosphatases (shoot); phospholipase C (root); glycerophosphodiester phosphodiesterases, thylakoid soluble phosphoprotein, and histone deacetylase (shoot and root); cysteine protease inhibitor 1, metallocarboxypeptidase inhibitor, and proteinase inhibitor type-2 (root and stolon); and methylketone synthase Ib (root, shoot and stolon). The tartrate-resistant acid phosphatases (TRAcPs) and glycerophosphodiester phosphodiesterases were highly down-regulated under N stress in shoots. The TRAcPs, known as purple acid phosphatases (PAPs), have been identified in plants, animals and fungi, and particularly in Arabidopsis thaliana under phosphate starvation48. Glycerophosphodiester phosphodiesterase plays an important role in lipid metabolism, which releases inorganic phosphorous from phospholipids during P starvation49,50. One GDSL esterase/lipase gene was down-regulated genes under N stress in potato roots. A Gossypium hirsutum GDSL esterase/lipase functions as hydrolytic enzyme and involved in ovule and fibre development and plays a key role in seed development in Arabidopsis thaliana51. The phospholipase C gene was up-regulated under N stress in potato roots, and a previous study has characterized isoforms of phosphoinositide-specific phospholipase C in potato52. The methylketone synthase Ib and thylakoid soluble phosphoprotein genes were down-regulated in shoots and the protein phosphatase-2c was highly up-regulated in roots under N stress. The plant protein serine/threonine phosphatases (PP1/PP2A and PP2C) are the major families found in plants and animals, which function in kinase-associated protein phosphatase signalling pathways. The role of ABI1/ABI2, PP2C enzymes, has been envisaged in signal transduction of abscisic acid pathway in Arabidopsis thaliana; and alfalfa PP2C family functions as a negative regulator in plant mitogen-activated protein kinase pathways53. The inhibitors like cysteine protease inhibitor 1, metallocarboxypeptidase inhibitor and proteinase inhibitor type-2 P303.51 were up-regulated under N stress in stolons. Role of cysteine proteinase inhibitor has been envisaged in plants for growth and developmental, defense and stress response54. A complex network of proteases and protease inhibitors has been observed during multiple functions related to time and space in various biological pathways in plant life cycle55. Metallocarboxypeptidase inhibitors have been found in cDNA libraries of the wounded tissues of leaves treated with abscisic acid56. Histone deacetylase play key roles in regulation of gene expression of histone deacetylation in response to biotic and abiotic stresses in plants57. Thus, our study unveils many gene networks playing very crucial roles in potato during N stress tolerance metabolism.

Conclusions

A large number of potential genes were discovered in this study in potato in response to N deficiency. Some of these genes showed very high up-regulation or down-regulation under N stress such as glutaredoxins, TFs (Myb, WRKY, BTB/POZ domain, LOB domain and F-box), high-affinity nitrate transporter/nitrate transporter, sodium/proline symporter, glucose-6-phosphate/phosphate translocator 2, cysteine protease inhibitor 1, metallocarboxypeptidase inhibitor, dehydraton responsive protein RD22, hydroxyproline-rich glycoprotein, ACC oxidase, and sesquiterpene synthase so on. In future, more investigations would be required in field-grown plants to observe transcripts changes in response to contrasting N regimes. Furthermore, functional validation of these potential genes via transgenic manipulation would also be essential. Taken together, the present study has greatly improved our knowledge on enrichment of gene networks and regulatory elements involved in N metabolism pathways in potato, which would necessarily strengthen research on N metabolism in potato in future.

Methods

Plant materials and nitrogen treatments

A potato variety Kufri Jyoti (KJ), N-inefficient2,15, was used in this study. In vitro plants were maintained at Indian Council of Agricultural Research-Central Potato Research Institute, Shimla (31.1048°N, 77.1734°E 2,276 m above mean sea level), Himachal Pradesh, India. Plants were grown in three replicates in aeroponic culture with two N treatments (low N: 0.2 mM N; and high N: 4 mM N, control) and other micro- and micro-nutrients remain same as described by Tiwari et al.15. N was supplied as NO3− and NH4+ forms in both treatments. Plant tissues (shoots, roots and stolons) from both low N and high N were collected in three replicates after 30 days (roots, and shoots i.e. leaves) and 40 days (stolons, when stolon initiation/ tuberization were started) of planting in aeroponic culture. The collected six tissues (KJ-HighN Shoot, KJ-LowN Shoot, KJ-HighN Root, KJ-LowN Root, KJ-HighN Stolon, and KJ-LowN Stolon) were snap-frozen in liquid nitrogen and stored at −80 °C for transcriptomes sequencing. The tissue from three biological replicates was pooled for single sample for transcriptomes sequencing.

Total RNA isolation and Illumina NextSeq500 PE library preparation

Total RNA was isolated from plant samples (shoots, roots and stolons of both low N and high N) using modified c-TAB and Lithium Chloride method58. The qualities and quantities of the isolated RNA samples were checked on 1% denaturing RNA agarose gel and NanoDrop (ThermoFisher Scientific, Wilmington, Delaware USA), respectively. RNA-seq paired end sequencing libraries were prepared from the QC passed RNA samples using Illumina TruSeq Stranded mRNA sample prep kit following the manufacturer instructions (Illumina, San Diego, CA, USA). Briefly, mRNA was enriched from the total RNA using poly-T attached magnetic beads, followed by enzymatic fragmentation, 1st strand cDNA conversion using SuperScript II and Act-D mix to facilitate RNA dependent synthesis. The 1st strand cDNA was then synthesized to second strand using second strand mix. The dscDNA was then purified using AMPure XP beads followed by A-tailing, adapter ligation and then enriched by limited no of PCR cycles. The PCR enriched libraries were analyzed on 4200 Tape Station system (Agilent Technologies, Santa Clara, CA, USA) using high sensitivity D1000 Screen tape as per the manufacturer instructions. The brief bioinformatics work flow is depicted in Fig. S1.

Cluster generation and sequencing

After obtaining the Qubit 3.0 (ThermoFisher Scientific, Waltham, Massachusetts, USA) concentrations for the libraries and the mean peak sizes from Agilent Tape Station profiles, the PE illumina libraries were loaded onto NextSeq500 for cluster generation and sequencing. Paired-end sequencing allows the template fragments to be sequenced in both the forward and reverse directions on NextSeq500. The kit reagents were used in binding of samples to complementary adapter oligos on paired-end flow cell. The adapters were designed to allow selective cleavage of the forward strands after re-synthesis of the reverse strand during sequencing. The copied reverse strands were then used to sequence from the opposite end of the fragment.

RNA sequencing and high quality read statistics

The raw data were processed to obtain high quality clean reads using Trimmomatic v0.38 to remove adapter sequences, ambiguous reads (reads with unknown nucleotides “N” larger than 5%), and low-quality sequences (reads with more than 10% quality threshold (QV) < 20 phred score). A minimum length of 50 nt (nucleotide) after trimming was applied. After removing the adapter and low quality sequences from the raw data, high quality reads were obtained. These high quality (QV > 20), paired-end reads were used for reference-based reads mapping. Parameters considered for filtration were as follows: i) SLIDINGWINDOW: Sliding window trimming of 10 bp, cutting once the average quality within the window falls below a threshold of 20, ii) LEADING: Cut bases off the start of a read, if below a threshold quality of 20, and iii) TRAILING: Cut bases off the end of a read, if below a threshold quality of 20.

Reads mapping to the reference potato genome

The reference genome of Solanum tuberosum Group Phureja DM1-3, with a genome size of 840 Mb and the associated annotations were downloaded from the Spud DB database of the potato genome sequence16. The download links for genome was http://solanaceae.plantbiology.msu.edu/data/potato_dm_v404_all_pm_un.fasta.zip and for annotation was http://solanaceae.plantbiology.msu.edu/data/PGSC_DM_V403_genes.gff.zip. The high quality reads were mapped to the reference genome using TopHat v2.1.1 with default parameters59.

Differential gene expression analysis

Cufflinks v2.2.1 program assembles transcriptomes from RNA-seq data and quantifies their expression. The individual gtf files of the transcriptomes were used for differential gene expression analysis using cuffdiff by blind dispersion method60. There are total of 39,028 protein coding genes present in annotation file of Solanum tuberosum Group Phureja DM1-3. Differential gene expression analysis was performed using cuffdiff version 2.2.1 between low N and high N (control) for shoots, roots and stolons. FPKM values were used to calculate the log fold change as log2 (FPKM_experimental/FPKM_control). Log2 FC values greater than zero were considered up-regulated whereas less than zero were down-regulated along with P-value threshold of 0.05 for statistically significant results. Venn diagrams of up-regulated and down-regulated DEGs were prepared using Venny 2.1 tool to analyze tissues-specific and tissues-independent genes for all three tissues61.

Heat map

An average linkage hierarchical cluster analysis was performed with the top 50 DEGs using multiple experiments viewer (MeV v4.9.0)62. The heatmap shows level of gene abundance. Levels of expression are represented as log2 ratio of gene abundance between low N and high N (control) samples. Heatmaps were constructed using the log-transformed and normalized value of genes based on Pearson uncentered distance and average linkage method. In heatmap, each horizontal line refers to a gene. The color represents the logarithmic intensity of the expressed genes. Relatively higher gene expression values are shown in red colour and lower expression are shown in green colour.

Scatter plot

The Eurofins proprietary R scripts were used to depict graphically the representing expression of genes in two distinct conditions of each sample combination i.e. low N and high N (control). It helps to identify DEGs in one sample with respect to another and also allows comparing two values associated with genes. In scatter plot, each dot represents a gene. The vertical position of each gene represents its expression level in the control (high N) samples while the horizontal position represents its expression level in low N samples. Thus, genes that fall above the diagonal are over-expressed and genes that fall below the diagonal are under-expressed as compared to their median expression level in experimental grouping of the experiment.

Volcano plot

The Eurofins proprietary R scripts were used to depict the graphical representation and distribution of DEGs which were found in low N as well as high N samples. The volcano plot arranges DEGs along dimensions of biological as well as statistical significance. The red block on the right side of zero represents the up regulated genes whereas green block on the left side of zero represents significant down regulated genes. The Y-axis represents the negative log of p-value (p < 0.05) of the performed statistical test where data points with low p-values (highly significant) appearing towards the top of the plot. Grey block shows the non-DEGs.

GO analysis

The GO annotations of the DEGs were obtained from the Ensembl Plants database for Solanum tuberosum. The GO annotations were categorized into up-regulated, down-regulated, expressed in both and exclusive DEGs. The information on number of genes was assigned into three main GO domains (biological process, cellular component, and molecular function). The bar plots, depicting the GO distribution, were obtained through the WEGO portal (http://wego.genomics.org.cn/cgi-bin/wego/index.pl)63.

KEGG pathways analysis

The functional annotations of the DEGs were carried out against the curated KEGG GENES database using KAAS (KEGG Automatic Annotation Server-. (http://www.genome.jp/kegg/ko.html)64. The KEGG orthology (KO) database of Nightshade family was used as the reference for pathways mapping. The result contains KO assignments and automatically generated KEGG pathways using KAAS BBH (bidirectional best hit) method against available database.

Conserved motifs analysis

Selected highly up-regulated ( ≥ 3 log2 fold change) and down-regulated ( ≤ −3 log2 fold change) DEGs with known function in shoots, roots and stolons were analysed using MEME (version 5.1.0) software65.

Real time – quantitative polymerase chain reaction (RT-qPCR) analysis

Selected twelve differentially expressed genes (two genes from each shoots, roots and stolons from both low N and high N) were validated by RT-qPCR analysis. The coding sequences of the selected genes were downloaded from the potato genome sequence database (http://solanaceae.plantbiology.msu.edu/pgsc_download.shtml), and RT-qPCR primers were designed using IDT PrimerQuest Tool (https://eu.idtdna.com/Primerquest/Home/Index) following default parameters, as summarised in Table S6. The same samples were used RNA isolation using RNeasy Mini Kit (Qiagen, Venlo, Limburg, Netherlands) and cDNA synthesis using TaqMan Reverse Transcription Reagent (Applied Biosystems, New Jersey, USA). RT-PCR analysis was performed using Power SYBR Green PCR Master Mix in ABI PRISM HT7900 (Applied Biosystems Warrington, UK) following temperature/time profile 50 °C/2 min; 95 °C/10 min; and 40 cycles of 95 °C/15 s, 60 °C/1 min, and 72 °C/30 s with an internal standard potato ubiquitin-ribosomal protein gene (ubi3; L22576) as described in Tiwari et al.66.

Plant biomass, total chlorophyll and N content analysis

Phenotypic traits, plant biomass, total chlorophyll and total N contents were measured as procedures described by Tiwari et al.15. Plants were grown for full crop cycle till 90 days crop duration to measure these traits. Observations were recorded on per pant basis in three replications for plant height (cm), leaf length (cm), leaf width (cm), tuber yield (g/plant), tuber number, total N content, total chlorophyll content (mg/g fresh weight), root dry weight (g) and shoot dry weight (g). Data were analysed by one way analysis of variance (ANOVA) using XLSTAT 2018.5 at least level of significance (p ≤ 0.05).