Sequence comparison of the GmPHO1 and GsPHO1 cDNAs
Fourteen putative PHO1 homologous genes were previously predicted [45], but 12 PHO1-like cDNAs were ultimately obtained from either the cultivated (prefix Gm) or wild (prefix Gs) soybeans (Additional file 1: Figures S1, S2). They were divided into Class Ι and Class II (Fig. 1, Additional file 1: Figure S2). When the Arabidopsis PHO1 genes were included, we found that the functionally characterized AtPHO1;H4 (SHB1) gene belonged to Class I, whereas AtPHO1;H1 and AtPHO1 were classified under Class II (Additional file 1: Figure S2), indicating functional divergence between Classes I and II. Moreover, each PHO1 gene from cultivated soybean Suinong 14 (SN14) and wild soybean ZYD00006 (ZYD6) were clustered together to form an orthologous pair (Fig. 1, Additional file 1: Figure S2), suggesting that these originated from a common progenitor. The paralogs underwent distinct diversification during evolution, with differences ranging from 26.9 to 95.2% (Additional file 2: Table S1), suggesting functional divergence, whereas the sequence identities of the 12 orthologous pairs ranged from 99.3 to 100%. Two members (PHO1;H9 and PHO1;H10) from Class I and two members (PHO1;H12 and PHO1;H14) from Class II were identical between SN14 and ZYD6, and the remaining pairs contained few amino acid substitutions and indels (Fig. 1, Additional file 2: Table S2). We then estimated the potential effects of these sequence variations between the orthologous pairs by PROVEAN and SNAP analyses, which indicated that all amino acid substitutions and indels were neutral (Additional file 2: Table S2), suggesting that the functions of the corresponding orthologous gene pairs might be not altered. Each PHO1 orthologous pair from the cultivated and wild soybeans may thus likely have similar functions.
Organ-specific expression of soybean PHO1-like genes
Assessment of expression diversity has been utilized to establish functional divergence among genes. To investigate the tissue-specific expression patterns of the PHO1-like genes, the total RNAs from the roots, leaves, stems, flowers, and developing fruits (1, 3, 5, and 15 days after pollination, shortened as fruit1, fruit3, fruit5, and fruit15) in SN14 and ZYD6 were subjected to qRT-PCR. Due to high sequence identity, common but specific primers were designed for paralogous gene pairs PHO1;H1/H4, PHO1;H9/H10, and PHO1;H12/H14. The results showed that the PHO1 genes belonging to Classes I and II exhibited different organ-specific expression patterns (Additional file 1: Figure S3). Most soybean PHO1 genes in Class Ι were differentially expressed among the vegetable and reproductive organs (Additional file 1: Figure S3a-e). A few extreme variations in expression were observed. PHO1;H7 was active in most developing organs but hardly detected in roots (Additional file 1: Figure S3d), whereas PHO1;H9/H10 were predominantly expressed in the roots (Additional file 1: Figure S3e). However, the expression patterns for genes in Class II were strikingly different from those in Class I. Soybean genes in Class II had similar expression patterns and were predominantly expressed in vegetable tissues, particularly the roots and stems (Additional file 1: Figure S3f-i). Moreover, the expression of these genes in fruits was relatively low (Additional file 1: Figure S3f-i).
Comparison of the orthologous Gm-GsPHO1 pairs indicated that Class I genes displayed diverse expression patterns in various organs between SN14 and ZYD6. GmPHO1;H7 had higher expression than GsPHO1;H7 in all assessed organs (Additional file 1: Fig. S3d), whereas GsPHO1;H3 had a relatively higher expression in flowers and early developing fruits (fruit1) than GmPHO1;H3 in SN14 (Additional file 1: Figure S3b). Nevertheless, some genes such as PHO1;H2 and PHO1;H6 showed complex expression patterns in different tissues between SN14 and ZYD6. GmPHO1;H2 exhibited lower expression levels in the roots, leaves, and flowers, but higher expression in other organs than GsPHO1;H2 (Additional file 1: Figure S3a). Similarly, PHO1;H6 showed lower expression levels in the roots and leaves but higher expression in the developing fruits in SN14 than in ZYD6 (Additional file 1: Figure S3c). Despite the complex expression patterns of genes in Class I, these showed relatively higher expression levels in developing fruits (fruit3-fruit15) in SN14 than in ZYD6. However, Class II genes were highly expressed in the vegetative plant parts in ZYD6 than in SN14 as well as in the roots, stem, or leaves (Additional file 1: Figure S3f-i).
The above observations between ZYD6 and SN14 were further verified by qRT-PCR analyses using additional four wild (Y1, Y2, Y3, and Y4) and four cultivated (Hefeng48, Nenfeng16, Heinong35, and Dongnong53) soybeans (Additional file 1: Figure S4). Despite few fluctuations, these reflected similar expression patterns in different tissues in ZYD6 and SN14. Taken together, the integrated analysis and main findings indicated that both Classes I and II exhibited conserved tissue-specific expression patterns either in wild or cultivated soybeans: Class I had a broader expression domains than Class II, in which most genes were predominantly expressed in the roots (Fig. 2). Interestingly, genes that had differential expression patterns in the roots between wild and cultivated soybeans had higher expression in the roots of wild soybeans than cultivated ones, whereas genes that were differentially expressed in fruit15 between G. max and G. soja had higher expression in fruit15 of cultivated soybeans than wild ones (Fig. 2), suggesting functional divergence of Classes I and II PHO1 genes in fruit (seed) development and stress responses.
Soybean PHO1-like gene expression under various stresses
Salinity stresses
To further illustrate the diversity of soybean PHO1-like genes in response to abiotic stresses, two-week-old seedlings of SN14 and ZYD6 were treated with different levels of salinity, and total RNA from the roots were subjected to qRT-PCR analysis (Fig. 3). The Class I genes PHO1;H2, PHO1;H7, and PHO1;H9/H10 and the Class II genes PHO1;H5, PHO1;H8, and PHO1;H12/H14 showed similar response patterns to salt stress between SN14 and ZYD6 (Fig. 3a, d-e g-i), whereas the remaining genes, including PHO1;H3 and PHO1;H6 in Class I and PHO1;H1/H4 from Class II, displayed relatively different expression patterns under salt stress between SN14 and ZYD6 (Fig. 3b-c, f). The Gm-GsPHO1;H12/H14 orthologous pairs were all consistently induced under all the concentration of salt stress with the stronger response in ZYD6. Although some orthologous pairs showed similar expression patterns under salt stress, one or two genes of the pairs only responded to certain concentrations such as PHO1;H5 and PHO1;H8 in Class II and the Class I members of PHO1;H2 and PHO1;H9/H10 (Fig. 3a, e, g-h). The expression of GmPHO1;H5 decreased only at 150 mM and 200 mM NaCl, but GsPHO1;H5 expression decrease at all concentrations. GmPHO1;H8 expression was induced by salt stress, whereas that of GsPHO1;H8 only increased with 200 mM and 250 mM NaCl treatment. Furthermore, paralogs of the Gs-GmPHO1 genes showed variations in both tendency and magnitude of gene expression. These findings indicate that the soybean PHO1 genes may have diverse roles in salt stress responses.
Low-Pi treatments
We also investigated the expression of the PHO1 genes in soybean roots and leaves under Pi-starvation (Fig. 4). Soybean PHO1 genes in Class I exhibited diverse expression profiles under normal Pi conditions (solid lines in Fig. 4a-e) and these showed complex responses to Pi-deficiency in the roots and leaves (dashed lines in Fig. 4a-e). However, the orthologous gene pairs of PHO1;H2, PHO1;H3, PHO1;H6, and PHO1;H7 showed similar responses to low-Pi treatments in both roots and leaves between ZYD6 and SN14, whereas PHO1;H9/10 depicted a similar expression pattern only in the leaves between ZYD6 and SN14. Similar PHO1;H9/10 expression patterns were observed in the roots of ZYD6 and SN14 in the first 2 weeks after low-Pi treatments; however, upregulation was observed in ZYD6 under Pi stress at the third week (Fig. 4e). In Class I, only PHO1;H9/10 was significantly and continuously upregulated in the roots of ZYD6 under Pi-starvation, suggesting its roles in response to low-Pi conditions.
Class II PHO1 genes (PHO1;H1/4, PHO1;H5, PHO1;H8, and PHO1;H12/14) showed similar expression profiles during root and leaf development under normal conditions (red and blue solid lines, Fig. 4f-i). With Pi-deficiency, these were similarly upregulated in the roots between ZYD6 and SN14 and reached a maximum at 3 weeks (red dashed lines in Fig. 4f-i). These genes were upregulated in the leaves of ZYD6, whereas their expression in the leaves of SN14 did not change (dashed blue lines in Fig. 4f-i).
These findings suggest that the PHO1 gene family may play roles in soybean adaptation and morphological divergence. Moreover, the expression variation is likely decisive evolutionary event between the orthologous pairs of PHO1 genes, which occurred during the divergence of wild and cultivated soybeans. We therefore functionally tested these assumptions by overexpressing the cDNA of these genes from cultivated soybean in transgenic Arabidopsis.
Overexpressing GmPHO1 genes in transgenic Arabidopsis
Due to the highly identical sequences within paralog pairs GmPHO1;H9/H10 and GmPHO1;H12/H14 (Additional file 2: Table S1), the GmPHO1;H10 and GmPHO1;H14 were excluded, and thus altogether 10 GmPHO1 genes from the SN14 (cultivar) were included in the transgenic Arabidopsis assays. Three T3 lines of each transgene that were confirmed via RT-PCR (Additional file 1: Figure S5) were further investigated, and all obtained GmPHO1 transgenic plants showed no obvious phenotypic variations compared to wild-type (WT) Arabidopsis, including flowering time, seed size, plant height, and germination rate under normal growth conditions (Additional file 2: Table S3), indicating that overexpressing GmPHO1 did not affect plant development.
GmPHO1 altered the salt tolerance of transgenic Arabidopsis plants
Because alterations in the expression of soybean PHO1 genes occur in response to salinity stresses, we first evaluated the behavior of the GmPHO1 transgenic plants under salinity stresses (Fig. 5, Additional file 1: Figure S6). The seed germination and seedling green rate were evaluated under salt stresses (125 and 175 mM NaCl, respectively). The germination rates of transgenic lines were all higher than those of WT under 125 mM NaCl for 2 days; however, all included plant lines could germinate normally with increasing treatment time (Additional file 2: Table S4). Similar patterns were observed with 175 mM NaCl treatment (Additional file 2: Table S4). These results indicate that transgenes could enhance seed germination speed.
No significant difference in the green rate of seedlings was observed under 125 mM NaCl conditions; however, differences were observed under 175 mM NaCl conditions (Additional file 1: Figure S6a). Only 35S::GmPHO1;H9 from Class I transgenic plants showed significantly higher seedling green rate than WT under NaCl stress (Fig. 5a, Additional file 1: Figure S6a). However, the 35S::GmPHO1;H1, 35S::GmPHO1;H5, 35S::GmPHO1;H8, and 35S::GmPHO1;H12 transgenic lines of Class II exhibited a higher seedling green rate than WT under 175 mM NaCl stress (Additional file 1: Figure S6a). These observations show that GmPHO1 enhances the salt tolerance, however, the function is obviously diverged between Class I and Class II since only one gene in Class I and many in Class II could confer salt tolerance to transgenic Arabidopsis plants.
Transgenic plants of 35S::GmPHO1;H6, 35S::GmPHO1;H8, and 35S::GmPHO1;H9 were selected for the evaluation of salt tolerance in the soil. The four-week old seedlings were supplied with 250 mM NaCl to induce salt stress. After 15 days of treatment, the WT and transgenic lines harboring GmPHO1;H6 showed extensive bleaching and wilting, whereas the 35S::GmPHO1;H8 and 35S::GmPHO1;H9 transgenic lines all exhibited intense green color and were more vigorous than the WT (Additional file 1: Figure S6b-e). Furthermore, 35S::GmPHO1;H8 and 35S::GmPHO1;H9 transgenic plants all showed significantly higher chlorophyll content than the WT and the 35S::GmPHO1;H6 transgenic lines (Additional file 1: Figure S6f). Moreover, the aboveground biomass of 35S::GmPHO1;H8 and 35S::GmPHO1;H9 transgenic plants after one-month treatment were much higher than the WT and 35S::GmPHO1;H6 lines (Additional file 1: Figure S6 g). These observations indicate that overexpressing soybean PHO1 genes (GmPHO1;H1, GmPHO1;H5, GmPHO1;H8, GmPHO1;H9, and GmPHO1;H12) improves salt tolerance in Arabidopsis.
GmPHO1 affects the salt tolerance pathways in transgenic Arabidopsis
The SOS2, SOS3, ADH, P5CS1, and FRY1 genes are essentially involved in salt tolerance pathways in plants [21,22,23,24,25,26]. To elucidate the possible connection of GmPHO1 with salt stress pathways, we checked the expression of these salt tolerance pathway genes in transgenic plants. The seven-day-old seedlings of WT, 35S::GmPHO1;H8, 35S::GmPHO1;H9, and 35S::GmPHO1;H6 transgenic Arabidopsis plants were transferred to a medium with 175 mM NaCl for 24 h, and total RNA from the whole seedlings was subjected to qRT-PCR. Under normal growth (0 mM NaCl) conditions, we found that each transgene could differentially affect the expression of these marker genes compared to the WT (highlighted by the black stars above the open color columns in Fig. 5b-d). Furthermore, we observed that ADH, P5CS1, SOS3, and FRY1 were upregulated by NaCl treatment, whereas SOS2 expression was not significantly affected in the WT background (highlighted by the black star above the blocked gray columns compared to the corresponding open gray columns in Fig. 5b-d). Moreover, these salt pathway marker genes were differentially expressed in the transgenic plants compared to the WT under salt stress conditions (highlighted by the red stars above the blocked color columns relative to the blocked gray column in Fig. 5b-d).
In the 35S::GmPHO1;H8 transgenic plants, except for the P5CS1 downregulation, the SOS2, SOS3, ADH, and FRY1 genes were upregulated under 0 mM NaCl. However, compared to the WT under salt treatments, only the expression of ADH sharply increased, whereas the other genes were downregulated (Fig. 5b). Nonetheless, relative to the corresponding transgenic lines under 0 mM NaCl, both ADH and P5CS1 were upregulated, and the other genes were downregulated under salinity stress (Fig. 5b). Compared to the WT, in 35S::GmPHO1;H9 transgenic Arabidopsis, SOS2, SOS3, ADH and P5CS1 were downregulated, whereas FRY1 was upregulated under 0 mM NaCl conditions (Fig. 5c). However, when treated with NaCl, the SOS3 and ADH genes were upregulated in these transgenic lines, and SOS2, FRY and P5CS1 were downregulated (Fig. 5c). Relative to the corresponding transgenic lines without NaCl supplementation, SOS3, ADH and P5CS1 were upregulated and the other genes were downregulated with NaCl treatments (Fig. 5c). Interestingly, in the 35S::GmPHO1;H6 transgenic lines, changes in the expression of the salt pathway marker genes were observed under both normal and salt stress conditions (Fig. 5d); however, no changes in salt tolerance were detected in the transgenic plants. Compared to the changes in gene expression patterns in the 35S::GmPHO1;H6 transgenic lines (Fig. 5d), the upregulation of the ADH gene in both 35S::GmPHO1;H9 and 35S::GmPHO1;H8 (Fig. 5b, c) might play an essential role in enhancing salt tolerance in the two transgenic analyses. Nonetheless, changes in the expression of the marker genes in the salt-tolerance pathways in overexpressing GmPHO1;H6 are suggestive of the inherent nature of PHO1 genes in salt response.
GmPHO1 is involved in the responses to Pi-starvation in transgenic Arabidopsis
We further investigated the responses of GmPHO1 transgenic plants under Pi-deficiency. The growth of primary and lateral roots, which are considered as diagnostic traits of plants in response to phosphate starvation, were shown in the presence of 35S::GmPHO1;H4 and 35S::GmPHO1;H6 (Fig. 6a-d). In response to Pi-deficiency, not only WT Arabidopsis showed a reduction of primary root length and an increase of lateral root number, but also all transgenic Arabidopsis lines did (Additional file 1: Figure S7). However, we found that the transgenic plants harboring GmPHO1;H4 and GmPHO1;H8 showed different responses to phosphate starvation compared to the WT, whereas the transgenic Arabidopsis lines harboring other GmPHO1 genes, including 35S::GmPHO1;H6 plants, did not exhibit alterations in the sensitivity to Pi-deficiency (Fig. 6a-d, Additional file 1: Figure S7). Overexpression GmPHO1;H4 in Arabidopsis was associated with longer primary roots and a higher number lateral roots than the WT (Fig. 6a, b, Additional file 1: Figure S7), thereby enhancing tolerance to Pi-starvation, whereas the GmPHO1;H8 transgenic lines displayed significantly shorter primary roots only under Pi-deficiency (Additional file 1: Figure S7), hence increasing the sensitivity of root to insufficient Pi.
GmPHO1 affects the Pi signaling pathway in transgenic Arabidopsis
IPS1, AT4, PHT1;4, and PHO1;H1 are upregulated in Arabidopsis during Pi-deficiency [4,5,6] and are considered to be marker genes in the Pi signal pathway. To further evaluate the possible roles of the GmPHO1 genes in the Pi pathway, we investigated the expression levels of these marker genes in transgenic Arabidopsis lines. The 35S::GmPHO1;H8, 35S::GmPHO1;H4, and 35S::GmPHO1;H6 plants were evaluated and compared to the WT. Under normal conditions, IPS1 and PHT1;4 were upregulated in the 35S::GmPHO1;H4 transgenic lines, whereas PHO;H1 was downregulated, and AT4 expression did not change (Fig. 6e-h). PHT1;4 expression did not change, whereas PHO;H1, IPS1, and AT4 were downregulated in the 35S::GmPHO1;H8 transgenic lines compared to the WT (Fig. 6i-l). The variations in the expression of these marker genes in the 35S::GmPHO1;H6 plants was similar to that in the 35S::GmPHO1;H8 plants (Fig. 6m-p). These results indicate that overexpression of GmPHO1 could affect the expression of genes that are involved in the Pi pathway.
However, under Pi-starvation conditions, we found that IPS1, AT4, and PHT1;4 were indeed upregulated by Pi-deficiency in WT Arabidopsis, whereas PHO1;H1 was significantly downregulated (comparisons between blue columns under Pi normal and deficiency conditions in Fig. 6e-p). Moreover, the expression levels of these marker genes in the 35S::GmPHO1;H4 transgenic lines decreased compared to the WT, indicating their insensitivity to Pi-deficiency (Fig. 6e-h), whereas these were all upregulated in the 35S::GmPHO1;H8 transgenic lines, showing hypersensitivity to Pi-deficiency (Fig. 6i-l). However, changes in the expression of these marker genes in the 35S::GmPHO1;H6 transgenic plants in response to Pi-starvation were similar to those of the WT (Fig. 6m-p). Therefore, in the GmPHO1 family, GmPHO1;H4 or GmPHO1;H8 might participate in the Pi pathway by affecting the expression of genes that are homologous to Arabidopsis PHO;H1, IPS1, AT4, and PHT1;4.
Responses of SN14 and ZYD6 to salinity stresses and Pi-starvation
Class II GmPHO1 genes apparently could affect the response of transgenic Arabidopsis plants to stresses. Whether do these genes contribute to the divergence of soybeans in response to stresses? The native expression of PHO1 genes in response to salinity and low-Pi stresses was diverse in the roots of SN14 and ZYD6 (Additional file 2: Table S5). However, if the expression variation pattern was compared, we found that most PHO1 genes, particularly Class II genes, responded similarly to salinity stresses, and half of these genes had sharp contrasts in response to Pi-starvation (Additional file 2: Table S6), hinting the differential response capacities to stresses between SN14 and ZYD6. We therefore evaluated the responses of SN14 and ZYD6 to salinity stress and Pi-starvation. Using various NaCl concentrations to evaluate soybean salinity tolerance [46, 47] in SN14 and ZYD6, soybean seedlings of both accessions showed similar phenotypic variations (Additional file 1: Figure S8), indicating that the two accessions have no significant differences in salinity tolerance. However, when we treated these soybeans in low-Pi conditions, differential responses, particularly in relation to root development, were observed between SN14 and ZYD6 (Fig. 7, Additional file 1: Figure S9). We found that the primary root length of ZYD6 was longer than that of SN14 during the three-week development after germination, whereas no distinct changes in response to Pi-deficiency were observed (Fig. 7a, b). However, the development of lateral roots in SN14 was apparently promoted by Pi-deficiency, but not in ZYD6 (Fig. 7a, c). The increase in the number of lateral roots in SN14 in response to Pi-starvation was accompanied by an increase both in total root length and biomass (Fig. 7a, c-f). However, Pi-deficiency induced a significant reduction in SN14 aboveground biomass (Fig. 7g, h), but not in ZYD6, indicating that wild soybean ZYD6 is insensitive to Pi-starvation.
Based on the results of transgenic Arabidopsis analysis, we then checked the expression of soybean genes that were homologous to the Arabidopsis Pi pathway marker genes. Searching in the NCBI and Phytozome databases did not identify any homologs of Arabidopsis IPS1 and AT4 in soybean, whereas soybean PHT1;4-like genes were found. Phylogenetic reconstruction of the PHT1;4 genes depicted the putative orthology of the soybean genes (Glyma.03G162800, Glyma.10G036800, and Glyma.19G164300) and Arabidopsis PHT1;4 (Additional file 1: Figure S10). The total RNAs from the roots were subjected to qRT-PCR analysis, which revealed that soybean PHT1;4-like genes showed similar expression profiles during root development in both ZYD6 and SN14, but had different expression levels in roots. In particular, Glyma.19G164300 showed trace levels of expression under normal conditions (Fig. 7i-k). Moreover, these were all upregulated in ZYD6 and SN14 under Pi-deficiency (Fig. 7i-k). Nevertheless, the increase in the level of gene expression in response to Pi-starvation in ZYD6 was significantly higher than that in SN14.
These results suggest the differential response capacities of wild and cultivated soybeans to salinity and Pi-starvation at both phenotypic and molecular levels, which was linked with the PHO1 gene family since the salt and Pi signaling pathways could be affected by the expression of some GmPHO1 genes.