Differences in growth characteristics under low N stress
Under the normal N supply (CK) treatment, the plant height of the three N-efficient genotypes were 36.73, 41.43, and 67.20 cm, respectively, while under low nitrogen (LN) treatment, the plant heights of the three genotypes were 33.44, 37.32, and 62.62 cm, respectively. Compared with CK treatment, the average plant height decreased by 8.56% under LN treatment after 40 days. For the N-inefficient genotypes, the average plant height decreased by 10.39%(Fig. 1A). The average ground diameter (GDn) of N-efficient genotypes decreased by 2.90%, while the GDn of the N-inefficient genotypes decreased by 5.94% (Fig. 1B). Similarly, under LN treatment, the average fresh weight of the stem (SFW) and dry weight of the stem (SDW) of N-efficient plants decreased by 12.89 and 15.74%, respectively, while the SFW and SDW of the N-inefficient plants decreased by 15.22 and 23.89%, respectively (Fig. 1C). Notably, the fresh weight of the roots (RFW) of the three N-efficient genotypes were 1.435, 1.525, and 3.455 g under CK treatment, respectively, while in the LN treatment, the RFW of the three genotypes increased slightly, to 1.440, 1.535, and 3.515 g, respectively. While the RFW of the N-inefficient genotypes decreased significantly (p < 0.05). By contrast, the dry weight of the roots (RDW) of the N-efficient and N-inefficient genotypes (except the C-1 clone) was reduced significantly under LN treatment (p < 0.05, Fig. 1D, Fig. S1 and S2). Moreover, the leaf fresh weight (LFW) and leaf dry weight (LDW) of a single leaf of the N-efficient and N-inefficient genotypes decreased under LN treatment, and the values of LFW and LDW of the N-efficient genotypes were higher than those of the N-inefficient genotypes under CK or LN treatment (p < 0.05, Fig. 1E, Fig. S1 and S2). As shown in Table S1, the low N adaptation coefficient (LNAC) of GDN, SFW, SDW, RFW, LFW, and LDW of the N-efficient genotypes was significantly higher than that of the N-inefficient genotypes (p < 0.05).
Differences in leaf morphology and chlorophyll content under low N stress
The contents of chlorophyll a (Chl a), chlorophyll b (Chl b), carotenoids (Car), and total chlorophyll (Chl) of plants (N-efficient and N-inefficient genotypes) under LN treatment were significantly lower than those under CK treatment (p < 0.05, Fig. 2A and Fig. S1). Besides, the chlorophyll contents (Chl a, Chl b, Car, and Chl) of the N-efficient genotypes were higher than those of the N-inefficient genotypes under LN or CK treatment (p < 0.05, Fig. 2A and Fig. S1). For example, under LN treatment, the Chl of three N-efficient genotypes were 4.22, 5.19, and 5.99 mg·g− 1, respectively, while the Chl of three N-inefficient genotypes were 3.26, 3.33, and 3.63 mg·g− 1, respectively (Fig. 2A). Compared with CK treatment, leaf morphological traits [leaf length (LL), leaf width (LW), and leaf area (LA)] were reduced under LN treatment (Fig. 2B, Fig. S1 and S2). The LL, LW, and LA of the N-efficient genotypes decreased by 9.43, 8.99, and 15.50%, respectively, while those of the N-inefficient genotypes decreased by 15.22, 14.64, and 30.58%, respectively. In particular, under LN treatment, the LA of the three N-efficient genotypes were 53.42, 54.49, and 83.10 cm2, which were higher than those of the three N-inefficient genotypes (p < 0.05, Fig. 2B). The average LNAC values of Chl a, Chl b, Car, Chl, LL, LW, and LA of the N-efficient genotypes were higher than those of the N-inefficient genotypes,; however, the differences were not significant (Table S1).
Differences in the response of enzyme activities under low N stress
To study the differences between the two contrasting genotypes in response to low N stress, we measured the enzyme activities (nitrate reductase (NR), glutamine synthetase (GS), glutamate dehydrogenase (GDH), and glutamine oxoglutarate aminotransferase (GOGAT)), total amino acid contents (AAs), and soluble sugar contents (SSs) of mixed leaf samples at seven time points (T0, T1, T2, T3, T4, T5, and T6 represent 0, 3, 5, 10, 20, 30, and 40 days of N treatment, respectively). As shown in Fig. S3, the enzyme activities in leaves of the N-efficient (A) and N-inefficient (C) genotypes were inhibited under LN treatment and the AA content decreased significantly after 5 days (T2) of LN treatment (p < 0.05). At T0, T2, T3, T5, and T6 of LN treatment, the NR activity of genotype A was higher than that of genotype C (p < 0.05), and the NR activity reached the lowest at T4 (Fig. S3, Table S2). With increasing treatment time, the activities of GS and GOGAT of genotypes A and C increased first and then decreased, with the turning point of change mostly occurring at T2 (except for the change of GS activity of genotype A). At T5 and T6, the activities of GS and GOGAT of genotypes A and C were decreased significantly under LN treatment (p < 0.05, Fig. S3, Table S2). The activity of GS in genotype A was higher than that in genotype C before LN treatment (p < 0.05). However, after 5 days (T2), there was no significant difference in GS activities between genotypes A and C, and the GS activity of genotype A was significantly lower than that of genotype C after 40 days (T6) of LN treatment (p < 0.05, Table S2). Besides, the GDH activity of genotype A decreased first and then increased, while that of genotype C decreased continuously (Fig. S3). Under CK treatment, there was no significant difference in GDH activity between genotypes A and C (except at T0). After 20 days (T4) of LN treatment, the GDH activity of genotype A was significantly higher than that of genotype C (p < 0.05, Table S2). The AA content of genotypes A and C decreased continuously in response to low N stress (Fig. S3). The AA content of genotypes A and C showed no significant difference at the same time points under CK treatment (except T2); however, the AA content of genotype A was higher than that of genotype C at T3, T4, and T6 under LN treatment (p < 0.05, Table S2). In particular, at the beginning of low N stress, the SS content of genotype A decreased and then increased gradually, while genotype C showed the opposite response (Fig. S3).
Evaluation of RNA sequencing data
According to the change trends of enzyme activities, AAs, and SSs in the process of LN treatment, the mixed leaf samples of N-efficient (A) and N-inefficient (C) genotypes at 0 (T0), 5 (T2), 20 (T4), and 40 (T6) days after LN treatment, with three biological replicates (T0-A-1, T0-A-2, T0-A-3, T0-C-1, T0-C-2, T0-C-3, T2-LN-A-1, T2-LN-A-2, T2-LN-A-3, T2-LN-C-1, T2-LN-C-2, T2-LN-C-3, T4-LN-A-1, T4-LN-A-2, T4-LN-A-3, T4-LN-C-1, T4-LN-C-2, T4-LN-C-3, T6-LN-A-1, T6-LN-A-2, T6-LN-A-3, T6-LN-C-1, T6-LN-C-2, and T6-LN-C-3) were selected for transcriptome sequencing analysis. A total of 36,866,796–61,136,254 clean reads were generated, of which approximately 76.89–87.63% were uniquely mapped to the genome of P. deltoides (Table S3). The Pearson correlation coefficient between the biological replicates ranged from 0.9810 to 0.9977 (Fig. S4), implying that the RNA-seq data were highly reliable.
DEGs between N-efficient and N-inefficient genotypes under low N stress
To study the differences of response to low N stress of the two genotypes (A and C), we detected the DEGs between A and C comparison groups at four-time points (T0-C vs. T0-A, T2-LN-C vs. T2-LN-A, T4-LN-C vs. T4-LN-A, and T6-LN-C vs. T6-LN-A), and 2383; 1243; 1292; and 3055 DEGs were detected, respectively. Compared with genotype C, 1149 DEGs were upregulated and 1906 genes were downregulated in genotype A at T6 (Fig. S5). As shown in the Venn diagram (Fig. S5), there were 397 common DEGs at T0, T2, T4, and T6 between genotypes A and C, and 972, 225, 224, and 1578 specific DEGs at T0, T2, T4, and T6, respectively.
At T0, T2, T4, and T6, the specific DEGs in genotype A were enriched for 40, 32, 34, and 42 Gene Ontology (GO) terms, and more DEGs were mainly enriched in the metabolic process, catalytic activity, and cell. Some DEGs were related to stimulus response in genotype A at different time points (Fig. S6). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of specific and common genes at different time points between genotypes A and C showed that the specific DEGs at T0 were mainly enriched in the biosynthesis of secondary metabolites, glutathione metabolism, starch and sucrose metabolism, phenylalanine metabolism, flavonoid biosynthesis, and thiamine metabolism (Fig. 3). The 225 specific DEGs at T2 were specifically enriched for ascorbate and aldarate metabolism, and one environmental information processing category (ABC transporters). The distinct DEGs at T6 were specifically enriched for metabolism (metabolic pathways, phenylpropanoid biosynthesis, tyrosine metabolism, fatty acid metabolism, pentose and glucuronate interconversions, and indole alkaloid biosynthesis) and plant hormone signal transduction. The common DEGs were mainly enriched in photosynthesis, oxidative phosphorylation, DNA replication, and homologous recombination pathways. These results suggested that genotypes A and C have different response mechanisms to low N stress.
Eight DEGs related to porphyrin and chlorophyll metabolism, and five DEGs related to carotenoid biosynthesis were detected between genotypes A and C, their expression levels in genotype A were higher than those in genotype C after 20 days (T4) of LN treatment (Fig. 4A, B). During the LN treatment, the expression of eight DEGs involved in carbon fixation in photosynthetic organisms in genotype A were higher than those in genotype C (Fig. 4C). After 20 days of low N stress, 22 of 28 DEGs related to biosynthesis of amino acids were expressed at higher levels in genotype A than in genotype C (Fig. 4D). The expression patterns of genes related to plant hormone signal transduction were different between genotype A and genotype C in response to low N stress: 11 of 17 DEGs were expressed at lower levels in genotype A than in genotype C at 40 days (T6) of low N stress (Fig. 4E). Among the seven DEGs related to nitrogen metabolism, four genes were upregulated and three genes were downregulated in genotype A; however, their expression levels were higher than those in genotype C (Fig. 4F).
DEGs during low N stress of N-efficient and N-inefficient genotypes
The DEGs of genotypes A and C at different time points under low N stress were determined to study the differences of response mechanisms of the two genotypes. Three different comparison groups of genotype A at the early stage (T0-A vs. T2-LN-A), the middle stage (T2-LN-A vs. T4-LN-A), and the late stage (T4-LN-A vs. T6-LN-A) were constructed, and 445 (272 upregulated and 173 downregulated), 2337 (1624 upregulated and 713 downregulated), and 692 (219 upregulated and 473 downregulated) DEGs were detected, respectively (Fig. S7). Three corresponding comparison groups of genotype C (early stage: T0-C vs. T2-LN-C, middle stage: T2-LN-C vs. T4-LN-C. and late-stage: T4-LN-C vs. T6-LN-C) were constructed, and 1635 (1231 upregulated and 404 downregulated), 598 (153 upregulated and 445 downregulated) and 226 (132 upregulated and 94 downregulated) DEGs were identified, respectively (Fig. S7). The results showed that the gene response of genotype A was mainly in the middle stage of low N stress, while that of genotype C was mainly in the early stage. There were 191, 319, and 50 common DEGs between genotypes A and C in the early, middle, and late stages of LN treatment, respectively. Notably, there were 926 common DEGs between genotypes A and C in response to low N stress; however, there were more specific DEGs in genotype A (2,062) than in genotype C (1,118; Fig. S7).
The specific DEGs in genotypes A were associated with protein kinase activity, protein phosphorylation, negative regulation of development process, and catalytic activity. By contrast, the DEGs specifically expressed in genotype C were enriched in microtubule-based process, intrinsic component of membrane, and membrane part (Fig. S8). Moreover, during the response of genotypes A and C to low N stress, the expression levels of genes related to metabolic pathways and biosynthesis of secondary metabolites changed significantly. The difference was that in genotype A, many DEGs were enriched in glutathione metabolism, phenylalanine metabolism, phenylpropanoid biosynthesis, ABC transporters, ribosome biogenesis in eukaryotes, and plant hormone signal transduction, whereas in genotype C, many DEGs were enriched in flavonoid biosynthesis, pyruvate metabolism, starch and sucrose metabolism, and steroid biosynthesis (Fig. 5).
Trend analysis of DEGs
For genotypes A and C, DEGs that were detected in response to low N stress at four-time points were clustered into 20 profiles, among which 2988 DEGs of genotype A were significantly enriched in six profiles (profile 0, 6, 9, 10, 12, and 13), among which profile 12 contained the largest number of genes (1168), which was similar to the previous research results, i.e., the genes in genotype A mainly responded in the middle stage of low N stress, and less in the early and late stages (Fig. S9). For genotype C, 2044 DEGs were significantly enriched in four profiles (profile 2, 14, 15, and 17), among which profile 17 contained the largest number of genes (629), which showed that genes mainly responded at the early stage of treatment, and fewer genes participated in the middle and late stages (Fig. S9).
Genes in profile 12 of genotype A (A-profile 12) were mainly involved in kinase activity, transport, localization and binding of various substances, as well as the defense response, whereas the genes in profile 17 of genotype C (C-profile 17) were mainly related to membrane, microtubule, cytoskeletal part, and cytokinesis (Fig. S10). Many genes in A-profile 12 and C-profile 17 were enriched in the pathway of biosynthesis of secondary metabolites. Genes in A-profile 12 were related to the biological processes of N absorption, transformation, transport, and assimilation (glutathione metabolism, phenylalanine metabolism, cysteine, and methionine metabolism, biosynthesis of amino acids, and ABC transporters), while the genes in C-profile 17 were significantly enriched in brassinosteroid biosynthesis, amino sugar, and nucleotide sugar metabolism, fatty acid elongation, and N metabolism (Fig. 6).
Identification of WGCNA modules associated with special traits
A total of 24 leaf samples of genotypes A and C at four-time points (three biological repeats) during low N treatment were used to carry out WGCNA. Based on the results (Fig. 7A), 4106 genes were divided into 19 modules. Except for the grey module with only one gene, the other module sizes ranged from 59 (‘royal blue’) to 539 (‘brown’). Through correlation analysis between modules and experimental traits (NR, GS, GDH, GOGAT, AAs and SSs; Fig. 7B), it was found that the correlation coefficients between GS and ‘magenta’, ‘blue’, ‘light green’, and ‘tan’ were − 0.79, − 0.81, − 0.78, and 0.77, respectively. Moreover, the ‘magenta’ module was not only significantly related to GS, but also correlated highly with GOGAT and AAs (|r| > 0.75, p < 0.05). Therefore, we speculated that genes in the ‘magenta’ module might play a key role in the plant response to low N stress.
Gene expression trends and function analysis of a specific module
The expression trends of genes in the ‘magenta’ module were similar in genotypes A and C, and most of the genes had low expression at T0, after which expression increased gradually with increasing treatment time of low N stress; however, the expression trends of a few other genes were the opposite. Moreover, the expression of these genes in genotype A was more obvious than that in genotype C (Fig. S11). An analysis of functional annotations (Gene Ontology terms) revealed that genes in the ‘magenta’ module are related to the membrane, catalytic activity, enzymatic activity, response to stress, and negative regulation of various biological processes (Fig. 8).
Co-expression network construction and key Gene Mining in specific modules
Based on the correlation of 171 genes in the ‘magenta’ module (Table S4), the top 10 genes in terms of connectivity (within the module) were selected as hub genes, and the gene pairs related to the hub genes with the top 150 weight values were selected to construct the network graph (Fig. 9). The 10 genes (except one TF gene, Podel.18G019200) in the center of the network were hub genes, including genes encoding domain-containing proteins [DUF668 domain-containing protein/DUF3475 domain-containing protein (Podel.02G021400), MACPF domain-containing protein At4g24290-like (Podel.19G001200), ACT domain-containing protein ACR4 isoform X1 (Podel.08G129500), C2 and GRAM domain-containing protein At1g03370-like (Podel.06G003100)], echinoderm microtubule-associated protein-like 6 (Podel.19G035300), ABC transporter C family member 3-like (Podel.03G215800), probable serine/threonine-protein kinase clkA (Podel.03G153500), plant intracellular Ras-group-related LRR protein 3-like (Podel.10G040100), UDP-glucoronosyl/UDP-glucosyl transferase family protein UDP (Podel.04G076900), and random slug protein 5-like (Podel.05G137500). Among them, Podel.03G153500 had the strongest connectivity with all other genes and participates in the plant response to stress (Table S5).
The expression levels of most hub genes in the leaves of genotypes A and C decreased slightly in the early stage of low N stress (from T0 to T2), and increased gradually from T2 to T6. Moreover, the expression levels of hub genes in genotype A were lower than those in genotype C at T0, but it was the opposite at T6, indicating more obvious changes in the expressions of hub genes in genotype A (Fig. 9, Table S4).
To further explore the key genes in the ‘magenta’ module, TF annotation was carried out. Finally, 12 TFs were annotated (Table S6), including WRKY TF family genes [PodelWRKY41 (Podel.01G096500), PodelWRKY75 (Podel.15G104200), and PodelWRKY18 (Podel.18G019200)], bHLHTF family genes [PodelBHLH25 (Podel.01G303700), PodelBHLH30 (Podel.02G115300), and PodelBHLH (Podel.06G162000)], LOB TF family genes [PodelLBD37 (Podel.07G061600) and PodelLBD1 (Podel.08G051900)], ARR-B TF family genes [PodelWER (Podel.12G083500)], bZIP TF family genes [PodelHY5 (Podel.06G265700)], SBP TF family genes [PodelSPL4 (Podel.01G421600)], and zf-HD TF family genes [PodelZHD4 (Podel.12G042700)]. Among them, the connectivity between PodelWRKY18 and other genes in the ‘magenta’ module was second only to the hub genes, and the changing trend of its expression with the duration of low N stress was consistent with that of the hub genes (Fig. 9, Table S6).
Difference in gene expression related to N metabolism under N stress
The expression levels of genes related to N uptake, transport, and assimilation in genotypes A and C during LN treatment were detected using quantitative real-time reverse transcription PCR (qRT-PCR), and the accuracy of transcriptome data was verified by comparing the results of qRT-PCR and RNA sequencing (Fig. S12). Under LN treatment, the expression levels of nitrate (NO3−) transport related genes [NRT1;1 (Podel.03G118100) and NRT1;2 (Podel.12G073800)] were inhibited, and the change trends in genotypes A and C were similar. The expression of AMT1;6 (Podel.09G046000) was inhibited in genotype A and promoted in genotype C, while the expression of AMT2;1 (Podel.06G111800) showed the reverse results. In genotype A, the expression levels of NR (Podel.05G183500) and NiR (Podel.04G146000) decreased at first (from T0 to T4) and then increased (from T4 to T6). In genotype C, these genes were inhibited all the time, but promoted in at T2 and T6. In genotypes A and C, the expression of GS2 (Podel.10G024300) was upregulated at T2 and gradually decreased from T2 to T6, to a point that was lower than that at T0, and the inhibition in genotype A was stronger than that in genotype C. The expression of GDH2 (Podel.15G117500) decreased from T0 to T2 and then gradually increased from T2 to T6 in genotype A, with the expression at T6 being higher than that at T0. In genotype C, the expression decreased from T0 to T6.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.