Next Article in Journal
Atractylodes macrocephala Root Rot Affects Microbial Communities in Various Root-Associated Niches
Previous Article in Journal
Detection of Rice Leaf Folder in Paddy Fields Based on Unmanned Aerial Vehicle-Based Hyperspectral Images
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Agronomy
Volume 14
Issue 11
10.3390/agronomy14112661
agronomy-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Integrative Transcriptomic and Metabolomic Analysis Reveals Quinoa Leaf Response Mechanisms to Different Phosphorus Concentrations During Filling Stage

by
Hongxin Wang
1,†,
Hanxue Li
1,†,
Xiaorong Li
2,†,
Qianchao Wang
1,
Junna Liu
1,
Ping Zhang
1,
Heng Xie
1,
Li Li
1 and
Peng Qin
1,*
1
College of Agronomy and Biotechnology, Yunnan Agricultural University, Kunming 650201, China
2
Chuxiong Academy of Agricultural Sciences, Chuxiong 675000, China
*
Author to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Agronomy 2024, 14(11), 2661; https://doi.org/10.3390/agronomy14112661
Submission received: 25 August 2024 / Revised: 23 September 2024 / Accepted: 8 November 2024 / Published: 12 November 2024
Browse Figures
Versions Notes

Abstract

:
Quinoa is an annual self-pollinating plant rich in polyphenols, flavonoids, saponins, and amino acids; its protein balance closely aligns with the ideal recommendation set by the Food and Agriculture Organization. Therefore, quinoa is considered the most suitable “all-nutrient food”. Phosphorus fertilization plays an important role in restricting the growth and development of quinoa; however, the effects of phosphorus fertilizer on quinoa growth remain unclear. Therefore, we conducted metabolome and transcriptome analyses on quinoa leaves during the filling stage, subjecting plants to different doses of phosphorus fertilizer. Overall, phosphorus treatment exerted a significant impact on the phenotypic characteristics of quinoa. Specifically, through a combined analysis of ultra-performance liquid chromatography–tandem mass spectrometry and transcriptome analysis, we identified the alteration and regulation of specific metabolites and genes within flavonoid biosynthesis pathways; this comprehensive evaluation helped elucidate the response mechanism of quinoa leaves during the grouting stage under various phosphorus conditions. Ultimately, the results of this study provide a reference for the selection of quinoa cultivars that exhibit tolerance to low- or high-phosphorus stress; additionally, we offer a theoretical basis for the rational application of phosphorus fertilizer and the enhancement of phosphorus utilization efficiency.

1. Introduction

Quinoa (Chenopodium quinoa Willd.) is an annual, dicotyledonous, self-pollinating, herbaceous plant in the subfamily Chenopodium of the Amaranthaceae family, with alternating single leaves. It is native to the Andes Mountains in South America and has a planting history of over 5000–7000 years. Quinoa, often referred to as a “super grain” [1,2,3,4,5], possesses a strong ability to adapt to the surrounding environment, with significant resistance to cold, drought, saline-alkali, and barren conditions; nonetheless, it is considered a C3 crop that prefers cold and high altitude. Furthermore, quinoa has been confirmed to possess resistance against plant pathogens [4,6]. Quinoa is considered a non-genetically modified nutritional food and stands out as the only single crop that provides all the essential nutrients required by humans [7,8,9]. In particular, quinoa is a whole-protein food, with an exceptionally high protein content, making it readily absorbed by the human body. Further, quinoa seeds are abundant in nine essential amino acids required for human physiological functions; notably, they also contain lysine, an amino acid typically absent in ordinary grains. Quinoa also has antibacterial, anti-inflammatory, blood pressure-lowering, and antioxidant effects [10,11,12,13]. Therefore, it is currently considered a popular health food among consumers.
Phosphorus is an abundant nutrient that affects crop growth, development, and yield. In particular, it promotes crop root development, enhances cold and drought resistance, and promotes plant maturity [14]. Further, it has been demonstrated to increase the yield of spike and plump grains, participate in various plant metabolic activities, and affect plant physiology and morphology. Phosphorus also plays an important role in photosynthesis, carbohydrate synthesis, transportation, nitrogen metabolism, and fat synthesis in plants [15]. As the seedling period is extended into the subsequent rapid growth phase, the inhibitory effects of low-phosphorus stress on plants become more pronounced [16]. The supply level of phosphorus directly governs the accumulation of phosphorus in cucumber and tomato seedlings, which in turn indirectly affects the content of other mineral elements and chlorophyll; ultimately, this alters the growth and development process of seedlings and influences the formation of their morphological characteristics [17]. Additionally, phosphorus deficiency or excessive phosphorus can inhibit the growth of soybean seedlings [18].
Functional genomics can be used to study the overall genetic expression and metabolic responses of plants when growing under limited- or excessive-phosphorus conditions; for example, a prior study demonstrated that stress-related metabolites, such as polyols, accumulate in phosphorus-deficient roots alongside sugars, which play a critical role in the induction of phosphorus stress genes [19]. In response to phosphorus deficiency, plants must first detect the occurrence of nutritional stress and transmit the appropriate signals to promote adaptation. The change in organic acid content observed in phosphorus-stressed roots is attributed to a notable reduction in the expression of TC 1864 isocitric acid dehydrogenase, a key regulatory enzyme in the tricarboxylic acid cycle [20]. Further, phosphorus stress induces oxidative reactions in soybeans, characterized by increased lipid peroxidation, elevated peroxide levels, and increased catalase and peroxidase activities [19]. Under phosphorus-deficient conditions, plants have also been observed to accumulate sugar and starch in their leaves; this accumulation leads to an increase in the threonine load in the phloem, facilitating the relocation of carbon resources to the roots and increasing their size relative to the tender shoots [21].
At present, various studies have explored the effects of different fertilizer compositions on the growth and development of quinoa; however, there remain limited reports on the effects of phosphorus fertilizer content on the growth and development of quinoa during the filling stage. During the filling stage, quinoa exhibits the highest rate of dry matter accumulation and leaf development. In the present study, ultra-high-performance liquid chromatography–tandem mass spectrometry (UPLC-MS/MS) and transcriptome analyses were used to comprehensively evaluate the effects of different phosphorus conditions on the differential metabolite and gene expression of quinoa at the grouting stage. Additionally, by assessing the morphological indicators, the response mechanisms of quinoa seedlings to different phosphorus levels were elucidated.

2. Materials and Methods

2.1. Material and Sample Preparation

Dianli-2912, a high-grade cultivar of white quinoa independently selected by Yunnan Agricultural University, was planted in the modern agricultural education and research base of Yunnan Agricultural University in Xundian County, Kunming (E 102°41′, N 25°20′). Uniform and consistent seeds were selected and evenly seeded in a pot (117 × 39 × 65 cm); the substrate used contained 2.75 g/kg, 1.66 g/kg, and 1.18 g/kg of CH4N2O, P2O5, and K2O, respectively. We then applied P2O5 at doses of 0 kg/hm2 (W0), 112.5 kg/hm2 (W1), 225.0 kg/hm2 (W2), 337.5 kg/hm2 (W3), 450.0 kg/hm2 (W4), and 562.5 kg/hm2 (W5) to these quinoa seeds. During this period, the fertilization rates of CH4N2O and K2O were maintained at 112.5 kg/hm2, ensuring approximately 50 seedlings in each pot. During the early stage, conventional cultivation and management techniques were followed, including an average temperature of 25.6 °C, approximately 10 h of sunlight, and a sowing depth of 2–3 cm, and using a 1:1 soil ratio of humus to red soil. Initial fertilization was conducted during the two-leaf stage of quinoa development; then, the second fertilization was conducted 30 d after the initial fertilization. In order to avoid errors, plants were marked on the day of flowering in this experiment, and we started sampling at 10 a.m. 15 days after the flowering with an average temperature of 25.6 °C and a rainfall at 0.0 mm on the sampling day. Immediately after sampling, we placed them in liquid nitrogen for rapid freezing and transfered them to −80 °C for storage. Quinoa exhibited the greatest phenotypic differences at the grouting stage, with stiff seedlings at a P2O5 dosage of 0 kg/hm2, vigorous growth at a P2O5 dosage of 450.0 kg/hm2, and a decreasing growth trend at a P2O5 dosage of 562.5 kg/hm2. Metabolomic and transcriptome analyses were performed on the leaves, with 6 gradients of P2O5 applied and 3 biological replicates, for a total of 18 samples (Table 1).

2.2. Morphological Data Collection

Samples were taken from the filling stage leaves of the quinoa plants 15 days after the flowering (repeated three times) to determine the plant height, leaf area, stem thickness, and other morphological indicators. Next, the quinoa plant height and stem thickness (distance from the base to the top of the unfolded leaf tip) were measured using vernier calipers. Plants were then killed at 110 °C for 30 min and dried to a constant weight at 85 °C to determine their dry weight. The leaf area of the seedlings was measured using a TPYX-A crop leaf morphology measuring instrument. Data are presented as the mean ± standard deviation. For all analyses, SPSS software (version 27.0) was used to determine significant differences via one-way analysis of variance (ANOVA) testing. A p-value < 0.05 was considered statistically significant.

