Populus tomentosa clones (TC1521) were grown in culture on a half-strength Murashige–Skoog (MS) medium (Murashige and Skoog, 1962) (pH = 6.2) containing 20 g L^–1 sucrose (Sigma-Aldrich, St. Louis, MO, United States) and 0.4 mg L^–1 indole-3-butyric acid (IBA) (Sigma-Aldrich) at 25°C under a 16/8 h (day/night) photoperiod. Sixty-day-old plants were transferred into a hydroponic solution with sufficient N level for 5 days, which was changed for fresh solution every 2 days. The plants were then transferred to a solution with or without sufficient N as the control and treatment groups for 3 days as described previously (Ren et al., 2015). Briefly, plants were grown in modified half-strength mass spectrometry (MS) liquid medium (pH = 6.2) with 2 mM NH₄NO₃ (Sigma-Aldrich) and 1 mM KNO₃ (Sigma-Aldrich) as sufficient N conditions (KK) (control) or with.01 mM NH₄NO₃ and 1 mM KCl (Sigma-Aldrich) instead of KNO₃ for low-N treatment (DN). Whole P. tomentosa plants were harvested in the midmorning, immediately frozen in liquid N, and stored at −80°C.

Samples were taken from six 60-day-old P. tomentosa plants with or without low-N treatment and ground in liquid N, and 50 ± 2.5 mg materials were transferred to 1.5-mL tubes, followed by the addition of 1 mL of 100% methanol (precooled to −20°C) and 10 mL of phenylalanine (10 μg/mL) as an internal standard, and centrifuged for 10 s. The tubes were preheated, ultrasonicated for 15 min at 60°C, and then centrifuged for 10 min. The supernatants (0.4 mL) were then transferred to 0.2 mL of acetonitrile precooled to 0°C, and 0.4 mL ultrapure water was added to the new tubes and then centrifuged for 15 min. Then, aliquots of 200 μL of the supernatants were transferred to glass bottles and dried under a gentle stream of N₂ gas. Methoxyamine pyridine hydrochloride at a concentration of 20 mg/L (30 μL) was added to the bottles and shaken for 30 s. The oximation reaction proceeded at 37°C for 15 min. Finally, 30 μL of N,O-bistrifluoroacetamide (containing 1% trimethylchlorosilane) derivatization reagent was added and allowed to react for 1 h at 70°C. After these reactions, the samples were analyzed for their metabolite contents.

Metabolites were detected by gas chromatography-mass spectrometry (GC-MS) (7890A/5975C GC-MS system; Agilent Technologies, Santa Clara, CA, United States) at Shanghai Sensichip Infotech Co. Ltd. (Shanghai, China). The Restek capillary column was an HP-5 ms (30 m × .25 mm × .25 μm) (Agilent Technologies). The parameters were as follows: injection port temperature, 280°C; EI ion source temperature, 230°C; quadrupole rod temperature, 150°C; carrier gas, high-purity helium (99.99%); splitless injection; and sample size, 1 μL. The temperature program consisted of an initial temperature of 70°C for 2 min, 10°C/min up to 320°C, and was put on hold for 6 min. GC-MS was performed by the full-scan method with a range from 50 to 550 mass-to-charge. The XC/MS software was used for metabolomics data preprocessing in the R software package (R Development Core Team, Vienna, Austria) and then compiled to remove impurity peaks due to losses from the column and the sample preparation process. The results were then organized as a two-dimensional (2D) matrix, including observation values (samples), variables (retention time/mass-to-charge ratio), and peak strength. Finally, each sample was normalized relative to the total mass using the internal standard, and the normalized data were input into SIMCA-P (ver. 11.0) for principal component analysis (PCA) using the PLS-DA model with variable importance in projection (VIP) values > 1, combined with Student’s t-test (p ≤ 0.05) to identify the differentially expressed metabolites, and searched for metabolites in the National Institute of Standards and Technology¹ and Kyoto Encyclopedia of Genes and Genomes (KEGG)² database.

Before the data analysis, all data were standardized for mean-centering and unit-variance scaling using SIMCA-P with the default parameters (ver. 11.5³). Hierarchical clustering analysis (HCA) and PCA models were tested using all samples. Significant differences among metabolites between DN and KK were examined using the t-test (p ≤ 0.05). A heatmap was built using Pearson’s test and hierarchical clustering, performed with MATLAB 7.5 (MathWorks, Inc., Natick, MA, United States). Metabolic pathways were constructed with Metaboanalyst, and the Arabidopsis metabolic pathway database was used as a reference for the global algorithm. The enrichment pathways of metabolites were analyzed based on the KEGG database with a p-value ≤ 0.05 established as the false discovery rate (FDR) for multiple tests. The interactions among different metabolites were determined using KEGGSOAP⁴ and metabolic pathway networks were constructed using Cytoscape⁵.

