Screening of Key Indices and the Gene Transcriptional Regulation Analysis Related to Salt Tolerance in Salix matsudana Seedlings

: Pot experiments were performed to comparatively study the differences in 16 salt tolerance indices between the seedlings of six Salix matsudana clones under the stress of various concentrations of NaCl (0, 0.1%, 0.3%, 0.5%, and 0.7%), including the salt injury index, shoot fresh weight, root fresh weight, leaf water content, relative conductivity, malondialdehyde content, and antioxidant enzyme activity. The salt-tolerant clones and key indices of salt tolerance were selected. Transcriptome sequencing analysis was performed on the selected salt-tolerant and salt-sensitive clones under salt stress, and the links between the physiological indices of salt tolerance and gene expression were analyzed. Results: (1) Superoxide dismutase (SOD), peroxidase (POD), chlorophyll, and net photosynthetic rate were closely related to the salt tolerance of Salix matsudana at the seedling stage. The regression equation was constructed as follows: salt tolerance index (y) = 0.224x 10 + 0.216x 11 + 0.127x 12 + 0.191x 7 − 0.187 (x 10 = chlorophyll, x 11 = SOD, x 12 = POD, x 7 = net photosynthetic rate). (2) The number of differentially expressed genes between the seedlings of salt-tolerant and salt-sensitive clones varied with the time of exposure (0 h, 4 h, 12 h, and 24 h) to 200 mmol · L − 1 NaCl stress. The most differentially expressed genes in Sm172 were detected upon 24 h vs. 4 h of salt treatment, while the most in Sm6 were in the 24 h vs. 0 h comparison. Gene Ontology analysis and Kyoto Encyclopedia of Genes and Genomes analysis showed that several differentially expressed genes were involved in carotenoid biosynthesis and plant mitogen-activated protein kinase signaling pathways. The nine highly expressed transcription factor genes ( Sm172-f2p30-2392 , Sm172-f2p28-2386 , Sm6-f8p60-2372 , Sm6-f2p39-2263 , Sm6-f16p60-2374 , Sm6-f3p60-931 , Sm6-f2p60-1067 , Sm172-f3p54-1980 , and Sm172-f3p54-1980 ) were closely correlated with the four key indices of salt tolerance. These genes could become genetic resources for salt tolerance breeding of Salix matsudana . transduction pathways al-tered in Salix matsudana under salt stress. More differentially expressed genes were detected after 12 h and 24 h of stress treatment, and more metabolic pathways were involved at these times, such as amino acid biosynthesis, plant hormone signal transduction, amino acid and nucleotide sugar metabolism, cysteine and methionine metabolism, mitogen-activated protein kinase (MAPK) signaling pathway, and carbon metabolism. These ﬁndings indicate that these metabolic pathways might play an important role in the salt stress response of Salix matsudana


Introduction
The coastal saline-alkali land of the Yellow River Delta suffers frequent natural disasters, soil salinization, and great difficulty in afforestation, which restricts the growth of many tree species [1][2][3]. The selection of salt-tolerant afforestation tree species is a key measure to improving the outcome of afforestation. Willow (Salix matsudana) is a native arbor species naturally distributed in this area that has a certain tolerance to salt. Selection and breeding of improved willow varieties could greatly help with the construction of shelter forests in the Yellow River Delta.
Sciences, as well as three unapproved clones with good traits (preliminarily named Jinan 1, Jinan 2, and Binzhou 1) in the experimental forest of Jinan, Shandong Province. The widely popular Willow No. 172 was used as the control. We numbered Luliu 2, Luliu 6, Jinan 1, Jinan 2, Binzhou 1, and Willow No. 172 as Sm2, Sm6, SmA, SmB, SmC, and Sm172 for this paper, respectively. The branches used in this study were collected from the experimental forest of Xicang Village, Changqing District, Jinan City, Shandong Province.

