Identification of WRKY transcription factor family genes in Pinus massoniana Lamb. and their expression patterns and functions in response to drought stress

Background Pinus massoniana Lamb. is the timber species with the widest distribution and the largest afforestation area in China, providing a large amount of timber, turpentine and ecological products. Seasonal drought caused by climate warming severely constrains the quality and growth of P. massoniana forests. WRKY transcription factors play an important role in plant responses to abiotic stress. In this study, the molecular mechanisms by which P. massoniana responds to drought stress were analysed based on the P. massoniana WRKY (PmWRKY) family of genes. Results Forty-three PmWRKYs are divided into three major families, 7 sub-families, and the conserved motifs are essentially the same. Among these 43 PmWRKYs express under drought stress but with different expression patterns in response to stress. PmWRKYs respond to drought stress induced by exogenous hormones of SA, ABA, and MeJA. The expression of PmWRKY6, PmWRKY10, and PmWRKY30 up-regulate in different families and tissues under drought stress, while PmWRKY22 down-regulate. Transgenetic tobaccos of PmWRKY31 are with lower malondialdehyde (MDA) content and higher proline (Pro) content than wild type (WT) tobaccos. In transgenic tobaccos of PmWRKY31, expression levels of related genes significantly improve, and drought tolerance enhance. Conclusions This study analysed the molecular biological characteristics of PmWRKYs and investigated the expression patterns and functions of PmWRKYs in response to drought stress in P. massoniana. The results of this study provide a basis for in-depth research of the molecular functions of PmWRKYs in response to drought stress. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03802-7.

highly conserved domain composed of 60 amino acids that contains a heptapeptide structure with a highly conserved WRKYGQK domain at its N-terminus, which is related to DNA binding activity [8]. Although the heptapeptide structure is highly conserved, WRKYGQK still has variations, such as WRKYGKK, WKKYGQK, and WRKYGRK [7]. In addition, gene function changes with the variation in conserved structure [9]. The C2H2 (CX 4-5 CX 22-23 HXH) and C2HC (CX 7 CX 23 HXC) zinc finger structures at the C-terminus participate in proteinprotein interactions (PPIs) and assist in DNA binding [8,10]. The WRKY genes can bind to the cis-acting element of the downstream gene promoter, which contains a very conserved (C/T)TGAC(T/C) sequence called the W-box [11]. WRKY genes bind to this element, thereby activating downstream gene expression. In addition, WRKY genes can also function by combining with elements such as W-box, PRE4, PRE2, G-box, WK box and SURE [12][13][14].
Drought is an important abiotic stress factor that hinders plant growth and development and is a frequent cause of damage. Plants can regulate their resistance to drought stress through the WRKY regulatory network and hormone signal transduction. For example, HbWRKY1 expression in Hevea brasiliensis [15] and HrWRKY21 expression in Hippophae rhamnoides [16] are upregulated in response to drought stress. GhWRKY41 and ThWRKY2 enhance the drought tolerance of cotton and Tamarix hispida [17,18]. EjWRKY17 from Eriobotrya japonica enhances drought tolerance in transgenic Arabidopsis [19]. WRKY genes can also improve the drought tolerance of plants by increasing osmotic regulatory substances, such as soluble sugar and proline (Pro) and other contents, and reducing the content of malondialdehyde (MDA) [20,21].
Pinus massoniana Lamb., a major plantation timber species in southern China, has excellent characteristics, such as barren tolerance, strong adaptability, fast growth, and a wide range of applications [22,23]. However, the uneven rainfall distribution and seasonal drought in southern China [24], severely restrict the growth of P. massoniana. Global warming and environmental deterioration further exacerbate this problem. Therefore, research on the drought tolerance of P. massoniana is of great importance. There are relatively few studies on the drought tolerance mechanisms of P. massoniana at the molecular level. Whole-genome sequencing of the P. massoniana genome has not yet been carried out due to its large size (expected size of approximately 24 Gb, gigabase pairs, chromosome-level assembly). Therefore, the identification of P. massoniana WRKY (PmWRKY) based on transcriptome analysis is the top choice. In this study, we identified the genes in the PmWRKY family in P. massoniana families with different drought tolerance based on transcriptome data and analysed their physicochemical properties, phylogenetic evolution, and conserved motifs. We also analysed the expression patterns of 4 PmWRKYs under hormone treatment and drought stress and the drought tolerance function of PmWRKY31 to provide a reference for the further in-depth study of the functions of PmWRKYs and their mechanisms of response to drought stress.

