Transcriptome Analysis of Populus euphratica under Salt Treatment and PeERF1 Gene Enhances Salt Tolerance in Transgenic Populus alba × Populus glandulosa

Populus euphratica is mainly distributed in desert environments with dry and hot climate in summer and cold in winter. Compared with other poplars, P. euphratica is more resistant to salt stress. It is critical to investigate the transcriptome and molecular basis of salt tolerance in order to uncover stress-related genes. In this study, salt-tolerant treatment of P. euphratica resulted in an increase in osmo-regulatory substances and recovery of antioxidant enzymes. To improve the mining efficiency of candidate genes, the analysis combining both the transcriptome WGCNA and the former GWAS results was selected, and a range of key regulatory factors with salt resistance were found. The PeERF1 gene was highly connected in the turquoise modules with significant differences in salt stress traits, and the expression levels were significantly different in each treatment. For further functional verification of PeERF1, we obtained stable overexpression and dominant suppression transgenic lines by transforming into Populus alba × Populus glandulosa. The growth and physiological characteristics of the PeERF1 overexpressed plants were better than that of the wild type under salt stress. Transcriptome analysis of leaves of transgenic lines and WT revealed that highly enriched GO terms in DEGs were associated with stress responses, including abiotic stimuli responses, chemical responses, and oxidative stress responses. The result is helpful for in-depth analysis of the salt tolerance mechanism of poplar. This work provides important genes for poplar breeding with salt tolerance.


Introduction
Soil salinity is considered to be a principal problem for international agricultural manufacturing and sustainability [1]. Salt stress now severely limits the growth, development, yield, and quality of plants worldwide [2]. Salt stress is projected to affect more than 6% of the world's land surface and over 20% of irrigated land, restricting the productivity and geographic dispersion of farmed crops.
For adaptation to saline environments, plants have evolved intricate methods for sensing environmental cues and coordinating metabolic processes and morphological characteristics through gene expression regulation [3]. Certain physiological structures of plants undergo adaptive changes to cope with unfavorable survival conditions, and such physiological changes are commonly seen in saline plants, where the outer skin cell walls thicken and leave flesh out to conserve water and accumulate salt [4]. Significant research has also been conducted around the changes in gene expression that occur in response to salt stress [2,5]. Salt tolerance genes in plants are associated with many pathways, GWAS results, and found that some genes, such as POPULUS_EUPHRATICA_07301 (ERF1), POPULUS_EUPHRATICA_13633 (NAC21/22), POPULUS_EUPHRATICA_28205 (ITN1), etc., are involved in the osmotic stress response of plants. The PeERF1 gene was selected to transform with poplar for functional validation. RNA sequencing of PeERF1 transgenic P. alba × P. glandulosa screened many stress-related genes and found that overexpressed transgenic P. alba × P. glandulosa functioned mainly after salt stress. Our results can help to screen excellent resistance gene resources in poplar, and combined with genetic engineering, it is important to breed resistant and high-quality tree species suitable for our needs.

Validation of Salt Stress Treatment in P. euphratica
To explore the salt tolerance mechanism of P. euphratica, Plants at four-time points (0, 12, 24, and 48 h) and two treatment concentrations (A, 150 mM NaCl; B, 300 mM NaCl) under NaCl were employed in physiological research. Before the harvest season, different treatments were started at the same time, and all samples were collected at the same time [35]. To evaluate the physiological responses to salt stress, we assessed superoxide dismutase (SOD), peroxidase (POD), and malondialdehyde (MDA).
Peroxidase is an important endogenous scavenger of reactive oxygen species, which is closely related to plant stress resistance. POD can reduce plant damage and remove excessive reactive oxygen species. Under 150 mM and 300 mM treatments, there was a significant upward trend from 0 h to 12 h and a significant decrease from 12 h to 48 h ( Figure 1). Interestingly, SOD increased first and then decreased under salt stress. MDA decreased first and then increased at 150 mM and 300 mM. RNA sequencing at four timepoints with two salt treatments. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed on DEG to identify the different metabolic pathways involved in the salt stress response. We com pared the overlapped genes between WGCNA modules and former GWAS results, and found that some genes, such as POPULUS_EUPHRATICA_07301 (ERF1), POPULUS_EU PHRATICA_13633 (NAC21/22), POPULUS_EUPHRATICA_28205 (ITN1), etc., are in volved in the osmotic stress response of plants. The PeERF1 gene was selected to trans form with poplar for functional validation. RNA sequencing of PeERF1 transgenic P. alb × P. glandulosa screened many stress-related genes and found that overexpressed trans genic P. alba × P. glandulosa functioned mainly after salt stress. Our results can help to screen excellent resistance gene resources in poplar, and combined with genetic engineer ing, it is important to breed resistant and high-quality tree species suitable for our needs

