3.1 Plant growth and root morphology
Seedling height of Pi-starved Chinese fir decreased while root elongation increased compared with that of Pi-sufficient Chinese fir seedlings (Fig. 1 A-C). After 30days of growth under the experimental conditions, the height growth of Chinese fir seedlings under Pi deficiency condition was significantly lower than that of Chinese fir seedlings grown under high Pi supply. The root growth of Chinese fir under Pi deficiency condition was significantly higher than that under normal Pi supply condition. The root growth increased with the increase of Pi deficiency stress intensity and stress time. The root/shoot ratio of Chinese fir seedlings under Pi deficiency was higher than that under normal condition. Under the experimental conditions, Pi deficiency inhibited seedling height growth and promoted root elongation compared with Pi sufficient conditions.
After 30 and 60 days of growth under different treatment conditions, phosphorus content in roots decreased significantly in Pi-deficient seedlings compared with Pi-sufficient seedlings (Fig. 1-D). The phosphorus concentrations in the shoot and root of plants grown without phosphorus at 30 days of stress were only 78% and 70% of those grown under adequate phosphorus condition, respectively, and further decreased to 57% and 67% at 60 days.
Analysis of other nutrient elements revealed that Pi deficiency inhibited the accumulation of Mg in shoot, increased the accumulation of Ca and Mn, but no significant change in the accumulation of K, Mg, Cu and Zn(Fig. 1-E).After Pi deficiency treatment, the changes of element contents in roots were different from those in shoot. Pi deficiency increased the concentrations of K and Cu in roots, decreased the concentrations of Fe, and the contents of Ca, Mg and Mn had no significant difference. Under Pi deficiency, the ratios of root/shoot element concentrations of P, K, Mg and Cu increased while that of Ca, Fe and Mn decreased.
3.2 Functional annotation and classification
Transcriptome sequencing data were deposited in the NCBI SRA database under the accession number PRJNA681878. The assembled Unigene sequences were compared with the Nr, SwissProt, KOG and KEGG databases using blastx (Fig. S1-A). A total of 50,808 Unigenes were annotated out of 126,413 Unigenes, and 45,445, 39,381, 34,957 and 21,222 Unigenes were annotated out of the four databases, accounting for 35.95%, 31.15%, 27.65% and 16.79% of total Unigenes respectively. Using blastx and Nr databases to do species distribution statistics (Fig. S1-B), the results showed that, there were high homology between Chinese fir and Anthuriumamnicola (7.66%), Amborellatrichopoda (7.64%), Nelumbonucifera (5.32%), Piceasitchensis (2.68%) and Elaeisguineensis (2.04%).
According to the Nr annotation information, the Gene Ontology (GO) annotation information of Unigenes was obtained by Blast2go, and all Unigenes were classified by WEGO, and the functional distribution characteristics of Unigenes were obtained (Fig. S1-C).The most represented GO subcategories within the biological process main term were ‘cellular and metabolic processes’ and ‘cellular process’, and the cellular component main term were mainly represented by the term ‘cell or cell part’ and ‘membrane or membrane part’. As for molecular function main term, 6,196 transcripts corresponded to ‘catalytic activity’ and 4,655 sequences to ‘binding category’. KOG classifies gene sequences into homology (Fig. S1-D). The KOG database contains 34,957 annotation results, grouped into 25 different biological functions. The ‘General function prediction only’ term contained the largest number of annotated Unigenes(10,894, or 18.96%), followed by ‘Posttranslational modification, protein turnover, chaperones’ term(7,431, or 12.93%), and ‘Signal transduction mechanisms’ term(6,762, or 11.77%).
In the KEGG annotated results (Fig. 2), 21,222 Unigenes were annotated into 135 different metabolic pathways in five main categories and 19 subcategories. Unigenes were most involved in ‘Metabolic pathways’ (4,221 Unigenes), accounting for 19.89% of KEGG annotated Unigenes, and 2,246 Unigenes (10.58%) were annotated into ‘Biosynthesis of secondary metabolites’ term. Sucrose signaling, hormone synthesis, organic acid metabolism and phosphorus transport induced by Pi deficiency were closely related to starch and sucrose metabolism, plant hormone signal transduction, tricarboxylic acid cycle (TCA cycle) and ABC transporter in KEGG metabolic pathway, and 369, 307, 228 and 182 Unigenes were involved in these terms, respectively.
3.3 Screening of differentially expressed genes
Significant differentially expressed genes were analyzed among the subgroups (Fig. S2).After 3 days of stress treatment, 204 differentially expressed genes were up-regulated and 318 differentially expressed genes were down-regulated in low Pi group compared with normal Pi supply group. Compared with high Pi group, 86 differential expressed genes were up-regulated and 920 differential genes were down-regulated in Pi starvation group while 82 differentially expressed genes up-regulated and 711 differentially expressed genes down-regulated compared with low Pi group. After 30 days of Pi deficiency treatment, 283 differentially expressed genes were up-regulated and 1,665 differentially expressed genes were down-regulated in low Pi group compared with the normal Pi supply group, and 471 differentially expressed genes were up-regulated in Pi deficiency group compared with normal phosphorus supply group. There were 1095 differentially expressed genes (838 up-regulated and 508 down-regulated) in Pi starvation group and low Pi group.
From different subgroups, the number of down-regulated genes was more than up-regulated genes under Pi deficiency. Compared with HP treatment, HP-VS-NP subgroup had more differentially expressed genes than HP-VS-LP subgroup under the same treatment time. The number of differentially expressed genes in 30-day treatment was higher than that in 3-day treatment. The results of gene ontology categories with significantly differentially expressed gene among groups showed that the root system responded to Pi deficiency by changing its cell components, regulating its catalytic activity and adapting to the stress environment through various metabolic pathways. The KEGG pathway annotation of differentially expressed genes under Pi deficiency showed that the accumulation of oxidative phosphorylation, transmembrane transporters, osmotic stress response, hormone signal transduction, pyruvate metabolism, tricarboxylic acid cycle, acid phosphatase synthesis and other metabolic pathways were closely related to the responses of Chinese fir to Pi deficiency.
The gene expression patterns were clustered according to the characteristics of three consecutive samples, and the genes with certain biological characteristics were selected from the results of the clustering. The trends of differentially expressed genes in the 3-day HP-LP-NP subgroup and 30-day HP-LP-NP subgroup were analyzed. There were 1,829 tendency genes in 3-day HP-LP-NP subgroup (Fig. 3-A). Among them, 42 genes showed an upward trend of gene expression with the increase of Pi deficiency (Fig. 3-A, profile7), and 4 Unigenes obtained annotation in GO and KEGG. The results showed that these 4 up-regulated genes were mainly involved in the phosphorylated carbohydrate metabolism of Chinese fir. There were 315 genes (Fig. 3-A, profile0) whose expression decreased with the increase of Pi deficiency. The GO annotation results were distributed in 20 subcategories. In the molecular function category, ‘catalytic activity’ term and ‘binding’ term enriched the most genes, and in the biological process category, ‘cellular’ process and ‘metabolic’ process accounted for the most. The top five pathways for KEGG annotation were microbial metabolism, lysine degradation, tricarboxylic acid cycle (TCA) cycle, carbon metabolism and ribosome in different environments. There were 3,257 tendency genes in 30-day HP-LP-NP subgroup (Fig. 3-B). Among them, there were 259 genes whose expression increased with the increase of Pi deficiency (Fig. 3-B, profile7).There were 32 genes that got the GO annotation, participating in ‘metabolic process’ term, ‘cellular process’ term, ‘cell’ term, ‘cell part’ term, ‘catalytic activity’ term and ‘binding term’ respectively. KEGG annotated 33 genes of ascending tendency, and the pathway that enriched the most genes was biosynthesis and metabolic pathway of secondary metabolites. The trend analysis showed that there were 318 genes with the characteristic of decreasing expression with increased degree of Pi deficiency (Fig. 3-B, profile0). There were 36 genes involved in the ‘metabolic process’, ‘cellular’ term, ‘cell’ term, ‘cell part’ term, ‘membrane’ term, ‘catalytic activity’ term. KEGG annotated 36 genes, including lipoic acid metabolism, phosphatidylinositol metabolism, starch and sucrose metabolism, tricarboxylic acid cycle, RNA degradation, protein export, RNA transport and other metabolic pathways.
3.4 Genes putatively involved in Pi deficiency response in Chinese fir root
When plants grow in phosphorus deficient conditions, the roots usually activate insoluble phosphate in the soil by secreting organic acids, and utilize the organic phosphorus around the rhizosphere by secreting acid phosphatase (ACP). Purple acid phosphatase (PAP) is a special class of acid phosphatase, and its members have different responses to Pi deficiency and their functions are regulated by transcription factors and SPX protein. SPX is an important regulator of phosphorus signaling network, which can interact with PHR1 to regulate the expression of phosphate transporter gene PHT1, thus promoting phosphorus uptake by roots under Pi deficiency. A total of 18 ACP genes, 70 PAP genes, 10 SPX structure and 10 protein genes were screened from the transcriptome data based on all Unigenes annotation in GO and KEGG databases. Among them, 17 PAP genes, three ACP genes and three SPX domain proteins were differentially expressed under Pi deficiency. The differentially expressed PAP, ACP, and SPX genes were analyzed by expression thermography (Fig. S3-A), and the results showed that the expression of Unigene074730 in SPX gene increased and the expression of Unigene074730 in SPX gene decreased. Except ACP gene Unigene0100264 and PAP gene Unigene0013971 and Unigene104552, the expression of other ACP and PAP genes increased with the increase of Pi deficiency degree and the extension of stress time. The results showed that Chinese fir had the physiological response of synthesizing ACP to activate the organic phosphorus in its rhizosphere under Pi deficiency.
According to the annotation of GO and KEGG database, 20 differentially expressed genes related to organic acid synthesis and metabolism were screened. The gene expression of these 20 genes under different Pi deficiency was analyzed by thermography (Fig. S3-B), and the results showed that the differentially expressed genes in the process of organic acid synthesis and metabolism of Chinese fir under Pi deficiency were up-regulated upon 30days of exposure to Pi deficiency. However, the up-regulation of genes related to organic acid synthesis and metabolism was not obvious under 3-day Pi deficiency treatment. This might be related to Chinese fir root activating soil phosphorus by synthesizing organic acids after30-dayPi deficiency, but the functions of genes need to be further studied.
Based on all Unigene annotation results, a total of 86 Unigene entries annotated as phosphate transporter family and phosphate protein carrier were retrieved. A total of 35 Unigene entries were confirmed to be involved in inorganic phosphate transport by sequence signature analysis compared to BLAST analysis. The other 51 genes may be involved in phosphorus transport, the function of which is unknown. Among the 35 phosphorus transporter genes, six genes, which had high homology with PHO1, were annotated as PHO1 family genes by Nr and Swissprot and classified as small molecule transporter by KOG. We found 16 genes with high homology with PHT1 family, inorganic phosphorus transporters; one carrier family with high homology to the PHT2 gene, sodium-dependent phosphate transporter; 12 carrier proteins annotated as mitochondrial phosphate, which were highly homologous and predicted to PHT3 family. Based on the analysis of the characteristics of the PHT family genes, and the prediction of the coding sequence (CDS) and the characteristics of the PHT gene multiple sequence alignment with Arabidopsis thaliana, rice and poplar, four genes with high homology to PHT4 gene were screened out, which were highly similar to the known PHT family genes. The Unigene annotation messages contained three probably anion transporters, whose pathway were annotated as perfuse of the major facilitator superfamily and one possibly organic anion transporter (Unigene 0074757). And two PHT5 genes were annotated with membrane proteins containing SPX domains. A total of 41 genes involved in inorganic phosphate transport were obtained.
Thermographic analysis of seven Unigenes with significant differential expression of the PHT gene family (Fig. 4-A) showed that Unigene0035030 (PREDICTED: PHO) was overexpressed after 30 days of growth in Pi deficient condition, but the expression level decreased with the increase of the extent of Pi deficiency. Unigene0006988 (PREDICTED: PHT1) was enhanced after 30 days of Pi deficiency. The expression of Unigene0010395 (PREDICTED: PHT1) increased with the increase of Pi deficiency time, and decreased with the increase of Pi deficiency. The expression of Unigene0048799 (PREDICTED: PHT1) increased with the increase of Pi deficiency time. The Unigene0093392 (PREDICTED: PHT3) increased with the increase of Pi deficiency time while the expression of Unigene0044771 (PREDICTED: PHT3) was enhanced after 30 days of Pi deficiency. The expression of Unigene0067125 (PREDICTED: PHT4) and Unigene0074757 (PREDICTED: PHT4) increased with the decrease of phosphorus supply, the expression level was higher after 3 days than after 30 days.
In order to better understand the role of PHT gene in phosphate uptake by Chinese fir root, and clarify the expression pattern of PHT gene in Chinese fir root under different phosphorus supply treatments, seven known up-regulated PHT gene families expressed in phosphate starvation response were selected and verified by qRT-PCR (Fig. 4-B) and fluorescence quantitative analysis. The expression of Unigene0006988 (PREDICTED: PHT1) increased with the increase of stress time, and significantly increased under the condition of Pi deficiency for 30 days. Unigene0010395 (PREDICTED: PHT1) gene expression increased with the increase of Pi stress time, but the expression of Pi starvation group was lower than low Pi group at the same P deficiency time. The expression of Unigene0048799 (PREDICTED: PHT1) increased with the increase of Pi deficiency time, and was significantly higher after 30 days than after 3 days. At the same time, the higher the stress degree, the higher the expression amount, which was different from the result of gene expression heat map obtained by transcriptome sequencing. The expression level of Unigene0093392 (PREDICTED: PHT3) in different phosphorus supply groups was not significantly different after 3 days, but increased with the time of stress. The expression of high Pi treatment group was lower than low Pi and Pi starvation groups after 30-dayof growth. The expression of Unigene0044771 (PREDICTED: PHT3) increased with the increase of Pi deficiency time and increased under the condition of Pi deficiency after 30 days. The expression of Unigene0067125 (PREDICTED: PHT4) increased with the increase of Pi deficiency time, and increased significantly after 30-day treatment. The expression level of Unigene0074757 (PREDICTED: PHT4) was higher under the condition of Pi deficiency than that of normal phosphate treatment, and the expression level of Unigene0074757 increased noticeably after 30 days of Pi deficiency. The expression pattern of the rest of the genes was similar to that of RNA-seq.
3.5 Citric acid and glyoxylate cycle pathway analysis associated with root-released organic acids
Citric acid and glyoxylate cycle play an important role in the synthesis of organic acids in plant tissues, and the transcriptome sequencing of Chinese fir roots under phosphorus starvation revealed that the main metabolic pathway for enrichment of differentially expressed genes was the biosynthetic pathway of organic acids. The annotated transcripts related to citric acid and glyoxylate cycle were mapped, and28 genes with significant differences in enzyme coding presumably involved citric acid and glyoxylate cycle acid were identified under Pi starvation (Fig.5-A). Differential gene expression thermograms showed that 11 genes were up-regulated, 12 genes were down-regulated and 5 genes had no specific expression pattern (Fig.5-B). In general, the metabolism of pyruvic acid of Chinese fir was enhanced under Pi deficiency, and 5 pyruvic dehydrogenase genes were up-regulated to promote the synthesis of citric acid, at the same time, the expression of isocitratelyase Unigene0025798 in glyoxylic acid synthesis pathway was also increased to promote glyoxylic acid synthesis. The expression of malate synthase Unigene003021 was up-regulated to promote malic acid synthesis under Pi deficiency, therefore, the synthesis of citric acid, malic acid, glyoxylic acid and other organic acids could be promoted under Pi deficiency.
In addition, qRT-PCR analysis was used to study the expression of transcripts in citric acid and glyoxylate cycle pathway, including pyruvate dehydrogenase, citrate synthase, aconitase, isocitrate dehydrogenase and isocitrate lyase (Fig.5-C). The results of expression verification were consistent with the main expression characteristics of heat map. As a whole, Unigene0051978, Unigene0080699, Unigene0007996 were up-expressed and Unigene0074595 were down-expressed with the enhancement of Pi deficiency and the extension of time. Thus, the up-regulation of pyruvate dehydrogenase promoted the synthesis of acetyl-CoA by pyruvate, and the up-regulation of citric acid synthase promoted the synthesis of citric acid, and the down regulation of aconitase reduced the decomposition of citric acid, the up-regulation of isocitrate dehydrogenase enhanced the synthesis and metabolism of organic acids in the root system of Chinese fir by promoting the synthesis of α-ketoglutaric acid.
3.6 The contents of organic acids and the activity of related enzymes
The activities of NAD-malate dehydrogenase, citrate synthase, pyruvate dehydrogenase and phosphoenolpyruvate carboxylase in the root of Chinese fir seedlings after 0, 3 and 30 days of different phosphorus treatments were determined (Fig. 6-A). There was no significant difference in activity of NAD-malate dehydrogenase under different phosphorus treatments after 3 days. However, the activity of NAD-malate dehydrogenase under Pi deficiency conditionafter30 days was significantly lower than that under normal Pi supply. Compared with the normal Pi supply, the activity of citric acid synthase increased under the Pi deficiency condition after 3 and 30 days. The activity of pyruvate dehydrogenase increased significantly under Pi deficiency after 3 days, and increased with the increase of stress intensity, but decreased under Pi deficiency after 30 days with no significant difference.
The contents of malic acid, citric acid and pyruvic acid in Chinese fir seedlings were measured after 0, 3 and 30 days of different phosphorus treatments (Fig. 6-B). The contents of malic acid in the roots of Chinese fir were higher in Pi deficient condition than that under normal Pi supply after 30 days, and the higher the stress intensity, the higher the contents of malic acid. The content of citric acid in the roots of Chinese fir seedlings increased significantly after 3 days without Pi supply than that of the roots grown under normal Pi supply, and the stronger the breaking strength, the higher the content of citric acid. Under Pi deficiency, pyruvate content was lower than that of normal Pi supply after3 days, but increased after30 days.