2.3. Qualitative and Quantitative Analysis of Metabolite Extracts

2.3.1. Sample Preparation and Metabolite Extraction

We accurately weighed 200 mg (±1%) of the sample and subjected it to vacuum freeze drying (Scientz-100F). We accurately added 0.6 mL of methanol (including internal standard) (Table S1), added 100 mg of glass beads, and vortexed for 1 min. Grinding was performed at 50 Hz for 60 s in a high-throughput tissue grinder and repeated twice. Room temperature ultrasound was conducted for 15 min. We centrifuged at 12,000 rpm and 4 °C for 10 min, took a 200 µL supernatant, and added it to the detector bottle. The detector bottle was a 2 mL brown injection bottle, which was used to contain the extract to be measured and injected into the liquid phase injection tray. Quality control (QC) samples were used to correct for deviations in the analytical results of the mixed samples and for errors due to the analytical instrument itself. In total, 20 µL of each sample to be tested was mixed into a QC sample. Chromatographic separation was used with an ACQUITY UPLC® HSS T3 (150 × 2.1 mm, 1.8 µm, Waters, Milford, MA, USA) column maintained at 40 °C. The temperature of the autosampler was 8 °C. The gradient elution of analytes was carried out with 5 mM ammonium formate (Brand: sigma, Darmstadt, Germany) in water (A) and acetonitrile (B) at a flow rate of 0.25 mL/min or 0.1% formic acid (Brand: TCI) in water (C) and 0.1% formic acid in acetonitrile (Brand: Thermo Fisher, Waltham, MA, USA) (D). The injection of 2 μL of each sample was performed after equilibration. An increasing linear gradient of solvent B (v/v) was used as follows: 0–1 min, 2% B/D; 1–9 min, 2–50% B/D; 9–12 min, 50–98% B/D; 12–13.5 min, 98% B/D; 13.5–14 min, 98–2% B/D; and 14–20 min, 2% D-positive model (14–17 min, 2% B-negative model). The mass spectrometer instrument used an electrospray ionization source (ESI) in positive and negative ion ionization modes with 3.50 kV for positive ion sprays and 2.50 kV for negative ion sprays (the correction fluid model: thermoA39239). Sheath gas and auxiliary gas were set at 30 and 10 arbitrary units, respectively. The capillary temperature was 325 °C. The analyzer scanned over a mass range of m/z 81–1000 for a full scan at a mass resolution of 60,000. Data-dependent acquisition (DDA) MS/MS experiments were performed with a High-Energy Collisional Dissociation (HCD) scan. HCD breaks up molecules by means of high-energy collisions, resulting in the production of fragment ions, which are subsequently detected and analyzed by a mass spectrometer. The normalized collision energy was 30 eV. Dynamic exclusion was implemented to remove some unnecessary information in MS/MS spectra.

2.3.2. Qualitative and Quantitative Analysis of Metabolites

