Transcriptome analysis for Fraxinus mandshurica seedlings from different carbon fixation and growth provenances in response to nitrogen deficiency

doi:10.21203/rs.2.11242/v1

Background

The carbon fixation characteristic of F. mandshurica seedlings from 20 provenances has been evaluated by our research group. In order to explore whether growth and carbon sequestration are related to nutritional adaptability, the foliage and roots of the annual Wuchang (WC) seedlings with high carbon content and Hailin (HL) seedlings with low carbon content which were grown in nitrogen-deficient nutrition and total nutrition were used for RNA-seq determination.

Results

Eight transcriptome libraries by high-throughput sequence were analysis. 1,235,174,984 clean reads and 88,655 unigenes with N50 length of 1,259 bp were obtained. Under normal nitrogen condition, 783 differentially expressed genes (DEGs) between WC and HL were identified, the number in foliage (669) of DEGs between HL and WC was more than that in roots (149). The number of transcription factors (TFs), hormone, and Protein kinase (PK) genes was significantly more in WC than that in HL. Compared to the normal nitrogen, 8173 DEGs related to nitrogen deficiency were identified and the number of DEGs in roots (6999) was more than that in foliage (1616). Several nitrogen deficiency-related metabolic pathways and many DEGs involved in nitrogen absorption and assimilation, carbon metabolism, hormones, transcription factors and kinases were identified. The numbers of DEGs encoding nitrate transporters, PK, TFs and hormone in WC were less than that in HL, which indicated the response of HL provenance seedlings to nitrogen deficiency was stronger than that of WC provenance seedlings. However, under nitrogen deficiency, the number of up-regulated DEGs in the WC provenance seedlings is more than that in the HL, which indicated that the WC provenance seedlings were more tolerant to nitrogen deficiency than HL.

Conclusions

The data showed that high carbon content and high growth rate of WC provenance seedlings are mainly attributed to the high transcriptional expression of many metabolic genes in foliage. The response of HL provenance seedlings to nitrogen deficiency is significantly greater than that of WC provenances seedlings, but WC provenance seedlings were more tolerant to nitrogen deficiency than HL. Many genes related to nitrogen deficiency were identified, which will expand our current understanding of nitrogen responses.

Plant Molecular Biology and Genetics

Fraxinus mandshurica

nitrogen deficiency

provenance

RNA sequencing

Provenances are the geographic locations from which seed sources originate. The tree provenances experiment is an experiment in which seedlings of a specie of different provenances are planted together. The provenance experiment can study the genetic diversity within the species, the interaction between the provenance and the ecological environment and the adaptability of the provenance. The aim of the provenance experiment is to select an excellent provenance of the trait of interest and then plant the selected elite provenance in a suitable ecological zone. Numerous studies indicated that many traits have extensive genetic differences between different provenances, such as the density fluctuations [ 1], the growth and stem formation[ 2], radial growth[ 3], terpenoid patterns in needles[ 4]. What’s more, Brynjar Skulason et al. found that subalpine fir (Abies lasiocarpa) provenances showed significant differences for all measured traits[ 5].

In the past, studies on the differences in provenances of forest trees were mainly focused on the growth traits (tree height, DBH)[ 6, 7], physiological and biochemical traits[ 8], fitting growth model[ 9, 10] and identification of molecular markers related to provenance[ 11]. However, there were few reports on the differences in RNA-Seq between provenances. Some scholars believe that the differences between provenances are mainly determined by genetic factors[ 12, 13], while others believe that differences between provenances are mainly caused by the environment, such as soil[ 14], climate[ 1, 15, 16], temperature[ 17-19], etc. In addition, the effects of nutrition on provenance growth have also been extensively studied. Xabier Santiso et al. found that there is only a small difference among Arbutus unedo L. provenances under nitrogen nutrition stress[ 20]. But significant differences in growth between different provenances after treatment with nutrition was detected for Picea abies Karst., Fagus sylvatica L[ 7].

Nitrogen is one of the essential nutrients in plant growth and development[ 21, 22]. It is an important component of many important macromolecules, including proteins, enzymes and several important plant hormones. Plants often suffer from nitrogen deficiency. Fertilization will not only increase the cost of production, but also cause environmental pollution. Therefore, it is very important to develop or breed low-nitrogen-tolerant varieties. Nitrogen in the soil is first absorbed by nitrate transporters (NRTs) and ammonium transporters (AMTs) in an active transport manner in plant roots. Nitrate is reduced to ammonium by nitrate reductase (NR) and nitrite reductase (NiR) and then ammonium is transferred to the carbon skeleton to form glutamic acid by glutamine synthetase (GS) and glutamate synthetase (GOGAT). Nitrogen deficiency not only affects nitrogen metabolism but also carbon metabolism, such as sucrose phosphate synthase (SPS), sucrose synthase (SS), and soluble acid invertase (AI). In addition, phytohormones also play an important role in nitrogen deficiency, especially cytokinins. The more we understand the metabolism of nitrogen, the more we have tools to increase tolerance to nitrogen deficiency in plants. Many researchers have conducted extensive and in-depth research on nitrogen deficiency in many species from growth, yield, physiology and biochemistry, genes, and transcriptional regulation. With the increased accuracy and reduced cost of high-throughput sequencing technologies, RNA-Seq has been used in a growing number of species. RNA-Seq can reveal the mechanism of nitrogen deficiency tolerance from the transcriptional expression of genes, and also identify many nitrogen deficiency-related genes[ 23, 24], which will increase our understanding of the mechanism of nitrogen metabolism. However, RNA-Seq analysis of the differences between Fraxinus mandshurica Rupr. different provenances with nitrogen deficiency has not been reported yet.

Fraxinus mandshurica Rupr. is distributed in the northwest and northeast of China, the Russian Far East, the northern part of the Korean Peninsula and northern Japan. It is a valuable hardwood species in the forest areas of northeast China. Because this species is distributed widely, the differences of the climate, soil and topography in the distribution area are significant and result in geographical and reproductive isolation. Under long-term natural selection, Fraxinus mandshurica Rupr. differentiates to form different provenances. The early research on the provenance of Fraxinus mandshurica Rupr. is mainly carried out by Jilin Forestry Research Institute, and recently it was mainly carried out by the research group of Northeast Forestry University. These studies were focused on good provenance selection based on growth traits (tree height), but studies on the selection of good provenances for carbon fixation have not been reported. The characteristic of carbon fixation of 20 Fraxinus mandshurica provenances seedlings has been evaluated by our research group and has been classified high, middle and low carbon fixation groups. Especially, the carbon content and carbon fixation quality of the Wuchang (WC) provenance seedlings is significantly higher than that of the Hailin (HL). In this study, in order to explore the mechanism of high carbon fixed provenance formation, especially the response mechanism under nitrogen deficiency stress, the annual WC seedlings with high carbon content and HL seedlings with low carbon content were selected for RNA-seq determination. What’s more, to identify genes associated with nitrogen deficiency and explore the differences in response to nitrogen deficiency, the two provenances seedlings treated with nitrogen deficiency were also used for RNA-seq determination.

RNA-seq and de novo transcriptome assembly

To obtain an overall view of the nitrogen deficiency transcriptome in the two provenances, RNA samples were prepared from the roots and foliage of both provenances on 30 days after nitrogen deficiency treatment. RIN values indicate that RNA quality can be used for cDNA library construction (Additional file 2). Gene expression profiles of the high carbon fixation Fraxinus mandshurica provenance (WC) and low carbon fixation Fraxinus mandshurica provenance(HL) seedlings roots and foliage under both normal nitrogen and nitrogen deficiency conditions were analyzed. For the eight samples, three biological replicates were performed in sequencing. In total, 1,235,174,984 clean reads originated from the 24 Fraxinus mandshurica cDNA libraries were obtained following removal of adaptors and low-quality reads (Additional file 3). The percent G+ C and fraction of bases with quality scores of Q30 for the 24 libraries averaged 45.2% and 94.2%, respectively. These numbers showed that the data was of high quality. Trinity software was then used to assemble the clean reads into 88,655 unigenes with N50 length of 1259 bp. The relative length of the assembled sequences is one standard by which a transcriptome assembly can be assessed and the summary statistics are shown in Fig.1. In the final assemblies there were 48,754 unigenes ranging from 200 to 500 bp in length, which accounted for 55% of the total. There were 19,220 (21.68%) and 20,681 (23.33%) unigenes from 500 to 1000 bp and > 1000 bp in length, respectively, and assembled sequences < 200 bp were not taken into account.

Annotation and functional classification of unigenes

The assembled sequences were used as queries in BLAST searches (E ≤ 10⁻⁵) against the NR, KOG, Swissprot and KEGG databases, and a total of 55,537 were annotated (Fig. 2a). The majority of the unigenes (53,603; 60.5%) were annotated from the NR database. In contrast, only 24,959 (28.2%), 35,879 (40.5%) and 42,650 (48.1%) annotated unigenes were matched to the KEGG, KOG and Swissprot database, respectively. As shown in Fig.2a, a total of 21,084 unigenes were annotated with all four databases. To further utilize the transcriptome data, the KOG, GO, and KEGG databases were used to identify the molecular processes that occur during nitrogen deficiency of Fraxinus mandshurica seedlings roots and folium. The annotated unigenes were initially classified into 25 KOG categories with the largest group being “General function prediction only” (10,345; 28.83%) (Fig.2b). The Gene Ontology (GO), which is an internationally standardized classification system, was used to assign 15,827 unigenes to the three principal GO domains: “cellular component”, “molecular function” and “biological process” (Fig.2c). In the “cellular component” category, “cell” and “cell part” were the most abundant terms. For the category of “molecular function”, “catalytic activity” and “binding” were the most prominent. “Metabolic process” and “cellular process” were the most abundant terms in the “biological process” category.

