1 Introduction

The symbiotic nitrogen fixation system utilized by legumes is the most efficient known biological nitrogen fixation system. In fact, it is responsible for generating the greatest amount of symbiotically fixed N2, totaling roughly 20 million tons each year (Herridge et al. 2008). Soybean root nodule nitrogen fixation is a symbiotic process, whereby Rhizobium colonize soybean roots to form nodular structures that carry out biological nitrogen fixation (Kijne et al. 1990; Yan et al. 2014; Nanjareddy et al. 2017; Dennis et al. 2015; Downie et al. 1990). In this process, root nodules capture N2 from the air and then the nitrogen is incorporated into nitrogen-containing compounds, which serve as an important nitrogen source. Thus, nitrogen fixation efficiency is a key determinant of soybean growth and yield. Indeed, root nodule-associated nitrogen fixation accounts for over 60% of total nitrogen accumulation in soybean plants during the various developmental stages (Downie et al. 1990; Zhao et al. 2017; Tate 1995).

Phosphorus (P) is the second most important element associated with restricted crop growth as 40% of global output of crops grown on arable land is limited by low soil phosphorus levels (Vance 2001; Richardson and Simpson 2011; Zhou et al. 2019; Ao, et al. 2014). Phosphorus plays an important role in soybean nitrogen fixation, since N2 fixation within root nodules demands abundant P for energy transfer and consumption, as well as for protein and lipid synthesis (Mullen et al. 1988; Gentili and Huss-Danell 2003). Consequently, inorganic phosphate (Pi) deficiency directly impairs many cellular functions, not only by reducing levels of available ATP necessary for enzymatic activity, membrane transport, utilization of inorganic N sources, and photosynthesis, but also by limiting biosynthesis of macromolecules (proteins, lipids, and nucleic acids) (Bosse and Kock 1998). By contrast, adequate phosphorus levels promote soybean morphogenesis and growth by supporting activities of glutamine synthetase, nitrogenase, and nitrate reductase, which are essential enzymes involved in amino acid and protein synthesis (Miao et al. 2006, 2007; Wu et al. 1999; Zarrin et al. 2006; Kutama et al. 2008). Proper application of P can enhance the capacity of nitrogen-fixing root nodules-promoting nitrogen by phosphorus (Li et al. 2021). Therefore, P is pivotal to nodule morphogenesis and efficient N2 fixation.

The process by which soybean plants sense soil P is important for symbiotic nitrogen fixation (Cassman et al. 1981; Miao and Liu 2010). The overall impact of P on nodule nitrogen fixation mainly stems from its effects on host plant growth, nodule formation, nodule growth, nodule function, and nodule metabolism (Robson and O’Hara 1981; Jakobsen 1985; Israel 1987, 1993; Sa and Israel 1991; Cassman et al. 1981; Martin et al. 2017; Tsvetkova and Georgiev 2003; Chen et al. 2011). The amount and rate of nodule nitrogen fixation can be influenced by either high or low soil P level and varies significantly among the various soybean growth stages (Yao et al. 2012, 2018; Wu 2012). Nevertheless, despite previous research efforts, the mechanism by which phosphorus affects soybean nodule-associated nitrogen fixation efficiency remains unknown.

P has dual functions as a critical mineral element in cells and as a signal when P is deficient in nodules (Gentili and Huss-Danell 2003). The PHT1 genes encode a family of plasma membrane phosphate transporters (Nussaume et al. 2011), GmPHT1;1 (Chen et al. 2019) and GmPHT1;4 (Qin et al. 2012) which are expressed in non-fixing regions of nodules, and play a role in nodulation and Pi homeostasis in nodules. High levels of GmPHT1;11 enhanced Pi accumulation in nodules, increased nodule size, and furthermore nitrogenase activity was enhanced by GmPHT1;11 (Lu et al. 2020). Stress-related metabolites such as polyols accumulate in P-deficient roots as well as sugars, which are known to be essential for P stress gene induction (Hernandez et al. 2007).

Recently, research advances have greatly enhanced our understanding of effects of phosphorus on soybean nitrogen fixation (Li et al. 2015). In this work, we analyzed effects of phosphorus levels on differential expression of metabolites and metabolic pathways during the nitrogen fixation process within root nodules. The results provided clues that reveal the key factors related to the physiological mechanism whereby phosphorus levels influence nitrogen fixation within soybean root nodules.

2 Materials and methods

2.1 Plant materials and sampling

