Tissue-cultured ‘M9-T337’ seedlings were transplanted to soil after rooting and grown for an additional 5 weeks in a greenhouse at 25°C (day) and 20°C (night), relative humidity 65–70%, light intensity 400 μM m⁻² s⁻¹, and a photoperiod of at 16 (light):8 h (dark), using full spectrum LED lamps. During this period, plants were watered weekly with modified 1/2-strength Hoagland nutrient solution [2.5 mM (NH4)₂Suc, 0.1 mM KH₂PO₄, 5.9 mM KCl, 1 mM CaCl₂·2H₂O, 2 mM MgSO₄, 0.1 mM Fe-EDTANa₂, 0.05 mM H₃BO₃, 0.001 mM ZnSO₄, 0.001 mM CuSO₄·5H₂O, 0.012 mM MnSO₄·H₂O, and 0.0002 mM Na₂M_OO₄·7H₂O) to ensure normal growth. Then, 180 plants of similar size with 8–9 leaves were selected to pre-culture in deionized water for about 10 d, allowing them to adapt to the hydroponic conditions and fully consume their stored nutrients. Next, the plants were divided randomly into four groups and separately exposed to (1) N+P-sufficient (NNNP, 5 mM KNO₃ and 1 mM KH₂PO₄), (2) P-deficient (NNLP, 5 mM KNO₃ and 0.001 mM KH₂PO₄), (3) N-deficient (LNNP, 0.2 mM KNO₃ and 1 mM KH₂PO₄), and (4) N+P-deficient (LNLP, 0.2 mM KNO₃ and 0.001 mM KH₂PO₄) conditions. The four nutrient solutions each contained 1 mM CaCl₂·2H₂O, 2 mM MgSO₄, 0.1 mM Fe-EDTANa₂, 0.05 mM H₃BO₃, 0.001 mM ZnSO₄, 0.001 mM CuSO₄·5H₂O, 0.012 mM MnSO₄·H₂O, and 0.0002 mM Na₂M_OO₄·7H₂O, and the pH was adjusted to 5.9. In the NNLP, LNHP, and LNLP treatments, KCl was used to replace KH₂PO₄ or KNO₃ to avoid K deficiency. In this experiment, 15 hydroponic boxes (1 L volume) were set up for each treatment. Each box contained 1 L of nutrient solution and three seedlings. The nutritional treatments were continued for 8 weeks and refreshed every 7 d. A submersible pump (power: 3W) was used to supply oxygen to the nutrient solution in each hydroponic box, and 30 minutes every two hours. For ¹⁵N determination, nine seedlings from each treatment were selected to treat with nutrient solutions that replaced potassium nitrate with ¹⁵N-KNO₃ (10.14%, Shanghai Research Institute of Chemical Industry, China).

After treatment, the roots were dispersed in water to separate them, and photographs were taken using a photo scanner (Epson, Japan). The images were analyzed using WinRHIZO Pro2009 (Regent Instruments, Canada) to determine the root parameters. Three biological replicates were assessed for each treatment.

The roots, stems, and leaves of three experimental replications from each treatment were harvested and the fresh weights measured. The components were then dried in an oven individually at 105°C for 30 min and then at 70°C to a constant weight. The dry weight was then measured.

The roots, stems, and leaves of ‘M9-T337’ seedlings labelled with ¹⁵N were washed three times with distilled water, fixed at 105°C for 30 min, and dried at 70°C to a constant weight before being ground into a powder. Then, 0.3 g of powdered material from each tissue type was digested using H₂SO₄-H₂O₂ and fixed in distilled water to a final volume of 50 mL. Five milliliters of the clear supernatant were removed to determine the total P content using the vanadium molybdate yellow colorimetric method (Zhang et al., 2020). The remaining 0.1 g of powdered material from each tissue type was used to determine the total N content and ¹⁵N abundance using a Flash 2000HT elemental analyzer coupled to a Finnigan DELTA V Advantage isotope ratio mass spectrometer (Thermo Fisher Scientific, Germany). The total N contents or ¹⁵N contents of the leaf, stem, and root were calculated by multiplying the N concentration or ¹⁵N concentration by the dry weight, respectively. ¹⁵N utilization efficiency (¹⁵NUE) was calculated as the ratio of the total ¹⁵N content in the seedling to the total ¹⁵N in the fertilizer.

On Day 7 of the experimental period, three seedlings from each treatment were used for ion flux measurements. The net H⁺ and NO 3   − flux in the rhizosphere was measured with a noninvasive micro test technique (NMT 100 Series, Younger USA LLC, USA). Ionic flux data were calculated using Mage Flux (http://xuyue.net/mageflux).

On Day 30 of the experimental period, the roots from each treatment were harvested and divided into two portions, one for widely targeted metabolomics determination and the other for transcriptome sequencing. Extraction, derivatization, detection, quantification, and data analyses were carried out by Wuhan MetWare Biotechnology Co., Ltd., China. Three biological replicates were assessed for each treatment, and one biological replicate was used for every nine seedling root mixes. The metabolites were extracted with a 70% methanol solution from 0.1 g of the lyophilized powdered samples, and, subsequently, the sample extracts were analyzed using an UPLC-ESI-MS/MS system (Applied Biosystems 4500 Q TRAP, Germany; SHIMADZU Nexera X2, MS, Japan).

The metabolites in the different samples were analyzed based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) compound database, MetWare database (MWDB), and multiple reaction monitoring (MRM). To study metabolite accession-specific accumulation, metabolites from 12 samples were used for hierarchical clustering analysis (HCA) and orthogonal partial least squares discriminant analysis (OPLS-DA) after being log₂-transformed and normalized. The metabolites were selected on the basis of the combination of a statistically significant threshold of variable influence in projection (VIP) values obtained from the OPLA-DA model and P values from a two-tailed Student’s t-test on the normalized peak areas from different groups, where metabolites with VIP > 1.0 and P< 0.05 were considered differentially accumulated metabolites (DAMs). Venn diagrams were used to illustrate the number of DAMs among different treatments. In order to study the trend of relative metabolite content under different N and P supply conditions, the relative content of all DAMs was standardized to z-core, and then Kmeans cluster analysis was conducted based on z-core. Z-core was calculated by equation (1). x represents the quantitative value of a specific metabolite, μ represents the average quantitative value of all metabolites, and σ represents the standard deviation.

