Transcriptome and Metabolome Analysis Reveals Salt-Tolerance Pathways in the Leaves and Roots of ZM-4 (Malus zumi) in the Early Stages of Salt Stress

The breeding of salt-tolerant rootstock relies heavily on the availability of salt-tolerant Malus germplasm resources. The first step in developing salt-tolerant resources is to learn their molecular and metabolic underpinnings. Hydroponic seedlings of both ZM-4 (salt-tolerant resource) and M9T337 (salt-sensitive rootstock) were treated with a solution of 75 mM salinity. ZM-4’s fresh weight increased, then decreased, and then increased again after being treated with NaCl, whereas M9T337′s fresh weight continued to decrease. The results of transcriptome and metabolome after 0 h (CK) and 24 h of NaCl treatment showed that the leaves of ZM-4 had a higher content of flavonoids (phloretinm, naringenin-7-O-glucoside, kaempferol-3-O-galactoside, epiafzelechin, etc.) and the genes (CHI, CYP, FLS, LAR, and ANR) related to the flavonoid synthesis pathway showed up-regulation, suggesting a high antioxidant capacity. In addition to the high polyphenol content (L-phenylalanine, 5-O-p-coumaroyl quinic acid) and the high related gene expression (4CLL9 and SAT), the roots of ZM-4 exhibited a high osmotic adjustment ability. Under normal growing conditions, the roots of ZM-4 contained a higher content of some amino acids (L-proline, tran-4-hydroxy-L-prolin, L-glutamine, etc.) and sugars (D−fructose 6−phosphate, D−glucose 6−phosphate, etc.), and the genes (GLT1, BAM7, INV1, etc.) related to these two pathways were highly expressed. Furthermore, some amino acids (S-(methyl) glutathione, N-methyl-trans-4-hydroxy-L-proline, etc.) and sugars (D-sucrose, maltotriose, etc.) increased and genes (ALD1, BCAT1, AMY1.1, etc.) related to the pathways showed up-regulation under salt stress. This research provided theoretical support for the application of breeding salt-tolerant rootstocks by elucidating the molecular and metabolic mechanisms of salt tolerance during the early stages of salt treatment for ZM-4.


Introduction
"The fruit trees go up the hill and down the beach, do not compete for land with grain and cotton" has become a policy that must be implemented over the long term. The utilization of the "four wastelands" (barren hills, waste valleys, barren hillocks, and desolated beaches) is an important strategic measure for the development of new orchards. This could improve the efficiency of land utilization while simultaneously safeguarding the

Comparison of Salt Tolerance between ZM-4 and M9T337 in Response to Salt Stress
The fresh weights of the hydroponic seedlings of ZM-4 and M9T337 were obtained at 0 h, 24 h, 48 h, 72 h, 96 h, 120 h, and 144 h after the NaCl treatments. The fresh weight of ZM-4 increased from 0 h to 24 h after the NaCl treatments, decreased until 96 h, and then increased again from 96 h to 144 h (13.40 g, 14.26 g, 13.41 g, 13.11 g, 12.87 g, 13.01 g, and 13.05 g at 0 h, 24 h, 48 h, 72 h, 96 h, 120 h, and 144 h, respectively). However, after being exposed to the NaCl solution, M9T337's fresh weight decreased steadily from 0 h to 144 h (14.88 g, 114.46 g, 13.53 g, 12.90 g, 12.37 g, 12.41 g, and 11.88 g at 0 h, 24 h, 48 h, 3 of 28 72 h, 96 h, 120 h, and 144 h, respectively). The fresh weight of ZM-4 was heavier than that of M9T337 after 48 h of NaCl treatments, despite the fact that the initial fresh weight was lower (Figure 1a).

Comparison of Salt Tolerance between ZM-4 and M9T337 in Response to Salt Stress
The fresh weights of the hydroponic seedlings of ZM-4 and M9T337 were obtained at 0 h, 24 h, 48 h, 72 h, 96 h, 120 h, and 144 h after the NaCl treatments. The fresh weight of ZM-4 increased from 0 h to 24 h after the NaCl treatments, decreased until 96 h, and then increased again from 96 h to 144 h (13.40 g, 14.26 g, 13.41 g, 13.11 g, 12.87 g, 13.01 g, and 13.05 g at 0 h, 24 h, 48 h, 72 h, 96 h, 120 h, and 144 h, respectively). However, after being exposed to the NaCl solution, M9T337's fresh weight decreased steadily from 0 h to 144 h (14.88 g, 114.46 g, 13.53 g, 12.90 g, 12.37 g, 12.41 g, and 11.88 g at 0 h, 24 h, 48 h, 72 h, 96 h, 120 h, and 144 h, respectively). The fresh weight of ZM-4 was heavier than that of M9T337 after 48 h of NaCl treatments, despite the fact that the initial fresh weight was lower (Figure 1a).   (Table S1, Table S2). The Pearson correlation coefficient was used to evaluate sample repeatability, and the results revealed that the repeatability among the biological replicates was reliable ( Figure S1).