Validation of Salt Stress Treatment in P. euphratica
To explore the salt tolerance mechanism of P. euphratica, Plants at four-time points (0 12, 24, and 48 h) and two treatment concentrations (A, 150 mM NaCl; B, 300 mM NaCl under NaCl were employed in physiological research. Before the harvest season, differen treatments were started at the same time, and all samples were collected at the same tim [35]. To evaluate the physiological responses to salt stress, we assessed superoxide dis mutase (SOD), peroxidase (POD), and malondialdehyde (MDA).
Peroxidase is an important endogenous scavenger of reactive oxygen species, which is closely related to plant stress resistance. POD can reduce plant damage and remov excessive reactive oxygen species. Under 150 mM and 300 mM treatments, there was a significant upward trend from 0 h to 12 h and a significant decrease from 12 h to 48 h ( Figure 1). Interestingly, SOD increased first and then decreased under salt stress. MDA decreased first and then increased at 150 mM and 300 mM. , and malonaldehyde (MDA) contents between different times and concentrations in P. euphratica; the control is water. Three biological replicate and three technical replicates were used in the studies. The student's t-test was used to assess th data. p < 0.05 of three-time intervals (12,24, and 48 h) versus the control (0 h) and two treatmen concentrations (150 mM and 300 mM) versus the control (0 mM), respectively. Abbreviations: PeuC P. euphratica control; A, 150 mM NaCl salt-stressed; B, 300 mM NaCl salt-stressed; 12, salt stress fo 12 h; 24, salt stress for 24 h; 48, salt stress for 48 h (* p-value < 0.05).

Analysis of RNA Sequences and Identification of Differentially Expressed Genes in P. euphratica under Salt Stress
To further elucidate the molecular basis of salt response in P. euphratica, RNA se quencing (RNA-seq) was used for transcriptome analysis. The cDNA libraries were cre ated separately from the leaves taken at two concentrations (150 mM and 300 mM) and , and malonaldehyde (MDA) contents between t different times and concentrations in P. euphratica; the control is water. Three biological replicates and three technical replicates were used in the studies. The student's t-test was used to assess the data. p < 0.05 of three-time intervals (12,24, and 48 h) versus the control (0 h) and two treatment concentrations (150 mM and 300 mM) versus the control (0 mM), respectively. Abbreviations: PeuC, P. euphratica control; A, 150 mM NaCl salt-stressed; B, 300 mM NaCl salt-stressed; 12, salt stress for 12 h; 24, salt stress for 24 h; 48, salt stress for 48 h (* p-value < 0.05).

Analysis of RNA Sequences and Identification of Differentially Expressed Genes in P. euphratica under Salt Stress
To further elucidate the molecular basis of salt response in P. euphratica, RNA sequencing (RNA-seq) was used for transcriptome analysis. The cDNA libraries were created separately from the leaves taken at two concentrations (150 mM and 300 mM) and fourtime points (0 h, 12 h, 24 h, and 48 h) following salt stress, with three duplicates at each time point, and then sequenced on the Illumina Hiseq 2500 platform. After the initial sequencing readings have been checked for quality, 915,771,706 clean readings (137.37 Gb) were retained, and the paired-end was 125 bp for further analysis. The content of GC was 43.96%, and that of Q30 was 93.7% (Supplementary Table S1).
The FPKM values were determined based on the uniquely mapped reads for both control and salt stress conditions to determine the relative expression levels of genes. Using comparison analysis, the data of RNA-seq revealed a total of 18,002 differently expressed genes in reaction to salt stress. Out of the total DEGs, 8624 (47.9%) of those DEGs were up-regulated, while 9378 (52.1%) of them were down-regulated (Supplementary Table S2). A larger amount of variation in gene expression in response to salt stress are indicated by a broader dispersion.
DEGs were discovered in several comparisons and grouped into four groups: Group I for various treatment timepoints under 150 mM salt stress compare to control (Peu12A vs. PeuC, Peu24A vs. PeuC, and Peu48A vs. PeuC) (Figure 2A), Group II for different treatment timepoints under 300 mM salt stress compare to control (Peu12B vs. PeuC, Peu24B vs. PeuC, and Peu48B vs. PeuC) ( Figure 2B), Group III for comparison with 150 mM salt stress at different times (Peu24A vs. Peu12A and Peu48A vs. Peu24A) ( Figure 2C) and Group IV for comparison with 300 mM salt stress at different times (Peu24B vs. Peu12B and Peu48B vs. Peu24B) ( Figure 2D). In the comparison of "Peu12B vs. PeuC," (a total of 3771 DEGs including 1932 up-regulated genes and 1839 down-regulated genes) the biggest DEG set was discovered, which indicated that the young leaves of P. euphratica rapidly reprogrammed the cellular response at the transcriptome level at salt stress for 12 h. The lowest DEG set was discovered in the comparison "Peu48B vs. Peu24B" (588 DEGs including 324 up-regulated genes and 264 down-regulated genes), followed by comparison "Peu48A vs. Peu24A" (762 DEGs including 221 up-regulated genes and 541 down-regulated genes) ( Figure 2E). In three comparisons, the overlapped DEGs of Group II (823 common genes, 4.57% of total DEGs) were bigger than those of the other groups, according to comparative analysis (Figure 2A,D). The most extensive subset of particular DEGs was "Peu12B vs. Peu12A" (6086 DEGs in Group III) and "Peu12A vs. PeuC" (3771 DEGs in Group II). Furthermore, in the fourth group, there were few overlapping of differential genes ( Figure 2D).

GO and KEGG Enrichment Analysis of DEGs
The identified DEGs were annotated with 18 biological processes (BP), 9 cellular components (CC), and 9 molecular functions (MF) in GO categories (Supplementary Figure S1), and significantly enriched (p ≤ 0.001) into 122 GO terms (Supplementary Table S3). The terms "cellular process", "metabolic process", "single-organism process", "biological regulation", "localization", "regulation of biological process" and "response to stimulus", were the dominant groups in the biological processes; "cell", "cell part", "membrane", "membrane-part", "organelle", and "macr-molecular complex" were of the representative groups in the cellular components. Among the molecular functions, a great number of DEGs were focused on categories of "binding", "catalytic activity", and "transporter activity".
To elucidate the molecular pathways driving salt tolerance in P. euphratica, the DEGs were characterized using the GO knowledgebase (http://geneontology.org/ (accessed on 20 December 2020)). The enhanced functions for all DEGs owing to salt stress are shown in GO enrichment scatterplots (Supplementary Figure S2 and Supplementary Table S4). DEGs were enriched including "cellular carbohydrate metabolic process" and "photosynthesis" that were prominent in Group I comparisons (Peu12A vs. PeuC, Peu24A vs. PeuC and Peu48A vs. PeuC), but did not enrich in Group II comparisons (Peu12B vs. PeuC, Peu24B vs. PeuC and Peu48B vs. PeuC). "Cell wall organization or biogenesis" and "aminoglycan catabolic process" were enriched in Group III comparisons (Peu24A vs. Peu12A and Peu48A vs. Peu24A), but did not enrich in Group IV (Peu24B vs. Peu12B and Peu48B vs. Peu24B). GO terms under the MF category include "carboxylic ester hydrolase activity" was enriched in Group I, "heme binding" and "tetrapyrrole binding" were enriched in Group II. In the CC category, DEGs were significantly enriched in GO terms of "extracellular region", "external encapsulating structure", "cell wall", "cell periphery", and "apoplast" in all groups.   To further understand the molecular response to salt of P. euphratica at the pathway level, we used KEGG pathways method to perform enrichment analysis. In this approach, a total of 20 routes were greatly enhanced (Supplementary Figure S3 and Supplementary  Table S5). In particular, "Plant hormone signal transduction", "Galactose metabolism", "MAPK signaling pathway-plant", and "Cutin, suberine and wax biosynthesis" were significantly enriched in most comparisons. The pathway of "Phenylpropanoid biosynthesis" was highly enriched in the early and middle phases (Peu24A vs. PeuC and Peu24A vs. Peu12A); the pathway of "Cysteine" and "Amino sugar and nucleotide sugar metabolism" was significantly enriched in medium and late stages (Peu48B vs. PeuC and Peu48A vs. Peu24A). Moreover, the pathways of "Starch and sucrose metabolism", "Nitrogen metabolism' and "Flavonoid biosynthesis" were specifically enriched in different stages in Group II and Group III comparisons.

Identification of WGCNA Modules Related to Phenotypic Traits
We used WGCNA to construct a co-expression network to study the interrelationships among salt responsive genes. The DEGs state that their expression patterns of expression were split into seventeen distinct modules. Seventeen separate modules were recognized and labeled with various colors in a dendrogram ( Figure 3A). These results indicated that different degrees of salt stress development were under phenotypic traits control. The turquoise module of 9438 DEGs, green of 763 DEGs and blue of 4624 DEGs were significantly and stably correlated with phenotypic traits including POD, SOD, and MDA (r > 0.8 or r < −0.8). The correlation network of the turquoise, green and blue module genes with WGCNA edge weight >0.3 were shown in Figure 3B-D. After analyzing the overlapped genes between WGCNA modules and prior GWAS results [11], we found 23, 7, and 2 candidate genes for GWAS intervals in the turquoise, blue and green modules, respectively (Supplementary Table S6). Analysis of major transcription factor (TF) families presented in different modules revealed that they were differentially sensitive to salt response effects (Supplementary Figure S6). Overall, TFs were strongly represented in turquoise module (217 TFs, 46.5% of DEGs), which including 55 MYB, 26 AP2, 22 b-ZIP, 18 B3, 17 WRKY, 15 ARF, 15 GRAS, 11 GATA, 10 TCP, 7 SBP, 7 FAR1, 6 HSF, 3 E2F/DP, 2 SAP, 2 ZF-HD, and 1 HD-ZIP (Supplementary Table S7). We observed that some genes (ERF1, NAC21/22, MYB48, WRKY6 etc.) are also involved in plants' osmotic stress response. The ethylene response factor (ERF) family is one of the largest plant-specific transcription factor families, playing an important role. Among them, the PeERF1 gene was highly connected in the turquoise module, and the traits were significantly different, and the expression levels of each treatment stage were significantly different. We selected PeERF1 for functional validation of transgenic poplar.
For quantitative Real-Time PCR (qRT-PCR) validation, eight genes were chosen from the DEG list (Supplementary Table S8

Cloning of PeERF1 and Protein Sequence Comparisons
To verify the reliability of WGCNA, we selected PeERF1 to conduct functional verification. P. euphratica provided the cDNA fragment for PeERF1. It comprises a 663 bp coding sequence that encodes 221 amino acids with a molecular weight of 53.814 k Da (Supplementary Figure S5).

PeERF1 Subcellular Location and Transcriptional Activity Analysis
To explore the distribution of PeERF1 transcription factor expression products in cells, we constructed expression vectors containing transcription factors using the GFP

Cloning of PeERF1 and Protein Sequence Comparisons
To verify the reliability of WGCNA, we selected PeERF1 to conduct functional verification. P. euphratica provided the cDNA fragment for PeERF1. It comprises a 663 bp coding sequence that encodes 221 amino acids with a molecular weight of 53.814 k Da (Supplementary Figure S5).

PeERF1 Subcellular Location and Transcriptional Activity Analysis
To explore the distribution of PeERF1 transcription factor expression products in cells, we constructed expression vectors containing transcription factors using the GFP fusion protein reporter, and transiently introduced the successfully constructed expression vectors into Nicotiana tabacum leaves [36]. As shown, the fluorescent signal of the 35S: PeERF1-YFP protein was exclusively observed in the nucleus. These results show that PeERF1 is a nuclear gene localized (Supplementary Figure S7A,B).
To investigate PeERF1 transcriptional activity, we created the pGBKT7-PeERF1 fusion vector and converted into Y2HGold yeast cells. Negative control yeast cells were Pgadt7 co-transformed yeast cells. Positive control yeast cells were pGADT7-T co-transformed yeast cells. As shown in Figure 5A, it has self-activating activity in yeast.

Promoter Activity and Relevant Transient Acting Elements Analysis of PeERF1
To evaluate the tissue-specific expression of the PeERF1 gene, various DNA segments of 1524 bp from the promoter of the PeERF1 promoter were inserted into poplar genome to induce the production of the β-glucuronidase (GUS). GUS histochemical staining revealed that GUS was expressed mostly in the top young leaves of poplar, but not in the stems or roots (Supplementary Figure S7I).
In order to investigate the role of PeERF1 gene and elements, the GCC box, DRE motif, TTG1 motif, and TTG2 motif, which can bind to ERF subfamily proteins, were selected with reference to the relevant literature [32,[37][38][39]. To determine whether PeERF1 can also bind to these elements, we performed a yeast one-hybrid assay. The constructed pHIS2

Promoter Activity and Relevant Transient Acting Elements Analysis of PeERF1
To evaluate the tissue-specific expression of the PeERF1 gene, various DNA segments of 1524 bp from the promoter of the PeERF1 promoter were inserted into poplar genome to induce the production of the β-glucuronidase (GUS). GUS histochemical staining revealed that GUS was expressed mostly in the top young leaves of poplar, but not in the stems or roots (Supplementary Figure S7).
In order to investigate the role of PeERF1 gene and elements, the GCC box, DRE motif, TTG1 motif, and TTG2 motif, which can bind to ERF subfamily proteins, were selected with reference to the relevant literature [32,[37][38][39]. To determine whether PeERF1 can also bind to these elements, we performed a yeast one-hybrid assay. The constructed pHIS2 reporter vector containing the GCC-box, DRE motif, TTG1 motif, and TTG2 motif and the pGADT7-Rec2 effector vector carrying PeERF1 and the pGADT7-Rec2 effector vector carrying PeERF1 were co-transformed into Y187 yeast cells. The yeast cells of the two co-transformed reporter and effector vectors were able to grow on DDO and TDO/3-AT (75 mM) medium, demonstrating that PeERF1 can bind to the GCC-box, DRE motif, TTG1 motif, and TTG2 motif, just as the positive control ( Figure 5B).

Morphological Changes of Transgenic P. alba × P. glandulosa under Salt Stress
We selected the high transformation efficiency transgenic P. alba × P. glandulosa healing culture method for transformation to obtain transgenic plants [40]. Four transgenic P. alba × P. glandulosa lines and the wild type (WT) were given different treatments of 0, 50, 75, and 100 mM NaCl. Under normal circumstances, the plant height, root length, and fresh weight of 35S: PeERF1 transgenic poplar were 1.06, 1.05, and 1.08 times more than those of WT, respectively. 35S: SRDX-PeERF1 transgenic poplar was 0.91, 0.99, and 0.91 times more than WT. Additionally, plant height, root length, and fresh weight of 35S: PeERF1 transgenic poplars were 1.36, 1.

Morphological Changes of Transgenic P. alba × P. glandulosa under Salt Stress
We selected the high transformation efficiency transgenic P. alba × P. glandulosa healing culture method for transformation to obtain transgenic plants [40].  To compare transgenic and WT under salt stress, one-month-old plants were treated with 150 mM salt concentration for 7 days. Plant height and fresh weight were determined. Transgenic poplars and WT grew normally in the control environment without obvious symptoms. However, under salt stress, WT started to wilt slightly and leaves fell off then. The growth state of 35S:PeERF1 (OE) was better than that of WT, while the To compare transgenic and WT under salt stress, one-month-old plants were treated with 150 mM salt concentration for 7 days. Plant height and fresh weight were determined. Transgenic poplars and WT grew normally in the control environment without obvious symptoms. However, under salt stress, WT started to wilt slightly and leaves fell off then. The growth state of 35S:PeERF1 (OE) was better than that of WT, while the growth state of 35S: SRDX-PeERF1 (SE) transgenic lines was worse than that of WT ( Figure 7A). The plant height and fresh weight of all OE plants were significantly higher than those of WT, while SE showed the opposite trend ( Figure 7B,C).

Analysis of Physiological Characteristics of PeERF1 Transgenic P. alba × P. glandulosa
To investigate the role of PeERF1 in osmotic stress, we measured physiological indicators of transgenic lines and WT under salt treatment. The accumulation of ROS in PeERF1 transgenic is lower than in WT after the treatments. After salt treatment, the transgenic plants' peroxidase (POD), soluble sugar content, and superoxide dismutase (SOD) levels were greater than the WT plants ( Figure 8A,B,D), although the transgenic plants' malondialdehyde (MDA) levels were lower, indicating that scavenging activity was increased ( Figure 8C). All of these findings suggest that PeERF1 overexpression improves salt tolerance.

Analysis of Physiological Characteristics of PeERF1 Transgenic P. alba × P. glandulosa
To investigate the role of PeERF1 in osmotic stress, we measured physiological indicators of transgenic lines and WT under salt treatment. The accumulation of ROS in PeERF1 transgenic is lower than in WT after the treatments. After salt treatment, the transgenic plants' peroxidase (POD), soluble sugar content, and superoxide dismutase (SOD) levels were greater than the WT plants ( Figure 8A,B,D), although the transgenic plants' malondialdehyde (MDA) levels were lower, indicating that scavenging activity was increased ( Figure 8C). All of these findings suggest that PeERF1 overexpression improves salt tolerance.
The electrolyte leakage rate (LR) and leaf chlorophyll content (CC) were measured to assess the plants' growth stage. The LR and CC of transgenic and WT lines were not substantially different under normal growing conditions ( Figure 8E,F). Under salt stress, the OE lines grew faster than the WT, with considerably greater LR and CC than the WT. The SE lines grew in the opposite direction, with much lower LR and CC than the WT.
The contents of hydrogen peroxide and superoxide were detected using Nitrotetrazolium blue chloride (NBT) and 3,3 -diaminobenzidine (DAB), respectively. Between transgenic and control lines, there was no obvious difference in NBT staining. After salinization, the staining area of 35S: PeERF1 transgenic was substantially less than that of WT after salinization. DAB staining has a similar trend (Figure 9). The staining area of transgenic 35S: SRDX-PeERF1 that suppressed expression was substantially bigger than that of WT. DAB staining has a similar trend.
PeERF1 transgenic is lower than in WT after the treatments. After salt treatment, the transgenic plants' peroxidase (POD), soluble sugar content, and superoxide dismutase (SOD) levels were greater than the WT plants ( Figure 8A,B,D), although the transgenic plants' malondialdehyde (MDA) levels were lower, indicating that scavenging activity was increased ( Figure 8C). All of these findings suggest that PeERF1 overexpression improves salt tolerance.  The electrolyte leakage rate (LR) and leaf chlorophyll content (CC) were measured to assess the plants' growth stage. The LR and CC of transgenic and WT lines were not substantially different under normal growing conditions ( Figure 8E,F). Under salt stress, the OE lines grew faster than the WT, with considerably greater LR and CC than the WT. The SE lines grew in the opposite direction, with much lower LR and CC than the WT.
The contents of hydrogen peroxide and superoxide were detected using Nitrotetrazolium blue chloride (NBT) and 3,3′-diaminobenzidine (DAB), respectively. Between transgenic and control lines, there was no obvious difference in NBT staining. After salinization, the staining area of 35S: PeERF1 transgenic was substantially less than that of WT after salinization. DAB staining has a similar trend (Figure 9). The staining area of transgenic 35S: SRDX-PeERF1 that suppressed expression was substantially bigger than that of WT. DAB staining has a similar trend.

Expression Analysis of Downstream Genes Regulated of PeERF1
To further investigate the regulation of those downstream genes to improve salt tolerance in plants. Based on the results of the analysis of variance, we identified genes as substantially different with FDR (p-value) < 0.05 and |log2fc| > 1. 126 differentially expressed genes were identified in the control from transgenic P. alba × P. glandulosa overexpressing the PeERF1 gene, of which 109 genes were up-regulated and 17 genes were down-regulated ( Figure 10). A total of 132 differentially expressed genes were identified in transgenic poplars, of which 96 genes were up-regulated and 36 genes were downregulated. In the treatment group, 935 differentially expressed genes were identified from transgenic poplars overexpressing PeERF1, of which 207 genes were up-regulated and 728 genes were down-regulated. A total of 44 differentially expressed genes were identified in transgenic poplars, of which 38 genes were up-regulated and 6 genes were down-regulated.

Expression Analysis of Downstream Genes Regulated of PeERF1
To further investigate the regulation of those downstream genes to improve salt tolerance in plants. Based on the results of the analysis of variance, we identified genes as substantially different with FDR (p-value) < 0.05 and |log2fc| > 1. 126 differentially expressed genes were identified in the control from transgenic P. alba × P. glandulosa overexpressing the PeERF1 gene, of which 109 genes were up-regulated and 17 genes were down-regulated ( Figure 10). A total of 132 differentially expressed genes were identified in transgenic poplars, of which 96 genes were up-regulated and 36 genes were down-regulated. In the treatment group, 935 differentially expressed genes were identified from transgenic poplars overexpressing PeERF1, of which 207 genes were up-regulated and 728 genes were down-regulated. A total of 44 differentially expressed genes were identified in transgenic poplars, of which 38 genes were up-regulated and 6 genes were down-regulated. To further investigate the function of these genes, overexpressed and WT differential genes under salt stress were classified into three categories: biological processes, cellular components and molecular functions. These genes are associated with response to acidic chemicals, seed dormancy process, dormancy process, oxidoreductase activity, negative regulation of the developmental process, response to abiotic stimulus, response to chemical, response to oxidative stress, response to oxygen-containing compound, and response to reactive oxygen species (Figure 11). Overexpression of transgenic PeERF1 was compared between salt stress treatments and controls using a q-value of ≤0.05. There were 88 significantly enriched pathways with To further investigate the function of these genes, overexpressed and WT differential genes under salt stress were classified into three categories: biological processes, cellular components and molecular functions. These genes are associated with response to acidic chemicals, seed dormancy process, dormancy process, oxidoreductase activity, negative regulation of the developmental process, response to abiotic stimulus, response to chemical, response to oxidative stress, response to oxygen-containing compound, and response to reactive oxygen species (Figure 11). To further investigate the function of these genes, overexpressed and WT differential genes under salt stress were classified into three categories: biological processes, cellular components and molecular functions. These genes are associated with response to acidic chemicals, seed dormancy process, dormancy process, oxidoreductase activity, negative regulation of the developmental process, response to abiotic stimulus, response to chemical, response to oxidative stress, response to oxygen-containing compound, and response to reactive oxygen species (Figure 11). Overexpression of transgenic PeERF1 was compared between salt stress treatments and controls using a q-value of ≤0.05. There were 88 significantly enriched pathways with Overexpression of transgenic PeERF1 was compared between salt stress treatments and controls using a q-value of ≤0.05. There were 88 significantly enriched pathways with 621 differentially significant genes involved in the salt stress treatment. The top 20 pathways by the number of differential genes are shown in Figure 11.
To verify the accuracy of the transcriptome sequencing results, we selected 16 stressrelated genes, designed primers in their conserved domains for RT-PCR experiments, and then compared them with the transcriptome results. The results presented were nearly identical to the transcriptome sequencing results (Figure 12). 621 differentially significant genes involved in the salt stress treatment. The top 20 pathways by the number of differential genes are shown in Figure 11.
To verify the accuracy of the transcriptome sequencing results, we selected 16 stressrelated genes, designed primers in their conserved domains for RT-PCR experiments, and then compared them with the transcriptome results. The results presented were nearly identical to the transcriptome sequencing results (Figure 12).

Discussion
The response signal pathway of abiotic stress is related to the physiological response of plants [41]. Treatment of P. euphratica with different concentrations of salt (50, 100, and 150 mM) reduced stomatal opening and leaf photosynthetic capacity, and the activities of POD and SOD were significantly higher than those of the control group This is similar to the results of physiological indicators of P. euphratica treated at different times in our study (Figure 1). Some stress-related genes of P. euphratica have also been studied in plants. For example, Compared with the control, OE-PeCBF4a exhibited higher SOD activity and significantly reduced MDA content, indicating that transgenic OE-PeCBF4a poplars were more tolerant to salt stress [42]. The activities of antioxidant enzymes such as SOD, POD, and CAT in tobacco plants overexpressing PeREM1.3 under salt stress significantly increased and decreased ROS levels, maintain ROS homeostasis, and improved the salt tolerance of plants [43]. In this study, the changes of POD, SOD, and MDA activities, soluble sugar content, chlorophyll content and relative electrical conductivity of transgenic poplar showed that the PeERF1 gene could significantly improve its salt tolerance ( Figure 9).

Discussion
The response signal pathway of abiotic stress is related to the physiological response of plants [41]. Treatment of P. euphratica with different concentrations of salt (50, 100, and 150 mM) reduced stomatal opening and leaf photosynthetic capacity, and the activities of POD and SOD were significantly higher than those of the control group. This is similar to the results of physiological indicators of P. euphratica treated at different times in our study (Figure 1). Some stress-related genes of P. euphratica have also been studied in plants. For example, Compared with the control, OE-PeCBF4a exhibited higher SOD activity and significantly reduced MDA content, indicating that transgenic OE-PeCBF4a poplars were more tolerant to salt stress [42]. The activities of antioxidant enzymes such as SOD, POD, and CAT in tobacco plants overexpressing PeREM1.3 under salt stress significantly increased and decreased ROS levels, maintain ROS homeostasis, and improved the salt tolerance of plants [43]. In this study, the changes of POD, SOD, and MDA activities, soluble sugar content, chlorophyll content and relative electrical conductivity of transgenic poplar showed that the PeERF1 gene could significantly improve its salt tolerance (Figure 9).
Based on the association analysis technology between genes and traits, WGCNA can reduce the dimensionality of complex data and simplify it into several modules. Weiss proposed Gene Module Association Study (GMAS) as a supplement to GWAS analysis results [44]. By integrating the results of GWAS and WGCNA analysis, Farber found that this method can significantly improve the mining efficiency of the micro-effect sites of the yellow seed phenotype, and found a shortcut for GWAS to solve its own limitations [45].
GWAS on the salt tolerance index of seed germination rate in four different branches of P. euphratica in Northwest China. A total of 38 SNP were found to be associated with seed salt tolerance, located within or near 82 genes [11]. Comparing the results with the GWAS, it was found that there were 41 overlapping genes in the 17 modules in WGCNA. qRT-PCR analysis of eight genes showed that the results were consistent with the transcriptome data ( Figure 4). For example, NAC21/22 responds to salt stress by indirectly regulating two metabolic pathways, phytohormone signaling and mitogen-activated protein kinase signaling [46]. The AtITN1 mutant negatively regulates salt stress by mediating the BOHC and RBOHD genes [47].
In order to verify the efficiency and reliability of WGCNA data screening for salt tolerance genes, PeERF1 gene was selected for salt tolerance analysis, and the positive response of PeERF1 gene to plant salt stress was demonstrated by functional verification. In this study, we proved that PeERF1 is a nuclear localization protein with transcriptional autoactivation activity. Overexpression of PeERF1 gene can enhance the salt tolerance of transgenic poplars, and dominant suppression of transgenic lines can inhibit the salt tolerance of transgenic plants. This result provides new progress and evidence for the involvement of PeERF1 gene in the study of plant salt tolerance.
To obtain more information about molecular pathways underlying PeERF1's response to salt stress, transcriptome sequencing of PeERF1 (OE), SRDX-PeERF1 (SE) strains and WT was carried out before and after salt stress treatments. In the analysis of GO enrichment of over differential genes after stress, it was found that most GOs were mainly associated with response to acid chemical, seed dormancy process, dormancy process, oxidoreductase activity, negative regulation of developmental process, response to abiotic stimulus, response to chemical, response to oxidative stress, response to oxygen-containing compound, response to reactive oxygen species and regulation of gibberellic acid mediated signaling pathway were found to be associated. The stress-related genes studied are as follows. PtrERF2 plays an important role in the growth and stress resistance of plants among the homologous genes of other plants [48]. Overexpression of miR160 could cause degradation of the ARF18 target gene transcript and increased salt tolerance in miR160 overexpress transgenic plants [49]. PtrWRKY70 enhanced tolerance to osmotic stress among its homologous genes in other plants [50]. PtrCAT1 belongs to the CAT enzyme family, which in other plant homologous genes can eliminate excess reactive oxygen species during stress resistance [51][52][53]. Thus, genes downstream of PeERF1 may respond to abiotic stresses through different biological processes or pathways, thereby enhancing plant stress tolerance.
Except for PeERF1, which was functionally validated in our investigation, the other genes may also play a role in plant salt tolerance. For example, compared with wild type plants, ZmPP2C-A1, ZmPP2C-A2, and ZmPP2C-A6 overexpression plants displayed greater germination rates after ABA and NaCl treatments than wild type plants in A. thaliana [54]. Although there have been no reports of CLH2 function in salt tolerance, it has been confirmed that AtPKL-1 mutant caused the up-regulation of chlorophyll degradation gene CLH2 gene under cold stress in A. thaliana [55]. More research is needed to determine whether these genes contribute to salt tolerance in P. euphratica.

Plant Growth and NaCl Treatment
All materials used in this study were seedlings because the phenotype of poplar from seedling to mature stage is more prominent than that of older trees, and the seedling stage can eliminate environmental differences. Freshly matured seeds were harvested from the P. euphratica sown, germinated in the culturing room, and the seedlings were moved to tissue culture bottles for clone propagation. Under the aseptic condition, the growth of the plant was the same on the rooting medium. After the fibrous roots grew, they were transplanted into a pot culture containing peat soil and vermiculite (the ratio of peat soil and vermiculite was 3:1). Under the condition of long sunshine, it grows in a growth chamber underneath 26-28 • C, light for 16 h, and dark for 8 h. Plants have been subjected to 150 mM and 300 mM NaCl to simulate four-time variables of salt stress (0 h, 12 h, 24 h, and 48 h). At the cease of every time point, young leaves have been harvested, straight away frozen in liquid nitrogen, and saved in a −80 • C freezer for similar evaluation and sequencing. We made three biological replicates. (For the convenience of marking the following P. euphratica as Peu, 150 mM as A, and 300 mM as B).

Measurement of the Enzyme Activity of POD, MDA, and SOD
We selected 3-month-old plants and treated them with 150 mM NaCl and 300 mM NaCl for 12 h, 24 h, and 48 h. Plants of 0.5 g were harvested, crushed in liquid nitrogen, and enzymes were extracted with extraction buffer. Subsequently, the extraction was centrifuged at 12,000× g for 10 min, and the supernatant was used for the enzyme activity measurement. Cell superoxide dismutase (SOD), malondialdehyde (MDA), and peroxidase (POD) were measured (Solarbio, Beijing, China).

RNA Extraction and Illumina Deep Sequencing
A Plant Total RNA Extraction Kit was used to extract total RNA (TIANGEN, Beijing, China). The purity and integrity of RNA were checked by nanodrop, agarose gel electrophoresis, and Agilent Bioanalyzer 2100 system (Agilent Technologies Co., Ltd., Beijing, China). The Illumina Hiseq platform was used to sequence the library (LC Sciences, Houston, TX, USA), and a 150 bp double-terminal analysis used to be generated. Clean readings had been mapped to the reference P. euphratica version 2 (v 2.0) genome [56] the usage of HISAT2 [57].

Sequencing Data Processing and Analysis
Gene expression levels were calculated as reads per kilobase of transcript sequence per million base pairs sequenced (FPKM) using Cufflinks [58]. Differentially expressed genes (DEGs) were identified by DESeq [59]. A threshold of padj ≤ 0.05 was used to retrieve DEGs.
To investigate the biological function of DEGs, Blast2GO was used for gene set enrichment and Gene Ontology (GO) terminology [60]. The KEGG (Kyoto Encyclopedia of Genes and Genomes) tool was used to investigate metabolic pathways [61]. The co-expression network was examined using the R package WGCNA [62]. The graphical network created by Cytoscape 3.9.0 software [63].

qRT-PCR Analysis of Illumina Sequencing Results
We selected eight candidate DEGs for quantitative Real-Time PCR, which was used to confirm the accuracy of the RNA-Seq analysis. RNA was extracted and cDNA was reverse transcribed with HiScript qRT Super Mix (+gDNA wiper) (TIANGEN, Beijing, China). An amount of 1 µ cDNA was used as a template for quantitative Real-Time PCR (qRT-PCR), amplified using a Roche Light Cyclerreg 96 Real-Time PCR detection system and SYBR Premix Ex Taq™ (Takara, Dalian, China). The initial denaturation was performed at 95 • C for 5 min, followed by 35 denaturation cycles of 30 s at 95 • C, 30 s of annealing, and 20 s of extension at 72 • C. The primers were created using Primer Premier 5 software (PREMIER Biosoft, Palo Alto, CA, USA) and are listed in Supplementary Table S8. For each gene, three biological replicates were employed. The relative expression levels of the DEGs were calculated using the2 −∆∆ct method and normalized to the expression of the internal reference gene actin [64].

Plasmid Construction and P. alba × P. glandulosa Transformation
To test the functioning of the candidate genes, we used gene transformation in P. alba × P. glandulosa. The coding sequences of the candidate genes were cloned into pDONR222.1 and sequenced. Under the direction of the CaM35S, the coding sequence of pMDC32 was cloned by the Gateway system (Invitrogen, Waltham, MA, USA), and 35S: PeERF1 and 35S: SRDX-PeERF1 plant expression vectors were obtained. P. alba × P. glandulosa were genetically transformed using the callus transformation technique. More than 15 transgenic lines were developed after screening on media containing 50 mg/L hygromycin (Invitrogen, Waltham, MA, USA). Finally, three distinct transgenic lines with high levels of candidate genes were employed in further investigations.

Subcellular Localization and Self-Activation Experiment of PeERF1
The coding region of PeERF1 lacking a termination codon was fused to the 5'terminus of the coding region of GFP to form PeERF1-GFP fusion protein. The proper coding sequence was cloned into pEarley Gate 101(C-YFP) under the control of the CaM35S utilizing the Gateway system (Invitrogen, Waltham, MA, USA), and cloned into pDONR222.1 for sequencing. The primers are shown in Supplementary Table S8. Agrobacterium tumefaciens with fusion vectors 35S: PeERF1-YFP was immersed in tobacco leaves by injection and cultured in darkness at room temperature for 36 h. The epidermises were then stained in the dark for 15 min with 100 ng/mL DAPI. The epidermis was dyed, rinsed three times with phosphate-buffered saline, and the fluorescence was analyzed using a confocal laser scanning microscope (LSM 700, Zeiss, Germany).
To verify whether PeERF1 has self-activation activity, we carried out a yeast two-hybrid test. we used the Yeast maker™ Yeast Transformation System 2 (Takara, Dalian, China) and performed yeast transformation according to the instructions in the user manual. The effector vector and reporter vector were tested in the Y2H Gold yeast strain using the SD/-Trp solid medium. As a negative control, pGBKT7-Rec2 transformed yeast cells were used. To examine the interactions, the successfully transformed yeast cells were cultured in SD/-Trp liquid medium and then coated on SD/-Trp and SD/-Trp and SD/-Trp/X-a-gal (20 mM) solid medium.

Promoter Assay of PeERF1
To investigate the tissue localization of PeERF1 gene expression, we transformed P. alba × P. glandulosa using PMDC164 plasmid vector containing the PeERF1 promoter. We use about 1524 bp upstream sequence of PeERF1 gene in P. euphratica genome for primer design. The primers are in Supplementary Table S8, and the Gateway method was used to construct the GUS promoter gene reporter vector. Transgenic callus of poplar was produced using the previously published Agrobacterium-mediated technique. The activity of the GUS enzyme was examined using a type microscope.

One-Hybrid Yeast Experiment
To verify whether PeERF1 can bind to the GCC-box (AGCCGCC), DRE motif (TACCGA-CAT), TTG1 motif (TTGTTTTGTT), TTG2 motif (TTTTTTTGTT), and their mutant elements GCC-m1 (AGTTGCC), GCC-m2 (ATCCTCC), GCC-m3 (TTTTTTTT), DRE-m1 (TATTGACAT), DRE-m2 (TACCTTCAT), TTG1-m1 (TTGCCTTGTT), TTG1-m2 (TTGTTTTGCC), TTG2-m1 (TTTTTCCGTT), and TTG2-m2 (TTTTTTTGCC), we performed yeast one-hybrid assay. The core sequences of the above four components and their corresponding base mutation sequences were concatenated in three repeats, and EcoR I and Sac I restriction sites were introduced into the upstream and downstream of each concatenated fragment to design primers. We performed yeast transformation using the Yeast maker™ (Takara, Dalian, China). The effector and reporter vectors were introduced into the Y187 yeast strain and tested on DDO solid medium. The pHIS2-p53/pGADT7-Rec2-p53 as positive control and pHIS2-p53/pGADT7Rec2-PeERF1 as negative control co-transformed yeast cells were utilized. To examine the interactions, the successfully transformed yeast cells were grown in DDO liquid medium and then coated on DDO and TDO/3-AT (75 mM) solid medium, and then tested for spots. 4.11. Analysis of Physiological Characteristics of Transgenic P. alba × P. glandulosa Plants On 1/2 MS media with 0, 50, 75, and 100 mM NaCl, transgenic and wild type (WT) were cultivated. The root length, plant height, and fresh weight of plants at one month old were measured using ten biological duplicates. In the greenhouse, one-month-old WT and transgenic plants were placed in soil containers for a month. They were subsequently given a treatment of 150 mM NaCl for 7 days, with water as a control. We assessed the POD, SOD, MDA, electrolyte leakage rate (LR), and leaf chlorophyll content (CC) after the stress treatments. The LR was determined using a DDSJ-308A conductivity meter and following the manufacturer's recommendations. According to Li's approach, the CC of several plants was calculated [29]. The activities of SOD, POD, MDA, and Soluble sugar content were measured by following the manufacturer's instructions (Beijing Solarbio Technology Co., Ltd., Beijing, China). Kumar and Huang's approach was used for staining of histochemical using 3,3 -diaminobenzidine (DAB) and Nitrotetrazolium blue chloride (NBT) [65]. Total RNA was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Eukaryotic mRNA was enriched by Oligo(dT) beads (Epicentre, Madison, WI, USA). The purified double-stranded cDNA fragments were prepared using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB #7530, New England Biolabs, Ipswich, MA, USA), which included end repair, A base addition, and ligation to Illumina sequencing adapters. The AMPure XP Beads (1.0×) were used to purify the ligation process. Ligated fragments were size-selected using agarose gel electrophoresis and amplified using polymerase chain reaction (PCR). Gene Denovo Biotechnology Co. sequenced the generated cDNA library on an Illumina Novaseq6000 (Guangzhou, China).

Gene Expression Characterization Using RNA-Seq
The short read alignment tool Bowtie2 [59,66] (version 2.2.8) was used to map reads to the ribosomal RNA (rRNA) database. A reference genome index was constructed, and paired-end clean reads were mapped to it with HISAT2.2.4 [57] with the "-rna-strandness RF" and other parameters set to default. StringTie v1.3.1 was used to create the mapped reads of each sample in a reference-based technique [67,68]. For each transcription area, an FPKM value was computed using the RSEM program to assess the amount and variability of its expression [69].

Conclusions
In this study, we obtained the DEGs in P. euphratica leaves with different salt treatments, which were involved in biological processes and pathways, and screened eight potential genes based on the co-expression network and GWAS results. Under salt stress, we found the accumulation of osmotic adjustment chemicals and an increase in antioxidant enzyme activity in P. euphratica. Many DEGs were shown to be enriched in stress-related activities such as transcriptional regulation, stress response, cell death, and so on. PeERF1 has been shown to improve salt tolerance in transgenic plants by overexpression and dominant inhibitory expression in P. alba × P. glandulosa. Our results will contribute to resolve the salt tolerance mechanism of P. euphratica, help to identify and breed poplar salt tolerance candidate genes, provide material for future related studies, and also lay the foundations for plant stress tolerance.