Identification of differentially expressed genes (DEGs)

The transcriptional levels were normalized using the FPKM method. Meanwhile, FDR<0.05 was used as screening thresholds to test the significance of difference in transcript abundance. From eight comparisons, including CP1 (WCKF vs. HCKF), CP2 (WCKR vs. HCKR), CP3 (WCKF vs. WCTF), CP4 (WCKR vs. WTR), CP5 (HCKR vs. HTR), CP6 (HCKF vs. HTF), CP7 (WTF vs. HTF), CP8 (WTR vs. HTR), a large number of differentially expressed genes (DEGs) were identified (table1). The numbers of DEGs detected were as follows: CP1 669 (285 up- and 384 down-regulated), CP2 149 (77 and 72), CP3 305 (165 and 140), CP4 4225 (1907 and 2318), CP5 4531 (2386 and 2145), CP6 1414 (991 and 423), CP7 141 (67 and 74), CP8 178 (44 and 134). More down-regulated DEGs were found in CP1 and CP8, which indicated that the numbers of DEGs up-regulated in WC were 34.74% and 204.55% greater than that in HL under normal nitrogen for foliage and nitrogen deficiency for roots. In order to study the difference in DEGs between two provenances, the DEGs in CP1, CP2 and CP7, CP8 were further analyzed together. Consequently, 1057 DEGs were identified using pair-wise comparison of each accession between WC and HL under normal condition and nitrogen deficiency (Additional file 4). 783 DEGs were identified between WC and HL under normal condition, including both up-regulated (345) and down-regulated (438) genes (fig.3a and fig.3b). Interestingly, DEGs in the CP1 (669) were nearly 4.5 times as much as those in the CP2 (149) (table1), which indicated that the number of DEGs in foliage was significantly higher than that in roots under normal nitrogen. Similarly, 300 DEGs were identified between WC and HL under nitrogen deficiency, including both up-regulated (102) and down-regulated (198) genes (fig.3a and fig.3b). And down-regulated DEGs in the CP8 (134) were nearly twice as much as those in the CP7 (74) (table1), which indicated that the number of DEGs in roots was significantly higher than that in foliage under nitrogen deficiency.

Moreover, analysis of DEGs in CP3, CP4 and CP5, CP6 can be used to identify genes that respond to nitrogen deficiency in WC and HL, and then the differences in response to nitrogen deficiency between the two provenances can be further analyzed. Consequently, 8173 DEGs under nitrogen deficiency were identified using pair-wise comparison of each accession between normal and nitrogen deficiency in WC and HL (Additional file 4: Figure S2). However, there were 1959 DEGs commonly found in both HL and WC. 4447 DEGs were identified between normal nitrogen and nitrogen deficiency in WC, including both up-regulated (2019) and down-regulated (2431) genes (fig.3c and fig.3d). Interestingly, DEGs in the CP4 (4225) were nearly 14 times as much as those in the CP3 (305), which indicated that the number of DEGs in roots was significantly higher than that in foliage in WC. Similarly, 5685 DEGs were identified between normal and nitrogen deficiency in HL, including both up-regulated (3210) and down-regulated (2518) genes (fig.3c and fig.3d). And DEGs in the CP5 (4531) were nearly 3.2 times as much as those in the CP6 (1414), which indicated that the number of DEGs in roots was significantly higher than that in foliage in HL. In addition, the transcriptional expression of DEGs in the two provenances foliage was significantly different in response to nitrogen deficiency. The two provenances displayed dissimilar expression patterns, in which the mount of total (1414), up-regulated (991) and down-regulated (423) DEGs in the CP6 were 4.6 times, 6.0 times and 3.0 times as much as those (305,165 and 140) in the CP3(table1), respectively, which indicated that the foliage of HL provenance responded to nitrogen deficiency more strongly than that of WC provenance.

To validate the gene expression results derived from the transcriptome data, nine DEGs related to transcription factors in CP1 and six DEGs involved in the phenylpropanoid pathway for the biosynthesis of lignin in CP4 and CP5 were selected for qRT-PCR. The results of experiment showed a strong correlation with the RNA-Seq data (Fig.4), which indicates that the results of the transcriptomic analysis are reliable.

GO annotation for DEGs

Gene annotation of DEGs was using GO database. Annotated genes were grouped into three major functional categories: cellular component, molecular function and biological process, and then were divided into subcategories (Fig.5; Additional file 5). For ease of presentation, data was expressed as the average of genes in each subcategory in WCKF vs. the average in HCKF. Within the category of biological process in WCKF vs. HCKF (CP1), genes involved in single-organism process, cellular process, metabolic process, localization, biological regulation, response to stimulus and developmental process were 11.32%,11.80%,13.24%, 3.03%, 3.03%, 2.87% and 1.28% of the total, respectively. In the cellular component category, 5.10%, 4.47%, 3.83%, 3.83%, 1.28% and 1.12% genes were related to cell, membrane, organelle, membrane part, organelle part and cell junction, respectively. In the molecular function category, 14.99%, 9.73% and 1.75% genes were involved in catalytic activity, binding and transporter activity, respectively (Fig.5). Annotation for DEGs in CP2, CP3, CP4, CP5, CP6, CP7 and CP8 pairwise comparisons was also carried out. In all three functional categories, genes belonging to the above three groupings showed a similar distribution pattern to those of CP1. However, macromolecular complex showed higher levels in CP3, CP4, CP5 and CP6, which suggested that macromolecular complex might play an important role in responding to nitrogen deficiency (Fig.5; Additional file 5).

KEGG pathway analysis for DEGs

Transcriptome sequencing can also be taken to identify a variety of metabolic processes. In this study, a total of 1057 DEGs were confirmed in CP1, CP2, CP7 and CP8 and 8173 DEGs between normal and nitrogen deficiency were identified in CP3, CP4, CP5 and CP6. After identifying DEGs, KEGG annotation was used for DEGs from the above comparisons. In CP1, CP2, CP3, CP4, CP5, CP6, CP7 and CP8,126, 19, 63, 1298, 1165, 232, 21 and 30 DEGs were assigned to 73, 24, 47, 111, 117, 90, 20 and 30 pathways of the KEGG database, respectively (Additional file 6). The most reliable significantly (p<0.05) enriched pathways of DEGs were represented as Additional file 3. Integrate CP3, CP4, CP5 and CP6 comparisons, the major pathways involved in the mechanism for responding to nitrogen deficiency were “Phenylpropanoid biosynthesis”, “Nitrogen metabolism”, “Diterpenoid biosynthesis”, “Isoquinoline alkaloid biosynthesis”. In addition, in the roots that responded most to nitrogen deficiency, “Ribosome” pathway had the highest number of DEGs in both provenances. “Oxidative phosphorylation” and “Starch and sucrose metabolism” pathways were specifically enriched in the roots of WC and HL, respectively. Similarly，because only 2 DEGs in CP2 are annotated by KEGG, we used the DEGs in CP1 to analyze the difference between two provenances. In CP1, there were 18 pathways which could be divided into eight types involving in the growth mechanism: “Carbohydrate metabolism”, “Energy metabolism”, “Lipid metabolism”, “Amino acid metabolism”, “Metabolism of other amino acids”, “Metabolism of terpenoids and polyketides”, “Biosynthesis of other secondary metabolites” and “Signal transduction” (Additional file 6).

DEGs involved in carbon and energy metabolism

The expression levels of many genes involved in carbohydrate metabolism, lipid metabolism and energy metabolism were differentially between HL and WC. For example, six genes involved in “starch and sucrose metabolism” pathway in CP1 showed increased transcript abundance in WC (Additional file 6, Additional file8). These genes encoded polygalacturonase (Unigene0056877), beta-glucosidase (Unigene0064953, Unigene0055864, Unigene0054701, Unigene0040728) and trehalose-phosphate phosphatase J(Unigene0037044) were identified and their expressions were increased by 2.3 to 42-fold. Five DEGs involved in “Nitrogen metabolism” pathway in CP1 showed increased transcript abundance in WC (Additional file 6). These genes encoded carbonic anhydrase (Unigene0063117, Unigene0054409, Unigene0063118), glutamine synthetase (Unigene0053381), high affinity nitrate transporter 2.5 (Unigene0050395) were identified and their expressions were increased by 3.4 to 16.7-fold. Nine DEGs involved in Lipid metabolism in CP1 showed increased transcript abundance in WC (Additional file 6). These genes encoded lipoxygenase (Unigene0062168, Unigene0059365), 3-ketoacyl-CoA synthase (Unigene0041866, Unigene0041868, Unigene0030198), 4-coumarate--CoA ligase (4CL6, Unigene0050299), 3-ketoacyl-CoA thiolase 2 (Unigene0043651), Jasmonic acid carboxyl methyltransferase (Unigene0030475), salicylic acid methyltransferase (SAMT, Unigene0060674) were identified and their expressions were increased by 2.2 to 13.0-fold.