Comparison of Transcriptional Profiling between ZM-4 and M9T337 in Response to Salt Stress
Genes were identified as novel if they were not already known from comparison to the reference genome. Overall, we were able to pinpoint 41,575 unique genes across all of our samples. The average numbers of genes found in M-CKL, M-CKR, M-TL, M-TR, Z-CKL, Z-CKR, Z-TL, and Z-TR were 33,692, 35,569, 33,633, 35,687, 33,901, 35,716, 33,640, and 35,441, respectively (Table S3). A total of 2307 novel genes were identified from the samples, 126 of which were annotated in the KEGG pathway and 240, 305, and 314 of which were annotated in GO in terms of the cellular components, molecular functions, and biological processes, respectively (Table S4).   (Tables S1 and S2). The Pearson correlation coefficient was used to evaluate sample repeatability, and the results revealed that the repeatability among the biological replicates was reliable ( Figure S1).

Comparison of Transcriptional Profiling between ZM-4 and M9T337 in Response to Salt Stress
Genes were identified as novel if they were not already known from comparison to the reference genome. Overall, we were able to pinpoint 41,575 unique genes across all of our samples. The average numbers of genes found in M-CKL, M-CKR, M-TL, M-TR,  Z-CKL, Z-CKR, Z-TL, and Z-TR were 33,692, 35,569, 33,633, 35,687, 33,901, 35,716, 33,640,  and 35,441, respectively (Table S3). A total of 2307 novel genes were identified from the samples, 126 of which were annotated in the KEGG pathway and 240, 305, and 314 of which were annotated in GO in terms of the cellular components, molecular functions, and biological processes, respectively (Table S4).

DEG Analysis in Leaves in Response to Salt Stress for ZM-4 and M9T337
DEGs were used to identify the key genes involved in salt tolerance. There were 1217, 862, 4711, and 4550 DEGs for M-CKL vs. M-TL, Z-CKL vs. Z-TL, M-CKL vs. Z-CKL, and M-TL vs. Z-TL, respectively (Figure 1b, Table S5). The number of DEGs between ZM-4 and M9T337 with the same treatments was significantly higher than that of the different treatments for the same samples.
In comparison with the DEGs at 0 h and 24 h after the NaCl treatment for the GO term in leaves, the top 20 GO terms for the DEGs in M-CKL vs. M-TL were biological processes and molecular functions, with five and fifteen GO terms for them, respectively. The top three terms with a high rich factor were the ionotropic glutamate receptor activity (GO: 0004970), extracellular ligand-gated ion channel activity (GO: 0005230), and glutamate receptor activity (GO: 0008066) (Figure 2a and Figure S2). While the top 20 GO terms for the DEGs in Z-CKL vs. Z-TL differed from those of M9T337, including biological processes, molecular functions, and cellular components, the top three terms with a high rich factor were the nucleosome (GO: 0000786), FANCM-MHF complex (GO: 0071821), and protein heterodimerization activity (GO: 0046982) (Figure 2b and Figure S3). The 20 top GO terms for the DEGs in M-TL vs. Z-TL were biological processes and molecular functions, with the ionotropic glutamate receptor activity (GO: 0004970), extracellular ligand-gated ion channel activity (GO: 0005230), and ADP binding (GO: 0043531) being the top three terms with a high rich factor (Figure 2c and Figure S4).
The top 20 DEG pathways in M-CKL vs. M-TL, according to the KEGG pathway analysis, were metabolism, environmental information processing, and organismal systems. The top three pathways with a high rich factor were the synthesis and degradation of ketone bodies (ko00072), monoterpenoid biosynthesis (ko00902), and galactose metabolism (ko00052) (Figure 3a and Figure S5). The top 20 pathways for the DEGs in Z-CKL vs. Z-TL were different from those of M9T337. These were metabolism, genetic information processing, cellular processes, and organismal systems. The top three pathways with a high rich factor were DNA replication (ko03030), linoleic acid metabolism (ko00591), and sesquiterpenoid and triterpenoid biosynthesis (ko00909) (Figure 3b and Figure S6). The top 20 pathways for the DEGs in M-TL vs. Z-TL were metabolism, cellular processes, and organismal systems. The top two pathways with a high rich factor in M-TL vs. Z-TL were aflatoxin biosynthesis (ko00254) and linoleic acid metabolism (ko00591) (Figure 3c and Figure S7).
A Venn diagram among M-TL vs. Z-TL, M-CKL vs. M-TL, Z-CKL vs. Z-TL, and M-CKL vs. Z-CKL was constructed. The results showed that 136 DEGs were coexpressed in the first three combinations, which could be a response to salt stress for both ZM-4 and M9T337. Of these 136 DEGs, two DEGs (MD02G1120200 (LOX1.5) and MD03G1138500 (XA21)) were enriched in the MAPK signaling pathway, and nine DEGs were enriched in the biosynthesis of secondary metabolites, such as starch and sucrose metabolism, brassinosteroid biosynthesis, phenylpropanoid biosynthesis, and so on ( Figure 4a, Table  S6). While 95 DEGs were identified in response to salt stress 24 h after the NaCl treatment of ZM-4, the same was not true for M9T337. Among the 95 DEGs, 12 DEGs were enriched in the biosynthesis of secondary metabolites. In addition to the above pathways, flavonoid biosynthesis (MD15G1024100 (DFR), MD11G1059500 (CYP), and MD10G1311100 (ANR)) and ascorbate and aldarate metabolism (MD15G1201000 (MIOX2)) are unique pathways that may be the main genes for the salt tolerance of the leaves of ZM-4.