Soybean (Suinong 14, SN14) plants were grown in sand medium in pots. The nutrient composition and concentration (mg/L) of the nutrient solution were as follows: (NH4)2SO4, 235.80; MgSO4, 240.00; CaCl2, 220.00; Na2MoO4·H2O, 0.03; CuSO4·5H2O, 0.08; ZnSO4·7H2O, 0.22; MnCl2·4H2O, 4.90; H3BO3, 2.86; FeSO4·7H2O, 5.57; and Na2EDTA, 7.45; P levels: P1 (P stress), KH2PO4, 4.39, K2SO4, 42.00, KCl, 36.00; P11 (P stress), KH2PO4, 48.26, K2SO4, 28.00, KCl, 24.00; P31 (normal P), KH2PO4, 136.00; P61 (high P), KH2PO4, 136.00, NaH2PO4·2H2O, 150.96. Nitrogen source is 15 N labeled ammonium sulfate with an abundance of 3.38%.

Before the vegetative cotyledon stage (VC, unfolded cotyledons), only 500 ml distilled water was supplied to plants once per day. From VC to V3 (third trifoliate leafs) stages, P31 nutrient solution was supplied. From V3 stage and thereafter, different P treatments were started for 15 days (R1 stage) supplied in 500 ml nutrient solution once a day.

Rhizobium inoculation was performed when opposite true leaves completely opened as follows: field soybean nodules from plants grown in Haerbin Heilongjiang province were collected during the previous year and stored in a refrigerator. They were ground and added to the nutrient solution, at a rate of approximately 5 g of nodules per liter. Inoculation of soybean plants was performed on five consecutive days to assure that each plant was well inoculated.

Samples (4 replicates for each treatment) were taken from 8:00 to 10:00 am, and the plants were divided into two parts above-ground and under-ground part at the cotyledon scars. After determining the fresh weight of the above-ground part, it was dried it at 65 °C and the dry weight measured. The nitrogen content and 15 N abundance of the dried sample enabled the calculation of nitrogen accumulation, nodule nitrogen accumulation and nodule nitrogen fixation rates. The underground part was removed completely, rinsed with distilled water and dried on filter paper to remove excess water, and then the nitrogenase activity of root nodules was measured. After the measurement, the nodules were removed, and dried at 65 °C for determining the dry weight.

Nodules were collected (6 replicate samples for each treatment) and frozen immediately in liquid nitrogen and stored at − 80 °C, and then used for metabolomics analysis.

2.2 Measurement and calculation methods

An acetylene reduction assay was used to measure the nitrogenase activity in root nodules as described by Gremaud and Harper (1989) and Chen et al. (2013). The washed soybean roots with nodules were placed in a 500 ml brown glass bottle and sealed with a rubber plug (with a hole). 50 ml air was extracted and 50 ml acetylene (99.9%) was added with a syringe, which reacted for 2 h. A GC7900 gas chromatograph analyzer (Model GC7900, Shanghai Techcomp Scientific Instrument Co., Ltd., China) was used to measure the ethylene content. Specific nitrogenase activity (SNA) was expressed as µmoles of ethylene formed per gram dry weight of nodules per hour. Acetylene reduction activity (ARA) was expressed as µmoles of ethylene formed per plant per hour.

For the determination of nitrogen content in plants, K2SO4 and CuSO4 were used as catalysts and concentrated sulfuric acid (H2SO4) for digestion, and the nitrogen content was measured by an automatic Kjeldahl nitrogen analyzer.

15N abundance determination: The 15N abundance was measured by a mass spectrometer (Thermo-Fisher Delta V Advantage IRMS, Waltham, MA, USA) in dual inlet (DI) mode.

$$\begin{array}{l}\mathrm{Total}\;\mathrm{Nitrogen}\;\mathrm{Accumulation}\;(\mathrm{TNA})=\mathrm{Dry}\;\mathrm{matter}\times\mathrm{Nitrogen}\;\mathrm{content}\\\mathrm{Fertilizer}\;\mathrm{Nitrogen}\;\mathrm{Accumulation}\;(\mathrm{FNA}\;)=\mathrm{TNA}\times({}^{15}\mathrm N\;\mathrm{abundance}\;\mathrm{of}\;\mathrm{sample})/{}^{15}\mathrm N\;\mathrm{abundance}\;\mathrm{of}\;\mathrm{fertilizer}\\\mathrm{Nodule}\;\mathrm{Nitrogen}\;\mathrm{Fixation}\;\mathrm{Accumulation}\;(\mathrm{NNFA})=\mathrm{TNA}-\mathrm{FNA}\\\mathrm{Ratio}\;\mathrm{of}\;\mathrm{Nodule}\;\mathrm{Nitrogen}\;\mathrm{Fixation}\;(\mathrm{RNNF})=(\mathrm{TNA}-\mathrm{FNA})/\mathrm{TNA}\end{array}$$