Experimental Methods
(1) Pot experiment: In April 2017, potted seedlings were cultivated with sieved soil. The soil was taken from the Forestry Experimental Station of Shandong Agricultural University. Before loading, it was screened to remove impurities and homogenized. The soil is sandy loam with a pH of 7.05 ± 0.03, a total nitrogen content of 1.27 ± 0.14 g·kg −1 , the hydrolyzable nitrogen content of 94.36 ± 13.89 mg·kg −1 , the available phosphorus content of 32.01 ± 4.10 mg·kg −1 , the available potassium content of 56.72 ± 7.83 mg·kg −1 , and the organic matter content of 16.27 ± 2.39 g·kg −1 . The cuttings of Salix matsudana clones Sm6 and Sm172 were selected and planted into greenhouse flowerpots at the Forestry Experimental Station of Shandong Agricultural University. The upper diameter of the flowerpots was 30 cm, the lower diameter was 20 cm, and the height was 25 cm. The weight of the soil in each pot was 10 kg. After planting, each cutting protruded approximately 2 cm above the soil surface. Regular watering management was conducted during the study. After the cuttage seedlings survived, the weak seedlings were eliminated in July 2017, while the seedlings that grew well and had no obvious diseases and pests were selected for experiments. The greenhouse was equipped with a water curtain, fan, and other temperature control equipment. During the experiment, the day and night temperatures in the greenhouse were controlled at (25 ± 2) • C and (20 ± 2) • C respectively, and the relative humidity was 65-70%.
Soil salinity was determined by the gravimetric method. A soil salinity gradient (0.1%, 0.3%, 0.5%, and 0.7%) was prepared, and deionized water was used as the blank control. First, the mass of NaCl needed for each stress level was calculated. NaCl solutions with the corresponding concentrations were prepared in deionized water and added three times in equal amounts (irrigation once every 7 days). A tray was placed below each pot to prevent the loss of salt. The physiological and biochemical indices were measured 15 days after the last addition of the salt.
(2) Transcriptome analysis: Transcriptome sequencing was performed using a combination of second-generation and third-generation sequencing technologies. Normally growing Sm6 and Sm172 adult plants (one each) were selected. Several branches with similar growth were cut at a length of approximately 15 cm. The cuttings were placed in beakers filled with deionized water and cultivated in a constant-temperature light incubator. After the branches took root and grew new leaves, various groups of branches were soaked in 200 mmol·L −1 NaCl solution for 0, 4, 12, and 24 h. The leaves were collected from all groups immediately after the soaking treatment. During sampling, leaves were collected from 3 plants in each group and mixed. Immediately after collection, the samples were placed in liquid nitrogen and then stored in a −80 • C freezer. The samples were saved for procedures such as RNA extraction.

Determination of Salt Tolerance Indices and Data Analysis
We grouped the indices related to the salt tolerance of Salix matsudana into four groups, salt injury indices, physiological and biochemical indices, biomass indices, and photosynthetic indices (Table 1), and used them to establish the standard equation of salt tolerance index. The root fresh weight (RFW), shoot fresh weight (SFW), relative leaf water content (RWC), relative seedling height (RH), and relative diameter (RD) were determined according to the method developed by Chen et al. [33]. The chlorophyll content was measured following Wang et al. [34]. Using a portable photosynthesis system (CIRAS-2), the net photosynthetic rate (Pn), stomatal conductance (GS), and intercellular CO 2 concentration (Ci) were measured between 9:00 a.m. and 11:00 a.m. at a saturated light intensity of 1200 µmol·m −2 s −1 and a leaf chamber temperature of 25 • C. The relative conductivity (RC) and malondialdehyde (MDA) content were determined as described by Qin et al. [35]. The activities of peroxidase (POD) and superoxide dismutase (SOD) were determined following Li et al., and the contents of soluble sugar (SS) and starch (SST) were determined according to Wang et al. [36,37]. All statistical analyses were conducted using Statistical Product Service Solutions for Windows 22.0 (SPSS, Chicago, IL, USA).

Sequencing Method and Analysis
Leaf RNA was extracted, and a cDNA library was constructed. Second-generation sequencing was performed on the Illumina HiSeq high-throughput sequencing platform. The third-generation full-length transcriptome sequencing was performed by Anoroad Gene Technology Co., Ltd. (Beijing, China) with the PacBio Sequel sequencer. Sequencing quality was assessed using FastQC. Functional annotation was performed with Trinotate. GO classification was performed with BLAST2GO. COG classification was performed using eggNOG. Metabolic pathway analysis was performed using the KEGG. Differential expression analysis was conducted with DESeq. The GO classification of the differentially expressed Unigene and Pathway enrichment analysis were performed using GOseq.