Total RNA was isolated from three 60-day-old P. tomentosa plants with or without low-N treatment using TRIzol reagent (Invitrogen, Carlsbad, CA, United States). The quantity and quantity of total RNA were determined using 1% agarose gel (Sigma-Aldrich) electrophoresis and an Agilent 2100 Bioanalyzer (Agilent Technologies). Complementary DNA (cDNA) libraries were constructed as described previously (Zhang et al., 2012). Briefly, total RNA was first treated with DNase I, and mRNAs were then enriched with oligo(dT) magnetic beads mixed with fragmentation buffer (Ambion, Austin, TX, United States). The fragmented messenger RNAs (mRNAs) were used to synthesize the random hexamer-primed cDNA, which was subjected to size selection and further PCR amplification. The quantity and quantity of cDNA libraries were determined with the Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, United States). Finally, the cDNA libraries were sequenced using the Illumina HiSeq 2000 platform (Illumina, San Diego, CA, United States), and raw sequencing reads were processed to remove the dirty reads, i.e., reads with adapters, >5% unknown nucleotides, and low-quality reads, based on the National Institutes of Health Sequence Read Archive database (accession number: SRP063920). The obtained clean reads were applied to de novo assembly with Trinity⁶, which contains three independent software modules: Inchworm, Chrysalis, and Butterfly (Supplementary Figure S1). Briefly, the programs in Trinity were applied sequentially to assemble the clean reads into unique full-length transcripts, map reads into contigs, and then assemble unigenes. Finally, unigene sequences were aligned by BLASTn to the National Center for Biotechnology Information (NCBI) non-redundant nucleotide database (NT)⁷ with an e-value cutoff < 10^–5. Furthermore, unigenes were aligned with the NCBI non-redundant protein database (NR) (see footnote 7) and Swiss-Prot⁸ protein database with an e-value cutoff < 10^–5. To predict and classify possible functions, unigenes were also searched against the Cluster of Orthologous Groups (COG) database⁹ by BLASTx with an e-value cutoff < 10^–5. Unigenes not aligned to any database were further scanned by ESTScan to obtain the amino sequences of hypothetical proteins. By blasting against these databases, four parts of the analysis of unigenes were performed: SSR analysis, unigene expression annotation, single nucleotide polymorphisms (SNP) analysis, and unigene function annotation. Further, PCA analysis and unigene expression difference analysis were conducted based on unigene expression annotation. The fragments per kilobase per million reads method was used to assess the differential expression of unigenes, as described previously (Chen et al., 2013; Hou et al., 2016), and unigenes with FDR ≤ 0.001 and |log2Ratio| ≥ 1 were regarded as differentially expressed genes (DEGs) between the DN and KK groups.

Unigenes were aligned with Kyoto Encyclopedia of Genes and Genomes (KEGG)² to predict the metabolic pathways, with p ≤ 0.05 and q ≤ 0.05 taken to indicate significant enrichment (Hou et al., 2016; Wang X. et al., 2016), and the functions of gene products (using BLASTx¹⁰) with an e-value cutoff < 10^–5. Furthermore, Gene Ontology (GO) analysis was performed with NR annotation to annotate the functions of unigenes with the Blast2GO program¹¹, with p ≤ 0.05 taken to indicate significance (Hou et al., 2016; Wang X. et al., 2016).

Samples were taken from three 60-day-old P. tomentosa plants with or without low-N treatment and analyzed for the concentrations of indole-3-acetic acid (IAA), abscisic acid (ABA), and gibberellic acid (GA). Extraction and purification were performed as described previously (Pan et al., 2010). Briefly, the samples were ground into a powder with a mortar and pestle, and 50-mg samples were transferred to precooled 2-mL tubes and kept in liquid N. Then, 500 μL of 2-propanol/H₂O/concentrated HCl (2:1:0.002, vol/vol/vol) extraction solvent was added to each tube, and various volumes of internal standard solutions were added. The tubes were centrifuged at 100 rpm for 30 min on a shaker at 4°C. Then, 1 mL of dichloromethane was added to each sample and shaken for 30 min at 4°C. The tubes were further centrifuged at 13,000 × g for 5 min at 4°C. Then, 900 μL of solvent from the lower phase was transferred into screw-capped vials and concentrated using an N evaporator. The samples were redissolved in 100 μL of methanol. Then, 50 μL of sample solution was injected into the reverse-phase C18 Gemini high-performance liquid chromatography (HPLC) column for high-performance liquid chromatography/electrospray ionization tandem mass spectrometry (HPLC-ESI-MS/MS) analysis. Quantitative analysis of each plant hormone was performed as described previously (Pan et al., 2010).