Statistical analyses were conducted using the SPSS 19.0 (SPSS Inc., Chicago IL) procedure. The multiple comparison was conducted by Duncan’s multiple range test at the 0.05 probability level.

2.3 Metabolites extraction

The freeze-dried samples were crushed with a mixer mill for 30 s at 45 Hz. Then 50 mg aliquot of each individual sample was precisely weighed and transferred to an Eppendorf tube. After addition of 1 ml of extract solution (methanol/water = 3:1, precooled at -40 °C, containing internal standard). After 30 s vortex, the samples were homogenized at 35 Hz for 4 min and sonicated for 15 min in ice-water bath. Then the samples were extracted over night at 4 °C on a shaker. Then centrifuged at 12,000 rpm for 15 min at 4 °C. The supernatant was carefully filtered through a 0.22 microporous membrane, then the resulting supernatants were diluted 10 times with methanol/ water mixture (v:v = 3:1, containing internal standard) and vortexed for 30 s and transferred to 2 mL glass vials. Then 30 ul was taken from each sample and pooled as QC sample. The sample was stored at -80 °C until the UHPLC-MS (Ultra high performance liquid chromatog-mass spectrometry) analysis.

2.4 UHPLC-MS analysis

The UHPLC separation was carried out using an EXIONLC System (SCIEX). The mobile phase A was 0.1% formic acid in water, and the mobile phase B was acetonitrile. The column temperature was set at 40 °C. The auto-sampler temperature was set at 4 °C and the injection volume was 2 µL.

A Sciex QTrap 6500+ (Sciex Technologies), was used for assay development. Typical ion source parameters were: IonSpray Voltage: +5500/-4500 V, Curtain Gas: 35 psi, Temperature: 400 °C, Ion Source Gas 1:60 psi, Ion Source Gas 2: 60 psi, DP: ± 100 V.

2.5 Data preprocessing and annotation

SCIEX Analyst Work Station Software (Version 1.6.3) was employed for MRM (Multiple reaction monitoring) data acquisition and processing (Hair 2010; Dunn et al. 2011). MS raw data (.wiff) files were converted to the TXT format using MSconventer. In-house R program and database were applied to peak detection and annotation.

3 Results

3.1 Effect of P levels on nitrogen accumulation, nodules and nitrogenase activity

With the increase of phosphorus levels, total nitrogen accumulation and nodule nitrogen fixation accumulation significantly increased, P stress significantly decreased the rate of nodule nitrogen fixation but the difference was not significant in P61 compared to P31 (Table 1). P stress significantly decreased the nodule number and weight but it was significantly increased in P61, according to single nodule weight. P stress significantly decreased nodule weight by reducing single nodule size, however the main reason of nodule weight increased in P61 was the increased nodule number (Table 2). The specific nitrogenase activity and acetylene reduction activity of P31 and P61 treatments was significantly higher than in the P stress treatments. The P61 treatment significantly increased acetylene reduction activity and had no significant impact on specific nitrogenase activity compared to P31 (Table 3).

The data are represented as mean values ± standard error (with four replicates). Values with the same letters are not significantly different at the 5% level according to Duncan’s multiple range test.

Table 1 Change in nitrogen accumulation vertical comparison
Full size table
Table 2 Change in the number and weight of root nodules vertical comparison
Full size table
Table 3 Change in specific nitrogenase activity (SNA) and acetylene reduction activity (ARA) of root nodules vertical comparison
Full size table

3.2 Nodule metabolic responses to P

There were significant changes in nitrogen accumulation, nodules and nitrogenase activity when different phosphorus levels were used as treatments for 15 days, as revealed by further metabolomics analysis.

Analysis of metabolic characteristics of nodules in plants grown at four concentrations of phosphorus revealed that there was a total of 1267 expressed metabolites. These included alkaloids, flavonoids, phenols, terpenoids, amino acids and amino acid derivatives, steroids and steroid derivatives, phenylpropanoids, fatty acyls, coumarins, organooxygen compounds, lignans, xanthones, nucleotides and nucleotide derivatives, phytohormone, and other compounds. The identified metabolites are listed in Table S1. The identified metabolites were analyzed using principal component analysis (Jolliffe 2002) to evaluate overall differences between P1, P11, P31, and P61 groups (Fig. 1). According to the first and the second principal component, the four P treatments were assigned to four different regions. Two principal components explained 18.4% of the total variance. The first principal component PC1 and the second principal component PC2 contributed to 11.3% and 7.1% of the interpretation rate, respectively. Thus, significant differences in variance were observed among the four P treatments.

