Abstract

   Nitrogen (N) is critical to the growth and productivity of crops. To
   understand the molecular mechanisms influenced by N stress, we used
   RNA-Sequencing (RNA-Seq) to analyze differentially expressed genes
   (DEGs) in root and leaf tissues of spinach. N stress negatively
   influenced photosynthesis, biomass accumulation, amino acid profiles,
   and partitioning of N across tissues. RNA-seq analysis revealed that N
   stress caused most transcriptomic changes in roots, identifying 1,346
   DEGs. High-affinity nitrate transporters (NRT2.1, NRT2.5) and glutamine
   amidotransferase (GAT1) genes were strongly induced in roots in
   response to N deplete and replete conditions, respectively. GO and KEGG
   analyses revealed that the functions associated with metabolic pathways
   and nutrient reservoir activity were enriched due to N stress. Whereas
   KEGG pathway enrichment analysis indicated the upregulation of DEGs
   associated with DNA replication, pyrimidine, and purine metabolism in
   the presence of high N in leaf tissue. A subset of transcription
   factors comprising bHLH, MYB, WRKY, and AP2/ERF family members was
   over-represented in both tissues in response to N perturbation.
   Interesting DEGs associated with N uptake, amino acid metabolism,
   hormonal pathway, carbon metabolism, along with transcription factors,
   were highlighted. The results provide valuable information about the
   underlying molecular processes in response to N stress in spinach and;
   could serve as a resource for functional analysis of candidate
   genes/pathways and enhancement of nitrogen use efficiency.

Introduction

   Nitrogen (N), being a constituent of amino acids, nucleic acids,
   cofactors, and secondary metabolites, is central to most biological
   activities in plants. Effects of N on plant morphology rely on its
   concentration, which is a determining factor in several metabolic
   processes, resource allocation, growth, and development of plants
   [[29]1]. N moves along a complex branching and merging pathway which
   interacts at numerous sites with the carbon (C) flow, ion, and
   assimilates at the cell and whole plant level [[30]2]. N is typically
   transported in the form of nitrate (NO[3]^−), dissolved ammonia, and
   amino acids. Most of the NO[3]^− is reduced in the shoots, while
   ammonia incorporates into amino acids for long-distance translocation.
   Among various external factors that affect plant growth, N
   partitioning, in particular, defines the productivity by altering the
   ratios of roots to shoots [[31]3, [32]4]. For example, in spinach
   (Spinacia oleracea), suboptimal N nutrition enhanced root growth at the
   expense of shoot biomass [[33]5, [34]6].

   Spinach is one of the popular nutritious green leafy vegetables grown
   in most parts of the world. Although understanding the mechanisms of N
   metabolism in spinach is of considerable significance, most of the
   previous studies have mainly focused on physiological and agronomical
   aspects. Most plants can consume only half of the applied N, losing it
   in the form of NO[3]^−, which leads to groundwater contamination,
   health, and environmental hazards [[35]7, [36]8]. In commercial spinach
   production, it is estimated that about 60% of N is lost through
   leaching [[37]9], most likely due to its shallow root system [[38]10]
   and short production cycle. Although spinach needs significant amounts
   of N to sustain rapid growth, it is relatively weak in NO[3]^−
   reduction [[39]11–[40]13]. One way to tackle this crisis without
   compromising productivity is to improve nitrogen use efficiency (NUE)
   [[41]14].

   N perturbation induces series of molecular changes associated with
   intrinsic processes such as chlorophyll synthesis [[42]15], root
   architecture [[43]16], and metabolism of sugars and sugar phosphates
   [[44]17]. In particular, N starvation results in diminished levels of N
   containing metabolites such as amino acids glutamate (Glu) and
   glutamine (Gln) [[45]18]. As plants take up NO[3]^− and ammonium
   (NH[4]^+) from the soil and assimilate into amino acids in
   roots/shoots, improving the amino acid partitioning from source to sink
   would be an effective strategy to optimize N utilization [[46]19]. To
   date, quite a few studies have been conducted to investigate
   genome-wide expression analysis in response to N stress in various
   plants such as Arabidopsis [[47]20–[48]22], rice [[49]23–[50]25], maize
   [[51]26, [52]27], wheat [[53]28, [54]29], and Brassica species
   [[55]30]. Although this information should help in understanding
   general molecular processes, each species has its unique and
   non-overlapping mechanisms to withstand N stress. Besides few studies
   focusing on mechanisms of N uptake [[56]31, [57]32] and genotypic
   variation [[58]6, [59]33, [60]34], not much is known about the
   transcriptomic changes in the vegetative or non-vegetative tissues of
   spinach in response to N perturbations.

   Here we performed RNA-sequencing (RNA-Seq) of leaf and root tissues of
   spinach combined with physiological and metabolic analyses using two N
   treatments: high N (200 ppm) and low N (50 ppm). A total of 1346 and
   1136 differentially expressed genes (DEGs) were identified in response
   to N in root and leaf tissues, respectively. The changes in the
   expression of the identified genes and the impacted pathways in
   response to N would help in advancing our understanding of
   tissue-specific molecular and physiological processes in spinach. The
   novel genes/pathways identified in this study have provided new targets
   for genetic manipulation or introgression breeding to improve NUE in
   spinach.

Materials and methods

Plant material and growth conditions

   The spinach plants were grown in an environmentally controlled growth
   chamber at the Texas A&M AgriLife Research and Extension Center,
   Uvalde, Texas. The seeds of spinach variety ‘Banjo’ were grown in
   plastic containers containing Turface ® based growing medium under 200
   μmol·mPPPPP^-2PPPPP·s PPPPP^-1 PPPPPlight intensity PPPPP PPPPP (12
   hours each light and dark period) at 23°C and 60–70% relative humidity.
   Two treatments, high N and low N, where N was added in the form of
   Ca(NO[3])[2], were used to maintain the concentration of 200 ppm for
   high N and 50 ppm for low N. CaCl[2] was compensated to keep the same
   concentration of calcium in both treatments. The concentrations of
   remaining macro- and micro-elements were kept constant using N-free
   Hoagland’s nutrient solution (No. 2 Basal Salt Mixture, Caisson Labs,
   Smithfield, UT, USA) comprising 2.86 mg·LPPPPP^−1PPPPP H[3]BO[3],
   554.90 mg·LPPPPP^−1PPPPP CaCl[2]R, 0.045 mg·LPPPPP^−1PPPPP CuCl[2],
   33.0 mg·LPPPPP^−1PPPPP C[10]RRH[12]N[2] NaFeO[8]·3H[2]O, 240.33
   mg·LPPPPP^−1PPPPP MgSO[4], 1.81 mg·LPPPPP^−1PPPPP MnCl[2]R·4H[2]O,
   0.025 mg·LPPPPP^−1PPPPP Na[2]MoO[4]·2H[2]O, 372.70 mg·LPPPPP^−1PPPPP
   KCl, 136.025 mg·LPPPPP^−1PPPPP KH[2]PO[4], and 0.1 mg·LPPPPP^−1PPPPP
   ZnCl[2]R. Samples collected from 6-week old plants were frozen in
   liquid nitrogen and stored at −80°C until subsequent analyses. Three
   independent plants were used for metabolic analysis, total RNA
   extractions, physiological and biochemical analysis.

Determination of physiological traits, N, and free amino acids

   Leaf photosynthetic rate (Pn), transpiration rate (Tr), stomatal
   conductance (Cs), and intercellular CO[2] concentration were measured
   with LICOR-6400-XT Photosynthetic system (Lincoln, USA) from the fully
   expanded leaves of 6-week old plants. The chlorophyll content was
   measured using a portable chlorophyll meter (SPAD-502, Konica Minolta,
   Tokyo, Japan). The plant samples were analyzed for total N (TKN),
   NO[3]^−, NH[4]^+ Rusing Easy Chem Plus analyzer (Chinchilla Scientific,
   Oak Brook, IL). The amino acid analysis was performed using WatersPPPPP
   ^RPPPPP Acquity H-class UPLC system coupled to WaterPPPPP^RPPPPP’s Xevo
   TQs MS-MS mass spectrometer with electrospray ionization (ESI) probe
   following pre-established protocol [[61]35]. Data integration and
   quantitation were performed using the Water’s TargetLynx^™ software.
   Statistical analysis was performed using JMP Pro 14 (SAS Institute,
   Cary, NC, USA).

RNA extraction, cDNA library preparation, and sequencing

Sample collection and preparation

   The flash-frozen plant samples in liquid nitrogen were ground to a fine
   powder using a paint shaker (Harbil, Wheeling, IL, USA) and
   3-mm-diameter steel balls (Abbott Ball, West Hartford, CT, USA). Total
   RNA was extracted using an RNeasy® Plant Mini Kit (QIAGEN Sciences,
   Germantown, MD, USA) as per the manufacturer’s protocol and treated
   with DNase1 (QIAGEN Sciences, Germantown, MD, USA). The purity of the
   RNA was analyzed using the NanoPhotometer® spectrophotometer (IMPLEN,
   CA, USA). RNA integrity and quantitation were assessed using the RNA
   Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent
   Technologies, CA, USA).

Library preparation for transcriptome sequencing

   One μg total RNA per sample was used for the RNA library preparations.
   Sequencing libraries were generated using NEBNext® Ultra^™ RNA Library
   Prep Kit for Illumina® (NEB, USA) following the manufacturer’s
   recommendations, and index codes were added to attribute sequences to
   each sample. Library concentration was first quantified using a Qubit
   2.0 fluorometer (Life Technologies), diluted to 1 ng/μl before checking
   the insert size on an Agilent Bioanalyzer 2100 system and quantified to
   greater accuracy by quantitative PCR (Q-PCR) (library activity >2 nM).