Samples were taken from three 60-day-old P. tomentosa plants with or without low-N treatment for 0, 1, 3, or 5 days as described above. Total RNA was extracted using an RNA prep Pure Plant Kit (Tiangen, Beijing, China) and reverse transcribed using a FastQuant RT Kit (With gDNase) (Tiangen). To verify RNA sequencing (RNA-Seq) profiles in this investigation, Real-Time Quantitative Reverse Transcription PCR (qRT-PCR) was performed on an Applied Biosystems 7500 Fast Real-Time PCR System using SYBR Premix Ex Taq™ (Tli RNaseH Plus; Takara, Shiga, Japan) in accordance with the manufacturer’s instructions. The primers of 14 randomly selected DEGs are listed in Supplementary Table S1. Reactions were performed in a volume of 10 μL containing 1 μL of cDNA, 5 μL of SYBR Green, 0.2 μL of forward primer, 0.2 μL of reverse primer, 0.2 μL of ROXII, and 3.4 μL of distilled water. The thermocycling conditions consisted of an initial denaturation step at 95°C for 10 min followed by 40 cycles of 95°C for 30 s, 95°C for 5 s, and 60°C for 30 s. All reactions were performed in triplicate for each gene. The 2^–ΔΔCT relative quantification method was used to evaluate and calculate variations (Livak and Schmittgen, 2001), with 18S rRNA used as an internal reference. The correlations of gene expression between RNA-Seq and qRT-PCR were analyzed by Pearson’s test with p ≤ 0.01 taken to indicate significance.

To determine the metabolomic regulatory mechanisms of the response of poplar to N deficiency, GC-MS was performed to analyze the changes in metabolites between P. tomentosa plants grown under conditions of low and sufficient N (DN and KK, respectively). The total ion chromatogram (TIC) is shown in Supplementary Figure S2. A total of 1,131 metabolites were finally identified by GC-MS (Supplementary Table S2), and partial least squares-discriminant analysis (PLS-DA) and PCA were performed to determine the accuracy and significance of differences between the KK and DN samples (Supplementary Figure S3). The KEGG enrichment analysis assigned the detected metabolites to 18 metabolic pathways, including amino acid- and sugar-related metabolism. Metabolites clustered in C fixation, starch and sucrose, and fructose and mannose metabolism were also enriched (Supplementary Table S3).

Among the 1,131 metabolites, 30 with significantly differential expression were identified by PLS-DA with VIP > 1 and p ≤ 0.05 (Figure 1). Most of these metabolites (70%) were significantly suppressed under conditions of low-N stress. Marked reductions in the levels of metabolites with amino groups (such as valine, L-isoleucine, L-alanine, cadaverine, ethylamine, and ethanolamine) were assumed to result in the cessation of the de novo synthesis of the free amino acids because of N deficiency (Table 1). Furthermore, the levels of alcohols were reduced under low-nitrate stress, as shown by the a-hydroxycyclohexene, cyclohexanol, ethylene glycol, inositol, and xylitol contents compared with the controls. However, N deficiency led to the accumulation of some soluble sugars, including D-fructose, D-galactose, and D-glucose, in P. tomentosa plants.

Based on the pathways of low-N-responsive metabolites, we further constructed metabolic networks to analyze their interactions (Figure 2). As shown in the metabolic networks, the induced sugars, including D-glucose, produced by the Calvin cycle influenced the synthesis of diverse downstream metabolites, including amino acids and putrescine belonging to the citrate cycle. Interestingly, in addition to the induced D-glucose, D-fructose, and D-galactose, other sugar-related metabolites, such as galactonic acid and oxalic acid, were enhanced under conditions of N deficiency in this study. The induced level of sugar may have been due to reduced carbohydrate metabolism during degradation and utilization by downstream metabolites of amino acids derived from N assimilation, in turn, due to low-N stress, consistent with the pattern of N assimilation-related gene regulation revealed by transcriptomics analysis, which also revealed the connection of C and N metabolism.

To further determine the genes involved in regulating low-N signaling in poplar under conditions of N deficiency, RNA-Seq was conducted in DN and KK P. tomentosa plants. After discarding the contaminated raw data, 104,843,476 clean reads (52,903,032 reads from DN, 51,940,444 reads from KK) containing 943,591,2840 nt were obtained (Supplementary Table S4). Based on these clean reads, 154,857 and 146,823 contigs of the treatment and control, respectively, were further assembled to 71,801 (DN) and 82,908 (KK) unigenes with mean lengths of 582 and 508 nt, respectively (Supplementary Table S5).

The unigene sequences were then aligned to the NCBI NR protein database, Swiss-Prot, KEGG, and COG by BLASTx with an e-value cutoff < 10^–5, and to the NCBI NT nucleotide database by BLASTn with an e-value cutoff < 10^–5, to retrieve proteins with the highest sequence similarity to the given unigenes along with their protein functional annotations (Supplementary Figure S1). For function annotation analysis, we obtained 52,816, 54,584, 32,342, 29,222, 18,431, and 42,264 unigenes, which were annotated to the NR, NT, Swiss-Prot, KEGG, COG, and GO databases, respectively, and a final total of 59,125 unigenes were annotated (Supplementary Tables S6, S7).

Then, we searched the unigene sequences against the COG database to predict the possible functions and understand the global gene function distributions of the species. A total of 18,431 sequences from 59,125 unigenes were mapped to the COG database; among the 25 categories, the largest group of unigenes (n = 5,411, 29.36%) were annotated as “General function prediction only,” followed by “Translation, ribosomal structure and biogenesis” (n = 3,384, 18.36%) and “Transcription” (n = 2,840, 15.41%), while genes functioning in “Nuclear structure” (n = 2, 0.01%) and “Extracellular structures” (n = 17, 0.09%) represented the smallest categories (Figure 3).

To identify the genes involved in the response of P. tomentosa to low-N stress, we further screened unigenes with significant changes based on the assembled data. A total of 1,561 downregulated and 1,101 upregulated DEGs were identified with a false discovery rate (FDR) ≤ 0.001 and | log2Ratio| ≥ 1 (Supplementary Table S8a and Supplementary Figure S4). To exclude statistical error, more rigorous criteria (fragments per kilobase per million reads of unigene in KK and DN > 1) were applied, and 1,017 DEGs were finally identified under conditions of N deficiency in P. tomentosa (Supplementary Table S8b). To understand the functions of these genes in the N stress response, we further aligned DEGs to the GO database. The results revealed 21, 15, and 12 classes of DEGs involved in biological process, cellular component, and molecular function, respectively (Figure 4). In the biological process, most DEGs fell into the cluster of “cellular process” (889, 1.50%), “metabolic process” (874, 1.48%), “single-organism process” (454,.77%), and “response to stimulus” (373,.63%). Three categories of DEGs, “organelle” (n = 931, 1.57%), “cell” (n = 1,072, 1.81%), and “cell part” (n = 1,072, 1.81%), were dominant in the cellular component, whereas “binding” (n = 604, 1.02%), “catalytic activity” (n = 464, 0.78%), and “structural molecule activity” (n = 342, 0.58%) represented the main DEG groups in molecular function (Supplementary Table S9).

To further understand the molecular pathways in which these P. tomentosa genes were involved, we scanned the DEGs into the KEGG database. As a result, 1,329 genes were finally annotated to 104 KEGG pathways, among which 16 pathways were significantly enriched at p ≤ 0.05, including “genetic information processing,” “cellular processes,” “metabolism,” “environmental information processing,” and “organismal systems” under conditions of low-N stress in P. tomentosa (Table 2). Intriguingly, five pathways related to amino acid metabolism were significantly enriched and 20 DEGs were involved in N metabolism, indicating that P. tomentosa shows a marked response to low-N stress. On the other hand, four pathways related to carbohydrate metabolism were also enriched, i.e., amino sugar and nucleotide sugar metabolism, inositol phosphate metabolism, pentose and glucuronate interconversion, and propanoate metabolism, suggesting that carbohydrate regulation is involved in the response to low-N stress (Table 2). In addition, many DEGs could be clustered into several metabolic pathways, including N metabolism, auxin-related pathway, or function as some transporters and kinases, etc. (Supplementary Table S10).

Plants take up N from two sources in the soil, ammonium, and nitrate. Nitrate is reduced to nitrite (NO₂^–) and then ammonium, which is incorporated into glutamate and glutamine with a net supply of 2-oxoglutarate (provided by isocitrate dehydrogenase, IDH), catalyzed by the GS/glutamate synthase (GOGAT) cycle. Other amino acids, such as aspartate and asparagine, are then generated by aspartate aminotransferase (AspAT) and asparagine synthetase (AS) (Suárez et al., 2002). In the nitrate assimilation process, N is incorporated into N-containing compounds, such as other amino acids, chlorophylls, and nucleic acids. Analysis of the transcriptome data revealed downregulation of N acquisition genes, such as putative nitrate reductase (NR; twofold downregulated) and GS (1.60-fold downregulated). The genes that participate in the biogenesis of various amino acids and other N-containing compounds from glutamine/glutamate were shown to be downregulated: putative IDH (11.60-fold downregulated), AspAT (1.17-fold downregulated), and AS (1.37-fold downregulated) (Supplementary Table S10 and Table 3).