The preparation, extract analysis, and qualitative and quantitative analyses of metabolites from the 18 quinoa seedling samples were performed according to the standard procedures described by BioNovoGene (Suzhou, China). UPLC-MS/MS can be used to accurately perform qualitative and quantitative analyses of the metabolites in samples. To obtain reliable and high-quality metabolomics data, quality control (QC) is required for this procedure [22]. In the present study, QC sample aggregation with good repeatability indicated that the system was stable [23]. A principal component analysis (PCA) score graph was used to observe the degree of aggregation and dispersion of the samples. The metabolite content data were normalized by unit variance scaling (UV). All samples and related data were calculated using a distance matrix; subsequently, hierarchical clustering of these samples was conducted to form a tree graph that illustrated the similarity between them. A hierarchical clustering analysis was conducted to evaluate the relative levels of different metabolites under varied experimental conditions; these results were visualized using heat maps. Ultimately, by linking metabolites with the same or similar metabolic patterns into classes, the biological functions of known and unknown metabolites could be inferred. Using the OPLS-DA model, we analyzed these metabolomic data; then, by producing score maps and permutations for each group, we further demonstrated the differences between each group [24]. To identify differential metabolites, we conducted a screening process with specific criteria. The relevant screening conditions for the differentially expressed metabolites were as follows: p-value ≤ 0.05 and variable importance in projection (VIP) ≥ 1 [25]; p-value ≤ 0.05 and fold change ≥ 1.5 or ≤0.667 [26]; and p-value ≤ 0.05 for multiple groups [27]. The identification of metabolites was first verified based on their precise molecular weights (with a molecular weight error < 15 ppm). Then, we used Metlin (http://metlin.scripps.edu) and MoNA (https://mona.fiehnlab.ucdavis.edu/, accessed on 3 December 2021) to analyze the specific MS/MS fragmentation pattern. Additionally, we cross-referenced and annotated the metabolite results with the BioNovoGene self-built standard product database. Finally, to assess the consistency of the changing trends among the identified metabolites, we conducted a correlation analysis of the differential metabolites. The correlation between each metabolite was analyzed by calculating the Pearson correlation coefficients for all pairs of metabolites.

2.4. Transcriptome Sequencing and Data Analysis

2.4.1. RNA Extraction, Quantification, Sequencing, and Data Analysis

After RNA extraction, purification, and library construction, the 18 quinoa leaf samples were subjected to next-generation sequencing using an Illumina sequencing platform for paired sequencing of these libraries. The upgraded version of HISAT2 software (http://ccb.jhu.edu/software/hisat2/index.shtml, accessed on 3 December 2021) from TopHat2 was used to compare the filtered reads to the reference genome. For this mapping, we established that if the reference genome selection is appropriate and there is no contamination in the relevant experiments, the mapping ratio of the sequences should typically be >70%. HTSeq was used to statistically compare the read count value of each gene with the original gene expression level. To ensure comparability between the gene expression levels of different genes and samples, we standardized expression using fragments per kilobase per million fragments (FPKM). In the reference transcriptome, a gene is considered to be expressed, with an accurate expression level, with an FPKM > 1. Additionally, the correlation of gene expression levels between samples served as an important indicator of experiment reliability and the appropriateness of sample selection. The correlation of gene expression levels between samples was assessed using Pearson’s correlation coefficient.
A Principal Component Analysis (PCA) reduces high-dimensional data to two or three dimensions by linear transformation while maintaining the features that contribute the most to the variance of each party, i.e., reducing the complexity of the data. When there are multiple samples, we used the DESeq package in R to perform PCA (Principal Component Analysis) on each sample based on expression. PCA was used to link similar samples together; specifically, samples with a shorter distance between them were considered to possess a higher similarity. We used DESeq [28] to conduct a differential analysis of gene expression with the following screening conditions: an expression difference multiple log2 fold change ≥ 1 and a p-value < 0.05. Based on the results of the differential gene expression analysis, a Venn diagram and upset diagram were produced; ultimately, these could be used to indicate the number of common unique differential genes among the comparison groups, perform a hierarchical clustering analysis, and produce a corresponding clustering heat map for each group. According to the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment results, the degree of enrichment was determined by the rich factor, FDR value, and the number of genes enriched in the pathway.

2.4.2. RT-qPCR

To verify the reliability of the transcriptome sequencing results, differential genes were selected for RT-qPCR validation. Additionally, TUB-6 was selected as an internal reference gene. Primers for the differential genes assessed using an RT-qPCR analysis were designed using Beacon Designer 7.9 (Table 2). The total volume of the reaction system was 20 µL; the PCR thermal cycle used in this study was as follows: 94 °C (30 s), 94 °C (5 s), and 60 °C (30 s) for 40 cycles. RTqPCR was conducted according to the instructions provided by PerfectStart SYBR qPCR Supermix (TransGen Biotech, Beijing, China). Finally, the relative gene expression levels were calculated using the 2−∆∆CT method [29].

2.5. Combined Transcriptome and Metabolome Analysis

A transcriptomic and metabolomic integrated analysis was used to normalize and statistically analyze data at different biomolecular levels to establish relationships between these molecules. Pearson correlation was used to calculate the correlation between the two omics’ data and their p-values, and heat maps were plotted based on the above data, which showed that there was a significant positive correlation between metabolites and genes, or there was a significant negative correlation. Based on the metabolic pathway maps in the KEGG database (www.kegg.jp, accessed on 20 December 2021), genes and metabolites involved in the same metabolic pathway were analyzed by association; the genes and metabolites were shown by the different fills of the logo module.

3. Results

3.1. Changes in Agronomic Traits of Quinoa Plants During the Grouting Stage Under Various Phosphorus Concentrations

The amount of P2O5 applied to different treatments was 0 kg/hm2 (W0), 112.5 kg/hm2 (W1), 225.0 kg/hm2 (W2), 337.5 kg/hm2 (W3), 450.0 kg/hm2 (W4), and 562.5 kg/hm2 (W5), respectively. Samples were taken from the filling stage leaves of the quinoa plants 15 days after the flowering (repeated three times) to determine the plant height, leaf area, stem thickness, and other morphological indicators. At this time, the most significant morphological and phenotypic differences in the quinoa plants were evident. Among them, W0 exhibited the shortest plant height, while W4 possessed the tallest height; however, W5 exhibited a lower plant height than W4. Overall, this indicated that as the phosphorus content increased, the height of the quinoa plants increased; however, excessive-phosphorus application led to a downward trend in quinoa plant height. This same trend was observed for leaf area and stem thickness. Notably, the leaf area, plant height, and stem diameter of quinoa under different phosphorus fertilizer treatments were higher in the W4 group than those in all other groups, indicating that the W4 phosphorus fertilizer concentration may be the most suitable phosphorus fertilizer concentration for quinoa growth (Figure 1).

3.2. Metabolomic Response Mechanism of Quinoa Leaves to Phosphorus Fertilizer During Grouting Stage

3.2.1. Qualitative Analysis of Differential Metabolites

In the metabolomic analysis, we utilized three biological replicates for each phosphorus treatment group; specifically, by using the UPLC-MS/MS detection platform and a self-built database, we conducted qualitative and quantitative analyses across the different sample groups. A total of 179 metabolites were detected across the 18 quinoa leaf samples, most of which were carboxylic acids and derivatives (26.72%), Fatty Acyls (18.97%), and Organooxygen compounds (10.35) (Table S2; Figure 2C). Next, we overlapped the total ion flow diagrams (TIC diagrams) of different samples for the mass spectrometry detection and analysis; overall, the ion peak overlap of different substances in the repeated samples was determined to be relatively high (Figure 2A,B), indicating that the experimental process used in this study was stable, and the detection results were reliable.

3.2.2. PCA and Repeatability Evaluation of Samples

In the PCA score plot, the contribution rate of the first principal component was determined to be 14.84% and that of the second principal component was 12.72%, indicating that the two principal components reflected the main feature information of the test sample (Figure 3A). Overall, the 18 quinoa leaf samples exhibited clear separation in the two-dimensional graph, indicating that the data processing of each sample was reliable and that there were significant differences among the samples. Correlations between all metabolites were analyzed by calculating their corresponding Pearson or Spearman rank correlation coefficients. This correlation analysis of these quinoa samples revealed a strong correlation between the different phosphorus treatments. Specifically, the differential metabolite clustering heat maps (Figure 3B–F) clearly show the differences in metabolite expression between the quinoa seedlings treated with different phosphorus fertilizer concentrations.

3.2.3. Screening and Enrichment Analysis of Differential Metabolites

According to the differential metabolite volcano map (Figure 4A–E), eight differential metabolites were observed in W0 and W1, with two upregulated and six downregulated metabolites in W1. Further, 18 differential metabolites were identified between W2 and W0, with 11 upregulated and 7 downregulated metabolites in W2. Additionally, 25 differential metabolites were detected between W3 and W0, of which 16 were upregulated and 9 were downregulated in W3. Moreover, 39 differential metabolites were identified between W4 and W0, with 22 upregulated and 17 downregulated metabolites in W4. Finally, 40 differential metabolites were observed between W5 and W0, of which 13 were upregulated and 27 were downregulated in W5 (Tables S3–S7). The upregulation of sphingosine content in W2, W3, W4, and W5 indicates that sphingosine plays an important role in the excessive application of phosphorus fertilizer. Meanwhile, we used W1 as the control group; eight differential metabolites were identified between W2 and W1, all of which were upregulated in W2. Additionally, 17 differential metabolites were detected between W3 and W1, of which 15 were upregulated and 2 were downregulated in W3. Moreover, 35 differential metabolites were identified between W4 and W1, with 27 upregulated and 8 downregulated metabolites in W4. Finally, 29 differential metabolites were observed between W5 and W1, of which 16 were upregulated and 13 were downregulated in W5. Among these metabolites, tyrosol, meso-2,6-diaminoheptanedioate, kaempferol 3-O-beta-D-glucosyl-(1→2)-beta-D-glucoside, guanosine 2,3-cyclic phosphate, methylmalonic acid, and cyanin appeared simultaneously in low-phosphate (LP; W0) and HP (W3, W4, and W5) groups. Specifically, tyrosol and meso-2,6-diaminoheptanedioate were downregulated in W0 and upregulated in W5; kaempferol 3-O-beta-D-glucose-(1→2)–beta-D-glucoside was downregulated in W0 and upregulated in W3; guanosine 2,3-cyclic phosphate was downregulated in W0 and upregulated in W4 and W5; methylmalonic acid was upregulated in W0 and downregulated in W3, W4, and W5; and cyanin was upregulated in W0 and downregulated in W5 (Tables S8 and S9). As both tyrosol and meso-2,6-diaminoheptanedioate were determined to be downregulated in W0 and upregulated in W5, these two substances were implicated in the response of quinoa leaves to phosphorus fertilizer.
KEGG is the primary public database used to link differentially expressed factors to specific pathways. Using a KEGG analysis, differential metabolites in the W1 vs. W0 group were determined to be enriched in 18 metabolic pathways. Alternatively, differential metabolites in the W2 vs. W0 group were enriched in 24 metabolic pathways, with significant enrichment in the sphingolipid metabolism, purine metabolism, lysine biosynthesis, α-linolenic acid metabolism, and arachidonic acid metabolism pathway. Additionally, differential metabolites in the W3 vs. W0 group were significantly enriched in 32 metabolic pathways, including purine metabolism, sphingolipid metabolism, ABC transporters, lysine biosynthesis, α-linolenic acid metabolism, monobactam biosynthesis, and glycine, serine, and threonine metabolism. Further, differential metabolites in the W4 vs. W0 group were enriched in 46 metabolic pathways, with significant enrichment observed in lysine biosynthesis, α-linoleic acid metabolism, propanoate metabolism, flavone and flavonol biosynthesis, cysteine and methionine metabolism, pyrimidine metabolism, pyruvate metabolism, and arachidonic acid metabolism. Finally, differential metabolites in the W5 vs. W0 group were enriched in 43 metabolic pathways, with significant enrichment in ABC transporters, zeatin biosynthesis, oxidative phosphorylation, linoleic acid metabolism, and thiamine metabolism pathways. The degree of enrichment of these differential metabolites in each pathway is shown in the KEGG enrichment bar chart (Figure 5A–E) and KEGG enrichment bubble chart (Figure 5F–J). Most differentially expressed metabolites detected across all groups were enriched in phenylalanine metabolism, ABC transporters, purine metabolism, thiamine metabolism, and linoleic acid metabolism.

3.3. Transcriptomic Response of Quinoa Leaves to Phosphorus Fertilizer at Grouting Stage

3.3.1. Quality Analysis of Transcriptome Sequencing Data

After the filtering and quality evaluation of the raw data, the clean data of each sample reached 6 Gb; further, the corresponding base percentage score of Q30 was determined to be ≥91% for these data (Table S10). Then, the filtered reads were compared to the reference genome using HISAT2 software. The distribution of reads on all expressed genes showed a uniform distribution; additionally, the mapping ratio of sequences was >80%, indicating that the reference genome selection was appropriate and that there was no contamination in these experiments. Ultimately, this indicated that the sequencing results in this study have high accuracy and could be effectively used for a further analysis. The density plot produced indicated that gene abundance changed with differing expression levels within the 18 quinoa leaf samples; specifically, this demonstrated that the gene expression level (FPKM) in these samples was concentrated between log10−2 and log104 (Figure 6A). A correlation analysis was conducted between each group of replicates; the corresponding box plot (Figure 6B) and correlation heatmap (Figure 6C) provided evidence for good biological repeatability within each group and indicated notable differences between the groups. According to the PCA diagram (Figure 6D), the contribution rates of the first and second principal components were 75% and 12%, respectively. The similarity among samples W0, W1, and W4 was determined to be high. Therefore, W4 was used as the high-phosphorus (HP) group, which was used to preliminarily evaluate the response mechanisms of quinoa seedlings to phosphorus fertilizer. Functional annotation was performed on the genes detected in this experiment using databases, such as eggNOG, cellular components (CCs), biological processes (BPs), Gene Ontology (GO), Swissprot, molecular functions (MFs), Non-Redundant Protein Sequence Database (NR), and KEGG. Overall, eggNOG annotated 41,797 genes, CC annotated 20,773 genes, BP annotated 21,343 genes, GO annotated 22,615 genes, Swissprot annotated 35,186 genes, MF annotated 19,806 genes, NR annotated 49,168 genes, and KEGG annotated 15,677 genes.

3.3.2. Differential Gene Expression Analysis

Next, we used DESeq to evaluate the differential expression of genes. The corresponding statistical chart of the differentially expressed gene (DEG) analysis indicated the presence of DEGs between the different comparison groups. Among them, 1152 DEGs were identified in the W0 vs. W1 comparison group, of which 515 were downregulated and 637 were upregulated; 3481 DEGs were identified in the W0 vs. W2 comparison group, of which 1651 were downregulated and 1830 were upregulated; 5562 DEGs were identified in the W0 vs. W3 comparison group, of which 2247 were downregulated and 3315 were upregulated; 2789 DEGs were identified in the W0 vs. W4 comparison group, of which 972 were downregulated and 1817 were upregulated; and 2717 DEGs were identified in the W1 vs. W5 comparison group, of which 1192 were downregulated and 1525 were upregulated (Figure 7A). Further, the volcanic map of these DEGs indicated that the number of upregulated DEGs in W1, W2, W3, W4, and W5, when compared with those in W0, was higher than the number of downregulated genes; this indicated that these DEGs are responsive differential genes that may play an important role in the response of quinoa leaves to different phosphorus fertilizers (Figure 7B–F).
To study the expression patterns of different sample genes under phosphorus-associated stress, we used the R language Pheatmap software package to perform a two-way cluster analysis (Figure 7G). This analysis included assessing the connections between different genes and samples from all comparison groups. Overall, we centralized and standardized the FPKM of genes and divided these genes into different clusters, grouping genes with similar gene expression trends together. The corresponding trend charts (Figure 7H) indicated that the DEGs could be categorized into nine distinct clusters; notably, genes in the eighth and ninth clusters exhibited significant overall changes in expression in the leaves across different phosphorus levels during the filling stage. These findings suggest that these DEGs could serve as potential markers for phosphorus regulation in quinoa.

3.3.3. Differential Gene KEGG Annotation and Enrichment Analysis

Mapping of DEGs across various metabolic pathways to the KEGG database revealed that these DEGs were enriched in metabolic pathways related to phosphorus metabolism. According to the KEGG enrichment bubble diagram, the metabolic pathways with more DEG enrichment in the W0 vs. W1 comparison group were plant hormone signal transduction, phenylpropanoid biosynthesis, flavonoid biosynthesis, plant circadian rhythms, and valine, leucine, and isoleucine degradation. Alternatively, the metabolic pathways with more DEG enrichment in the W0 vs. W2 comparison group were determined by plant–pathogen interaction; photosynthesis antenna proteins; phenylpropanoid biosynthesis; starch and sucrose metabolism; and cutin, suberine, and wax biosynthesis. The metabolic pathways with more DEG enrichment in the W0 vs. W3 comparison group were established as plant–pathogen interactions, plant hormone signal transduction, MAPK signaling pathway—plant, alpha-linolenic acid metabolism, and pentose and glucuronate interconversions. Additionally, the metabolic pathways with more DEG enrichment in the W0 vs. W4 comparison group included plant–pathogen interactions, MAPK signaling pathway—plant, starch and sucrose metabolism, zeatin biosynthesis, and alpha-linolenic acid metabolism. Finally, the metabolic pathways with more DEG enrichment in the W0 vs. W5 comparison group were determined to be plant–pathogen interaction, pentose and glucuronate interconversions, MAPK signaling pathway—plant, phenylpropanoid biosynthesis, and circadian rhythm—plant (Figure 8A–E).

3.3.4. GO Annotation and Enrichment Analysis

Next, all DEGs were subjected to a GO enrichment analysis using topGO. The GO enrichment analysis results of the DEGs were categorized according to MF, BP, and CC. The top 10 GO term entries with the most significant enrichment were identified in each GO classification (Figure 9A–E). In the W0 vs. W1 comparison group, 3174 significantly enriched GO functions were obtained, including 252 CCs, 741 MFs, and 2181 BPs. In the W0 vs. W2 comparison group, 5597 significantly enriched GO functions were obtained, including 473 CCs, 1287 MFs, and 3837 BPs. Further, in the W0 vs. W3 comparison group, 6973 significantly enriched GO functions were obtained, including 618 CCs, 1533 MFs, and 4822 BPs. Additionally, in the W0 vs. W4 comparison group, 4832 significantly enriched GO functions were obtained, including 392 CCs, 1088 MFs, and 3352 BPs. Finally, in the W0 vs. W5 comparison group, 4579 significantly enriched GO functions were obtained, including 363 CCs, 1075 MFs, and 3141 BPs.
According to the GO enrichment column chart of DEGs, the BP with the most significant enrichment in the W0 vs. W1 comparative group was mainly enriched in responses to stimuli, responses to chemical substances, responses to oxygen-containing compounds, responses to abiotic stimuli, and responses to organic substances. In the W0 vs. W2 comparative group, DEGs were primarily enriched in responses to stimuli, defense response, response to chitin, response to stress, and cell wall organization or biogenesis. In the W0 vs. W3 comparative group, DEGs were mainly enriched in responses to stimuli, cell communication, signal transduction, signaling, and defense response. In the W0 vs. W4 comparative group, DEGs were primarily enriched in responses to stimuli, responses to chitin, responses to organonitrogen compounds, responses to oxygen-containing compounds, and response to stress. Finally, in the W0 vs. W5 comparative group, DEGs were predominantly enriched in responses to chitin, responses to organonitrogen compounds, defense response, response to stress, and response to an external stimulus.

3.3.5. Transcription Factor Family Analysis

In this study, we also identified the transcription factors of all DEGs; we established that the response of quinoa leaves to phosphorus fertilizer at the grouting stage was regulated by various transcription factor families, such as bHLH, NAC, MYB-related, FAR1, ERF, C2H2, GRAS, WRKY, M-type-MADS, bZIP, B3, ZF-HD, MYB, C3H, and G2-like families. Specifically, 878, 620, and 572 DEGs were observed to be regulated by the bHLH, NAC, and MYB transcription factor families, respectively (Figure 10). Moreover, the expression levels of the genes regulated by these transcription factors were determined to be relatively high. Considering the important roles that these transcription families play in the response to environmental stress and internal physiological responses in quinoa plants, these factors should be investigated further.

3.3.6. RT-qPCR Validation

To determine the authenticity and reliability of the transcriptome data and the differential expression levels of candidate genes, we randomly selected 10 genes for RT-qPCR verification. The expression of genes gene-LOC110707704, gene-LOC110700414, gene-LOC110727076, gene-LOC110682410, gene-LOC110727927, and gene-LOC110685810 was consistent with the RNA-seq data, indicating that the transcriptome sequencing results in this study were reliable (Figure 11).

3.4. The Integrative Metabolomic and Transcriptomic Analysis of the Response Mechanism of Quinoa Leaves to Different Phosphorus Concentrations at the Grouting Stage

After phenotypic, metabolomic, and transcriptomic analyses, all data of the plants under W5 treatment were established to be significantly different from those of the normal plants. Therefore, we used W5 as a high-phosphorus treatment and W0 as a phosphorus-deficiency treatment; we then used the transcriptomic and metabolomic data to jointly analyze the response mechanism of quinoa leaves to phosphorus fertilizer at the grouting stage. The association heat map shows whether the relationship between metabolites and genes is positive or negative (Figure 12). Simultaneous mapping of metabolites and genes to the KEGG pathway revealed that both genes and metabolites were enriched for flavonoid biosynthesis in both the W0 vs. W1 group and the W5 vs. W1 group.
In the flavonoid biosynthesis pathway (Figure 13), the genes LOC110727184, LOC110724462, and LOC110727183 in chalcone synthesis (CHS); LOC110694697 and LOC110724781 in F3H; and LOC110714529 and LOC110725334 in FLS were upregulated under LP conditions and HP conditions. The gene LOC110724467 in chalcone synthesis (CHS) was downregulated under LP conditions and HP conditions. This indicates that these eight genes play a crucial role in the response of quinoa to high- or low-phosphorus stress. The accumulation upregulated kaempferol, Luteolin, Dihydrokaempferol, and Cyanidin under low-phosphorus treatment, and interestingly, the downregulation of kaempferol, Luteolin, Dihydrokaempferol, and Cyanidin under high-phosphorus conditions suggests that high phosphorus inhibits the synthesis of substances in the flavonoid biosynthetic pathway such as kaempferol.

4. Discussion

Phosphorus is one of the “three elements” of crop nutrition, and plays an important role in crop growth and development, as well as in the formation of quality and yield. Phosphorus is a component of ATP (adenosine triphosphate), which constitutes a part of nucleic acid molecules such as DNA and RNA, and all the physiological activities in the process of crop growth and development are carried out under the expression of the genetic material DNA and RNA information, and phosphorus has an extremely important role in the respiration, photosynthesis, and energy conversion of plants [30]. The phosphorus concentration in the soil usually does not reach the minimum requirement of phosphorus for plants, so the application of phosphorus fertilizers is necessary to increase crop yields. With adequate supply of phosphorus, the ability of root mass production and deep penetration into the soil can be improved, which will enable plants to absorb more water from deeper soil and alleviate drought. However, phosphorus in the soil can easily react with some cations in the soil and produce precipitation or be adsorbed and fixed by the soil, thus reducing the utilization rate of phosphorus fertilizer [31].
The present study combined transcriptomic and metabolomic data to analyze the response mechanism of quinoa leaves to phosphorus fertilizer at the grouting stage, there by providing a theoretical basis for improving the yield and quality of quinoa and establishing sustainable use of soil [32]. Previously, we conducted a preliminary analysis of the response mechanism of quinoa to phosphorus fertilizers. Specifically, we assessed three concentrations of phosphorus fertilizer, namely 0 kg/hm2, 112.5 kg/hm2, and 337.5 kg/hm2. This prior study established that quinoa seedlings grew better under the 337.5 kg/hm2 phosphorus fertilizer treatment, with plants being shorter under LP conditions [32]. Therefore, the present study utilized a gradient of phosphorus fertilizer, from 0 kg/hm2 phosphorus treatment to treatment containing five times the standard phosphorus treatment (562.5 kg/hm2). Overall, we established that quinoa plants grew more effectively when treated with 450.0 kg/hm2 phosphorus fertilizer; however, treatment with 562.5 kg/hm2 phosphorus fertilizer resulted in delayed seedling development and a significant inhibition of growth. Therefore, we used W5 as a representative treatment for high-phosphorus stress in quinoa.
In recent years, many researchers have conducted relevant research on plant responses to phosphorus stress using omics data. For example, Ding et al. determined that the synthesis of flavonoids and phosphorylated metabolites was significantly reduced in the absence and presence of excess phosphorus [33]. Consistent with the results of the present study, transcriptomic and metabolomic correlation analyses revealed that metabolite synthesis in the flavonoid biosynthesis pathway was inhibited in the presence of a phosphorus fertilization excess. Further, Wang et al. conducted a metabolomic analysis of quinoa, identifying the differential expression of 177 lipid metabolites, 172 flavonoid metabolites, and 166 phe.
Zhu et al. determined that plants can release phosphatidylinositol from membrane phospholipids in response to phosphorus-deficiency stress [34]. Phosphorus deficiency reduces glutamate levels and increases glutamine content; in contrast, excessive phosphorus leads to a decrease in glutamine content [33]. Muhammad et al. found that the differential metabolites observed between low-phosphorus stress and normal fertilization were primarily sugars, amino acids, and organic acids [35]. Additionally, Sun et al. determined that apple DEGs were mainly enriched in flavonoid biosynthesis, phenylpropane biosynthesis, and ATP-binding cassette transporters under low-phosphorus stress conditions [36]. The metabolites and enrichment pathways identified in the above studies were similar to those identified in the present study. Specifically, our findings revealed that the differential metabolites across the different phosphorus treatments primarily included lipids, amino acids, carbohydrates, and nucleotides; further, we determined that these differential metabolites were mainly enriched in pathways such as phenylalanine metabolism, α-linoleic acid metabolism, and flavonoid biosynthesis.
Plants cannot directly absorb large amounts of organic phosphorus but instead absorb inorganic phosphate from the soil through their roots [14]. Wang et al. proposed that quinoa seedlings respond to varying phosphorus environments by regulating changes in metabolites and gene expression within pathways related to glycerides and phospholipids [31]. Wu et al. combined metabolome and transcriptome analyses, revealing that phosphorus deficiency can significantly inhibit plant growth and development. Additionally, consistent with the results of the current study, Wu et al. determined that phenylpropane biosynthesis pathways are closely related to plant responses to phosphorus-deficiency stress [37]. Under phosphorus-deficient conditions, anthocyanin accumulation can be induced in plant leaves. Additionally, under phosphorus-deficient conditions in the absence of inositol polyphosphate, PAP1 is released from SPX, leading to the activation of anthocyanin biosynthesis [38,39]. Our current study demonstrated that cyanin accumulates under low-phosphorus conditions, while its expression is downregulated under high-phosphorus conditions. In our metabolomics study, we noted that cyanin (Cyanidin) accumulates under low-phosphorus conditions but downregulates its expression under high-phosphorus conditions. We postulate that cyanin has a certain effect on the response of resveratrol to low-phosphorus stress. But according to reports, quinoa cannot synthesize anthocyanins and they are replaced by betalains, with betalains and anthocyanins being mutually exclusive. Betalains are special water column segments that exist only in Caryophylllaceae plants. Compared to anthocyanins and carotenoids, the synthesis and regulation mechanism of betalains is simpler [40].
Peroxidases, being oxidoreductases, are involved in various cellular processes, including lipid peroxidation, and are often found in peroxisomes. Nonetheless, peroxisomes also contain enzymes related to phospholipid synthesis and, therefore, also participate in lipid synthesis. Cocozza et al. evaluated the potential effects of phenylpropanes and peroxidase on the adaptability of Artemisia donovaria under limited- and excessive-phosphorus conditions [41].
Previous studies have shown that under phosphorus stress, gene families, such as MYB, WRKY, ARF, and BHLH, are induced, which aligns with the findings of the present study. Our results indicated that the primary transcription factor families associated with the response of quinoa leaves to phosphorus fertilizer during the grouting stage include bHLH, NAC, WRKY, bZIP, B3, ZF-HD, and MYB, alongside other gene families [42,43,44,45]. WRKY transcription factors are a family of regulatory genes linked to various non-biotic stresses in plants; however, their role in phosphorus deficiency has rarely been reported. Nonetheless, Liu et al. previously demonstrated that soybean-related genes belonging to the WRKY family are involved in the regulation of soybean phosphorus-deficiency tolerance [46]. Therefore, our results further demonstrated that the WRKY transcription family is involved in the regulation of plant responses to phosphorus stress.
This study analyzed the response of quinoa (Dianli-2912) to phosphorus fertilizer based on transcriptomics and metabolomics, which is not representative of all quinoa lines, and the molecular mechanism of quinoa’s response to phosphorus fertilizer in the future still needs to be further explored. In agricultural production, in order to improve crop yield, the application of large amounts of nitrogen and phosphorus fertilizers is mainly relied on to achieve increased yields, but too much phosphorus fertilizer will cause soil acidification, water eutrophication, and other ecosystems to be damaged, and is not conducive to the sustainable development of agriculture [47,48]. The application of phosphorus fertilizer should be rationally combined with the soil demand, and at the same time, the development of new fertilization technology based on the principle of the inter-root regulation of phosphorus fertilizer should be carried out to further improve the efficiency of phosphorus fertilizer utilization and crop yield, which is conducive to the realization of the green and sustainable development of agriculture.

5. Conclusions

In this study, we conducted metabolomic and transcriptomic analyses of quinoa leaves at the grouting stage under different phosphorus fertilizer levels. Overall, different phosphorus conditions had a significant impact on the phenotypic characteristics of quinoa. Additionally, we identified the changes in the regulation of metabolites and genes in flavonoid biosynthesis pathways, which explained the corresponding response mechanisms of these quinoa leaves under varied phosphorus conditions. Different phosphorus conditions had specific impacts on the differentially expressed metabolites and genes of quinoa. Ultimately, we established that these pathways can adapt and respond to different phosphorus fertilizer environments. This study provides theoretical references for improving quinoa yield and achieving sustainable use of phosphorus fertilizer.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/agronomy14112661/s1, Figure S1: Appraisal chart of daughter ions (mz); Table S1: Information on internal standards. Table S2: Supplementary table for the determined metabolites. Table S3: DAMs in W0 vs.W1. Table S4: DAMs in W0 vs.W2. Table S5: DAMs in W0 vs.W3. Table S6: DAMs in W0 vs.W4. Table S7: DAMs in W0 vs.W5. Table S8: DEMs in the group of LP vs. CK(W0 VS W1). Table S9: DEMs in the group of HP vs. CK(W5 VS W1). Table S10: Statistics of disembarking data.

Author Contributions

Writing—original draft preparation, H.W., Q.W., J.L. and P.Z.; writing—review and editing, H.L., X.L., L.L. and H.X.; methodology, Q.W. and H.X.; validation, H.L., X.L. and H.W.; formal analysis, J.L.; investigation, H.L. and L.L.; data curation, J.L. and P.Z.; supervision, H.L. and H.X.; project administration, P.Q.; funding acquisition, P.Q.; resources, P.Q.; conceptualization, P.Q. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Yunnan Province Academician Workstation (202405AF140012), the Yunnan Province’s “Xing Dian Talent Support Program” (XDYC-CYCX-2022-0031) and the Yunnan Province Expert Workstation (202205AF150001).

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Materials; further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Bhargava, A.; Shukla, S.; Ohri, D. Chenopodium quinoa—An Indian perspective. Ind. Crops Prod. 2006, 23, 73–87. [Google Scholar] [CrossRef]
  2. Roman, V.J.; Toom, L.; Gamiz, C.C.; Pijl, N.; Linden, C. Differential responses to salt stress in ion dynamics, growth and seed yield of European quinoa varieties. Environ. Exp. Bot. 2020, 177, 104146. [Google Scholar] [CrossRef]
  3. Adolf, V.I.; Shabala, S.; Andersen, M.N.; Razzaghi, F.; Jacobsen, S.E. Varietal differences of quinoa’s tolerance to saline conditions. Plant Soil 2012, 357, 117–129. [Google Scholar] [CrossRef]
  4. Hariadi, Y.; Marandon, K.; Tian, Y.; Jacobsen, S.E.; Shabala, S. Ionic and osmotic relations in quinoa (Chenopodium quinoa Willd.) plants grown at various salinity levels. J. Exp. Bot. 2011, 62, 185–193. [Google Scholar] [CrossRef] [PubMed]
  5. Vega-Gálvez, A.; Miranda, M.; Vergara, J.; Uribe, E.; Puente, L.; Martínez, E.A. Nutrition facts and functional potential of quinoa (Chenopodium quinoa willd.), an ancient Andean grain: A review. J. Sci. Food Agric. 2010, 90, 2541–2547. [Google Scholar] [CrossRef]
  6. Jacobsen, S.E. The worldwide potential for quinoa (Chenopodium quinoa Willd.). Food Rev. Int. 2003, 19, 167–177. [Google Scholar] [CrossRef]
  7. Escuredo, O.; González Martín, M.I.; Wells Moncada, G.; Fischer, S.; Hernández Hierro, J.M. Amino acid profile of the quinoa (Chenopodium quinoa Willd.) using near infrared spectroscopy and chemometric techniques. J. Cereal Sci. 2014, 60, 67–74. [Google Scholar] [CrossRef]
  8. FAO. Food and Agriculture Organization of the United States/World Health Organization/United Nations University, Energy and Protein Requirements. Report of a Joint FAO/WHO/UNU Meeting; World Health Organization: Geneva, Switzerland, 1985. [Google Scholar]
  9. Repo-Carrasco, R.; Espinoza, C.; Jacobsen, S.-E. Nutritional Value and Use of the Andean Crops Quinoa (Chenopodium quinoa) and kañiwa (Chenopodium pallidicaule). Food Rev. Int. 2003, 19, 179–189. [Google Scholar] [CrossRef]
  10. Dakhili, S.; Abdolalizadeh, L.; Hosseini, S.M.; Shojaee-Aliabadi, S.; Mirmoghtadaie, L. Quinoa protein: Composition, structure and functional properties. Food Chem. 2019, 299, 125–161, Corrigendum in Food Chem. 2019, 310, 125318. [Google Scholar] [CrossRef]
  11. Yang, X.; Hou, Z.; Xue, P. Antibacterial activity and mechanism of action saponins from chenopodium quinoa willd. husks against foodborne pathogenic bacteria. Ind. Crops Prod. 2020, 149, 112350. [Google Scholar] [CrossRef]
  12. Jin, H.M.; Wei, P. Anti-Fatigue Properties of Tartary Buckwheat Extracts in Mice. Int. J. Mol. Sci. 2011, 12, 4770–4780. [Google Scholar] [CrossRef] [PubMed]
  13. Guo, H.; Hao, Y.; Richel, A.; Everaert, N.; Ren, G. Antihypertensive effect of quinoa protein under simulated gastrointestinal digestion and peptide characterization. J. Sci. Food Agric. 2020, 100, 5569–5576. [Google Scholar] [CrossRef] [PubMed]
  14. Mimura, T. Regulation of Phosphate Transport and Homeostasis in Plant Cells. Int. Rev. Cytol. 1999, 191, 149–200. [Google Scholar] [CrossRef]
  15. Shen, J.; Yuan, L.; Zhang, J.; Li, H.; Bai, Z.; Chen, X. Phosphorus dynamics: From soil to plant. Plant Physiol. 2011, 156, 997–1005. [Google Scholar] [CrossRef]
  16. Zhou, J.; Shi, W.; Pan, K.; Ying, Y.; Sun, C. Effect of low phosphorus stress on growth and nutrient physiology of Phyllostachys edulis seedlings. J. Zhejiang AF Univ. 2022, 39, 1010–1017. [Google Scholar] [CrossRef]
  17. Naureen, Z.; Sham, A.; Al Ashram, H.; Gilani, S.A.; Al Gheilani, S.; Mabood, F.; Hussain, J.; Al Harrasi, A.; AbuQamar, S.F. Effect of phosphate nutrition on growth, physiology and phosphate transporter expression of cucumber seedlings. Plant Physiol. Biochem. 2018, 127, 211–222. [Google Scholar] [CrossRef]
  18. Liao, Z.X.; Li, X.H.; Xue, Y.B.; Yang, N.D.; Wu, Z.W.; Liu, Y. Effects of phosphorus on growth and development of soybean seedlings. Key Eng. Mater. 2022, 905, 353–358. [Google Scholar] [CrossRef]
  19. Juszczuk, I.; Malusà, E.; Rychter, A.M. Oxidative stress during phosphate deficiency in roots of bean plants (phaseolus vulgaris L.). J. Plant Physiol. 2001, 158, 1299–1305. [Google Scholar] [CrossRef]
  20. Herná, G.; Ramírez, M.; Valdés-López, O.; Tesfaye, M.; Graham, M.A.; Czechowski, T.; Schlereth, A.; Wandrey, M.; Erban, A.; Cheung, F.; et al. Focus issue on legume biology:phosphorus stress in common bean:root transcript and metabolic responses. Plant Physiol. 2007, 144, 752–767. [Google Scholar] [CrossRef]
  21. Verbruggen, N. How do plants respond to nutrient shortage by biomass allocation? Trends Plant Sci. 2010, 11, 610–617. [Google Scholar] [CrossRef]
  22. Dunn, W.; Broadhurst, D.; Begley, P.; Zelena, E.; Francis-McIntyre, S.; Anderson, N.; Brown, M.; Knowles, J.D.; Halsall, A.; Haselden, J.N.; et al. Procedures for large-scale metabolic profiling of serum and plasma using gas chromatography and liquid chromatography coupled to mass spectrometry. Nat. Protoc. 2011, 6, 1060–1083. [Google Scholar] [CrossRef] [PubMed]
  23. Want, E.; Wilson, I.; Gika, H.; Theodoridis, G.; Plumb, R.S.; Shockcor, J.; Holmes, E.; Nicholson, J.K. Global metabolic profiling procedures for urine using UPLC-MS. Nat. Protoc. 2010, 5, 1005–1018. [Google Scholar] [CrossRef] [PubMed]
  24. Thévenot, E.A.; Roux, A.; Xu, Y.; Ezan, E.; Junot, C. Analysis of the human adult urinary metabolome variations with age, body mass index, and gender by implementing a comprehensive workflow for univariate and opls statistical analyses. J. Proteome Res. 2015, 14, 3322–3335. [Google Scholar] [CrossRef] [PubMed]
  25. Trygg, J.; Wold, S. Orthogonal projections to latent structures (O-PLS). J. Chemom. 2002, 16, 119–128. [Google Scholar] [CrossRef]
  26. Kieffer, D.A.; Piccolo, B.D.; Vaziri, N.D.; Liu, S.; Lau, W.L.; Khazaeli, M. Resistant starch alters gut microbiome and metabolomic profiles concurrent with amelioration of chronic kidney disease in rats. Am. J. Physiol. Renal Physiol. 2016, 310, F857–F871. [Google Scholar] [CrossRef]
  27. Saude, E.J.; Skappak, C.D.; Regush, S.; Cook, K.; Ben-Zvi, A.; Becker, A.; Saude, R.E.J.; Skappak, C.D.; Regush, S.; Cook, K.; et al. Metabolomic profiling of asthma: Diagnostic utility of urine nuclear magnetic resonance spectroscopy. J. Allergy Clin. Immunol. 2011, 127, 757–764. [Google Scholar] [CrossRef]
  28. Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139. [Google Scholar] [CrossRef]
  29. Livak, K.J.; Schmittgen, T. Analysis of relative gene expression data using real-time quantitative pcr and the 2-∆∆ct method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
  30. Sikma, M.; Ikawa, H.; Heusinkveld, B.G.; Yoshimoto, M.; Hasegawa, T.; Haar, L.T.G.; Anten, N.P.R.; Nakamura, H.; de Arellano, J.V.-G.; Sakai, H.; et al. Quantifying the feedback between rice architecture, physiology, and microclimate under current and future CO2 conditions. J. Geophys. Res. Biogeosci. 2020, 125, e2019JG005452. [Google Scholar] [CrossRef]
  31. Parihar, M.; Meena, V.S.; Mishra, P.K.; Rakshit, A.; Choudhary, M.; Yadav, R.P.; Rana, K.; Bisht, J.K. Arbuscular mycorrhiza: A viable strategy for soil nutrient loss reduction. Arch Microbiol. 2019, 201, 723–735. [Google Scholar] [CrossRef]
  32. Wang, Q.; Guo, Y.; Huang, T.; Zhang, X.; Zhang, P.; Xie, H.; Liu, J.; Li, L.; Kong, Z.; Qin, P. Transcriptome and Metabolome Analyses Revealed the Response Mechanism of Quinoa Seedlings to Different Phosphorus Stresses. Int. J. Mol. Sci. 2022, 23, 4704. [Google Scholar] [CrossRef] [PubMed]
  33. Ding, Z.; Jia, S.; Wang, Y.; Xiao, J.; Zhang, Y. Phosphate stresses affect ionome and metabolome in tea plants. Plant Physiol. Biochem. 2017, 129, 30. [Google Scholar] [CrossRef] [PubMed]
  34. Zhu, S.; Liang, C.; Tian, J.; Xue, Y. Advances in Plant Lipid Metabolism Responses to Phosphate Scarcity. Plants 2022, 11, 2238. [Google Scholar] [CrossRef]
  35. Muhammad, I.I.; Nor, S.; Abdullah, A.; Saud, H.M.; Gori, A. The dynamic responses of oil palm leaf and root metabolome to phosphorus deficiency. Metabolites. 2021, 11, 217. [Google Scholar] [CrossRef]
  36. Sun, T.; Zhang, J.; Zhang, Q.; Li, X.; Li, M.; Yang, Y.; Zhou, J.; Wei, Q.; Zhou, B. Transcriptome and metabolome analyses revealed the response mechanism of apple to different phosphorus stresses. Plant Physiol. Biochem. 2021, 167, 639–650. [Google Scholar] [CrossRef]
  37. Wu, Q.; Yang, L.; Liang, H.; Yin, L.; Chen, D.; Shen, P. Integrated analyses reveal the response of peanut to phosphorus deficiency on phenotype, transcriptome and metabolome. BMC Plant Biol. 2022, 22, 524. [Google Scholar] [CrossRef] [PubMed]
  38. Peng, Z.; Tian, J.; Luo, R.; Kang, Y.; Lu, Y.; Hu, Y.; Liu, N.; Zhang, J.; Cheng, H.; Niu, S.; et al. Mir399d and epigenetic modification comodulate anthocyanin accumulation in malus leaves suffering from phosphorus deficiency. Plant Cell Environ. 2020, 43, 1148–1159. [Google Scholar] [CrossRef]
  39. He, Y.; Zhang, X.; Li, L.; Sun, Z.; Li, J.; Chen, X.; Hong, G. Spx4 interacts with both phr1 and pap1 to regulate critical steps in phosphorus-status-dependent anthocyanin biosynthesis. New Phytol. 2021, 230, 205–217. [Google Scholar] [CrossRef]
  40. Zhao, X.; Zhang, Y.; Long, T.; Wang, S.; Yang, J. Regulation Mechanism of Plant Pigments Biosynthesis: Anthocyanins, Carotenoids, and Betalains. Metabolites 2022, 12, 871. [Google Scholar] [CrossRef]
  41. Cocozza, C.; Bartolini, P.; Brunetti, C.; Miozzi, L.; Pignattelli, S.; Podda, A.; Scippa, G.S.; Trupiano, D.; Rotunno, S.; Brilli, F.; et al. Modulation of class III peroxidase pathways and phenylpropanoids in Arundo donax under salt and phosphorus stress. Plant Physiol. Biochem. 2022, 183, 151–159. [Google Scholar] [CrossRef]
  42. Dai, X.; Wang, Y.; Wen-Hao, Z. Oswrky74, a wrky transcription factor, modulates tolerance to phosphate starvation in rice. J. Exp. Bot. 2016, 67, 947–960. [Google Scholar] [CrossRef] [PubMed]
  43. Nilsson, L.; Müller, R.; Nielsen, T.H. Increased expression of the myb-related transcription factor, phr1, leads to enhanced phosphate uptake in arabidopsis thaliana. Plant Cell Environ. 2010, 30, 1499–1512. [Google Scholar] [CrossRef] [PubMed]
  44. Filiz, E.; Vatansever, R.; Ozyigit, I.I. Dissecting a co-expression network of basic helix-loop-helix (bhlh) genes from phosphate (pi)-starved soybean (Glycine max). Plant Gene 2016, 9, 19–25. [Google Scholar] [CrossRef]
  45. Wang, S.; Zhang, S.; Sun, C.; Xu, Y.; Chen, Y.; Yu, C.; Qian, Q.; Jiang, D.-A.; Qi, Y. Auxin response factor (OsARF12), a novel regulator for phosphate homeostasis in rice (Oryza sativa). New Phytol. 2014, 201, 91–103. [Google Scholar] [CrossRef]
  46. Liu, X.; Yang, Y.; Wang, R.; Cui, R.; Xu, H.; Sun, C.; Wang, J.; Zhang, H.; Chen, H.; Zhang, D. GmWRKY46, a WRKY transcription factor, negatively regulates phosphorus tolerance primarily through modifying root morphology in soybean. Plant Sci. 2021, 315, 111148. [Google Scholar] [CrossRef]
  47. Le Noë, J.; Garnier, J.; Billen, G. Phosphorus management in cropping systems of the Paris Basin: From farm to regional scale. J. Environ. Manag. 2018, 205, 18–28. [Google Scholar] [CrossRef]
  48. Mahapatra, D.M.; Satapathy, K.C.; Panda, B. Biofertilizers and nanofertilizers for sustainable agriculture: Phycoprospects and challenges. Sci. Total Environ. 2022, 803, 14999. [Google Scholar] [CrossRef]
Agronomy 14 02661 g001
Figure 1. Agronomic characteristics of quinoa under treatment with different doses of P2O5.
Figure 1. Agronomic characteristics of quinoa under treatment with different doses of P2O5.
Agronomy 14 02661 g001
Agronomy 14 02661 g002
Figure 2. (A,B) Base peak chromatogram (BPC) of typical quinoa leaf samples from each group treated with different doses of P2O5. (C) Differential metabolite classification.
Figure 2. (A,B) Base peak chromatogram (BPC) of typical quinoa leaf samples from each group treated with different doses of P2O5. (C) Differential metabolite classification.
Agronomy 14 02661 g002
Agronomy 14 02661 g003
Figure 3. (A) PCA score chart of quality spectrum data from each group of quinoa samples; (B–F) clustering heat maps of differentially expressed metabolites in these quinoa leaves treated with different doses of P2O5. (B) W0 vs. W1; (C) W0 vs. W2; (D) W0 vs. W3; (E) W0 vs. W4; (F) W0 vs. W5.
Figure 3. (A) PCA score chart of quality spectrum data from each group of quinoa samples; (B–F) clustering heat maps of differentially expressed metabolites in these quinoa leaves treated with different doses of P2O5. (B) W0 vs. W1; (C) W0 vs. W2; (D) W0 vs. W3; (E) W0 vs. W4; (F) W0 vs. W5.
Agronomy 14 02661 g003
Agronomy 14 02661 g004aAgronomy 14 02661 g004b
Figure 4. Volcanic maps of the differential metabolites in quinoa leaves treated with different doses of P2O5. (A) W0 vs. W1; (B) W0 vs. W2; (C) W0 vs. W3; (D) W0 vs. W4; (E) W0 vs. W5.
Figure 4. Volcanic maps of the differential metabolites in quinoa leaves treated with different doses of P2O5. (A) W0 vs. W1; (B) W0 vs. W2; (C) W0 vs. W3; (D) W0 vs. W4; (E) W0 vs. W5.
Agronomy 14 02661 g004aAgronomy 14 02661 g004b
Agronomy 14 02661 g005aAgronomy 14 02661 g005b
Figure 5. (A–E) Bubble charts and (F–J) bar charts of metabolic pathway-influencing factors in quinoa leaves treated with different doses of P2O5. (A) W0 vs. W1; (B) W0 vs. W2; (C) W0 vs. W3; (D) W0 vs. W4; (E) W0 vs. W5; (F) W0 vs. W1; (G) W0 vs. W2; (H) W0 vs. W3; (I) W0 vs. W4; (J) W0 vs. W5.
Figure 5. (A–E) Bubble charts and (F–J) bar charts of metabolic pathway-influencing factors in quinoa leaves treated with different doses of P2O5. (A) W0 vs. W1; (B) W0 vs. W2; (C) W0 vs. W3; (D) W0 vs. W4; (E) W0 vs. W5; (F) W0 vs. W1; (G) W0 vs. W2; (H) W0 vs. W3; (I) W0 vs. W4; (J) W0 vs. W5.
Agronomy 14 02661 g005aAgronomy 14 02661 g005b
Agronomy 14 02661 g006
Figure 6. (A) Distribution of FPKM value expression density; (B) violin diagram of FPKM value expression; (C) sample correlation graph; and (D) PCA score chart of quinoa leaves treated with different doses of P2O5.
Figure 6. (A) Distribution of FPKM value expression density; (B) violin diagram of FPKM value expression; (C) sample correlation graph; and (D) PCA score chart of quinoa leaves treated with different doses of P2O5.
Agronomy 14 02661 g006
Agronomy 14 02661 g007aAgronomy 14 02661 g007bAgronomy 14 02661 g007c
Figure 7. The differential gene expression analysis of quinoa leaves treated with different doses of P2O5. (A) The statistical map of the expression difference analysis; (B–F) volcano maps of differentially expressed genes; (G) the clustering heat map of differentially expressed genes; and (H) the trend analysis chart. In the trend analysis chart, the gray line indicates the expression patterns of genes in each cluster and the blue line represents the average expression levels of all genes in the cluster. (B) W0 vs. W1; (C) W0 vs. W2; (D) W0 vs. W3; (E) W0 vs. W4; (F) W0 vs. W5.
Figure 7. The differential gene expression analysis of quinoa leaves treated with different doses of P2O5. (A) The statistical map of the expression difference analysis; (B–F) volcano maps of differentially expressed genes; (G) the clustering heat map of differentially expressed genes; and (H) the trend analysis chart. In the trend analysis chart, the gray line indicates the expression patterns of genes in each cluster and the blue line represents the average expression levels of all genes in the cluster. (B) W0 vs. W1; (C) W0 vs. W2; (D) W0 vs. W3; (E) W0 vs. W4; (F) W0 vs. W5.
Agronomy 14 02661 g007aAgronomy 14 02661 g007bAgronomy 14 02661 g007c
Agronomy 14 02661 g008aAgronomy 14 02661 g008b
Figure 8. Bubble diagrams illustrating KEGG enrichment analysis of quinoa leaf samples from plants treated with different doses of P2O5. (A) W0 vs. W1; (B) W0 vs. W2; (C) W0 vs. W3; (D) W0 vs. W4; (E) W0 vs. W5.
Figure 8. Bubble diagrams illustrating KEGG enrichment analysis of quinoa leaf samples from plants treated with different doses of P2O5. (A) W0 vs. W1; (B) W0 vs. W2; (C) W0 vs. W3; (D) W0 vs. W4; (E) W0 vs. W5.
Agronomy 14 02661 g008aAgronomy 14 02661 g008b
Agronomy 14 02661 g009aAgronomy 14 02661 g009b
Figure 9. Histograms illustrating GO enrichment analysis of quinoa leaves treated with different doses of P2O5. (A) W0 vs. W1; (B) W0 vs. W2; (C) W0 vs. W3; (D) W0 vs. W4; (E) W0 vs. W5.
Figure 9. Histograms illustrating GO enrichment analysis of quinoa leaves treated with different doses of P2O5. (A) W0 vs. W1; (B) W0 vs. W2; (C) W0 vs. W3; (D) W0 vs. W4; (E) W0 vs. W5.
Agronomy 14 02661 g009aAgronomy 14 02661 g009b
Agronomy 14 02661 g010
Figure 10. Transcription factor family analysis of quinoa leaves treated with different doses of P2O5.
Figure 10. Transcription factor family analysis of quinoa leaves treated with different doses of P2O5.
Agronomy 14 02661 g010
Agronomy 14 02661 g011aAgronomy 14 02661 g011b
Figure 11. The validation of transcriptome data using qRT-PCR: the line graph indicates the transcriptome sequencing-detected transcript abundance FPKM expression with differential multiplicity. The bar graph shows the relative expression detected by qRT-PCR. Data are averages of three replicates and three biological replicates.
Figure 11. The validation of transcriptome data using qRT-PCR: the line graph indicates the transcriptome sequencing-detected transcript abundance FPKM expression with differential multiplicity. The bar graph shows the relative expression detected by qRT-PCR. Data are averages of three replicates and three biological replicates.
Agronomy 14 02661 g011aAgronomy 14 02661 g011b
Agronomy 14 02661 g012
Figure 12. Associate heat maps. The horizontal axis represents the name of the gene, and the vertical axis represents the name of the metabolite. Each square represents the correlation and significance between the gene and the metabolite. The magnitude of correlation is represented by different colors and shades, with red representing positive correlation and blue representing negative correlation. The darker the color, the higher the correlation, and the lighter the color, the lower the correlation. The significance p-value is represented by *; 0.01 < p < 0.05, *; 0.001 < p < 0.01, **; p < 0.001, ***. (A) CK vs. LP; (B) CK vs. HP.
Figure 12. Associate heat maps. The horizontal axis represents the name of the gene, and the vertical axis represents the name of the metabolite. Each square represents the correlation and significance between the gene and the metabolite. The magnitude of correlation is represented by different colors and shades, with red representing positive correlation and blue representing negative correlation. The darker the color, the higher the correlation, and the lighter the color, the lower the correlation. The significance p-value is represented by *; 0.01 < p < 0.05, *; 0.001 < p < 0.01, **; p < 0.001, ***. (A) CK vs. LP; (B) CK vs. HP.
Agronomy 14 02661 g012
Agronomy 14 02661 g013
Figure 13. Expression profiles of genes and metabolites related to flavonoid biosynthesis in quinoa under different phosphate fertilizer dosages. The boxes in the pathway represent DEGs or DAMs. Red and blue represent upregulated and downregulated genes/metabolites, respectively.
Figure 13. Expression profiles of genes and metabolites related to flavonoid biosynthesis in quinoa under different phosphate fertilizer dosages. The boxes in the pathway represent DEGs or DAMs. Red and blue represent upregulated and downregulated genes/metabolites, respectively.
Agronomy 14 02661 g013
Table 1. Sample Information.
Table 1. Sample Information.
GroupRepetitionPhosphorus Fertilizer Usage
W0W0-10 kg/hm2
W0-2
W0-3
W1W1-1112.5 kg/hm2
W1-2
W1-3
W2W2-1225.0 kg/hm2
W2-2
W2-3
W3W3-1337.5 kg/hm2
W3-2
W3-3
W4W4-1450.0 kg/hm2
W4-2
W4-3
W5W5-1562.5 kg/hm2
W5-2
W5-3
Table 2. Gene primer sequences in RT-qPCR.
Table 2. Gene primer sequences in RT-qPCR.
IDPrimerForward/Reverse PrimerPrimer Sequences (5′ to 3′)
1geneLOC110707704Forward primerGTCGTCATCATCATCATC
Reverse primerATATTCCAACCACTGTCT
2geneLOC110689796Forward primerGTTCGTGTTAATGTCTATG
Reverse primerTTCTATCCTCAAGTTCTTC
3geneLOC110702485Forward primerATTGAAGGCTATGAATCC
Reverse primerGTCTGAGGTTGATATGTC
4geneLOC110689796Forward primerGTTCGTGTTAATGTCTATG
Reverse primerTTCTATCCTCAAGTTCTTC
5geneLOC11000414Forward primerGTAAGGAGCACTATAACAA
Reverse primerGGACATATCAGAGACAAC
6geneLOC110735654Forward primerTTGCTAATCACTTCTCTG
Reverse primerGCTATCCACTTCACTATC
7TUB6Forward primerTGAGAACGCAGATGAGTGTATG
Reverse primerGAAACGAAGACAGCAAGTGACA
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Wang, H.; Li, H.; Li, X.; Wang, Q.; Liu, J.; Zhang, P.; Xie, H.; Li, L.; Qin, P. Integrative Transcriptomic and Metabolomic Analysis Reveals Quinoa Leaf Response Mechanisms to Different Phosphorus Concentrations During Filling Stage. Agronomy 2024, 14, 2661. https://doi.org/10.3390/agronomy14112661

AMA Style

Wang H, Li H, Li X, Wang Q, Liu J, Zhang P, Xie H, Li L, Qin P. Integrative Transcriptomic and Metabolomic Analysis Reveals Quinoa Leaf Response Mechanisms to Different Phosphorus Concentrations During Filling Stage. Agronomy. 2024; 14(11):2661. https://doi.org/10.3390/agronomy14112661

Chicago/Turabian Style

Wang, Hongxin, Hanxue Li, Xiaorong Li, Qianchao Wang, Junna Liu, Ping Zhang, Heng Xie, Li Li, and Peng Qin. 2024. "Integrative Transcriptomic and Metabolomic Analysis Reveals Quinoa Leaf Response Mechanisms to Different Phosphorus Concentrations During Filling Stage" Agronomy 14, no. 11: 2661. https://doi.org/10.3390/agronomy14112661

APA Style

Wang, H., Li, H., Li, X., Wang, Q., Liu, J., Zhang, P., Xie, H., Li, L., & Qin, P. (2024). Integrative Transcriptomic and Metabolomic Analysis Reveals Quinoa Leaf Response Mechanisms to Different Phosphorus Concentrations During Filling Stage. Agronomy, 14(11), 2661. https://doi.org/10.3390/agronomy14112661

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop