Characterization of the Gene Expression Profile Response to Drought Stress in Populus ussuriensis Using PacBio SMRT and Illumina Sequencing

In this study, we characterized the gene expression profile in the roots of Populus ussuriensis at 0, 6, 12, 24, 48 and 120 h after the start of polyethylene glycol (PEG)-induced drought stress using PacBio single-molecule real-time sequencing (SMRT-seq) and Illumina RNA sequencing. Compared to the control, 2244 differentially expressed genes (DEGs) were identified, and many of these DEGs were associated with the signal transduction, antioxidant system, ion accumulation and drought-inducing proteins. Changes in certain physiological and biochemical indexes, such as antioxidant activity and the contents of Ca2+, proline, and total soluble sugars, were further confirmed in P. ussuriensis roots. Furthermore, most of the differentially expressed transcription factors were members of the AP2/ERF, C2H2, MYB, NAC, C2C2 and WRKY families. Additionally, based on PacBio SMRT-seq results, 5955 long non-coding RNAs and 700 alternative splicing events were identified. Our results provide a global view of the gene expression profile that contributes to drought resistance in P. ussuriensis and meaningful information for genetic engineering research in the future.


Introduction
Drought is one of the major natural disasters facing plants. In recent years, due to the influence of the global climate, the frequency of droughts has been increasing [1]. In plants, the root is the first organ to sense drought stress, and it plays a vital role in the response of plants to drought stress [2]. Hence, the roots are important in maintaining a good growth in plants, especially in soils with insufficient moisture or nutrients. During water deficit, the root perceives the soil water status, promotes biosynthesis and signal transduction of abscisic acid (ABA), and promotes changes in the expression level of drought-related genes [3,4].
Plants are sessile organisms with complex mechanisms to cope with drought stress. Strategies for dealing with drought stress in plants include (1) changes in the structure of the root system at the phenotypic level, (2) changes in osmotic regulation and antioxidant metabolism at the physiological level, and (3) regulation of drought signaling, transcription factors (TFs), and cellular defense at the molecular level [5][6][7]. Studies have exhaustively investigated molecular responses to drought stress in plants, such as Populus davidiana [8], Gossypium hirsutum and Gossypium tomentosum [9], Arabidopsis thaliana [10], and Oryza sativa [11]. TFs play an important role in the molecular response mechanism under drought stress, such as ABA-responsive element (ABRE), myeloblastosis (MYB), ethylene-responsive factor (ERF) and WRKY [12][13][14]. Additionally, synthesis of certain specific proteins, such as heat shock protein (HSP) and late embryogenesis abundant (LEA) protein, enhances the drought tolerance of plants [15,16]. However, the gene network involved in the response For Illumina RNA-seq, a total of 816,789,418 clean reads were obtained, with a clean Q30 base rate > 93% for all samples. Then, the clean reads were successfully mapped back to the full-length transcripts of PacBio SMRT-seq. As shown in Table S1, the average mapping ratio is 69.3%, the average multimap ratio is 28.8%.

SSR and lncRNA Prediction
In this study, 3121 SSRs were identified using the MISA software (Table S3). As shown in Figure 1C, the type distribution included 1252 dinucleotides (40.1%), 1673 trinucleotides (53.6%), and 196 tetranucleotides. As shown in Figure 1D, Coding-non-Coding Index (CNCI), Coding Potential Assessment Tool (CPAT), and Coding Potential Calculator (CPC) software predicted 11,484,11,294, and 8481 lncRNAs, respectively. For Illumina RNA-seq, a total of 816,789,418 clean reads were obtained, with a clean Q30 base rate > 93% for all samples. Then, the clean reads were successfully mapped back to the full-length transcripts of PacBio SMRT-seq. As shown in Table S1, the average mapping ratio is 69.3%, the average multimap ratio is 28.8%.

Identification of Alternative Splicing
In total, 700 AS events were identified using the AStalavista tool (Table S3). As shown in Figure 2A, the majority of AS events were of retained intron (RI), alternative 3' splice site (A3SS), or alternative 5' splice site (A5SS) types. Gene ontology (GO) enrichment analysis was performed on AS genes; the top 10 functional groups of biological process (BP), molecular function (MF), and cellular component (CC) categories are

