Comparative Transcriptome Analysis of Sophora japonica (L.) Roots Reveals Key Pathways and Genes in Response to PEG-Induced Drought Stress under Different Nitrogen Conditions

: Sophora japonica is a native leguminous tree species in China. The high stress tolerance contributes to its long lifespan of thousands of years. The lack of genomic resources greatly limits genetic studies on the stress responses of S. japonica . In this study, RNA-seq was conducted for S. japonica roots grown under short-term 20% polyethylene glycol (PEG) 6000-induced drought stress under normal N and N starvation conditions (1 and 0 mM NH 4 NO 3 , respectively). In each of the libraries, we generated more than 25 million clean reads, which were then de novo assembled to 46,852 unigenes with an average length of 1310.49 bp. In the differential expression analyses, more differentially expressed genes (DEGs) were found under drought with N starvation than under single stresses. The number of transcripts identiﬁed under N starvation and drought in S. japonica was nearly the same, but more upregulated genes were induced by drought, while more downregulated genes were induced by N starvation. Genes involved in “phenylpropanoid biosynthesis” and “biosynthesis of amino acids” pathways were upregulated according to KEGG enrichment analyses, irrespective of the stress treatments. Additionally, upregulated N metabolism genes were enriched upon drought, and downregulated photosynthesis genes were enriched under N starvation. We found 4,372 and 5,430 drought-responsive DEGs under normal N and N starvation conditions, respectively. N starvation may aggravate drought by downregulating transcripts in the “carbon metabolism”, “ribosome”, “arginine biosynthesis pathway”, “oxidative phosphorylation” and “aminoacyl-tRNA biosynthesis” pathways. We identiﬁed 78 genes related to N uptake and assimilation, 38 of which exhibited differential expression under stress. A total of 395 DEGs were categorized as transcription factors, of which AR2/ERF-ERF, WRKY, NAC, MYB, bHLH, C3H and C2C2-Dof families played key roles in drought and N starvation stresses. The transcriptome data obtained, and the genes identiﬁed facilitate our understanding of the mechanisms of S. japonica responses to drought and N starvation stresses and provide a molecular foundation for understanding the mechanisms of its long lifespan for breeding resistant varieties for greening. partitioning, responses amino acids


Introduction
Sophora japonica (Chinese scholar tree) is a native deciduous tree species in China belonging to the Fabaceae family. Because of its well-developed root system, large crown and dense foliage, S. japonica is widely used for urban and rural greening in northern China. S. japonica has a long lifespan and many individuals in China have lived for thousands of years, which may be attributed to its strong adaptability and resistance to various abiotic in the amino acid metabolism, transport and stress categories were upregulated, whereas genes in hormone metabolism and redox pathways were downregulated in the roots of Arabidopsis thaliana under N starvation [32]. Lu et al. [33] studied the mitigating effect of N supply on drought and found that genes related to the antioxidative system and secondary metabolism were promoted by N addition in poplar. The mechanisms underlying how plants respond to single drought or N deficiency stresses have been widely studied, whereas the effect of N starvation on drought stress has rarely been studied. In recent years, transcriptomics has become an effective strategy to study changes in gene expression or metabolic processes when an organism is under stress, especially in non-modal plants [31,34,35]. There are only two genetic studies on S. japonica, which have only used the aboveground plant parts as sequencing materials [36,37]. Thus, the mechanisms of how S. japonica roots respond to drought and N deficiency at the molecular level remain unknown.
In this study, S. japonica seedlings were exposed to short-term treatments with N starvation and 20% polyethylene glycol(PEG)-6000-induced drought stresses. The leaf chlorophyll content and the MDA and proline contents in both leaves and roots were measured. Simultaneously, twelve fine root samples were sequenced using an Illumina HiSeq 2000 platform. The main objectives of this study were as follows: (1) to understand the molecular regulation of drought and N starvation stresses in S. japonica; (2) to preliminarily explore the effect of N starvation on drought stress; and (3) to identify candidate genes participating in N uptake and metabolism. The S. japonica transcriptome profile and metabolic pathway analysis under N deficiency and drought stress can provide the molecular basis for gene expression, genomics, and functional studies on the abiotic stress response of S. japonica, as well as other leguminous tree species.