DEGs related to hormone signaling

A total of 71 DEGs were involved in several plant hormone signal transduction pathways in our study, such as jasmonic acid(JA), auxin, cytikinine(CK), ethylene, salicylic acid(SA), gibberellin(GA), brassinosteroid and abscisic acid(ABA) (Additional file 7). Among them, auxin, JA and SA pathways were predominant. According to the above analysis, many genes involved in plant hormone signal transduction pathway showed increased transcript abundance in WC. Ten genes involved in “Plant hormone signal transduction” and encoded Zinc-finger protein (ZIM) precursor (TIFY) (Unigene0007822, Unigene0058809, Unigene0052515, Unigene0052516, Unigene0058810), ethylene-responsive transcription factor 1B-like (Unigene0033087), auxin-induced protein 15A-like (Unigene0043970),auxin response factor 9-like (Unigene0064310), Auxin-responsive GH3 family protein (Unigene0056132), histidine-containing phosphotransfer protein 4-like (Unigene0081003) were identified and their expressions were increased by 2.3 to 13.9-fold (Additional file 7).

Other differentially regulated genes

Many genes involved in pathways of “Phenylpropanoid biosynthesis”, “Glucosinolate biosynthesis” “Terpenoid backbone biosynthesis”, “Monoterpenoid biosynthesis”, “Diterpenoid biosynthesis” and “Cyanoamino acid metabolism” showed increased transcript abundance in WC under normal nitrogen in CP1 (Additional file 6). For example, seven genes involved in “Phenylpropanoid biosynthesis” pathway which encoded beta-glucosidase (Unigene0055864, Unigene0040728, Unigene0064953 and Unigene0054701), cinnamyl alcohol dehydrogenase 6 (Unigene0045855), peroxidase 19 (Unigene0066438), caffeic acid 3-O-methyltransferase (Unigene0061140) were identified and their expressions were increased by 2.3 to 40.8-fold (Additional file 6). Nine genes involved in “Terpenoid backbone biosynthesis”, “Monoterpenoid biosynthesis” and “Diterpenoid biosynthesis” pathways which encoded cytochrome P450 CYP82U4 (Unigene0059965),terpene synthase (Unigene0055252, Unigene0055251, Unigene0062513, Unigene0009791), geraniol synthase (Unigene0059790),1-deoxy-D-xylulose-5-phosphate reductoisomerase (Unigene0042530), heterodimeric geranylgeranyl pyrophosphate synthase small subunit (Unigene0040100), 1-deoxy-D-xylulose-5-phosphate synthase 2 (Unigene0033472) were identified and their expressions were increased by 2.4 to 19.7-fold (Additional file 7). Five genes involved in “Cyanoamino acid metabolism” pathways which encoded beta-glucosidase (Unigene0054701, Unigene0064953, Unigene0040728, Unigene0055864), isoleucine N-monooxygenase 1-like (Unigene0035112) and their expressions were increased by 2.3 to 40.8-fold (Additional file 6, Additional file 8).

In addition, protein kinases and transcription factors (TFs) under normal nitrogen and nitrogen deficiency were also analyzed (Additional file 9). In CP1, CP2, CP7 and CP8, a total of 23, 0, 8 and 0 DEGs encoding protein kinases and a total of 28, 2, 2 and 2 DEGs encoding transcription factors were found, respectively (Additional file 10 and Additional file 11). This indicated that the transcription levels of protein kinases and transcription factors in foliage were different between the two provenances under both nitrogen deficiency conditions and normal nitrogen. Therefore, the transcriptional expression of protein kinases and transcription factors in foliage (CP1) was further analyzed. A total of 23(3up- and 20down-regulated) and 28(11up- and 17down-regulated) DEGs encoding protein kinases and transcription factors were found in CP1, respectively (Additional file 10 and Additional file 11). Among DEGs encoding protein kinases, 16 (2 up- and 14 down-regulated) receptor-like serine/threonine-protein kinase (RLKs), 2 (0 and 2) CBL-interacting serine/threonine-protein kinase 14-like (CIPKs) and 3 (1 and 2) Serine/threonine-protein kinase (STPKs) were identified, respectively (Additional file 10). Among DEGs encoding transcription factors, 7 (2 up- and 5 down-regulated) AP2/ERF, 9 (3 and 6) bHLH, 3 (3 and 0) MYB, 2 (2 and 0) NAC, 2 (0 and 2) TCP and 2(0 and 2) WRKY were identified, respectively (Additional file 11, Additional file 8). The results indicated that the protein kinases and transcription factors were more active in WC rather than in HL.

DEGs between normal nitrogen and nitrogen deficiency in Fraxinus mandshurica

In order to study the two Fraxinus mandshurica provenances transcriptional expression difference in response to nitrogen deficiency, DEGs between normal condition and nitrogen deficiency from roots and foliage of WC and HL were obtained, respectively. KEGG enrichment analysis showed that the transcriptional expression of foliage from two provenances was significantly different in response to nitrogen deficiency. For example, the number of KEGG pathways which were significantly enriched was 8 in WC foliage (CP3), less than that (14) in HL foliage (CP6). The pathways with a larger number of genes enriched in WC foliage (CP3) were “Sulfur metabolism” and “Nitrogen metabolism”, however, that in HL foliage (CP6) were “Plant-pathogen interaction”, “Plant hormone signal transduction”, “Protein processing in endoplasmic reticulum”, “Phenylpropanoid biosynthesis” and “Starch and sucrose metabolism” (Fig6d). This indicated that the response patterns of WC and HL to nitrogen deficiency were different in the foliage. Similarly, the DEGs of the roots were analyzed. For example, the number of KEGG pathways which were significantly enriched was 8 in WC roots (CP4), less than that (17) in HL roots (CP5). The pathways with a larger number of genes enriched in WC roots (CP4) were “Ribosome”, “Oxidative phosphorylation”, “Phenylpropanoid biosynthesis”, “Cyanoamino acid metabolism” and “Diterpenoid biosynthesis”, however, that in HL roots (CP5) were “Ribosome”, “Starch and sucrose metabolism”, “Oxidative phosphorylation” , “Plant hormone signal transduction”, “Pentose and glucuronate interconversions”, “Tyrosine metabolism” and “Phenylpropanoid biosynthesis” (Fig.6c). This indicated the response patterns of WC and HL to nitrogen deficiency were also different in the roots.

Protein kinases and transcription factors (TFs) responding to nitrogen deficiency

A total of 14, 92, 132 and 152 DEGs encoding protein kinases were found in CP3, CP4, CP5 and CP6, respectively (Additional file 10, Additional file 8). Among them, 14 (11 up- and 3 down-regulated), 55 (52 and 3), 86 (80 and 6) and 118 (115 and 3) receptor-like kinases (RLKs) were identified. And there were 0 (0 up- and 0 down-regulated), 3 (3 and 0), 5 (4 and 1) and 5 (5 and 0) genes of CP3, CP4, CP5 and CP6 involved in the mitogen-activated protein kinases cascade. We found 0, 1, 2 and 2 DEGs of calcium-dependent protein kinases family, 0, 1, 2 and 1 DEGs of CBL-interacting protein kinases (CIPK) family and 0, 1, 5 and 2 protein phosphatase 2C were all up regulated in CP3, CP4, CP5 and CP6. Moreover, we identified 0, 12, 18 and 10 DEGs of Serine/threonine-protein kinase (STPK) in CP3, CP4, CP5 and CP6. The details of these PK and protein phosphatase genes are shown in Additional file 10. These results indicated that protein phosphorylation may be induced by nitrogen deficiency in Fraxinus mandshurica. In addition, it was interesting that there were more up-regulated PK in HL(CP5) than WC(CP4) (Fig.7a and Fig.7b).

In this study, totally 4 (1 up- and 3 down-regulated), 52 (34 and 18), 84 (63 and 21) and 60 (52 and 8) transcription factors (TF) genes were identified in CP3, CP4, CP5 and CP6 (Additional file 11). Obviously, the number of up-regulated expressed transcription TFs in CP5 (CP6) was greater than that of CP4 (CP3) (Fig.7c-f), which indicated HL is stronger than WC in response to nitrogen deficiency. Generally, the major TF families were AP2/ERF, WRKY, bHLH, MYB, MYB-related, NAC and bZIP. In the CP4, the MYB family was the largest group (21%), followed by the AP/ERF (17%) and bHLH (13%) family (Additional file 9: Figure S4); for CP5, the top three family were AP/ERF (21%), bHLH (21%) and MYB (19%) (Additional file 9: Figure S5); and for CP6, the top three were WRKY (35%), bHLH (27%) and AP/ERF (17%) (Additional file 9: Figure S6); however, for CP3, only 4 differentially expressed transcription factor genes were identified (Fig.7c). Of AP2/ERF genes, 3, 14 and 10 were up-regulated in CP4, CP5 and CP6, respectively; of bHLH genes, 6, 17 and 13; of MYB genes, 6, 10 and 2; of WRKY genes, 2, 4 and 20. The results indicated that in the response to nitrogen deficiency, bHLH, MYB and AP2/ERF play a major role in roots, and WRKY plays a major role in foliage (Fig.7c-f).

DEGs related to hormone signaling responding to nitrogen deficiency

