Metabolomics and Transcriptomics Integration of Early Response of Populus tomentosa to Reduced Nitrogen Availability

Nitrogen (N) is one of the most crucial elements for plant growth and development. However, little is known about the metabolic regulation of trees under conditions of N deficiency. In this investigation, gas chromatography-mass spectrometry (GC-MS) was used to determine global changes in metabolites and regulatory pathways in Populus tomentosa. Thirty metabolites were found to be changed significantly under conditions of low-N stress. N deficiency resulted in increased levels of carbohydrates and decreases in amino acids and some alcohols, as well as some secondary metabolites. Furthermore, an RNA-sequencing (RNA-Seq) analysis was performed to characterize the transcriptomic profiles, and 1,662 differentially expressed genes were identified in P. tomentosa. Intriguingly, four pathways related to carbohydrate metabolism were enriched. Genes involved in the gibberellic acid and indole-3-acetic acid pathways were found to be responsive to low-N stress, and the contents of hormones were then validated by high-performance liquid chromatography/electrospray ionization tandem mass spectrometry (HPLC-ESI-MS/MS). Coordinated metabolomics and transcriptomics analysis revealed a pattern of co-expression of five pairs of metabolites and unigenes. Overall, our investigation showed that metabolism directly related to N deficiency was depressed, while some components of energy metabolism were increased. These observations provided insights into the metabolic and molecular mechanisms underlying the interactions of N and carbon in poplar.


