Fruits of two bayberry cultivars (‘Dongkui’ and ‘Shuijing’) were selected as experimental materials, collected from LinGe Family Ecological Farm in Lin’an District, Hangzhou City, Zhejiang Province. Fruit samples were collected at five different developmental stages (designated S1-S5), corresponding to 51, 58, 65, 72, and 80 days after flowering, respectively. The experiment was conducted with three biological replicates, each containing at least five fruits. All collected samples were immediately frozen in liquid nitrogen and stored at -80°C for subsequent analysis.

To screen the WRKY genes of bayberry, the genome sequences and annotation information of bayberry was downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_003952965.2/). The WRKY sequences of Arabidopsis were obtained from The Arabidopsis Information Resource (TAIR) (https://www.arabidopsis.org/). The hidden Markov model (HMM) for the WRKY conserved domain (PF03106) was retrieved from the Pfam protein family database. Using this HMM model and the Arabidopsis WRKY protein sequences, candidate WRKY genes in bayberry was identified through HMMER and BLASTP tools. Initially, the Arabidopsis sequences were used as query sequences to perform a BLASTP search against the bayberry proteome (e-value < 1e-10, identity >= 40). The WRKY proteins in the bayberry genome was identified using the HMM search function in TBtools, and redundant sequences were removed from the results. Molecular weight (Da), isoelectric points, and other physicochemical properties of WRKY proteins were determined using the online ExPASy tool (http://web.expasy.org/protparam/).

A Clustal X2 sequence alignment was performed on the full-length protein sequences of Arabidopsis and bayberry WRKY proteins. The phylogenetic tree was constructed using the MEGA X software with the Neighbor-Joining method, and bootstrap resampling was set to 1000 repetitions. The phylogenetic tree was subsequently polished using the ITOL online tool (https://itol.embl.de/).

The exon and intron information of bayberry WRKY genes was retrieved from the bayberry genome annotation data. The gene structure diagram was generated using the Visualiza Gene Structure function in TBtools. For motif analysis, the MEME suite (http://meme-suite.org/tools/meme) was employed to identify six conserved motifs within the bayberry WRKY protein sequences. The distribution of these motifs was visualized using the Gene Structure View function in TBtools.

The chromosome positions and tandem repeats of MrWRKY genes were visualized through the BLAST and MCScanX functions in TBtools software based on the MrWRKY protein sequences. Based on the bayberry genome and annotation files, whole-genome sequence alignment was conducted through the MCScanX and gene position extract functions in TBtools, yielding syntenic relationship files of MrWRKYs within the bayberry species. The chromosome length file was obtained using the fasta stats function, while the gene density file was generated via the fasta stats table function. Finally, by integrating the syntenic relationship file, chromosome length file, and gene density file, the homologous gene diagram for MrWRKYs were constructed using the Advanced Circos function zone of TBtools software.

The cis-acting elements in the 2000 bp promoter regions of each bayberry WRKY gene family members were analyzed using the PlantCARE website (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/). And the analysis results were visualized using the TBtools software.

Total RNA was extracted from various cultivars of bayberry during different fruit development stages using the CTAB method (Liu et al., 2024). RNA quantity and quality were determined using a NanoDrop One spectrophotometer, and RNA integrity was assessed by 1% agarose gel electrophoresis. High-quality cDNA libraries were constructed, and gene expression analysis was performed. A transcriptome-wide expression profile of bayberry WRKY genes was generated using strand-specific RNA-seq, with transcript abundance calculated as reads per kilobase of transcript per million mapped reads (RPKM). Sequencing data used in the current study are available in the NCBI Sequence Read Archive database (project number: PRJNA1105392). Expression patterns of bayberry WRKY genes were visualized as a heatmap using the Heatmap illustrator tool in TBtools. The relative expression levels of genes were calculated using the 2^−ΔΔCt method, with Mr-actin serving as the reference gene for normalization. All experiments were conducted with three biological and technical replicates to ensure reproducibility.

The promoter sequence was inserted upstream of the firefly luciferase (LUC) reporter gene vector, while the Renilla luciferase (REN) gene, driven by the 35S promoter, served as an internal control on the same vector. The coding sequence (CDS) of MrWRKY was cloned into the pGreenII 0029 62-SK vector as the effector gene, with the empty vector serving as the negative control. For transient expression assays, Agrobacterium tumefaciens strains harboring either the effector or reporter constructs were co-infiltrated into the abaxial side of Nicotiana benthamiana leaves using Agrobacterium-mediated transformation. After 48 hours of incubation, luciferase activity was measured using a dual-luciferase reporter assay system. The relative LUC/REN ratio was calculated according to the instructions of the dual luciferase reporter gene assay kit.

Data were processed, statistically analyzed and calculated using Microsoft Excel software. Data visualization was performed with OriginPro 2023. Independent two-sample t-tests were performed using IBM SPSS Statistics 27.0, with P<0.05 considered statistically significant. Gene expression profiles were compared using Pearson correlation analysis.

By integrating the hidden Markov model (HMM) and BLAST search, a total of 55 WRKY gene family members were identified from the bayberry genome data, and were systematically designated as MrWRKY01 to MrWRKY55 based on their chromosomal locations. The MrWRKY proteins exhibited significant variations in their physicochemical properties, including amino acid length and theoretical isoelectric point values ( Supplementary Table 1 ). Specifically, the amino acid numbers of MrWRKYs varied from 131aa (MrWRKY02) to 990 aa (MrWRKY45), the molecular weights ranged from 14879.22 Da (MrWRKY02) to 108143.45 Da (MrWRKY02). And the pI values fluctuated between 4.94 (MrWRKY30) and 9.89 (MrWRKY48). In addition, 42 (76.36%) members of the MrWRKYs were localized in the nuclear, 8 (14.54%) members were localized in the cytoplasmic, 3 (5.45%) members were localized in the chloroplast.

To comprehensively understand the evolutionary relationships and functional conservation of the WRKY genes, multiple sequence alignments were performed on WRKY protein sequences from bayberry and Arabidopsis. As shown in Figure 1 , a phylogenetic tree was constructed using the Neighbor-Joining (NJ) method. The WRKY genes of bayberry were systematically classified following the established classification scheme of Arabidopsis WRKY genes. The results showed that bayberry WRKY genes could be divided into five subfamilies: WRKY-Ia, WRKY-Ib, WRKY-IIa, WRKY-IIb, and WRKY-III. Among them, WRKY-Ia subfamily comprised 14 genes, representing the largest subgroup within the MrWRKY gene family, suggesting its highest level of diversity and complexity. In contrast, the WRKY-IIa and WRKY-IIb subfamilies each contained 11 genes, reflecting distinct evolutionary trajectories and adaptive capabilities among WRKY subfamilies. The WRKY-Ib and WRKY-III subfamilies, each consisting of 10 genes, constitute the smallest subgroups, and these quantitative variations likely reflect differential selection pressures and retention mechanisms during evolutionary processes.

The chromosome location analysis of bayberry WRKY genes were conducted by the Tbtools tool. As shown in Figure 2 , 51 MrWRKY genes were mapped across the eight chromosomes of bayberry, while four MrWRKY genes (MrWRKY52, MrWRKY53, MrWRKY54, and MrWRKY55) were localized on the scaffold. These 51 MrWRKY genes were unevenly distributed across eight chromosomes. Notably, Chromosome 2 contained the highest number of MrWRKY genes (14 members), while chromosomes 1, 3, 4, 6, 7, and 8 harbored 2, 7, 5, 9, 8, and 5 genes, respectively. Strikingly, only a single MrWRKY gene was found on chromosome 5. This uneven genomic distribution pattern suggested that the chromosomal allocation of WRKY genes may be associated with structural and functional characteristics of the respective chromosomes.

Interestingly, chromosomes 2 and 5 exhibited the most abundant distribution of WRKY gene subfamilies, encompassing WRKY-Ia, WRKY-Ib, WRKY-IIa, WRKY-IIb, and WRKY-III subfamilies. Considering the pivotal role of chromosome scaffolds in gene expression and regulation, these scaffold-associated WRKY genes may confer distinctive biological functions.

To gain deeper insights into the gene structure and phylogenetic relationships of bayberry WRKY genes, we conducted a detailed analysis of the structural features and conserved motifs of the MrWRKY family members. The evolutionary relationship of the bayberry WRKY genes was presented in Figure 3A . As shown in Figure 3B , six motifs were identified among the 55 bayberry WRKY protein sequences. Specifically, 71.43% of the WRKY-Ia subfamily contained four motifs, while approximately 60.00% of the WRKY-Ib subfamily had three motifs. Similarly, 6 (54.54%) members of the WRKY-IIa subfamily possessed five motifs, whereas 17 (80.95%) members of the WRKY-IIb and WRKY-III subfamilies (excluding MrWRKY15, MrWRKY27, and MrWRKY34) had only two motifs.

Domain analysis clearly revealed that each WRKY gene contained at least one WRKY domain ( Figure 3C ). As a hallmark feature of the WRKY gene family, this domain played a pivotal role in regulating various biological processes and was essential for gene expression and functional activity. By comparing structural differences in domains among different genes, we could further elucidate their evolutionary relationships and functional divergence.

Analysis of the MrWRKY gene family structure revealed that 17 (30.90%) MrWRKY members contained untranslated regions (UTRs) ( Figure 3D ). Specifically, 17 members possessed both 5’ and 3’ UTR exons, with exon numbers ranging from 1 to 13. Notably, MrWRKY43 contained the highest number of exons (13), while MrWRKY26 had only one exon and lacks introns. Within the WRKY-III subfamily, members with three exons accounted for 90.0% of the total. The MrWRKY gene family exhibited significant variation in gene length and intron-exon numbers, while maintaining relatively conserved gene structures within each subfamily. The similarity in gene structures and motif features among subfamilies further validated the reliability of the classification.

To further investigate the potential functions and regulatory mechanisms of bayberry WRKY genes, sequences from the upstream 2000 bp promoter regions of WRKY genes were extracted, and an extensive cis-acting element analysis was conducted using the PlantCARE online tool ( Figure 4 ). Analyzing the promoters of bayberry WRKY genes could help identify cis-regulatory elements and understand their regulatory networks, while also uncovering how WRKY genes interact with other signaling pathways in the plant. The analysis revealed that the promoter regions of two WRKY genes, MrWRKY52 (a member of the WRKY-IIa subfamily) and MrWRKY53 (a member of the WRKY-IIb subfamily), showed no detectable cis-acting elements.

Through a systematic analysis of the bayberry WRKY promoters, this study identified and classified various cis-regulatory elements into three major categories: hormone response elements, light response elements, and abiotic stress response elements. Notably, light response elements were widely distributed across the WRKY gene family. Specifically, 186 light response elements were detected in the WRKY-Ia subfamily, accounting for 44.39% of the total cis-acting elements in this subfamily. Similarly, 139 light response elements were found in the WRKY-Ib subfamily, representing 44.13% of its total cis-acting elements. Additionally, the WRKY-IIa, WRKY-IIb, and WRKY-III subfamilies contained 127, 103, and 116 light response elements, respectively, which accounted for 44.41%, 39.77%, and 39.73% of their respective total cis-acting elements.

In terms of hormone response elements, the WRKY-Ia subfamily exhibited the highest abundance, containing 62 methyl jasmonate (MeJA) response elements (14.79%) and 50 abscisic acid (ABA) response elements (11.93%). Analysis of stress-responsive elements revealed that the WRKY-III subfamily possessed the highest quantity, containing 60 elements in total, accounting for 20.55% of its total cis-acting elements repertoire. The investigation further demonstrated that the majority of WRKY promoters contained multiple cis-acting elements, with MrWRKY54 exhibiting the maximum count of 49 cis-acting elements. The distribution characteristics of these cis- regulatory elements suggested that WRKY genes likely participate in responding to multiple hormone signals and environmental stresses, indicating that their expression undergoes sophisticated multi-tiered regulation during plant growth and development.

To investigate the evolutionary history of bayberry WRKY genes, collinearity analyses were conducted within and across species. As shown in Figure 5 , the results revealed six segmental duplication events involving 12 WRKY genes, accounting for 21.81% of the entire WRKY gene family. Notably, five pairs of duplicated genes were distributed across different chromosomes, indicating that the MrWRKY gene family members experienced chromosomal segment duplication during evolution. Furthermore, a comparative collinearity analysis was performed between bayberry and other species, including Arabidopsis, apple, mandarin orange, peach, grape, and kiwifruit. As shown in Figure 6 , 30 orthologous gene pairs were identified between bayberry and Arabidopsis, while apple, mandarin orange, peach, grape, and kiwifruit shared 81, 54, 52, 50, and 88 orthologous gene pairs, respectively. These findings suggested a closer evolutionary relationship between the WRKY gene families of bayberry and kiwifruit.

To investigate the expression patterns and biological functions of WRKY genes in bayberry fruits, two cultivars of bayberry fruits at different developmental stages were used as experimental materials ( Figure 7 ). The results revealed diverse expression patterns among the bayberry WRKY subfamilies. In the WRKY-Ia subfamily, MrWRKY28 was specifically highly expressed in ‘Shuijing’ at stages S1 and S2, while MrWRKY13 was predominantly expressed during the early fruit development stage in ‘Dongkui’. MrWRKY01 and MrWRKY24 exhibited peak expression levels in ‘Dongkui’ during the fruit ripening stage. MrWRKY14 exhibited higher expression levels in ‘Shuijing’ fruits compared to ‘Dongkui’. Specifically, it showed markedly elevated expression during early developmental stages (S1-S3) in ‘Shuijing’, with expression abundance progressively increasing, followed by a decline during fruit maturation (S4-S5). In contrast, the expression level in ‘Dongkui’ displayed a rapid decrease throughout the entire fruit development process.

In the WRKY-Ib subfamily, four genes showed high expression during the early fruit development stage. Among these, MrWRKY30 and MrWRKY32 were highly expressed in ‘Shuijing’ during fruit ripening, whereas MrWRKY54 displayed the highest expression levels in ‘Dongkui’ at this stage. For the WRKY-IIa subfamily, MrWRKY50 and MrWRKY53 were both highly expressed in ‘Shuijing’ at stages S1 and S2. In the WRKY-IIb subfamily, four genes exhibited relatively high expression levels in ‘Shuijing’ during the early fruit development stage. Notably, MrWRKY47 reached its highest expression at stage S4 in ‘Shuijing’, while the expression levels of MrWRKY23 gradually decreased as the fruit developed.

Within the WRKY II subfamily, MrWRKY18 and MrWRKY27 showed high expression in ‘Dongkui’ at stage S1. In ‘Shuijing’, MrWRKY06 and MrWRKY07 exhibited increasing expression levels as the fruit developed. These findings demonstrate that the bayberry WRKY gene family exhibits a wide variety of expression patterns during fruit development and ripening.

To further investigate the biological functions of the WRKY gene family members during bayberry fruit development, molecular regulatory mechanisms were studied. Based on transcriptome sequencing data, a key WRKY transcription factor member, MrWRKY14, which negatively regulated sugar accumulation in the ‘Dongkui’ fruit, was identified. Bioinformatics analysis revealed multiple potential WRKY transcription factor binding sites in the SWEET1 promoter region, suggesting a possible regulatory relationship between them ( Figure 8A ). Further analysis of the expression patterns of MrWRKY14 and MrSWEET1 during bayberry fruit development showed that both genes exhibited significantly downregulated expression as the fruit matured ( Figures 8B, C ).

Pearson correlation analysis demonstrated a highly significant positive correlation between the expression levels of MrSWEET1 and MrWRKY14 (R = 0.913, P < 0.05) ( Figure 8D ). To validate their regulatory relationship, dual-luciferase reporter assays confirmed that MrWRKY14 significantly enhanced the transcriptional activity of the SWEET1 promoter, with a statistically significant difference compared to the empty vector control (P < 0.05) ( Figure 8E ). Therefore, this study demonstrated that MrWRKY14 negatively regulated fruit sugar accumulation by activating the transcriptional activity of MrSWEET1.