Plant Materials and Treatments
Seeds of Chinese scholar trees were purchased from Shaanxi Sanyuan Gold Seed Industry Technology Co., Ltd., Yangling Branch. After the seeds germinated, they were then planted in pots (10 L) filled with fine sand. Seedlings were cultivated in a greenhouse (natural light, day/night 28/20 • C, 80% relative humidity), and irrigated every three days with 50 mL of modified Hoagland nutrient solution. After 8 weeks, the sand attached to the seedling roots was washed away with tap water, and uniform seedlings were transferred and acclimated to hydroponic culture for 4 weeks. The nutrient solution used for hydroponics and irrigation was reported in a previous study [18]. The hydroponic solution was renewed every 3 days and aerated continuously by air pumps.
To determine the appropriate time and PEG-6000 concentration for drought treatment, a preliminary experiment was conducted. The plants described above were treated with four PEG concentrations (0%, 10%, 20%, and 30%) for 24 h. Roots were harvested across four time points (0, 6, 12 and 24 h) for POD and SOD activity assays. As the results showed ( Figure S1), compared with the control, both activity levels changed rapidly at 0-6 h under all three concentrations of PEG-6000. Three hours was selected as the experimental time for the drought treatment, considering that changes at the transcriptomic level are much more rapid than those at the physiological level, as suggested by a previous study in Arabidopsis [38]. Additionally, the two enzyme activity levels changed relatively slowly under the 10% PEG treatment. The 30% PEG concentration was so high that the cell osmotic balance may have been disrupted at an earlier time point, as SOD activity under 30% PEG decreased continuously at 6-24 h. Therefor, 20% PEG was used to induce moderate drought stress in this study.
In addition to the preliminary experiment, hydroponic-adapted plants were also used for the formal experiment. Briefly, we designed four treatments in this study as described in Figure 1: (a) NN (1 mM NH 4 NO 3 applied for 15 h); (b) NL (1 mM NH 4 NO 3 applied for 12 h, then 20% PEG-6000 applied for 3 h); (c) LN (0 mM NH 4 NO 3 applied for 15 h); (d) LL (0 mM NH 4 NO 3 applied for 12 h, then 20% PEG-6000 applied for 3 h). Leaves and fine roots of the four treatment groups were harvested simultaneously but kept separate.
Three independent biological replicates, each containing a pool of three different plants, were sampled. Roots were used for high-throughput transcriptomic profiling and qRT-PCR verification. Both roots and leaves were used for physiological index measurements. The samples were rapidly frozen in liquid nitrogen, and then stored at −80 • C.

Physiological Index Measurements
POD and SOD activities were measured using physiological assay kits (Nanjing Jiancheng Bioengineering Institute, Nanjing, China) according to the manufacturers' recommendations. The MDA content was measured using the thiobarbituric acid (TBA) method [39]. The proline content was determined according to the method described by Chołuj et al. [40]. We analyzed the variables using one-way ANOVA in SPSS 23 (SPSS Inc., Chicago, IL, USA) by Tukey's test at the p < 0.05 probability level.

Transcriptome Sequencing
Total RNA was extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) with RIN > 7.5 in all samples (Table S1). cDNA libraries were generated using the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) according to the manufacturer's protocols and were sequenced by Biomarker Technologies Co., Ltd. (Beijing, China) on an Illumina HiSeq 2000 platform. Reads with high quality were necessary for the downstream analyses. Clean data were obtained by removing reads containing adapters, reads containing poly-N and low-quality reads from the raw data. The Q30 (the proportion of bases with a quality value greater than 30) and GC-content (the percentage of nitrogenous bases in a DNA or RNA molecule that are either guanine or cytosine) of the clean data were calculated at the same time. Then, transcriptome assembly was accomplished using Trinity software [41] with min_kmer_cov set to 2 by default. The assembled sequences were designated as "unigenes" for further annotation.