Verification of the Differentially Expressed Genes
RNAs were extracted from Sm6 and Sm172 that had received salt treatment for 0, 4, 12, or 24 h. The RNAs were reverse-transcribed, and 12 randomly selected salt-tolerancerelated genes of Salix matsudana were subjected to real-time fluorescence quantitative verification. The polymerase chain reaction (PCR) system (20 µL) was set up as follows: SYBR Green 10 µL, forward primer (10 mM) 1 µL, reverse primer (10 mM) 1 µL, cDNA 2 µL, and ddH 2 O to a total volume of 20 µL. The PCR conditions were as follows: predenaturation at 95 • C for 30 s and 40 cycles of denaturation at 95 • C for 5 s and optimal annealing temperature 60 • C for 30 s.  (Table 2) showed that all indexes showed significant or extremely significant differences under different NaCl concentrations and among different clones. Analysis of variance showed that the data of the indicators selected in this test were reliable. Correlation analysis of the 16 indices (Table 3) showed that there were various degrees of correlation between the indices in Salix matsudana exposed to salt stress. The salt injury index was positively correlated with the intercellular concentrations of CO 2 and malondialdehyde and negatively correlated with stomatal conductance, net photosynthetic rate, chlorophyll, biomass, and growth weight. The range of change varied between the indices, indicating that the results might be one-sided if single indices were used to evaluate the salt tolerance of Salix matsudana. The salt tolerance of Salix matsudana is a complex comprehensive trait.

Grey Relational Cluster Analysis of the Indices of Various Clones
The grey correlation degree was calculated from the membership function values of the relative values of the 16 indices of the six varieties (Table 4), and the indices were subjected to cluster analysis. At a fixed critical value r ∈ (0,1), the index xi and the index x j were considered to be the same type of index when ξ ij ≥ r (i = j). In this experiment, the critical value r = 0.57 was used to screen out two types of index clusters: C 1 = {x 4 , x 8 , x 11 , x 12 , x 14 , x 15 , x 16 } and C 2 = {x 1 , x 2 , x 3 , x 5 , x 6 , x 7 , x 9 , x 10 , x 13 }. In the above classification, leaf water content, intercellular CO 2 concentration, SOD, POD, electrical conductivity, soluble sugar, and soluble starch were classified into one category, while the salt injury index, chlorophyll content, shoot fresh weight, root fresh weight, seedling height growth, ground diameter growth, stomatal conductance, net photosynthetic rate, and malondialdehyde were classified into another category. Table 5 shows that these two types of indices contributed different amounts to the comprehensive evaluation of salt tolerance.  Since the cumulative variance contribution rates of the two types of indices were both greater than 85% (85.94% and 87.15%, respectively), these indices could be subjected to principal component analysis (Table 6). To evaluate the salt tolerance of the examined varieties, the comprehensive salt tolerance value (S) of each variety was calculated from the scores of the two types of indices obtained by principal component analysis (S = F 1 × W 1 + F 2 × W 2 ). The scoring formula for the C 1 comprehensive indices was F 1 = hc 1 × yc 1 , while the scoring formula for C 2 comprehensive indices was F 2 = hc 2 × yc 2 . The weight coefficient of principal components (h) was generally expressed as the variance contribution rate (namely, hc 1 = 0.8594, hc 2 = 0.8715). The linear combination of principal components of the C 1 and C 2 indices was obtained through principal component analysis: yc 1 = 0.113x 3 + 0.116x 10 + 0.124x 5 + 0.123x 1 + 0.118x 9 + 0.121x 6 + 0.123x 2 + 0.128x 7 + 0.112x 13 ; yc 2 = −0.15x 8 − 0.21x 4 + 0.391x 11 + 0.013x 14 − 0.262x 16 − 0.26x 15 + 0.344x 12 . The weight coefficients of the index clusters (w 1 = 0.561 and w 2 = 0.439) were calculated through the analytic hierarchy process. Notes: F 1 : Score of class C 1 index; F 2 : Score of class C 2 index; S: Consolidated values.
The S values are given in Table 6. On the basis of a cut-off value of 0.7, Sm6, Sm2, and SmA of the six Salix matsudana varieties were salt-tolerant, whereas Sm172, SmC, and SmB were salt-sensitive. Using the comprehensive salt tolerance value of the six varieties and the salt-stressrelated indices, a regression equation was established to identify the salt tolerance of Salix matsudana at the seedling stage. The equation was used to screen out the identification indices of salt tolerance at the seedling stage. An optimal regression equation was established through stepwise regression analysis using the comprehensive S value as the dependent variable and the relative values of the indices as the independent variables: y = 0.224x 10 + 0.216x 11 + 0.127x 12 + 0.191x 7 − 0.187 (r = 0.949, R 2 = 0.90), where y is the S index, x 10 is chlorophyll, x 11 is SOD, x 12 is POD, and x 7 is Pn. The regressors x 10 , x 11 , and x 12 were extremely significantly correlated with y (p = 0.009, 0.000, and 0.002, respectively), while x 7 was significantly correlated with y (p = 0.034). According to the equation, chlorophyll, SOD, POD, and Pn could be used as key indices to identify the salt tolerance of Salix matsudana at the seedling stage.