DEG Analysis of Roots in Response to Salt Stress for ZM-4 and M9T337
A total of 1259, 565, 7035, and 6594 DEGs for M-CKR vs. M-TR, Z-CKR vs. Z-TR, M-CKR vs. Z-CKR, and M-TR vs. Z-TR were identified, respectively (Figure 1b, Table S5). However, the number of DEGs in roots was greater than that in leaves.
The DEGs at 0 h and 24 h after the NaCl treatment were compared for the GO term in roots. The top 20 GO terms for the DEGs in M-CKR vs. M-TR were biological processes, molecular functions, and cellular components, and the top three terms with a high rich factor were the pectin catabolic process (GO:0045490), pectin metabolic process (GO: 00045488), and cellular hormone metabolic process (GO: 0034754) (Figure 2d and Figure S8). While the top 20 GO terms for the DEGs in Z-CKR vs. Z-TR were biological processes and molecular functions, the top three GO terms for the DEGs in Z-CKR vs. Z-TR were the hyaluronan metabolic process (GO: 0030212), mucopolysaccharide metabolic process (GO: 1903510), and threonine-phosphate decarboxylase activity (GO: 0048472) (Figure 2e and Figure S9). The top 20 GO terms for the DEGs in M-TR vs. Z-TR were biological processes and molecular functions, and the top three GO terms for the DEGs in M-TR vs. Z-TR were ADP binding (GO: 0043531), peptide:proton symporter activity (GO: 0015333), and proton-dependent peptide secondary activity (GO: 0022897) (Figure 2f and Figure S10).  processing, cellular processes, and organismal systems. The top three pathways with a high rich factor were DNA replication (ko03030), linoleic acid metabolism (ko00591), and sesquiterpenoid and triterpenoid biosynthesis (ko00909) (Figures 3b and S6). The top 20 pathways for the DEGs in M-TL vs. Z-TL were metabolism, cellular processes, and organismal systems. The top two pathways with a high rich factor in M-TL vs. Z-TL were aflatoxin biosynthesis (ko00254) and linoleic acid metabolism (ko00591) (Figures 3c  and S7).  The KEGG pathway analysis showed that the top 20 pathways for the DEGs in M-CKR vs. M-TR were metabolism, environmental information processing, and organismal systems. The top two pathways with a high rich factor were flavone and flavonol biosynthesis (ko00944) and zeatin biosynthesis (ko00908) (Figure 3d and Figure S11). The top 20 pathways for the DEGs in Z-CKR vs. Z-TR were different from those of M9T337. These were metabolism and environmental information processing. The top three pathways with a high rich factor were glucosinolate biosynthesis (ko00966), cyanoamino acid metabolism (ko00460), and glycosaminoglycan degradation (ko00531) (Figure 3e and Figure S12). The top 20 GO terms for the DEGs in M-TR vs. Z-TR were metabolism, environmental information processing, and organismal systems, and the top three pathways with a high rich factor in M-TR vs. Z-TR were zeatin biosynthesis (ko00908), diterpenoid biosynthesis (ko00904), and carotenoid biosynthesis (ko00906) (Figure 3f and Figure S13).
( Figure 4a, Table S6). While 95 DEGs were identified in response to salt stress 24 h after the NaCl treatment of ZM-4, the same was not true for M9T337. Among the 95 DEGs, 12 DEGs were enriched in the biosynthesis of secondary metabolites. In addition to the above pathways, flavonoid biosynthesis (MD15G1024100 (DFR), MD11G1059500 (CYP), and MD10G1311100 (ANR)) and ascorbate and aldarate metabolism (MD15G1201000(MIOX2)) are unique pathways that may be the main genes for the salt tolerance of the leaves of ZM-4.  Table  S5). However, the number of DEGs in roots was greater than that in leaves.
The DEGs at 0 h and 24 h after the NaCl treatment were compared for the GO term in roots. The top 20 GO terms for the DEGs in M-CKR vs. M-TR were biological processes, molecular functions, and cellular components, and the top three terms with a high rich factor were the pectin catabolic process (GO:0045490), pectin metabolic process (GO: 00045488), and cellular hormone metabolic process (GO: 0034754) (Figures 2d and  S8). While the top 20 GO terms for the DEGs in Z-CKR vs. Z-TR were biological processes and molecular functions, the top three GO terms for the DEGs in Z-CKR vs. Z-TR were the hyaluronan metabolic process (GO: 0030212), mucopolysaccharide metabolic process (GO: 1903510), and threonine-phosphate decarboxylase activity (GO: 0048472) (Figures 2e and S9). The top 20 GO terms for the DEGs in M-TR vs. Z-TR were biological processes and molecular functions, and the top three GO terms for the DEGs in M-TR vs. Z-TR were ADP binding (GO: 0043531), peptide:proton symporter activity (GO: 0015333), and proton-dependent peptide secondary activity (GO: 0022897) (Figures 2f and S10).
The KEGG pathway analysis showed that the top 20 pathways for the DEGs in M-CKR vs. M-TR were metabolism, environmental information processing, and organismal systems. The top two pathways with a high rich factor were flavone and flavonol biosynthesis (ko00944) and zeatin biosynthesis (ko00908) (Figures 3d and S11). The top The Venn diagram results for roots showed that 63 DEGs were coexpressed in M-TR vs. Z-TR, M-CKR vs. M-TR, and Z-CKR vs. Z-TR, which may be a response to salt stress for both ZM-4 and M9T337. Five DEGs were enriched in the biosynthesis of secondary metabolites, including phenylpropanoid biosynthesis, cysteine and methionine metabolism, and tyrosine metabolism. One DEGs (MD15G1007000 (P4H3)) was found to be enriched in arginine and proline metabolism. While 113 DEGs were identified as a response to salt stress 24 h after the NaCl treatment for ZM-4, this was not true for M9T337. Two DEGs (MD01G1158500 (PYL4) and MD07G1227100 (PYL4)) were enriched in the MAPK signaling pathway, one was enriched in plant hormone signal transduction, and seven DEGs were enriched in the biosynthesis of secondary metabolites, including zeatin biosynthesis, riboflavin metabolism, linoleic acid metabolism, and so on, which may be the main genes for the salt tolerance of the roots of ZM-4 ( Figure 4b, Table S6).