Identification of Alternative Splicing
In total, 700 AS events were identified using the AStalavista tool (Table S3). As shown in Figure 2A, the majority of AS events were of retained intron (RI), alternative 3' splice site (A3SS), or alternative 5' splice site (A5SS) types. Gene ontology (GO) enrichment analysis was performed on AS genes; the top 10 functional groups of biological process (BP), molecular function (MF), and cellular component (CC) categories are shown in Figure  S1. The top three GO enrichment results for BP category were "transcription", "regulation of transcription" and "mRNA processing", and the top three for MF category were "ATP binding", "metal ion binding" and "DNA binding". As shown in Figure 2B, 281 AS genes were mapped to the 37 pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.
shown in Figure S1. The top three GO enrichment results for BP category were "transcription", "regulation of transcription" and "mRNA processing", and the top three for MF category were "ATP binding", "metal ion binding" and "DNA binding". As shown in Figure 2B, 281 AS genes were mapped to the 37 pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Transcript Integrity Analysis and Gene Annotation
The HQ isoforms generated by PacBio SMRT-seq were used to open reading frames (ORFs) identification after integrity assessment using single-copy orthologs (BUSCO). A total of 34,522 ORFs (297-6624 b.p. long) were identified using TransDecoder. The average length of ORFs is 1134 b.p., and the N50 of ORFs is 1473 b.p. ( Figure 2C). After removing the redundancy of ORFs, 14,136 genes were annotated using Trinotate (Table  S5). For further functional prediction and categorization, 8967 genes were annotated using the euKaryotic Orthologous Groups (KOG) method ( Figure S2). For GO enrichment analysis, 10,692 ORFs were successfully enriched into BP (9211), CC (9383), and MF (9290) categories ( Figure S3).

Identification and Analysis of DEGs
The differentially expressed genes (DEGs) were determined using the DEGseq software, which is based on the read count values of each transcript. In total, 2244 DEGs (1328 DEGs were up-regulated and 959 DEGs were down-regulated) were identified under PEG-induced drought stress ( Figure 3B, Table S6). Our results demonstrated that 871 (65.59%) and 786 (59.19%) DEGs were up-regulated at 6 or 12 h, respectively, while 679 (70.8%) and 429 (44.73%) DEGs were down-regulated at 6 or 12 h, respectively (Fig-

Transcript Integrity Analysis and Gene Annotation
The HQ isoforms generated by PacBio SMRT-seq were used to open reading frames (ORFs) identification after integrity assessment using single-copy orthologs (BUSCO). A total of 34,522 ORFs (297-6624 b.p. long) were identified using TransDecoder. The average length of ORFs is 1134 b.p., and the N50 of ORFs is 1473 b.p. ( Figure 2C). After removing the redundancy of ORFs, 14,136 genes were annotated using Trinotate (Table S5). For further functional prediction and categorization, 8967 genes were annotated using the euKaryotic Orthologous Groups (KOG) method ( Figure S2). For GO enrichment analysis, 10,692 ORFs were successfully enriched into BP (9211), CC (9383), and MF (9290) categories ( Figure S3).

Identification and Analysis of DEGs
The differentially expressed genes (DEGs) were determined using the DEGseq software, which is based on the read count values of each transcript. In total, 2244 DEGs (1328 DEGs were up-regulated and 959 DEGs were down-regulated) were identified under PEG-induced drought stress ( Figure 3B, Table S6). Our results demonstrated that 871 (65.59%) and 786 (59.19%) DEGs were up-regulated at 6 or 12 h, respectively, while 679 (70.8%) and 429 (44.73%) DEGs were down-regulated at 6 or 12 h, respectively ( Figure 3B). In the down-regulated category, only 29 genes showed decreased expression at all time points, and in the up-regulated category, only 89 genes showed decreased expression at all time points. A heatmap, which can visualize the changing pattern of DEGs under different time points of PEG treatment was constructed based on the read count values ( Figure 3A). ure 3B). In the down-regulated category, only 29 genes showed decreased expression at all time points, and in the up-regulated category, only 89 genes showed decreased expression at all time points. A heatmap, which can visualize the changing pattern of DEGs under different time points of PEG treatment was constructed based on the read count values ( Figure 3A).  In GO enrichment analysis of all DEGs, 1607, 1604, and 1520 DEGs were enriched in the BP, MF, and CC categories, respectively ( Figure S4). A total of 662 DEGs were mapped to 134 KEGG pathways, and the metabolic pathways regulated by PEG-induced drought stress were screened out (p < 0.05; Figure 4). Our results demonstrated that metabolic pathways related to drought stress were more significantly enriched in the early stages (6, 12, and 24 h) than in the late stages (48 and 120 h). The phenylpropanoid biosynthesis pathway was significantly enriched at all time points, and the MAPK signaling, flavonoid biosynthesis, and amino acid metabolism pathways were significantly enriched at most of the time points.
In GO enrichment analysis of all DEGs, 1607, 1604, and 1520 DEGs were enriched in the BP, MF, and CC categories, respectively ( Figure S4). A total of 662 DEGs were mapped to 134 KEGG pathways, and the metabolic pathways regulated by PEG-induced drought stress were screened out (p < 0.05; Figure 4). Our results demonstrated that metabolic pathways related to drought stress were more significantly enriched in the early stages (6, 12, and 24 h) than in the late stages (48 and 120 h). The phenylpropanoid biosynthesis pathway was significantly enriched at all time points, and the MAPK signaling, flavonoid biosynthesis, and amino acid metabolism pathways were significantly enriched at most of the time points.