In addition to the basic roles in growth and development, phytohormones are also involved in various environmental responses, such as light, salt and drought. It has been proposed that some hormones coordinate demand and acquisition of nitrogen. In this study, the number of up-regulated DEGs related to hormone signaling was greater than that of down-regulated DEGs in respond to nitrogen deficiency. Totally 0 (0 up- and 0 down-regulated), 22 (17 and 5), 40 (34 and 6) and 15 (14 and 1) DEGs were identified in CP3, CP4, CP5 and CP6 (Additional file 7). Obviously, the number of up-regulated DEGs related to hormone signaling in CP5 (CP6) was greater than that of CP4 (CP3), which indicated the response to nitrogen deficiency, HL is stronger than WC. Among them, auxin, ABA, SA, JA and CK pathways were predominatly induced. In the auxin signaling pathway, there were 8, 19 and 2 genes up-regulated in CP4, CP5 and CP6; in the SA pathway, the number of up-regulated genes was 5, 5 and 1; in the ABA pathway, the number of up-regulated genes was 0, 5 and 3; in the JA pathway, the number of up-regulated genes was 0, 1 and 7; in the CK pathway, the number of up-regulated genes was 2, 5 and 0 (Additional file 7). In addition, the response patterns of DEGs related to hormones under nitrogen deficiency in roots and foliage were different. For example, the hormones that predominantly played a role in the foliage were JA and ABA but in the roots were auxin, SA, CK and ABA.

DEGs related to nitrogen assimilation responding to nitrogen deficiency

Many genes involved in nitrogen absorption and assimilation were differentially expressed under nitrogen deficiency relative to the normal condition. In the current study, 10 DEGs encoding nitrate transporters were detected (Additional file 12, Additional file 8). Compared to normal nitrogen, 7 DEGs were up-regulated and 3 DEGs were down-regulated under nitrogen deficiency. Whereas, the number of these DEGs in WC was less than that in HL, and there were 2 up-regulated DEGs (Unigene0034112 and Unigene0050395) which encoded NRT3.2 and NRT2.5 and 2 down-regulated DEGs (Unigene0046155 and Unigene0047301) which encode NRT2.7 and NRT in both WC and HL. In addition, four DEGs (Unigene0018686, Unigene0022313, Unigene0026134, Unigene0050965) which encoded NRT3.1 and NRT2.1 were up-regulated and one DEG (Unigene0039451) which encoded NRT3.1 was down-regulated only in HL, and one DEG (Unigene0071953) which encoded NRT was up-regulated in WC, respectively. Four nitrogen deficiency responsive genes encoding ammonium transporters (Unigene0034277, Unigene0045039, Unigene0045040, Unigene0054529) were identified (Additional file 12), in which, Unigene0054529 was responsive to nitrogen deficiency only in HL. Moreover, fourteen, six and four DEGs encoding amino acid, lysine histidine and oligopeptide transporters were found (Additional file 12), respectively. There were also three, six, two and two DEGs encoding glutamine synthetase (GS), glutamate dehydrogenase (GDH), nitrite reductases (NiR) and nitrate reductase (NR) which were the key enzymes in nitrate assimilation and their transcription expression were all down-regulated except two GS unigenes (Unigene0009471, Unigene0053381).

In this study, there were also one, one, three and one DEGs encoding Dof zinc finger protein DOF1.5-like (CDF4), Pyruvate kinase (PYK1), Ribulose-1,5-bisphosphate carboxylase (RuBP) and neutral invertase A (INVA) which were the key enzymes in carbon assimilation and their transcription expression were all down-regulated except one INVA unigene (Additional file 12). Meanwhile, expression of many genes associated with absorption or translocation of other nutrients changed under nitrogen deficiency, such as phosphate (9), potassium (3), sulfate (5), Zinc (4), iron (3) and molybdate transporter (2), indicating that uptake of these nutrients in Fraxinus mandshurica is affected by nitrogen metabolism under cross-talking regulation.

DEGs between WC and HL under normal nitrogen

Under normal conditions, compared with HL provenance, there were 669 DEGs in the foliage of WC provenance, of which 285 were down-regulated and 384 were up-regulated; 149 DEGs were in roots, of which 77 were down-regulated, and 72 were up-regulated. It indicated that the difference in growth rate of the two provenances under normal nitrogen is mainly caused by the DEGs in the foliage. Through comparative analysis, it was found that the number of up-regulated genes in WC provenance foliage was more than that in HL. Through KEGG annotation analysis, it was found that the DEGs were significantly enriched in 8 metabolic processes such as biosynthesis of other secondary metabolites, signal transduction, metabolism of terpenoids and polyketides, carbohydrate metabolism, energy metabolism, lipid metabolism, amino acid metabolism and metabolism of other amino acids, that is, 17 metabolic pathways (Additional file 6: CP1). Further analysis of DEGs on these significant enrichment pathways revealed that transcription factor, hormone, and kinase-related genes were significantly more in WC provenance than in HL, but among the DEGs of carbon metabolism, nitrogen metabolism and lipid metabolism, only the number of lipid-related genes in WC provenance was significantly higher than that in HL. Five of the nine lipid metabolism-related differentially expressed genes are in the alpha-Linolenic acid metabolism pathway, and the transcription level of four jasmonic acid synthesis-related genes (4CL5, lipoxygenase 2.1, 3-ketoacyl-CoA thiolase 2, SAMT) was significantly higher in WC than that in HL. Correspondingly, the JAZ protein transcription level in the jasmonic acid signaling pathway was higher in the WC than that in HL (TIFY10A, TIFY10B, TIFY11B, TIFY3B).

Studies have shown that in addition to MYC2 transcription factors[ 25, 26], JAZ proteins can bind directly or indirectly to transcription factors such as BHLH[ 27], MYB[ 28], AP2/ERF[ 29] and WRKY [ 30]. In Fraxinus mandshurica, there are three BHLH, three ERF, three WRKY transcription factors whose expression were consistent with JAZ expression. In addition to the regulation of transcription factors and hormones on the plants, the kinase also has a cascade of amplification functions in the growth and signaling process of the plant. In the significantly enriched KEGG metabolic pathway, twelve kinases whose transcription level in the WC provenance is significantly higher than that in HL provenance. Six of the twelve kinases belong to the serine/threonine protein kinase, two are leucine-rich protein kinases, and two are receptor-associated protein kinases on the cell wall.

Since growth and carbon storage of WC provenance seedlings are higher than that of HL provenance seedlings, we wanted to find genes with high expression of carbon metabolism-related genes in WC provenance seedlings. Three genes related to carbon metabolism, INVA, β-glucosidase and trehalose phosphate phosphatase J (TPPJ) were found. Invertase, was also known as sucrose invertase. The sucrose is broken down into fructose and glucose by invertase. Studies have shown that highly active invertases in plants are associated with rapid growth of tissues, such as in seedlings, young foliage, radicles and young fruits. The analysis showed that the transcription expression of sucrose invertase from WC provenance was significantly higher than that of HL provenance, indicating that the utilization rate of sucrose in WC was higher. Studies have shown that trehalose 6-phosphate (T6P) is a signal substance that can reflect the adequacy of carbon sources. In this study, it was found that the transcription expression level of T6P phosphorylase TPPJ in WC provenance was 5 times that of HL provenance, but there was no significant difference in the transcription expression of T6P metabolic enzyme TPS, indicating that trehalose content in WC provenance seedlings may be higher than HL provenance. It is speculated that the carbon use efficiency of WC provenance was higher than that of HL provenance.

In summary, plants produce a large amount of carbohydrates through photosynthesis. Sucrose, is the main form of carbon transport, whose content increase causes plants to produce higher levels of T6P. In the WC provenance, on the one hand, the high concentration of sucrose is decomposed by invertase, which is fully used for plant growth and development; on the other hand, the content of T6P is decreased to regulate growth by increasing the transcriptional expression of TPPJ. Due to the low utilization rate of sucrose, coupled with the inhibition of carbon source utilization by high concentration of T6P, the growth of HL provenance seedlings is slower than that of WC seedlings.

Identification of DEGs related to nitrogen deficiency

DEGs related to nitrogen deficiency have been identified in many species, such as duckweed[ 31], Magnaporthe oryzae[ 32], rice (Oryza sativa L.)[ 33], Tibetan wild barley[ 34], Maize[ 35], soybean[ 24], Arabidopsis[ 23, 36], etc. However, little is known about the DEGs related to nitrogen deficiency in Fraxinus mandshurica. In current study, analysis of DEGs in CP3, CP4 and CP5, CP6 can be used to identify genes that respond to nitrogen deficiency in Fraxinus mandshurica, and then the differences in response to nitrogen deficiency between two provenances can be further analyzed. Consequently, 8173 DEGs under nitrogen deficiency were identified using pair-wise comparison of each accession between normal and nitrogen deficiency in WC and HL. Different parts of plants respond differently to nutrient stress and the number of DEGs in roots was significantly higher than that in foliage under nitrogen deficiency. For example, DEGs in the CP4 (4225) were nearly 14 times as much as that in the CP3 (305) and DEGs in the CP5 (4531) were nearly 3.2 times as much as that in the CP6 (1414).