Gene Functional Annotation
For functional annotation, all unigenes were searched against public databases including the NR (NCBI nonredundant protein sequences), KOG/COG/eggNOG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KEGG (Kyoto Encyclopedia of Genes and Genomes) and GO (Gene Ontology) databases by using Blast software [42] with an E-value ≤ 10 −5 . Pfam (Protein family) annotation was performed using HMMER software [43] with an E-value less than 10 −10 . The plant TF database (http://planttfdb.cbi.pku.edu.cn/ accessed on 25 Mar 2020) was used to search for putative transcription factors in S. japonica.

Identification of Differentially Expressed Genes (DEGs)
The generated clean data were mapped back to the assembled transcriptome. Gene expression levels were estimated by RSEM [44] and normalized to FPKM (fragments per kilobase per million mapped reads) values. DEGs were identified using the DESeq package [45] with a false discovery rate (FDR) < 0.05 and fold change ≥2 used as thresholds (positive or negative for either over-or underexpression, respectively). In the downstream DEG analyses, GO and KEGG databases were used to classify DEGs, and KEGG enrichment was conducted by KOBAS software [46] to test the statistical enrichment of the DEGs. A heat map showing the expression of selected DEGs was drawn using TBtools [47]. Phylogenetic trees of the identified NRT and AMT proteins with associated proteins in Arabidopsis were constructed via the neighbor-joining (NL, 1000 bootstrap replicates) method with MEGA 7.0 software (Center for Evolutionary Medicine and Informatics, Tempe, AZ, USA). The protein sequences are provided in Table S2.

qRT-PCR Validation
To validate the reliability of the transcriptome sequencing results, we randomly selected 15 genes and subjected them to quantitative real-time PCR. RNA samples for transcriptome sequencing were also used for RT-qPCR. First-strand cDNA synthesis was performed using EasyScript ® One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China) according to the manufacturer's instructions using an anchored oligo (dT) 18 primer and 500 ng of total RNA. Reactions were performed on a LightCycler ® 96 real-time PCR system (Roche, Mannheim, Germany) using 2 × SYBR real-time PCR premixture (Bioteke, Beijing, China) in a 20 µL reaction under the following conditions: one cycle of 150 s at 95 • C, followed by 45 cycles at 95 • C for 10 s, 58 • C for 20 s and 72 • C for 30 s. The primers used are listed in Table S3. Helicase was used as the reference gene. For each of the selected genes, three biological replicates with three technical replicates were assayed. The relative expression levels were calculated using the 2 −∆∆CT method [48].

Physiological Characteristics Affected by Drought and N Starvation
After 15 h of exposure to N starvation and 3 h of exposure to 20% PEG, the leaf chlorophyll content and root MDA content did not change under all treatments (Figure 2A,C). The MDA content in leaves increased significantly only under drought with N starvation (LL) when compared with that of the control (NN) ( Figure 2B). The proline content in leaves increased significantly only under drought ( Figure 2D), and the proline content in roots increased upon drought but decreased under N starvation, whereas it was not affected by drought with N starvation ( Figure 2E).