INTRODUCTION
Plant growth is perturbed by various biotic and abiotic stresses (Rejeb et al., 2014). Among abiotic stresses, nitrogen (N) stress has a major effect on plant physiological activity. Many biological molecules, including nucleic acids, amino acids, proteins, chlorophyll, lipids, and a variety of other metabolites, contain N, which is required for their synthesis (Kusano et al., 2011). N is actively transported or taken up by the plant root system. Organisms utilize N from three sources. The first source is N 2 from the air, which can be assimilated by leguminous rhizobia (Zahran, 1999), and the second is organic N in the soil, which can be taken up by plants in specific environments (Jones et al., 2005). However, the main sources of N, which are taken up by most higher plants via transporters, are ammonium (NH 4 + ) and nitrate (NO 3 − ) from the soil (Jackson et al., 2008). The absorption and utilization of N (NH 4 + and NO 3 − ) are highly regulated in plants (Patterson et al., 2010). A wide range of physiological activities of plants is disrupted by N deficiency, including photosynthesis, signal transduction, and the synthesis of phospholipids, endogenous hormones, and many secondary metabolites (Shao et al., 2020;de Bang et al., 2021;Mu and Chen, 2021).
Nitrogen (N) deficiency usually reduces amino acid synthesis (Albinsky et al., 2010), while metabolite profiling analyses of Arabidopsis and maize have shown that growth under low-N conditions causes increases in the levels of many amino acids (North et al., 2009;Broyart et al., 2010;Trachsel et al., 2013). This may be due to varying experimental conditions. For example, gamma-aminobutyric acid (GABA) was unaltered under conditions of nitrate deficiency in high-light conditions but was induced in low-light conditions in tomatoes (Urbanczyk-Wochniak and Fernie, 2005). Furthermore, N deficiency was shown to affect the biosynthesis of some carbohydrates, as the fundamental processes of carbon (C) and N metabolism are tightly coordinated (Cao et al., 2019). The C skeleton and energy provided by carbohydrates are required for photosynthesis and N uptake (Gutieìrrez et al., 2005;Zheng, 2009). It has been demonstrated that N deficiency suppresses the levels of carbohydrate and major soluble sugars, but stimulates the accumulation of starch (Scheible et al., 1997). A wide range of genes is involved in the low-N stress response in plants (Shi et al., 2016). Some N assimilation process-related genes, especially those involved in the ornithine urea cycle (OUC) and tricarboxylic acid (TCA) cycle, were identified in Aureococcus anophagefferens by transcriptomics analysis . High levels of expression of genes involved in longchain fatty acid and hydrocarbon biosynthesis were also found in Botryococcus braunii (Chlorophyta) under conditions of N deprivation (Fang et al., 2015). Moreover, transcription factors with driving roles in N × C interactions were shown to be associated with N stress in maize .
It is well known that forest plantations of poplar species have large effects on C mitigation, the pulp industry, and biomass production (Studer et al., 2011). In recent years, some fastgrowing tree species, such as Populus spp., have been widely planted worldwide (Rennenberg et al., 2010). Unlike annual plants, the yields of which are highly dependent on the addition of high-N fertilizer, perennial plants, such as trees, achieve their yields with minimal N input (20-50% less) because of their remobilization of resources, such as N in bark storage proteins (Karp and Shield, 2008). This difference exerts a large influence on life cycle studies of bioenergy chains, considering the energy consumption required to produce N fertilizer. Elucidating the genetic regulation underlying N use efficiency (NUE) , and identifying the important genes via genomic and other "-omics" approaches (Tuskan et al., 2018), will facilitate progress in genetically modified tree breeding for sustainable and efficient supply of biomass plants in the future. These investigations help maintain environmental and financial sustainability (Karp and Shield, 2008). Considering the increasingly recognized importance of forestry in ecological balance and the accelerating exhaustion of mineral resources for fertilizer, it is necessary to investigate the regulatory mechanisms and genes involved in the NUE of trees for tree improvement.
Populus tomentosa, also known as Chinese white poplar, is one of the fastest-growing poplar species; it is widely distributed in northern China and is of great economic and ecological importance (Du et al., 2012). P. tomentosa is considered a model system to explore and understand the morphological changes and molecular mechanisms of tree growth and development, as well as responses to the environment. Previously, our laboratory focused on systematically characterizing the molecular responses of P. tomentosa under conditions of N deficiency, including the global profiling of microRNAs (miRNAs) (Ren et al., 2015), the degradome (Chen M. et al., 2015), and long non-coding RNAs (lncRNAs) . It was reported that total C content, reactive oxygen species (ROS), ATP, peroxidase, superoxide dismutase (SOD), and glutamine synthetase (GS) were increased in two contrasting poplar clones Nanlin 1388 and Nanlin 895 . In addition, the results of transcriptomics analyses of N signaling, metabolism, and storage in poplar shoot growth and development have been reported. Transcriptomics studies showed that N starvation suppressed the expression of genes encoding most nitrate transporters (NRTs) and ammonium transporters (AMTs) in poplar leaves and genes involved in N assimilation in both roots and leaves (Luo et al., 2013). N starvation treatment was also shown to increase the fine root length and surface area, foliar starch concentration, and transcript abundance of several AMTs (AMT1;2) and NRTs (NRT1;2 and NRT2;4B) in the roots of slow-growing species (P. popularis) and a fastgrowing species (P. alba × P. glandulosa) during acclimation to limiting N supply (Luo et al., 2013). Global transcriptomic reprogramming was shown to play a critical role underlying the physiological and morphological response of poplar leaves and roots to N starvation and excess (Luo et al., 2015). Similarly, global transcriptome reprogramming and activation of root growth were also revealed in poplar (Populus tremula × Populus alba) to low-N supply (Wei et al., 2013a). PtaNAC1-centered subnetwork was further revealed to be involved in increasing root biomass, which was helpful in the dynamic adjustment of poplar root architecture to low-N availability (Wei et al., 2013a,b). Further investigation demonstrated that PtaNAC1 and PtaRAP2.11 encoding transcription factors, F-box proteinencoding gene similar to Hawaiian Skirt (PtaHWS) had a markable influence on root development of poplar under low N (Dash et al., 2015).
Most studies to date have focused on morphological, physiological, and transcriptional changes in poplar, and few have examined the global metabolic changes in poplar combined with transcriptomics profiles under low-N conditions. To identify the metabolites and the corresponding regulatory pathways and genes involved in low-N signaling in trees, we identified genes and metabolites produced responsive to low-N stress through metabolomics and transcriptomics profile analyses. Intriguingly, we detected alterations in metabolites and transcriptional reprogramming, providing insights into the physiological and metabolic changes involved in growth and development, and obtained information to improve NUE in plantations in both agriculture and forestry.

Plant Materials and Treatment
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 4 NO 3 (Sigma-Aldrich) and 1 mM KNO 3 (Sigma-Aldrich) as sufficient N conditions (KK) (control) or with.01 mM NH 4 NO 3 and 1 mM KCl (Sigma-Aldrich) instead of KNO 3 for low-N treatment (DN). Whole P. tomentosa plants were harvested in the midmorning, immediately frozen in liquid N, and stored at −80 • C.

Metabolite Extraction
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 2 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.

Gas Chromatography-Mass Spectrometry and Metabolite Profile Analysis
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/massto-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 1 and Kyoto Encyclopedia of Genes and Genomes (KEGG) 2 database.

Metabolomics Data Analysis and Metabolic Pathway Construction
Before the data analysis, all data were standardized for meancentering and unit-variance scaling using SIMCA-P with the default parameters (ver. 11.5 3 ). 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 4 and metabolic pathway networks were constructed using Cytoscape 5 .