Identification and physiochemical characteristics of PmWRKYs
Based on previous full-length transcriptome data of the second and third generations, the gene sequences annotated WRKY were obtained. A total of 496 sequences were obtained using HMMER software and the BLAST method. After removing redundant sequences, 281 candidate sequences were obtained. The conserved WRKY domains of candidate protein sequences were analysed using NCBI CD search, SMART, and Pfam. After the sequences that contain incomplete WRKY domains or do not contain WRKY domains were removed, 43 PmWRKYs were finally identified and named based on homologous sequence data (Supplementary Table 1).
Analysis of the physicochemical properties and functional structure of proteins revealed that the PmWRKY proteins encode 120 (PmWRKY35) to 908 (PmWRKY1) amino acids; the molecular weights of the PmWRKY proteins range from 13.48 (PmWRKY35) to 98.15 (PmWRKY1) kDa. The theoretical isoelectric point (pI) values for PmWRKYs range from 4.68 (PmWRKY10) to 10.23 (PmWRKY40), with an average pI of 7.93. Among the PmWRKYs, 15 have a pI < 7, which is relatively acidic, and 28 have a pI > 7, which is alkaline. The instability coefficient of PmWRKY40 protein is 36.07, and the instability coefficients of 42 PmWRKY proteins are all greater than 40, indicating that they are unstable proteins. Among the 43 PmWRKY proteins, only PmWRKY25 (-0.459), PmWRKY37 (-0.320), and PmWRKY39 (-0.307) have an average hydropathicity > -0.5, indicating that they are hydrophilic; the rest are hydrophobic. The subcellular localization of all proteins was predicted to be in the nucleus, with no transmembrane structure (Supplementary Table 1).

Analysis of the phylogenetic, conserved motifs and multiple sequence alignments of PmWRKYs
The phylogenetic tree showed that the WRKY family can be classified into families I, II and III. Family I, which contains 2 WRKY domains and a C2H2 zinc finger structure, is further categorized into the IC and IN subfamilies. Family II, which contains a WRKY domain and a C2H2 zinc finger structure, is categorized into subfamilies IIa, IIb, IIc, IId and IIe. Family III contains a WRKY domain and a C2HC zinc finger structure [25,26]. Families I, II and III include 13 PmWRKYs, 29 PmWRKYs,  and 1 PmWRKY (PmWRKY6), respectively. The IIa, IIb,  IIc, IId and IIe subfamilies includes 10, 4, 6, 6, and 3 PmWRKYs, respectively. In addition, P. massoniana has a close genetic relationship with Pinus taeda and Picea glauca (both in the Gymnosperm Pinaceae) (Fig. 1).
Analyses of conserved gene motifs revealed 10 conserved motifs with a length of 15-57 amino acids. Motif1 and Motif3 are conserved heptapeptide domains of WRKY. Motif1 and Motif2 distribute throughout all subfamilies, while Motif3 only distributes in subfamily I. Motif4, Motif5, and Motif6 distribute in subfamilies IIa and IIb. Motif7 distributes in family I and subfamilies IIa, IIb and IIc. Motif8 and Motif10 distribute in family I and subfamily IIa. Motif9 distributes in family I and subfamily IIe. The evolution of WRKY genes is highly conserved, and they may have similar functions. However, there are also inconsistent conserved motifs contained in the same protein subfamily; for example, PmWRKY10 in subfamily IIe contains more Motif9 than does other members (Fig. 2a, b). Multiple sequence alignment of PmWRKY protein sequences revealed that the protein sequences of the same subfamily genes have high similarity, further indicating that they may have similar gene functions. Most of the PmWRKY protein sequences contain the conserved heptapeptide WRKYGQK at the N-terminus, and only 5 proteins in subfamilies IN, IId, and IIc have amino acid variations, i.e., WTKYGKR, WRKYGQN, WRKYGKK, and WRKYGRK. There is a variable zinc finger structure at the C-terminus of the sequence, such as the deletion of the zinc finger structure in the C-terminus conserved domain of PmWRKY36 and the presence of zinc finger structure variation in PmWRKY39 and PmWRKY43 proteins ( Supplementary Fig. 1).

Protein-protein interactions of PmWRKYs
The PPI network of PmWRKYs (Fig. 3) shows that most PmWRKYs can interact with more than 1 protein, e.g., PPIs between PmWRKY40 and PmWRKY30 and PPIs between PmWRKY6 and PmWRKY20. PmWRKY40 (AtWRKY33) plays an important role in the entire PPI network. AtWRKY33 is involved in the response to various abiotic stresses. Hence, PmWRKY40 may have similar functions [27].

Expression analysis of PmWRKYs by RNA-seq
An expression heatmap of 43 PmWRKYs under drought conditions was generated based on the transcriptome data ( Fig. 4). All the genes were expressed during drought stress. The expression of 4 genes peaked on the 7th day of drought (D1), expression of 4 genes peaked on the 7th day of drought (D2), that of 15 genes peaked on the 8th day of drought (D3), and that of the remaining genes all peaked on C1-C3. With FDR ≤ 0.001 and |log2FC|≥ 2 as significant difference gene screening standards, at the same time combined with the phylogenetic tree, we eventually chose four differentially expressed genes, PmWRKY6, PmWRKY10, PmWRKY22, and PmWRKY30, for research of expression patterns under drought stress after treatments of hormones.

Analysis of the expression of PmWRKYs under hormone treatments and drought stress
In needles, PmWRKYs had same expression pattern in different drought-resistant families, but also had different expression patterns. For example, expression of PmWRKY6 was earlier induced by MeJA and ABA in DT families, while it was induced by MD stress in DS families. The expression of PmWRKY6 was induced by SA, MeJA and ABA in DT family in response to drought stress. PmWRKY10 and PmWRKY22 showed higher expression levels in DT family than in DS family Fig. 3 Prediction of the protein-protein interaction (PPI) network between PmWRKYs and Arabidopsis proteins. The PPI network of PmWRKYs was constructed using STRING and Cytoscape software and was analysed and predicted based on the WRKY proteins of Arabidopsis at different stress stages. However, after re-watering, the expression levels of PmWRKY10 and PmWRKY22 in DS families were significantly higher than those in DT families. Combined with the phenotypic data, it was suggested that DS families need rapid up-regulation of these genes to cope with drought stress, while DT families showed more durable drought resistance. Under SD stress, the expression of PmWRKY30 increased sharply in DT families, but its expression under hormone treatments was lower than that of CK2, while its expression of hormone treatments was higher than that of CK2 in DS families (Fig. 5).
In stems, compared with CK1, expression patterns of 4 PmWRKYs were basically same in the two families under drought stress. Expression of PmWRKY6, PmWRKY10 and PmWRKY30 was up-regulated, while the expression of PmWRKY22 was down-regulated. Expression of PmWRKY6 in stems of DT was much higher than that in DS stems. Compared with CK2, there were different expression patterns for these 4 genes in different drought resistant families. Expression of PmWRKY6 was induced by SA, MeJA and ABA in DT, and expression of 4 genes was significantly induced by ABA in DS. After re-watering, compared with CK1, only expression of PmWRKY30 was down-regulated in drought-resistant families CK2, while the other treatments showed that expression levels of 4 genes in CK2 were significantly higher than those in CK1 and higher than those in hormone treatments (Fig. 6).
In roots, compared with CK1, expression of PmWRKY30 was up-regulated in the two families. Expression levels of PmWRKY6, PmWRKY10 and PmWRKY22 were higher than those of CK1 during light drought stress in DS, while all three genes were down-regulated in DT. Compared with CK2, PmWRKY6 was up-regulated by MeJA and ABA in DS, and was Heatmap of PmWRKY expression using RNA-seq data. Blue represents low expression levels, and red represents high expression levels. C and D represent the control group (normal growth, C) and the treatment group (drought stress, D), respectively; 1, 2, and 3 represent the control and drought treatments on the 7th day, 7 h after rewatering on the 8th day, and the 8th day, respectively. TBtools was used to generate the gene expression heatmap up-regulated by SA and MeJA in DT. PmWRKY10 and PmWRKY22 were up-regulated by MeJA in DS. In DT, PmWRKY10 was up-regulated by MeJA and ABA, and PmWRKY22 was up-regulated by ABA. PmWRKY30 was up-regulated by three hormones in DT, and only up-regulated by ABA in DS. After re-watering, PmWRKY6 and PmWRKY30 were up-regulated by ABA in drought-sensitive families, while these two genes were induced by SA and MeJA in drought-resistant families (Fig. 7).
These 4 PmWRKYs were expressed in roots, stems, and leaves, according to results of their expression in response to hormone and drought stress in each tissue of the two families. Different PmWRKYs responded differently to hormone signals. There might be different hormone pathways in response to drought stress in different drought-resistant families. It was speculated that the expression of PmWRKY6, PmWRKY10 and PmWRKY30 could be involved in response to drought stress.

Phylogenetic analysis of PmWRKY31 gene information and overexpression in tobacco under drought stress
Our research group previously analysed the insect resistance mechanism of PmWRKY31 (PmWRKY30 in this paper) [28] and found, through transcriptome data, that PmWRKY31 was upregulated under drought conditions. Therefore, the drought tolerance function of PmWRKY31 was verified by overexpressing PmWRKY31 in tobacco (Fig. 8). After the simulated drought treatment, the phenotypic changes in water loss and wilting were obvious in both transgenic tobacco and WT tobacco (Fig. 8a). Under normal culture conditions, there was no difference in the growth status of transgenic tobacco and wild-type (WT) tobacco (Fig. 8b). On the 2nd day of the stress treatment, both transgenic and WT tobacco plants exhibited significant changes. The leaves of the WT tobacco plants were severely wilted, and the upper leaves of the transgenic tobacco plants of the W and W3 lines grew normally. On the 4th day of stress treatment, WT tobacco plants were severely wilted, with significant water loss in the stems, while the transgenic plants only exhibited wilting leaves. The results indicate that PmWRKY31 improved the drought tolerance of tobacco plants (Fig. 8a).
The phylogenetic tree showed that P. massoniana and Pinus kesiya had the closest genetic relationship. The similarity between PmWRKY31 and the amino acid sequences of other plants was between 41.54% and 99.53%, and the homology between PmWRKY31 and Pinus kesiya was highest, which is consistent with the phylogenetic tree, indicating significant differences in the protein structures of the same gene family between gymnosperms and angiosperms (Fig. 8c).

Changes in MDA and Pro contents in tobacco plants overexpressing PmWRKY31 under drought stress
Before treatment (0 d), the MDA content in the W3 line was the highest. During the drought stress process, the MDA content of the WT was significantly increased and was significantly higher than that of the transgenic tobacco (P < 0.05). This result indicates that PmWRKY31 further increased the drought tolerance of tobacco by reducing MDA accumulation (Fig. 8d). The Pro content in the transgenic tobacco was higher than that in WT tobacco during the drought stress process. The results indicated that under drought stress, PmWRKY31 improved drought tolerance by promoting Pro accumulation (Fig. 8d).

Analysis of the gene expression patterns in tobacco plants overexpressing PmWRKY31 under drought stress
During drought stress, no PmWRKY31 expression was detected in WT tobacco (Fig. 9a). Under drought stress, PmWRKY31 expression first increased and then decreased in the W and W3 lines. PmWRKY31 expression in the W3 line was higher than that in the W line during the stress process (Fig. 9a).
The expression of the NtAPX, NtCBL and NtCAT genes in tobacco was upregulated overall, but the expression patterns were different. The expression of NtAPX was upregulated on the 4th day of treatment, the expression of NtCBL was upregulated during the entire stress process, and the expression of NtCAT decreased after the initial increase but was still higher than that of WT tobacco (Fig. 9b). These results indicated that PmWRKY31 might enhance the drought tolerance of transgenic tobacco by regulating downstream drought tolerance-related genes.
In addition, under drought stress, the expression of NtWRKYs in WT tobacco, except for 2 genes, decreased on the 4th day of treatment, and the expression of other genes was upregulated during the entire stress process and was higher than their expression in W and W3. These results indicate that PmWRKY31 inhibits the expression of NtWRKYs in tobacco (Fig. 9c).

Discussion
With the development of omics technology, more species of TF families have been identified. A large number of studies have shown that the WRKY TF family is involved in plant growth and development and plant response to abiotic stress [29,30]. In this study, 43 PmWRKYs were identified based on transcriptome data. The number of identified PmWRKYs was lower than the 72 WRKYs in A. thaliana [5], 104 WRKYs in P. Trichocarpa [7], 61 WRKYs in Taxus chinensis (also a gymnosperm) [31] but was higher than the reported 25 [23] or 31 WRKYs [32] in P. massoniana, indicating that this study is more comprehensive. This study revealed that PmWRKYs had 4 types of variations in the conserved heptapeptide structure. The variations included structural variations, which may cause the genes to be unable to specifically bind to downstream genes, thereby losing gene function. For example, in Crystalline soybean, GmWRKY6 and GmWRKY21 with variations in the heptapeptide structure cannot bind to the downstream gene W-box [33]. In addition, zinc finger deletion (PmWRKY36) and mutation (PmWRKY39 and PmWRKY43) occurred in PmWRKY proteins, but whether the function of the zinc finger structure (such as PPI and assisting DNA binding) is lost in these 3 genes still needs to be further studied [8]. Conserved motif analysis showed that the conserved motifs of genes in the same subfamily were basically the same, suggesting that genes in the same subfamily might have similar functions [25]. In this study, the PmWRKY gene family was classified into families I, II and III. Family II was categorized into 5 subfamilies. The phylogenetic tree shows that subfamilies IIa and IIb have a close genetic relationship, that subfamilies IId and IIe are clustered on the same branch, and that subfamily IIc is far from the above 4 subfamilies and has evolutionary characteristics similar to those of family II of A. thaliana [34]. Furthermore, the subfamilies in P. massoniana had closer relationships with gymnosperms, indicating that the evolutionary relationship between gymnosperms was close. In summary, the study results indicate that the WRKY gene family of P. massoniana is highly conserved with regard to gene structure and evolutionary level.
In this study, under drought stress and hormone induction, the 4 PmWRKYs were constitutively expressed under drought treatment, PmWRKY6, PmWRKY10 and PmWRKY30 positively regulated drought stress, and PmWRKY22 negatively regulated drought stress. It was reported that MaWRKYs enhance the cold tolerance of bananas by combining the W-box elements of MaNCED1 and MaNCED2 [35]. Whether each PmWRKY acts alone or multiple PmWRKYs act together in response to drought stress needs to be further investigated. Combined with the expression results for PmWRKYs in families with different drought tolerance, it is speculated that PmWRKY6 and PmWRKY30 are the key genes that enhance the response of the drought-tolerant family to drought.
Hormones play a vital role in the growth and development and abiotic stress response of plants. For example, SA, ABA, and JA signals can independently interact with each other to improve the drought tolerance of plants [36,37]. For example, TFs can improve plant stress resistance by participating in SA signalling, AtWRKY39 overexpression in Arabidopsis can improve its heat tolerance, and AtWRKY39 expression is induced by SA and MeJA to improve plant heat tolerance [38]. ZmWRKY79 enhances drought tolerance in plants by regulating the synthesis of ABA [39]. The above findings indicate that TFs can play a role by participating in hormone synthesis or by participating in the hormone signal transduction network. In this study, under drought stress, PmWRKYs responded to exogenous SA, ABA, and MeJA, and each PmWRKY responded to at least 1 hormone, indicating that PmWRKYs may improve drought tolerance in P. massoniana. JrWRKY6 and JrWRKY53 in Juglans regia may enhance its resistance to salinity, drought and heat stress through the ABA signalling pathway and W-box binding activity [39]. CsWRKY2 in Camellia sinensis is involved in the ABA downstream signalling pathway and improvements in drought tolerance [40]. In this study, PmWRKY6 and PmWRKY30 expression was upregulated by treatment with 3 hormones. It is speculated that PmWRKY6 and PmWRKY30 may be involved in a hormone signalling pathway and improve the drought tolerance of P. massoniana. However, the specific mechanisms need to studied further.
Interestingly, although the expression of PmWRKY22 was induced by hormones, its expression was lower than that of CK1. Studies have shown that AtWRKY18, AtWRKY40, and AtWRKY60 are negative regulators of ABA signalling and that these genes, as the upstream genes of the ABA signalling pathway, participate in the ABA signalling regulatory network by inhibiting ABA signalling pathway-related genes such as MYB2 and ABI5 [41,42]. Therefore, it is speculated that PmWRKY22 may negatively regulate hormonal signals in response to drought. The expression of PmWRKYs induced by hormone signals was different in the 2 families. This result was similar to the expression pattern of CsPRXs in 2 varieties of C. sinensis [43], results that may be caused by the different genetic backgrounds of the 2 families with different drought tolerance. In addition, it is reasonable to believe that the high expression of PmWRKYs increases the drought tolerance of the drought-tolerant family of P. massoniana.
Currently, in plants with no established genetic transformation system, the effects of genes on plant phenotypes, physiology, and related genes are analysed through functional validation in model plants to explore the functions of genes in improving plant resistance [30,44]. In this study, PmWRKY31-overexpressing transgenic tobacco was more drought tolerant than WT tobacco; its MDA content was significantly lower than that of WT tobacco, and its Pro accumulation was higher than that of WT tobacco. The results indicate that transgenic tobacco could better alleviate membrane damage and maintain osmotic regulation and thereby has better drought tolerance [17,20].
TFs mainly regulate downstream genes, activate or inhibit the expression of downstream genes, and interact with drought-related genes to resist drought stress [17,45]. Pbr-WRKY53 improves drought tolerance of tobacco by increasing the expression of drought tolerance-related genes such as NtSOD and NtCAT [46]. In this study, the expression of NtCBL, NtCAT , and NtAPX in transgenic tobacco was upregulated, indicating that PmWRKY31 might increase the drought tolerance of tobacco by increasing the expression of drought stress-related genes. Interestingly, PmWRKY31 inhibited the expression of NtWRKYs in PmWRKY31overexpressing transgenic tobacco to improve drought tolerance, but the specific regulatory mechanism needs to be further studied.