Fig. 1
figure 1

Score scatter plot for PCA model TOTAL with QC

Full size image

A supervised orthogonal partial least squares discriminant analysis model was used to compare metabolic characteristics of P1, P11, P31, and P61. The OPLS-DA scatter scores revealed significant differences among pairwise comparison groups (Fig. 2) (Rygg and Wold 2002). Moreover, values of R2Y and Q2, along with results of 200 permutation tests confirmed that the quality of each supervised model was acceptable and that overfitting had not occurred (Table S2).

Fig. 2
figure 2

Score scatter plot of OPLS-DA model

Full size image

Analysis of results obtained using a volcano plot based visualization revealed differences in metabolites between different P treatments, with each point representing a metabolite and the point size representing a variable importance parameter (VIP) value based on the OPLS-DA model (Fig. 3). Ultimately, 200, 199, and 107 different metabolites were detected for P1 vs. P31, P11 vs. P31, and P61 vs. P31, comparisons, respectively. Regarding these metabolites, abundance of 110, 100, and 52 metabolites increased, while abundance of 90, 99, and 55 metabolites decreased, respectively (relative to P31 levels), with final values obtained for comparisons of VIP > 1.0 and p < 0.05 (Table S3). We found an interesting trend whereby the greater the difference in P levels between compared groups, the greater the numbers of differential metabolites were detected through pairwise comparisons.

Fig. 3
figure 3

Volcano plot. a P1 vs. P31; b P11 vs. P31; c P61 vs. P31

Full size image

Each point in the middle represents a metabolite, the abscissa represents the multiple change of the group compared to each substance (take the logarithm based on 2), and the ordinate represents the P-value of the Student’s t test (take the logarithm based on 10 Negative number), the size of the scatter point represents the VIP value of the OPLS-DA model, the larger the scatter point, the greater the VIP value. The scattered colors represent the final screening results. Metabolites that are significantly up-regulated are shown in red, metabolites that are significantly down-regulated are shown in blue, and metabolites that are not significantly different are in gray.

As a result of the large number of different metabolites screened among the various P treatments in this study, we applied the combination cutoff of VIP > 1.0 together with a fold change > 2 or < 0.5 as the threshold used for decision-making as to whether or not to conduct further analysis of a given differential metabolites. Ultimately, 273 and 168 differential metabolites were obtained from all groups (P1, P11, P31, P61) based on analysis of variance (ANOVA) significance cutoffs of < 0.05 and < 0.01, respectively. In Fig. 4, a heat map displays chemical classes and abundance of differential metabolites (with ANOVA significance values ≤ 0.01) within nodules of the four groups.

Fig. 4
figure 4

Heatmap of hierarchical clustering analysis

Full size image

Overall, 84 differential metabolites were detected in the P1 vs. P31 groups (50 with increased abundance and 34 with decreased abundance relative to P31), 81 differential metabolites in the P11 vs. P31 comparison group (43 with increased abundance and 38 with decreased abundance), and 20 differential metabolites in the P61 vs. P31 group (11 with increased abundance and 9 with decreased abundance relative to P31). These differential metabolites were then further analyzed for associations with metabolic pathways as determined using KEGG analysis.

3.3 Nodule metabolite changes associated with different P levels

Pairwise comparisons between P1, P11 or P61 with P31, revealed differences in differential metabolites based on chemical classification and regulation mode (Table S4). In P1 vs. P31, 84 differential metabolites were identified that belonged to 28 classes that mainly included categories for alkaloids (17.86%), amino acids and their derivatives (16.67%), flavonoids (10.71%), nucleotides and their derivatives (5.95%), and organooxygen compounds (5.95%). Regarding abundance of specific metabolites within each of these categories, 80% were increased, 78.57% were increased, 55.56% were decreased, 60% were decreased, and 100% were decreased, respectively, in P1 relative to P31. In P11 vs. P31, 81 differential metabolites were identified that belonged to 27 classes that mainly included alkaloids (20.99%), amino acids and their derivatives (12.35%), flavonoids (11.11%), and nucleotides and their derivatives (7.41%). Again, regarding abundance of metabolites within these categories, 76.47% were increased, 80% were increased, 55.56% were decreased, and 66.67% were decreased, respectively, in P11 relative to P31. In P61 vs. P31, 20 differential metabolites were identified that belonged to 12 classes that mainly included categories of alkaloids (15%), flavonoids (15%), organooxygen compounds (15%), amino acid and their derivatives (10%), and terpenoids (10%), With respect to abundance of metabolites in these categories, 33.33%, 66.67%, 66.67%, 50%, 50% were increased, respectively, in P61 relative to P31. Thus, these metabolites may be influenced by phosphorus-level effects on nitrogen fixation processes within soybean nodules.

Based on the abovementioned results, we further compared and analyzed the different metabolite compounds within each class and generated a heatmap of the changes of differential metabolites (VIP > 1.0 and fold change > 2 or < 0.5). The pairwise comparisons based on P treatment level are shown in Fig. S1. For each pairwise comparison group, corresponding ratios of quantitative values of corresponding differential metabolites were calculated, with the first 15 metabolites with most markedly increased or decreased abundance (relative to P31) is displayed using matchstick illustrations (Fig. S2). For each differential metabolite, values of correlation coefficients were calculated using the Pearson method and then were displayed in the heatmap (Fig. S3).

3.4 KEGG annotation and metabolic pathway analysis of differential metabolites

We identified 57, 59, and 34 KEGG pathways associated with differential metabolites detected for pairwise comparisons P1 vs. P31, P11 vs. P31, and P61 vs. P31, that included 413, 394, and 178 compounds, respectively (Table 4). KEGG metabolic pathway terms associated with highest numbers of differential metabolites included Metabolic pathways (gmx01100), Biosynthesis of secondary metabolites (gmx01110), Biosynthesis of amino acids (gmx01230), ABC transporters (gmx02010), Aminoacyl-tRNA biosynthesis (gmx00970), with numbers of differential metabolites differing among the three comparisons is provided (Table S5). Therefore, nitrogen fixation in soybean root nodules may respond to phosphorus stress through activities associated with these metabolic pathways. Phosphate and amino acid transporters, phosphate, glutamate, glutamine and aspartate were decreased in P1 vs. P31.

Table 4 Pairwise comparison of KEGG pathways and compound differences
Full size table

After KEGG annotation analysis revealed all the pathways associated with differential metabolites, we performed enrichment analysis and topological analysis on these pathways to reveal key pathways associated with different P levels. Ultimately, the results revealed detection of 50, 55, 27 pathways for P1 vs. P31, P11 vs. P31 and P61 vs. P31, respectively (Fig. 5). The numbers of compounds, differential metabolites and chemical classes differed among the pairwise comparisons (Table 5). Aminoacyl-tRNA biosynthesis, beta-alanine metabolism, alanine, aspartate and glutamate metabolism, galactose metabolism, glycine, serine and threonine metabolism, isoquinoline alkaloid biosynthesis, indole alkaloid biosynthesis, arginine and proline metabolism, ascorbate and aldarate metabolism, and glutathione metabolism were the main pathways identified for P1 vs. P31. Flavone and flavonol biosynthesis, nicotinate and nicotinamide metabolism, glycine, serine and threonine metabolism, phenylalanine metabolism, beta-alanine metabolism, isoquinoline alkaloid biosynthesis, amino sugar and nucleotide sugar metabolism, and glutathione metabolism were the main pathways identified for P11 vs. P31. Alanine, aspartate and glutamate metabolism, beta-alanine metabolism, flavone and flavonol biosynthesis, galactose metabolism, glutathione metabolism, isoquinoline alkaloid biosynthesis, pantothenate and CoA biosynthesis, and tyrosine metabolism were the main pathways identified for P61 vs. P31. Ultimately, these metabolic pathways may provide important clues for use in uncovering mechanisms by which phosphorus levels affect nodule nitrogen fixation. Notably, compounds associated with these metabolic pathways for P1 vs. P31 belonged to 6 chemical classes (amino acids and their derivatives, nucleotide and their derivatives, vitamins, pyrimidine nucleotides, carbohydrates and alkaloids). For P11 vs. P31 the pathways belonged to 7 chemical classes (amino acids and their derivatives, nucleotides and their derivatives, pyrimidine nucleotides, carbohydrates, organooxygen compounds, flavonoids, and benzene and substituted derivatives), and for P61 vs. P31 they belonged to 7 chemical classes (amino acid and derivatives, alkaloids, flavonoids, organooxygen compounds, hydroxycinnamoyl derivatives, phenylpropanoids and others). For each pathway identified for each P level-based pairwise comparison, associated metabolite compounds and their relative abundance changes (increased or decreased relative to P31) were different and shown in (Table S6).

Fig. 5
figure 5

Bubble plot display pathway analysis results. a P1 vs. P31; b P11 vs. P31; c P61 vs. P31

Full size image

Each bubble in the bubble diagram represents a metabolic pathway. The abscissa of the bubble and the size of the bubble indicate the influence factor of the pathway in the topology analysis. The larger the size, the greater the influence factor; the ordinate and color of the bubble indicates the enrichment analysis. P value (take the negative natural logarithm, namely -ln(p)), the darker the color, the smaller the P value, the more significant is the degree of enrichment.

