Introduction

Camellia huana, the only native species of Camellia in Guizhou Province, China, is a perennial evergreen shrub belonging to the section Chrysantha Chang and the family Theaceae1. Its flowers and leaves are rich in tea polyphenols2,3, ascorbic acid4, polysaccharides5 and other physiologically active ingredients. Therefore, it is often used as a novel healthcare material and is widely bred in China. Notably, C. huana prefers shade and high humidity, and the soil environment should have high phosphorus (P) levels, intermediate potassium (K) levels and low nitrogen (N) levels6. Among them, as one of the essential elements of C. huana, P is crucial for its successful planting and cultivation in the future7. However, the P content in the soil of many areas in which it is introduced often cannot meet the planting requirements. Therefore, exploring the adaptability of C. huana in low-P environments is necessary. At present, studies of C. huana have focused mainly on individual selection and breeding8 or on specific functional genes9. Few systematic studies on its transcriptome and metabolites under low-P conditions have been conducted. Therefore, in this study, the Illumina HiSeq platform was used to perform transcriptome sequencing of C. huana under low-P stress, and the results were used to analyze the genes involved in the response to low P levels and elucidate the molecular mechanisms that might be involved in the adaptation to stress. Additionally, qRT‒PCR was performed to verify the sequencing results and changes in the related physiological and biochemical indices. The results of this study increase our understanding of the mechanism by which C. huana adapts to the soil environment and provide a theoretical reference for subsequent research on the introduction, protection, development and utilization of this plant resource.

Materials and methods

Plant material and experimental conditions

The experimental site was located in the medicinal botanical garden of Guiyang Health Vocational University (106°61′06″E, 26°61′97″N). The seedlings used in this study were obtained from Changtang town, Nanning, Guangxi Province, China (108.69°42′98′′E, 22.80°77′62′′N). The town is renowned for its proficient techniques and advanced management methods for growing C. huana, which could minimize individual variations caused by factors such as planting methods. Individuals of the same clones with similar heights, uniform growth and strong growth potential were selected as the experimental materials, and total nutritional Hoagland nutrient solution (Ca(NO3)2, 946 mg L− 1; KNO3, 506 mg L− 1; NH4NO3, 80 mg L− 1; KH₂PO₄, 136 mg L− 1; MgSO4, 493 mg L− 1; iron salt solution, 2.5 mL; and trace elements, 5 mL) was used for culture10. The plants were cultured in Hoagland nutrient solution for 30 days, after which they were watered every 4 days. The Hoagland nutrient mixture was then replaced with phosphorus-free Hoagland nutrient mixture (without KH₂PO₄), and the stress treatment was performed using the same experimental methods. Healthy leaves with intact shapes and no obvious lesions were sampled after 0 days (control check group, CK), 15 days (15 D) and 30 days (30 D) for subsequent experiments.

RNA extraction and transcriptome sequencing

Total RNA was extracted from different treatment (CK, 15 D, and 30 D) samples, with three biological replicates for each treatment, using a Plant RNA Isolation Kit (Majorbio, Shanghai, China). Sequencing libraries were constructed using the RNA Library Prep Kit for Illumina (NEB, Boston, MA, USA). Moreover, the trimmed mean M-value (TMM) method was used to normalize the data to eliminate the difference in sequencing depth between different samples. The libraries were subsequently sequenced using the HiSeq 4000 platform (Illumina, San Diego, CA, USA), and 150 bp (bp) paired-end reads were generated. Sequences with lengths less than 30 bp, N ratios greater than 10% and without inserted fragments for reasons such as connector self-connection and adapter sequences were removed using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) to obtain clean reads.

Transcriptome assembly, functional annotation and identification of differentially expressed genes (DEGs)