Data processing and analysis

   Clustering, sequencing, and quality control. The clustering of the
   index-coded samples was performed on a cBot Cluster Generation System
   using PE Cluster Kit cBot-HS (Illumina) according to the manufacturer’s
   instructions. After cluster generation, the libraries were sequenced on
   an Illumina Hiseq platform, and 150 bp paired-end reads were generated.
   Raw reads of fastq format were processed to obtain clean reads by
   removing the adapter, reads containing ploy N, and low-quality reads
   from raw data. At the same time, Q20, Q30, and GC content, the clean
   data were calculated.

   Reads mapping to the reference genome. Reference genome and gene model
   annotation files were downloaded from SpinachBase
   ([62]http://spinachbase.org/). Index of the reference genome was built
   using Bowtie v2.2.3, and paired-end clean reads were aligned to the
   reference genome using TopHat v2.0.12.

   Gene expression quantification and DEG analysis:

   HTSeq v0.6.1 was used to count the reads mapped to each gene. FPKM
   [[63]36] of each gene was calculated based on the length of the gene
   and reads count mapped to this gene. Differential expression analysis
   of high N and low N conditions (three biological replicates per tissue
   per treatment) was performed using the DESeq R package (1.18.0)
   [[64]37]. Genes with P-value < 0.05 found by DESeq were assigned as
   differentially expressed.

   GO and KEGG enrichment analysis. Gene Ontology (GO) [[65]38] enrichment
   analysis of differentially expressed genes was implemented by the GOseq
   R package, in which gene length bias was corrected. GO terms with
   P-value less than 0.05 were considered significantly enriched by DEGs.
   KOBAS software was used to test the statistical enrichment of
   differential expression genes in the Kyoto encyclopedia of genes and
   genome pathways database (KEGG; [66]http://www.genome.jp/kegg/)
   [[67]39].

Validation by quantitative real-time PCR

   To validate the RNA-seq data, quantitative real-time PCR (RT-qPCR) was
   conducted to examine the expression pattern of twenty DEGs. Primer
   Premier 3.0 software was used to design gene-specific primers on the
   basis of the selected unigene sequences ([68]S1 Table). Total RNA was
   extracted with the Quick-RNA^™ Miniprep Kit (Zymo Research Corporation,
   Irvine, CA) treated with DNase1 (Zymo Research Corporation, Irvine,
   CA), and subjected to reverse transcription using iScript RT Supermix
   (Bio-Rad Laboratories, Inc, Hercules, USA). The quality and quantity of
   the RNA were analyzed by a Denovix DS-11+ spectrophotometer
   (Wilmington, Delaware, USA). Gene expression analysis via reverse
   transcription-qPCR was performed using the BioRad CFX96 qPCR instrument
   using SsoAdv Univer SYBR GRN Master Kit (Bio-Rad Laboratories, Inc,
   Hercules, USA). The expression levels of selected DEGs were normalized
   by comparing with an internal reference gene, 18SrRNA [[69]40]. The
   relative expression levels (Cq values) for each gene were normalized to
   that of reference genes by taking an average of three biological
   replicates. The relative expression levels were calculated using the
   2^−ΔΔCt method [[70]41]. All RT-qPCR were repeated in three biological
   and three technical replications.

Results and discussions

Validation of N stress for expression profiling in spinach

   N availability affects its acquisition and metabolism in plants. Hence
   a detailed phenotypic characterization preceded sampling of material
   for RNA-Seq analysis. A generalized response to N limitation was
   confirmed by evaluating responses of gas exchange parameters in the
   leaves of 6-week old spinach plants. N stress significantly reduced the
   net photosynthetic rate, stomatal conductance, and transpiration rate,
   implying compromised photosynthetic performance ([71]Table 1). Further,
   elemental N and NO[3]^− [was] significantly reduced in leaf, petiole,
   and root tissues validating treatment effect ([72]Fig 1). NH[4]^+
   content was also lower in leaf tissues due to N stress. As most N in
   plants is converted to organic form as free amino acids, we compared
   the amounts of six highly abundant amino acids; Asparagine (Asn),
   Aspartate (Asp), Gln, Glu, Serine (Ser) and Alanine (Ala) in four
   tissue types of spinach. The contents of all the amino acids showed a
   significant reduction in the roots, while foremost N carrying amino
   acids (Asp, Asn, Glu, Gln) were significantly decreased in the leaf
   tissue ([73]Fig 2). Unlike total N and content of NO[3]^−, no
   significant changes were observed in amino acids in petiole due to N
   stress.

Table 1. Responses of gas exchange parameters to N stress in spinach.

   Treatment P[n] (μmol CO[2] m^-2s^-1) G[s] (mol[2] H[2]O m^-2 s^-1) C[i]
   (μmol of CO[2] / mol) E (mmol H[2]O m^-2 s^-1) SPAD
   high N 9.79 ± 0.42 0.58 ± 0.04 347.1 ± 2.12 5.56 ± 0.27 43.7 ± 2.19
   low N 6.01 ± 0.46 0.26 ± 0.05 341.2 ± 5.13 3.79 ± 0.5 35.8 ± 1.15
   p-Value 0.00[74]^* 0.01[75]^* 0.34 0.02[76]^* 0.03[77]^*
   [78]Open in a new tab

   P[n]—net photosynthetic rate, E- transpiration rate, G[s]—stomatal
   conductance, C[i]—intercellular CO[2] concentration, SPAD (Soil Plant
   Analysis Development) units * significant at P<0.05.

Fig 1. Percent changes in N content in spinach tissues due to N stress.

   [79]Fig 1
   [80]Open in a new tab

   (n = 3, Mean ± SE) * significant at P<0.05.

Fig 2. Changes in free amino acids in spinach due to N stress.

   [81]Fig 2
   [82]Open in a new tab

   Asparagine (Asn), Aspartate (Asp), Glutamine (Gln), Glutamate (Glu),
   Serine (Ser), and Alanine (Ala) (n = 3, Mean ± SE) * significant at
   P<0.05.

Transcriptome sequencing and assembly of sequencing data

   The transcriptomic changes induced by N stress in a spinach leaf and
   root tissues were analyzed by RNA-Seq. Total twelve independent
   libraries; six from leaves (HNL1, HNL2, HNL3 for high N and LNL1, LNL2,
   LNL3 for low N) and six from roots (HNR1, HNR2, HNR3 for high N and
   LNR1, LNR2, LNR3 for low N) were sequenced using the Illumina HiSeq
   platform. On average, 45.53 to 43.47 million raw reads were generated
   for leaf and root tissues in both N treatments ([83]Table 2). Across
   all reads, the Q20 and Q30 percentage was more than 97 and 93%,
   respectively (sequencing error rate was less than 0.03%), and GC
   content for the libraries was more than 43%. Additional significant
   characteristics of these libraries are summarized in [84]S2 Table.
   Among all the libraries, the ratio of total mapped reads and multiple
   mapped reads in both tissue types was above 85% and 2.5%, respectively.
   Average 85% reads were uniquely mapped to the reference genome in both
   tissue types under both N treatments. The data generated from all
   libraries provided a foundation for quality analyses. The RNA-Seq
   dataset is accessible through GEO Series accession number [85]GSE145943
   ([86]https://www.ncbi.nlm.nih.gov/geo/).

Table 2. Summary of sequencing data quality of spinach samples in both tissue
types and N treatments.

   Treatment Tissue Sample name Raw reads Clean reads Error rate (%) Q20
   (%) Q30 (%) GC content (%)
   high N Leaf HNL1 41597938 40201618 0.03 97.76 93.65 44.15
   HNL2 47471760 45594984 0.03 97.72 93.64 44.48
   HNL3 47246618 45862032 0.03 97.53 93.07 44.20
   Root HNR1 43542944 42695366 0.03 97.66 93.34 43.31
   HNR2 42288790 41132378 0.03 97.52 93.08 43.21
   HNR3 42159430 39940286 0.03 97.60 93.30 43.22
   low N Leaf LNL1 46979108 42951112 0.03 97.75 93.62 44.14
   LNL2 43156630 41716986 0.03 97.85 93.81 44.62
   LNL3 46758928 45391850 0.03 97.52 93.10 44.15
   Root LNR1 45260090 44242952 0.03 97.71 93.49 43.29
   LNR2 46795700 45438476 0.03 97.72 93.49 43.57
   LNR3 40777766 39803014 0.03 97.90 93.95 43.33
   [87]Open in a new tab

   The correlations among biological replicates were assessed using the
   Pearson correlation coefficient ([88]S1 Fig) to validate the
   reliability of RNA-seq data. The libraries for the same treatment
   (i.e., biological replicates) were highly correlated. The weak
   correlation between tissue types within treatments suggests a
   significant impact of N stress on gene expression profiles. Distinct
   gene expression levels under different experiment conditions also
   suggest that root tissue is more sensitive to N perturbation than leaf
   tissue ([89]S2 Fig).

Analysis of differentially expressed genes (DEGs)

   The cluster analysis confirmed that a large number of genes were
   differentially expressed in leaf and root tissues in both treatments
   (low N and high N) ([90]S3 Fig). The DEG analysis of roots revealed
   that 1346 transcripts were significantly (p < 0.05) altered when plants
   were grown in high N. These included 726 upregulated and 620
   downregulated transcripts ([91]Fig 3, [92]S3 Table, [93]S4 Table). A
   total of 1136 genes were differentially expressed in leaf tissue, of
   which 550 genes were upregulated ([94]Fig 3, [95]S5 Table), while 586
   were downregulated in high N ([96]Fig 3, [97]S6 Table). Within
   treatments, the number of DEGs was higher in high N treatment while
   within tissue types, it was more abundant in root tissue in both N
   treatments ([98]Fig 4). Between treatments, the DEG analysis revealed
   that expression of 222 genes was significantly (p < 0.05) affected by N
   stress, of which 88 genes were upregulated and 134 genes downregulated
   in the presence of high N availability ([99]S4 Fig). Together, all the
   DEGs in root and leaf tissues under both treatments are presented as a
   single Venn diagram ([100]Fig 5). Total 200 and 181 DEGs were uniquely
   expressed in root and leaf tissue, respectively.

Fig 3.

   [101]Fig 3
   [102]Open in a new tab

   The volcano maps showing the number of differentially expressed genes
   (DEGs) in the root (A) and leaf (B) tissues in high N. Red dots
   represent up-regulated genes, and green dots represent down-regulated
   genes (p < 0.05).

Fig 4.

   [103]Fig 4
   [104]Open in a new tab

   The volcano maps showing the number of differentially expressed genes
   (DEGs) in the leaf of spinach plants grown under low N (A) and high N
   (B) treatments. Red dots represent up-regulated genes, and green dots
   represent down-regulated genes (p < 0.05).

Fig 5. The Venn diagram presenting the number of differentially expressed
genes (DEGs) within the treatment.

   [105]Fig 5
   [106]Open in a new tab

   The sum of the number in the circle presents the total number of DEGs,
   and the overlap presents the genes in common.

Functional annotation and GO enrichment analysis of DEGs

   Among the DEGs that were significantly upregulated in leaf under Low N,
   functional annotation showed that 15 biological processes, 9 cell
   component metabolic pathways, and 6 molecular functions were
   significantly over-represented (p-value < 0.05) ([107]Fig 6).
   Similarly, DEGs showing down-regulation in leaf affected 21 molecular
   functions, followed by 8 biological processes under low N availability
   ([108]Fig 6). In leaves, the low N availability mainly upregulated DEGs
   associated with biological processes, metabolic and organo-nitrogen
   metabolic processes, while the membrane and its integral components,
   and protein modification processes were over-represented in DEGs that
   were down-regulated. In the case of high N treatment, among the
   up-regulated DEGs, functional annotation showed 15 biological
   processes, 12 cellular components, and 3 molecular functions were
   significantly affected ([109]Fig 7). While among the downregulated
   DEGs, 7 biological processes and cellular components, and 16 molecular
   functions were over-represented. Similar to the upregulated DEGs
   associated with biosynthetic processes in leaves under low N, high N
   treatment clearly represented processes related to the biosynthesis and
   metabolism of nitrogenous compounds and peptides. Among the
   downregulated DEGs, functions associated with membrane modifications,
   response to a stimulus, protein modifications, and transferase activity
   were significantly over-represented in leaf tissue in high N treatment.

Fig 6.

   [110]Fig 6
   [111]Open in a new tab

   Gene ontology (GO) enrichment bar chart of DEGs showing most enriched
   G) terms of (A) upregulated and (B) downregulated DEGs in leaf under
   Low-N. The y-axis is the enriched GO term; the x-axis is the number of
   DEGs enriched in the listed term. Colors represent different GO types:
   biological process, cellular component, and molecular function. *
   significantly enriched term (p-value < 0.05).

Fig 7.

   [112]Fig 7
   [113]Open in a new tab

   Gene ontology (GO) enrichment bar chart of DEGs showing most enriched
   G) terms of (A) upregulated and (B) downregulated DEGs in leaf tissues
   under High-N. The y-axis is the enriched GO term, and the x-axis is the
   number of DEGs enriched in the listed term. Colors represent different
   GO types: biological process, cellular component, and molecular
   function. * significantly enriched term (p-value < 0.05).

   Overall, in root tissue, the GO terms such as membrane and cellular
   components, biological processes associated with a single organism,
   response to stimulus and transferase activity were over-represented in
   the up-regulated DEGs while nutrient reservoir activity related
   function was represented in the down-regulated DEGs due to N stress
   ([114]S5 Fig). While in the case of leaf tissue, only the functions
   associated with catalytic, transferase, and kinase activities were
   represented in the up-regulated DEGs during N stress ([115]S6 Fig).
   Taken together, N stress over-represented GO terms such as biological
   processes, molecular functions like catalytic and transferase activity,
   and cellular locations like a membrane or its components in the
   up-regulated DEGs ([116]S7 Fig).

KEGG pathways associated with differentially expressed genes

   The enrichment of KEGG pathways was similar in leaf and root tissues in
   both N treatments ([117]Fig 8), showing over-representation of DEGs
   associated with metabolic pathways, especially in amino acid metabolism
   and biosynthesis of secondary metabolites. DEGs related to amino acid
   metabolic pathways such as Ala, Asp, and Glu synthesis, glutathione
   metabolism, phenylalanine metabolism were enriched in both N
   treatments. KEGG pathways of DEGs identified in both N treatments in
   leaf Vs. root tishsue is shown in [118]S8 Fig. Like tissue types, the
   analysis showed that the number of DEGs in both N treatments were also
   associated with metabolic and biosynthesis of secondary metabolites
   pathways.

Fig 8.

   [119]Fig 8
   [120]Open in a new tab

   KEGG pathway enrichment scatters plot analysis of DEGs in leaf (A) and
   root (B) tissues between the N treatments. The Rich factor is the ratio
   of DEGs in this pathway term to all gene numbers annotated in this
   pathway term. A q value is the corrected p-value ranging from 0 to 1,
   and a lower value indicates higher pathway enrichment. The pathway
   names are shown on the vertical axis, rich factor on the horizontal
   axis, the size of the point represents the number of DEGs, and the
   color of the dot represents the q value.

Quantitative real-time-PCR validation of DEGs from RNA-Seq

   To validate RNA-seq results, we selected twenty DEGs to confirm their
   expression patterns of the by quantitative real-time PCR. The
   fold-changes of the selected genes using RT-qPCR were consistent
   ([121]S9 Fig) with the results obtained with RNA-Seq analysis (R^2 =
   0.085, r^2 = 0.93 for roots and R^2 = 0.080, r^2 = 0.91 for leaf),
   indicating reproducibility and credibility of the RNA-seq data to
   evaluate the expression of nitrogen-responsive genes in spinach.
   Histograms were produced by comparing the Log2 fold-changes by
   transcriptome analysis and RT-qPCR ([122]S10 Fig).

Key genes responding to N perturbations

   A complex network of developmental, physiological, biochemical, and
   molecular responses regulates uptake and partitioning of N. We selected
   a subset of DEGs that were directly associated with N metabolisms, such
   as N uptake, assimilation, and transport to understand their role in
   leaf or root tissues based on N availability.

Nitrate transporters

   Among the four families of transporters, Nitrate Transporter 1 (NRT1)
   [[123]42] and Nitrate Transporter 2 (NRT2) [[124]43] participate in
   NO[3]^− uptake, distribution, or storage. As per the external NO[3]^−
   concentration, plants use two specialized transport systems [[125]44],
   a constitutive, low-affinity transport system (LATS) and an inducible,
   high-affinity transport system (HATS) to maximize N uptake efficiency
   [[126]45–[127]47]. The LATS allows transport during high external
   NO[3]^− while low concentrations activate HATS [[128]44], although it
   also contributes to fulfilling NO[3]^− demand at higher concentrations
   [[129]48]. We first looked at differentially expressed NO[3]^−
   transporters under HATS induced in roots in response to N stress. We
   found at least three NRT2.1 (Spo03988, Spo09968, Spo03990) members from
   HATS were induced in roots due to N stress. In particular, the
   expression of Spo03988 (NRT2.5) and Spo09968 (NRT2.1) were upregulated
   4.3 and 2.6-fold higher in N stressed roots. The affinity of NRT2
   transporters for now NO[3]^− influx in the root has already been
   demonstrated for NRT2.1 and NRT2.5 [[130]43, [131]49, [132]50]. NRT2.5
   is a part of HATS that enables roots to absorb NO[3]^− under limited
   availability, contributes to phloem loading in shoots, and activates
   nitrate-inducible genes [[133]49, [134]51]. Both; NRT2.1 and NRT2.5 in
   Arabidopsis have been characterized as NO[3]^− responsive genes
   [[135]52]. Additionally, expression of Spo15883 (NRG2, Nitrate
   regulatory gene) was also induced in roots due to limited N
   availability. These genes are required for NO[3]^− signaling and
   regulate expression of the nitrate-responsive genes, including NRT2.1
   [[136]53]. Under limited N, we also observed the induction of three
   ammonium transporters (Spo24676, Spo16971, Spo24677), implying
   alternative strategies deployed by roots to accomplish N requirement.
   Under high N ([137]S3 Table), the expression of NRT1 transporters
   (Spo05548, Spo25744, Spo07601) was significantly upregulated in roots.
   Protein NRT1/ PTR FAMILY 6.2 has been shown to function as a
   low-affinity proton-dependent NO[3]^− transporter [[138]54], justifying
   its induction in spinach roots. While NPF2.13 (NRT1/ PTR FAMILY 2.13)
   is involved in phloem loading and NO[3]^− remobilization [[139]55]. It
   has been proposed that as the higher NO[3]^− concentrations in soil are
   first encountered by the advancing root tip, root tip specific NRT1
   transporters would best capture the available NO[3]^− [[140]56].
   Intriguingly, in leaf tissue, expression of six members of the NRT1
   family (Spo17357, Spo25690, Spo01109, Spo19480, Spo20248, and Spo26413)
   under low N and two members (Spo24247 and Spo25782) under high N were
   significantly upregulated. It has been demonstrated that several
   NO[3]^− transporters identified in leaves are closely correlated with,
   e.g., stomatal opening [[141]57], NO[3]^− reductase activity [[142]52],
   accumulation and remobilization of NO[3]^− [[143]55].

Amino acid metabolism

   In plants, inorganic N (NO[3]^− and NH[4]^+) is incorporated into Glu
   and Gln in roots before converting to other amino acids and nitrogenous
   compounds. Amino acids are later transported to the sink organs by the
   long-distance transport systems of the plant. We found that the
   expression of a putative glutamine amidotransferase GAT1 (Spo00216) was
   almost 5 -fold (Log2) up in roots as well as 4 -fold up in leaf under
   high N. Root specific GAT1 gene has been shown to be highly responsive
   to N status in Arabidopsis [[144]58]. It is plausible to assume that
   under high N conditions, GAT1 facilitates the conversion of Gln to Glu
   for channeling C to the TCA cycle. We found two Glutamate receptor
   (GLR) genes (Spo15921 and Spo01837) up-regulated in response to high N
   in roots. GLR genes are implicated in C: N signaling, hypocotyl
   elongation, root growth, stress responses, and general ion transport
   [[145]59–[146]62]. Asn and Gln serve as the major N transport and
   storage compounds in most non-leguminous plants [[147]63]. The
   expression of Asparagine synthetases (AS) (Spo05295 and Spo08792) were
   significantly up-regulated in the roots, and leaf under high N. AS
   activity is regulated by N status [[148]60, [149]64] as it catalyzes
   the conversion of Asp and Gln to Asn and Glu in an ATP-dependent
   reaction. Aminotransferases have been suggested to play essential roles
   in redirecting N resources to different pathways [[150]65]. The
   expression of the branched-chain amino acid (BCAA) aminotransferase
   genes BCAT-2 (Spo14673 and Spo13997) regulating conversion of
   leucine/valine/isoleucine to Glu were significantly upregulated in
   roots. Among the most enriched pathway terms, several DEGs associated
   with degradation of branched-chain amino acids were significantly
   upregulated in the root (33 out of 48 reference genes annotated to the
   BCAA degradation pathway) as well as leaf (23 out of 48) tissues.
   Although BCAA catabolism likely has many functions in plants
   [[151]66–[152]68], its involvement in N redistribution needs further
   functional characterization of the genes involved in BCAA catabolism.
   Additionally, we found expression of serine decarboxylase (SDC;
   Spo13254) significantly upregulated in response to high N in roots.
   Most recently, SDC like gene showing remarkably high expression in
   young roots under sufficient N has been characterized in tea plants
   [[153]69]. We performed a BLASTP homology search of the NCBI database.
   We found that the deduced protein sequence of the spinach SDC gene
   shared high similarity with characterized Tea Camellia sinensis (74.6%
   identity) SDC as well as other plant species such as Arabidopsis
   (77.26% identity), and rice (74.7% identity) ([154]S11 Fig).
   Decarboxylation of Ser is the major source of ethanolamine production
   required for plant growth [[155]70]. Moreover, it has been proposed
   that Ser could serve as a source of N in non-photosynthetic tissues
   like roots [[156]71]. Nevertheless, additional functional
   characterization of this gene will be required to understand its role
   in N metabolism. Under high N, expression of putative amino acid
   transporters such as WAT-1 related (UMAMIT34, Usually Multiple Acids
   Move In and out Transporters; Spo14333), vacuolar amino acid
   transporters (Spo21518 and Spo21518) and ABC transporter C family
   member (Spo05865) were significantly upregulated (≥2 -fold (Log2) in
   leaf implying increased N assimilation.

Nucleic acid metabolism

   The development and proliferation of living cells require DNA
   replication, which is responsible for genome duplication in plants
   [[157]72]. DNA replication is initiated and facilitated by the
   formation of the replication fork comprising the Minichromosome
   Maintenance MCM(2–7) protein helicase complex [[158]73]. The MCMs are
   found to be highly expressed in dividing tissues such as shoot apex and
   root tips of several plants [[159]74], and their expression is
   coordinated during plant development, possibly at the level of
   transcription [[160]75]. Consistent with these reports, nearly half of
   the DEGs (24 out of 50; P < 0.00029) associated with DNA replication
   fork assembly were up-regulated in leaf tissue under high N. Since DNA
   replication is linked to cell cycle progression and DNA repair
   processes, upregulation of genes associated with MCM complex, RFAs,
   genes regulating DNA polymerase complexes (α primase, δ, and ε), DNA
   ligases (Lig 1, Fen1) and helicase (Dna2) were not unexpected ([161]S12
   Fig). Besides replication fork, upregulation of several DEGs associated
   with pyrimidine (34 out 116 reference genes annotated for the
   pyrimidine pathway; P <0.004; [162]S13 Fig) and purine metabolism (41
   out of 158 reference genes annotated for the purine pathway; P <0.006;
   [163]S14 Fig) also supported increased nucleic acid synthesis and plant
   growth in response to N availability.

Hormonal pathways

   Several studies have demonstrated that NO[3]^− can regulate
   biosynthesis, degradation, transport, and signaling of phytohormones
   [[164]76, [165]77]. The auxin transport and signaling components are
   critical for altering root architecture in response to N availability
   in plants [[166]78, [167]79]. A relationship between the availability
   of NO[3]^− and auxin metabolism has been validated in Arabidopsis,
   maize, and rice [[168]80–[169]82]. Among the DEGs associated with
   hormonal pathways, 11 genes coding Auxin-binding proteins (ABP) with
   putative auxin receptor function were up-regulated in roots under high
   N treatment. Although these proteins share homology with many ABP19 and
   ABP20 proteins from plants, a detailed functional characterization
   would be required to understand their significance during excess N
   availability. The Auxin Binding Protein 1 (ABP1) has been extensively
   studied as a candidate auxin receptor and has been shown to be a key
   regulator for auxin-mediated responses [[170]83, [171]84]. Due to N
   stress, two GA2oxs; Gibberellin 2-beta-dioxygenase (Spo11903 and
   Spo03697) and a Gibberellin-regulated (Spo07806) genes were
   up-regulated in roots. GA2oxs regulate plant growth by inactivating
   endogenous bioactive gibberellins. The gibberellin's impact on
   leaf‐growth and enhanced apical dominance is well known; however,
   little information is available about root-specific changes in
   gibberellin in response to N stress. Nonetheless, reduction in
   endogenous GA concentrations in roots due to the down-regulation of GA
   biosynthesis at the transcriptional level has been shown to regulate
   root growth and NO[3]^− uptake in cucumber plants [[172]85].

Carbon metabolism

   In high N treatment, several genes associated with C metabolism were
   differentially expressed in leaf tissue. In particular, the expression
   of Phosphoenolpyruvate carboxykinase (PEPCK; Spo20411) was up-regulated
   by 4-fold (Log2 scale). PEPCK is involved in multiple C
   metabolism-related functions such as CO[2] assimilation, relocating the
   consumed TCA intermediates in the biosynthesis, and assimilation of N
   [[173]86]. Similarly, expression of Isocitrate lyase (ICL; Spo13898)
   and Malate synthase (MS; Spo16696) involved in the glyoxylate pathway
   were up-regulated (>2-fold Log2). It has been proposed that cellular
   lipids can be metabolized via the glyoxylate cycle. Sucrose can be
   synthesized from the four C products of the glyoxylate cycle, which
   converts TCA cycle components to sucrose via the process of
   gluconeogenesis. Although detailed enzyme kinetics and tracer
   experiments will be required to understand the role of these genes in
   spinach, unusually coordinated expression of PEPCK with MS and ICL
   poses a possibility of the presence of gluconeogenic pathway in leaf
   tissue. In C4 plants, like maize, gluconeogenesis has been shown to
   regulate the import/mobilization of nitrogenous compounds.
   Up-regulation of Expansins (Spo03956 and Spo11774) and Sugar
   transporter proteins (MST4; Spo11108 and Spo11108) and Beta-amylase 2
   (Spo10422) also supports cell expansion and enhanced influx of
   photo-assimilates in leaf tissue under high N treatment, respectively.

Differentially expressed transcription factors (TFs) responding to N stress

   We used iTAK [[174]87] to perform the transcription factor analysis. A
   total of 786 transcription factor related genes were identified from
   spinach and the 56 TF families that expressed genes differentially
   between leaf Vs. root tissues under low N or high N are shown in
   [175]Fig 9. The largest members of the TFs belonged to the bHLH family
   (48), followed by MYB (33), AP2-EREBP (31), WRKY (29), NAC (19), in
   case of high N treatment. While for low N, the bHLH family (45),
   followed by MYB (36), AP2-EREBP (27), WRKY (28), NAC (17) members of
   the families, were differentially regulated. The distribution of the
   differentially expressed TFs between leaf vs. root tissues under the N
   stress was separately investigated. A list of differentially expressed
   TFs showing at least ≥ 2 Log2 fold significant (p-value < 0.05) change
   in the expression in leaf and root tissues due to N treatments ([176]S7
   Table). It will be valuable to functionally characterize the identified
   TFs for enhancing plant performance under limited N availability and
   thereby providing novel targets for genetic manipulation.

Fig 9. Transcription factor analysis.

   [177]Fig 9
   [178]Open in a new tab

   The graph shows the number of transcription factors (vertical axis)
   from 56 TF families that expressed genes differentially (change in the
   expression <2 Log2-fold) in leaf tissue under Low-N or High-N. iTAK is
   used to perform the transcription factor analysis of plants (p-value
   < 0.05). The details of TFs families are published in Zheng et al.,
   2016.

   In roots, low N up-regulated expression (≥ 2-fold log2) of three
   ERFs—TINY (Spo16233), ERF17 (Spo20470), and ERF12 (Spo10826); and two
   NAC TFs (NAC43; Spo25926 and NAC86; Spo15445). The interaction between
   ethylene and availability of N affects numerous physiological
   processes, including roots architecture, leaf, reproductive organ
   development, as well as the synthesis of amino acids, proteins, and
   enzymes [[179]88]. N starvation increases the number or affinity of
   root receptors, allowing roots to respond to lower concentrations of
   ethylene than those found in unstressed roots. The AP2/ERF-TF family
   has been involved in signaling processes and the responses to
   environmental stresses [[180]89], such as micronutrient deficiency
   [[181]90, [182]91]. Similar to our observations, under N starvation,
   seven ERF genes in cucumber seedlings were also differentially
   regulated [[183]92]. Under high N, the expression of
   Ethylene-responsive transcription factor 13 (ERF13; Spo11940) and two
   of the NAC29 TFs was induced (≥ 2-fold Log2). We found that the
   expression of 1-aminocyclopropane-1-carboxylate oxidase1 (ACO; Spo03782
   and Spo03777) genes involved in ethylene synthesis were also
   upregulated in roots. Induction of expression and enzymatic activities
   of ethylene biosynthetic genes, ACO, and ACC synthase (ACS) in response
   to high N has been demonstrated in other plants [[184]93–[185]95].

Conclusions

   In conclusion, comparative analysis of the leaf and root transcriptomes
   indicated a coordinated regulation of multiple genes and pathways and
   provided a catalog of transcriptomic variation between the spinach
   tissues in response to N. Although transcriptomic responses to N
   perturbation have been studied in Arabidopsis and few other crops, to
   the best of our knowledge, this is the first study that provides
   information about the partitioning of leaf and root transcriptomic
   responses in spinach. The genomic resources generated from this RNA-Seq
   experiment can be used to characterize genes associated with NUE in
   spinach breeding materials and to develop informative markers
   specifically to select parental lines. RNA-seq analysis revealed that N
   stress caused most transcriptomic changes in roots, identifying 1,346
   DEGs compared to 1136 DEGs in leaf. Significant upregulation of NRT2.1
   (Spo03988, Spo09968, Spo03990) in roots validated its role in N
   acquisition under limited availability. The root-specific genetic
   manipulation to alter the expression of NRT2.1 could be a way forward
   to enhance NUE in spinach. Concomitant activation of putative glutamine
   amidotransferase GAT1 (Spo00216) and Asparagine synthetases (AS)
   (Spo05295 and Spo08792) along with Glutamate receptor (GLR) (Spo15921
   and Spo01837) in response to excess N highlights the functional
   significance of Gln reduction in N assimilation. Consistent with the
   effects of higher N availability on plant growth, expression of genes
   associated with nucleic acid synthesis and DNA replication were
   significantly upregulated. The RNA-Seq data generated in this study
   will be a valuable resource for the identification of key genes related
   to N metabolism and identifying potential targets for genetic
   manipulation to enhance N uptake and utilization. The functional
   characterization of the differentially expressed genes during N stress
   would provide new insights into developing N use efficient varieties.

Supporting information

   S1 Fig. Heat maps of the correlation coefficient between samples.

   The scatter diagrams demonstrate the correlation coefficient between
   all samples, R2, the square of the Pearson coefficient. LNR: low N
   root; LNL: low N leaf; HNR: high N root and HNL: high N leaf; The
   numbers 1,2,3 are independent replications.

   (TIF)
   [186]Click here for additional data file.^ (280.8KB, TIF)
   S2 Fig. Different gene expression levels in N treatments.

   FPKM violin Plot, the x-axis shows the sample names where T2 is high N
   and T1 low N; LF = Leaf and RT = Root, and the y-axis shows the
   log10(FPKM+1). Each violin has five statistical magnitudes (max value,
   upper quartile, median, lower quartile, and min value). The violin
   width shows the gene density.

   (TIF)
   [187]Click here for additional data file.^ (170.1KB, TIF)
   S3 Fig. Hierarchical cluster analysis of differential gene expression
   in the leaf and root of spinach plants at high N and low N.

   (TIF)
   [188]Click here for additional data file.^ (132.9KB, TIF)
   S4 Fig. The volcano maps showing the number of differentially expressed
   genes in spinach plants grown under High-N.

   Red dots represent up-regulated genes, and green dots represent
   down-regulated genes (p < 0.05).

   (TIF)
   [189]Click here for additional data file.^ (120.4KB, TIF)
   S5 Fig

   Gene ontology (GO) enrichment bar chart of DEGs showing most enriched
   GO terms of (A) upregulated and (B) downregulated DEGs in root tissue
   under Low-N. The y-axis is the enriched GO term, and the x-axis is the
   number of DEGs enriched in the listed term. Colors represent different
   GO types: biological process, cellular component, and molecular
   function. * significantly enriched term (p-value < 0.05).

   (TIF)
   [190]Click here for additional data file.^ (222.4KB, TIF)
   S6 Fig

   Gene ontology (GO) enrichment bar chart of DEGs showing most enriched
   GO terms of (A) upregulated and (B) downregulated DEGs in leaf tissue
   under low N. The y-axis is the enriched GO term, the x-axis is the
   number of DEGs enriched in the listed term. Colors represent different
   GO types: biological process, cellular component, and molecular
   function. * significantly enriched term (p-value < 0.05).

   (TIF)
   [191]Click here for additional data file.^ (235KB, TIF)
   S7 Fig. Gene ontology (GO) enrichment bar chart of DEGs showing most
   enriched GO terms of (A) upregulated and (B) downregulated DEGs in
   Low-N in both tissues.

   The y-axis is the enriched GO term, and the x-axis is the number of
   DEGs enriched in the listed term. Colors represent different GO types:
   biological process, cellular component, and molecular function. *
   significantly enriched term (p-value < 0.05)

   (TIF)
   [192]Click here for additional data file.^ (230.3KB, TIF)
   S8 Fig

   KEGG pathway enrichment scatter-plot analysis of DEGs in low N (A) and
   high N (B) in leaf Vs. root tissue. The Rich factor is the ratio of
   DEGs in this pathway term to all gene numbers annotated in this pathway
   term. A q value is the corrected p-value ranging from 0 to 1, and a
   lower value indicates greater pathway enrichment. The pathway names are
   shown on the vertical axis, rich factor on the horizontal axis, the
   size of the point represents the number of DEGs, and the color of the
   dot represents the q value

   (TIF)
   [193]Click here for additional data file.^ (210.7KB, TIF)
   S9 Fig. Linear regressions involving the RNA sequencing data and the
   RT-qPCR validation data expressed in terms of Log2 fold-change (FC).

   FC was calculated as the ratio between the drought-stressed and control
   plants. (A, B) indicate roots and leaves, respectively. * Significant
   Pearson’s correlation coefficient P < 0.01. (R2 = 0.85 and 0.80 for
   root and leaf tissues respectively).

   (TIF)
   [194]Click here for additional data file.^ (169.5KB, TIF)
   S10 Fig. qRT–PCR validation of RNA-seq results.

   Twenty randomly selected DEGs identified by RNA-seq (red bars) in
   roots, and leaf tissue of spinach were selected to analyze by qRT-PCR
   (black bars). Comparison of the fold change of RNA-seq and qRT–PCR
   shows a correlation coefficient of 0.94, indicating the reliability of
   RNA-seq results. Error bars represent the SEM.

   (TIF)
   [195]Click here for additional data file.^ (143.4KB, TIF)
   S11 Fig. CLUSTAL W (1.83) multiple sequence alignment of spinach SDC.

   Multiple alignments of the protein sequence of spinach SDC (Spo13254)
   with selected SDCs from other species. Identical amino acids are shown
   with asterisk *. The NCBI GeneBank IDs used to Serine decarboxylases
   from plants are as follows: Arabidopsis ([196]NP_175036.1), Camellia
   sinensis ([197]XP_028123129.1), Solanum lycopersicum
   ([198]XP_004237774.1), Oryza sativa Japonica Group ([199]BAS79103.1),
   Zea mays ([200]PWZ22363.1), Spinacia oleracea (Spo13254)

   (TIFF)
   [201]Click here for additional data file.^ (1.2MB, tiff)
   S12 Fig. DNA replication fork (KEGG: 03030).

   DNA replication pathway diagram highlighting DEGs up-regulated (green
   boxes with red border) in leaf tissue under high N.

   (TIF)
   [202]Click here for additional data file.^ (115.3KB, TIF)
   S13 Fig. Pyrimidine metabolism (KEGG: 00240).

   Diagram highlighting DEGs up-regulated (green boxes with red border) in
   leaf tissue under high N.

   (TIF)
   [203]Click here for additional data file.^ (297.4KB, TIF)
   S14 Fig. Purine metabolism (KEGG: 00230).

   Diagram highlighting DEGs up-regulated (green boxes with red border) in
   leaf tissue under high N.

   (TIF)
   [204]Click here for additional data file.^ (353.2KB, TIF)
   S1 Table. Primers used for RT-qPCR.

   (XLSX)
   [205]Click here for additional data file.^ (18.3KB, xlsx)
   S2 Table. Detail statistics of sequencing results.

   Total mapped—number of reads that can be mapped to the reference
   genome; Multiple mapped—number of reads mapped to multiple sites in the
   reference genome. Uniquely mapped—number of reads mapped to the
   reference genome; the number of reads mapped to the positive strand (+)
   or the minus strand (-); Splice reads mapped to two exons (junction
   reads), and non-splice reads are mapped entirely to a single exon.HN =
   high N; LN = low N, L1, L2 L3 = Leaf replications 1, 2,3; R1, R2, R3 =
   Root Replications 1,2,3.

   (XLSX)
   [206]Click here for additional data file.^ (18.9KB, xlsx)
   S3 Table. List of differentially expressed genes upregulated in roots
   under high N treatment.

   (XLSX)
   [207]Click here for additional data file.^ (63.2KB, xlsx)
   S4 Table. List of differentially expressed genes upregulated in roots
   under low N treatment.

   (XLSX)
   [208]Click here for additional data file.^ (59.6KB, xlsx)
   S5 Table. List of differentially expressed genes upregulated in leaf
   under high N treatment.

   (XLSX)
   [209]Click here for additional data file.^ (61KB, xlsx)
   S6 Table. List of differentially expressed genes upregulated in leaf
   under low N treatment.

   (XLSX)
   [210]Click here for additional data file.^ (59KB, xlsx)
   S7 Table. Differentially expressed transcription factors (TFs).

   TFs with <2-fold change (Log2) in the expression due to N stress roots
   and leaf tissues (P-value < 0.05). The details of TFs families are
   published [[211]87].

   (XLSX)
   [212]Click here for additional data file.^ (11.3KB, xlsx)

Acknowledgments