De Novo Transcriptome Assembly
Transcriptome analysis using the Illumina HiSeq 2000 sequencing platform was performed. We obtained a total of 109.04 Gb clean data. For each cDNA library, 25.72-34.61 million clean (high-quality) reads were generated (Table S4). The Q30 quality scores were over 92.11%, and the GC content was approximately 45% for all samples. The transcriptome was then de novo assembled to 86,672 transcripts (total length of 1.25 × 10 8 bp), which were grouped into 46,852 trinity 'unigenes' with an average length = 1310.49 bp (Table S5). Approximately 48.1% of the unigenes had a length greater than 1,000 bp ( Figure S2).
Sequence annotation against public databases resulted in the successful annotation of 37,183 (79.36%) unigenes in at least one of the databases (Table S6 and S7). In brief, a total of 35,524 unigenes were annotated against the NR database. Additionally, 24,070 transcripts were assigned to 5,786 GO terms and categorized into 50 functional groups ( Figure S3; Table S8), and 18,607 putative proteins were assigned to five main categories that included 131 KEGG pathways (Table S9).

Differentially Expressed Gene Analyses
Principal component analysis of the gene expression results showed that the main difference between the samples was derived from drought stress, which explained 52.7% of the total detected variation. N status (0 and 1 mM) contributed to another 27.6% of the variability ( Figure 3A). The three replicates of each sampling group clustered well together, with the exception of sample LL-1, which clustered with the NL samples. In addition, LL-1 had a low correlation with LL-2 and LL-3 (Table S10). After excluding the outlier LL-1, the correlation among biological replicates was higher than that across samples. A total of 11 samples were used for downstream DEG analyses.

DEGs Induced by Drought and N Starvation
Under drought and N starvation, we identified 11,808 DEGs from fine roots of S. japonica in five pairwise comparisons. When compared with the control (NN), more DEGs were found under drought with N starvation than under a single stress, with 7,389 genes being differentially expressed. The number of transcripts identified between N starvation and drought in S. japonica varied little, but more upregulated genes were induced by drought, while more downregulated genes were induced by N starvation ( Figure 3B). Among all these DEGs, 520 were commonly upregulated by drought, N starvation and drought with N starvation, whereas 390 transcripts were commonly downregulated ( Figure 3C). To understand the functions of these genes, the GO database, and KEGG pathway enrichment were used.
When the DEGs were mapped to terms in the GO database and compared to the whole transcriptome background, they were categorized into small functional groups in three main categories (biological process, cellular component, and molecular function) of the GO classification. Under all three stresses, metabolic process, cellular process and singleorganism process were dominant within the biological process category. The cell, cell part and membrane were the most frequent groups in the cellular component category. Catalytic activity and binding were major parts of the molecular function category (Table S11).
To further understand the biochemical pathways of DEGs, KEGG analysis was conducted. The numbers of KEGG pathways were individually determined, and the significant pathways were further filtered using a cut-off p-value of 0.05. Genes in the "biosynthesis of amino acids", "phenylpropanoid biosynthesis" and "plant-pathogen interaction" pathways were upregulated, irrespective of the stress treatments, while "nitrogen metabolism" pathway was enriched under drought and drought with N starvation ( Figure 4A,C,E). Remarkably, genes in the "plant hormone signal transduction" and "starch and sucrose metabolism" pathways were downregulated upon drought ( Figure 4B), while "photosynthesis" pathway was significantly enriched in downregulated DEGs under N starvation and drought with N starvation (Figure 4D,F). In addition, the number of DEGs in the "ribosome" pathway was much larger under both N starvation and drought with N starvation than under drought ( Figure S4).