DEGs Responding to Drought Stress in P. ussuriensis
DEGs associated with the response to drought stress in plants are mainly involved in drought signal transduction and defense system activation. All DEGs responding to drought stress are shown in Tables S7-S11, and some of these are listed in Table 2. The perception and transmission of drought signals by P. ussuriensis plays a crucial role in the response to drought stress. In this study, we identified 38 DEGs involved in the transmission of drought signals in P. ussuriensis (Table S7). Of these, 14 DEGs, including two 9-cis-epoxycarotenoid dioxygenase (NCED) homologous genes, four PYR1-like (PYL) homologous genes, and two protein phosphatase 2C (PP2C) homologous genes, were involved in ABA synthesis and binding. Of the total 38 DEGs, 20 DEGs, including calmodulin-like (CML), calmodulin-binding protein, and calcium-dependent protein kinase (CPK/CDPK) homologous genes, were significantly differentially expressed under drought stress. In addition, we found that four mitogen-activated protein kinase kinase kinase (MAPKKK) homologous genes were significantly up-regulated under PEG-induced drought stress.

DEGs Responding to Drought Stress in P. ussuriensis
DEGs associated with the response to drought stress in plants are mainly involved in drought signal transduction and defense system activation. All DEGs responding to drought stress are shown in Tables S7-S11, and some of these are listed in Table 2. The perception and transmission of drought signals by P. ussuriensis plays a crucial role in the response to drought stress. In this study, we identified 38 DEGs involved in the transmission of drought signals in P. ussuriensis (Table S7). Of these, 14 DEGs, including two 9-cis-epoxycarotenoid dioxygenase (NCED) homologous genes, four PYR1-like (PYL) homologous genes, and two protein phosphatase 2C (PP2C) homologous genes, were involved in ABA synthesis and binding. Of the total 38 DEGs, 20 DEGs, including calmodulin-like (CML), calmodulin-binding protein, and calcium-dependent protein kinase (CPK/CDPK) homologous genes, were significantly differentially expressed under drought stress. In addition, we found that four mitogen-activated protein kinase kinase kinase (MAPKKK) homologous genes were significantly up-regulated under PEG-induced drought stress.
Furthermore, 56 DEGs related to reactive oxygen species (ROS) scavenging were identified under PEG-induced drought stress, including those encoding catalase (CAT), glutathione S-transferase (GST), peroxidase (POD) and superoxide dismutase (SOD) ( Table S8). Of these, most of the GST homologous genes were significantly up-regulated after 6, 12, and 24 h of PEG-induced drought stress. In P. ussuriensis, a total of 28 POD homologous genes were also found to be significantly differentially expressed under PEG-induced drought stress. After PEG treatment, 13 DEGs related to ion transport were found in the roots of P. ussuriensis. The expression levels of two genes (f7p60_2856_19557 and f2p60_2781_19668) encoding potassium transporter (KT) homologous genes, one gene (f2p37_2396_20662) encoding cyclic nucleotide-gated ion channel (CNGC) homologous gene, two genes (f2p55_3201_1361 and f2p60_3228_18926) encoding glutamate receptor (GLR) homologous genes, two genes (f3p60_1980_22126 and f4p60_1851_22830) encoding cation exchanger (CAX) homologous genes, and two genes (f2p57_4009_433 and f2p20_3494_18605) encoding calcium-transporting ATPase (ACA) homologous genes were significantly different.

Transcription Factors Responding to Drought Stress in P. ussuriensis
According to reports, many TFs are affected by drought conditions. During PEGinduced drought stress, 122 TFs were differentially expressed, which belonged to 18 different families (Table S12). Of these, 10 TF families accounted for 86% of stress response, including AP2/ERF, C2H2, MYB, NAC, C2C2, WRKY, HB, bHLH, AUX/IAA, and C3H ( Figure 5A). A heatmap was constructed for visualizing the change of TFs at all time points of the PEG-induced drought stress treatment ( Figure 5B). The number of TFs was greater at 6 and 12 h than that at 48 and 120 h, probably because many TFs were activated in the early stages of drought stress. As shown in Figure 5C, we selected several TFs (e.g., TFs encoding ERF017, ZAT11, HB7, MYB58, NAC35, and WRKY75 homologous genes) to conduct quantitative reverse transcription-polymerase chain reaction (qRT-PCR) analysis.

Validation of Gene Expression Levels
The results of qRT-PCR analysis of 12 genes were consistent with that of Illumina RNA-seq data (Figure 6), indicating the reliability of the Illumina RNA-seq data.

Validation of Gene Expression Levels
The results of qRT-PCR analysis of 12 genes were consistent with that of Illumina RNA-seq data ( Figure 6), indicating the reliability of the Illumina RNA-seq data.

