Genome-wide identification and expression analysis of the WRKY gene family in Rhododendron henanense subsp. lingbaoense

Prediction of physicochemical properties

The physicochemical properties of the RhlWRKY protein were predicted using the compute PI/MW tool of ExPASy database (http://web.expasy.org/compute_pi/), and subcellular localization was predicted using the WoLFPSORT website (https://www.genscript.com/wolf-psort.html).

Gene identification and phylogenetic tree construction

Firstly, InterProScan (v5.50-84.0) was used to annotate the RhlWRKY protein domain based on InterPro protein database, and the annotation results of pfam gene family were extracted. Other transcription factor families not annotated by InterPro are annotated with PlantTFDB (https://planttfdb.gao-lab.org/). A total of 65 putative WRKY genes were identified, and identify the candidate RhlWRKY from the genome of Rhl (SAR: PRJNA656593) (Zhou et al., 2022). Next, MEGA7.0 was used to conduct cluster analysis of the 90 Arabidopsis protein sequences and the 65 putative RhlWRKY proteins to compare their evolutionary with 1,000 replicates bootstrap analysis for statistical reliability. We further performed a NJ phylogenetic tree of the whole RhlWRKY protein sequences relationships.

Gene structure and conserved motif analysis

A phylogenetic tree was constructed using MEGA7.0, and then the RhlWRKY Protein sequence was used to find conserved motifs by using the Multiple Em for MEME online website (https://meme-suite.org/meme/) with zero and one per sequence, maximum number of motifs sets at 12 and optimum width of motif from 6 to 200. The motif with an e-value less than 1e−10 was retained for further analysis. The genetic structure was determined for the RhlWRKY genes based on the genome Gff3 annotation file. The Gene Structure View (Advanced) software of the TBtools (v2.038) was used to create a phylogenetic tree, motif, and genetic structure diagrams.

Cis-acting elements

The 2,000 bp upstream sequences of the start codons of RhlWRKY genes were collected to obtain the promoter sequences using the GXF Sequences Program of the TBtools. Extract of the and analyzed using the PlantCARE database (Lescot et al., 2002) to predict potential cis-acting elements.

Chromosomal location and collinearity analysis

A chromosomal location map was generated using TBtools and collinearity analysis was performed using McScanX software. The BLAS-BLAST GUI Wrapper algorithm was used for Rhl protein internal self-comparison, the Text Marge for MCscanX for regularizing the gff format file, and the Quick Run MCScanX Wrapper algorithm for self-comparison. Repeat events and collinearity relationships of the RhlWRKY genes were analyzed. The Ka/Ks values for duplicate RhlWRKY gene pairs were calculated using KaKs_Calculator2.0.

Gene ontology functional annotation

Gene Ontology (GO) annotations for each RhlWRKY gene were obtained using Blast2GO software by aligning between the WRKY protein sequences and the NCBI non redundant protein (Conesa et al., 2005). Subsequently, WEGO software (https://wego.genomics.cn/) was used for GO functional classification and to determine gene function distribution (Zhang & Wang, 2005).

Stress treatments and transcriptomic analysis

Seeds of Rhl were collected from Xiaoqinling, Lingbao City, Henan Province. The seeds were decontaminated, sterilized and cultured on MS medium at 25 °C (16 h light/8 h dark). Two-month-old seedlings were subjected to treatments with low-temperature (4 °C), 20% PEG, NaCl (200 mM), and MeJA (200 µM) solutions. Samples were collected at 0, 1, 3, 6, and 12 h with light after treatment and stored at −80 °C. Transcriptome analysis was conducted in May 2023, using material from the hypocotyls of 14-day-old seedlings, as well as roots, stems, leaves, and fully bloomed flowers from three-month-old tissue-cultured plants. Samples were stored at −80 °C for transcriptome sequencing, which was performed in triplicate (three biological and three technical replicates; sequencing services provided by Wuhan Frasergen Bioinformatics Co., Ltd., Wuhan, China). Using SOAPnuke software (v2.1.0) to filter Raw reads, high-quality Clean reads are obtained by removing paired reads containing connectors, reads with a N ratio greater than 0.5%, and reads with a mass value Q ≤ 20 and base number accounting for more than 50% of the entire read. Clean reads were compared with the reference genome using HISAT (v2.2.1) software. The prediction of new transcriptomes was performed using StringTie v1.3.4d software, using the Fragments Per Kilobase of exon model per million mapped fragments values (FPKM values) as an indicator of gene expression level. The distribution of reads on the reference genome generally shows only the distribution of the first 25 longest chromosomes, or Scaold (Fig. S1). DESeq v1.22.2 was used to calculate FPKM ratio (FC, Fold change) and False discovery rate (FDR, False discovery rate) between the difference comparison samples. |log2FC| ≥ 1 and FDR < 0.05 were used as DEGs screening criteria. The statistical power of this experimental design, calculated in RNA SeqPower was shown in Table S1. A transcriptional pattern heatmap of the RhlWRKY gene family was generated using the HeatMap feature in TBtools software (with log2(FPKM) Col Scale settings). The RNA-seq data could be searched from the Biosample with accession numbers SAMN37367964–37367966 (root), SAMN37367967–37367969 (stem), SAMN37367970–37367972 (leaf), SAMN37367973–37367975 (hypocotyl), SAMN37367976–37367978 (flower).

qRT-PCR treatments

RNA was extracted using a magnetic bead total RNA kit from Shanghai Lingjun Biotechnology Co., Ltd., Shanghai, China. The concentration of total RNA and OD260/OD280 values were determined using an ultramicro spectrophotometer (NanoDrop One, Thermo Scientific, Waltham, MA, USA). Electrophoresis was used to identify the integrity of RNA. This RNA was used to synthesize cDNA by reverse transcription. Primers for five genes (see Table S2) were designed using GenScript (https://www.genscript.com) software and synthesized by Sangon Biotech (Shanghai) Co., Ltd., Shanghai, China. The EF1α gene was used as a reference gene. The cDNA from Rhl plants treated with NaCl, PEG, MeJA, and low-temperature treatments was used in qPCR reactions. Relative expression levels were calculated using the 2^−ΔΔCt method. Each response had a negative control group and contained three biological and three technical replicates. Actin gene (EF1α) was used as internal reference gene. Results of NTC was not detected. The obtained results were inputted into GraphPadPrism (9.3) software (grouped –mean set-summary data-separated bar graph) to plot expression model diagram, and one-way ANOVA test was used to assess differences by SPSS (13.0) (Analyze–Compare Means-One-way ANOVA), and the P-value of pairwise comparison was obtained by LSD test.

Characteristics analysis

The results, as shown in Table 1, indicate that the 65 RhlWRKY genes (RhlWRKY_1- RhlWRKY_65) exhibit significant variations in the number of encoded amino acids, ranging from 112 amino acids (RhlWRKY_18) to 1,037 amino acids (RhlWRKY_41). The molecular weight of RhlWRKY proteins ranges from 13.0801 kDa (RhlWRKY_41) to 114.6188 kDa (RhlWRKY_18), and their isoelectric points (pI) range from 4.39 (RhlWRKY_34) to 9.87 (RhlWRKY_41), with an average pI of 6.96. Among these, 26 are alkaline (pI > 7) and 39 are acidic (pI < 7). The hydrophilic coefficients are <0 for all RhlWRKY proteins, indicating all RhlWRKY proteins are hydrophilic proteins. Subcellular localization analysis reveals that the majority of these proteins are predicted to be localized in the cell nucleus, with only a few (RhlWRKY_13, _18, _19, _56, _64, and _65) predicted to be in the chloroplast and cytoplasm. The results revealed that these newly-identified RhlWRKY proteins in Rhododendron henanense exhibited a various subcellular distribution, which may be associated with functional diversification in abiotic stress responses. The similar observations for WRKY proteins had been reported in maize (Hu et al., 2021).

Conservation domain and phylogenetic analysis of the WRKY gene family

A multiple sequence alignment was performed for the conserved domains of the RhlWRKY family, and a SeqLogo was generated (see Fig. S2). The majority of RhlWRKY transcription factors exhibit conservation in the N-terminal heptapeptide domain (WRKYGQK) and the C-terminal zinc finger (C2HH/C). These domains are similar to the amino acid variation patterns observed in WRKY proteins from rice and other species (Zhang & Wang, 2005). In addition to the conserved WRKYGQK sequence, three variants were identified: WKKYGEK in Group I (RhlWRKY36), WRKYGKK in Group II-c (RhlWRKY64), and WRKYGRK in Group II-d (RhlWRKY65).

A phylogenetic tree (see Fig. 1) was constructed by aggregating sequences of WRKY genes from Rhl (65) and Arabidopsis (90). WRKY transcription factors were categorized into three subgroups, with Group II further divided into five subclasses: Group II-a, -b, -c, -d, and -e. The three subgroups are in the following order, from largest to smallest: Group II, Group I, and Group III, with a total of 92, 37, and 27 WRKY genes, respectively.

Gene structure and conserved domains of RhlWRKY

Figure 2 shows that RhlWRKY genes contain 0–9 introns and 1–10 exons, and genes within the same subtribe exhibit similar exon-intron structures. This indicates that the 65 RhlWRKY genes have relatively conserved genetic structures. A total of 12 conserved motifs were identified within the RhlWRKY genes, and are designated as motifs 1–12. Motifs 1, 2, 8, and 4 are found in the majority of Rhl genes, indicating a high level of conservation in these four motifs. Among these, motif1 and motif3 contain the WRKYGQK conserved domain. Motif3, with the presence of CC (cysteine), and motif5, with HH/HC, together form the complete C2HH/C2HC zinc finger. Motif1 (containing C), motif2 (containing C), and motif8, or nearby sequences (containing HH), form a complete C2HH zinc finger. As shown in Fig. 2, nearly all genes contain conserved domains and zinc fingers. In Group I, there are two WRKY domains and two C2H2 zinc fingers. In Group II, there is one WRKY domain and one C2H2 structure, with Group II-a and -b containing a unique motif6 and Group II-d and -e mostly containing motif10, suggesting genes in these groups may have similar functions. In Group III, there is one WRKY domain and one C2HC zinc finger. Within the same subtribe, conserved motifs exhibit similarities, indicating that genes within the same subgroup likely share similar biological functions.

Analysis of cis-acting elements

A total of 84 elements were identified in the upstream promoter regions of RhlWRKY (Fig. 3), and the results revealed the presence of various cis-elements in the flanking regions associated with stress, hormones, transcription and development etc. The majority of these elements fall into the stress response category, including elements associated with low-temperature, injury, and drought responses. Among the genes, 63 contain low-temperature response elements (MYB), 33 sequences exhibit wound response elements (WUM-motif), and all 65 genes contain drought response elements (MYC). Additionally, 28 genes have TC-rich repeat sequences, known to protect plants from damage in adverse conditions (Sun et al., 2016), and 18 genes feature light response elements (MRE). A total of 37 genes contain W-box elements, participating not only in wound and pathogen stress responses (Jiang et al., 2016) but also in self-regulation and cross-regulation with other genes (Chi et al., 2013). In terms of hormone response elements, 52 genes contain MeJA-induced response elements (CGTCA-motif), and 55 genes have abscisic acid response elements (ABRE). All 65 genes contain transport cis-acting elements (CAAT-box).

Chromosome localization, collinearity analysis and Ka/Ks analysis

The results of chromosomal location analysis (as shown in Fig. 4) indicate that out of the 65 RhlWRKY gene family members, only two genes, RhlWRKY_46 and RhlWRKY_51, were not mapped to definitive chromosomal locations. It may be due to the fact that the scaffold could not be assigned to a chromosome. The remaining RhlWRKY family genes are distributed across 12 chromosomes, with the exception of chromosome 9, indicating a relatively small preference for the chromosomal distribution of RhlWRKY genes. The highest numbers of genes mapped to Chr3, Chr8, and Chr12, each with eight genes, while the fewest genes mapped to Chr1, Chr10, and Chr11, each with three genes. Two genes separated by five or fewer genes and located in a chromosomal segment less than 100 kb in length are considered tandem duplicate genes (Liu et al., 2020). Based on these criteria, a total of eight pairs of tandem duplicate sequences were identified (see Fig. 4, and Table S3), including four pairs on chromosome 5, two pairs on chromosome 8, and one pair on chromosome 6. Collinearity analysis revealed the presence of 97 segmental duplication relationships among the 65 RhlWRKY genes across the 12 chromosomes (see Fig. 5, and Table S3). To study evolutionary factors, Ka/Ks values were calculated for 105 duplicate gene pairs (see Table S3). A Ka/Ks ratio less than 1 indicates that gene pairs are in a state of negative selection or purifying selection, a Ka/Ks ratio greater than 1 indicates positive selection, and a Ka/Ks ratio equal to 1 suggests neutral selection (Wang et al., 2018). The homologous gene pair of RhlWRKY_18 and RhlWRKY_60 has a Ka/Ks ratio greater than 1, indicating positive selection, but the remaining 104 homologous pairs have Ka/Ks values less than 1, suggesting that the vast majority of RhlWRKY genes underwent purification due to environmental pressures.

GO functional annotation analysis

GO annotation was used to classify gene functions into three categories: molecular functions, cellular components, and biological processes. As shown in Fig. 6, WRKY proteins are almost entirely absent in the cellular component category. Instead, these genes are primarily associated with molecular functions related to binding and transcription regulation activity. In terms of biological processes, WRKY proteins are involved in metabolic processes, biological regulation, regulation of biological processes, and cellular processes.

Transcriptional abundance analysis of the RhlWRKY gene family

Transcriptional patterns were determined based on RNA-seq data from five different tissues of Rhl: roots, stems, leaves, flowers, and hypocotyls (see Fig. 7). RhlWRKY_56 and _65 were not detected in any of the five tissues, suggesting these might be pseudogenes. There were substantial differences in the expression profiles of various RhlWRKY gene family members across different tissues. Among the 43 RhlWRKY genes, 66.15% were expressed in all five tissues, with 24 of them showing relatively high expression levels (36.92%, FPKM > 4). Combining the RhlWRKY cluster analysis, it is evident that, in Group I, apart from RhlWRKY_7, the rest of the genes exhibited relatively high expression levels in all five tissues, with RhlWRKY_14 exhibiting the highest expression level. In Group II, RhlWRKY_5, _6, _18, _21, _37, _38, _42, _44, and _60 showed relatively high expression levels in all five tissues. RhlWRKY_4, _9, _30, _34, _41, _57, and _64 exhibited lower expression levels or were undetectable in all five tissues. In Group III, RhlWRKY_54, _59, and _62 showed relatively high expression levels in all five tissues, whereas RhlWRKY_17, _25, and _48 had lower expression levels. As shown in Fig. 7, RhlWRKY_14 is highly expressed in all tissues (with FPKM values of 236.62, 390.34, 360.22, 255.24, and 322.50 in roots, stems, leaves, flowers, and hypocotyls, respectively). Additionally, tissue-specific expression patterns were observed, with RhlWRKY_37 showing specificity in stems, leaves, and hypocotyls. Genes with relatively high expression levels in stems included RhlWRKY_6 and _36, and RhlWRKY_60, _42, _21, and _6 exhibited higher expression levels in hypocotyls. Furthermore, RhlWRKY_18 displayed specific expression in flowers, and RhlWRKY_42 and _37 showed tissue-specific expression in roots. These genes likely play key roles in their respective tissues.

Analysis of expression patterns under abiotic stress

Expression analysis was performed on five genes belonging to three subtribes of the WRKY gene family. The five selected genes (RhlWRKY_17,_19,_37,_42,_45) show relatively high expression levels in the selected tissues and were orthologous to Arabidopsis (Fig. 8). Under MeJA treatment, WRKY_17 exhibited significant upregulation at 12 h, and WRKY_37 showed significant down-regulation at 3, 6, and 12 h. WRKY_42 showed significant up-regulation at 1 and 3 h, and WRKY_45 had significant up-regulation at 1 and 12 h with respect to 0 h. Under low-temperature stress, WRKY_17 showed significant up-regulation at 1 and 3 h, WRKY_19 was up-regulated at 3 and 12 h, WRKY_37 was up-regulated at 3 h and further at 12 h, WRKY_42 showed increased expression at 12 h, and WRKY_45 showed higher expression at 3 and 12 h with respect to 0 h. In response to drought stress, WRKY_17 exhibited upregulation at 3 h, WRKY_19 showed significant downregulation at 1 h, and WRKY_37 showed downregulation at 3 and 12 h, with a more significant decrease at 6 h. Additionally, WRKY_42 showed significant upregulation at 1 h, and WRKY_45 had significant upregulation at 3 h and 12 h with respect to 0 h. Under high-salinity stress, WRKY_19 had upregulation at 6 h and 12 h, while WRKY_42 showed a more significant down-regulation at 3 h and 6 h with respect to 0 h. WRKY_17 exhibited a 10-fold upregulation at 12 h under MeJA treatment with respect to 0 h, suggesting this is a major gene that responds to MeJA treatment. WRKY_19 showed an 8.7-fold upregulation at 12 h under low-temperature treatment with respect to 0 h, indicating it as the major gene to respond to low-temperature stress. WRKY_42 showed a 4-fold upregulation at 1 h under drought treatment with respect to 0 h, indicating it as a major gene involved in drought response. Interestingly, WRKY_42 exhibited relatively high expression levels across various tissues (roots, stems, leaves, flowers, and hypocotyls). WRKY_19 showed a 4.5-fold upregulation at 6 h under high-salinity treatment with respect to 0 h, indicating that this gene likely acts to counter high-salinity stress.