3.1 Appearance characteristics of Lotus japonicus under low-phosphorus and drought stress
Lotus japonicus was grown for one month under low-phosphorus and drought stresses. Compared with the control, the growth of the aboveground parts of the pulse roots became weaker, the stem became shorter and the number of leaves decreased (Fig. 1a, 1b). The growth potential of the underground root system of Baimai roots became weaker, the main roots became shorter, and the number of lateral roots decreased significantly (Fig. 1C, 1D).
3.2 Overview of transcriptome sequencing
Through the transcriptome analysis of 9 Lotus japonicus samples subjected to the control, low-phosphorus stress and drought stress treatments, a total of 57.09 GB clean data were obtained, the amount of clean data from each sample reached 6.05 GB, and the q30 base percentage was 94.67% or more (Supplementary Table 1). In addition, we mapped the sequence base content and sequencing quality values to determine the sequencing quality (Supplementary Fig. 1). The clean reads of each sample were sequenced to obtain the genomes of pulse roots, and the alignment efficiency ranged from 82.44–84.02% (Supplementary Table 2). Moreover, the coverage depth distribution map of the reads of each sample mapped against the reference genome (Supplementary Fig. 2) and the distribution map of the reads of each sample mapped against different regions of the genome (Supplementary Fig. 3) were generated. The degree of random mRNA fragmentation was tested by simulating the results of mRNA fragmentation by assessing the distribution of mapped reads in each mRNA transcript (Supplementary Fig. 4). To test the effect of magnetic bead purification during library preparation, an insert fragment length test was performed (Supplementary Fig. 5). The number of genes sequenced from each sample was determined based on saturation, and the performance trend was the same (Supplementary Fig. 6). The gene expression distribution intervals in the transcriptomes of the 9 pulse root samples were highly conserved (Fig. 2A). The assessment of correlations (squared Pearson correlation coefficient) of expression between different samples showed that the correlations among the control, drought and low-phosphorus treatments were strong (Fig. 2B). Principal component analysis (PCA) showed that the gene expression patterns differed greatly under the control, drought and low-phosphorus treatments, indicating that the sample treatment results were reliable (Fig. 2C).
3.3 Differentially expressed genes and enrichment analysis
Differentially expressed genes were determined by pairwise analyses between the control, low-phosphorus, and drought treatments using the DESeq2 tool (Supplement: Fig. 7). A total of 2176 differentially expressed genes were screened in zl1 vs. zl2(zl12), including 994 upregulated genes and 1182 downregulated genes (Supplement: Table 3). The screening of zl1 vs. zl3 (zl13) revealed 3026 differentially expressed genes, including 1990 upregulated genes and 1036 downregulated genes (Supplement: Table 4). A total of 2980 differentially expressed genes were screened in zl2 vs. zl3(zl23), among which 2115 were upregulated, and 865 were downregulated (Supplement: Table 5). In addition, Venn analysis of the three groups of differentially expressed genes was performed (Fig. 3A), and the numbers of common differentially expressed genes among different groups were recorded. zl12 and zl13 shared 700 differentially expressed genes, zl12 and zl23 shared 593 differentially expressed genes, zl13 and zl23 shared 1121 differentially expressed genes, and zl12, zl13 and zl23 shared 168 differentially expressed genes (Fig. 3A).
To explore the function of differentially expressed genes between different treatments, we annotated the GO functions of differentially expressed genes between different treatments and displayed the top 20 GO functions with significant enrichment (Fig. 3B). According to the differential expression gene enrichment functions of control zl1 vs. zl2(zl12), control zl1 vs. zl3(zl13) and control zl2 vs. zl3 (zl23), the genes identified in the three comparisons were significantly enriched in cell, membrane and plant type functions in the cellular component category. In addition, three kinds of photosystem functions were enriched under low-phosphorus–drought. In terms of molecular functions, the genes identified among the three comparison groups were significantly enriched in various enzyme activities and ion binding functions. In addition, the low-phosphorus control was enriched in transcription factor functions. Among biological processes, the genes of the three groups were mainly concentrated in various metabolic and transport-related functions. In addition, the GO enrichment results of the common differentially expressed genes of zl12, zl13 and zl23 were all in accord with the above results (Supplement, Fig. 8).
To understand the KEGG pathway enrichment of differentially expressed genes, we focused on the KEGG pathway annotation of differentially expressed genes among the low-phosphorus control, drought control and low-phosphorus–drought treatments and identified the top 20 pathways with significant enrichment (Fig. 3C). The differentially expressed genes in zl12 were enriched in the plant hormone signal transduction, oxidative phosphorylation, peroxisome, and phosphate management pathways. The differentially expressed genes in zl13 were enriched in plant hormone signal transduction, MAPK signaling pathway-plant and flavonoid biosynthesis, which are pathways that regulate drought stress in plants. The differentially expressed genes in zl23 were enriched in the phosphatidylinositol signaling system and flavonoid biosynthesis, which regulate phosphoric acid and drought stress pathways. In addition, the differentially expressed genes common of zl12, zl13 and zl23 were enriched in MAPK signaling pathway-plant and plant hormone signal transduction (Supplement, Fig. 8).
3.4 Osmotic stress
The osmotic stress response is the core of the plant response to drought stress. Therefore, we mapped the differentially expressed genes identified in drought-related controls to osmotic stress regulatory networks (Supplement, Table 6). The results showed a total of 6 differentially expressed genes in the ABA signaling pathway, 2 of which were related to ABA receptor PYR/PYLs, and their expression levels were increased. One differentially expressed gene corresponded to RALF, and two differentially expressed genes corresponded to fer. One differentially expressed gene corresponded to GEF, and its expression was reduced. A total of 9 differentially expressed genes were included in the osmotic stress signaling pathway, 6 of which encoded MAPKs, and their expression was increased. The expression of three differentially expressed genes corresponding to RBOHB/C was increased. In addition, 72 differentially expressed genes encoded transcription factors. Most of the 22 WRKY genes were highly expressed under drought stress, and some of these genes were also highly expressed under low-phosphorus stress. The expression of 7 of 10 bhlhs was increased. The expression of 10 of the 11 NACs was increased. The expression of 18 of the 24 mybs was increased. The expression of 2 of the 5 AP2s was increased. It is worth noting that POD and HKT, related to nonenzymatic antioxidants, were also ranked as ninth and first, respectively. The expression of POD increased and that of HKT decreased (Fig. 4).
Note ABA, abscisic acid; PYR/PYLs, pyrabactin resistance 1-like protein; PP2C, type 2C protein phosphatases; SnRKs, Snf1 (sucrose nonfermenting-1)-related protein kinases; RBOHD/F, respiratory burst oxidase homolog D; SLAC1, slow anion channel 1; CPK/CBL, calcium-regulated phosphorylation systems; CAT, catalase; RALF, rapid alkalinization factor; FER, Feronia; GEF, guanine nucleotide exchange factor; MAPK, mitogen-activated protein kinase; ABFs, ABF, transcription factor; WRKY, WRKY transcription factor; bHlH, bHlH transcription factor; NAC, NAC transcription factor; MYB, MYB transcription factor; AP2, AP2/ERF transcription factor, GRAS, GRAS transcription factor; POD, peroxidase; GST, glutathione S-transferase zeta; HKTs, cation transporters.
3.5 Plant hormone signal transduction pathway response gene analysis
We selected the FPKM values of differentially expressed genes enriched in the plant hormone signal transduction pathway from the zl1 vs. zl2(zl12), zl1 vs. zl3(zl13) and zl2 vs. zl3(zl23) comparisons for analysis. The results showed that among the zl12 differentially expressed genes, the auxin pathway was enriched in two gene types, and two GH3 genes and one SAUR gene were upregulated. The gibberellin pathway was enriched in 4 TF genes, and 3 of these genes were upregulated. Abscisic acid was enriched in two PYR/PYL genes, and their expression was upregulated. Only one ERF1/2 gene was enriched in the ethylene pathway, and its expression was downregulated. Nine JAZ genes were enriched in the jasmonic acid pathway, and their expression was upregulated. The salicylic acid pathway was enriched in 2 NPR1 genes, 3 TGA genes and 3 PR-1 genes, all of which were downregulated (Fig. 5A). Among the zl13 differentially expressed genes, the auxin pathway was enriched in one TIR1 gene, whose expression was downregulated; one AUX/IAA gene, which was upregulated; six GH3 genes, 4 of which were upregulated; and three SAUR genes, one of which was upregulated. One AHP and one B-ARR gene were upregulated in the cytokinin pathway. The gibberellin pathway was enriched in two TF genes, one of which was upregulated. The abscisic acid pathway was enriched in two PYR/PYL genes, whose expression was upregulated. In the ethylene pathway, two EBF1/2 and two ERF1/2 genes were enriched and upregulated. The jasmonic acid pathway was enriched in 10 JAZ genes, and their expression was upregulated. Two TGA genes were annotated in the salicylic acid pathway, and their expression was downregulated. Two PR-1 genes were also identified, one of which was upregulated (Fig. 5B). Among the zl23 differentially expressed genes, three AUX/IAA genes were enriched in the auxin pathway, and one of the AUX/IAA genes was upregulated. Five GH3 genes were identified, four of which were upregulated. Three AHP genes were enriched in the cytokinin pathway, and two were upregulated. The gibberellin pathway was enriched in one TF gene, which was downregulated. The ABF gene was enriched in the abscisic acid pathway and was upregulated. The ethylene pathway enriched in two ERF1/2 genes, and their expression was upregulated. The jasmonic acid pathway was enriched in 6 JAZ genes, and their expression was upregulated. The salicylic acid pathway was annotated with 2 PR-1 genes, both of which were upregulated (Fig. 5C).
3.6 Transcription factors
We selected bZIP (67), PP2C (91), LEA (70), PLATZ (14) and ACX (3) transcription factors related to growth and development and stress and selected SPS (8), ALA (20) and PIP5K (17) gene families related to the phosphorus response. The measurement of gene expression differences showed that the bZIP and PLATZ families presented little difference in the number of members with significant expression under control, low-phosphorus and drought conditions. Most PP2C and LEA genes were significantly expressed under drought conditions. Two ACX2 genes were significantly expressed under drought treatment. Most SPS genes were significantly expressed under low-phosphorus conditions. Most ALA genes were significantly expressed under drought conditions. Most PIP5K genes were significantly expressed under control conditions (Fig. 6).
3.7 qRT - PCR testing of differentially expressed genes
To verify the reliability of the gene expression data obtained by transcriptome sequencing, 10 genes were randomly selected from 168 differentially expressed genes in the control group, low-phosphorus group and drought group for qRT-PCR verification. The qRT-PCR expression levels and expression trends of the differentially expressed genes were consistent with the transcriptome sequencing results, indicating high reliability of the transcriptional data (Fig. 7).
3.8 Differential metabolite analysis
We analyzed 730 metabolites from 18 samples with 6 replicates under control, low-phosphorus and drought stresses(Fig. 8A). There were 390 different metabolites identified in the control group, among which 76 were upregulated, and 314 were downregulated. There were 397 differential metabolites identified in the drought control, among which 130 were upregulated, and 267 were downregulated. There were 332 different metabolites identified under low-phosphorus–drought stress, among which 209 were upregulated, and 123 were downregulated. A total of 620 differential metabolites were obtained (Supplement, Table 7), which were enriched in 59 KEGG pathways, including 55 metabolic-type pathways, and 41 metabolites were enriched in metabolic pathways. The biosynthesis of secondary metabolites was enriched in 26 kinds of metabolites (Fig. 8B). Among the top 20 pathways showing significant enrichment, the citrate cycle (TCA cycle), sucrose metabolism, arginine and proline metabolism and pentose phosphate pathway and other pathways were involved in plant responses to drought stress and phosphorylation (Fig. 8C). Notably, 74 kinds of common metabolites were screened, which were enriched in 18 KEGG pathways, including 17 metabolism-type pathways. Metabolic pathways and the biosynthesis of secondary metabolites were enriched in 7 kinds of metabolites (Fig. 8D). Furthermore, various flavonoid biosynthesis, arginine and proline metabolism, and glycerophospholipid metabolism pathways participate in the metabolic and phosphorylation reactions of substances under drought stress (Fig. 8E). In addition, heat map analysis of the quantitative results for 74 kinds of metabolites showed that 44 kinds of differential metabolites presented high contents under control conditions, 6 kinds of differential metabolites presented high contents under low-phosphorus conditions, and 24 kinds of differential metabolites presented high contents under drought stress (Fig. 8F).
3.9 Differentially expressed protein analysis
A total of 470 protein-encoding genes were identified by the proteomic analysis of nine samples from the control, low-phosphorus and drought stress treatments (Supplement: Table 9). We analyzed the differentially expressed proteins among different treatments, and 71 differentially expressed proteins were identified in total (Fig. 9A), among which 7 differentially expressed proteins were found in the control group, including 3 upregulated and 4 downregulated proteins. There were 29 differentially expressed proteins related to drought control, including 16 upregulated and 13 downregulated proteins. There were 56 differentially expressed proteins identified under the low-phosphorus–drought treatment, 20 of which were upregulated and 36 of which were downregulated. In addition, 71 differentially expressed proteins were analyzed for enrichment. GO enrichment analysis showed that 41 differentially expressed proteins were related to metabolic processes, while 28 and 24 differentially expressed proteins were related to cellular processes and single organism processes, respectively (Fig. 9B). Neither immune system processes nor multiorganism processes were enriched in differentially expressed proteins. Among cellular components, the cell and cell part categories were enriched in 26 differentially expressed proteins, and the organelle category was enriched in 22 differentially expressed proteins. The membrane-enclosed lumen, nucleoid, virion and virion part categories were not enriched in differentially expressed proteins. Among molecular functions, binding and catalytic activities were enriched in 33 and 32 differentially expressed proteins, respectively. The antioxidant activity, molecular activity and protein tag categories were not enriched in differentially expressed proteins. The results of the KEGG pathway enrichment analysis showed that the ribosome pathway was enriched in 9 differentially expressed proteins, and the other 19 pathways all included 1 or 2 differentially expressed proteins. Sucrose metabolism and the phosphoric acid reaction pathway of carbon metabolism, starch and sucrose metabolism and pyruvate metabolism (Fig. 9C). In addition, heat maps were drawn based on the quantitative values of 71 differentially expressed proteins, 13 of which were high in the control group, and 22 of which were high in the low-phosphorus and drought groups. The quantitative values of most differentially expressed proteins under drought conditions were lower than those under the control and low-phosphorus treatments (Fig. 9D).