Consistent with these observations, several genes encoding putative amino acid synthesis-related proteins were suppressed under low-N conditions, such as proline-rich protein (CL2231.Contig3_All), aspartate aminotransferase (CL8615.Contig2_All), cationic amino acid transporter (Unigene38715_All), and GS (Unigene27220_All, Unigene22953_All) (Supplementary Table S10). These results indicated the decreased level of N assimilation and further utilization under conditions of N deficiency stress in P. tomentosa.

C assimilation and N absorption control the status of plant growth and development. Extensive studies have shown the interdependence and interaction of C and N metabolism in photosynthesis, amino acid metabolism, lipid metabolism, and carbohydrate metabolism, including the TCA cycle, Calvin cycle, glycolysis, etc. (Sinha et al., 2015; Goel et al., 2016). Briefly, C metabolism is mostly dependent on photosynthesis, in which N participates in the biosynthesis of chlorophyll and other related proteins. Correspondingly, N assimilation demands energy and C skeletons from C metabolism.

Obviously, a low-N environment affects the functions of the chloroplasts and chlorophyll, and putative related proteins were downregulated in P. tomentosa, such as chlorophyll a–b binding protein of light-harvesting complex (LHC) II type I (2.73-fold downregulated), chlorophyll a–b binding protein CP26 (1.21-fold downregulated), photosystem I reaction center subunit psaK (1.37-fold downregulated), oxygen-evolving enhancer protein 2 of photosystem II (1.89-fold downregulated), ferredoxin–NR (1.05-fold downregulated), and chloroplast processing peptidase, indicating depressed photosynthesis capacity under low-N conditions (Supplementary Table S10).

On the other hand, our transcriptomics data revealed significant changes in the enzymes involved in the TCA cycle, glycolysis, and lipid degradation. Particularly, three enzymes involved in glycolysis, i.e., 6-phosphofructokinase 4, fructose-bisphosphate aldolase, and phosphoglycerate mutase, were decreased under low-N conditions. Moreover, five enzymes in the TCA cycle were altered in P. tomentosa under conditions of low-N stress, and genes encoding putative citrate synthase, succinyl-CoA synthetase, and succinate dehydrogenase were upregulated, while genes encoding putative aconitate hydratase, IDH (time-limiting enzyme), and PEPCase (flux regulator of TCA) were downregulated. Conversely, genes encoding five putative enzymes involved in lipid degradation were induced, i.e., stearoyl-CoA desaturase, 3-hydroxyacyl-CoA dehydrogenase, acetyl-CoA acyl-transferase, fatty acyl-CoA reductase, and long-chain-aldehyde dehydrogenase (Supplementary Table S10). Our data indicated decreased energy from glucose and increased energy and C skeletons from lipid, while the downregulation of putative phosphoenolpyruvate carboxylase (PEPCase) and IDH depressed the flux of oxaloacetate into the TCA cycle and lowered the production of ATP and nicotinamide adenine dinucleotide (NADH) in plants, respectively. Taken together, these results suggest reduced carbohydrate metabolism and energy production, particularly decreased energy generated from glucose in P. tomentosa, under conditions of low-N stress.

Auxin and ABA were reported to be closely related to N signaling (Sakakibara, 2003; Mockaitis and Estelle, 2008; Lu et al., 2015). The plant hormones IAA, ABA, and GA in P. tomentosa were quantified by HPLC-ESI-MS/MS in this study. The concentration of IAA was increased under low-N conditions (1.52) compared with normal conditions (1.16) (Figure 5). The accumulation of IAA was consistent with the increased expression of genes encoding auxin response genes. For example, putative ARF6 (1.50-fold increase) and IAA10 (2.23-fold increase) were markedly induced (Table 4).

On the other hand, the concentrations of ABA and GA were generally reduced, especially ABA, with more significant repression under conditions of low-N stress (Figure 5). ABA and GA are biosynthesized in the methylerythritol-4-phosphate (MEP) pathway. In the network, the conversion between isopentenyl pyrophosphate (IPP) and dimethylallyl diphosphate (DMAPP) is catalyzed by putative isopentenyl diphosphate isomerase (IDI; 12.09-fold). Next, isoprenyl pyrophosphate synthase (IPPS) catalyzes IPP to generate isopentenyl AMP (iAMP), the precursor of cytokinin (CTK). Meanwhile, IPP and DMAPP can be transformed to geranylgeranyl diphosphate (GGDP), further generating GA and other carotenoids. In addition, GGDP can be reduced to phytyl diphosphate by putative geranyl-geranyl reductase (1.12-fold), thus providing phytol (0.507) for chlorophyll synthesis (Table 4). The concentrations of GA and ABA in this pathway showed small reductions, whereas another route of phytol utilization was increased under conditions of N deficiency.