DAM Analysis in Roots in
The KEGG pathway analysis of M-CKR vs. M-TR showed that DAMs were mainly enriched in the metabolic pathways, biosynthesis of secondary metabolites, biosynthesis of antibiotics, biosynthesis of amino acids, and so on ( Figure S18). The top three VIP values of metabolites were L-arginine, chlorogenic acid, and cryptochlorogenic acid. LysoPC 18:3 and lysoPC 18:1 were induced to accumulate by salt stress in M9T337 (Figure 7d, Table S9). The KEGG pathway analysis of Z-CKR vs. Z-TR showed that DAMs were enriched in the metabolic pathways, biosynthesis of secondary metabolites, microbial metabolism in diverse environments, biosynthesis of antibiotics, and so on ( Figure  S19). The top three VIPs of metabolites were 3-hydroxyphloretin, L-norleucine, and L-leucine. L-cysteine and sanguisorbigenin were induced to accumulate by salt stress in ZM-4 ( Figure 7e, Table S9). The KEGG pathway analysis of M-TL vs. Z-TL showed that DAMs were enriched in the biosynthesis of secondary metabolites, biosynthesis of amino acids, flavonoid biosynthesis, flavone and flavonol biosynthesis, and so on ( Figure  S20). The top three VIPs of metabolites were 2α,3α,19α,23-tetrahydroxy-12-ursen-28-oic acid, quercetin-3-O-glucoside, and 1-O-feruloyl-β-d-glucose. Sanguisorbigenin only accumulated in Z-TR (Figure 7f, Table S9).

Integrated Analysis of the Transcriptome and Metabolome of ZM-4 and M9T337 Responsive to Salt Stress
A coexpression network analysis of the transcriptome and metabolome was performed to further exploit the relationship between DEGs and DAMs in the leaves and roots of ZM-4 and M9T337 under salt stress. The pathway function model and bidirectional orthogonal projections to latent structures model (O2PLS) were carried out to screen the associated genes and metabolites that had an influence on the sample grouping and analysis of the association characteristics. For the leaves, we used 227 DAMs and 6980 DEGs in the association analysis, and for the roots, we used 254 DAMs and 9436 DEGs (Table S11).

Integrated Analysis of the Transcriptome and Metabolome of Leaves Responsive to Salt Stress
The KEGG pathways shared by DAMs and DEGs were examined. A total of 40 pathways were obtained for M-CKL vs. M-TL, including 229 DEGs and 41 DAMs (Tables S12 and S13). The first three pathways were phenylpropanoid biosynthesis, the biosynthesis of secondary metabolites, and metabolic pathways; both the Gene_Pvalue and Gene_Qvalue were less than 0.05, but only the Metabolite_Pvalue of the biosynthesis of secondary metabolites was less than 0.05, making it a candidate pathway for further investigation. A total of 30 pathways were obtained for Z-CKL vs. Z-TL, including 151 DEGs and 41 DAMs. The first three pathways were the same as M-CKL vs. M-TL, while the pathway for which both the Gene_Pvalue and Metabolite_Pvalue were less than 0.05 was flavonoid biosynthesis. A total of 54 pathways were obtained for M-TL vs. Z-TL, including 793 DEGs and 54 DAMs. The first three pathways were the biosynthesis of secondary metabolites, flavonoid biosynthesis, and metabolic pathways. The Gene_Pvalue, Gene_Qvalue, Metabolite_Pvalue, and Metabolite_Qvalue of flavonoid biosynthesis were all less than 0.05, suggesting that this pathway was the major characteristic in the response to salt stress for ZM-4, which included 20 DEGs and 11 DAMs.

