Initial data preparation involved downloading the genome sequence and annotation file of R. chingii from the Rosaceae genome database (https://www.rosaceae.org/), obtaining the WRKY domain file with the Pfam database number PF00403 (http://pfam-legacy.xfam.org/), and downloading all 71 AtWRKY protein sequences from The Arabidopsis Information Resource (TAIR; https://www.arabidopsis.org/). Subsequently, WRKY domain and all AtWRKY sequences were compared using the “Simple HMM Search” and “Blast Several Sequences to a Big Database” programs of the TBtools software, and the intersection of the comparison results was uploaded to the InterProScan database (https://www.ebi.ac.uk/interpro/result/interprosca/). Finally, all members of the RcWRKY family were identified (Chen et al., 2023).Simple Hmm Search results determine members based on an E-value less than 1e-10. The parameter setting is an E-value less than 1e-5, and the blast program is executed.

To better understand the basic information of the family members, this study analyzed the key physicochemical properties, such as amino acid quantity, molecular weight, isoelectric point, and instability coefficient of RcWRKY, using the Expasy ProtParam online tool (https://web.expasy.org/protparam/). Position information of the family members was extracted using the “GXF Gene Position & Info. Extract” program of the TBtools software (Chen et al., 2023; Su et al., 2023). Finally, specific locations of the family member proteins in cells were predicted using the Bologna Unified Subcellular Component Annotator database (https://busca.biocomp.unibo.it/).

To better understand the grouping of the family members, we constructed a phylogenetic tree of AtWRKY and RcWRKY. We aligned all protein sequences using the MAFFT v7.471 program, constructed a phylogenetic tree using the MEGA7 software, and displayed the phylogenetic tree on the Chiplot website (https://www.chiplot.online/). Based on the classification of the 71 AtWRKY protein sequences in TAIR database (https://www.arabidopsis.org/), RcWRKY proteins were grouped together.

MEME online tool (https://meme-suite.org/) was used to analyze the motifs in the RcWRKY protein sequence, resulting in the identification of 10 conserved motif sequences. These conserved motif sequences were uploaded to the InterProScan database (https://www.ebi.ac.uk/interpro/result/interprosca/) for functional annotation. Finally, Gene Structure View of the TBtools software was used to visualize the motifs, conserved domains, and gene structures (Chen et al., 2023).

Using the Fasta Extract program of the TBtools software, the 2000-bp promoter sequence preceding the ATG start codon of all RcWRKY genes was extracted. This sequence was uploaded to the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/) for cis-regulatory element prediction, and the number of such elements was counted. Finally, the results were visualized as a heat map using the HeatMap program of the TBtools software (Chen et al., 2023).

Next, collinearity of the RcWRKY family members was analyzed using the MCScanX software, and a collinearity map was drawn using the Advanced Circos function of the TBtools software. Replication types of the RcWRKY family were analyzed using DupGenfinder. Replication types were classified into five categories: Whole-genome (WGD), transposed, dispersed, tandem (TD), and proximal duplication (Chen et al., 2023; Qiao et al., 2019; Wang et al., 2012). The non-synonymous substitution rate/synonymous substitution rate (Ka/Ks) of gene pairs with different replication types was calculated using the Simple Ka/Ks Calculator of the TBtools software (Chen et al., 2023).

The genome and annotation files of maize and grape were downloaded from the Chinese National Genome Database (https://ngdc.cncb.ac.cn/), those of rice and Arabidopsis thaliana from the Ensembl Plants database (https://plants.ensembl.org/), and those of Fragaria vesca, Malus domestica, Pyrus communis, Prunus persica, R. idaeus, and R. occidentalis from the Rosaceae Genome Database (https://www.rosaceae.org/). MCScanX software was used to conduct collinearity analysis of the RcWRKY family members of different species. Finally, the results were visualized as a heat map using the HeatMap program of the TBtools software (Chen et al., 2023; Wang et al., 2012).

Transcriptome data of R. chingii at different stages (MG, GY, YO, and RE) were obtained from the National Center for Biotechnology Information (bioproject number: PRJNA671545). The data were generated using the Illumina sequencing platform, and clean reads obtained from the National Center for Biotechnology Information fastq files were aligned with the R. chingii genome using the Hisat2 v2.0.5 software. The reads of each gene were calculated using the FeatureCounts tool (http://subread.sourceforge.net/), and fragments per kilobase pair per million reads values were calculated based on the transcript length.

Members of the flavonoid synthesis gene family were identified using the same method described in section 1.1. The 10 flavonoid synthesis-related protein structure domain files with Pfam database numbers (Supplementary Table 1) and protein sequences of 23 genes in 10 flavonoid synthesis-related gene families of A. thaliana were downloaded from TAIR database (https://www.arabidopsis.org/). Then, 10 flavonoid biosynthetic gene family members in R. chingii were identified, and a phylogenetic tree was constructed. Expression levels of the RcWRKY and flavonoid synthesis gene family members were statistically analyzed, and the average values for each period were calculated. HeatMap program of the TBtools software was used to visualize the heat map (Chen et al., 2023).

Correlations between members of the RcWRKY and flavonoid synthesis gene families during fruit development were analyzed using the Chiplot website (https://www.chiplot.online/), and a correlation heatmap was constructed (Xie et al., 2023). Simultaneously, W-box sites in the promoters of all flavonoid synthesis gene family members were analyzed, as described in section 1.3. Finally, the binding sites of promoters were drawn using the Simple BioSequence Viewer program of the TBtools software (Chen et al., 2023).

In order to detect the expression of genes at different stages of fruit ripening, R. chingii fruits were harvested at 30 days (MG), 40 days (GY), 50 days (YO), and 60 days (RE) after flowering. The fruits were rapidly frozen with liquid nitrogen and stored at -80°C in the refrigerator. Each replicate consisted of a minimum of 10 fruits collected from 5 to 6 R. chingii trees, meticulously mixed to ensure equal distribution. A minimum of six biological replicates (≥60 fruits in total) were established for each stages of R. chingii, amounting to a total of 240 fruits collected for this study. RNA was extracted from young leaves using the Plant RNA Extraction Kit (Takara, Beijing, China), and RNA purity was determined using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA was reverse-transcribed into cDNA using the PrimeScriptTM RT Reagent Kit with gDNA Eraser (Takara).

RT-qPCR analysis of 12 samples (MG, GY, YO, and RE; three replicates for each treatment, totaling 12 samples) was performed using the SYBR PrimeScript RT-PCR Kit (Takara) and ABI 7500 Real-Time PCR system (ABI 7500; Thermo Fisher, Singapore). Relative gene expression was calculated using the 2^-ΔΔCT method (Livak and Schmittgen, 2001). Fluorescent qPCR primers were designed using the Batch q-PCR Primer Design tool in the TBtools software. All primers are listed in Supplementary Table 2; actin was used as the internal reference gene in R. chingii (Chen et al., 2023). Finally, a bar chart was constructed using GraphPad Prism v.10.

Using the seamless cloning technique, open reading frame sequences of the removed stop codons of RcWRKY34 and RcWRKY37 were ligated into the pAN580 vector, and green fluorescent protein was connected at the N-terminus. These recombinant vectors were introduced into tobacco protoplasts via polyethylene glycol-mediated transformation. Green fluorescent protein fluorescence signals were observed under 470 nm excitation light using a laser confocal microscope to analyze the localization of RcWRKY34 and RcWRKY37 in cells (Yoo et al., 2007).

Open reading frame sequences of RcWRKY34 and RcWRKY37 were ligated into the pGreen II 62-SK effect vector, and the 2000-bp promoter sequence before the ATG start codon of LG07.48 was connected to LUC at the upper end of the pGreen II 0800-LUC vector. Both recombinant vectors were transformed into the Agrobacterium GV3101 strain of root-knot fungus, and tobacco leaves were infected with Agrobacterium. Three leaves were injected from each tobacco plant, and each leaf was injected with 9 areas (1 mL of bacterial solution could cover 9 areas). At least 3 tobacco plants were injected in the same batch to ensure the reliability of the results. Infected tobacco was cultured in the dark for 24 h and transferred to a 25°C constant-temperature incubator under 16 h alternating light and dark conditions for two days. LUC fluorescence signals were observed using a chemiluminescence imaging system (Tanon 5200; Tanon, Shanghai, China). LUC and REN fluorescence activities were measured using a dual-luciferase reporter gene detection kit (DL101-01; Vazyme, Nanjing, China), and the relative LUC/REN ratio was calculated.

In this study, we identified 52 RcWRKY family members in R. chingii. These 52 RcWRKY members were unevenly distributed across six chromosomes, with seven on chromosome 2, eight on chromosome 3, six on chromosome 4, seven on chromosome 5, 14 on chromosome 6, and 10 on chromosome 7 (Supplementary Table 3). Analysis of physicochemical properties revealed that the number of amino acids in these 52 proteins was 123–1448, molecular weight was 13520.87–159535 Da, isoelectric point was 4.53–9.96, and instability coefficient was 33.37–72.11, with an average hydrophilicity coefficient ranging from –1.189 to –0.354 (Supplementary Table 3). Subcellular localization prediction results showed that all 52 RcWRKY proteins were localized in the nucleus (Supplementary Table 3). By constructing a phylogenetic tree of WRKY proteins from A. thaliana and R. chingii and classifying them according to AtWRKY proteins in A. thaliana, this study divided RcWRKY proteins into groups I, II-a–e, and III. Group I contained nine RcWRKY and 14 AtWRKY proteins, group II-a contained three RcWRKY and three AtWRKY proteins, group II-b contained seven RcWRKY and eight AtWRKY proteins, group II-c contained eight RcWRKY and 18 AtWRKY proteins, group II-d contained six RcWRKY and seven AtWRKY proteins, group II-e contained seven RcWRKY and eight AtWRKY proteins, and group III contained 12 RcWRKY and 13 AtWRKY proteins (Figure 1).

We assessed the distribution of the top 10 conserved motifs in the RcWRKY family. Notably, motif sequences in the gene family members of different groups exhibited a high degree of consistency, with some motifs being specific to certain groups, such as motif 5, which was present only in group I, and motif 6, which was present only in groups II-a and II-b (Figures 2A, B). Additionally, motifs 1/2/5/6/10 were conserved in the WRKY domain (Supplementary Table 4). Further analysis of the gene structures of the 52 RcWRKY family members revealed that the gene structures of members within the same group were highly consistent. For example, members of group II-c had 2–4 exons, whereas those of group II-d all had three exons (Figure 2C).

Promoters of the RcWRKY family members were rich in various cis-acting elements, including response elements related to growth, development, hormones, and stress (Figure 3). Among these elements, those related to growth and development, especially photoreceptor elements, were the most abundant. The most common hormone response elements were abscisic acid and jasmonic acid methyl ester response elements. The most abundant stress response elements were anaerobic induction and drought response elements. The heatmap clearly showed that, except photoreceptor elements, the most abundant cis-acting elements in the RcWRKY family were jasmonic acid methyl ester response elements.

Collinearity analysis was conducted within the R. chingii genome, and 20 collinear gene pairs were identified (Figure 4A). Calculations of their Ka/Ks ratio revealed that, except for three gene pairs showing high divergence (RcWRKY1/RcWRKY35, RcWRKY9/RcWRKY36, and RcWRKY25/RcWRKY45), Ka/Ks ratios of other gene pairs were all less than 1 (Supplementary Table 5).The three gene pairs might be due to a low sequence similarity (such as a nucleotide similarity of less than 50%), and it is considered to be excluded from the analysis because the Ka/Ks ratio might be unreliable. Further analysis revealed eight gene pairs of dispersed duplication, three gene pairs of transposed duplication, two gene pairs of TD, and 20 gene pairs of WGD in the RcWRKY family (Supplementary Table 6). Compared to its role in the entire genome, WGD played a much greater role in RcWRKY family expansion (Figure 4B).

Co-linearity analysis of the RcWRKY family members with other species revealed 37 collinear gene pairs with rice, 33 pairs with corn, 60 pairs with A. thaliana, 75 pairs with grape, 83 pairs with F. vesca, 166 pairs with M. domestica, 163 pairs with Pyrus communis, 85 pairs with Prunus persica, 90 pairs with R. Idaeus, and 86 pairs with R. occidentalis (Figure 5). Further analysis revealed that only a few RcWRKY members were collinear with monocotyledonous plants, whereas most were collinear with dicotyledonous plants. For example, members of group II-b were collinear only with dicotyledonous plants. Notably, RcWRKY2 and RcWRKY48 did not exhibit any collinear genes in other species (Figure 5).

Transcriptome data analysis revealed that, among the 52 members, five were not expressed, whereas the other 47 members were highly expressed in fruits at different time points. Specifically, RcWRKY33, RcWRKY6, RcWRKY41, RcWRKY2, RcWRKY23, RcWRKY50, RcWRKY26, RcWRKY38, RcWRKY37, RcWRKY34, and RcWRKY11 were highly expressed during the MG period, RcWRKY13, RcWRKY5, RcWRKY9, RcWRKY12, RcWRKY4, RcWRKY47, RcWRKY42, RcWRKY52, RcWRKY21, RcWRKY31, RcWRKY32, RcWRKY19, RcWRKY24, RcWRKY49, RcWRKY43, RcWRKY39, and RcWRKY18 were highly expressed during the GY period, RcWRKY14, RcWRKY16, RcWRKY8, RcWRKY35, and RcWRKY20 were highly expressed during the YO period, and RcWRKY15, RcWRKY22, RcWRKY45, RcWRKY46, RcWRKY10, RcWRKY3, RcWRKY44, RcWRKY36, RcWRKY40, RcWRKY1, RcWRKY28, RcWRKY51, RcWRKY25, and RcWRKY48 were highly expressed during the RE period (Figure 6).

In this study, we identified 30 flavonoid synthesis-related genes belonging to 10 gene families in R. chingii. These included 12 members of the 4CL family, one member of the C4H family, three members of the F3’H family, two members of the PAL family, two members of the DFR family, three members of the CHI family, two members of the CHS family, two members of the F3H family, one member of the ANS family, and two members of the FLS family (Figure 7A). Further analysis of these 30 flavonoid synthesis-related genes in R. chingii revealed that most members exhibited highest expression levels during the MG period (Figure 7B).

Correlation analysis was conducted between 19 flavonoid synthesis genes highly expressed during the MG period and 11 corresponding WRKY genes. Notably, RcWRKY33, RcWRKY41, RcWRKY38, RcWRKY37, RcWRKY34, and RcWRKY11 were significantly positively correlated with most flavonoid synthesis genes (Figure 8A). Therefore, these six WAKY genes were subsequently selected for quantitative verification. Moreover, among the 19 flavonoid synthesis genes highly expressed during the MG period, 11 contained W-box cis-acting elements in their promoter regions (Figure 8B). From these 11 genes, the genes that contain more than two W-box elements and are significantly positively correlated with the 11 highly expressed waky genes during the same period were selected for subsequent quantitative verification and functional verification.

To verify the reliability of the transcriptome data, we conducted RT-qPCR validation of RcWRKY33, RcWRKY41, RcWRKY38, RcWRKY37, RcWRKY34, and RcWRKY11. RcWRKY37 and RcWRKY34 were highly expressed during the MG period (Figure 9). Simultaneously, we selected LG07.48 from the 4CL gene family, LG02.725 from the F3’H gene family, and LG02.181 from the DFR gene family for validation. All three genes were significantly highly expressed during the MG period (Figure 9).

Subcellular localization studies confirmed that both RcWRKY34 and RcWRKY37 were located in the nucleus, consistent with the predicted results (Figure 10A). Based on Figure 10B, we constructed expression vectors and selected RcWRKY34 and RcWRKY37 as effector genes and proLG07.48 as the reporter promoter to verify flavonoid synthesis regulation by the RcWRKY family in R. chingii. Fluorescence intensity was significantly higher in the tobacco leaves co-transfected with 35S::RcWRKY34 and 35S::RcWRKY37 than those transfected only with proLG07.48::LUC (Figure 10C). Quantitative analysis also confirmed higher luciferase activity in the leaves co-transfected with 35S::RcWRKY34 and 35S::RcWRKY37, suggesting that both RcWRKY34 and RcWRKY37 specifically bind to proLG07.48 and exert positive regulatory effects (Figure 10D).