To validate the expression profiles of the identified DEGs and determine their possible dynamic responses to low-N stress at different treatment stages, the levels of expression of 14 DEGs were investigated by qRT-PCR after 0, 1, 3, and 5 days of low-N treatment (Figure 6). Except for Unigene24651, Unigene3952, CL5328.Contig5, and CL4139.Contig1, the expression changes of most (10) DEGs at 3 days under low-N treatment conditions quantified by qRT-PCR were consistent with the abundance determined by RNA-Seq analysis. The inconsistencies in the expression of the four genes quantified by these two methods may have been due to insufficient coverage sequencing depth to reflect the true distribution of these genes or to differences in the data normalization criteria for the two methods. Further expression correlation analyses of these ten DEGs demonstrated that the expression determined by RNA-Seq was positively correlated with that revealed by qRT-PCR (Supplementary Figure S5).

After further examining the expression changes by qRT-PCR, we found three different expression patterns (Figure 6). First, the expression abundance increased gradually after 1 and 3 days of low-N treatment, peaked at 3 days, and then declined steadily at 5 days for Unigene 22453, Unigene19959, Unigene25149, and Unigene24651. It is worth noting that Unigene24651 was markedly repressed at 5 days relative to 0 days, in contrast to Unigene 22453, Unigene19959, and Unigene25149. Second, the expression levels of Unigene3952, CL5328.Contig5, Unigene24078, Unigene3286, CL933.Contig 1, and CL9290.Contig2 peaked after 1 day of low-N treatment and decreased after 3 days. The expression levels of Unigene3952, CL5328.Contig5, and Unigene24078 were increased, while those of CL933.Contig 1 and CL9290.Contig2 were reduced, after 5 days. Both of these patterns showed the dynamic fluctuation of miRNA expression in P. tomentosa in response to low-N stress at different stages. Third, Unigene18608, CL6853.Contig2, and CL4139.Contig1 were all induced, while Unigene13209 was repressed, at all stages of low-N stress.

Based on the global metabolome, we found an array of primary metabolites with differential changes under DN treatment, meanwhile, the transcriptome data showed accompanied changes of genes, which were involved in biosynthesis or degradation of these key metabolites. Thereafter we compared the expression of inducement or suppression of genes with these corresponding metabolites. Coordinated metabolomics and transcriptomics analyses revealed the relative expression of five pairs of metabolites and unigenes (Table 5). Our results revealed five genes encoding enzymes involved in the biogenesis or degradation of associated carbohydrates, and genes and metabolites showed a co-expression pattern. The data indicated that the accumulation of D-glucose and galactonic acid could be attributed to their increased biosynthesis, catalyzed by dTDP-D-glucose 4,6-dehydratase and β-galactosidase, respectively, while the accumulation of D-fructose and D-galactose may have been the result of a reduction of fructose-2,6-bisphosphatase and galactose mutarotase expression.

The results of the present study showed that N deficiency leads to the differential expression of many genes, depresses N assimilation, reduces amino acid biosynthesis, and reduces photosynthetic and energy production capacity. Similar to previous studies, a large number of DEGs were identified in P. tomentosa in response to low-N conditions, including several genes directly related to N metabolism, encoding putative nitrate reductase (NR), glutamate dehydrogenase (GDH), glutamine synthetase (GS), among others (Supplementary Table S10). Generally, plants can directly use ammonium, but not nitrate (Jackson et al., 2008). Nitrate first must be reduced by NR (Bowsher et al., 1988), which is then reduced to ammonium by the concerted action of NR (Foyer and Mullineaux, 1998). With the action of two enzymes, GOGAT and GS, ammonium can form the amino acids that combine with the organic acids produced by photosynthesis.

This study showed that five transcripts of putative NR and two genes encoding putative GS were downregulated under conditions of low-N stress because of N source limitations (Supplementary Table S10). However, the expression levels of NR and GS in Arabidopsis and rice were elevated under low-N conditions, possibly indicating differences in the mechanisms regulating N uptake and accumulation among P. tomentosa, Arabidopsis, and rice. A previous investigation revealed direct coupling of N assimilation and photosynthesis in chloroplasts, in a mechanism designated as nitrate photo assimilation (Bot et al., 2009). Several studies have revealed universal physiological and metabolic changes involved in the responses of plants to N deficiency, including significant reductions of chlorophyll content, growth, and photosynthesis, protein, starch, etc. (Mesnard and Ratcliffe, 2005). We detected four decreased transcripts of putative chlorophyll a–b binding proteins, indicating downregulated photosynthesis under low-N conditions. Similarly, previous investigations showed that genes involved in photosynthesis were downregulated in the leaves of Arabidopsis (Liu et al., 2017) and roots of rice (Cai et al., 2012).