Since no reference genome for C. huana is available, the method of de novo assembly was used in this study. Trinity version 2.5.0 (https://github.com/trinityrnaseq/trinityrnaseq/wiki) was used to reassemble the clean reads, perform the classification and finally obtain the transcripts. The software was run with the default parameters. TransRate version 1.0.3 (http://hibberdlab.com/transrate/) and BUSCO version 3.0.2 (Benchmarking Universal Single-Copy Orthologs) (http://busco.ezlab.org) were used for the initial filtering, optimization and evaluation of the assembled sequence.

The sequences were subjected to a BLAST search against the Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.genome.jp/kegg/) public database to functionally annotate the genes. The KEGG enrichment analysis of the DEGs was performed with KOBAS version 3.0. In this software, false-positivity was assessed via the Benjamini–Hochberg procedure with a false discovery rate (FDR) correction cutoff E value of 10− 6. When the adjusted p value of a KEGG pathway was less than 0.05, we considered this KEGG pathway significantly enriched.

The gene expression level for each sample was determined according to the transcripts per million reads (TPM) using RSEM version 1.2.31 (Univ. of Wisconsin, Madison, WI, USA) with all parameters set to their defaults, in which the read counts were normalized using the edgeR package with the trimmed mean of M values method, and then the length of the gene (L) and the normalized counts (N) were used to calculate the TPM. The DEGseq R package was used to analyze the differential expression of genes in different samples. Based on the Benjamini and Hochberg (BH) correction for multiple tests, the significant differential expression thresholds were set as |log2(fold change)|>1 and P (P-adjusted) < 0.05. The heatmap and differential expression of genes among samples were analyzed using the pheatmap R package (Massachusetts General Hospital, Boston, MA, USA). The fold-change cutoffs for differential expression were consistent with the RNA-seq criteria.

qRT‒PCR validation of DEGs

Nine significantly differentially expressed genes were randomly selected from key pathways showing significant up- or downregulation (based on KEGG and DEG analyses) for qRT‒PCR verification of the accuracy of the transcriptome data, with Actin used as the reference gene11. Three independent biological replicates and three technical replicates for each biological replicate were run. Each reaction mixture consisted of 10 µL of 2×SGExcel FastSYBR qPCR mix (Sangon Biotech, Shanghai, China), 1 µL of forward and reverse primers (10 µM), 2 µL of complementary DNA (cDNA), and 6.0 µL of ddH2O. The gene-specific primers used are shown in Table 1. The reaction procedure consisted of the following steps: 95 °C for 3 min (predegeneration); 35 cycles of 95 °C for 5 s, 55 °C for 15 s, and 72 °C for 15 s (amplification); followed by 95 °C for 15 s, 60 °C for 1 min, 95 °C for 15 s, and 60 °C for 15 s (melting curve). The accuracy of the qRT‒PCR results was further validated by analyzing the melting curves. Gene expression levels were calculated using the 2−∆∆Ct method, with all the data from the qRT‒PCR experiments compared with those obtained from transcriptome sequencing.

Table 1 Primers used in this study.
Full size table

Validation of physiological and biochemical parameters

The levels of ATP, superoxide dismutase (SOD; EC 1.15.1.1), catalase (CAT; EC 1.11. 1.6), soluble sugar (SS), starch and plant chlorophyll were measured using the methods described in the Experimental Principles and Techniques of Plant Physiology and Biochemistry (3rd Edition)12. The levels of the four hormones abscisic acid (ABA), ethylene (ETH), methyl jasmonate (MeJA) and indoleacetic acid (IAA) were determined according to the cellulose absorption index (CAI) method13. Glutathione S-transferase (GST) activity was determined using the methods reported by Materon14. All of the above measurements were biologically repeated three times for each treatment.

Results

Transcriptome sequencing and de novo assembly

The raw transcriptome data can be found in the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/) Sequence Read Archive (SRA) database (accession number: PRJNA1101880). A total of 138,000,514 raw reads, 68,414,648 clean reads and 20,489,133,698 clean bases were obtained by filtering the raw data, with an average error rate of 0.01216% (Table 2). Moreover, 17,703,312, 18,203,342 and 17,456,012 reads were matched in the three groups, respectively, with matching rates exceeding 77%, indicating a high degree of matching. The contents of Q20 and Q30 were more than 98% and 95%, respectively. According to industry standards in plant transcriptomics, high-quality sequencing results require Q20 to be greater than 90% and Q30 to be greater than 85%. Our results are significantly higher than these indicators. The duplicate sequences were all relatively low in different groups, indicating that the sequencing quality reached a high standard. The results of the quality control analysis revealed that the sequencing quality was high and that the error rate was low, which met the requirements for subsequent analysis.