Total RNA Extraction and Illumina Sequencing Analysis
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 . 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 lowquality 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 6 , 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) 7 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 8 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 9 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.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analysis
Unigenes were aligned with Kyoto Encyclopedia of Genes and Genomes (KEGG) 2 to predict the metabolic pathways, with p ≤ 0.05 and q ≤ 0.05 taken to indicate significant enrichment Wang X. et al., 2016), and the functions of gene products (using BLASTx 10 ) 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 11 , with p ≤ 0.05 taken to indicate significance Wang X. et al., 2016).

Hormone Quantification by High-Performance Liquid Chromatography/Electrospray Ionization Tandem Mass Spectrometry
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 2 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).

Real-Time Quantitative Reverse Transcription PCR to Detect the Transcripts of Differentially Expressed Genes in Response to Low-N Treatment (LN)
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 TM (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.

Metabolite Profiling of Populus tomentosa Under N Deficiency
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, Lisoleucine, 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.

Network Construction of Responsive Metabolites in Populus tomentosa Under N Deficiency
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, Dfructose, 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.

RNA-Seq of Populus tomentosa Under N Deficiency and Gene Function Annotations
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).

Global Transcriptomic Changes of Populus tomentosa Under N Deficiency
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).

Genes Related to N Acquisition, Allocation, and Assimilation in Response to Low-N Stress
Plants take up N from two sources in the soil, ammonium, and nitrate. Nitrate is reduced to nitrite (NO 2 − ) 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.

Interaction of C and N Metabolism Under Conditions of N Deficiency
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.21fold downregulated), photosystem I reaction center subunit psaK (1.37-fold downregulated), oxygen-evolving enhancer protein   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, fructosebisphosphate 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 longchain-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.

Hormone Signaling-Related Gene Profiles and Hormone Quantification
Auxin and ABA were reported to be closely related to N signaling (Sakakibara, 2003;Mockaitis and Estelle, 2008;. 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.23fold 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.

Detection of Dynamic Expression of Differentially Expressed Genes by Real-Time Quantitative Reverse Transcription PCR
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

Correlation of Gene Expression and Metabolic Changes Under 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-Dglucose 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.

DISCUSSION
Due to their ability to perform high-throughput genome-wide characterization, transcriptomics analyses have been extensively applied to systemically investigate the genes responsible for the development and adaptation of plants to different environmental stresses (Vogel et al., 2016;Wang Y. N. et al., 2016;Woo et al., 2016). It is well known that N, one of the most important nutrients, functions as a component of various key cell molecules, such as proteins (amino acids), nucleic acids, chlorophyll, and secondary metabolites. Consequently, many studies have focused on the mechanisms underlying the roles of N in plants, and DEGs in response to N deficiency. Several genes involved in N assimilation and utilization have been identified in a variety  of plants, including cucumber (Zhao et al., 2015), rice , maize , Medicago truncatula (Bonneau et al., 2013), and Arabidopsis (Günther et al., 2012).

Global Changes in the Transcriptomic Network of Populus tomentosa Under N Deficiency
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).

Responses of Hormone and Hormone-Related Genes to N Deficiency
Among plant hormones, auxin and ABA have been reported to be closely related to N signaling (Sakakibara, 2003;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-miR167cf, 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-Nglucosyltransferase (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.

Global Changes of Metabolic Networks in Response to N Deficiency
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), Lalanine (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 Dglucose (.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).

CONCLUSION
In summary, this study revealed the global changes of transcripts and metabolites occurring in response to N deprivation in P. tomentosa. Specifically, 2,662 DEGs, 30 significantly changed metabolites, and three altered hormones were detected in our study. Combining transcriptomic and metabolic profiles, our investigation revealed a general depression of molecular and physiological metabolism under low N stress, including N absorption and assimilation, photosynthesis, glycolysis, and the TCA cycle. This extends our understanding of the relationships between responsive genes and downstream metabolic responses and provides a basis for the integrated and comprehensive analysis of molecular responses to N deficiency, which may improve understanding of the molecular and metabolic mechanisms underlying the phenotype of P. tomentosa, and the interaction of C and N under conditions of N deficiency.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: National Center for Biotechnology Information (NCBI) BioProject database under accession number PRJNA296440.

AUTHOR CONTRIBUTIONS
YW designed the experiment and edited the manuscript. MC and YY performed the research and wrote the manuscript. XY, LZ, TF, and XH contributed to the analytical tools and reagents. All the authors read and approved the final manuscript.