Many genes involved in nitrogen absorption and assimilation were differentially expressed under nitrogen deficiency relative to the normal condition. In plants, nitrogen is first actively absorbed by nitrate transporters in roots. In the current study, 10 DEGs encoding nitrate transporters were detected (Additional file 12), such as Unigene0034112, Unigene0050395, Unigene0046155, Unigene0047301, Unigene0018686, Unigene0022313, Unigene0026134, Unigene0050965, Unigene0039451 and Unigene0071953. Four, three, six, two and two DEGs encoding AMT, GS, GDH, NiR and NR were identified (Additional file 12), which were key enzymes in nitrogen assimilation. The transcription expression of genes encoding GS, NR and GOGAT was also down-regulated under nitrogen starvation in duckweed[ 31]. What’s more, macromolecular complex showed higher levels in CP3, CP4, CP5 and CP6, which suggested that macromolecular complex might play an important role in responding to nitrogen deficiency.

Nitrogen deficiency also has an important impact on carbon metabolism[ 31, 37] and leads to increased accumulation of carbohydrates. A total of 6 DEGs encoding CDF4, PYK1, RuBP and INVA were identified under nitrogen deficiency. Phenylalanine (Phe) plays an important role in the interconnection of plant primary metabolism and secondary metabolism. The metabolism of Phe plays a central role in the channeling of carbon from photosynthesis to the biosynthesis of phenylpropanoids and lignin is one of Phe-derived compounds. Lignin is a cross-linked polymer of three monolignols: p-hydroxyphenyl (H), guaiacyl (G), and syringyl (S) units[ 38]. Genes involved in monolignol biosynthesis and wood-related pathways are L-phenylalanine ammonia-lyase (PAL), cinnamate 4-hydroxylase (C4H), 4CL, cinnamyl alcohol dehydrogenase (CAD), hydroxycinnamoyl CoA:shikimate hydroxycinnamoyl transferase (HCT), 4-coumarate-3-hydroxylase (C3H), caffeoyl CoA 3-O-thyltransferase (CCoAOMT), Ferulate/coniferaldehyde 5-hydroxylase (F5H) and caffeic acid/5-hydroxyconiferaldehyde 3-O-methyltransferase (COMT)[ 38]. In the current study, the DEGs in the phenylpropanoid pathway for the biosynthesis of lignin (Fig9) were also significantly up-regulated such as the 4CL gene and lignin content is increased but the increase is not significant in roots under nitrogen deficiency (Additional file 13: Figure S7). Suppression of 4CL gene expression exerted the biggest impact on lignin production of all of the genetic manipulations of phenylpropanoid related genes in conifers[ 39] and 4CL RNAi lines exhibited a reduced lignin content of approximately 50%[ 40]. The lignin content of T 0 and T 1 generation plants over-expressed antisense 4CL gene was also reduced by 45.77% and 31.97% in Alfalfa (Medicago truncatula), respectively[ 41]. Induction of lignin biosynthesis is an adaptive response of plants subjected to many abiotic stresses[ 42]. What’s more, in maize shoots, N luxury significantly reduces lignin total content as well as the generation of S, H and G monomers of the lignin[ 43]. Transcriptional and physiological analyses identify a regulatory role for hydrogen peroxide in the lignin biosynthesis of copper-stressed rice roots[ 44]. Peroxidases are the main enzymes which are involved in the process of lignin biosynthesis[ 45]. What’s more, gibberellin (GA) is a positive regulator of lignification[ 38]. The key genes for lignin synthesis, 4CL and CAD1, are significantly induced by GA, which in turn changes the composition of the cell wall[ 46]. In this study, the DEGs in the terpenoids and polyketides metabolism pathway for the biosynthesis of abscisic acid (ABA) and gibberellin A4 (GA4) (Fig8) were also significantly up-regulated in roots under nitrogen deficiency. Compared with the control, the endogenous ABA content and GA content were significantly reduced and increased in the Fraxinus mandshurica roots by the nitrogen deficiency, respectively(Fig10). Meanwhile, transcription expression of many genes associated with absorption or translocation of other nutrients changed under nitrogen deficiency, such as phosphate (9), potassium (3), sulfate (5), Zinc (4), iron (3) and molybdate transporter (2), indicating that uptake of these nutrients in Fraxinus mandshurica is affected by nitrogen metabolism under cross-talking regulation.

Hormones, transcription factors and kinases are the three major regulators in plants. Studies have shown that mitogen-activated protein kinase kinase 9 plays an important role in plant nitrogen tolerance[ 47]. Protein phosphorylation may be induced by nitrogen deficiency in Fraxinus mandshurica. For example, a total of 14, 92, 132 and 152 DEGs encoding protein kinases were found in CP3, CP4, CP5 and CP6, respectively (Additional file 10). Transcriptional expression of nitrogen metabolism genes may be activated by transcription factors under nitrogen deficiency[ 21]. The results indicated that in the response to nitrogen deficiency, bHLH, MYB and AP2/ERF play a major role in roots, and WRKY plays a major role in foliage. The hormones homeostasis influences plant growth by regulating nutrients metabolism, such as proteins, nucleic acids and soluble carbohydrates[ 48-50]. Nitrogen deficiency caused a significant decrease in GA and IAA content but increase in ABA content in rice foliage[ 51]. In this study, Hormone (auxin, ABA, SA, JA and CK) pathways were predominatly induced by nitrogen deficiency, especially ABA and GA pathways (Fig.8). Compared to control, endogenous ABA content was reduced by nitrogen deficiency, but endogenous ZR, GA and IAA content was increased. Therefore, Fraxinus mandshurica is different from other species in response to nitrogen deficiency in endogenous hormones, and the reason for this difference needs further research. In addition, the response patterns of DEGs related to hormones under nitrogen deficiency in roots and foliage were different. For example, the hormones that predominantly played a role in the foliage were JA and ABA but in the roots were auxin, SA, CK and ABA.

Several nitrogen deficiency-related metabolic pathways were also identified by KEGG analysis[ 31]. Integrate CP3, CP4, CP5 and CP6 comparisons, the major pathways involved in the mechanism for responding to nitrogen deficiency were “Phenylpropanoid biosynthesis”, “Nitrogen metabolism”, “Diterpenoid biosynthesis”, “Isoquinoline alkaloid biosynthesis”. In addition, the pathway which had the most DEGs in the roots that responded most to nitrogen deficiency in both provenances was the “Ribosome” pathway. Studies on Magnaporthe oryzae proteomics have shown that nitrogen deficiency induces the synthesis of many extracellular proteins, and protein degradation and translation have also undergone extensive changes[ 32].

The transcriptional expression difference of two provenances in response to nitrogen deficiency

The number of DEGs of HL provenance seedlings was less under nitrogen stress (CP8) but significantly more in response to nitrogen deficiency than that of WC provenance seedlings, respectively. For example, the foliage of HL provenance responded to nitrogen deficiency more strongly than that of WC provenance. The numbers of KEGG pathways which were significantly enriched both in WC foliage (CP3) and roots (CP4) were less than that both in HL foliage (CP6) and roots (CP5), respectively. KEGG enrichment analysis showed that the response patterns of WC and HL to nitrogen deficiency were different both in the foliage and roots. The pathways with a larger number of genes enriched in WC foliage (CP3) were “Sulfur metabolism” and “Nitrogen metabolism”, however, that in HL foliage (CP6) were “Plant-pathogen interaction”, “Plant hormone signal transduction”, “Protein processing in endoplasmic reticulum”, “Phenylpropanoid biosynthesis” and “Starch and sucrose metabolism” (Fig6d). Similarly, the pathways with a larger number of genes enriched in WC roots (CP4) were “Ribosome”, “Oxidative phosphorylation”, “Phenylpropanoid biosynthesis”, “Cyanoamino acid metabolism” and “Diterpenoid biosynthesis”, however, that in HL roots (CP5) were “Ribosome”, “Starch and sucrose metabolism”, “Oxidative phosphorylation”, “Plant hormone signal transduction”, “Pentose and glucuronate interconversions”, “Tyrosine metabolism” and “Phenylpropanoid biosynthesis” (Fig.6c), which indicated that “Oxidative phosphorylation” and “Starch and sucrose metabolism” pathways were specifically enriched in the roots of WC and HL, respectively. The numbers of DEGs encoding NRTs, PK, TFs and hormone in HL were more than that in WC. The higher expression of NRTs could be attributed to its more nitrogen uptake and higher low nitrogen tolerance[ 34]. Studies have shown that trehalose can increase plant growth by increasing the expression of nitrogen metabolism genes under nitrogen deficiency conditions[ 52]. In plants, UDP-glucose and glucose-6-phosphate are catalyzed by trehalose 6-phosphate synthases (TPSs) to synthesize T6P, which is catalyzed by trehalose 6-phosphate phosphatases (TPPs) to form trehalose[ 53, 54]. In this study, one and four TPS genes transcription expression was up-regulated and down-regulated significantly under nitrogen deficiency in HL provenance seedlings roots, respectively. The endogenous ZR content in HL foliage was significantly higher than that in WC foliage under nitrogen deficiency (Fig10). The root lignin content in HL was higher than that in WC under nitrogen deficiency but not significant (Additional file 13: Figure S7). In summary, the response of HL provenance seedlings to nitrogen deficiency was stronger than that of WC provenance seedlings. However, under nitrogen deficiency, the number of up-regulated DEGs in the WC provenance seedlings is more than that in the HL (CP7 and CP8), especially in the roots (CP8), which indicated that the WC provenance seedlings were more tolerant to nitrogen deficiency than HL.