Classification of the Salix matsudana Genes That Were Differentially Expressed in Response to Salt Stress
The second-generation sequencing data were compared and quantified using the full-length transcripts as a reference. The obtained read count was compared between groups, yielding intergroup ratios. The genes with |log 2 Ratio| ≥ 1 and q < 0.05 were considered significantly differentially expressed genes between two groups. The number of differential genes varied with the duration of salt treatment ( Table 7). The differential genes were roughly divided into three major categories on the basis of the clustering results: 12 h salt treatment, 24 h salt treatment, and 0-4 h salt treatment ( Figure 1). The results of NT, NR, and BLASTX annotation showed that a greater number of differential genes were annotated at 12 h vs. 0 h, 12 h vs. 4 h, 24 h vs. 0 h, and 24 h vs. 4 h (4362 genes on average). According to the annotation results of the major databases, more than 98.33% of genes in the NT and NR databases were annotated (Table 8).

Pathway Analysis of the Differentially Expressed Genes
All the metabolic pathways in the KEGG were subjected to enrichment analysis using the hypergeometric test. The 20 KEGG metabolic pathways with the most significant enrichment of the differentially expressed genes were identified (Table 9), which gave us an intuitive understanding of the metabolic processes and signal transduction pathways altered in Salix matsudana under salt stress. More differentially expressed genes were detected after 12 h and 24 h of stress treatment, and more metabolic pathways were involved at these times, such as amino acid biosynthesis, plant hormone signal transduction, amino acid and nucleotide sugar metabolism, cysteine and methionine metabolism, mitogen-activated protein kinase (MAPK) signaling pathway, and carbon metabolism. These findings indicate that these metabolic pathways might play an important role in the salt stress response of Salix matsudana.

Salix matsudana Genes Differentially Expressed in Response to Salt Stress and Expression Verification
Analysis of the transcriptome data revealed that Salix matsudana adapted to and resisted salt stress and conducted self-regulation mainly through participating in carotenoid biosynthesis, MAPK signaling, phytohormone signal transduction, flavonoid biosynthesis, starch metabolism, and sucrose metabolism. The differential genes were screened out using a multiple of difference >2 and a corrected p < 0.01. Among the nine differential genes with the most significant gene expression difference, Sm172-f2p30-2392, Sm172-f2p28-2386, Sm6-f8p60-2372, and Sm172-f2p39-2386 were involved in the synthesis of 9-cis-epoxycarotenoid dioxygenase (NCED) and the regulation of NCED and abscisic acid (ABA) contents. Sm6-f2p39-2263 participated in the synthesis of asparagine synthase and positively regulated the response of Salix matsudana to salt stress. Sm6-f16p60-2374 was involved in circadian rhythm and photosensitive processes. Sm172-f3p54-1980 was involved in ABA-activated signaling pathways and was closely related to leaf senescence. Sm6-f3p60-931, Sm6-f2p60-1067, and Sm172-f2p30-1863 were involved in MAPK metabolic pathway. Sm6-f2p60-1067 was related to the PYL gene. It promoted the ABA signal response by inducing PYL overexpression. Sm172-f2p30-1863 was related to the synthesis of protein phosphatase 2C (PP2C). Silencing of Sm172-f2p30-1863 promoted the signal response of ABA. As a result of these changes, Salix matsudana became more sensitive to ABA and more promptly regulated its own response mechanism in the face of stress.

Discussion
This study analyzed the chosen indices using a combination of the analytic hierarchy process and principal component analysis, which was more reasonable than directly doing principal component analysis on all indices. SOD, POD, chlorophyll, and net photosynthetic rate may be used as key indices to identify the salt tolerance of Salix matsudana at the seedling stage [38][39][40]. Upon salt stress, the antioxidant enzymes such as SOD and POD were activated in all clones of Salix matsudana, which led to the elimination of the toxicity of free radicals. SOD prevents the formation of free radicals and the damage to cells mainly by regulating the content of O 2− . It is the first line of defense against cell damage [39,41]. POD plays an important role in protecting cells from H 2 O 2 stress [42]. In response to salt stress, photosynthesis is reduced in Salix matsudana. The accumulation of photosynthetic products is inhibited, which hinders the growth and development of plants and even causes plant death [43]. As an important photosynthetic pigment, chlorophyll plays an important role in absorption, transmission, and transformation of light energy [3,44,45]. Therefore, the change in chlorophyll content is an important index of plant growth [46]. Salt stress leads to ion toxicity, hinders the absorption of Mg 2+ by plants, blocks the synthesis of chlorophyll, reduces chlorophyll content, and thus seriously hinders photosynthesis [47][48][49][50][51]. Such findings are consistent with the research results of Gong et al. [52,53].
Under salt stress, reactive oxygen species activate a variety of MAPK signaling cascades, inducing specific cellular responses. This study found that in Salix matsudana exposed to salt stress, the most strongly expressed MAPK signaling pathway was the mitogen-activated protein kinase (MPK). This study showed that H 2 O 2 phosphorylated mitogen-activated protein kinase kinase 4/5 (MKK4/5) by activating ANP1 upstream of the MPK cascade, thereby phosphorylating MPK3. This result is consistent with the findings of Kim et al. [54] and Kovtun et al. [55]. H 2 O 2 also directly activated MPK3 by inducing the formation of nucleoside diphosphate kinase 2 (NDPK2), in line with the findings of Moon et al. [56]. The serine-threonine protein kinase OXI1 (oxidative-signal-inducible 1) was also activated by peroxide and responded through the downstream MPK3. In the process of cell death caused by reactive oxygen species through the above three pathways, MPK3 is significantly upregulated, and the level of reactive oxygen species (such as H 2 O 2 ) is increased [57], which eventually leads to cell death. When Salix matsudana is exposed to low salt stress or short-term salt stress, the upregulation of MPK3 and the increase in H 2 O 2 lead to elevated activities of SOD and POD in the antioxidant system, which allows the resistance to the damage of reactive oxygen species and maintenance of normal physiological function of cells. Therefore, improving the activity of antioxidant enzymes in plants and increasing the level of plant antioxidant metabolism are important ways to enhance the salt tolerance of plants. Alleviation of the harm of salt stress to Salix matsudana could be achieved by inhibiting the expression of the MPK3-related gene Sm6-f3p60-931, which suppressed the activity of MPK3 and inhibited the production of reactive oxygen species. This biological process played a certain positive role in the salt stress resistance of Salix matsudana and provided evidence that SOD and POD are key indices of salt tolerance in Salix matsudana.
In addition, this study showed that MPK3 of the MAPK pathway cascade was significantly upregulated under salt stress, which reduced the function of SPEECHLESS (SPCH) and inhibited stomatal development. The basic helix-loop-helix (bHLH)-type transcription factors that regulate stomatal development include three members: SPCH, MUTE, and FAMA. These proteins are involved in the regulation of the initiation of the M cell lineage, the formation of guard mother cells, and the symmetrical division of guard mother cells [58,59]. SPCH has a region between the bHLH domain and the C-terminal region that is rich in the acidic amino acids serine and threonine. It mediates the direct phosphorylation of SPCH by MPK3 [60][61][62], which hinders stomatal development and directly hinders photosynthesis. Therefore, physiological study of Salix matsudana showed that, with more concentrated or longer salt treatment, the net photosynthetic rate was reduced in each Salix matsudana clone. This finding provides evidence that the net photosynthetic rate is a key index of salt tolerance in Salix matsudana.
Moreover, genes with large differences in expression levels are involved in the synthesis of NCED. According to the metabolic pathway map of carotenoid biosynthesis, carotenoid biosynthesis is related to ABA synthesis [54,63]. ABA plays a central role in plant responses to stressful environments [64,65]. In the process of ABA synthesis, zeaxanthin epoxidase catalyzes the conversion of carotene to violaxanthin. NCED further catalyzes the conversion of 9-cis-flavin to xanthotoxin, which further leads to the production of abscisic aldehyde. Abscisic aldehyde is eventually converted into ABA under the action of an amine oxidase. These findings are consistent with the results of Li et al. [23]. NCED is a member of the carotenoid cleavage dioxygenase family. It is the most important rate-limiting enzyme in the process of ABA biosynthesis in higher plants and belongs to a class of key genes involved in the regulation of ABA biosynthesis [66,67]. NCED expression level is positively correlated with ABA content in plants [68][69][70]. Therefore, the goal of plant physiological self-regulation might be achieved through the control of the expression of NCED synthesis-related genes such as Sm172-f2p28-2386, Sm172-f2p30-2392, Sm172-f2p39-2386, and Sm6-f8p60-2372, which further control the expression of NCED and regulate the content of ABA in plants. In future research on the salt tolerance of Salix matsudana, ABA might become a key index.
These data provide a scientific basis for future research on the mechanism of Salix matsudana salt tolerance, deliver a simpler method for the identification and breeding of salt-tolerant Salix matsudana varieties, and can guide the selection of salt-tolerant tree species in the coastal saline-alkali land of the Yellow River Delta.

Conflicts of Interest:
The authors declare no conflict of interest.