Changes in Physiological Indexes in Response to Drought Stress
We measured the H2O2 content and the activities of GST, SOD, and POD for evaluating the drought-induced phenotypes. After PEG-induced drought stress, from 6 to 48 h, the H2O2 content levels in P. ussuriensis were significantly increased; however, the levels dropped significantly at 120 h ( Figure 7A). POD, SOD, CAT, and GST activities showed a similar trend after PEG-induced drought stress by significantly increasing at 6, 12, 24, 48, and 120 h and reaching the highest at 24 h ( Figure 7B-E).
After PEG-induced drought stress, the proline content levels in P. ussuriensis significantly increased at 24, 48, and 120 h ( Figure 7F). The total soluble sugar content exhib-

Changes in Physiological Indexes in Response to Drought Stress
We measured the H 2 O 2 content and the activities of GST, SOD, and POD for evaluating the drought-induced phenotypes. After PEG-induced drought stress, from 6 to 48 h, the H 2 O 2 content levels in P. ussuriensis were significantly increased; however, the levels dropped significantly at 120 h ( Figure 7A). POD, SOD, CAT, and GST activities showed a similar trend after PEG-induced drought stress by significantly increasing at 6, 12, 24, 48, and 120 h and reaching the highest at 24 h ( Figure 7B-E). tained at a high level at 12, 24, and 48 h ( Figure 7G). The malondialdehyde (MDA) content in P. ussuriensis gradually increased from 0 to 120 h and peaked at 120 h ( Figure  7H).

Measurements of Minerals Content
In this study, we found that some genes related to ion transport were significantly different in expression after PEG treatment. Therefore, the mineral content (Ca, K, and Na) in the roots of P. ussuriensis was measured. The potassium (K), sodium (Na), and calcium (Ca) contents in the roots of P. ussuriensis were measured using inductively coupled plasma optical emission spectrometry (ICP-OES). The K and Ca contents in the roots of P. ussuriensis were significantly increased after PEG-induced drought stress and maintained at a high level at 120 h ( Figure 8A,C). However, during PEG-induced drought stress, the Na content in the roots of P. ussuriensis did not significantly increase at 6, 12, and 24 h and sharply increased at 120 h ( Figure 8B). In addition, the K and Ca levels in the roots of P. ussuriensis were much higher than Na levels. After PEG-induced drought stress, the proline content levels in P. ussuriensis significantly increased at 24, 48, and 120 h ( Figure 7F). The total soluble sugar content exhibited a significant increase during PEG-induced drought stress treatment and was maintained at a high level at 12, 24, and 48 h ( Figure 7G). The malondialdehyde (MDA) content in P. ussuriensis gradually increased from 0 to 120 h and peaked at 120 h ( Figure 7H).

Measurements of Minerals Content
In this study, we found that some genes related to ion transport were significantly different in expression after PEG treatment. Therefore, the mineral content (Ca, K, and Na) in the roots of P. ussuriensis was measured. The potassium (K), sodium (Na), and calcium (Ca) contents in the roots of P. ussuriensis were measured using inductively coupled plasma optical emission spectrometry (ICP-OES). The K and Ca contents in the roots of P. ussuriensis were significantly increased after PEG-induced drought stress and maintained at a high level at 120 h ( Figure 8A,C). However, during PEG-induced drought stress, the Na content in the roots of P. ussuriensis did not significantly increase at 6, 12, and 24 h and sharply increased at 120 h ( Figure 8B). In addition, the K and Ca levels in the roots of P. ussuriensis were much higher than Na levels. Confocal images of Fluo-4 AM (Biyuntian Biotechnology Co., Ltd., Shanghai, China) fluorescence showed an increase in the fluorescence intensity of meristematic cells in roots after PEG treatment, which suggested a higher free Ca 2+ accumulation in PEG treatment groups ( Figure 9A-F). In the control group, the Ca 2+ accumulation in the root tip cells was mainly distributed on the cell wall and nucleus ( Figure 9A). With the increase in time for PEG-induced drought stress treatment, the Ca 2+ accumulation in the root tip cells gradually and evenly filled the cytoplasm.  Confocal images of Fluo-4 AM (Biyuntian Biotechnology Co., Ltd., Shanghai, China) fluorescence showed an increase in the fluorescence intensity of meristematic cells in roots after PEG treatment, which suggested a higher free Ca 2+ accumulation in PEG treatment groups ( Figure 9A-F). In the control group, the Ca 2+ accumulation in the root tip cells was mainly distributed on the cell wall and nucleus ( Figure 9A). With the increase in time for PEG-induced drought stress treatment, the Ca 2+ accumulation in the root tip cells gradually and evenly filled the cytoplasm. Confocal images of Fluo-4 AM (Biyuntian Biotechnology Co., Ltd., Shanghai, China) fluorescence showed an increase in the fluorescence intensity of meristematic cells in roots after PEG treatment, which suggested a higher free Ca 2+ accumulation in PEG treatment groups ( Figure 9A-F). In the control group, the Ca 2+ accumulation in the root tip cells was mainly distributed on the cell wall and nucleus ( Figure 9A). With the increase in time for PEG-induced drought stress treatment, the Ca 2+ accumulation in the root tip cells gradually and evenly filled the cytoplasm.