Among plant hormones, auxin and ABA have been reported to be closely related to N signaling (Sakakibara, 2003; Lu et al., 2015; Yu et al., 2016). Our results demonstrated the interaction of GA signaling with N deprivation. Transcripts of putative GA-insensitive (GAI, a DELLA domain-containing GRAS family transcription factor) were induced (1.57-fold upregulation) and putative repressor of GA1-3-like 1 (RGL1) was suppressed (1.78-fold downregulation) under conditions of low-N stress. In cucumbers, suboptimal root zone temperatures were shown to suppress GA biogenesis and plant growth, while exogenous GA restored seedling biomass and enhanced N uptake (Bai et al., 2016). These results show the interdependent promotional role between GA biogenesis and N acquisition.

In Arabidopsis, N deprivation induced an increase in auxin content and expression of the auxin synthetic gene, TAR2, and the expression of auxin influx carriers (AUX/LAX family, such as AUX1, LAX1, LAX2, LAX3) and efflux transporter factors (PIN and PGP, e.g., PIN1, PIN2, and PIN4) could be regulated by the nitrogen/carbon (N/C) ratio (Gutiérrez et al., 2007). In addition, auxin signaling is related to N availability, and ARF8 could be induced by nitrate deficiency. miR160 has been shown to suppress ARF10, ARF16, and ARF17, while miR167 targets ARF6 and ARF8 (Rhoades et al., 2002). In maize, low N supply increased the levels of root auxin and NO, and then enhanced root elongation (Mi et al., 2008). Auxin was shown to trigger the accumulation of miR393, leading to the repression of TIR1/AFBs, depressed auxin perception and response, and ultimately, homeostasis (Chen et al., 2011). In addition, by binding to auxin and AUX/IAA, TIR1 could direct the ubiquitination and degradation of AUX/IAA proteins. On the other hand, binding to AUX/IAA promotes the activation of ARFs and the depression of other early genes in the auxin response pathway (Mockaitis and Estelle, 2008). Studies of Arabidopsis have revealed that the coordination of miRNA nodes directs the nitrate and auxin signaling involved in lateral root formation. Both miR393 and its target, AFB3, are known to be induced by nitrate and glutamine/glutamate, which affect the uptake of auxin (Vidal et al., 2010). However, miR167 is depressed under conditions of nitrate treatment, thus causing the accumulation of ARF8 and the downstream gene product, GH3, leading to a homeostatic level of auxin and modulation of the root architecture (Yang et al., 2006). Previous studies have validated the interaction of pto-miR160a and its target mRNA, pto-ARF16, which are involved in tree growth (Tian et al., 2016). Our previous investigation demonstrated the downregulation of pto-miR160f, pto-miR167a/b, pto-miR167c–f, pto-miR167g, and pto-miR390a–d and the upregulation of pto-miR393a/b under conditions of low-N stress in P. tomentosa, and thus the association of nitrate with auxin signaling (Ren et al., 2015). Furthermore, our transcriptomics data showed the corresponding induction of putative pto-IAA10 (CL5591.Contig2_All, 2.2-fold upregulation) and pto-ARF6 (Unigene21573_All, 1.5-fold upregulation), indicating that the nitrate-induced auxin signaling pathway plays a role in root architecture plasticity.

Previous studies provided evidence that the biosynthesis of CTK is dependent on N status to induce IPT genes, and that nitrate and glutamine act as inducers in Arabidopsis (Takei et al., 2004) and rice (Kamada-Nobusada et al., 2013). Similar to ABA, reduced CTK in a low-N environment was reported to enhance lateral root formation (Kiba et al., 2010). In addition, cytokinin-N-glucosyltransferase was shown to be involved in CTK homeostasis in Arabidopsis. Our study detected the downregulation of putative cytokinin-N-glucosyltransferase (1.77-fold downregulation), consistent with the results seen in Arabidopsis under conditions of N starvation (Bi et al., 2007) and rice (Cai et al., 2012). By catalyzing N-glucosylation of CTKs, CTKs then regulated CTK-response genes, such as CKX3, AHK2, AHK3, ARR1, and LOG2 (Wang et al., 2013), thus controlling the response to N deficiency in P. tomentosa. In addition, ethylene response regulator was increased in Botrytis cinerea-infected Solanum lycopersicum under nitrate-limiting conditions, suggesting that N status may be related to pathogen susceptibility (Vega et al., 2015). As a consequence, diverse plant hormones should be involved in the adaptation of P. tomentosa to conditions of low-N stress.