DEG Analyses about the Effect of N Starvation on Drought Stress
To explore the effect of N starvation on drought stress, DEGs between two pairwise comparisons were analyzed, NN vs. NL and LN vs. LL. Upon drought, 2,512 genes were coexpressed (1,786 upregulated and 718 downregulated) under the two N statuses ( Figure 3D). With a similar number of upregulated genes, approximately 1000 more DEGs were downregulated upon drought under N starvation than under normal N condition.
The GO enrichment analysis of the two pairwise comparisons did not reveal differences among treatments (Table S11). Then, the identified DEGs were mapped to pathways of the KEGG database to identify significantly enriched biochemical pathways. Interestingly, the results indicated that genes involved in "starch and sucrose metabolism" were only downregulated under normal N condition ( Figure 4B), whereas those related to "carbon metabolism", "alanine, aspartate and glutamate metabolism" and "RNA transport" were only downregulated under N starvation condition ( Figure 5B). Furthermore, the "plant hormone signal transduction" pathway was significantly enriched in both upand downregulated DEGs under normal N condition but not enriched under N starvation ( Figure 4A,B and Figure 5A,B). Additionally, the KEGG map indicated that the pathways of "ribosome" (Figure S5), "arginine biosynthesis pathway" (Figure S6), "oxidative phosphorylation" ( Figure S7) and "aminoacyl-tRNA biosynthesis" ( Figure S8) were almost downregulated under N deficiency but upregulated under normal N condition.

Candidate Genes Related to N Uptake and Metabolism
According to the NR annotation, 78 unigenes were involved in N uptake and metabolism. Among them, 38 were differentially expressed in at least one comparison ( Figure 6; Table S12). In response to drought (NN vs. NL), most predicted NRTs were downregulated (except for c1253196.graph_c0 and c1254472.graph_c0), while two AMTs and all DEGs involved in N assimilation were upregulated. Similarly, most of the DEGs subjected to drought combined with N starvation (NN vs. LL) had consistent expression changes with those under drought, and the differences included a lack of an increase in NRT transcripts and the downregulation of one GS. Conversely, under N starvation (NN vs. LN), all three differentially expressed NRTs had increased transcript abundances, and the expression of only one NR was altered but was downregulated. However, other differentially expressed GOGATs and GDHs were all upregulated, in keeping with under drought treatment and drought with N starvation. The expression levels of DEGs predicted to participate in N transport and metabolism varied in response to drought under different N statuses. When N was deficient, the majority of NRTs were also downregulated upon drought, but 11 of 15 NRTs had a different trend compared with those under normal N condition. The expression of all NRs and NiR was enhanced regardless of N concentration, whereas two GS genes and two GOGATs were downregulated and all GDHs were expressed normally under N deficiency condition.

Transcription Factors
Transcription factors play important roles in upstream regulation when plants suffer from abiotic and biotic stresses. In the present study, we identified 1,051 TFs, which were categorized into 60 TF families. C2H2 was the most abundant family, followed by the C3H and bZIP families. Among them, 395 TFs belonging to 52 families showed different expression in at least one comparison (Figure 7; Table S13). The number of differentially expressed transcription factors (DETs) was higher under drought with N starvation than under either drought or N starvation alone. Most of the AP2/ERF-ERF family DETs were upregulated, whereas bHLH family DETs were downregulated by drought and drought with N starvation. Under N starvation, the enriched upregulated TF groups were members of the WRKY and NAC families, while C2C2-Dof was the largest downregulated family upon N starvation. Furthermore, drought stress increased the expression of genes involved in the AP2/ERF-ERF, WRKY and NAC families but reduced the transcription of genes in the bHLH family, regardless of N status.

Validation of RNA-Seq by qRT-PCR Validation
To verify the RNA-seq data, we validated the expression of 15 genes by qRT-qPCR using the same samples as for RNA-seq. Although the magnitude of the change was not completely equal, the qRT-PCR results of the tested DEGs were generally consistent with the RNA-sequencing data (Figure 8), which means our transcriptome is reliable for further studies.