Table 5 Number of compounds, differential metabolites and classifications
Full size table

3.5 Analysis of regulatory network of differential metabolites

A metabolic network is a network composed of metabolic reactions and regulatory mechanisms that regulates reactions associated with cellular metabolic and physiological processes. After obtaining corresponding differential metabolite information for each P level-based pairwise comparison group, we performed pathway searches and regulatory interaction network analysis using the soybean KEGG database. For the P1 vs. P31 comparison group, we identified three pathway matches that included alanine, aspartate and glutamate metabolism (gmx00250), glycine, serine and threonine metabolism (gmx00260), and beta-alanine metabolism (gmx00410), as well as 2 modules of serine biosynthesis, glycerate-3P = > serine (M00020) and polyamine biosynthesis, arginine = > ornithine (M00134), 27 enzymes, 158 reactions, and 110 compounds. For the P11 vs. P31 comparison group, we identified one pathway match, alanine, aspartate and glutamate metabolism (gmx00250), as well as 2 modules that included GABA (gamma-aminobutyrate) shunt (M00027) and polyamine biosynthesis, arginine = > ornithine (M00134), 24 enzymes, 163 reactions, and 110 compounds. For the P61 vs. P31 comparison group, there were 3 pathway matches that included galactose metabolism (gmx00052), tropane, piperidine and pyridine alkaloid biosynthesis (gmx00960), and betalain biosynthesis (gmx00965), as well as 1 module of biotin biosynthesis, pimeloyl-ACP/CoA = > biot (M00123), 48 enzymes, 187 reactions, and 61 compounds (Fig. S4 and Table S7). These networks may play important roles in processes associated with phosphorus-based regulation of nodule nitrogen fixation.

4 Discussion

Metabolomics is the study of metabolic phenotypes and their relationships with corresponding genes under specific conditions. When used to determine endpoints of biological systems regulation, metabolomics can provide important information associated with changes in organisms. The changes include biological phenotypes, which are observable and measurable changes (e.g., metabolite characteristics) that can be used to understand biological processes and the mechanisms, both intuitively and effectively. For this reason, the field of metabolomics has developed rapidly and has been incorporated into many research fields, including those focused on understanding molecular mechanisms that determine plant phenotypes (Liu et al. 2020; Wei et al. 2019; Tang et al. 2014; Dunn et al. 2011; Wiklund et al. 2008). Many studies have shown that phosphorus plays specific roles in nodule initiation, growth, and function, in addition to its involvement in host plant growth processes (Yao et al. 2018; Miao et al. 2007; Israel 1987, 1993; Gentili and Huss-Danell 2003; Sulieman et al. 2013; Thuynsma et al. 2014). In this study, we analyzed the effects of phosphorus levels on metabolite types and metabolite-associated pathways within soybean nodules.

4.1 Effects of phosphorus level on soybean metabolites and KEGG analysis results

Stress-related primary and secondary metabolites, such as polyols and sugars, are known to accumulate in root nodules of P-deficient bean plants (Hernandez et al. 2007). Subsequently, 10% of these compounds are released by root nodules into the surrounding environment (Tawaraya et al. 2014). Primary metabolites in roots of plants grown under P-deficiency have been shown previously to involve amino acids, organic acids, polyhydroxy acids, fatty acids, sugars, sugar phosphates, polyols, and other nitrogenous compounds (Hernandez et al. 2007). In the present study, we detected abundant metabolites within nodules of plants grown under low phosphorus (P1, P11) and high phosphorus (P61) conditions; and then identified changes in metabolite types and abundance in P1, P11 or P61 groups relative to corresponding P31 indicators (P1 vs. P31, P11 vs. P31, P61 vs. P31) (Table S1).

Our results demonstrated that relative to P31, greater numbers of differential metabolites were found in P1 and P11 than in P61 (Table S3). These metabolites mainly included the amino acids (L-tyrosine, L-arginine, L-asparagine, L-ornithine, L-citrulline, L-histidine) that were present in increased abundance in P1 and P11. This result conflicted with results of previous studies that demonstrated reduced biosynthetic pathway activities associated with amino acid synthesis (gmx01230) in P-deficient nodules (Hernandez, et al. 2009). The possible reasons are that the material used in the present study was soybean and not the common bean, and in addition the concentration, time and treatment methods of phosphorus stress treatment werte quite different. However no significant differences, associated with amino acid rmetabolite quantity or expression, were observed between P61 and P31. Although P stress inhibited nitrogen fixation ability of root nodules by reducing number, weight, specific nitrogenase activity and acetylene reduction activity, the root nodules responded to P stress by increasing the number of amino acids and their derivates. These differential amino acids were mainly up-regulated and involved ABC transporters.