To investigate the effects of low-N stress on the metabolic regulation of P. tomentosa and validate gene expression at the level of the transcriptome, GC-MS was further applied for identification and robust quantification of primary metabolites in plant samples, including sugars, sugar alcohols, amino acids, organic acids, and polyamines. These results should enable a comprehensive evaluation of metabolic changes in P. tomentosa exposed to low-N stress, and highlight fields of interest for future studies of the NUE of poplar. The contents of amino acids were changed markedly under conditions of low N supply, including valine (1.32-fold downregulation), L-alanine (0.69-fold downregulation), and L-isoleucine (0.76-fold downregulation), similar to the gene expression pattern of the transcriptome. These changes were somewhat different from those seen in previous studies related to N starvation. In tomatoes, the L-alanine level was reported to decrease to less than half that recorded at the start of the experiment, while L-isoleucine and valine did not change markedly (Urbanczyk-Wochniak and Fernie, 2005). Tschoep et al. (2009) found that a low-N growth regime caused shoot growth reduction and increased the levels of many amino acids. The levels of amino acids, including alanine, serine, glycine, valine, leucine, isoleucine, phenylalanine, and tyrosine, derived from pyruvate and phosphoenolpyruvate increased in Synechocystis under conditions of N deficiency (Tschoep et al., 2009). This may have been related to differences between the species or experimental conditions. Urbanczyk-Wochniak and Fernie (2005) reported that different environmental factors could cause different changes in amino acids, although the variation of amino acids was not linear under low-N conditions, with some amino acids showing a marginal decline initially but an increased level at the end of the experiment. The results of this investigation were also somewhat consistent with previous studies of the poplar transcriptome, which indicated that poplar species could slow down N assimilation under conditions of limited N supply (Luo et al., 2013).

In this study, the downregulation of genes participating in the glyoxylate cycle [isocitrate lyase, malate synthase, and phosphoenolpyruvate carboxykinase (ATP: oxaloacetate carboxy-lyase (transphosphorylation), EC 4.1.1.49] may have had adverse effects on the growth of P. tomentosa, indicating a depressed state of C metabolism providing less carbon (by synthesizing a variety of carbohydrates) and energy (providing NADPH and then ATP) for plant growth and development. On the other hand, N and C metabolism could be mutually affected due to the requirement for C skeletons for the production of N metabolites, such as amino acids. In dicots, nitrate levels could affect carbohydrate metabolism (Scheible et al., 1997). Phloem transport and C export are suppressed under low-N conditions (Nunes-Nesi et al., 2010). The significantly altered amino acid production seen in this study was related to pyruvate recycling. The accumulation of foliar starch was also found in Arabidopsis and maize exposed to low-N conditions due to the reduced demand of C skeletons for N compounds, including amino acids and proteins (Ikram et al., 2012; Schlüter et al., 2012). Concerning C metabolism, the expression of some sugars, i.e., D-fructose (.433-fold), D-galactose (.171-fold), and D-glucose (.229-fold), was induced under conditions of N deficiency in this investigation. Similarly, it was reported that glucose and fructose accumulated under low-N conditions in Synechocystis (Asayama et al., 2004). In contrast, fructose and galactose did not vary markedly under high-light conditions in tomatoes, whereas they accumulated under low-light conditions (Urbanczyk-Wochniak and Fernie, 2005). Glucose and fructose are signaling molecules and energy sources in plants under conditions of abiotic stress (Bonneau et al., 2013). Although foliar starch, glucose, and fructose accumulate under low-N conditions, N deficiency inhibits photosynthetic capacity and growth. Poplar roots actively forage for nutrients under low-N conditions, while the acquired N appears insufficient for the biosynthesis of photosynthetic enzymes and metabolic precursors, leading to decreased photosynthetic capacity (Sardans and Peñuelas, 2012). In this study, the level of phytol related to chlorophyll synthesis increased significantly, which would directly affect the chlorophyll content in photosynthesis. Moreover, there was a marked accumulation of oxalic acid, which is involved in in the photosynthetic C3 cycle and photorespiration. Accumulation of phosphoric acid was also observed; this is involved in cyclic photophosphorylation. Interestingly, in contrast to most previous studies, four amines (putrescine, ethanolamine, ethylamine, and cadaverine) showed significant changes in this study. As components of N assimilation, the synthesis of ethylamine, ethanolamine, and putrescine decreased. Consistently, previous studies have shown that polyamines play an important role in modulating the sensitivity/tolerance to N stress (Syed et al., 2011). In addition to carbohydrate metabolic processes, N deficiency also affects the metabolism of other nutrients, such as phosphate and sulfur transport, indicating alteration of the N:C/P/S ratios and potential interactions (Supplementary Table S7). Specifically, several sugar-related pathways are known to interconnect with N-responsive pathways, and C/N balance has particular significance for plant growth. Obviously, crosstalk between N and C/P/S metabolism is involved in nutrient utilization and signal transduction, and a limitation of one element could cause imbalances in the nutrient network and affect the uptake, accumulation, and utilization of other nutrients (Schachtman and Shin, 2007; Rennenberg et al., 2010).