Integrated Analysis of the Transcriptome and Metabolome of Roots Responsive to Salt Stress
The KEGG pathways shared by DAMs and DEGs were examined. A total of 48 pathways were obtained for M-CKR vs. M-TR, including 250 DEGs and 46 DAMs (Tables S12 and S13). The first three pathways were metabolic pathways, the biosynthesis of secondary metabolites, and pentose and glucuronate interconversions, with both the Gene_Pvalue and Gene_Qvalue being less than 0.05, but only the Metabo-lite_Pvalue and Metabolite_Qvalue of the biosynthesis of secondary metabolites being less than 0.05, so it was the candidate pathway for further investigation. A total of 36 pathways were obtained for Z-CKR vs. Z-TR, including 124 DEGs and 38 DAMs. The first three pathways were the biosynthesis of secondary metabolites, metabolic pathways, and cyanoamino acid metabolism. The Gene_Pvalue, Gene_Qvalue, Metabolite_Pvalue, and Metabolite_Qvalue were less than 0.05. A total of 46 pathways were obtained for M-TR vs. Z-TR, including 1092 DEGs and 42 DAMs. The first three pathways were the biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, and metabolic pathways. The Gene_Pvalue, Gene_Qvalue, Metabolite_Pvalue, and Metabolite_Qvalue of the first two pathways were all less than 0.05, indicating that these pathways are the most important in terms of ZM-4's response to salt stress.    (Table S14). The top 25 DAMs and DEGs were thought to have a high relevance according to the sum of the squares of the first two dimension loading values (Table S15, Figure 9d-f). The first three DAMs and DEGs of M-CKR vs. M-TR were fhloretin-4'-O-(6"-feruloyl)glucoside, L-leucine, L-histidine, MD01G1146100 (At2g39920), MD17G1029300 (LBD25), and MD17G1043400 (IP5P5), respectively. The first three DAMs and DEGs of Z-CKR vs. Z-TR were LysoPE 18:2, L-asparagine, 2-hydroxyoleanolic acid, MSTRG.9301, MD00G1015600 (fmdA), and MD14G1153600 (rps3), respectively. The first three DAMs and DEGs of M-TR vs. Z-TR were 3,4-dihydroxybenzoic acid, kaempferol-3-Oarabinoside, 2-hydroxycinnamic acid, MD07G1256300 (IRE), MD15G1179600 (RPV1), and MD08G1150100, respectively.

The Salt Tolerance Mechanisms Predicted in Leaves and Roots for ZM-4
The DAMs based on the VIP values (>1) and FC (≥1 or ≤ 0.5) of the metabolites for M-CKL vs. M-TL, Z-CKL vs. Z-TL, and M-TL vs. Z-TL were mainly classified as flavonoids (Table S9). The flavonoid biosynthesis pathway (ko00941) was also enriched by a DEG and DAM joint analysis (Tables S12 and S13). The heatmaps of expression for the DEGs and DAMs of this pathway in M-CKL, M-TL, Z-CKL, and Z-TL were drawn out ( Figure 10). Except for quercetin glycosides, kaempferol-3-O-sophoroside, and phloretin-2 -O-glucoside, the other DAMs were upregulated in Z-TL, and the genes related to the DAMs (CHI, CYP, FLS, LAR, and ANR) also showed high expression levels in Z-TL. As a result, the higher expression of the genes and components in the phenylpropanoid biosynthesis pathway was the primary salt-tolerance mechanism in the leaves of ZM-4. In addition to the flavonoid pathway, the phenylpropanoid biosynthesis pathway was enriched in the roots of M-CKR vs. M-TR, Z-CKR vs. Z-TR, and M-TR vs. Z-TR (Table S9, Table S12, Table S13). 5-O-caffeoyl shikimic acid, chlorogenic acid, and trans-5-O-(p-coumaroyl) shikimate were found in high concentrations in M-TR, whereas 5-O-p-coumaroyl quinic acid and L-phenylalanine were found in high concentrations in In addition to the flavonoid pathway, the phenylpropanoid biosynthesis pathway was enriched in the roots of M-CKR vs. M-TR, Z-CKR vs. Z-TR, and M-TR vs. Z-TR (Tables S9, S12 and Table S13). 5-O-caffeoyl shikimic acid, chlorogenic acid, and trans-5-O-(p-coumaroyl) shikimate were found in high concentrations in M-TR, whereas 5-O-pcoumaroyl quinic acid and L-phenylalanine were found in high concentrations in T-CKR ( Figure S21). The expression of SAT (MD10G1151400), HST (MD17G1224900), and 4CLL9 (MD14G1161200) were induced by salt stress in the roots of ZM-4. The expression of SAT (MD16G1108700 and MD13G1109000) was induced by salt stress in the roots of M9T337.
Furthermore, ZM-4's fresh weight continued to rise even after 24 h following the NaCl treatment, indicating that osmotic regulation played a crucial role in the response to salt stress. As a result, the osmotic regulating substances, such as amino acids and sugars, were analyzed using a heatmap ( Figure 11). There were two paths for maintaining the osmotic pressure. First, it was found that ZM-4, untreated with NaCl, contained a relatively high concentration of various compounds, such as L-proline, tran-4-hydroxy-Lprolin, L-glutamine, L-asparagine, cis-aconitic acid, glutaric acid, argininosuccinic acid, D−fructose 6−phosphate, D−glucose 6−phosphate, D−sorbitol, and so on (  Table S5, Table S8).

Discussion
High Na + concentrations in the soil generate hyperosmotic conditions that can severely impact plant nutrients and water uptake, as well as cause leaf withering and, eventually, plant death [41,42]. When exposed to saline soil, the first stress experienced by a plant is osmotic stress, which impedes plant growth and development [11,[43][44][45].

Discussion
High Na + concentrations in the soil generate hyperosmotic conditions that can severely impact plant nutrients and water uptake, as well as cause leaf withering and, eventually, plant death [41,42]. When exposed to saline soil, the first stress experienced by a plant is osmotic stress, which impedes plant growth and development [11,[43][44][45]. According to the findings of this study, the fresh weight of M9T337 fell following the NaCl treatment, whereas the fresh weight of ZM-4 increased from 0 h to 24 h, decreased until 96 h, and then climbed again from 120 h to 144 h. Although lower at 0 h, the fresh weight of ZM-4 was higher than that of M9T337 after 72 h. This phenotype indicates that ZM-4 was less influenced by the osmotic stress from 0 h to 24 h, even in the later stages. This phenotype demonstrates that ZM-4 was less influenced by the osmotic stress after 24 h. ZM-4 was shown to be more salt tolerant because it both naturally contained and was induced to produce more osmotic-regulating substances in response to salt stress.
The mapped reads ranged from 51.07% to 70.95% among the eight libraries, yielding a total of 41,575 genes and resulting in the identification of 2307 novel genes. The total number of genes was less than that of the GDDH13 (reference genome) [46] and the HFTH1 [47]. The ratio of mapped reads was slightly lower and more novel genes were discovered, which could be attributed to ZM-4's and M9T337's genetic link being distantly related to golden delicious (GDDH13), the level of which was lower in the leaves than in the roots. This was primarily because leaves have fewer genes than roots. It was preferable to use the Malus zumi genome de novo assembly as the reference genome.
ROS are produced as a result of high salinity stress and are largely detoxified by the enzymes SOD, CAT, APX, and POD [44,48,49]. This study discovered that the genes of APX3 (MD15G1242900, MD15G1243000, MD00G1173600, and MD05G1120800), APXT (MD00G1029300 and MD00G1029400), SOD (MD01G1164700), SODA (MD07G1232200), and CAT1 (MD06G1008600) were found to be highly expressed in the leaves of ZM-4. It could be one of the reasons for the salt tolerance of ZM-4. Along with these enzymes, flavonoids also perform a major function in ROS scavenging [50][51][52][53]. DEGs and DAMs were found to be more enriched in the flavonoid pathway in ZM-4 leaves. CHI, CYP, FLS, LAR, and ANR were up-regulated in the leaves of ZM-4, and high accumulations of phloretin, naringenin chalcoe, kaempferol-3-O-galactoside, kaempferol-3-O-glucoside, epiafzelechin, and epicatechin were found. As a result, the flavonoid pathway may play a key role in the response of ZM-4 leaves to salt stress.
The principal locations of salt stress response are the roots, and the initial salt reaction is derived primarily from the roots. Roots are severely damaged by stress-induced ROS excess [54,55]. The POD, SOD, and CAT activities were higher in Malus zumi roots than in the leaves [46]. Polyphenols were also crucial in reducing ROS, alongside antioxidant enzymes and flavonoids [56][57][58]. In the roots of ZM-4, the content of L-phenylalanine and 5-O-p-coumaroyl quinic acid accumulated, and MD17G1224900 (HST), MD10G1161400 (SAT), and MD14G1161200 (4CLL9) were highly expressed in the roots of ZM-4. This suggests that ZM-4's salt tolerance is in part due to its roots' antioxidant capacity.
Not only can reactive oxygen species (ROS) damage root function, but so does osmotic dysregulation. Roots are not able to absorb water if subjected to osmotic stress [59]. In this study, the fresh weight of ZM-4 increased before 24 h following the NaCl treatment, while that of M9T337 kept falling afterward. This suggests that M9T337 was subjected to higher degrees of osmotic stress than ZM-4, implying that ZM-4's roots contained more osmoticregulating substances. Soluble sugars, palmitoleic acid, D-arginine, pheophytin a, rutin, and vanillin were found to be crucial in promoting salt-stressed plant growth [60][61][62]. Under normal conditions, ZM-4 contained a high concentration of L-proline, tran-4-hydroxy-Lprolin, L-glutamine, lactobiose, D-glucosamine, and other compounds. It can maintain a high osmotic potential in the roots and absorb nutrients and water properly in the early phases of salt stress; this may be the one of the main reasons for its tolerance to salt stress. Furthermore, some other amino acids and sugars were induced to accumulate after salt stress, such as S-(methyl) glutathione, N-methyl-trans-4-hydroxy-L-proline, D-sucrose, planteose, maltotriose, and so on, and the genes (PKP1, PFK2, GLT1, SUS2, SUS3, SS4, etc.) related to these substances were up-regulated. This increased the roots' ability to regulate the osmotic pressure. This demonstrates that plant hormone signal transduction was enriched in the leaves, while starch and sucrose metabolism, as well as proline biosynthesis, were enriched in trifoliate orange roots [10]. Sucrose's reactions as a signaling molecule were crucial to Malus halliana in maintaining an osmotic equilibrium and eliminating ROS during salt stress. It directly mediated the accumulation of D-phenylalanine, tryptophan, and alkaloids (vindoline and ecgonine), as well as the expression of aspartate and glutamaterelated proteins (ASP3, ASN1, NIT4, and GLN1−1) [24]. Malus robusta's tolerance to salt, alkali, and salt-alkali stress was conferred by SNP182G on MdRGLG3, which converted a leucine to an arginine at the vWFA domain [32]. This could be the principal salt-tolerance mechanism in roots for ZM-4, because it had a high content of some amino acids and sugars and a high expression of related genes on its own. Additionally, salt stress induced the accumulation and expression of other amino acids, sugars, and related genes.

Plant Materials
The experiment was conducted in an artificial climate chamber at 26-28 • C with a photoperiod of 16 h light/8 h dark. Tissue culture seedlings of ZM-4 and M9T337 were transplanted to point trays after rooting for a month in the tissue culture bottles, and then the seedlings were transplanted to hydroponic tanks with Hoagland nutrient solution (Ca(NO 3

Phenotype of the Samples
The fresh weight was measured at 0 h, 24 h, 48 h, 72 h, 96 h, 120 h, and 144 h after the 75 mM NaCl treatment. Three seedlings were used to collect the fresh weight for the hydroponic seedlings of ZM-4 and M9T337. The average fresh weight of three seedlings was used as the phenotype to evaluate the salt tolerance of ZM-4 and M9T337.

RNA Quantification and Qualification
The total RNA was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. The RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase free agarose gel electrophoresis. After the total RNA was extracted, the rRNAs were removed to retain mRNAs. The enriched mRNAs were fragmented into short fragments by using a fragmentation buffer and reverse-transcribed into cDNA with random primers. Second-strand cDNA was synthesized with DNA polymerase I, RNase H, dNTP (dUTP instead of dTTP), and buffer. Next, the cDNA fragments were purified with a QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), the ends were repaired, poly (A) was added, and the fragments were ligated to Illumina sequencing adapters. Then, UNG (uracil-N-glycosylase) was used to digest the second-strand cDNA. The digested products were size-selected by agarose gel electrophoresis, amplified with PCR, and sequenced using Illumina HiSeqTM 4000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Filtering of Clean Reads
The reads obtained from the sequencing machines included raw reads containing adapters or low-quality bases, which would affect the following assembly and analysis. Thus, to obtain high-quality clean reads, the reads were further filtered by fastp (version 0.18.0) [63]. The parameters were as follows: (1) reads containing adapters were removed; (2) reads containing more than 10% of unknown nucleotides (N) were removed; and (3) low-quality reads containing more than 50% low-quality (Q-value ≤ 20) bases were removed.
The short-reads alignment tool Bowtie2 (version 2.2.8) [64] was used for mapping the reads to the ribosome RNA (rRNA) database. The rRNA-mapped reads were then removed. The remaining reads were further used in the assembly and analysis of the transcriptome.

Novel Transcript Identification and Annotation
An index of the reference genome was built, and paired-end clean reads were mapped to the reference genome (GDDH13 v 1.1) [46] using HISAT2 (version 2.1.0) [65] with "-rnastrandness RF" and other parameters set as a default.
The reconstruction of transcripts was carried out with the software Stringtie (version 1.3.4) [66,67], which, together with HISAT2, allows biologists to identify new genes and new splice variants of known ones. To identify the new transcripts, all of the reconstructed transcripts were aligned to the reference genome and divided into twelve categories by using Cuffcompare. Transcripts with one of the classcodes "u, i, j, x, c, e, o" were defined as novel transcripts. We used the following parameters to identify reliable novel genes: a transcript length of longer than 200 bp and an exon number of more than 1. Novel transcripts were then aligned to the Nr, KEGG, and GO databases to obtain the protein functional annotation.

Quantification of Transcript Abundance
The transcript abundances were quantified by the software StringTie in a referencebased approach. For each transcription region, an FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated to quantify its expression abundance and variations using the StringTie software.
The FPKM formula is shown as follows: FPKM = 10 6 C NL/10 3 where FPKM(A) is the expression of transcript A, C is the number of fragments mapped to transcript A, N is the total number of fragments that were mapped to reference genes, and L is the number of bases on transcript A. The FPKM method was able to eliminate the influence of different transcript lengths and the amount of sequencing data on the calculation of transcript expression. Therefore, the calculated transcript expression can be directly used for comparing the difference in transcript expression among the samples.

Differentially Expressed Transcript (DEG) Analysis
The differentially expressed transcripts of coding RNAs were analyzed. An RNA differential expression analysis was performed by the DESeq2 [68] software between two different groups (and by edgeR [69] between two samples). The genes/transcripts with a false discovery rate (FDR) parameter below 0.05 and an absolute fold change ≥ 2 were considered differentially expressed genes/transcripts. Differentially expressed coding RNAs were then subjected to an enrichment analysis of GO functions and KEGG pathways.

Chemicals and Reagents
All chemicals and reagents were of analytical grade. Methyl alcohol, acetonitrile, and ethyl alcohol were purchased from Merck Company, Germany. Milli-Q system (Millipore Corp., Bedford, MA, USA) ultrapure water was used throughout the study. Authentic standards were purchased from BioBioPha Co., Ltd.(Kunming, China) and Sigma-Aldrich (St. Louis, MO, USA).

Sample Preparation and Extraction
The freeze-dried samples were crushed using a mixer mill (MM 400, Retsch) with a zirconia bead for 1.5 min at 30 Hz. Then, 100 mg of powder was weighed and extracted overnight at 4 • C with 1.0 mL of 70% aqueous methanol containing 0.1 mg/L lidocaine for an internal standard. Following centrifugation at 10,000× g for 10 min, the supernatant was absorbed and filtrated (SCAA-104, 0.22 µm pore size; ANPEL, Shanghai, China, www.anpel.com.cn/ (accessed on 9 May 2022)) before the LC-MS/MS analysis. Quality control (QC) samples were mixed by all samples to assess the reproducibility of the whole experiment. The LIT and triple quadrupole (QQQ) scans were acquired on a triple quadrupolelinear ion trap mass spectrometer (Q TRAP), AB Sciex QTRAP6500 system, equipped with an ESI-Turbo Ion-Spray interface, operating in a positive ion mode and controlled by the Analyst 1.6.1 software (AB Sciex). The operation parameters were as follows: ESI source temperature, 500 • C; ion spray voltage (IS), 5500 V; curtain gas (CUR), 25psi; and the collision-activated dissociation (CAD) was set to the highest. The QQQ scans were acquired as MRM experiments with an optimized declustering potential (DP) and collision energy (CE) for each individual MRM transition. The m/z range was set between 50 and 1000.

Data Pre-Processing and Metabolite Identification
The data filtering, peak detection, alignment, and calculations were performed using the Analyst 1.6.1 software. To produce a matrix containing fewer biased and redundant data, peaks were checked manually for a signal/noise (s/n) > 10 and in-house software written in Perl was used to remove the redundant signals caused by different isotopes; in-source fragmentation; K + , Na + , and NH4 + adducts; and dimerization. To facilitate the identification/annotation of metabolites, an accurate m/z for each Q1 was obtained. The total ion chromatograms (TICs) and an extracted ion chromatogram (EICs or XICs) for the QC samples were exported to give an overview of the metabolite profiles of all samples. The area of each chromatographic peak was calculated. The peaks were aligned across the different samples based on the spectral pattern and retention time. The metabolites were identified by searching an internal database and public databases (MassBank, KNApSAcK, HMDB [71], MoTo DB, and METLIN [72]) and comparing the m/z values, the RT, and the fragmentation patterns with the standards.

Multivariate Statistical Analysis
For a preliminary visualization of the differences between different groups of samples, an unsupervised dimensionality reduction method principal component analysis (PCA) was applied to all samples using R package models.
The OPLS-DA model was further validated by cross-validation and 200 permutation tests [73]. For cross-validation, the data were partitioned into seven subsets, where each of the subsets was then used as a validation set. R 2 indicated the total variation in the data matrix that was explained by the model. The predictive ability (Q 2 ) values represented the most recognized diagnostic statistical parameter to validate the OPLS-DA model in metabolomics. An acceptable predictive model was considered for a Q 2 value greater than 0.4. A good predictive model was considered for a Q 2 value greater than 0.9. A permutation test randomly permuted class labels 200 times and then produced a distribution of R 2 values and Q 2 values. In essence, a reliable model should yield significantly larger R 2 and Q 2 values compared to the R 2 and Q 2 values generated from random models using the same dataset.
The loadings from (O) PLS were the directions of projection with respect to the original variables. Variables whose loadings were far away from the origin in a loadings plot might be inferred to have the greatest contribution to class separation.

Differential Metabolite Analysis
A variable importance in projection (VIP) score of the (O) PLS model was applied to rank the metabolites that best distinguished between the two groups. The threshold of the VIP was set to 1. In addition, a t-test was also used as a univariate analysis for screening differential metabolites. Those with a p-value of the t-test < 0.05 and a VIP ≥ 1 were considered differential metabolites between two groups.

KEGG Pathway Analysis
KEGG is the major public pathway-related database that includes not only genes, but also metabolites [74]. The metabolites were mapped to KEGG metabolic pathways for a pathway analysis and an enrichment analysis. A pathway enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways in differential metabolites by comparing with the whole background. The calculating formula is as follows: Here, N is the number of all metabolites with KEGG annotation, n is the number of differential metabolites in N, M is the number of all metabolites annotated to specific pathways, and m is number of differential metabolites in M. The calculated p-value was obtained through an FDR correction, taking FDR ≤ 0.05 as a threshold. Pathways meeting this condition were defined as significantly enriched pathways in differential metabolites. KEGG pathway maps are the linking of genomic or transcriptomic contents of genes to chemical structures of endogenous molecules, thus providing a method of performing an integration analysis of genes and metabolites. All differentially expressed genes and metabolites in this study were mapped to the KEGG pathway database to obtain their links in metabolic pathways.

O2PLS Model
In order to integrate the transcriptomic and metabolomic data, we performed a twoway orthogonal PLS (O2PLS) analysis [75]. This method decomposes the variation present in the two data matrices into three parts: the joint variation between the two datasets, the orthogonal variation that is unique to each dataset, and noise. The model assumes that some latent variables are responsible for the variation in the joint and orthogonal parts. O2PLS models were calculated using the OmicsPLS package of R. To determine the optimal number of components, the proposed alternative cross-validation procedure was utilized [76]. The best models were used for the integration analysis.

Conclusions
For this experiment, we tested the hydroponic seedlings of ZM-4 and M9T337 by treating them with a 75 mM NaCl solution for 144 h to see how they handled the salt stress. ZM-4 was found to be more tolerant to salt stress when tested using fresh weight. The tolerance of leaves and roots differed significantly. The tolerance of the leaves of ZM-4 was primarily due to their higher content of flavonoids, and they had a strong antioxidant capacity. Not only did the ZM-4 roots have a greater anti-oxidant capacity, but they also contained more osmotic regulators such as amino acids and sugars (L-proline, tran-4hydroxy-L-prolin, D−fructose 6−phosphate, D−glucose 6−phosphate, D−sorbitol, and so on) under normal conditions, and even more amino acids and sugars (N-methyl-trans-4-hydroxy-L-proline, 2-hydroxyethylphosphonic acid, D-sucrose, maltotriose, melezitose, and so on) induced by salt stress. Furthermore, the genes (GLT1, ALDC, AA5GT, BAM7, CM3, METK4, PKP1, AMY1.1, etc.) that responded in both directions were highly expressed. This research could serve as a theoretical foundation for the use of salt-tolerant ZM-4 in the breeding of salt-tolerant apple rootstocks.

Data Availability Statement:
The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (Genomics, Proteomics & Bioinformatics 2021) [77] in the National Genomics Data Center (Nucleic Acids Res 2022) [78], China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA009194, accessed on 10 December 2022), which is publicly accessible at https://ngdc.cncb.ac.cn/gsa.