Characteristics of PacBio SMRT-seq and Analysis of lncRNAs and AS Events
The PacBio SMRT-seq can generate long-read or full-length transcripts without splicing fragments. The full-length transcripts obtained from the PacBio SMRT-seq can be used for the identification and analysis of AS events and lncRNAs. The lncRNAs refer to RNA molecules whose transcripts are longer than 200 nt and do not encode proteins. The lncRNAs could function as competing endogenous RNAs (ceRNAs) to participate in response to drought stress in plants, and represent a new gene regulation layer [17,18]. In this study, we identified lncRNAs in P. ussuriensis under PEG-induced drought stress. However, their functions need to be further studied. Some mRNA precursors of eukaryotic genes produce different mRNA splicing isoforms through different splicing methods, which produce proteome and genetic diversity in eukaryotes [19,20]. There have been reports that AS mediates drought response, for example, maize plants overexpressing three different ZmCCA1 isoforms distinguished by AS events showed different tolerance to drought stress [21], and BrSR45a regulates drought stress response by influencing AS of target genes [22]. Although many AS events have been identified in this study, different functions of these AS events on drought stress need to be further investigated. Therefore, in the future, research on the functions of these lncRNAs and AS events will help to better understand the response of P. ussuriensis to drought stress at the molecular level.

Sensing and Transmission of Drought Signals in P. ussuriensis
The drought stress perception and signal transduction in plants is a complex process. Plant cells sense drought stress signals through changes in the cell surface or cell membrane. After sensing the stress signal, plant cells stimulate and induce changes in the concentration of Ca 2+ , and then use calcium effector protein to transmit the signals [9,23]. In this study, the Ca 2+ content in the roots of P. ussuriensis significantly increased after PEG-induced drought stress (Figures 8C and 9). At the gene expression level, the genes encoding Ca 2+ transporters (e.g., ACA, CAX, and CNGC2 homologous genes) were up-regulated after drought stress. The results were consistent with that of the previous study that reported the contribution of AtCNGC2 to Ca 2+ uptake in A. thaliana [24]. Furthermore, we found that the genes encoding CML, calmodulin-binding protein, and CPK/CDPK homologous genes in the roots of P. ussuriensis were differentially expressed after PEG-induced drought stress. Notably, as calcium effector protein, most of these genes were only significantly up-regulated at 6 h and 12 h of PEG-induced drought stress. These results indicated that the transmission of Ca 2+ signal mainly occurs in the early stages of drought stress.
Previous studies reported that ABA signaling plays an important role in plant acclimation to drought stress [9]. Zeaxanthin epoxidase (ZEP) is an important enzyme for the zeaxanthin epoxidation reaction in the ABA biosynthesis pathway, and overexpression of AtZEP enhanced the plant's tolerance to drought stress [25]. NCED is the rate-limiting enzyme of the ABA biosynthesis pathway, and the AtNCED3 gene is highly induced under drought stress [26]. In the roots of P. ussuriensis, one ZEP homologous gene and one NCED3 homologous gene were up-regulated under PEG-induced drought stress, and their expression levels reached the highest value at 12 h of PEG-induced drought stress. This result was consistent with results in Gossypium spp. and Vitis riparia that NCED homologous genes were up-regulated under drought stress, suggesting that ABA signaling plays a role in plant responses to drought stress [9,27]. The ABA receptor PYR/PYL/RCAR, the negative regulator PP2C and the positive regulator SnRK2 are the main members in the ABA signaling pathway [28]. During PEG-induced drought stress, four PYL homologous genes were down-regulated and two PP2C homologous genes were up-regulated in the roots of P. ussuriensis. However, the details of the regulation of PYL-PP2C-SnRK2 in P. ussuriensis under drought stress are unknown.

Antioxidant Defense against Mechanism in Response to Drought Stress
ROS is a byproduct of plant aerobic metabolism, and excessive accumulation of ROS can damage the structure and function of plant cells [29,30]. In this study, after PEG-induced drought stress, the H 2 O 2 content in the roots of P. ussuriensis was significantly increased (Figure 7A), and the count of MDA was also significantly increased ( Figure 7H). These results indicated that the excessive accumulation of ROS might be one of the factors that cause severe damage to the cell membrane of P. ussuriensis under drought stress. Under abiotic stress, plants control the overproduction of ROS through enzymatic components and non-enzymatic antioxidants [29,30]. During drought stress, relatively high activity of enzymes (SOD, POD, CAT, and GST), which decompose ROS to protect cells from oxidative damage, was observed in the roots of P. ussuriensis ( Figure 7B-E). The differential expression of ROS scavenging-related genes was observed in the roots of P. ussuriensis after PEG-induced drought stress (Table S7), which provides evidence at the molecular level. Among them, more ROS clearance-related genes expression levels were significantly changed in the early stages (6, 12 and 24 h) of PEG-induced drought stress. These results on ROS in P. ussuriensis under drought stress are consistent with that of research reports on Pinus halepensis, Malus pumila, and Hevea brasiliensis [31][32][33].

Role of TFs in P. ussuriensis in Response to Drought Stress
TFs play a role in regulating the activation or inactivation of drought-responsive genes, which greatly affects the response of plants to drought stress [34,35]. In this study, the number of significantly differentially expressed TFs was highest at 6 h of PEG-induced drought stress, and gradually decreased along with the progressive drought (Table S11). This result indicated that most TFs play a role in the early stages of drought stress. Similar results were obtained in A. thaliana, O. sativa, and P. davidiana, indicating that TFs are mainly involved in drought stress responses at 6 h and 12 h [8,36,37]. In this study, most of the differentially expressed TFs were members of the AP2/ERF, C2H2, MYB, NAC, C2C2, and WRKY TF families ( Figure 5A). A substantial amount of evidence has shown that these TF subfamilies were involved in plant drought resistance [12,[38][39][40].
Some TF family members play an important role in drought resistance through ABAmediated pathways [10,41]. In-depth study of the WRKY TF family in response to drought stress has been carried out in many species. For example, overexpression of GsWRKY20 in A. thaliana enhanced the drought tolerance by mediated ABA signaling, and overexpression of BdWRKY36 in tobacco enhanced the drought tolerance and ABA biosynthesis [10,41]. In the present study, 10 WRKY TF homologous genes were differentially expressed, suggesting a specific role for these WRKY TFs in drought tolerance. MYB is one of the largest plant TF families, which positively modulates drought stress through ABA-mediated pathways. The overexpression of SiMYB75 in A. Thaliana improved the plant's tolerance to drought stress. Further analysis showed that ABA content was increased, and ABA-related genes were differentially expressed in SiMYB75-overexpressing lines [42]. Fang et al. [43] observed that overexpression of PtrMYB94 in Populus enhanced the drought tolerance and increased the ABA content, suggesting that PtrMYB94 is involved in the regulation of ABA-dependent drought stress in Populus. Several MYB family genes were differentially expressed at 6 h of PEG-induced drought stress, suggesting a regulatory role of MYB in the early stages of drought stress in P. ussuriensis.
Furthermore, in this study, we identified some AP2/ERF and NAC TF family members whose expression levels were significantly altered by PEG-induced drought stress. Similar results were obtained with other plants, such as V. riparia, P. davidiana and L. kaempferi [8,27,38]. Furthermore, the genes encoding DRE1D, ERF5, NAC2, and NAC72/RD26 were upregulated under PEG-induced drought stress in this study. Additionally, these TFs have been shown to have a positive regulation on drought resistance in other plants [44][45][46][47]. Furthermore, the genes that were down-regulated during drought stress also contributed to the adaptation of plants to drought stress. For example, Dong and Liu [48] observed that AtRAP2.1 negatively regulated drought stress in A. thaliana. Moreover, overexpression of PeNAC034 in A. thaliana reduced drought tolerance, and PeNAC045 overexpression line of Populus showed a drought-sensitive behavior [49]. In this study, we found that three RAP homologous genes and five NAC homologous genes were specifically down-regulated during PEG-induced drought stress in the roots of P. ussuriensis (Table S11). Further exploration of these TFs may provide a greater understanding of drought resistance regulation in P. ussuriensis.

Ion Accumulation in P. ussuriensis
During drought stress, the osmotic balance in plant cells is destroyed, and water uptake is disturbed [50]. The changes in the ions in plant cells contribute to osmotic regulation. One unexpected finding of this study was the extent to which Ca 2+ content increased significantly in cells at 48 and 120 h of PEG-induced drought stress (Figure 9). In addition, one CNGC2 homologous gene was up-regulated at all time points of drought stress (Table 2), which mediates the absorption and accumulation of Ca 2+ [51]. These results indicated that, in addition to acting as a second messenger, Ca 2+ might also play an important role in osmotic regulation. Studies in A. thaliana have shown that Ca 2+ contributed to the conversion of water channels from the active to an inactive state, thereby reducing the water permeability of the cytoplasmic membrane of the A. thaliana roots [52].
Prior studies have noted the importance of K + in enhancing osmotic regulation in plants [53]. In this study, the K content in the roots of P. ussuriensis was significantly increased, and two HAK/KT homologous genes were significantly differentially expressed under PEG-induced drought stress (Table 2). Some transporters, namely HvHAK1, HvAKT2, OsHAK1, and OsHAK5, have been shown to be involved in the uptake of K + in the roots [53,54]. Expression analysis in Cajanus cajan revealed that several K + transport genes, that is, CcHAK6, CcHAK15, CcAKT2, and CcHKT1, were up-regulated in response to drought stress [55]. Notably, the K content in the roots of P. ussuriensis was significantly higher than that of Na and Ca under PEG-induced drought stress ( Figure 9A-C). The results showed that the increase in K + absorption might be the first choice for osmotic regulation in P. ussuriensis, under drought stress. Under PEG-induced drought stress, the Na content was significantly elevated in the roots of P. ussuriensis, but it was much lower than the K content levels. These results indicated that K + /Na + /Ca 2+ accumulates in the roots of P. ussuriensis, which may contribute to the osmotic adjustment of the roots and limit the dehydration caused by PEG-induced drought stress.

Proteins, Osmolytes and Lignin Response to Drought Stress in P. ussuriensis
Along with antioxidant and ion accumulation, accumulation of proteins, osmolytes (such as proline and soluble sugars), and lignin plays an important role in drought stress. LEA and HSP proteins are important drought-inducing proteins in plants. LEA proteins belong to a large protein family of hydrophilic proteins that are closely associated with resistance to drought stress [16,56]. In-depth studies of the LEA protein family in response to drought stress have been carried out in other plants. For example, transgenic rice lines overexpressing OsLEA3-1 have enhanced drought resistance compared to wild-type [57], and overexpression of IpLEA confers greater drought tolerance in A. thaliana [58]. In this study, we identified some of the LEA genes that were differentially expressed in the roots of P. ussuriensis under PEG-induced drought stress, especially five dehydrin genes that were up-regulated at all time points of PEG-induced drought stress (Table S10). Furthermore, most of the DEGs encoding HSPs were up-regulated in the roots of P. ussuriensis under PEG-induced drought stress. The result is consistent with previous studies on Capsicum annuum, O. sativa, and Solanum lycopersicum showing that HSPs positively regulate the tolerance to drought stress [11,15,59].
The essential role of osmolytes is to limit water loss in plants under drought stress. Proline is a common plant organic osmotic agent, which is closely related to drought stress [50]. Proline content increased under drought stress, which has been reported in A. thaliana, Populus euphratica, and Populus yunnanensis [60][61][62]. One of the most common permeates in plants is soluble sugar, and the accumulation of total soluble sugar contributes directly to drought tolerance. For example, drought significantly induced the accumulation of total soluble sugar in Michelia macclurei and Schima superba [63]. In addition, Regier et al. [64] observed that total soluble sugar content increased in the roots of water-limited Populus nigra compared to well-watered P. nigra. Our results clearly showed that the proline content and total soluble sugar content of P. ussuriensis were significantly increased in response to PEG-induced drought stress ( Figure 7F,G). In addition, the amino acid biosynthesis; cysteine and methionine metabolism; and glycine, serine, and threonine metabolism pathways were significantly enriched at 6 and 12 h of PEG-induced drought stress (Figure 4). The results suggested that proline and total soluble sugar levels positively regulated the drought stress tolerance in P. ussuriensis.
As the main component of the secondary cell wall, lignin plays an important role in plant resistance to drought stress [65,66]. Recently, it has been demonstrated that the increase in lignin enhanced the tolerance of P. ussuriensis to drought stress [67]. The lignin monomers were carried out through the phenylpropanoid biosynthesis pathway. Moreover, in the present study, we found that the phenylpropanoid biosynthesis pathway was significantly enriched under PEG-induced drought stress ( Figure 4). Increased expression of genes involved in lignin biosynthesis, such as PAL, 4CL, CAD, CCR, C4H, and COMT homologous genes, was observed in the roots of P. ussuriensis under PEG-induced drought stress (Table S11). Similar results were obtained in Paeonia ostii [68] and L. kaempferi [38], which suggested that plants might promote lignin biosynthesis resulting in enhanced tolerance to drought stress.

Plant Material and Drought Treatment
In this study, the tissue culture seedlings of P. ussuriensis were used as research material. The stem (3-cm long) of P. ussuriensis was cut and cultured in half-strength Murashige and Skoog (MS) medium for rooting and grown in a tissue-culture chamber (temperature 25 ± 2 • C, photoperiod of 16 h light and 8 h dark). After three weeks, the tissue culture seedlings of P. ussuriensis had a complete root development ( Figure S5). Three-week-old P. ussuriensis seedlings with uniform root growth were then transferred to liquid media containing 6% PEG 6000 (Sigma-Aldrich, St. Louis, MO, USA). After treatment of 0, 6, 12, 24, 48, and 120 h, the roots were harvested at the same time and frozen in liquid nitrogen. Each sample was from at least six different seedlings, and three biological replicates were collected at each treatment.

Preparation of Full-Length cDNA Library and Illumina RNA-seq Library and Sequencing
Total RNA was isolated as described previously [69], with a slight modification. Then, RNA integrity was assessed according to the requirement of 28S:18S ≥ 1.4, RNA Integrity Number (RIN) ≥ 8. The total RNA from 18 samples was combined in equal amounts for cDNA synthesis and SMRTbell library construction. The cDNA synthesis was performed using the Clontech SMARTer PCR cDNA synthesis kit (Takara Biotechnology, Dalian, China), and the SMRTbell library was constructed using the PacBio SMRTbell Template Prep Kit (PacBio, Menlo Park, CA, USA). In total, 2 µg of RNA was used to generate sequencing libraries using the NEBNext Ultra RNA Library Prep Kit (#E7530L, NEB, Ipswich, MA, USA). Libraries were sequenced in 150 b.p. paired-end mode, using an Illumina HiSeq X Ten system. The above sequencing was performed by Zhejiang Annoroad Biotechnology Co., Ltd. (Beijing, China).

Quality Control and Error Correction of PacBio SMRT-seq and RNA-seq
Quality control of the raw data was performed to avoid inaccuracy in subsequent information analysis. After data quality control, the data volume of subreads was calculated. The consensus sequence generated by multiple subread sequences are CCS reads, and the CCS reads are further corrected to obtain HQ isoforms (with accuracy greater than 0.99) and LQ isoforms. The isoforms with a large number of redundant sequences were clustered together using an isoform-level clustering algorithm known as iterative clustering for error correction (ICE) to obtain a new consensus isoform [70]. The consensus isoform obtained after ICE was used as the reference sequence, and the short-read sequences obtained from 18 Illumina RNA-seq libraries were compared to the reference sequence using RNA-seq by Expectation Maximization (RSEM) (v1.2.31) software (Li, Dewey, WI, USA) [71].

Transcript Integrity Analysis and Gene Annotation
In this study, the integrity of transcripts generated by PacBio SMRT-seq was analyzed using benchmarking universal single-copy orthologs (BUSCO) [77]. The ORFs were predicted using TransDecoder version 3.0.1 [78]. Then, the functional annotations of predicted ORFs were combined using Trinotate, which uses multiple well-referenced methods for functional annotation.

Analysis of Differentially Expressed Transcripts
Using the results from the PacBio SMRT-seq as a reference, the differential gene expression analysis between different treatment groups was performed using the DEGSeq software [79]. Genes satisfying the conditions of |log2 fold change rate| ≥ 1 and p-value < 0.05 were identified as DEGs. The visualization of the expression patterns of DEGs is achieved by building a heatmap using the pheatmap package (https://cran.r-project.org/ web/packages/pheatmap, accessed on 6 April 2021). GO enrichment analysis and KEGG pathway enrichment analysis of DEGs using GOseq R package and KOBAS software, respectively [80,81].

Determination of Physiological Indexes
The MDA content was measured using a plant MDA assay kit. The total soluble sugar content was measured using a plant soluble sugar content test kit. The proline content was measured using a proline assay kit. The CAT activity was measured using a CAT assay kit. The GST activity was measured using a GST assay kit. The SOD activity was measured using a total SOD assay kit. The POD activity was measured using a POD assay kit. All assay kits were purchased from Nanjing Jiancheng Bioengineering Institute, China. These data were averaged from three biological replicates.

Measurements of Minerals Content
The root samples were taken in triplicate at 0, 6, 12, 24, 48, and 120 h after PEG-induced drought stress. The whole roots were dried using an oven. The minerals content (Ca, K, and Na) in the roots of P. ussuriensis was assayed using ICP-OES, as described previously, with slight modifications [82]. These data were averaged from three biological replicates. In addition, according to the protocol outlined by Li et al. [14], the free Ca 2+ content in root tip cells was measured using calcium indicator Flou-4 AM (Biyuntian, Shanghai, China), and observed through a confocal laser scanning microscope imager 700 (Carl Zeiss, Docuval, Germany).

Statistical Analysis of Data
In this study, statistical analysis of data was performed using statistical package for social sciences (SPSS) version 20 (IBM, Chicago, IL, USA), and analysis of significance was performed using analysis of variance (ANOVA). In addition, correction for multiple testing using the Bonferroni approach. Significant differences in relative levels were represented by different lowercase letters, p < 0.05.

Conclusions
The purpose of the current study was to determine the responses of P. ussuriensis root to drought stress using physiochemical measurements, SMRT-seq and RNA-seq analysis. A hypothetical model of the drought tolerance mechanism for P. ussuriensis was established by integrating the physiological and molecular responses as well as plant growth under drought stress ( Figure S6). The findings clearly indicate that a series of potential droughtresponsive genes related to signal transduction, antioxidant, ion accumulation and osmotic regulation in the roots of P. ussuriensis were differentially expressed under drought stress. Combined with the physiological and biochemical indexes, we found the ion accumulation, the content of proline and soluble sugar, antioxidase activity (SOD, CAT, POD) and nonenzymatic antioxidant (GST) activity were all changed in P. ussuriensis root during PEGinduced drought stress. In summary, the combined analysis of temporal transcriptomic and physiological revealed complex genes and networks involved in the response of P. ussuriensis roots to different stages of drought stress.  Data Availability Statement: The raw sequence read from PacBio single-molecule real-time sequencing has been uploaded to the Sequence Read Archive database (SRA accession: PRJNA809518). The raw sequence read from Illumina RNA sequencing has been uploaded to the GEO database (GEO accession: GSE197096).

Conflicts of Interest:
There are no conflicts of interest to declare.