Abscisic acid and phosphorylcholine were decreased under P-deficient conditions (P1 and P11) but were increased in P61. ABA is a hormone that is involved in plant defense mechanisms against diverse stresses and it initiates the stress response in plants through transduction of stress signals. It also improves the ability of plants to resist effects of various types of stress factors (Chen et al. 2006; Gawronska et al. 2003). The reduction of ABA under P deficiency is a very effective mechanism for plant survival under P deficiency stress (Zhu et al. 2018). In this study, abscisic acid (ABA) was mainly associated with metabolic pathways (gmx01100), biosynthesis of secondary metabolites (gmx01110), and plant hormone signal transduction (gmx04075). ABA also inhibited carotenoid biosynthesis as an adaptation to P stress in root nodules. Another identified metabolite, phosphorylcholine, has been shown to be important in phospholipid metabolism. In rice it has been shown to be involved in controlling phosphorus balance and adaptive responses to phosphorus deficiency (Chen et al. 2006). It has also been shown to play a highly significant role in phosphorus deficiency-associated cellular phospholipid remodeling processes and adaptations in rice and other crops. Phosphorylcholine was down-regulated which affected glycerophospholipid metabolism in response to P stress in our study.

Abundance of D-glucose 6-phosphate and galactose 1-phosphate decreased under P1 and P11, but increased under P61 treatment compared to P31. D-glucose 6-phosphate is a key compound associated with sugar metabolism and is involved in glycolysis pathway, pentose phosphate pathway, starch and sucrose metabolism, etc. Phosphorus is associated with energy storage and utilization pathways, and our results suggest that D-glucose 6-phosphate was mainly associated with metabolic pathways (gmx01100), biosynthesis of secondary metabolites-other antibiotics (gmx00998), and starch and sucrose metabolism (gmx00500). ATP is necessary for nitrogen fixation in soybean nodules. The decrease in the abundance of D-glucose 6-phosphate affects energy metabolism, and ATP decreases, thereby inhibiting nitrogen fixation in soybean nodules.

Abundance of uridine 5’-diphospho-D-glucose, adenosine monophosphate and phosphoric acid decreased under P1 and P11. Purine metabolism, down-regulation affects the synthesis of nitrogen metabolism intermediates allantoate and allantoin, which are the main forms of nitrogen transport and storage in soybean and this results in a reduction in nitrogen fixation products (Zhang et al. 1988).

In the current study, changes in ABA, phosphorylcholine, D-glucose 6-phosphate and Pi content mainly involved carotenoid biosynthesis, glycerophospholipid metabolism and sugar metabolism and were associated with the effect of P stress on nitrogen fixation in soybean root nodules. This inhibited the process of nitrogen fixation and the synthesis of nitrogen fixation products.

4.2 Flavonoid differences in nodules

With regard to the symbiosis between soybean roots and rhizobia, the conduction of biochemical signals is a prerequisite for establishing a symbiotic system (Cooper 2004; Miao and Liu 2010). In the symbiotic process between rhizobia and legumes, flavonoids act as initial signaling molecules for nodulation induction by acting to chemoattract rhizobia, regulate synthesis of rhizobia nodulation factor, increase transcriptional expression of nodulation genes, and promote nodulation of rhizobia to enable nitrogen fixation (Chou and Wei 2010; Cooper 2007; Redmond and Batley 1986; Geurts et al. 2005; Walker et al. 2003). Flavonoids such as citrinin, saccharol, apigenin, quercetin, and others have been found in root secretions of legumes, such as soybean (Glycine max), alfalfa (Medicago Sativa), peas (Pisum sativum L.) and pigeon peas (Cajanus cajan). In this study we discovered differences in types and amounts of isoflavones in soybean nodules in plants grown under different phosphorus levels.

Numerous papers dealing with the effects of flavonoids on nodulation can be found in the literature. Several studies have demonstrated that digitoxin obtained from alfalfa can effectively induce expression of the nod gene of Rhizobium meliloti (Novak et al. 2002). In rapeseed, flavonoids have been shown to effectively increase the colonization rate of the root system by Azorhizobium caulinodans ORS571 (Hu et al. 1999). In soybean, the simultaneous addition of isoflavones daidzein and genistein has been shown to effectively stimulate nodulation (An et al. 2017). In bean plants, additions of naringenin, hesperetin, quercetin, and genistein all significantly affected faba bean rhizobial colony diameter, number, and absorbance (Wang et al. 2018). Another study demonstrated that when naringenin, hesperetin, and quercetin were added to a pea-based culture system in high concentrations (10–100 µg/mL), the flavonoids promoted curling and deformation of pea root hairs (Novak et al. 2002).

In this study, we obtained some new ideas on the basis of results that were consistent with previous reports. For example, 25, 21, 9 flavonoids were detected in P1 vs. P31, P11 vs. P31, and P61 vs. P31 groups, in which the abundance of 14, 9, 4 flavonoids (relative to P31 levels) were increased, respectively (Table S8). This indicated that, under P stress, more flavonoids were involved in metabolic processes in soybean root nodules, including nodule nitrogen fixation. Specifically, quercetin abundance was increased in nodules of all P-level based treatment groups. This leads to the speculation that under P stress, the decrease in the number of nodules was smaller than the decrease in the weight of a single nodule. Thus the role of quercetin in nitrogen fixation was mainly to promote the formation of root nodules. Hesperetin, dihydromyricetin, 3’,4’,7-trihydroxyisoflavone, linarin and homopterocarpin were increased, whereas daidzin, glycitin, irisflorentin, pelargonidin-3,5-O-diglucoside chloride, and hesperetin 7-O-glucoside were decreased in soybean root nodules to regulate the effect of P stress on nodule nitrogen fixation.

4.3 Metabolic pathways and networks in response to phosphorus level changes

Great progress has been made in elucidating the steps involved in nodule formation. The first step in this process involves a molecular dialogue between leguminous plants and rhizobial bacteria that is established through recognition of bacterial Nod factors by legume Nod factor receptors, such as NFR1 and NFR5 (Limpens and Bisseling 2003; Oldroyd and Downie 2004; Popp and Ott 2011). For example, this enables Lotus japonicus to recognize its bacterial microsymbiont Mesorhizobium loti, that is required for the initial physiological and cellular responses to the lipochitin-oligosaccharide signal (Radutoiu et al. 2003). Plant responses to P deficiency, as reflected by metabolites within root exudates, differ among metabolites and among growth periods (Tawaraya et al. 2014). Nevertheless, effects of phosphorus on metabolic pathways and metabolic networks within nodules still needs to be further clarified.

After completion of KEGG annotation-based identification of the pathways involved in differential production of metabolites, metabolic transmission analyses (enrichment analysis and mutation analysis) were identified and the various metabolites studied find key pathways with highest correlations with metabolite differences. Further screening, and network analysis of metabolic pathway regulation and interaction revealed that the numbers of differential metabolites detected for P1 vs. P31 and P11 vs. P31 were greater than numbers detected for P61 vs. P31. This indicated that metabolic pathway responses differed among comparison groups due to varying degrees of phosphorus stress (Fig. 5). Specifically, a total of 1155, 1292, and 788 metabolic pathway-associated compounds were detected for comparison groups P1 vs. P31, P11 vs. P31, and P61 vs. P31, respectively, with 125, 123, 48 of the compounds shown to be enriched relative to P31 in P1, P11, and P61, respectively (Table S6). For P1 vs. P31, 15 metabolites were associated with this pathway, of which 3 were differential metabolites. For P11 vs. P31, 15 metabolites were associated with the pathway, including 2 different metabolites. By contrast, no metabolites of this metabolic pathway were detected for P61 vs. P31. Several metabolic pathways were identified through KEGG analysis, such as glutathione metabolism, isoquinoline alkaloid biosynthesis, and aminoacyl-tRNA biosynthesis and these played key roles in plants under each treatment (May et al. 1998; Vitale et al. 1993; Chevalier and Rossignol 2011). However, metabolites involved in metabolic pathways and their regulatory effects differed among groups based on P levels. The pathway ascorbate and aldarate metabolism was closely related to P stress, which implies that this pathway and associated metabolites are involved in the adaptation of the symbiotic soybean root nodule system to P starvation.

5 Conclusion

Under low phosphorus stress conditions (P1, P11), the numbers of differential metabolites and expression abundance of associated metabolic pathways were significantly increased relative to corresponding to the P61 treatment. The metabolites with increased abundance mainly involved amino acids and various flavonoids, while abscisic acid, phosphorylcholine, D-glucose 6-phosphate, uridine 5’-diphospho-D-glucose, and Pi were decreased. The differential metabolites were mainly associated with carotenoid biosynthesis, glycerophospholipid metabolism, sugar metabolism, ascorbic acid and aldonic acid metabolism and these regulate nodule nitrogen fixation under P stress. The results of the present study will be used to analyze further the physiological mechanisms of nodulation and nitrogen fixation in soybean under phosphorus stress.