Conclusions
In this study, we successfully identified 43 PmWRKYs and analysed the physicochemical properties, phylogenetic evolution, and conserved motifs of PmWRKY genes family members. We found that the genes of this family have specific domains with highly conserved evolutionary levels. According to transcriptome data and gene expression analysis, PmWRKYs respond to drought stress and phytohormone treatments at different degrees and stress stage. It is also found that PmWRKYs in response to drought stress and phytohormones might have the same or different regulatory mechanisms. Furthermore, the transformation tobaccos of PmWRKY31 show that by regulating the expression of NtCBL, NtCAT and NtAPX, MDA content decreases, Pro content increases, and drought resistance is improved by alleviating membrane system damage and maintaining osmoregulation. These will contribute to further study on molecular mechanism of PmWRKYs in response to drought in P. massoniana.

Plant material and drought stress treatments
The drought-tolerant family 19-220 (DT) and the droughtsensitive family 19-214 (DS) of P. massoniana were provided by the Timber Forest Research Institute, Guangxi Forestry Research Institute (GFRI) Nanning, China. Two germplasm were selected from 132 excellent families by treatments of simulating natural drought stress, and then comprehensive evaluation of wilting degree, survival rate, seedling height, ground diameter and other indicators. We had got the permission to collect and use the two germplasm lines which were stored in GFRI. The experimental research and field studies on plants comply with relevant institutional, national, and international guidelines and legislation. Half-year-old seedlings of the 2 families with a good and basically consistent growth status were selected and placed in the drought shed of the nursery. Based on the results of previous experiments (Supplementary Table 2, Supplementary Fig. 2), 50 mg/L SA, 0.5 mmol/L MeJA, and 25 mg/L ABA were used for spraying treatments. Five groups were set up, i.e., CK1, normal watering; CK2, drought stress; 50 mg/L SA + drought stress; 0.5 mmol/L MeJA + drought stress; 25 mg/L ABA + drought stress, with 3 replicates for each treatment and 10 plants for each replicate. After 3 days of continuous watering, the above concentrations of SA, MeJA, ABA and distilled water (CK1 and CK2) were sprayed on the plants, with 10 mL sprayed on the aboveground part (stems and needles) and 10 mL sprayed on the underground part (roots). After spraying at 8 am every day for 4 consecutive days, all experimental seedlings were supplemented with 20 mL of water, which was recorded as the 0th day of drought stress, and the seedlings were subjected to continuous natural drought stress. CK1 was the normal watering control group, and the seedlings in this group were watered once every 2 days. The needles (2 cm from the top of the plant), stems (3.5 cm from the tip of the stem), and roots (3 cm from the tips of the main root and lateral roots) of P. massoniana were collected at light drought stress (LD, 55-70%, soil moisture content, MC %), moderate drought stress (MD, 45-55%), severe drought stress (SD, 30-45%) [47], and 48 h after re-watering (RW) and stored in a -80 °C freezer. The classification of drought grade and the determination of soil moisture content during sampling were determined according to the treatment phenotype of plants, as shown in Supplementary Fig. 2.
WT tobacco (Nicotiana benthamiana) and the F4 generation lines W and W3 of PmWRKY31-overexpressing transgenic tobacco were used as experimental materials. Seeds of WT tobacco were provided by Gold Tech Biotechnology Co., LTD (Nanning China). WT tobacco leaves were used as material to induce sterile tobacco seedlings, and pBI121_PmWRKY31 was constructed for Agrobacterium transformation. Hygromycin was used for prescreening. DNA of seedlings from transgenic tobaccos was extracted by CTAB method from leaves. According the sequence of PmWRKY31, PCR primers were designed for specific detection (PmWRKY31F:CCC TCT AGA ATG GAA GCA GTG GGG TTGAG, PmWRKY31R: TTT CCC GG GTC ACT GGA GCA ATT TAG CC). Transgenic Tobaccos with PmWRKY31 was cultured to F4 generation for further experiment [28]. They were cultured for 6 weeks in an artificial climate chamber (25 °C during the daytime/28 °C at night; relative humidity (RH) = 60%; light/ dark cycle = 16 h/8 h). Subsequently, polyethylene glycol (PEG) was used to simulate drought stress in both transgenic tobacco and WT tobacco. The specific method was as follows: the roots of each plant were watered with 15 mL of 20% PEG6000, and 1/3 of the 20% PEG6000 was poured into the bottom of the saucer tray to maintain the osmotic state [48]. The phenotypic differences in tobacco were observed. The leaves of transgenic tobacco and WT tobacco were sampled at the same position on the 0th, 2nd, and 4th days of drought stress to assess physiological indicators and the expression of drought tolerance-related genes.

Identification and physicochemical properties analysis of WRKY family members of P. massoniana
PmWRKY data were obtained from the third-generation full-length transcriptome and insect-resistant transcriptome of P. massoniana sequenced by our research group [49], a transcriptome for lateral branch differentiation [50], transcriptome of Aluminum stress [51], and a drought-tolerant, low-temperature stress, transcriptome (unpublished). Illumina HiSeq 4000 high-throughput sequencer was used for sequencing. The sequencing read length was 150 bp (Sequencing was performed by Novogene Bioinformatics Technology CO., LTD.). Clean Data of all samples was assembled by Trinity software [52]. Unigene sequences were compared with NR, SwisS-ProT, GO, COG, KOG and KEGG databases to obtain the annotation information. Bowtie method [53] and RSEM were used to estimate the expression levels of genes. Reference genomes were Norway spruce (Picea abies (L.) Karst.), loblolly pine (Pinus taeda L.) and white spruce (Picea glauca) as the main species. TransDecoder software was used to predict the coding DNA sequence (CDS). TBtools [54] and the ORF finder (https:// www. ncbi. nlm. nih. gov/ orffi nder) of the National Center for Biotechnology Information (NCBI) were used to predict the longest open reading frame (ORF). The 3 software programs were used together to obtain the nucleic acid sequences and amino acid sequences. Functional annotation was performed using the NR, SwissProt, and Pfam databases. Subsequently, the union of the 3 sets of annotation results was used to extract the nucleic acid and protein sequences using the Perl language. A hidden Markov model of the WRKY gene domain (PF03106) was downloaded from the Pfam database (https:// pfam. xfam. org/) [55]. WRKY genes were identified using HMMER software. The amino acid sequences of Arabidopsis WRKY genes downloaded from PlantTFDB (http:// plant tfdb. gao-lab. org/ index_ ext. php) [56] were compared with the amino acid sequences of P. massoniana WRKY genes to identify protein sequences. CD-HIT (http:// weizh onglab. ucsd. edu/ webMGA/ server/ cd-hit/) was used for clustering analysis, and redundant sequences were removed. NCBI CD Search (https:// www. ncbi. nlm. nih. gov/ Struc ture/ cdd/ wrpsb. cgi) [57], SMART (http:// smart. embl. de/) [58] and Pfam were used to detect the WRKY domains of the candidate proteins. Last, PmWRKYs of P. massoniana were identified by removing the sequences that contained incomplete domains or did not contain domains.

Phylogenetic analysis, conserved motif analysis and multiple sequence alignment of PmWRKYs
The WRKY gene protein sequences of Arabidopsis (72), Pinus taeda (20), and Picea glauca (16) were downloaded from PlantTFDB. MEGA 7 software [62] was used to perform multiple sequence alignment of the WRKY protein sequences of P. massoniana, A. thaliana, P. taeda, and P. glauca. The neighbour-joining (NJ) method was used to reconstruct phylogenetic trees using the P-distance model through 1000 bootstrap iterations. Subsequently, iTOL [63] was used to improve the phylogenetic trees. MEME (https:// meme-suite. org/ meme/) [64] was used to analyse the conserved motifs of the WRKY gene of P. massoniana with the following parameter settings: number of predicted motifs = 10, and length of the motifs = 6-60 AA. ClustalW software was used to perform multiple sequence alignments of WRKY proteins in each subfamily, and GeneDoc software was used to generate multiple sequence alignment diagrams.

Analysis of protein interaction among PmWRKYs
The STRING website (https:// string-db. org/ cgi) [65] and Cytoscape software [66] were used to predict the potential PPI network and biological functions of PmWRKY proteins in P. massoniana based on the WRKY protein analysis of A. thaliana.

RNA-Seq data analysis of PmWRKYs
An expression heatmap of PmWRKYs under drought stress was generated based on previous drought transcriptome data. Transcriptome data were sequenced by double-ended sequencing (PE150) using an Illumina HiSeq 4000 sequencer, and the resulting clean reads were assembled using Trinity software to obtain spliced transcripts (Unigene) for subsequent analysis. A treatment group (D) and control group (C) were set up. The treatment group was under continuous natural drought stress, and the control group was normally watered. The roots of the seedlings of the treatment group and the control group were collected on the 7th day (D1), 7 h after rewatering on the 7th day (D2), and 8th day (D3) of drought stress treatmen. For detailed processing and sampling, see Supplementary Fig. 3. Using the obtained transcriptome data, TBtools was used to plot the expression heatmap for the PmWRKY gene family under drought stress.

RNA extraction, cDNA synthesis, and RT-qPCR
RNA was extracted using a polyphenol polysaccharide plant RNA extraction kit from Tiangen (Beijing, China). M-MLV reverse transcriptase (Takara, Shanghai, China) was used for reverse transcription following the specific steps in the instruction manual, reverse transcription primer was Oligo(dT) 18 :5'-GGC CAC GCG TCG ACTAG TAC(T) 18 -3' . The concentration of all cDNA samples was adjusted to 50 ng/μL. Primer Premier 5.0 was used to design fluorescent quantitative PCR primers (Supplementary Table 3) with PmCYP as the internal reference gene [67], Nine N. benthamiana gene sequences were downloaded from NCBI and Adachi et al. [68]. A TB Green ® Premix Ex Taq ™ II (Tli RNaseH Plus) fluorescence quantification kit (Takara, Shanghai, China) was used to configure the reaction system. PCR analysis was performed using a quantitative PCR system (CFX96, Bio-Rad, California, USA). The measurement was biological replicated 3 times for each sample, and the relative expression of each gene was calculated using the 2 −ΔΔCt method [69].

Data analysis
One-way analysis of variance (ANOVA) in SPSS software was used to analyse the significance of the data. TBtools was used to generate the expression heatmap, and Graph-Pad Prism 8 was used to plot histograms.