Discussion
Drought and N deficiency negatively affect plant growth, and many physiological and biochemical parameters change to maintain the growth of plants under adverse conditions. In this study, POD and SOD activities rapidly changed in response to drought stress at 0-6 h in S. japonica ( Figure S1). Increased POD activity is beneficial to reduce oxidative damage in plant cells. SOD activity decreased first under drought at 6-12 h due to the sudden increase in reactive oxygen species but then increased at 6-12 h, which may be because drought induced SOD resynthesis to enhance drought resistance. The same trend was also observed in Camellia sinensis leaves under cold stress [49]. The chlorophyll content was negatively affected by all the stress treatments, although the change was not significant (Figure 2A), which may have been the result of the short stress duration. Interestingly, the MDA content in leaves noticeably increased under drought combined with N starvation but did not change in roots ( Figure 2B,C), which indicated that the cell membrane in leaves may be more sensitive than that in roots of S. japonica. Earlier leaf aging than root aging may be a survival strategy for perennials under adverse conditions [50]. Furthermore, the proline content in roots increased under drought but decreased under N starvation. Notably, the increases were even greater under adequate N condition than under N starvation in both leaves and roots ( Figure 2D,E). N is an essential element in proline synthesis. This result may suggest that N deficiency could aggravate drought stress to some extent. However, more efforts are required to explain the underlying mechanisms.
With the rapid development of high-throughput sequencing technology in recent years, transcriptomics provides opportunities for large-scale gene discovery and gene expression studies in nonmodel species without genomic resources. In this regard, we used an RNA-seq approach to explore the molecular mechanisms of the responses to drought and N deficiency stresses in the fine roots of S. japonica. A total of 46,852 unigenes were generated, with a higher N50 value of 1310 bp than two previous sequencing studies of S. japonica [36,37]. Higher quality transcriptome data can enrich the genetic information of S. japonica and provide a good foundation for the study of gene expression in this species.
Plants often activate rapid responses to biotic and abiotic stresses. For example, more than 16% of the genome exhibited altered expression levels in response to drought stress in cotton [51]. Similarly, we observed that 11,809 genes (approximately 25% of all the unigenes) were differentially expressed in fine roots of S. japonica under drought stress and N starvation. In addition, DEGs in different comparisons were categorized into large and varied metabolic pathways ( Figure S4). The quick and strong response indicated that S. japonica has a complex regulatory process in response to abiotic stress. Plants are not exposed to only single stresses in nature, as adverse conditions often occur simultaneously. The changes in genes were greater in S. japonica when drought and N deficiency were applied together than when they were applied alone ( Figure 3B,C). Whether seedlings were threatened by drought or N starvation, genes related to the "biosynthesis of amino acids" and "phenylpropanoid biosynthesis" were enriched in upregulated pathways in this study. The activation of the phenylpropanoid biosynthesis pathway contributes to the accumulation of various phenolic compounds that have the potential to scavenge harmful reactive oxygen species [52]. Amino acids are not only building blocks of different cellular components, but can also act as compatible osmolytes, regulate pH or act as a nitrogen or carbon reserve when plants are under stress [53]. The two pathways may play key roles in counteracting the adverse effects of stresses in S. japonica. It is worth noting that N starvation decreased the expression of genes in the photosynthesis pathway ( Figure 4D), consistent with previous study results, since N is an essential component of chlorophyll synthesis [54]. In addition, downregulated genes in the ribosome pathway accounted for the largest number of DEGs under both N starvation and drought with N starvation ( Figure S8B,C), indicating ribosome dysfunction caused by N deficiency, which would further hinder protein synthesis and may result in plant growth being restricted.
Many researchers have focused on the impact of nitrogen concentrations on the effects caused by drought. In S. japonica, more genes were upregulated than downregulated upon drought under normal N condition, which was in line with previous studies in various species [34,51]. However, with little difference in the number of upregulated genes, one thousand more genes were downregulated under N starvation than under normal N condition under drought. In particular, the KEGG pathway of "carbon metabolism" was only enriched in downregulated genes under N starvation, and genes in the KEGG pathways of "ribosome", "arginine biosynthesis pathway", "oxidative phosphorylation" and "aminoacyl-tRNA biosynthesis" were markedly downregulated under N starvation. C and N metabolism has been proven to interactively regulate many aspects of the physiology and development of plants [55]. N assimilation requires energy produced by carbon metabolism, and photosynthetic carbon assimilation requires a large amount of nitrogen [56]. Approximately 4-fold more genes involved in carbon metabolism were downregulated under N starvation than under normal N condition, indicating that the balance of carbon and nitrogen metabolism was disturbed, ultimately strengthening drought stress when N was limited. Ribosome and aminoacyl-tRNA underlie protein synthesis and support cell growth [57]. Arginine also plays a key role in stress responses and plant development, as it is a precursor to the synthesis of polyamines (including putrescine, spermidine and spermine) [58,59]. Oxidative phosphorylation in mitochondria affects the oxidative stress response in plants [60]. Decreased activity of COX, which is a terminal oxidase in oxidative phosphorylation pathways, leads to excessive reactive oxygen species. In summary, we speculated that N deficiency may aggravate drought stress by decreasing energy generation and storage, blocking protein synthesis, and reducing the contents of osmotic regulators. This was partly consistent with a previous study on poplars-adequate N could alleviate drought through enhancing expression of genes involved in plant hormone signal transduction and antioxidant defense, and through interactions between carbon metabolism and nitrogen assimilation [33]. However, these details need to be further studied in S. japonica.
Plants absorb and utilize inorganic N with the help of genes involved in N transport and metabolism. Our transcriptome analysis first identified major genes involved in N uptake and metabolism in S. japonica roots (Table S12). In the present study, most NRTs were downregulated under drought, while AMTs were upregulated. Similar conclusions were observed in poplar [61]. Absorption of more ammonium than nitrate may result in an increased ammonium to nitrate ratio, which could improve the drought resistance of plants [18,20]. The sharp upregulation of NR and NiR may also contribute to increasing the ammonium to nitrate ratio. Upon drought, GS, GOGAT and GDH were upregulated under normal N but downregulated or unchanged under N starvation, indicating that adequate N may promote the N assimilation process to provide substrates for stressresistant components. Under N starvation, the expression of three NRTs (c1261082.graph_c0, c704513.graph_c0 and c1215796.graph_c0) increased, but they were expressed normally under drought or drought with N starvation, suggesting that these genes could only be activated by N starvation but not drought. Phylogenetic analysis with NRTs in Arabidopsis showed that c1261082.graph_c0 and c1215796.graph_c0 were separately clustered together with AtNRT1.1 and AtNRT2.5 ( Figure S9A). In Arabidopsis, the dual-affinity nitrate transporter NRT1.1 and the high-affinity nitrate transporter NRT2.5 mediate nitrate uptake from soil [62]. AtNRT1.1 also participates in auxin transport to mediate lateral root development [63]. The two identified NRTs in S. japonica may exhibit a similar function when the N concentration was low. The N assimilation process was relatively slow upon N starvation, which may be explained by the fact that plants slowed their growth to require less N when N availability was limited [10]. One AMT gene (c1238615.graph_c0) was significantly induced under all three kinds of stresses and was grouped into the AtAMT1 family ( Figure S9B), which consists of high-affinity ammonium transporters, indicating its major role in stress resistance. These N-related genes showed various changes in expression under stress conditions, but their functions need to be further verified. The identification of these N-related genes has greatly promoted our understanding of the N utilization of S. japonica under stress conditions, and provides a molecular basis for further studies.
TFs are important regulators that have various functions in responses to abiotic stresses in plants [5]. In the present study, genes related to the AP2/ERF-ERF, WRKY and NAC families were most abundant among the upregulated TFs under stress. These families have been shown to play an important role in plant-enhanced resistance to stresses [64]. For instance, expression of CsERF1 enhanced cold tolerance in tobacco [65], and overexpression of WRKY11 and ANAC072 increased tolerance to drought and osmotic stresses in Arabidopsis [66,67]. Among the downregulated genes, members of the bHLH and C3H families accounted for the majority upon drought. In contrast to many previous studies indicating that overexpression of several members of the bHLH and C3H families could improve drought resistance [68,69], these TFs may have distinct roles in plants under drought. Furthermore, the C2C2-Dof family had a relatively high proportion of downregulated TFs under N starvation. It was proven that expression of Dof1 in transgenic Arabidopsis improved growth under low-nitrogen conditions [70], highlighting the great roles of C2C2-Dof in N starvation responses. These DEGs in various TF families tightly participate in the molecular regulation of the stress response, which may facilitate the acclimation of S. japonica to adverse environments.
In summary, the physiological and molecular changes play key roles in the resistance to drought and N starvation of S. japonica. At the physiological level, antioxidant enzymes and osmotic regulators help counteract the adverse effects of stress. In present study, S. japonica had high drought tolerance as SOD and POD activities as well as proline contents increased significantly in the roots. At the molecular level, the activation of "phenylpropanoid biosynthesis" and "biosynthesis of amino acids" pathways, the enhanced ability of N uptake and metabolism and the upregulated expressions of transcription factors in AP2/ERF-ERF, WRKY and NAC families may be beneficial to plant-enhanced resistance to either drought or N starvation. On the other hand, N deficiency could aggravate drought stress, as N starvation slowed down the processes of "photosynthesis","ribosome", "arginine biosynthesis", "oxidative phosphorylation" and "aminoacyl-tRNA biosynthesis", indicating that adequate N may be necessary for improving drought tolerance. These results suggest that Sophora japonica can activate stress pathways rapidly via modulating gene expressions in response to drought and N starvation. Thus, studying these genes could improve our knowledge of resistance mechanisms in S. japonica roots to short-term environmental stresses.