The numbers of assembled unigenes and transcripts were 149,012 and 243,443, with average lengths of 904.19 bp and 1031.99 bp, respectively, and N50 lengths of 1359 bp and 1634 bp, respectively. Clean reads were assembled using trinity software, and 45% of the total 200–500 bp sequences, 29% of the total 501–1000 bp sequences and 10% of the total 1001–1500 bp sequences were unigene sequences (Table 3). The assembly results were evaluated using transrate and BUSCO, and the percentages of unigenes were 0.28432 and 68.6%, respectively, whereas the percentages of transcripts were 0.40282 and 89.3%, respectively. Furthermore, the percentages of fragmented and missing BUSCO genes were low. Combining the N50 length and E90N50 length data, we believe that the sequencing results are relatively reliable and can be used for downstream analyses.

Table 2 Summary of sequencing data quality control.
Full size table
Table 3 Length distribution and software evaluation of unigenes and transcripts.
Full size table

Gene expression and KEGG enrichment analyses

A total of 5,030 genes were differentially expressed between the samples from the control and phosphorous-deficient plants. Among them, 1228, 1877 and 2670 genes were differentially expressed between the 15 D and CK groups, between the 30 D and CK groups, and between the 30 D and 15 D groups, respectively (Fig. 1A). Among all the DEGs, 426/802 and 861/1678 genes were up-/downregulated at 15 D and 30 D, respectively, compared with those in the CK group, and 1240/1430 were up-/downregulated at 30 D compared with those at 15 D (Fig. 1B). The KEGG enrichment analysis revealed that these genes were involved mainly in different pathways, including motor protein-related pathways; photosynthesis; starch and sucrose metabolism; plant hormone transduction; the MAPK signaling pathway; and flavonoid biosynthesis and membrane lipid metabolism (Fig. 1C, D and E).

Fig. 1
figure 1

Differential gene expression in C. huana under low-P stress. (A) Venn diagram of DEGs. (B) Statistical data of DEGs between different comparisons. Red and blue represent up- and downregulated expression, respectively. (C–E) KEGG enrichment analysis of DEGs between the 15 D and CK groups, the 30 D and CK groups, and the 30 D and 15 D groups. The adjusted P values are shown in Tables A1–A3 in Appendix A.

Full size image

qRT‒PCR validation

The abbreviations of all the genes are shown in Table A in Appendix A. The experimental results showed that the expression levels of AGP, APR, AUX2, LHCB3, LPIN, and PIP continued to decrease with increasing treatment time and reached their lowest values on day 30 (Fig. 2). The qRT‒PCR results for the six genes were very close to those of RNA-seq according to the TPM algorithm. However, the qRT‒PCR data for ACS and SnRK2 revealed that their expression first increased, peaked at 15 D, and then gradually decreased but remained higher than their expression in the CK group, whereas the expression level of TPM continued to increase. In contrast to the expression of ACS and SnRK2, the qRT‒PCR results showed that the expression level of PK decreased first but then increased, whereas the expression level of TPM continued to decrease. In general, the qRT‒PCR results were essentially consistent with the results for the key DEGs identified by transcriptome sequencing, indicating that the transcriptome sequencing results were reliable.

Fig. 2
figure 2

Changes in the expression of the nine randomly selected genes were determined via qRT‒PCR and sequencing. The x-axis represents different processing times, and the y-axis represents changes in gene expression under stress. The data show the fold change in the expression of each gene relative to the control conditions. Red indicates the RNA-sequencing results obtained using the TPM algorithm, and blue indicates the qRT‒PCR results. The related expression value of the RNA-seq and qRT‒PCR are shown in Table A4 in Appendix A.

Full size image

Validation of physiological and biochemical parameters

Figure 3A shows the energy metabolism pathway of C. huana in a low-P environment, which involves mainly glucose metabolism and the enzymes and other substances involved in this process. Figure 3B shows the relative changes in the expression of key genes in different metabolic pathways. In the glycometabolic pathway, the expression levels of β-amylase, SP and HK were significantly upregulated (p < 0.05, the same as below). However, the expression of AGP decreased gradually, and the expression of SSS, SBE and GBSS, the key genes involved in starch synthesis, did not change significantly (Fig. 3B). In the Embden–Meyerhof–Parnas (EMP) pathway, the expression of pyruvate kinase (PK), which catalyzes the transformation of phosphoenolpyruvate (PEP) to pyruvate, was significantly downregulated, which was most obvious in the 15 D group compared with the CK group. Notably, the expression levels of all the genes involved in photosynthesis tended to decrease to different degrees with increasing treatment time. In conclusion, stress induced by low-P conditions might have negative effects on the energy metabolism of C. huana.

Physiological and biochemical verification experiments revealed that with increasing treatment time, the total chlorophyll content first tended to increase but then tended to decrease, but a significant difference was not observed compared with the CK group (Fig. 4). The starch content continued to increase and reached its peak at 30 D (40.5823 mg/g), an increase of 133.5% compared with 17.3870 mg/g in the CK group, which was similar to the significantly upregulated trend of SP expression at 30 D (Fig. 3A). The level of soluble sugar continued to decrease from 26.3780 mg/g to 16.0560 mg/g at 30 D, which was closely related to the continuous downregulation of PK, a key enzyme in the EMP pathway (Fig. 3B). The ATP content also gradually decreased from 0.1637 µmol/g to 0.1395 µmol/g at 30 days, which was consistent with the systematic decrease in the activity of the photosynthesis-related genes psaA and psbA (Fig. 3B).

Fig. 3
figure 3

Effects of low-P stress on energy metabolism. (A) The main energy pathways in plants. The black and yellow arrows represent glycometabolism and the EMP pathway, respectively. The different enzymes are shown in red. (B) Fold changes in the expression of genes involved in metabolic pathways in response to stress. White indicates no change, red indicates upregulation, and blue indicates downregulation in each treatment group, as shown in the color bar for a log2 fold change scale. The abbreviations in the figure are shown in Tables A5 in Appendix A.

Full size image
Fig. 4
figure 4

Effects of low-P stress on different physiological and biochemical indices. Significant differences (p < 0.05) are indicated with lowercase letters. The specific data of the relevant indicators are shown in Table A6 in Appendix A.

Full size image

Cytomembrane lipid metabolism under stress

Figure 5 shows the general cytomembrane lipid metabolism pathway and changes in the expression of related genes in C. huana under low-P stress. The expression of genes involved in phospholipid metabolism was generally downregulated in response to stress. Among them, glycerol-3-phosphate acyltransferase (GPAT), which catalyzes the formation of phosphatidic acid from glycerol-3-phosphate; LPIN, which is required to form diacylglycerols; and PIP (phosphatidylinositol 4 phosphate 5-kinase), which is required to form phosphatidylinositols, exhibited the most obvious changes, and their expression levels were significantly different from those in the CK group after 30 days of treatment. Among the genes involved in galactolipid metabolism, the expression of PLA (phospholipase A) decreased significantly, while that of PLC increased significantly, and the expression of the remaining genes did not change significantly (Fig. 5B). In contrast to the genes involved in galactolipid metabolism, the genes involved in sulfur lipid metabolism were generally upregulated. Among them, the expression levels of UGDH, APS and APR were significantly higher in the 15 D and 30 D groups than in the CK group (Fig. 5B).

Fig. 5
figure 5

Effects of low-P stress on cytomembrane lipid metabolism. The black, orange and purple arrows represent the phospholipid, galactolipid and sulfur lipid pathways, respectively. (A) The different genes are presented in red font. (B) Changes in the expression of genes involved in metabolic pathways in response to stress. The abbreviations in the figure are shown in Tables A5 in Appendix A.

Full size image

Hormone metabolism and physiological and biochemical verification

Figure 6 shows the relationships between hormone metabolism and the expression of related genes in C. huana under low-P conditions. In general, among the four tested hormones, the IAA and MeJA synthesis pathways were generally inhibited, whereas the ETH and ABA synthesis pathways were activated. Among them, the expression levels of all the genes involved in IAA synthesis (green arrow in Fig. 6A), except cytochrome P450 monooxygenase (CYP), were downregulated to varying degrees and reached significant levels (Fig. 6B). The expression levels of all genes involved in the MeJA metabolic pathway were reduced. Among them, allene oxide cyclase (AOC) expression decreased most significantly (Fig. 6B), and the log2 fold changes in expression indicated 2.20- and 3.32-fold decreases in the 15 D and 30 D groups, respectively, compared with those in the CK group. On the other hand, the expression patterns of ETH- and ABA-related genes tended to increase under low-P stress. The expression of 9-cis-epoxycarotenoid dioxygenase (NCED) and dehydrogenase/reductase (DR), the key genes involved in ABA synthesis, increased continuously over time (Fig. 6B). Compared with the expression of genes in the ABA synthesis pathway, the expression of all genes in the ETH synthesis pathway also tended to increase, but the degree of increase was not as high as that of the ABA synthesis-related genes (Fig. 6B).

The physiological and biochemical results showed that the contents of IAA and MeJA were significantly decreased with increasing treatment time, reaching their lowest values at 30 D, with contents of 47.0286 ng g− 1 and 74.8373 ng g− 1, respectively. The differences between the treatment groups reached a significant level (Fig. 7). In contrast, the contents of ABA and ETH were significantly increased, peaking at 30 days (280.0712 ng g− 1 and 5.4879 µg g− 1, respectively). In addition, low-P stress significantly affected the activity of antioxidant enzymes. The activities of superoxide dismutase (SOD), glutathione S-transferase (GST), and catalase (CAT) increased significantly over time, peaking at 30 D, with contents of 671.20 U g− 1, 193.64 U g− 1, and 1813.84 U g− 1, respectively. The contents of SOD and GST were approximately 4.5 times higher than those in the CK group at 30 days, and the increase in CAT activity was the most significant, which was nearly 7 times higher than that in the CK group at 30 days.

Fig. 6
figure 6

Effects of low-P stress on hormone metabolism. (A) Hormone biosynthesis pathways. The orange, blue, green, and purple arrows represent the ABA, MeJA, IAA, and ETH pathways, respectively. The different genes are presented in red. (B) Changes in the expression of genes involved in metabolic pathways in response to stress. The abbreviations in the figure are shown in Tables A5 in Appendix A.

Full size image
Fig. 7
figure 7

Effects of low-P stress on different physiological and biochemical indices and hormone and antioxidant enzyme activities. Significant differences are indicated with lowercase letters. The specific data for the relevant indicators are shown in Table A7 in Appendix A.

Full size image

Discussion

Effects of low-P stress on the energy metabolism of C. huana

P plays a central role in most major metabolic processes in plants, especially photosynthesis and respiration. In terms of photosynthesis, the results of this study indicated that the expression levels of genes decreased to varying degrees with increasing stress, indicating that an insufficient P supply strongly affects the normal photosynthesis of C. huana. Studies have shown that P deficiency can inhibit the expression of psaA, psaO, and psbA, which are involved in ATP synthesis for photosynthesis15,16,17, and reduce the synthesis rate of ribulose-1,5-bisphosphate (RuBP), thereby reducing the ability of leaves to assimilate CO218,19. These two factors could lead to a decrease in the photosynthetic rate of the leaves, which strongly affects the photosynthetic characteristics of plants, including decreases in chlorophyll20 and soluble sugar21 contents. Our qRT‒PCR (Fig. 2) and physiological and biochemical test results (Fig. 4) confirmed this conclusion. On the other hand, as a key rate-limiting enzyme in the EMP pathway, PK was significantly downregulated in plants under stress (Fig. 3), suggesting a metabolic shift from glycolysis to alternative ATP generation pathways, possibly oxidative phosphorylation or lipid oxidation, to compensate for energy deficits under P deficiency. Similar metabolic adjustments have been observed in Cotton Genotypes22, Beta vulgaris23 and Solanum lycopersicum24. This phenomenon is essential for maintaining energy homeostasis and plant growth under stress.

Moreover, the starch content in C. huana gradually increased (Fig. 4). This finding was consistent with the results reported by Zhang et al.25, who studied changes in the starch content in maize (Zea mays L.) leaves under low-P stress. Previous studies have shown that plants maintain the exchange of metabolites between the cytoplasm and chloroplasts through the accumulation of starch in P-deficient environments26,27. A similar phenomenon was also observed in C. huana, where starch accumulation was suggested as an adaptive strategy to cope with P deficiency by serving as a temporary carbon reservoir and energy buffer28. In the starch synthesis pathway, three different starch-forming synthetases (GBSS, SSS and SBE) are key regulatory enzymes. In this study, the expression level of SBE was higher than that of GBSS and SSS under stress (Fig. 3). SBE generates amylopectin by specifically binding to starch granules. Due to the high expression of SBE, we suspected that amylopectin would tend to be synthesized in C. huana under low-P stress. However, compared with different starch-forming synthetases, the expression of SP was significantly higher (Fig. 3), which indicated that in a low-P environment, the generated starch would be quickly degraded and converted into sugars that could be directly used, such as glucose, to quickly supplement energy29,30,31. Yao et al.32 and Jiang et al.33 also found that low-P conditions led to the upregulation of SP expression when studying the starch synthesis mechanism of Tetraselmis subcorformis, similar to the results of this study.

Effects of low-P stress on lipid metabolism in the cell membrane of C. huana

In this study, most of the genes involved in phospholipid metabolism were downregulated, possibly due to an insufficient P content and the restriction of phosphate, a key raw material required for phospholipid synthesis34. Notably, for genes involved in the galactose lipid metabolism pathway, the expression of phospholipase C (PLC) in the phospholipase gene family was significantly increased, whereas that of phospholipase A (PLA) was significantly decreased (Fig. 5B). This pattern suggests a compensatory mechanism in which PLC activity increases to facilitate membrane remodeling and maintain cell integrity, despite a reduction in overall phospholipid synthesis35. The experimental results described by Zheng et al.36 for P deficiency stress in Triticum aestivum also showed that although the total phospholipid content decreased, the accumulation of diacylglycerol (DGDG) and quinoline diacylglycerol (SQDG) increased. These changes are beneficial to the function and membrane integrity of chloroplasts, jointly improving tolerance. Previous studies have shown that PLA is involved mainly in the hydrolysis of glycerol phospholipids37,38, and in a low-P environment, the level of glycerol phospholipids will be affected, resulting in the inhibition of PLA. On the other hand, PLC mainly catalyzes the synthesis of ester bonds between phospholipids and participates in the repair of cell membranes39. Due to the lack of P, the synthesis and maintenance of the cellular phospholipid bilayer were affected; therefore, the expression of PLC was significantly increased, which could protect the function of the membrane phospholipid bilayer in this environment. Little change in the expression of other related genes was observed in response to stress, indicating that the galactose lipid pathway may not be the main regulatory pathway involved in the response to low-P conditions. Because metabolic factors are multifaceted, the reason for the change in this mechanism remains to be verified. Most of the genes involved in sulfur lipid metabolism were upregulated. Sulfur has many functions in plants, including participating in signal transduction and stress resistance40. In the absence of P, sulfur-containing lipids such as sulfoquinoline diacylglycerol (SQDG) can partially replace phospholipids in chloroplast membranes, thus maintaining photosynthetic efficiency and membrane stability41,42. In addition, sulfur metabolism plays a key role in regulating antioxidant defense mechanisms, which could help plants mitigate oxidative stress caused by a phosphorus deficiency43. Therefore, the increase in sulfur lipid metabolism can not only compensate for the decrease in phospholipids but also increase the adaptability of plants in low-phosphorus environments by regulating cellular redox homeostasis and membrane lipid remodeling. The KEGG results for C. huana in this study also revealed the upregulation of its lipid metabolism pathway, indicating that it may adapt to a low-P environment through membrane lipid remodeling and P recovery mechanisms44,45.

Effects of low-P stress on hormone metabolism in C. huana

P is an essential element for plant growth, and its supply directly affects the synthesis, transport and response of plant hormones. The results showed that genes involved in the ETH and ABA synthesis pathways were generally upregulated (Fig. 6), consistent with the physiological and biochemical results (Fig. 7). Liu et al.46 and Wang et al.47 reported that P deficiency promoted the synthesis of ETH in plants, which inhibited the expression of the root hair suppressor factor GLABRA2 (GL2) by activating the core transcription factor ethylene-insensitive 3 (EIN3), thereby promoting ectopic root formation, including increased branching and root hair growth48,49,50. ETH improves P utilization by altering the root morphology and optimizing the rhizosphere microenvironment, maintaining growth and development in P-deficient environments. The results of the ABA analysis revealed that with increasing P deficiency stress duration, the expression of all the genes in the ABA synthesis pathway increased, especially that of the NCED gene. The catalysis of nine-cis-neoxanthin conversion to xanthoxin by NCED is a key step in the ABA metabolic pathway in plastids. Therefore, an increase in NCED expression promoted the expression of other genes in the ABA metabolic pathway, thereby increasing ABA accumulation (Fig. 7). One of the main functions of ABA in plants is to regulate stomatal opening. Under P-deficient conditions, the photosynthetic rate of leaves decreases, which leads to the shrinkage of mesophyll cells, resulting in a decrease in stomatal conductance on the leaf surface51. ABA regulates stomatal opening and closing to maintain water and nutrient balance in plants, thereby increasing plant stress resistance.

Regarding other hormones, the biosynthesis of IAA involves four pathways (Fig. 6A). In a low-P environment, genes involved in the tryptamine, IPA and IA pathways were inhibited, whereas CYP expression was significantly upregulated. These findings suggest that under stress conditions, C. huana primarily synthesizes IAA via the CYP pathway, likely as an adaptive response to P limitation. The downregulation of IAA inhibited excessive leaf growth and promoted root development, enabling plants to access resources better and reduce unnecessary energy consumption in P-limited environments52, similar to the upregulation of ABA. The expression levels of MeJA synthesis-related genes were generally reduced (Fig. 6B), and the content of MeJA also decreased with increasing stress (Fig. 7). Researchers have speculated that a reduction in MeJA levels is conducive to activating the defense system and increasing plant adaptability to low P-induced stress53. The changes in different hormones generally promote better resource acquisition and reduce unnecessary energy consumption for plants in P-limited environments54. The physiological and biochemical results indicated that the levels of SOD, peroxidase (POD; EC 1.11.1.7), CAT, glutathione S-transferase, ABA and ETH all increased gradually, consistent with the results reported by Wang et al.55, and confirmed the accuracy of the transcriptome data.

In this study of low-P stress in C. huana, the physiological verification and transcriptome data were generally consistent, but differences in the expression of some genes, including PK, ACS and SnRK2, were observed. These differences may be due to experimental errors in the implementation of the two measurement methods, such as differences caused by qRT‒PCR primer specificity, amplification efficiency, or a lag effect between mRNA accumulation and enzyme activity56,57. Overall, changes in enzyme activity reveal the complexity of plant adaptation mechanisms under low-P stress, suggesting that the regulation of enzyme activity is influenced by multiple factors.

Conclusions

The effect of decreased P levels on plants is known; however, little research has been conducted specifically on C. huana. In this study, we explored the molecular response of C. huana to low-P stress. The results showed that the genes whose expression changed significantly were involved mainly in sugar and energy metabolism, cell membrane lipid metabolism and hormone metabolism. The starch content increased through the actions of multiple genes. Moreover, the photosynthetic rate of the whole plant decreased significantly during the different periods. Significant changes in the expression of genes involved in lipid metabolism were observed at 30 D. Moreover, the levels of ABA and ETH gradually increased, whereas those of IAA and MeJA decreased. This report is the first to describe the high-throughput sequencing of P deficiency-tolerant C. huana, and this study provides novel insights into genes related to phosphorous deficiency tolerance and a valuable molecular basis for further studies of C. huana. At the breeding level, C. huana varieties or wild populations with high utilization of low P levels should be selected, and their low-P-tolerant characteristics should be evaluated using molecular markers or phenotypes. Low-P-tolerant parents were selected for hybridization, and generation breeding was conducted in a low-P environment to improve low-P adaptability. The survival ability of low-P-tolerant varieties in poor soil can be improved by breeding them. The ecological planting of C. huana should be realized, and the quality and market value of tea should be improved.