Total RNA of the samples was extracted with an RNAprep Pure Plant kit (DP441, Tiangen, China). The high-quality mRNA was randomly fragmented. First-strand cDNA was synthesized using the M-MuLV reverse transcriptase system. The RNA strand was then degraded by RNase H, and second-strand cDNA was synthesized using DNA polymerase. The double-stranded cDNAs were ligated to sequencing adapters. After amplification and purification, cDNA libraries were obtained and sequenced using a Novaseq6000 system (Illumina, USA).

Gene expression levels were determined using the RPKM (reads per kb per million reads) method. Differential expression analysis was performed using DESeq2 with a design formula [|log₂(FoldChange)|≥ 1 and false discovery rate (FDR)< 0.05] that took into account the contrast among samples from four treatments. The genes/transcripts with an adjusted P value of<0.05, using multi-testing error correction (Benjamini-Hochberg test) and corresponding to a false discovery rate (FDR) of 1%, were considered to be differentially expressed genes (DEGs). All DEGs were mapped to pathway terms in the KEGG database (http://www.genome.jp/kegg/). The FPKM values after centralization and standardization were used for DEGs hierarchical clustering analysis, and the clustering heatmap was drawn.

Quantitative values of genes and metabolites in all samples were used for correlation analysis using R (corrplot, Version 0.84). Correlation results with correlation coefficient greater than 0.8 and p value less than 0.05 were selected to draw nine-quadrant heatmap and correlation heatmap.

Purified RNA was reverse transcribed to first-strand cDNA with a PrimeScript™ RT Master Mix cDNA Reverse Transcription Kit (Takara, China). qRT-PCR was conducted with a ChamQ SYBR qPCR Master Mix kit (Vazyme, China) and a C1000 Touch™ Thermal Cycler system (Bio-Rad, USA). Relative transcript levels were calculated according to the 2^−ΔΔCp method using MdActin as a reference. Three biological and technical replications were performed. The primes sequences are shown in Table S1 .

The open reading frames of MdEXPA4 and MdEXLB1 were respectively inserted into the pRI 101-AN vector (TaKaRa, China), and the recombinant plasmids were respectively transformed into Agrobacterium strain GV3101. The recombinant plasmids were transformed into Col-0 using floral-dip transformation. T0 seeds were sown on 1/2 MS medium containing Kanamycin (50 mg/L) for transgenic selection. Two homozygous T3 transgenic A. thaliana lines (#1 and #2) were selected for further study. The MdEXLB1 plasmids were transfected into Solanum lycopersicum (cultivar Ailsa Craig) plants through leaf explant infection. Three third generation homozygous transgenic lines (#5, #8, and #11) were selected for further study.

Seeds of wild-type (WT) and MdEXLB1-OE S. lycopersicum were sown in MS medium. After sprouting, the seedlings were transplanted to soil and grown for an additional 2 weeks. Then, seedlings of similar size were randomly divided into four groups and separately exposed to N+P-sufficient, P-deficient, N-deficient, and N+P-deficient conditions, as in section 2.1. The nutritional treatments were continued for 8 weeks and refreshed every 7 d.

The seeds of wild-type and MdEXPA4-OE A. thaliana were sown in MS medium. After sprouting, the seedlings were grown on MS medium with different concentrations of NO 3   − and HPO 4   2 − : (1) N+P-sufficient (NNNP, 60 mM NO 3   − and 1.25 mM HPO 4   2 − ), (2) P-deficient (NNLP, 60 mM NO 3   − and 0.01 mM HPO 4   2 − ), (3) N-deficient (LNNP, 0.2 mM NO 3   − and 1.25 mM HPO 4   2 − ), and (4) N+P-deficient (LNLP, 0.2 mM NO 3   − and 0.01 mM HPO 4   2 − ) conditions.

All physiological data were graphed using GraphPad Prism v6.01 (GraphPad Software Inc., USA). Duncan’s test (P< 0.05) was used to analyze the statistical significance by SAS v9.0 (SAS Institute Inc. USA).

Compared with NNNP, the above-ground growth of ‘M9-T337’ seedlings was inhibited under NNLP, LNNP, and LNLP conditions ( Figure 1 ). Relative to NNNP, the leaf fresh weight decreased by 57.64%, 68.16%, and 81.24% under NNLP, LNNP, and LNLP conditions, respectively, and the leaf dry weight decreased by 48.10%, 46.32%, and 76.49%, respectively. In addition, the stem fresh weight decreased significantly only under LNLP, whereas the stem dry weight decreased by 29.46% and 58.63% under NNLP and LNLP conditions, respectively. However, the root fresh weight increased 1.78-fold under LNLP, and the root dry weight increased 1.35- and 1.94-fold under LNNP and LNLP conditions, respectively. Overall, the total plant fresh weight decreased in the order NNNP > NNLP > LNNP > LNLP, while the total plant dry weight decreased by 36.42% and 57.55% under NNLP and LNLP, respectively ( Figures 1B, C ). Based on the dry weights of the shoot and root, 1.86-fold, 2.14-fold, and 5.71-fold increases were found in the root-to-shoot ratio under NNLP, LNLP, and LNLP conditions, respectively ( Figure 1D ).

Taking into account the root parameters, NNLP, LNNP, and LNLP significantly affected the total root length, volume, surface area, average diameter, and number of tips. The total values for root length, volume, and surface area decreased among treatments in the following order: NNLP > LNLP > LNNP > NNNP, and the differences were significant ( Figures 2A–C ). The root average diameter increased 1.13 times under LNNP, whereas it reduced by 11.36% and 24.40% under NNLP and LNLP, respectively ( Figure 2D ). The greatest increase in root forks was found under the LNLP condition ( Figure 2E ). Additionally, 19.20%, 25.33%, and 39.18% increases in root tips were observed in NNLP, LNNP, and LNLP, respectively ( Figure 2F ).

Compared with that in NNNP, the total leaf N contents in NNLP, LNNP, and LNLP were significantly reduced by 18.87%, 77.67%, and 75.94%, respectively; the total stem N content was increased 1.15-fold in LNNP, but reduced by 15.87% in LNLP; and the total root N contents in NNLP, LNNP, and LNLP were increased 1.39-fold, 1.30-fold, and 1.72-fold. Based the total N content in all organs, the ranking for the total N content in whole plants under different treatments was NNNP > NNLP > LNNP ≌ LNLP ( Figure 3A ). Similarly, compared to those grown under NNNP conditions, the ¹⁵N contents under NNLP, LNNP, and LNLP conditions were reduced by 46.87%, 95.39%, and 91.75% in leaf, respectively. In stem and root, the ¹⁵N contents decreased significantly under LNNP and LNLP, but the greatest decreases were found in the LNNP treatment (75.54% and 62.57%); the decreases under NNLP were not significant ( Figure 3C ). Based on the ¹⁵N contents of different organs, 3.30-fold and 4.60-fold increases were found in the ¹⁵NUE under LNNP and LNLP, but no difference was found under NNLP ( Figure 3D ). NNLP caused a reduction of NO 3   − influx in comparison with NNNP, while there was much less NO 3   − influx under LNNP and LNLP in comparison with NNNP and NNLP ( Figures 3E, G ).

The total leaf P content by treatment decreased in the order: NNNP > NNLP > LNNP > LNLP, and the differences were significant. The greatest decrease in the total stem P content was observed in the LNLP treatment (64.03%). Relative to NNNP, the total root P content decreased by 31.80% under NNLP but increased by 29.90% and 16.83% under LNNP and LNLP, respectively. The total P contents in whole plants under NNLP, LNNP, and LNLP decreased, with the greatest decrease found in the LNLP treatment (62.08%) ( Figure 3B ). The maximum H⁺ influx was detected at 400 μm from the apical end of the root tip in roots of ‘M9-T337’ seedlings after growth under NNNP conditions for 7 d ( Figure 3F ). However, obvious H⁺ efflux was detected in the same region of the roots of ‘M9-T337’ seedlings under NNLP, LNNP, and LNLP treatments. The maximum H⁺ efflux was detected in the roots of ‘M9-T337’ seedlings under the LNNP treatment ( Figure 3H ).

A total of 1,069 metabolites were obtained, and these substrates were classified into 11 categories at Class 1 level; specifically: 216 flavonoids, 227 phenolic acids, 152 lipids, 96 amino acids and derivatives, 85 organic acids, 69 nucleotides and derivatives, 57 alkaloids, 41 triterpenes, 40 lignans and coumarins, 19 tannins, and 97 others ( Figure S1 ). To study the metabolite accumulation patterns, a K-means cluster analysis was performed. The results showed the metabolites were divided into five subclasses ( Figure 4A ). Metabolites in subclass 1 showed an increasing trend under LNNP, and the content of metabolites in subclass 3 tended to increase under NNLP and LNLP compared with NNNP, while subclass 4 metabolites showed a decreasing trend under LNLP, and those in subclasses 2 and 5 exhibited a decreasing trend under NNLP, LNNP, and LNLP in comparison to NNNP. Compared with NNNP, 264 DAMs were observed under NNLP, of which 21.97% were significantly up-regulated and 78.03% exhibited significant down-regulation. The 313 DAMs were observed under LNLP, of which 38.02% were up-regulated and 72.52% were down-regulated. In addition, 229 DAMs were observed under LNNP, of which 17.90% were up-regulated and 82.10% were down-regulated ( Figure 4B ). Venn analysis revealed 79 common DAMs in the comparisons of NNNP vs NNLP, NNNP vs LNNP, and NNNP vs LNLP. In addition, 78 DAMs were shared in comparisons of the NNNP vs NNLP and NNNP vs LNLP treatments, 22 were shared in the comparisons of NNNP vs NNLP and NNNP vs LNNP treatments, and 47 were shared in the comparisons of NNNP vs LNNP and NNNP vs LNLP treatments ( Figure 4C ).

To further analyze genome-wide changes in gene expression, the 12 samples were used to construct RNA libraries and perform RNA-seq. A total of 2025 DEGs were identified in the NNNP vs NNLP comparison, of which 52.40% were up-regulated and the rest were down-regulated; a total of 2851 genes were differentially expressed in the NNNP vs LNNP comparison, of which 61.24% were up-regulated, and the rest were down-regulated; a total of 1963 genes exhibited significant differences in expression in the NNNP vs LNLP comparison, of which 52.29% were up-regulated, and the rest were down-regulated ( Figure 5A ). Five hundred DEGs were common in the comparisons of NNNP vs NNLP, NNNP vs LNNP, and NNNP vs LNLP ( Figure 5B ). Compared with NNNP, the changes in gene transcription levels under NNLP and LNLP were greater than those under LNNP ( Figure 5C ).

KEGG enrichment analysis showed that DEGs and DAMs in the comparison of NNNP vs NNLP were association with 62 co-mapped pathways; those in the comparison of NNNP vs LNNP were associated with 49 co-mapped pathways; and those in the comparison of NNNP vs LNLP were associated with 63 co-mapped pathways. Metabolism-related assignments such as metabolic pathways, biosynthesis of secondary metabolites, biosynthesis of cofactors, biosynthesis of amino acids, flavonoid biosynthesis, 2-oxocarboxylic acid metabolism, phenylpropanoid biosynthesis, carbon metabolism, and purine metabolism, as well as environmental information processing-related pathways, including ABC transporters, were conspicuously enriched pathways under N and/or P deficient conditions ( Figure S2 ). According to the Pearson correlation coefficients (PCCs) for genes and metabolites detected in each sample, the metabolite changes may have been positively regulated by genes in the third and seventh quadrants, while the patterns of genes and metabolites from the first and ninth quadrants showed the opposite ( Figures 6A–C ). The DAMs and DEGs with PCCs greater than 0.8 are shown in Figures 6D–F . Respectively, 256 DAMs and 1,915 DEGs, 228 DAMs and 2,656 DEGs, and 306 DAMs and 1,893 DEGs were screened in the comparisons of NNNP vs. NNLP, NNNP vs. LNNP, and NNNP vs. LNLP. The DAMs under NNLP and LNLP conditions were distributed in the 11 groups of Class I, while the DAMs under LNNP were distributed among 10 groups. Notably, tannin substances such as procyanidin C2, strictinin, corilagin, and Sanguiin H1 were not screened under LNNP.

Association analysis based on the KEGG database showed that, in the comparison of NNNP vs. NNLP, the 68 DAMs and 413 DEGs with PCCs greater than 0.8 were mapped to 36 KEGG pathways; in the comparison of NNNP vs. LNNP, 59 DAMs and 521 DEGs with higher PCCs were mapped to 23 KEGG pathways in the comparison of NNNP vs. LNNP; and 80 DAMs and 360 DEGs with higher PCCs were mapped to 30 KEGG pathways in the comparison of NNNP vs. LNLP ( Figure S3 , Table S2 ). Among them, pentose and glucuronate interconversions, starch and sucrose metabolism, galactose metabolism, amino sugar and nucleotide sugar metabolism, and phenylpropanoid biosynthesis were conspicuously enriched pathways ( Figure S4 ).

The content of D-Glucoronic acid and D-Arabinose, which involved in pectin synthesis, were reduced under LNNP, in addition, the content of Glucose-1-phosphate, Uridine 5’-diphospho-D-glucose, and D-Galacturonic acid, which related to pectin metabolism, tended to decrease under NNLP, however, the metabolites in pentose and glucuronate interconversions pathway which related to pectin metabolism did not change significantly under LNLP. The metabolites in starch and sucrose metabolism pathway and galactose metabolism pathway, and involved in hemicellulose metabolism such as D-Glucose 6-phosphate, D-Fructose 6-Phosphate, Glucose-1-phosphate, Uridine 5’-diphospho-D-glucose, D-Glucose 1,6-bisphosphate, and Dihydroxyacetone phosphate were all down-regulated under NNLP, moreover, D-Glucose 6-phosphate, D-Fructose 6-Phosphate, and Raffinose were down-regulated under LNLP, furthermore, D-Mannose and Dulcitol were down-regulated under LNNP, but D-Glucose 1,6-bisphosphate was up-regulated under LNNP. The metabolites involved in lignin metabolism pathway such as Syringin and Caffeoylquinic acid were down-regulated under NNLP and LNLP, Ferulic acid, Sinapic acid, and Sinapyl alcohol were up-regulated under NNLP and LNLP, Coumarin and Cinnamaldehyde were down-regulated under LNLP. Under LNNP, Cinnamic acid and Sinapine were down-regulated significantly. The content of amino acid that made up extension, such as L-lysine showed a decreasing trend under LNLP, and the content of amino acid that made up expansin such as L-tryptophan was decreased under LNLP. The content of above amino acids did not change significantly under NNLP and LNNP ( Figure 7 ).

A total of 261 non-redundant DEGs encoding cell wall proteins or enzymes involved in cell wall metabolism were screened, of which 44 were involved in hemicellulose metabolism, 10 in cell wall synthesis, 129 in lignin metabolism, 44 in glycoprotein synthesis, and 41 in pectin metabolism ( Figure 8 ).

Overall, 25 (20 up-regulated and five down-regulated), 27 (13 up-regulated and 14 down-regulated), and 18 (12 up-regulated and six down-regulated) genes involved in cellulose and hemicellulose metabolism were identified in the comparisons of NNNP vs. NNLP, NNNP vs. LNNP, and NNNP vs. LNLP, respectively. Five, two, and four genes encoding xyloglucan endotransglucosylase/hydrolase (XTH) were up-regulated in comparisons of NNNP vs NNLP, NNNP vs LNNP, and NNNP vs LNLP, respectively. The expression of genes encoding type II glycosyltransferases, such as MD03G1126500 (IRX7), MD11G1146800 (IRX7) MD04G1103300 (IRX9), MD15G1287000 (IRX15), MD12G1123400 (IRX9), MD15G1287000 (IRX15), and MD02G1177100 (IRX15-L), showed an increasing trend in three comparisons, while novel.2662 (GUX2) was down-regulated in the comparison of NNNP vs LNNP. Seven of ten DEGs encoding type III polysaccharide synthases, including MD00G1061100 (CESA4), MD05G1189700 (CESA4), MD05G1296600 (CESA8), MD10G1276500 (CESA8), MD02G1095700 (CSLE2), MD15G1415100 (CSLG2), and MD15G1415200 (CSLG2), were up-regulated in the comparison of NNNP vs LNNP ( Figure 8A ).

A total of 58 (27 up-regulated and 31 down-regulated), 66 (29 up-regulated and 37 down-regulated), and 63 DEGs (28 up-regulated and 35 down-regulated) involved in lignin metabolism were screened in the comparisons of NNNP vs NNLP, NNNP vs LNNP, and NNNP vs LNLP, 13 DEGs, including MD03G1051400 (CYP73A5), MD17G1224900 (HCT), MD12G1211500 (BGLU), MD15G1374300 (BGLU), MD07G1161100 (COMT), MD07G1301000 (COMT), MD03G1288300 (COMT), MD17G1119700 (CAD), MD00G1088400 (PRX), MD02G1124700 (PRX), MD02G1138900 (PRX), MD03G1014200 (PRX), and MD13G1011000 (PRX), were down-regulated in all three comparisons, and nine DEGs, MD04G1016900 (HCT), MD07G1095200 (HCT), MD09G1169600 (HCT), MD09G1267700 (HCT), MD10G1161400 (HCT), MD10G1199500 (HCT), MD06G1038800 (PRX), MD13G1053000 (PRX), and MD15G1238800 (PRX), were up-regulated in all three comparisons ( Figure 8B ).

Six, two, and five expansin family genes were up-regulated in the comparisons of NNNP vs NNLP, NNNP vs LNNP, and NNNP vs LNLP, respectively. All DEGs in the comparison of NNNP vs NNLP were up-regulated; MD11G1054500 (EXPA4) and MD06G1041000 (EXLB1) were shared in three comparisons. The qRT-PCR results were generally consistent with the RNA-seq data ( Figure 8E ). A total of 18 extensin genes were differentially expressed, of which MD17G1257500 (PERK13) was up-regulated in the three comparisons. Three genes encoding arabinogalactan protein, MD03G1004000 (AGP14), MD01G1076000 (AGP16), and MD02G1072000 (AGP20), were up-regulated under NNLP and LNLP, but not differentially expressed under LNNP, while two genes, MD10G1293100 (AGP30) and MD02G1314100 (AGP23), were down-regulated under LNNP. Two genes encoding proline-rich proteins, MD07G1225500 (PRP1) and MD01G1156800 (PRP1), were down-regulated under LNNP ( Figure 8C ).

Three, seven, and three genes encoding pectinesterase (PME) were down-regulated in the comparisons of NNNP vs. NNLP, NNNP vs. LNNP, and NNNP vs. LNLP, respectively, of which MD12G1018900 (PME58) was the common DEG. Three, one, and two genes were up-regulated in the comparisons of NNNP vs NNLP, NNNP vs LNNP, and NNNP vs LNLP, respectively. A total of seven (four up-regulated and three down-regulated), 15 (three up-regulated and 12 down-regulated), and six (two up-regulated and four down-regulated) encoding polygalacturonase (PG) were identified in the three comparisons NNNP vs. NNLP, NNNP vs. LNNP, and NNNP vs. LNLP, of which MD12G1077600 (PG) was up-regulated and MD17G1279500 (PG) was down-regulated in the three comparisons. Compared with NNNP, the expression profile of a gene encoding pectate lyase (PL), MD13G1231000 (PL8), was up-regulated under NNLP, while three PL genes, MD01G1100600 (PL18), MD15G1149300 (PL18), novel.1309 (PL18), were down-regulated under LNNP, and the remaining two genes, MD05G1350300 (PL16) and MD09G1019900 (PL), were down-regulated under LNLP ( Figure 8D ).

Although MdEXPA4 transcription was not induced under NNNP treatment ( Figure 8E ), this gene was expressed under the other three treatments. To clarify the role of MdEXPA4 in the response to N and/or P deficiency, the phenotypes of wild-type (WT) and MdEXPA4-overexpressing A. thaliana subjected to NNNP, NNLP, LNNP, and LNLP conditions were analyzed. The results showed that MdEXPA4-overexpressing lines exhibited longer root length and more root forks than the WT when treated with NNNP ( Figures S8E, F ). Under NNLP conditions, the root lengths of transgenic lines #1 and #2 were 1.70-fold and 1.72-fold higher than that of the WT, respectively ( Figure S8E ). LNNP treatment led to 1.63-fold, 1.75-fold and 9.73-fold, 12.47-fold increases in the number of root forks and tips in transgenic lines #1 and #2, respectively, compared to the WT ( Figures S8F, G ). The LNLP condition caused the root lengths of transgenic lines #1 and #2 to increase by 1.34-fold and 1.54-fold compared to the WT, respectively, while the numbers of root tips also increased significantly ( Figures S8E, G ).

Another gene, MdEXLB1, was hypothesized to be an important candidate for involvement in the response to N and/or P deficiency ( Figure 8E , Figure S10 ). To evaluate the function of MdEXLB1, WT and MdEXLB1-overexpressing S. lycopersicum plants were cultivated hydroponically under different P and/or N supply conditions. The results showed that under the NNNP treatment, the fresh weights of transgenic lines #5, #8, and #11 were 1.78-fold, 1.79-fold, and 2.07-fold higher than that of WT, and the root fresh weights of these transgenic plants were 3.50-fold, 3.83-fold, and 4.19-fold higher than that of the WT, respectively. Under NNLP, LNNP, and LNLP conditions, the fresh weights of the transgenic plants were significantly higher than that of the WT; in particular, the roots exhibited a significant increase in fresh weight compared with those of the WT ( Figure 9 ). The total length, surface area, and volume of roots in MdEXLB1-overexpressing lines were significantly higher than those of the WT under NNNP. Under NNLP, LNNP, and LNLP, the root surface areas of the transgenic plants were significantly higher than that of the WT ( Figure S9 ). Moreover, NNNP treatment led to 1.26-fold, 1.53-fold, and 1.76-fold and 1.41-fold, 1.39-fold, and 1.28-fold increases in the total P and total N accumulation in transgenic lines #5, #8, and #11, respectively. Under the NNLP, LNNP, and LNLP treatments, the values for total P and total N accumulation were significantly higher than those in the WT ( Figure S9 ).

In field conditions, the availability of Pi and NO 3   − is one of the most important limiting factors for apple development (Sun et al., 2021a; Sun et al., 2021b; Chai et al., 2022). The availability of these nutrients impacts the phenome, ionome, transcriptome, and metabolome composition of the plants (Gan et al., 2016; Sun et al., 2020b; Lv et al., 2021; Mu and Chen, 2021; Nezamivand-Chegini et al., 2022; Wen et al., 2022). The paramount relevance of the rootstock is its key role at the interface between the apple plants and soil for nutrient mining (Chauhan et al., 2020). In the present study, we observed that the aboveground growth and development of apple was inhibited by N and/or P deficiency ( Figure 1 ), and the partitioning of carbohydrate and mineral nutrient resources to root biomass tended to increase ( Figure 3A ), thus leading to an increase in the root-to-shoot ratio of the hydroponically-grown dwarfing rootstock ‘M9-T337’ ( Figure 1D ), which was consistent with previous results (Sun et al., 2021a; Sun et al., 2021b). Generally, root development was enhanced by N and/or P deficiency, although P deficiency promoted root elongation, while N deficiency did not, and similarities existed between the root morphology under N+P deficit condition and that under P deficit and N deficit conditions, respectively, but the effect of N+P deficiency on root morphogenesis was not the accumulation of effects of P deficiency and N deficiency, and compared with P deficiency, root development was weakened under N+P deficiency ( Figures 1 , 2 ). These results suggested that the ‘M9-T337’ rootstock is highly plastic and responsive to the changing soil environment (Ma et al., 2020; Tahir et al., 2021a; Tahir et al., 2021b; Zhang et al., 2021b).

Previous studies have shown that N deficiency increased H⁺-ATPase activity and H⁺ efflux, therefore affecting lateral root (LR) development and NO 3   − uptake (Mlodzinska et al., 2015; Jiang et al., 2017; Lv et al., 2021). In a study of O. sativa, Wang et al. (2021b) found that H⁺ efflux was important for root elongation and P uptake under P deficient conditions. Moreover, Xu et al. (2020) confirmed the important role of the H⁺ pump in white lupin adaptation to P deficiency. To measure the concentration gradient of NO 3   − and H⁺, while preserving the integrity of the sample, NMT was used in this study (Huang et al., 2017; Lv et al., 2021). The ionome data indicated that, although the reduced net NO 3   − influx rate resulted in decreased total ¹⁵N accumulation in plants, N deficiency and N+P deficiency improved the ¹⁵NUE compared to a sufficient N+P supply ( Figure 3 ). We monitored obvious H⁺ efflux in the roots under N-deficient and N+P-deficient conditions ( Figure 3H ), which resulted in a decrease in the pH of the nutrient solution ( Figure S7 ). These results confirmed that H⁺ pumps stimulated by N deficiency play a crucial role in NO 3   − absorption of apple rootstock. The H⁺ efflux monitored under P deficiency implied that the H⁺ pumps were key players in the response to P deficiency.

N transporters (NRT) and P transporters (PHT) play important roles in plant uptake and transport of N and P, respectively (Zhang et al., 2021c; Chai et al., 2022). NRTs in higher plants include low affinity transporters encoded by some members of the NRT1/PTR (also named as NPF) family and high affinity transporters encoded by the NRT2 gene family (Tahir et al., 2021c; Misawa et al., 2022), as well as NRT1.1, a dual-affinity NO 3   − transporter (Ye et al., 2019). Four NTR2 transporters (NRT2.1, NRT2.2, NRT2.4, and NRT2.5) and two NRT1 transporters (NPF6.3 and NPF4.6) have root NO 3   − uptake activities (Lezhneva et al., 2014). In this study, 20 NRTs (NRT2.5, NRT2.7, and 18 NRT1s) were up-regulated under N deficiency, and 10 of them were also induced by N+P deficiency ( Figure S5 ), and therefore may be high affinity NO 3   − transporters. Lezhneva et al. (2014) reported that NRT2.5 plays an essential role in response to severe N deficiency. NRT2.7 is a seed-specific high-affinity nitrate transporter controlling NO 3   − content in A. thaliana (David et al., 2014; Zheng et al., 2016). In this study, it was found that NRT2.5 and NRT2.7 may be involved in NO 3   − uptake of ‘M9-T337’ root under N deficiency and P deficiency, respectively ( Figure S5 ). AtNRT1.5 drives root-to-shoot transport of NO 3   − and is also an indole-3-butyric acid transporter, involved in LR development (Watanabe et al., 2020; Chen et al., 2021). Gho and Jung (2019) found that AtNRT1.5 was induced by Pi deficiency in A. thaliana; they demonstrated that atnrt1.5 mutants displayed conspicuously longer PRs along with significantly reduced LR density under P deficient conditions. Therefore, NRT1.5 may also be involved in apple root architecture modification under N and/or P deficient conditions. Chai et al. (2022) found overexpression of MdNRT2.4 enhanced net H⁺ extrusion from the root surface and NO 3   − absorption. In this study, three MdNRT2.4 genes were only significantly upregulated under N+P deficient condition; thus, they may have important roles in N absorption under combined nutritional deficiencies. Additionally, 19 NRT1s were down-regulated in the root of ‘M9-T337’ under N deficiency, and 11 of them were also down-regulated under N+P-deficient conditions ( Figure S5 ), implying they may be low affinity transporters, which do not act at low NO 3   − concentrations (Sinha et al., 2020). Similarly, Hu et al. (2011) found that plants repress N assimilation and save more energy for Pi acquisition under P deficiency ( Figure S5 ), which is in agreement with our results. The low affinity NO 3   − transporters were down-regulated in response to P deficiency (Nezamivand-Chegini et al., 2022).

Seven PHT1s were upregulated in response to P deficiency and N+P deficiency ( Figure S5 ). The PHT1 subfamily contains high affinity Pi transporters, which function under low P concentrations in soil (Gho and Jung, 2019). Moreover, SPX1, SPX3, and SPX5 were up-regulated under P deficiency and N+P deficiency ( Figure S5 ). Studies have shown that most SPX subfamily members are induced by Pi deficiency in plants (Stefanovic et al., 2007; Zhang et al., 2016; Du et al., 2017; Wang et al., 2021a), and their major role is to interact with PHR1, thereby regulating the transcription of downstream PSIs (Lv et al., 2014; Dong et al., 2019). Furthermore, PHO1 was down-regulated under three treatments. According to previous studies, its homologs AtPHO1, AtPHO1;H1, and OsPHO1;2 are involved in regulating internal Pi homeostasis by transporting Pi from the root xylem to shoots (Secco and Poirier, 2010; Stefanovic et al., 2011), implying the translocation of Pi from the roots to shoot is inhibited by N and/or P deficiency. The PHTs exhibited a lower transcript level under N+P deficiency, and five PHTs (two PHT1 and three PHT4 genes) were down-regulated under N deficiency. The expression of SPXs was affected by N availability, it was speculated that NO 3   − deficiency may posted a negative effect on P absorption in apple dwarfing rootstock (Hu et al., 2019; Medici et al., 2019; Zhang et al., 2021e).

The number of DAMs specifically responding to N deficiency and P deficiency was fewer than the number responding to N+P deficiency ( Figure 4 ). The number of DEGs responding specifically to P deficiency and N+P deficiency was almost same, but fewer than the number that responded to N deficiency ( Figure 5 ). Thus, greater disturbance to metabolites occurred in response to combined nutritional deficiencies than a single nutrient deficiency, and the molecular mechanisms in response to combined nutritional deficiencies were rather complex. Several pathways involved in cell wall reprogramming were conspicuously enriched. However, ‘M9-T337’ adapts it root development differently under N and/or P deficiency, due to the variation of metabolites related to cell wall synthesis was different ( Figure S2 ). The plant cell wall is a dynamic network composed of cellulose, hemicellulose, pectin, lignin, and multiple types of structural proteins (Lampugnani et al., 2018). Under N deficiency, the plants altered their cell wall organization and adapted their root architecture by elongating the LRs to maximize N acquisition (Sun et al., 2017a). Kumar et al. (2021) revealed that cell wall reorganization, and associated activity-related gene up-regulation, is a strategy for tolerating P deficiency in O. sativa. Cell wall synthesis, composition, and reprogramming contribute to the meticulous modulation of root architecture and plant adaptations to biotic and abiotic stress (Ogden et al., 2018).

Cellulose is a load-bearing polymer present in the plant cell wall (Zhang et al., 2021a). Cellulose synthase complexes (CESAs) are involved in the formation of both primary and secondary cell walls (Endler and Persson, 2011). The roles of different subunits of CESAs vary, with CESA4, CESA7, and CESA8 involved in secondary cell wall formation in A. thaliana (Endler and Persson, 2011). In this study, CESA4 and CESA8, which are highly homologous to AtCESA4 and AtCESA8, were highly induced by N deficiency ( Figure 8 ). Therefore, N deficiency could increase cellulose formation in the root of ‘M9-T337’. A similar result was also reported by Nezamivand-Chegini et al. (2022). However, under P deficiency and N+P deficiency, no distinct cellulose synthesis was discerned. Hemicellulosic polysaccharides, including xylans, xyloglucans, mannans, glucomannans, and β-(1,3;1,4)-glucan, all harbor β-(1,4)-glycosyl-linked backbones with similar equatorial configurations (Zhang et al., 2021a). Xylans are the most abundant hemicellulosic polymers in vascular plants (Smith et al., 2017). IRX10 and its homolog IRX10L (members of the GT47 family), as well as IRX9, IRX9L, IRX14, and IRX14L (members of the GT43 family), are key components required for xylan backbone synthesis (Wu et al., 2010; Chen et al., 2013). IRX9, IRX10, and IRX14 are mainly involved in the extension of the backbone, while GUX1 and GUX2 are able to synthesize almost all side chains (Zhang et al., 2021a). In addition, IRX7 (a member of the GT47 family), IRX15, and its homolog IRX15L (members of the GT8 family), are also involved in xylan synthesis (Lee et al., 2009; Brown et al., 2011; Jensen et al., 2011). The transcriptome data showed that under all N deficiency, P deficiency, and N+P deficiency treatments, IRX7, IRX9, IRX10, IRX15, and its homolog IRX15L, were up-regulated, while the expression of GUX2 showed a decreasing trend under LNNP. These results suggested that xylan backbone synthesis in root of ‘M9-T337’ was enhanced by N and/or P deficiency, but side chain synthesis was inhibited by N deficiency. Xyloglucan is the most abundant primary-wall hemicellulose in all spermatophytes except grasses (Zhang et al., 2021a). XTHs are involved in the rearrangement of xyloglucans by cutting and rejoining them, causing reversible or irreversible loosening of cell walls to permit cell expansion (Nezamivand-Chegini et al., 2022). Zhang et al. (2022) reported that GmXTH38 is induced by P deficiency in both roots and leaves of soybean (Glycine max), and overexpression of GmXTH38 increases the sensitivity of LR formation under P deficiency. In this study, the transcription levels of XTHs in roots of ‘M9-T337’ were not only induced by P deficiency and N+P deficiency, but also by N deficiency ( Figures 7A , 8A ). These results showed that N and/or P deficiency had a role in the formation of xyloglucans in primary-wall hemicellulose of root cells, thereby mediating LR formation (Bourquin et al., 2002).

Pectins are a group of cell wall polysaccharides that possess α‐(1,4)‐linked galacturonic acids in their backbone. They play important roles in cell wall remodeling. Pectin deesterification stiffens the cell wall and prevents cell wall loosening (Hocq et al., 2017). Previous studies have shown that the enzymes involved in pectin synthesis are co-regulated together with NO 3   − carriers (Landi and Esposito, 2017). In addition, pectin has been demonstrated to contribute greatly to P remobilization from the cell wall. Elevated PME activity can facilitate the remobilization of P deposited in the cell wall (Zhu et al., 2016). Wu et al. (2022) found that OsPME14 overexpression increased PME activity, released more P from the root cell wall, and improved plant resistance to P deficiency. Zhu et al. (2018) reported that in O. sativa, NO 3   − deficiency caused higher levels of nitric oxide to accumulate in the roots and increased the activity of nitrate reductase. The nitric oxide stimulated cell wall pectin synthesis and demethylation of pectin, therefore increasing cell wall P release by increasing PME activity under P deficiency. Unlike O. sativa, N and/or P deficiency were unfavourable to the biosynthesis of cell wall pectin in the root of ‘M9-T337’ ( Figures 7A , 8A ), the effect of N deficiency was more significant since more PMEs, PGs, and PLs were significantly down-regulated under N deficiency ( Figure 8A ), suggesting that the response mechanism is unique among different species. Nevertheless, the expression of several PME genes showed an increasing trend under N and/or P deficiency ( Figure 7A ), thus the release of cell wall P may be increased.

Lignin is an unordered polymer composed of phenylalanine-derived aromatic monomer substances (Vanholme et al., 2019), and it is the primary structural component of thickened secondary cell walls in vascular plants (Bonawitz and Chapple, 2010). Lignin is composed of p-hydroxyphenyl units without methoxyl group (synthesized by the polymerization of p-coumaryl), guaiacyl units with one methoxyl group (synthesized by the polymerization of coniferyl), and syringyl units with two methoxyl groups (synthesized by sinapyl alcohol) (Xie et al., 2018; Gou et al., 2019). Qin et al. (2019) reported that N deficiency reduces the activity of most peroxidases, inhibits lignin biosynthesis, and lowers root solidity. Similar results were found in present study. Most PRX genes which encoded peroxidase and were involved in syringyl lignin biosynthesis were down-regulated under N deficiency and the content of coniferaldehyde and coniferyl alcohol tended to decrease ( Figures 8B , 7D ), therefore, N deficiency inhibited synthesis of guaiacyl units and syringyl units in the roots of apple dwarfing rootstock. However, P deficiency and N+P deficiency had no obvious inhibitory effect on guaiacyl synthesis. Additionally, in contrast to N deficiency, P deficiency and N+P deficiency posted a positive effect on the formation of syringyl units during lignin synthesis. Study has shown that overexpression of Gh4CL7 promoted lignin biosynthesis in A. thaliana, increased lignin content in plants, and promoted root elongation, thus improved drought tolerance of plants (Sun et al., 2020a). Therefore, the up-regulated 4CL genes and up-regulated syringyl units may be one of the important reasons that P deficiency and N+P deficiency promoted root elongation of ‘M9-T337’.

In summary, N and/or P deficiency enhanced the formation of hemicellulose xyloglucans, the extension of xylan backbone, and the loosening of cell wall. N deficiency increased cellulose formation, but inhibited the biosynthesis of hemicellulose xylan side chains, pectin and lignin. P deficiency and N+P deficiency posted a positive effect on the formation of syringyl units during lignin synthesis. These changes lead to different root morphogenesis under N and/or P deficient conditions.

Expansins are often regulated by nutrient stress to subsequently affect plant growth and nutrient uptake (Lee and Kim, 2013; Marowa et al., 2016; Kong et al., 2019). Liu et al. (2019) found that HvEXPA1 participates in root cell elongation and influences Al content by regulating root cell wall loosening when exposed to Al stress. Wu et al. (2022) screened 40 differentially expressed expansin genes in roots of oilseed rape (Brassica napus) in response to boron deficiency. Similar results were also reported by Nezamivand-Chegini et al. (2022). Through analyzing the similarities in the transcriptional responses in the roots under N and/or P deficit conditions, we observed up-regulation of six, two, and five expansin genes under N deficiency, P deficiency, and N+P deficiency, respectively. The content of L-tryptophan showed an increasing trend under N deficiency ( Figure 7C ). This amino acid is not only a precursor to auxin synthesis, but is also involved in expansin protein synthesis (Dai et al., 2013; Marowa et al., 2016). Lee and Kim (2013) reported that overexpression of EXPA17 in A. thaliana increased the density LRs treated with auxin. Under N deficiency, the recruitment of Thr101 phosphorylated NRT1.1 into the plasma membrane facilitates auxin transport, therefore root abundance increases under the action of expansins ( Figures 1 , 2 ) (Bouguyon et al., 2015; Zhang et al., 2019). Ding et al. (2021) found that P deficiency induced the expression of multiple expansin genes in roots, related to increased IAA content caused by P deficiency, thereby promoting root lengthening, which agrees with our results. In soybean, the protein product of GmPTF1, the expression of which is induced by Pi deficiency, directly binds to the E-box motif in the promoter region of the cell wall relaxation gene GmEXPB2, influencing LR formation and elongation, thereby improving P uptake efficiency (Yang et al., 2021).

Expansins break the noncovalent bonds between cellulose microfibrils and associated matrix polysaccharides and induce extension of the cell wall (Lee and Kim, 2013; Marowa et al., 2016). A drop in pH may be one of the reasons for increasing wall extension and the activity of the wall-modifying agent expansin ( Figure S7 ) (Arsuffi and Braybrook, 2018). Phylogenetic analysis showed that different expansins from various species falling within the same clade have similar effects on plant growth and development (Marowa et al., 2016). It has been reported that overexpressing RhEXPA4 increased the abundance of LRs in transgenic A. thaliana plants (Lü et al., 2013). In this study, we observed similar phenotypes; overexpression of MdEXPA4, which is highly homologous to RhEXPA4, promoted A. thaliana LR formation and development under N+P sufficiency ( Figure S5A ). Moreover, compared to the WT, the tolerance to P deficiency and N+P deficiency was enhanced in transgenic A. thaliana plants. Furthermore, under N deficiency, the enhanced root system resulted in a clear benefit to the aboveground growth of the transgenic plants. Boron et al. (2015) found that overexpression of AtEXLA promoted root elongation. Moreover, a recent study showed that GmEXLB1 expression was induced in soybean under P deficiency, and overexpression of this gene improved P acquisition by regulating root elongation and architecture in A. thaliana (Kong et al., 2019), which provided a possible direction for research on the function of this gene in apple. In this study, we confirmed that overexpression of MdEXLB1, which is highly homologous to GmEXLB1 and AtEXLA2, altered the root architecture of transgenic S. lycopersicum by increasing the number and length of LRs, thereby enhancing P acquisition. Additionally, we found that MdEXLB1 was induced by N deficiency. Overexpression of MdEXLB1 also enhanced N acquisition and biomass accumulation, and thus promoted the tolerance of transgenic S. lycopersicum to N deficiency ( Figures 9 , 10 ). Further understanding of the molecular mechanisms of expansin-based responses to N and/or P deficiency may lead to insights for improving nutrient utilization efficiency in apple rootstock.