Conclusions
In conclusion, we performed a combination of transcriptome sequencing and differential expression analyses to identify pathways and genes altered under drought, N starvation and drought with N starvation treatments. Several physiological traits showed various trends, and transcripts in various metabolic processes reacted rapidly and strongly in response to drought and N starvations stresses. Genes in pathways of phenylpropanoid biosynthesis and biosynthesis of amino acids were greatly activated, and they may play important roles in stress resistance in roots of S. japonica. N deficiency may aggravate drought stress by decreasing energy generation and storage, blocking protein synthesis, and reducing the contents of osmotic regulators. Additionally, many N-related genes and TFs were identified differentially expressed and played key roles in drought and N starvation stresses in this study, but the functions of these genes need to be further studied in the future. The rapid and strong transcriptomic changes of Sophora japonica roots to short-term drought and N starvation reflect an initial response to abiotic stresses, which is a basis of the adaptation to long-term stresses of ancient trees. Our research is the first to explore the molecular mechanism of -and provide valuable information on differentially expressed metabolic pathways and genes of S. japonica under abiotic stresses. This research also reveals the molecular basis for protecting ancient tree resources and breeding resistant varieties for greening.  Figure S9: Neighbor-joining phylogenetic analysis of (A) NTRs and (B) AMTs with related genes in Arabidopsis. The trees were constructed by neighbor-joining phylogeny test, and 1000bootstrao replicates. The protein sequences are provided in Table S3. The green triangles represent the genes proved to have functions to transports nitrate in Arabidopsis. The red stars represent the identified genes involved in N uptake and metabolism at present study, Table S1: RNA sample details including quality and concentration measurements, Table S2: Gene accession numbers of Arabidopsis thaliana used in phylogenetic analysis, Table S3: The primers used for qRT-PCR, Table S4: Summary of sequences analysis, Table S5: Length distribution of assembled unigenes, Table S6: Summary of the annotation results, Table S7: Information of unigene annotations to public databases, Table S8: GO annotation of unigenes, Table S9: KEGG pathway annotations of unigenes, Table S10: The correlation of three biological repeats for transcriptome data, Table S11: GO classifications of DEGs in paired comparisons of different treatments, Table S12: Candidate differentially expressed genes related to nitrogen uptake and metabolism, Table S13: