Transcriptome-wide identification of WRKY transcription factors and their expression profiles under different stress in Cynanchum thesioides

Identification of WRKY genes in Cynanchum thesioides

Due to the lack of genomic information, we used different transcriptome data to identify the WRKY gene family to ensure sequence accuracy. The gene sequences of Cynanchum thesioides were obtained from drought stress transcriptome data from previous studies (Zhang et al., 2019b) and unpublished transcriptome data of different anther development from the laboratory. The Hidden Markov Model (HMM) file (PF03106) for the WRKY domain was downloaded from the Pfam database (https://www.ebi.ac.uk/interpro) (Punta et al., 2004). In order to screen out WRKY candidate genes, BLASTp detection was performed on the Cynanchum thesioides transcriptome data using HMM with the e-value set to 0.01. Validation of candidate WRKY genes was performed using an online program (http://pfam.xfam.org/search) containing a conserved domain database search (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and an online tool (http://smart.embl-heidelberg.de/) to delete genes without WRKY domains. Sequences with 97% or more similarity between different transcriptome databases were manually removed. The search results were predicted by searching open reading frames (ORFs) (https://www.ncbi.nlm.nih.gov/gorf/gorf.html) and subsequent experiments were performed on genes with complete ORFs.

Sequence analyses of Cynanchum thesioides

The physical properties of WRKY members, such as length, molecular weight (MW), and isoelectric point (pI), were predicted using the online ExPasy program (http://www.expasy.org/tools/) (Wilkins et al., 1999).

We analysed the conserved motifs using the MEME website (https://meme-suite.org/meme/) and annotated with InterPro Scan (http://www.ebi.ac.uk/interpro/). The conserved motif of CtWRKY was detected by the online tool MEME (https://meme-suite.org/meme/) with parameters set to 0 or 1 occurrence per sequence, up to 10 subjects, and the range of the motif length was set to 5-50 aa (Bailey et al., 2015).

Phylogenetic analysis of Cynanchum thesioides

The amino acid sequence of CtWRKY with a complete ORF was compared with the WRKY gene of Arabidopsis thaliana by multiple sequence alignment using the Clustal Omega tool. Based on the comparison results, we constructed phylogenetic trees using the maximum likelihood method using MEGAX with 1,000 bootstrap replications of maximum likelihood (Tamura et al., 2011). The online ITOL (http://itol.embl.de/help.cgi) tool was used to embellish the phylogenetic tree.

Interaction network analysis of CtWRKY proteins

An interaction network was constructed for all CtWRKY. The homologs of each CtWRKY in Arabidopsis were identified using BLAST software (Chen et al., 2020a). After that, the Arabidopsis homologs corresponding to each CtWRKY gene were submitted to the STRING database (https://string-db.org/cgi/input.pl) (Szklarczyk et al., 2017) with default parameters to construct the interaction network.

Plant materials and treatments

The seeds of Cynanchum thesioides were collected at the Inner Mongolia Agricultural University, Huhot, Inner Mongolia Province, Northwest China (111°69′E, 40° 80′N), and seedlings were raised according to a previous method (Zhang et al., 2019b) and used for each assay in this paper.

Seedlings were treated with abiotic stress and hormones, including 4 °C (cold), 150 mM NaCl (salt), 0.1 mM abscisic acid (ABA), and 500 µM ethephon (ETH) at the ten-leaf stage. Leaves were collected at 0 h, 12 h, and 24 h after treatment and used to analyse the expression patterns of the different treatments. Then, Cynanchum thesioides roots, stems, leaves, flowers, fruits, and anthers at different developmental stages (according to the different lengths of the flower buds, it is named the T1-T4 stage) were collected to analyse the expression patterns of different tissues (Fig. 1). All samples were immediately frozen in liquid nitrogen and stored at −80 °C for subsequent extraction of total RNA. Three replicates were used for all samples.

RNA extraction and quantitative real-time PCR (qRT −PCR)

Total RNA was extracted from Cynanchum thesioides using an RNA kit (TIANGEN BIOTECH, Beijing, China) according to the manufacturer’s instructions. First-strand cDNA was synthesized using TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen Biotech, Beijing, China). Real-time fluorescent quantitative PCR (qRT −PCR) was performed using TB Green Premix Ex Taq II (RR420Q TaKaRa Biotechnology, Beijing, China) on an FTC-3000P (Funglyn Biotech, Toronto, Canada) system for real-time fluorescence quantitative PCR (qRT −PCR), and the primers are shown in Table S1. The reaction procedure was completed under the following program: 60 s of predenaturation at 95 °C, 40 cycles of 15 s at 95 °C and 15 s at 60 °C, and a final hold at 4 °C. All samples were assayed using three replicates. Relative expression levels were calculated using the 2^−ΔΔCT method (Livak & Schmittgen, 2001). By screening four reference genes of Cynanchum auriculatum, we found that ACT7 was the most stable under different abiotic stresses and hormone treatments and that GAPDH was the most stable in different tissues, so we selected ACT7 and GAPDH as internal control (Zhu et al., 2022).

Identification and characteristics of the WRKY gene family

WRKY structural domains were analysed by an HMMER search for candidate genes using SMART and NCBI CDD tools, confirming that these genes contain single or double WRKY structural domains and indicating that they are indeed members of the WRKY gene family. 96 WRKY genes were predicted by searching open reading frames (ORFs) in C. thesioides. Finally, 19 WRKY genes with complete ORFs were obtained. Since there is no standard nomenclature for WRKYs, the 19 WRKYs were designated CtWRKY1 to CtWRKY19 based on their gene IDs. A similar nomenclature is found in Pinus massoniana Lamb (Yao et al., 2020). All of these CtWKRY genes more conserved WRKYGQK structural domain and a C2H2 or C2HC-type zinc finger motif. CtWRKY ORFs range in length from 172aa (CtWRKY1) to 569aa (CtWRKY14), in molecular weight from 18.8 kDa (CtWRKY1) to 63.0 kDa (CtWRKY14), and in predicted isoelectric point from 4.81 (CtWRKY5) to 9.71 (CtWRKY10). Subcellular localization predictions showed that all CtWRKY genes were localized in the nucleus (Table 1).

Phylogenetic analysis and conserved motif detection of CtWRKY proteins

Based on previous studies, 69 AtWRKY proteins from different subgroups were randomly selected as representatives of Arabidopsis for comparison with Cynanchum thesioides (Eulgem et al., 2000; Yao et al., 2020). The most distinctive feature of WRKY proteins is the WRKY structural domain containing 60 amino acids, including the highly conserved feature “WRKYGQK” and a C2H2- or C2HC-type zinc finger motif (Li et al., 2016). The results showed that the structural domain sequences of CtWRKY all contained the conserved heptapeptide “WRKYGQK” (Fig. 2). Phylogenetic trees were constructed based on the WRKY protein sequences of Arabidopsis and Cynanchum thesioides. The results showed that the 19 CtWRKYs were classified into three categories (I-III) according to the taxonomic criteria of Arabidopsis, with seven CtWRKYs in group I, eight CtWRKYs in group II, and four CtWRKYs in group III. Group II had a further division into five categories, with one CtWRKY in IIa, two each in IIb and IId, none in IIc, and three in IIe (Fig. 3).

To study the conserved motifs of WRKY proteins in Cynanchum thesioides, we analysed the conserved motifs using the MEME website and annotated with InterPro Scan. A total of 10 predicted conserved motifs, called motifs 1–10, were identified with amino acid lengths ranging from 9–50, and the number of motifs in each protein sequence ranged from two to ten. The first four motifs were annotated as DNA-binding domain (Fig. 4). Among these motifs, motif 1 is an important structure of the WRKY structural domain present in all CtWRKY proteins. In addition to motif 1, motif 2 is also widely available in most members. The higher number of motifs in group I compared with other subgroups and the presence of motif 4 only in group I suggest that these genes may have multiple functions. Within a subpopulation, CtWRKY proteins have similar motifs and are relatively conserved. In addition, Although CtWRKY9 belongs to group II a, it shares a similar conserved motif with groups II b and II d.

Interaction network analysis of CtWRKY proteins

To better understand the potential interactions of CtWRKY proteins, interaction networks CtWRKY proteins were constructed using STRING software. These WRKY-interacting proteins mainly include VQ proteins (MKS1) involved in regulating plant defense responses, MPK3 and MPK4 proteins involved in plant responses to pathogens and stress, IKU1 proteins involved in regulating endosperm growth and seed size, and other proteins (Petersen et al., 2010; Berriri et al., 2012; Yan et al., 2021; Wang et al., 2010). Figure 5 shows the predicted interaction network between the 11 CtWRKY proteins VQ and MPK proteins. Among them, AtWRKY33 (CtWRKY2) had a strong interacts with 10 VQ proteins and 3 MPK proteins. In addition, AtWRKY53 (CtWRKY16), AtWRKY40 (CtWRKY9), and AtWRKY70 (CtWRKY6) interact with multiple VQ proteins to form a complex regulatory network. Previous studies suggest that there are interactions between WRKY proteins (Chi et al., 2013). Through interaction networks, we found that CtWRKYs are involved in biotic stress response and defence (AtWRKY33/CtWRKY2, AtWRKY40/CtWRKY9, AtWRKY53/CtWRKY16) (Zheng et al., 2006; Lozano-Durán et al., 2013; Wang et al., 2020), abiotic stress response (ZAP1/CtWRKY14) (Qiao, Li & Zhang, 2016), SA and JA signalling (AtWRKY22/CtWRKY11) (Kloth et al., 2016), salicylate signalling (AtWRKY70/CtWRKY6) (Ulker, Shahid Mukhtar & Somssich, 2007), pollen development (AtWRKY2/CtWRKY12) (Lei et al., 2017), and other biological processes. These results suggest that CtWRKYs may play important roles in various biological processes, which may provide important clues for understanding and validating the role of CtWRKYs in the response to various stresses.

Expression pattern of the CtWRKY genes in different tissues

To investigate the possible role of the CtWRKY gene in the growth and development of Cynanchum thesioides, we performed qRT–PCR expression analysis of five tissues of Cynanchum thesioides, root, stem, leaf, flower, and fruit, as well as anthers at different developmental stages (Fig. 6). The results showed that most CtWRKY genes were highly expressed in at least one tissue and only at low levels in other tissues; a result that may be related to the interaction of WRKY transcription factors with other genes or proteins during plant growth and development. There were 11, five, and eight WRKY genes highly expressed in the roots, stems, and flowers of Cynanchum thesioides. Respectively, CtWRKY5 and CtWRKY12 are the most highly expressed in roots and CtWRKY2, CtWRKY7, and CtWRKY18 are the most highly expressed in flowers. Only four genes were highly expressed in fruits, while only two were expressed in leaves. These 19 CtWRKY genes were expressed in different tissues, indicating that CtWRKY genes have a tissue-specific expression profile.

19 CtWRKY genes were analysed by qRT–PCR to determine their expression at different developmental stages of the anther (Fig. 7). The results showed that five genes (CtWRKY2, CtWRKY8, CtWRKY9, CtWRKY12, and CtWRKY17) were downregulated in expression more than five-fold at T2, and six genes (CtWRKY6, CtWRKY7, CtWRKY10, CtWRKY15, CtWRKY18, and CtWRKY19) were expressed at T3. The expression of two genes (CtWRKY1 and CtWRKY14) were highly upregulated by six-fold and five-fold at T4 (Fig. 7). The results suggest that the WRKY gene family plays an important role in the development of Cynanchum thesioides anthers.

Expression pattern of the CtWRKY genes under ETH and ABA treatments

It has been shown that ETH and ABA play important roles in regulating plant development, stress responses, and various physiological processes (Lin, Zhong & Grierson, 2009; Ku et al., 2018). Therefore, to investigate the role of the CtWRKY genes in response to ETH and ABA, we analysed the expression pattern of the CtWRKY genes under ETH and ABA treatment by qRT–PCR. The results showed that three genes (CtWRKY7, CtWRKY18, and CtWRKY19) were upregulated at 12 h and 24 h under ETH treatment, with expression up to 10-fold compared with 0 h. Two genes (CtWRKY2 and CtWRKY9) showed a trend of upregulation followed by downregulation, peaking at 12 h, which was three times higher than that at 0 h. In contrast, four genes (CtWRKY1, CtWRKY8, CtWRKY16, and CtWRKY17) showed a significant decreasing trend after ETH treatment (Fig. 8A).

In addition, we found that CtWRKY18 showed the same expression pattern under ABA treatment as under ETH treatment, and other CtWRKY genes showed significant dynamic changes when subjected to ABA treatment. Among them, four genes (CtWRKY9, CtWRKY11, CtWRKY14, and CtWRKY15) showed significantly upregulated expression, while two genes (CtWRKY6 and CtWRKY8) showed significantly downregulated expression. In addition, CtWRKY1 and CtWRKY19 were significantly upregulated at 12 h, whereas CtWRKY13 and CtWRKY16 expression reached the highest level at 12 h and then decreased, but it remained significantly higher than that at 0 h (Fig. 8B).

Expression pattern of the CtWRKY genes under salt and cold stress

Abiotic stresses severely affect plant growth and geographic distribution. To further investigate the potential role of the CtWRKY genes under abiotic stress, we analysed the expression pattern of the CtWRKY genes under low temperature and salt stress by qRT–PCR. The results showed that all CtWRKY genes responded positively under low-temperature stress treatment. Among them, five genes (CtWRKY2, CtWRKY5, CtWRKY6, CtWRKY7, and CtWRKY10) showed a gradual downregulation under low-temperature stress induction. CtWRKY8 and CtWRKY13 were significantly upregulated at 24 h and were four-fold higher than the control (0 h). In contrast, CtWRKY18 and CtWRKY19 were significantly upregulated at all time points and reached the highest expression at 24 h, which was 30-fold higher than the control (0 h) (Fig. 8C). In addition, CtWRKY18 maintained the same expression pattern under salt stress treatment. On the other hand, CtWRKY3, CtWRKY9, and CtWRKY19 were induced by salt stress and reached the highest expression level at 12 h of treatment, which was 20-fold higher than the control (0 h). Four genes (CtWRKY1, CtWRKY8, CtWRKY16, and CtWRKY17) were significantly downregulated at various time points of salt stress treatment (Fig. 8D).