Conclusion

The data showed that high carbon content and high growth rate of WC provenance seedlings are mainly attributed to the high transcriptional expression of many metabolic genes in foliage under normal condition. The response of HL provenance seedlings to nitrogen deficiency is significantly greater than that of WC provenances seedlings, but the WC provenance seedlings were more tolerant to nitrogen deficiency than HL. Many genes related to nitrogen deficiency were identified, which will expand our current understanding of nitrogen responses.

Plant cultivation, nitrogen treatment and harvesting

According to the results about different carbon fixation Fraxinus mandshurica provenances research, we collected high carbon fixation (Wu Chang, WC) and low carbon fixation (Hai Lin, HL) provenances of seeds. After obtaining approval from the owner of the orchard, the seeds were obtained from the Qingshan Forest Seed Orchard of the Weihe Forestry Bureau and the Hailin Forestry Bureau of Heilongjiang Province and were named as WC and HL respectively. At May 20, the seeds were transplanted in plastic pots (800mm × 400mm) filled with soil, soil composition and proportion is volume (nutritional soil): volume (meteorite) = 3:1, and adjust the incubator to maintain the temperature and humidity at (23 ± 2) °C and 60% ~ 80%, adjust the fluorescent lamp to ensure sunshine 16 hours in length. The seedlings were transplanted in plastic pots (100mm diameter) filled with fine sand on July 20, a pot a tree. The sand is soaked with a solution of 1 mol per liter of hydrochloric acid for 4 hours and then washed with distilled water until the pH of the aqueous solution of sand is 6.5. Seedlings of two provenances were divided into three groups with 10 seedlings for each group and were provided with modified LA solution (0.5mM KCl, 0.9mM CaCl₂, 0.3mM MgSO₄, 0.6 mM KH₂PO₄, 42μM KH₂PO₄,10 μM Fe-EDTA, 2 μM MnSO₄, 10 μM H₃BO₃, 7 μM Na₂MoO₄, 0.2 μM ZnSO₄, 0.05μM C_OSO₄, and 0.2 μM CuSO₄) containing 0 (nitrogen deficiency) or 1000 (normal nitrogen) μM NH₄NO₃, respectively, and the nutrient solution was adjusted to pH 5.5. The LA solution was added every 2 days with 1000ml until August 20. The roots and foliage samples used for cDNA library construction were collected on 30 days after treatment, with three biological replicates per group. All samples were frozen in liquid nitrogen immediately upon collection and stored at -80°C.

cDNA library construction and sequencing

