The genome and protein sequences of S. flavescens and Sophora moorcroftiana were obtained from previous studies (Qu et al., 2023; Yin et al., 2023). The HMM configuration file for WRKY domain (PF03106) was downloaded from the Pfam database (El-Gebali et al., 2019). Candidate WRKY members coding in the genomes of S. flavescens and S. moorcroftiana were identified using HMMER (v3.2.1) software, with an E-value threshold set at 10^-2 (Prakash et al., 2017). Only sequences containing the WRKY domain were considered as members of the WRKY family. To ensure the completeness of SfWRKY repertoire, we also examined the assembled novel transcripts obtained from the transcriptome assembly in section 2.5 for member identification. Furthermore, the WRKY genes in A. thaliana were derived from previous study (Abdullah-Zawawi et al., 2021).

The WRKY genes of S. flavescens were named according to their relative positions on the chromosomes. To investigate their protein properties, we utilized the ProtParam program (ExPASy tools, http://web.expasy.org/protparam/) to estimate the molecular weights (MWs) and theoretical isoelectric points (pI).

To construct the phylogenetic tree of the WRKY members, we utilized the ClustalW program in MEGA 7 (v7.0.26) software (Kumar et al., 2016) to perform multiple sequence alignments of the WRKY domain regions from above mentioned three species. Subsequently, we employed the Neighbor-Joining (NJ) method within MEGA 7 to build the phylogenetic tree, selecting Poisson model for amino acid substitution and applying pairwise deletion for gap treatment. To ensure the accuracy, we assessed the support for each relative branch through 1000 bootstrap replicates.

The DupGen_finder software (Qiao et al., 2019) was employed to conduct analysis of gene duplication patterns in S. flavescens. It includes whole genome duplication (WGD), tandem duplication (TD), transposed duplication (TRD), proximal duplication (PD), and dispersed duplication (DSD). The non synonymous substitution rate (Ka) and synonymous substitution rate (Ks) for TD gene pairs were calculated using KaKs_calculator3 (Zhang, 2022). The YN (Yang and Nielsen) model (Yang and Nielsen, 2000) was selected to compute the Ka/Ks ratio, which serves as an indicator of selective pressure on duplicated gene pairs. The gene density and intergenomic syntenic block analysis were conducted following the methods of Feng et al. (2024). Based on the genomic annotation information, TBtools software was utilized to obtain the upstream two kb genome sequence of SfWRKY genes from start codon. Subsequently, potential cis-regulatory elements were predicted and identified using default parameters via the PlantCARE website (Lescot, 2002).

To investigate the expression patterns of the SfWRKY genes in different tissues and growth stages, we collected four tissues (stem, flower, leaf, and root) of S. flavescens cultivated for over five years at Changzhi International Shennong Traditional Chinese Medicine Cultural Expo Park (Shangdang District, Changzhi City, Shanxi Province) on July 12, 2021. Additionally, we sampled pods at six different developmental stages between July 12 and August 6, 2021, with sampling conducted every five days. The roots of S. flavescens at eight distinct developmental stages, which were sown in September 2022, were collected on the 20th day of each month from April to November 2023 in Dianshang County (Lucheng District, Changzhi City, Shanxi Province). Each sample comprised three biological replicates collected in 50 mL centrifuge tubes and immediately frozen in liquid nitrogen before being stored at -80°C for further analysis. The experimental methods and analytical approaches for transcriptomics (stem, flower, leaf, and root; pod and root development) and metabolomics (pod and root development) were adapted from Zhong et al. (2024). The raw data from RNA-seq samples were archived in the NCBI database under accession number PRJNA1136989. The co-expression network of transcriptomes and metabolomes for pods and roots was constructed using R (v 4.2.2) with the WGCNA (v1.71) (Langfelder and Horvath, 2008) following the methods of Liu et al. (2024). The networks with inter-gene weight values greater than 0.3 were visualized using Cytoscape (v.3.8.2) (Otasek et al., 2019).

The roots of S. flavescens cultivated in Dianshang County were collected for real-time quantitative PCR (RT-qPCR) analysis. The plants included two different cultivation years, C and S represent the S. flavescens sowed in 2024 and 2022, respectively (sowing occurs every April), with samples collected in July 2024. For salt stress treatment, one-year-old seedlings of S. flavescens were irrigated with 250 mM NaCl solution (Salt), while the normal condition (NC) irrigated with equivalent distilled water. After treatment for 14 days, leaves were collected for RT-qPCR analysis. Each sample consisted of three biological replicates were collected into 50 mL centrifuge tubes and immediately frozen in liquid nitrogen before being stored at -80°C.

Total RNA was extracted from the tissues of S. flavescens using the polysaccharide polyphenol total RNA extraction kit (Beijing GeneBetter, China). The integrity of the RNA was confirmed with a Nanodrop 2000 spectrophotometer (Thermo Fisher, USA). Reverse transcription was performed using the HiScript II Q RT SuperMix kit (Takara, Takara Biomedical Technology (Beijing) Co., Ltd.), followed by qPCR utilizing the SYBR qPCR Master Mix kit (Vazyme, Nanjing, China). Specific primers were designed using Primer5 software (v5.00) ( Supplementary Table S1 ). EF-1α was utilized as the reference gene. All reactions were performed in triplicate, and the relative expression levels of genes were calculated using the 2^−△△Ct method. Statistical significance was assessed using Student’s t test.

We totally identified 69 SfWRKY genes based on genomic and transcriptomic data, which were named SfWRKY01-SfWRKY69 according to their positions on the chromosomes ( Figure 1 ; Supplementary Table S2 ). The 69 SfWRKY genes encoded proteins of varying sizes, ranging from the largest protein (SfWRKY12) with a molecular weight (MW) of 83.5 kD and composed of 757 amino acids, to the smallest protein (SfWRKY56) with an MW of 18.8 kD and containing 167 amino acids. The theoretical isoelectric point (pI) ranges from 4.8 (SfWRKY05) to 9.78 (SfWRKY18), indicating that different SfWRKY proteins perform various functions under different microenvironments. The 69 SfWRKY genes are unevenly distributed across nine chromosomes of S. flavescens ( Figure 1 ), with the majority located on chromosome 4 (Chr4, 15 genes) and Chr6 (12 genes), followed by Chr7 with 10 genes, and the least on Chr9 with only three genes.

Through correlation analysis between the WRKY gene numbers from the previous studies and their corresponding genome size, the former was not linearly correlated with the latter (Pearson’s correlation coefficient (r) = -0.0440, p > 0.05; Supplementary Figure S1 ). The number of SfWRKY genes was lower than those of A. thaliana (72), S. moorcroftiana (83), Medicago sativa (91), and Arachis hypogaea (158). However, it is relatively higher compared to other plants, such as Panicum miliaceum (32), Platycodon grandiflorus (42), and Dendrobioum catenatum (62). This indicates that there is no significant correlation between the size of a species’ genome and the size of the WRKY gene family.

To further understand the evolutionary diversity of the SfWRKY genes, we constructed an NJ tree for the WRKY family members of S. flavescens, S. moorcroftiana, and A. thaliana based on their WRKY domains ( Figure 2 ). Accordingly, the SfWRKY members were classified into three well defined groups. The 69 SfWRKY genes were unevenly distributed among the three groups, with 16 members in Group I, 44 members in Group II, and nine members in Group III ( Figure 2B ). The majority of S. flavescens and S. moorcroftiana WRKY members exhibit a one-to-one clustering pattern on the evolutionary tree, whereas the Arabidopsis WRKY members tend to form a cohesive cluster. Members in Group I have two WRKY domains located in N-terminal and C-terminal regions. Group II has largest member number, which is further subdivided into five subfamilies ( Figure 2C , IIa-IIe). Groups IIa and IIb tend to cluster into one branch, while groups IId and IIe tend to cluster together. In Group III, four closely located members (within 0.32 Mb) on the chromosome of S. moorcroftiana tend to cluster together, possibly originating from tandem duplication events. Among these groups or subfamilies, the member number in S. flavescens and S. moorcroftiana was similar ( Figures 2B, C ).

Subsequently, multiple sequence alignments were conducted for the WRKY domains of the three SfWRKY groups ( Figure 2D ). The results revealed that the WRKY domains were highly conserved across all categories, containing a heptapeptide domain (WRKYGQK) and a zinc-finger domain (C2H2 or C2HC). Additionally, there were six WRKY domain variants, five of which were WRKYGEK (SfWRKY56, and SfWRKY61). The other three variants included WRKYGKK (SfWRKY35 and SfWRKY47), WKKYAQT (SfWRKY29), and WRVKGQE (SfWRKY28). The gene structure and motif distribution analyses showed members from same group or subfamily showed similar characteristics ( Figure 3 ), which also provided evidence for the classification based on only phylogenetic analysis.

The DupGen_finder was used for analysis of SfWRKY gene duplication pattern, which included one pair of genes derived from WGD (such as SfWRKY04 and SfWRKY33), 8 pairs from TRD (such as SfWRKY20 and SfWRKY33), and 10 pairs of DSD (such as SfWRKY61 and SfWRKY20) events ( Figure 4A ). Among these TRD and DSD gene pairs, it is generally presenting one gene corresponding to multiple duplicate genes rather than one corresponding to one, which may be related to multiple rounds of duplication occurred in S. flavescens genome. For example, the TRD genes corresponding to SfWRKY33 are SfWRKY03, SfWRKY20, SfWRKY30, SfWRKY61, and SfWRKY68, and the DSD genes corresponding to SfWRKY53 are SfWRKY06, SfWRKY31, and SfWRKY48. We further calculated the Ka, Ks, and Ka/Ks values for the duplicated SfWRKY gene pairs, and the results showed that the Ka/Ks values ranged from 0.09 to 0.39 (< 1), indicating that these genes underwent strong purifying selection during evolution ( Table 1 ). This implies that gene duplications may contribute to the diversification and expansion of the SfWRKY gene family, while these duplicated genes are subject to strong functional constraints.

We also conducted a collinearity analysis among the three species, S. flavescens, S. moorcroftiana, and A. thaliana. The results showed that 16 SfWRKY genes from S. flavescens are present in collinear blocks with 11 counterparts from A. thaliana, and 24 SfWRKY genes from S. flavescens are present in collinear blocks with 12 counterparts from S. moorcroftiana ( Figure 4B ). It was observed that collinear genes tended to cluster into the same groups in phylogenetic tree ( Figure 2 ), implying that they have a common evolutionary origin.

To understand the potential roles of the SfWRKY family members in plant growth and development, response to plant hormones, and environmental stresses, we analyzed the distribution of cis-elements in the upstream promoter regions of the SfWRKY genes ( Figure 5 ). Totally, 18 plant growth and development related cis-elements were identified. The most abundant elements were the Box-4 (ATTAAT) and the G-box (TACGTG), which were involved in light responsiveness. They accounted for 30% and 20% of the total elements identified in this category. SfWRKY51 contains the highest number of such elements (13). Other elements were also observed, such as circadian control and tissue-specific motifs like the GT1-motif (photosynthetic reaction regulation), TCT-motif (light-responsive elements), and GATA-motif (plant development). Regarding the phytohormone responsive elements, the ABRE element involved in ABA responsiveness, the TCA-element involved in salicylic acid responsiveness, the TGACG-motif involved in MeJA-responsiveness and the ERE element involved in ethylene signal regulation were identified abundant in upstream regions of the SfWRKY genes. These four pattern elements represent more than 45% of the total hormone-responsive elements. In SfWRKY67, the number of ABRE elements is the highest (8), while in SfWRKY11, the number of ERE elements is the highest (7). In addition, other elements were also found, like the as-1 element involved in SA and oxidative stress responsiveness, and the CGTCA-motif involved in MeJA-responsiveness in this category.

In the abiotic and biotic stress category, different elements associated with stress responses, such as oxidation, defense, drought, wounding, heat, and low temperature, were observed. After our analysis in S. flavescens, the largest part of the elements belonging to the abiotic and biotic stress category corresponds to two general stress-responsive motifs, namely the MYB (CCAAT box) and MYC (CACATG box) binding sites, representing 30% and 20% of the total identified cis-elements, respectively. In addition, other stress-specific cis-elements were identified. Several of them were responsive to wounding and pathogens—including the WRKY-box (W-box), the TC-rich repeats, and the wound responsive motif (WUN-motif). Temperature-related elements, such as the stress responsive element (STRE) and the low-temperature responsive (LTR) motifs; drought related elements, such as the MYB BINDING SITE (MBS) and the Dehydration -responsive element (DRE)-core; and anaerobic conditions like the Anaerobic response element (ARE) motif ( Figure 5 ). Among them, SfWRKY37 contains the highest number of such elements, followed by SfWRKY60.

To investigate the expression patterns of SfWRKY family members among different tissues, we analyzed their expression levels in four different tissues, including leaves, flowers, roots, and stems, as well as different development stages of pods and roots. The results showed significant variations in the expression patterns of different SfWRKY genes across various tissues ( Figure 6 ). Notably, SfWRKY52 and SfWRKY66 exhibited consistently high expression levels throughout all examined tissues and developmental stages. We can classify the expression patterns of SfWRKY genes among different tissues mainly into three categories based on the tissue expression profiles: the first category includes genes with high expression levels, such as SfWRKY52, SfWRKY66, SfWRKY43, and SfWRKY52, which are highly expressed among the root, stem, leaf, and flower tissues, and across all the six developmental stages of the pod tissue. The second category includes low-expressed genes among many tissues, such as SfWRKY28, SfWRKY68, and SfWRKY47, among all the four tissues; SfWRKY12 and SfWRKY28 across the pod tissues; and SfWRKY12, SfWRKY28, and SfWRKY21 across the root tissues. The third category includes genes that do not express, such as SfWRKY61, SfWRKY63, and SfWRKY62, which are not expressed in any of the five tissues. It is noteworthy that the genes exhibiting high expression levels in one tissue or nearly ubiquitous high expression across all tissues, implying their predominant involvement in plant growth and development.

The gene coexpression networks were constructed based on transcriptome and metabolome data obtained from developmental stages of S. flavescens pods and roots. The results revealed that SfWRKY29 and SfWRKY41 were identified as hub-genes among the pod network ( Figure 6B ), indicating their pivotal roles in regulating the biosynthesis of secondary metabolites during pod development. IFR (2’-hydroxyisoflavanone reductase) and IFS (2-hydroxyisoflavanone synthase) emerged as key enzymes within the highly interconnected metabolic pathways. We further analyzed the coexpression network constructed based on the transcriptomics and metabolomics of S. flavescens root development, revealing SfWRKY29 as a central hub with the highest connectivity to other genes ( Figure 6C ), indicating its crucial involvement in regulating biosynthesis of secondary metabolites during root development. Notably, only I2’H (isoflavone 2’-hydroxylase) and I3’H (isoflavone 3’-hydroxylase) were identified as pivotal enzymes operating within this metabolic pathway.

As we know, the roots of S. flavescens can accumulate secondary metabolites such as flavonoids and alkaloids as they grow, and the accumulation increases with the growing years (Lei et al., 2021). Excessive accumulation often leads to autotoxicity, which inhibits plant growth or causes continuous monocropping obstacle. Given the outstanding performance of WRKY genes in resistance to adverse conditions, we used RT-qPCR technology to analyze the expression patterns of WRKY family members in the roots of S. flavescens with cultivated years. As shown in Figure 7 , among the 15 selected SfWRKY genes, all showed significantly higher expression in the SR (roots from S. flavescens sowed two years ago) than in the CR (roots from S. flavescens sowed in current years), with the greatest expression difference being observed in SfWRKY44 (upregulated by more than 335 times), followed by SfWRKY41 and SfWRKY39. Other SfWRKY genes also showed varying degrees of significant upregulation. This result implies that WRKY genes in S. flavescens play an important role in the response to the accumulation of secondary metabolites.

To investigate the expression patterns of WRKY family members in the leaf tissues of S. flavescens under salt stress, we used RT-qPCR to analyze the expression patterns of 9 SfWRKY genes ( Figure 8 ). The results showed that seven SfWRKY genes (SfWRKY18, SfWRKY20, SfWRKY35, SfWRKY52, SfWRKY55, SfWRKY66, and SfWRKY67) were downregulated under salt stress, while two genes (SfWRKY03 and SfWRKY50) were upregulated. The downregulation of SfWRKY55 was the most significant, while the SfWRKY03 was upregulated by 1.8 times when exposed to salt stress. Their significant expression differences under salt stress may suggest that they are involved in the response of S. flavescens to salt stress.