Total RNA was extracted by the Trizol method and Total RNA concentration was determined by Nanodrop 2000 (Thermo Fisher Scientific Nanodrop 2000). The RNA integrity number (RIN) by Agilent 2100 for assessing RNA quality. Total RNA was extracted from the 8 samples separately with 3 duplications, 24 RNA samples. After total RNA was extracted, mRNA was enriched by Oligo(dT) beads. Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcripted into cDNA with random primers. Second-strand cDNA were synthesized by DNA polymerase Ⅰ, RNase H, dNTP and buffer. Then the cDNA fragments were purified with QiaQuick PCR extraction kit, end repaired, poly(A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeq^TM 4000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

De novo cDNA assembly and functional annotation

To obtain clean reads, sequencing adaptors and low quality reads were removed from each library. Transcriptome de novo assembly was carried out with short reads assembling program – Trinity. Trinity is a modular method and software package which combines three components: Inchworm, Chrysalis and Butterfly. Firstly, Inchworm assembles reads by a greedy k-mer based approach, resulting in a collection of linear contigs. Next, Chrysalis clusters related contigs that correspond to portions of alternatively spliced transcripts or otherwise unique portions of paralogous genes, and then builds a de Bruijn graphs for each cluster of related contigs. Finally, Butterfly analyzes the paths taken by reads and read pairings in the context of the corresponding de Bruijn graph, and outputs one linear sequence for each alternatively spliced isoform and transcripts derived from paralogous genes. And the transcriptome feference database was obtained. All raw read data were deposited in the Genome Sequence Archive with the project ID. The unigene expression was calculated and normalized to RPKM (Reads Per kb per Million reads).

Basic annotation of unigenes includes protein functional annotation, pathway annotation, COG/KOG functional annotation and Gene Ontology (GO) annotation. To annotate the unigenes, we used BLASTx program (http://www.ncbi.nlm.nih.gov/BLAST/) with an E-value threshold of 1e−5 to NCBI non-redundant protein (Nr) database (http://www.ncbi.nlm.nih.gov), the Swiss-Prot protein database (http://www.expasy.ch/sprot), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg), and the COG/KOG database (http://www.ncbi.nlm.nih.gov/COG). Protein functional annotations could then be obtained according to the best alignment results.

Differential gene expression analysis

Each of the 24 Fraxinus mandshurica cDNA libraries was aligned separately to the transcriptome assemblies using Bowtie. The counting of alignments was estimated using the RSEM package. The Reads Per Kilobase of exon model per Million mapped reads (RPKM) method was then used to calculate the numbers of differentially expressed genes (DEGs) from pairwise comparisons of treatment group (nitrogen deficiency) Vs. control group (normal nitrogen) and WC Vs. HL. For the sake of brevity, WCKF, WCKR, HCKF, HCKR indicate cDNA libraries of the Wuchang provenance foliage, Wuchang provenance roots, Hailin provenance foliage, Hailin provenance roots under normal nitrogen, respectively. WTF, WTR, HTF, HTR indicate cDNA libraries of the Wuchang provenance foliage, Wuchang provenance roots, Hailin provenance foliage, Hailin provenance roots under nitrogen deficiency, respectively. From the biological replicates in this study, DESeq implemented in R package was used to analyze the differential expression between two groups. P values were adjusted using the Benjamini-Hochberg approach to control the false discovery rate (FDR). Genes with an adjusted P-value < 0.05 found by DESeq were usually identified as DEGs. However, genes with FDR < 0.01 and FC (Fold Change) ≥2 served as standards in the screening, and FC represents the specific value between two samples.

Quantitative real-time PCR analysis

Total RNA was extracted in triplicate from the eight samples of Fraxinus mandshurica (total of 24 samples) using Trizol® Reagent (Invitrogen, USA), and then reversed transcribed into cDNA using a PrimeScript® RT reagent Kit (Takara, Dalian, China). Nine DEGs associated with provenances seedlings of Fraxinus mandshurica and six DEGs between normal nitrogen and nitrogen deficiency were selected to validate the transcriptome data using quantitative real-time PCR (Additional file 1). The gene-specific primers for the nine unigene sequences were designed with Primer Premier 5 software, and the tubulin gene was used as an internal gene expression control: the gene was amplified with forward primer 5′- AGGACGCTGCCAACAACTTT -3′ and reverse primer 5′- TTGAGGGGAAGGGTAAATAGTG -3′. The amplifications were performed in 20 μl reactions consisting of 10 μl 2×TransStart® Green qPCR SuperMix (TransGen Biotech, Beijing, China), 0.8 μl of a mixture of the forward and reverse primers, 1.0 μl cDNA template, and 7.8 μl sterile dH₂O. The amplifications were performed on a ABI Prism7500 real-time PCR system using the following quantitative real-time PCR program: 95 °C for 30s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 34 s. Each amplification reaction was performed in triplicate.

Determination of phytohormone and lignin content

From each treatment, the fully expanded leaves were selected for the indoleacetic acid (IAA), gibberellin (GA), zeatinriboside (ZR) and abscisic acid (ABA) analyses. For the method of extracting and measuring phytohormone content, please refer to He[ 55]. The method of extracting and measuring lignin content is the same as phytohormone.

WC: Wuchang provenance; HL: Hailin provenance; DEGs: Differentially expressed genes; TFs: Transcription factors; PK: Protein kinase; NRTs: Nitrate transporters; AMTs: Ammonium transporters; NR: Nitrate reductase; NiR: Nitrite reductase; GS: Glutamine synthetase; GOGAT: Glutamate synthetase; SPS: sucrose phosphate synthase; SS: sucrose synthase; AI: acid invertase; RIN: RNA integrity number; RPKM: Reads Per kb per Million reads; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; WCKF: cDNA libraries of the WC foliage under normal nitrogen; WCKR: cDNA libraries of the WC roots under normal nitrogen; HCKF: cDNA libraries of the HL foliage under normal nitrogen; HCKR: cDNA libraries of the HL roots under normal nitrogen; WTF: cDNA libraries of the WC foliage under nitrogen deficiency; WTR: cDNA libraries of the WC roots under nitrogen deficiency; HTF: cDNA libraries of the HL foliage under nitrogen deficiency; HTR: cDNA libraries of the HL roots under nitrogen deficiency; FC: Fold Change; IAA: Indoleacetic acid; GA: Gibberellin; ZR: Zeatinriboside; ABA: Abscisic acid; CP1:WCKF vs. HCKF; CP2: WCKR vs. HCKR; CP3: WCKF vs. WCTF; CP4: WCKR vs. WTR; CP5: HCKR vs. HTR; CP6: HCKF vs. HTF; CP7: WTF vs. HTF; CP8: WTR vs. HTR; JA: Jasmonic acid; SA: Salicylic acid; ZIM: Zinc-finger protein; TIFY: ZIM precursor; RLKs: receptor-like serine/threonine-protein kinase; CIPKs: CBL-interacting serine/threonine-protein kinase 14-like; STPKs: Serine/threonine-protein kinase; CDF4: Dof zinc finger protein DOF1.5-like; PYK1: Pyruvate kinase; RuBP: Ribulose-1,5-bisphosphate carboxylase; INVA: Neutral invertase A; SAMT: salicylic acid methyltransferase; T6P: trehalose 6-phosphate; TPSs: trehalose 6-phosphate synthases; TPPs: trehalose 6-phosphate phosphatases; TPPJ: Trehalose phosphate phosphatase J; PAL: L-phenylalanine ammonia-lyase; C4H: cinnamate 4-hydroxylase; 4CL: 4-coumarate--CoA ligase; CAD: cinnamyl alcohol dehydrogenase; HCT: hydroxycinnamoyl CoA:shikimate hydroxycinnamoyl transferase; C3H: 4-coumarate-3-hydroxylase; CCoAOMT: caffeoyl CoA 3-O-thyltransferase; F5H: Ferulate/coniferaldehyde 5-hydroxylase; COMT: caffeic acid/5-hydroxyconiferaldehyde 3-O-methyltransferase.

Acknowledgements

We thank the OmicShare forum, a free online platform for data analysis (www.omicshare.com/tools).

Funding

This work was supported by Heilongjiang Province Applied Technology Research and Development Program (GA19B201). The funders had no role in study design, data collection, analysis and interpretation, or preparation of the manuscript.

Availability of data and materials

The data sets supporting the results of this article are included within the article and as Additional files.

Authors’ contribution

X T Z, Y Z and F Z: Contributed equally in objective formulation, experimental design and set up data collection, analysis and interpretation. X T Z, L Y, X Z and Z L: Contributed equally in all tissue collection and qRT-PCR analysis. T Z, H L, X W and Y X: Contributed equally in seed collection and dormancy release processing. All authors have read and approved the manuscript.

Competing interests

The authors declare that they have no competing interests.

Author details

¹National Key Laboratory of Tree Genetics and Breeding, Northeast Forestry University, Xiangfang Distract, Hexing road No.26, Harbin, Heilongjiang150040, People’s Republic of China.

² Life Science College, Northeast Forestry University, Xiangfang Distract, Hexing road No.26, Harbin, Heilongjiang150040, People’s Republic of China.

³Heilongjiang Forestry and Grassland Bureau, Xiangfang Distract, Hengshan road No.10, Harbin, Heilongjiang150090, People’s Republic of China.

⁴Shanhetun Forestry Bureau of Heilongjiang Province, Shanhetun Forestry Areas, Wuchang City, Harbin, Heilongjiang150232, People’s Republic of China.

Ethics approval and consent to participate

Not applicabole.

Consent for publication.

Not applicabole.

1. Nabais C, Hansen JK, David-Schwartz R, Klisz M, López R, Rozenberg P: The effect of climate on wood density: What provenance trials tell us? Forest Ecology and Management 2018, 408:148-156.

2. Nieto V, Giraldo-Charria D, Sarmiento M, Borralho N: Effects of provenance and genetic variation on the growth and stem formation of Eucalyptus pellita in Colombia. Journal of Tropical Forest Science 2016, 28(3):227-234.

3. Chhin S, Zalesny RS, Parker WC, Brissette J: Dendroclimatic analysis of white pine (Pinus strobus L.) using long-term provenance test sites across eastern North America. Forest Ecosystems 2018, 5(1).

4. Kleiber A, Duan Q, Jansen K, Verena Junker L, Kammerer B, Rennenberg H, Ensminger I, Gessler A, Kreuzwieser J: Drought effects on root and needle terpenoid content of a coastal and an interior Douglas fir provenance. Tree physiology 2017, 37(12):1648-1658.

5. Skulason B, Hansen OK, Nielsen UB: Provenance variation in Christmas tree characteristics in subalpine fir (Abies lasiocarpa) planted in Denmark and Iceland. Scandinavian Journal of Forest Research 2017, 33(3):222-232.

6. Václav B, Petr N, Jaroslav D: Results of Czech ash provenance experiment. Journal of Forest Science 2017, 63(No. 6):263-274.

7. Pan W, Sun J, Xia X, Liu R, Wu X, Li Y: Provenance variation in growth, stem-form, and wood basic density of 24-year-old Liriodendron. AUSTRIAN JOURNAL OF FOREST SCIENCE 2018, 135(4):343-362.

8. Zhou CB: Effect of Provenance and Climate on Xylem Anatomy of Haloxylon Ammodendron (C. A. Mey) Bunge in the Gurbantunggut Desert, China. Applied Ecology and Environmental Research 2017, 15(3):1309-1321.

9. Lu P, Parker WC, Colombo SJ, Man R, Baeten L: Restructuring tree provenance test data to conform to reciprocal transplant experiments for detecting local adaptation. Journal of Applied Ecology 2016, 53(4):1088-1097.

10. Thomson AM, Riddell CL, Parker WH: Boreal forest provenance tests used to predict optimal growth and response to climate change: 2. Black spruce. Canadian Journal of Forest Research 2009, 39(1):143-153.

11. Hebel I, Haas R, Dounavi A: Genetic Variation of Common Ash (Fraxinus excelsior L.) Populations From Provenance Regions in Southern Germany by Using Nuclear and Chloroplast Microsatellites. Silvae Genetica 2006, 55(1-6):38-44.

12. Santoni I, Callone E, Sandak A, Sandak J, Dire S: Solid state NMR and IR characterization of wood polymer structure in relation to tree provenance. Carbohydrate polymers 2015, 117:710-721.

13. Calvo L, Hernández V, Valbuena L, Taboada A: Provenance and seed mass determine seed tolerance to high temperatures associated to forest fires in Pinus pinaster. Annals of Forest Science 2015, 73(2):381-391.

14. Spinnler D, Egli P, Körner C: Provenance effects and allometry in beech and spruce under elevated CO2 and nitrogen on two different forest soils. Basic and Applied Ecology 2003, 4(5):467-478.

15. Klisz M, Ukalska J, Koprowski M, Tereba A, Puchałka R, Przybylski P, Jastrzębowski S, Nabais C: Effect of provenance and climate on intra-annual density fluctuations of Norway spruce Picea abies (L.) Karst. in Poland. Agricultural and Forest Meteorology 2019, 269-270:145-156.

16. Amusant N, Digeon A, Descroix L, BruneAu O, BezArD V, B eAuchène J: Planting rosewood for the sustainable essential oil production:Influence of surrounding forest and seed provenance on tree growth and essential oil yield. Bois et Forêts des Tropiques 2015(326):57-65.

17. Klisz M, Buras A, Sass-Klaassen U, Puchalka R, Koprowski M, Ukalska J: Limitations at the Limit? Diminishing of Genetic Effects in Norway Spruce Provenance Trials. Frontiers in plant science 2019, 10:306.

18. Montwe D, Isaac-Renton M, Hamann A, Spiecker H: Cold adaptation recorded in tree rings highlights risks associated with climate change and assisted migration. Nature communications 2018, 9(1):1574.

19. Thomson AM, Parker WH: Boreal forest provenance tests used to predict optimal growth and response to climate change. 1. Jack pine. Canadian Journal of Forest Research 2008, 38(1):157-170.

20. Santiso X, Retuerto R: Low among-provenance differences in structural and functional plasticity in response to nutrients in saplings of the circum-Mediterranean tree Arbutus unedo L. Tree physiology 2015, 35(10):1118-1128.

21. Jiang L, Ball G, Hodgman C, Coules A, Zhao H, Lu C: Analysis of Gene Regulatory Networks of Maize in Response to Nitrogen. Genes 2018, 9(3).

22. Wang Q, Nian J, Xie X, Yu H, Zhang J, Bai J, Dong G, Hu J, Bai B, Chen L et al: Genetic variations in ARE1 mediate grain yield by modulating nitrogen utilization in rice. Nature communications 2018, 9(1):735.

23. Bi YM, Wang RL, Zhu T, Rothstein SJ: Global transcription profiling reveals differential responses to chronic nitrogen stress and putative nitrogen regulatory components in Arabidopsis. BMC genomics 2007, 8:281.

24. Hao QN, Zhou XA, Ai HS, Wang C, Zhou R, Chen SL: Identification of genes associated with nitrogen-use efficiency by genome-wide transcriptional analysis of two soybean genotypes. BMC genomics 2011, 12(1):525.

25. Chini A, Boter M, Solano R: Plant oxylipins: COI1/JAZs/MYC2 as the core jasmonic acid-signalling module. The FEBS journal 2009, 276(17):4682-4692.

26. Lorenzo O, Chico JM, Sanchez-Serrano JJ, Solano R: JASMONATE-INSENSITIVE1 encodes a MYC transcription factor essential to discriminate between different jasmonate-regulated defense responses in Arabidopsis. The Plant cell 2004, 16(7):1938-1950.

27. Chung HS, Niu Y, Browse J, Howe GA: Top hits in contemporary JAZ: an update on jasmonate signaling. Phytochemistry 2009, 70(13-14):1547-1559.

28. Wang W, Liu G, Niu H, Timko MP, Zhang H: The F-box protein COI1 functions upstream of MYB305 to regulate primary carbohydrate metabolism in tobacco (Nicotiana tabacum L. cv. TN90). Journal of experimental botany 2014, 65(8):2147-2160.

29. Zhang X, Zhu Z, An F, Hao D, Li P, Song J, Yi C, Guo H: Jasmonate-activated MYC2 represses ETHYLENE INSENSITIVE3 activity to antagonize ethylene-promoted apical hook formation in Arabidopsis. The Plant cell 2014, 26(3):1105-1117.

30. Qiu D, Xiao J, Ding X, Xiong M, Cai M, Cao Y, Li X, Xu C, Wang S: OsWRKY13 Mediates Rice Disease Resistance by Regulating Defense-Related Genes in Salicylate- and Jasmonate-Dependent Signaling. Molecular Plant-Microbe Interactions 2007, 20(5):492-499.

31. Yu C, Zhao X, Qi G, Bai Z, Wang Y, Wang S, Ma Y, Liu Q, Hu R, Zhou G: Integrated analysis of transcriptome and metabolites reveals an essential role of metabolic flux in starch accumulation under nitrogen starvation in duckweed. Biotechnology for biofuels 2017, 10:167.

32. Oh Y, Robertson SL, Parker J, Muddiman DC, Dean RA: Comparative proteomic analysis between nitrogen supplemented and starved conditions in Magnaporthe oryzae. Proteome science 2017, 15:20.

33. Yang W, Yoon J, Choi H, Fan Y, Chen R, An G: Transcriptome analysis of nitrogen-starvation-responsive genes in rice. BMC plant biology 2015, 15:31.

34. Quan X, Zeng J, Ye L, Chen G, Han Z, Shah JM, Zhang G: Transcriptome profiling analysis for two Tibetan wild barley genotypes in responses to low nitrogen. BMC plant biology 2016, 16:30.

35. Bi YM, Meyer A, Downs GS, Shi X, EI-kereamy A, Lukens L, Rothstein SJ: High throughput RNA sequencing of a hybrid maize and its parents shows different mechanisms responsive to nitrogen limitation. BMC genomics 2014, 15(1):77.

36. Gutierrez RA, Lejay LV, Dean A, Chiaromonte F, Shasha DE, Coruzzi GM: Qualitative network models and genome-wide expression data define carbon/nitrogen-responsive molecular machines in Arabidopsis. Genome Biol 2007, 8(1):R7.

37. Ning P, Yang L, Li C, Fritschi FB: Post-silking carbon partitioning under nitrogen deficiency revealed sink limitation of grain yield in maize. Journal of experimental botany 2018, 69(7):1707-1719.

38. Zhao Q: Lignification: Flexibility, Biosynthesis and Regulation. Trends in plant science 2016, 21(8):713-721.

39. Pascual MB, El-Azaz J, de la Torre FN, Canas RA, Avila C, Canovas FM: Biosynthesis and Metabolic Fate of Phenylalanine in Conifers. Frontiers in plant science 2016, 7:1030.

40. Wagner A, Donaldson L, Kim H, Phillips L, Flint H, Steward D, Torr K, Koch G, Schmitt U, Ralph J: Suppression of 4-coumarate-CoA ligase in the coniferous gymnosperm Pinus radiata. Plant physiology 2009, 149(1):370-383.

41. Meng J, Li C, Zhao M, Wang C, Ru Y, Cui Z, Han Y: Lignin biosynthesis regulated by the antisense 4CL gene in alfalfa. Czech Journal of Genetics and Plant Breeding 2018, 54(No. 1):26-29.

42. Moura JC, Bonine CA, de Oliveira Fernandes Viana J, Dornelas MC, Mazzafera P: Abiotic and biotic stresses and changes in the lignin content and composition in plants. Journal of integrative plant biology 2010, 52(4):360-376.

43. Sun Q, Liu X, Yang J, Liu W, Du Q, Wang H, Fu C, Li WX: MicroRNA528 Affects Lodging Resistance of Maize by Regulating Lignin Biosynthesis under Nitrogen-Luxury Conditions. Molecular plant 2018, 11(6):806-814.

44. Liu Q, Zheng L, He F, Zhao F-J, Shen Z, Zheng L: Transcriptional and physiological analyses identify a regulatory role for hydrogen peroxide in the lignin biosynthesis of copper-stressed rice roots. Plant and Soil 2014, 387(1-2):323-336.

45. Herrero J, Esteban-Carrasco A, Zapata JM: Looking for Arabidopsis thaliana peroxidases involved in lignin biosynthesis. Plant physiology and biochemistry : PPB 2013, 67:77-86.

46. Garcia-Rojas M, Meneses M, Oviedo K, Carrasco C, Defilippi B, Gonzalez-Aguero M, Leon G, Hinrichsen P: Exogenous gibberellic acid application induces the overexpression of key genes for pedicel lignification and an increase in berry drop in table grape. Plant physiology and biochemistry : PPB 2018, 126:32-38.

47. Luo J, Wang X, Feng L, Li Y, He JX: The mitogen-activated protein kinase kinase 9 (MKK9) modulates nitrogen acquisition and anthocyanin accumulation under nitrogen-limiting condition in Arabidopsis. Biochemical and biophysical research communications 2017, 487(3):539-544.

48. Hisano H, Matsuura T, Mori IC, Yamane M, Sato K: Endogenous hormone levels affect the regeneration ability of callus derived from different organs in barley. Plant physiology and biochemistry : PPB 2016, 99:66-72.

49. Verbon EH, Liberman LM: Beneficial Microbes Affect Endogenous Mechanisms Controlling Root Development. Trends in plant science 2016, 21(3):218-229.

50. Albacete AA, Martinez-Andujar C, Perez-Alfocea F: Hormonal and metabolic regulation of source-sink relations under salinity and drought: from plant survival to crop yield stability. Biotechnology advances 2014, 32(1):12-30.

51. Xiong Q, Tang G, Zhong L, He H, Chen X: Response to Nitrogen Deficiency and Compensation on Physiological Characteristics, Yield Formation, and Nitrogen Utilization of Rice. Frontiers in plant science 2018, 9.

52. Lin Y, Zhang J, Gao W, Chen Y, Li H, Lawlor DW, Paul MJ, Pan W: Exogenous trehalose improves growth under limiting nitrogen through upregulation of nitrogen metabolism. BMC plant biology 2017, 17(1):247.

53. Lunn JE, Delorge I, Figueroa CM, Van Dijck P, Stitt M: Trehalose metabolism in plants. The Plant journal : for cell and molecular biology 2014, 79(4):544-567.

54. Goddijn OJM, van DK: Trehalose metabolism in plants. Trends in plant science 1999, 4(8):315-319.

55. He Z, Zhan Y, Zeng F, Zhao X, Wang X: Drought physiology and gene expression characteristics of Fraxinus interspecific hybrids. Plant Growth Regulation 2015, 78(2):179-193.

Table1 Differentially expressed genes (DEGs) identified in pairwise comparisons

Code	Stage comparison	Total DEGs	Up-regulated	Down-regulated
CP1	WCKF vs. HCKF	669	285	384
CP2	WCKR vs. HCKR	149	77	72
CP3	WCKF vs. WCTF	305	165	140
CP4	WCKR vs. WTR	4225	1907	2318
CP5	HCKR vs. HTR	4531	2386	2145
CP6	HCKF vs. HTF	1414	991	423
CP7	WTF vs. HTF	141	67	74
CP8	WTR vs. HTR	178	44	134

WCKF, WCKR, HCKF, HCKR indicate cDNA libraries of the Wuchang provenance foliage, Wuchang provenance roots, Hailin provenance foliage, Hailin provenance roots under normal nitrogen, respectively. WTF, WTR, HTF, HTR indicate cDNA libraries of the Wuchang provenance foliage, Wuchang provenance roots, Hailin provenance foliage, Hailin provenance roots under nitrogen deficiency, respectively.

Additional file1: Real-time quantitative PCR primers. (xlsx 11.0kb)

Additional file2: RNA quality assessment of 24 transcripts of Fraxinus mandshurica. (xlsx 11.1kb)

Additional file3: The most reliable significantly enriched pathways of DEGs. (xlsx 15.4kb)

Additional file4: Differentially expressed genes (DEGs) in HL and WC. (pdf 161kb)

Additional file5: GO annotation of DEGs. (xlsx 25.6kb)

Additional file6: KEGG pathway analysis for DEGs. (xlsx 342kb)

Additional file7: DEGs related to hormone signaling. (xlsx 44.7kb)

Additional file8: DEGs related to metabolism. (xlsx 115kb)

Additional file9: Differentially expressed transcription factors between two Fraxinus mandshurica provenances. (pdf 457kb)

Additional file10: DEGs related to RLKs. (xlsx 42.9kb)

Additional file11: DEGs related to transcription factors. (xlsx 34.1kb)

Additional file12: DEGs encoding protein transporters and nitrate assimilation enzymes. (xlsx 41.7kb)

Additional file13: The content of lignin in roots of HL and WC under normal nitrogen and nitrogen deficiency. (pdf 207kb)

Transcriptome analysis for Fraxinus mandshurica seedlings from different carbon fixation and growth provenances in response to nitrogen deficiency

Status:

Version 1

Abstract

Figures