Comprehensive identification of GSK3 gene family in celery and functional characterization of AgSK22 involvement in salt stress

Glycogen synthase kinases 3 (GSK3s), which are also known as the Shaggy-like protein kinases, are a kind of conserved serine/threonine protein kinases in eukaryotes and serve as important signaling mediators in a wide range of biological processes. In mammals, two GSK3 isoforms (GSK3α/β) are identified with essential roles in a variety of development and disease processes, including glycogen metabolism, determination of cell fate, embryogenesis, tumorigenesis, and regulation of transcription factor stability and activity^[1−4].

GSK3 appears to be more diverse and multiple homologs of GSK3 are identified in plants^[5]. There are ten and nine GSK3s in Arabidopsis and rice, respectively, implying that GSK3s may perform more diverse biological functions accompanied by gene duplication events in the plant kingdom^[6]. Plant GSK3s are also a class of multitasking kinases involved in a large number of developmental processes, including stem cambial activity, xylem and phloem differentiation, vascular tissue development, pollen development, stomatal development, root hair patterning, root elongation, chloroplast development, seed carbon partitioning, seed germination, and plant responses to various abiotic and biotic stresses^[7]. Further studies indicate that plant GSK3s serve as essential regulatory nodes to control plant development and responses through integrating different signal transduction pathways. They perform diverse roles through phosphorylating a large number of downstream substrates, including transcription factors and their cofactors, which could be affected on protein stability, localization, or activity upon GSK3s phosphorylation^[8−10].

In Arabidopsis, GSK3s are uniformly named AtSKs (Arabidopsis thaliana SHAGGY-like KINASES) and divided into four subgroups (I−IV). Subgroups I to IV include 3, 3, 2 and 2 members, respectively^[8]. AtSK21, the first-characterized and well-characterized GSK3 kinase in Arabidopsis from genetic screening, is also named BR-INSENSITIVE 2 (BIN2) due to the insensitivity of its gain-of-function mutant to brassinosteroids (BRs). AtBIN2, together with two other GSK3s, functions as negative regulators of BR signaling through phosphorylating and inactivating two key transcription factors, BRI1-EMS SUPPRESSOR 1 (BES1) and BRASSINAZOLE RESISTANT 1 (BZR1)^[11]. Moreover, subgroup II GSK3s mediate BR function in stomatal development by inhibiting a MAPK cascade^[12]. Subgroups I and II GSK3s regulate xylem cell differentiation in the downstream of TDIF-TDR (tracheary element differentiation inhibitory factor-TDIF receptor) signaling through phosphorylating the transcription factor BES1^[13]. AtBIN2 and AtBIL1/2 were recently demonstrated to play important roles in root development and seed germination by integrating antagonistic ABA and BR (abscisic acid) signaling^[14−19]. Cellulose synthase A1 (CESA1) phosphorylation by AtBIN2 results in a reduction of cellulose biosynthesis^[20]. In addition, suppression of AtSK32 leads to pollen defect indicating its important roles in pollen development^[21].

GSK3s also play important roles in the plant responses to abiotic and biotic stress. AtSK11 directly phosphorylates and activates glucose-6-phosphate dehydrogenase (G6PD), causing enhanced salt stress tolerance in the AtSK11-overexpressing plants^[10]. Drought stress risks many plants worldwide and greatly reduces crop yield. BIN2 positively regulates drought stress through phosphorylating a transcription factor RD26 (Responsive to Desiccation 26) in Arabidopsis, which would activate dehydration-responsive genes^[22]. BIN2 is also known to mediate plant cold tolerance by phosphorylating the transcriptional factor ICE1 (Inducer of CBF Expression1) and reducing its binding to the promoter regions of CBFs (C-repeat binding factors), which is known to be key regulators in plant cold stress response^[23,24]. Additionally, AtSK11 positively regulates the plant innate immunity in Arabidopsis, due to loss of AtSK11 causing increased susceptibility to Pseudomonas syringae infection^[25,26].

GSK3s have also been demonstrated to affect development and stress responses in many other plant species, including rice (Oryza sativa)^[27,28], soybean (Glycine max)^[29], tomato (Solanum lycopersicum)^[30], potato (Solanum tuberosum)^[31], tobacco (Nicotiana benthamiana)^[32], pepper (Capsicum annuum)^[33], cotton (Gossypium hirsutum)^[34], wheat (Triticum aestivum)^[35,36], and barley (Hordeum vulgare)^[37]. Celery, as an important Apiaceae vegetable crop, is widely cultivated worldwide and has significant economic value. It also possesses a variety of pharmacological effects, such as the prevention of constipation due to its richness in dietary fibrin. However, the production of celery worldwide is greatly affected by multiple external cues. It is noteworthy that although GSK3s have been demonstrated to perform versatile functions in plant growth and stress responses in many crop species, the functions of GSK3s in celery are barely known. The genomes information of three typical representative vegetable plants of Apiaceae, including celery^[38], coriander^[39] and carrot^[40], became available recently, which provided a solid foundation for a comprehensive study of key gene families in the whole-genome level for plant resistance breeding. In the current study, a comparative genomics study of GSK3s was conducted among three Apiaceae species, including celery, coriander, and carrot. Lettuce was also analyzed in this study due to its closer relationship with Apiaceae^[41]. Grape served as a reference due to possessing the most similar genome structure to the dicots’ ancestor^[42]. These data provided a basis for further characterization of diverse biological roles of GSK3s in celery, and for improvement of celery varieties through transgenic development or genetic breeding.

To identify the GSK3 family members in celery, the amino acid sequences of ten Arabidopsis GSK3s were used as queries to search the celery genome. In total, 11 putative GSK3s were identified in celery (Apium graveolens). The genomes of grape (Vitis vinifera), lettuce (Lactuca sativa), coriander (Coriandrum sativum), and carrot (Daucus carota) were also analyzed for the GSK3 family members as references, and a total of seven, nine, ten, and ten GSK3 genes were identified, respectively. Consistent with the nomenclature of this gene family in Arabidopsis, these genes were unified as SKs (Shaggy/GSK3-like kinases), and numbered sequentially according to the genomic locus within each subgroup in the phylogenetic tree (Supplemental Table S1).

To explore the phylogenetic relationship of the GSK3s, we constructed a phylogenetic tree based on the GSK3 amino acid sequences in the selected six species. All the GSK3s were categorized into four subgroups (I−IV) in the phylogenetic tree with Arabidopsis thaliana GSK3s as references. Further analysis indicated that the distribution of the GSK3s in different subgroups was not uniform. In celery, Subgroup I (AgSK11, AgSK12, AgSK13, and AgSK14) and Subgroup II (AgSK21, AgSK22, AgSK23, and AgSK24) contains four AgSKs, respectively. Subgroup III contains two AgSK proteins (AgSK31 and AgSK32). However, Subgroup IV has the smallest number of AgSK proteins with only one AgSK41. Moreover, GSK3s within each subgroup in the phylogenetic tree indicated a close genetic relationship among different species. The GSK3s of celery exhibited the closest genetic relationship with those of coriander and carrot, followed by those of lettuce, Arabidopsis, and grape, which is consistent with the closer phylogenetic relationship of the three Apiaceae species (Fig. 1).

Figure 1. 

Phylogenetic analysis of the GSK3 proteins in Arabidopsis (Arabidopsis thaliana, At), celery (Apium graveolens, Ag), coriander (Coriandrum sativum, Cs), carrot (Daucus carota, Dc), lettuce (Lactuca sativa, Ls) and grape (Vitis vinifera, Vv). Clades with colored branches referred to four subgroups I−IV.

Conserved motif analysis of GSK3 family members in celery

To further clarify the evolutionary history of GSK3 family genes, we examined conserved motifs and gene structure based on all the GSK3s identified in the selected six species. The number of both exons and introns of AgSK genes varied from eight to 14. Additionally, four AgSK genes didn’t possess UTR regions, which may be due to inaccurate predictions of some genes in the celery genome. Similar structures of exon-intron in AgSK genes were identified despite there are being some differences in the length of introns. We then performed conserved motifs analysis of AgSK proteins, and identified 15 conserved components (Supplemental Table S2). STKs_GSK3 domain corresponded to the motifs 1−7 in most of AgSKs, whereas this kinase domain contained much fewer motifs in AgSK11, AgSK13, AgSK21 and AgSK23. Whether the variation on the kinase domain was also caused by inaccurate predictions of some genes in the celery genome required to be further studied. Meanwhile, several motifs were identified to be subgroup-specific. For example, the motif 11 was specific to the subgroup III, while the motifs 10 and 13 were lacking in the subgroups II and III, respectively. The results indicated that AgSKs in different subgroups may have undergone functional divergence during the evolution process (Fig. 2).

Figure 2. 

Conserved motifs and gene structure analyses of the GSK3s in celery and Arabidopsis. (a) Phylogenetic tree of the GSK3 proteins. (b) Conserved motifs and domains of the GSK3 proteins. Different color boxes indicated different kinds of motifs. The position of the conserved domain was framed by rectangular boxes. (c) Structure prediction of the GSK3 genes. The coding sequence (CDS) and untranslated region (UTR) were displayed in yellow and blue boxes, respectively. The lines between boxes represent introns.

We extracted the STKs_GSK3 domain sequence from Arabidopsis and celery, and performed multiple sequence alignments. In Arabidopsis, the motifs CDFGSAK, SYICSR, TREE and M/LEYV are reported to be critical for modulating the GSK3 kinase activity. The three kinds of motifs were also present and highly conserved in the celery homologous of AtSKs. Additionally, the motif M/LEYV is the key residue for docking of bikinin, which is the chemical inhibitor of several AtSK proteins. The motif M/LEYV was found to be present in AgSK22, AgSK24, AgSK31, and AgSK32, which is consistent with the results that the activity of some AtSKs could not be inhibited by bikinin in Arabidopsis. The results showed the celery STKs_GSK3 domain contained most of the highly conserved sequences compared with this domain of Arabidopsis GSK3s, indicating highly conserved functions of celery GSK3s (Supplemental Fig. S1).

Collinearity analysis of GSK3 gene family

We investigated chromosome localizations of GSK3 genes in three Apiaceae species. All of the 11 celery GSK3 genes were unevenly distributed on six chromosomes, including chromosomes Ag2, Ag6, Ag7, Ag9, Ag10, and Ag11. Chromosome Ag6 harbored the most three AgSK genes, while Ag2, Ag7, and Ag9 harbored only one AgSK gene on each chromosome. Chromosome Ag10 and Ag11 harbored two AgSK genes, respectively. Moreover, two AgSK genes (AgSK21 and AgSK22) clustered at the beginning of chromosome Ag6 may be due to gene duplication.

In order to explore the evolutionary relationship of the GSK3 family genes, paralogous and orthologous GSK3 gene pairs were investigated among celery, coriander, and carrot. A total of four, six, and six paralogous GSK3s were detected within celery, coriander, and carrot, respectively. Notably, more orthologous GSK3s were identified between celery and coriander than that between celery and carrot, indicating a closer genetic relationship between celery and coriander (Fig. 3, Supplemental Tables S3 & S4).

Figure 3. 

Analysis of the GSK3 paralogs and orthologs. (a) Circos plot of the GSK3 paralogs within each Apiaceae species. (b) Circos plot of the GSK3 orthologs among the three Apiaceae species. The number on the scale represents the physical location in each chromosome. Duplicated GSK3 gene pairs are linked by lines with the corresponding colors. Chromosomal distributions are indicated in the outer layer.

To further clarify the possible duplication mechanisms leading to the expansion of GSK3 genes, gene duplication types were examined. The results showed that 45.45%, 40.0%, and 60.0% duplicated gene pairs resulted from whole-genome duplication (WGD) or segmental duplication in celery, coriander, and carrot, respectively. Additionally, there are 54.55% GSK3 duplicated gene pairs resulted from dispersed duplication in celery. Moreover, the percentages of collinear GSK3s in those three Apiaceae species were considerably higher than the average level within each genome, indicating that WGD brought important effects in the expansion of this gene family (Supplemental Tables S5−S7).

Cis-acting elements analysis of AgSK genes

To evaluate the transcriptional regulation of AgSK genes, putative cis-acting elements were screened in the promoter regions of the 11 AgSK genes. A total of 76 cis-acting elements were identified, and classified into seven categories, including light-responsive elements, promoter-related elements, environmental stress-related elements, hormone-responsive elements, development-related elements, site-binding related elements and other elements. In each category, we calculated the number of detected genes including and not-including replicates, and marked them as the numbers of 'hits' and 'genes', respectively. The most abundant elements were light-responsive elements with 71 hits and 11 genes, followed by the promoter-related elements (38 hits, 11 genes), environmental stress-related elements (37 hits, 11 genes), hormone-responsive elements (36 hits, 10 genes), and development-related elements (31 hits, 11 genes) (Supplemental Table S8). The elements associated with light, development, hormone, and stress were widely distributed among the AgSK genes. Notably, most of the AgSK genes were multiply responsive to several factors.

Since GSK3 genes were demonstrated to perform essential functions in various environmental stress responses in multiple plant species, the involvement of AgSK genes in the environmental stresses was investigated in detail to evaluate their possible functions on abiotic stress. The results showed that one AgSK gene could be induced by several kinds of environmental stresses, and each environmental stress also could induce multiple AgSK genes. For example, six AgSK genes were low-temperature responsive (AgSK11, AgSK13, AgSK21, AgSK23, AgSK24, and AgSK32), and three AgSK genes were drought responsive (AgSK11, AgSK24, and AgSK32). Hence, AgSK11, AgSK24, and AgSK32 were low-temperature and drought dual responsive genes (Supplemental Table S9).

In summary, the results from cis-acting elements analysis implied that AgSKs are also a class of multitasking kinases involved in a wide range of plant developmental processes, and plant responses to various abiotic and biotic stresses in celery.

Expression patterns of AgSK genes

Profiling expression patterns would provide a new perspective for understanding the biological functions of different gene family members. The expression patterns of AgSKs in different tissues were determined by transcriptome sequencing and quantitative real-time PCR (qRT-PCR) in celery (Fig. 4). Most of the AgSKs were expressed in a tissue-specific manner, and AgSKs subjected to different subgroups showed distinct expression patterns. For example, AgSKs in the subgroups I and II were highly expressed in the roots (R) and leaves (L), respectively, from both transcriptomic and qRT-PCR results. However, the expression patterns of AgSK24, AgSK31, AgSK32, and AgSK41 were not consistent in the transcriptomic and qRT-PCR results. AgSK41 was highly expressed in the roots and leaves from the transcriptomic and qRT-PCR results, respectively. AgSK24, AgSK31 and AgSK32 showed quite a low level of expression in all the detected tissues in the qRT-PCR experiment. The different expression patterns of AgSK homologous genes may indicate the functional diversity and differentiation among the family, and they may play different roles in different tissues in celery.

Figure 4. 

Tissues-specific expression patterns of the AgSK genes (root, petiole and leaf). (a) Hierarchical clustering of gene expression profiles of all the AgSK genes. The log₂-transformed FPKM values were utilized to generate the heatmap. (b) Validation by qRT-PCR of tissues-specific expression patterns of the AgSK genes.

Due to GSK3 genes performing important regulatory functions in various environmental stress responses in many plant species, the regulation of the AgSK genes expression was investigated to evaluate their possible functions on abiotic stress response under different environmental stresses. Overall, the AgSK genes were extensively affected by different environmental stresses, indicating the potential critical functions of them in the growth and development of celery. In the subgroup I, AgSK11 and AgSK13 were considerably sensitive to all the three kinds of stresses, while no obvious change in the expression of AgSK12 was observed. All the AgSK genes in the subgroup II were up-regulated under cold, salt, and drought stresses, with an expression peak at 12 h cold treatment, 6 h salt treatment, and 1 h drought treatment, respectively. The AgSKs in the subgroups III and IV (AgSK31, AgSK32 and AgSK41) were up-regulated under salt and drought stresses, not under cold stress. Moreover, the expression peak induced by salt stress occurred at 6 h for all the responsive AgSKs. The results demonstrated that the AgSK genes were extensively involved in the plant responses to these environmental stresses, implying their potentially critical functions in the growth and development of celery. Meanwhile, the different responsiveness of AgSK genes in different stresses indicated the diversity and differentiation of AgSK genes during evolution (Fig. 5).

Figure 5. 

The responses of the AgSKs to (a) cold, (b) salt, and (c) drought stresses were analyzed by qRT-PCR. Different lowercase letters among different timepoints for each gene in each stress indicate significant difference on the relative expression following Duncan’s test (p < 0.05), while the same letters represent no significant difference.

Functional analysis of AgSKs in Arabidopsis

GSK3s control a large number of processes involving plant growth and stress responses in Arabidopsis, and AtSK21 (BIN2) negatively regulates BR signaling through phosphorylating and inactivating the transcriptional factors AtBZR1/BES1. Considering the closest phylogenetic relationship between AgSK21/22/23 and AtSK21 (BIN2) in the subgroup II, we performed yeast two hybrid assays to examine whether the three AgSKs proteins can directly interact with AgBZR1 (a celery ortholog of AtBZR1) in celery. All the AgSK21/22/23 proteins strongly interacted with AgBZR1 in yeast, indicating conserved functions of AgSK21/22/23 in BR signaling in celery (Fig. 6a). We further examine the subcellular localization of AgSK21/22/23 in tobacco leaves. The results showed that all the three AgSK proteins were mainly localized in the nucleus and plasma membrane (Fig. 6b).

Figure 6. 

AgSKs directly interact with AgBZR1. (a) Yeast-two-hybrid assay showed AgSK interacted with AgBZR1 in yeast. The pGBKT7-AgSK and pGADT7-AgBZR1 fusion vectors were co-transformed into AH109 yeast cells. The transformants were screened on SD/-Leu-Trp-His (-LTH) medium. (b) Subcellular localization analysis of AgSK proteins. Tobacco (Nicotiana benthamiana) leaves transiently expressed AgSK-GFP fusion proteins were observed through the laser scanning confocal microscope. Scale bars represent 100 μm.

Previous studies have indicated that multiple GSKs were involved in coordinating salt stress and plant growth in Arabidopsis, especially, AtSK21(BIN2) served as a molecular switch between the salt stress response and plant growth recovery^[43]. Considering salt stress considerably regulated the expression of the subgroup II AgSKs (Fig. 5b), and AgSK22 shared the highest amino acid sequence identity with AtSK21(BIN2), the roles of AgSK22 in the salt stress response were explored using transgenic Arabidopsis plants over-expressing AgSK22-GFP (AgSK22-OE). A transgenic line, which didn’t show obvious phenotypes under normal growth condition, was selected to assess its salt stress response. The seed germination was monitored on 1/2 MS media supplemented with different concentrations of NaCl. No differences were observed in the germination of the transgenic plants without salt stress. However, the germination of the transgenic plants AgSK22-OE was delayed compared to Col-0 at 2 and 3 days on 1/2 MS media containing 150 mM NaCl, suggesting that AgSK22 negatively regulated salt stress response (Fig. 7). Additionally, the roles of AgSK23 in the salt stress response were also explored using transgenic Arabidopsis plants over-expressing AgSK23-GFP (AgSK23-OE), and the seed germination of the transgenic plants AgSK23-OE was also delayed by 150 mM NaCl (Supplemental Fig. S2).

Figure 7. 

Responses of AgSK22 overexpression transgenic Arabidopsis plants to NaCl. (a) AgSK22 protein was detected using an anti-GFP antibody. (b) AgSK22 overexpression transgenic lines showed more sensitive phenotypes when germinated on 1/2 MS media with 0 or 150 mM NaCl for 3 days. 1, 3, 5: Col; 2, 4, 6: AgSK22-OE. (c) & (d) Germination rate of Col-0 and AgSK22-OE overexpression transgenic lines at 2-, 3-, and 4-d on 1/2 MS media containing 0 or 150 mM NaCl. * denotes p < 0.05, as determined by a Student’s t-test.

GSK3s have been demonstrated to regulate a variety of biological processes associated with plant growth, development, and stress responses in many plant species^[44,45]. However, the roles of GSK3s in celery are barely studied. In this study, a genome-wide analysis of the GSK3 gene family was carried out in celery, including conserved motifs, phylogenetic classification, collinearity, cis-acting elements, subcellular localization, expression patterns, and the biological effects of AgSK22 to salt stress. The GSK3 gene family in coriander, carrot, lettuce, grape, and Arabidopsis species were also analyzed as references.

In Arabidopsis, there are ten GSK3 genes named AtSKs. We identified 11, ten, ten, and seven GSK3 genes in celery, coriander, carrot, and grape, and also rename them as corresponding SKs. The quantity of SKs within each Apiaceae species was similar to that in Arabidopsis, which is consistent with the previous report that both the three Apiaceae species and Arabidopsis underwent two rounds of WGD events after divergence from their common ancestor of dicotyledonous species^[38]. Grape didn’t undergo any polyploidization events and preserved the genome structure of the dicots’ ancestor. Hence, if no gene loss or translocation events occurred after polyploidization events, one grape gene was supposed to correspond to four celery orthologs. The number of GSK3 genes was supposed to be 28 since there are seven GSK3 genes in the grape genome. The actual number of GSK3 genes (11) is much fewer than expected (28), indicating extensive gene loss occurred in this gene family during evolution over time. WGD is demonstrated to be followed by loss of many duplicated genes over a long-time span, while the genes performing essential functions tend to be retained, such as transcription factors, kinases, phosphatases, and transporters^[46,47]. A similar quantity of GSK3s was retained among the selected species in the current study indicating its critical and conserved functions in the plant kingdom.

Furthermore, GSK3s in celery (AgSKs) clustered in the same clade in the phylogenetic tree had similar exon/intron organization, suggesting the organizational structure is conserved within each subgroup. In Arabidopsis, the sequence motifs CDFGSAK, SYICSR, and TREE are reported to be critical for modulating the GSK3 kinase activity^[48]. Additionally, the M/LEYV residues are critical for the docking of bikinin, which is the chemical inhibitor of several AtSK proteins^[49]. In this study, a set of motifs corresponding to the GSK3-kinase domain (motifs 1−7) were present in most of the AgSK proteins. The CDFGSAK, SYICSR, TREE, and M/LEYV motifs were found to be highly conserved and also present in the celery AgSKs. Notably, multiple subgroup-specific motifs were identified in the N-terminal extension of AgSK proteins. The results suggested celery GSK3s shared a conserved kinase domain but had variations in the N-terminal extension. These findings also indicated that GSK3 proteins have essential kinase activity among different species in plants.

Different expression patterns of homologous genes in different tissues indicate functional diversity and differentiation within the same gene family. In this study, the AgSKs were expressed in a tissue-specific manner and have different expression patterns in different subgroups, with the subgroup I members highly expressed in root, and subgroup II members highly expressed in leaf, suggesting functional differentiation among homologous genes in the same family. Additionally, AgSK24, AgSK31, and AgSK32 were found to be expressed at a considerably low level. It has been previously reported that duplicated genes may be retained with reduced expression level and promoter complexity after WGD events, and lowly expressed genes may undergo more rapid divergence^[50].

GSK3s were demonstrated to play essential roles in multiple abiotic stress responses in Arabidopsis and rice^[45]. Hence, the regulation of AgSK genes expression was determined under salt, cold, and drought conditions to predict the potential functions of AgSK genes under abiotic stresses. Most of the AgSKs undergo considerable changes in gene expression after being stimulated by those abiotic stress treatments. AgSK11 and AgSK13 in the subgroup I and all the AgSK genes in the subgroup II were sensitive to all the three kinds of stresses. The AgSKs in the subgroups III and IV (AgSK31, AgSK32, and AgSK41) were up-regulated under salt and drought, not under cold stresses. The results suggested that the expression of AgSKs could be regulated by more than one of the environmental stresses and each environmental stress could modulate the expression of more than one AgSKs. For example, AgSK11, AgSK24, and AgSK32 were responsive to both low-temperature and drought stresses, which is consistent with the results from both the cis-acting elements and qRT-PCR analysis. The extensive involvement of the AgSK genes in response to these environmental stresses implied their potential critical functions in the growth and development of celery. Meanwhile, the different responsiveness of the AgSK genes in different stresses indicated their functional diversity and differentiation during evolution, and they were involved in the responses to different stresses.

In Arabidopsis thaliana, multiple AtSKs were reported to be involved in the response to salt stress. AtSK13, AtSK22, AtSK31, and AtSK42 genes are up-regulated by salt treatment. Overexpression of AtSK11 or AtSK22 greatly induced salt tolerance in transgenic plants. By contrast, the gain-of-function mutant of AtSK21/AtBIN2 (bin2-1), exhibited decreased salt tolerance. The role of AgSK22 (a celery ortholog of AtSK21/AtBIN2) involved in salt stress was investigated in this study. One hundred and fifty mM NaCl treatment delayed the germination time of AgSK22 overexpression transgenic Arabidopsis plants, suggesting that AgSK22 negatively regulated salt stress response. The function of AgSK22 is consistent with AtSK21 but in contrast to AtSK11 in the response to salt stress, indicating that different GSK3 members may perform distinct functions in different plant species. The molecular mechanism by which AgSKs regulates stress requires to be further studied.

The genome-wide analysis of AgSKs indicated that the critical motifs modulating the GSK3 kinase activity were highly conserved across the whole family in celery. Moreover, AgSKs were extensively involved in the plant response to cold, salt, and drought environmental stresses. The transgenic Arabidopsis plants that overexpressing AgSK22 were found to be more sensitive to salt stress on seeds germination, suggesting that AgSK22 negatively regulated salt stress response.

Genome information of celery (Apium graveolens), coriander (Coriandrum sativum), carrot (Daucus carota), lettuce (Lactuca sativa), and grape ( Vitis vinifera) were retrieved from CGD (celery, http://celerydb.bio2db.com/), CGDB (coriander, http://cgdb.bio2db.com/), and phytozome (carrot, lettuce, and grape, https://phytozome-next.jgi.doe.gov/), respectively. To identify GSKs family members in the selected species, the amino acid sequences of GSK3 proteins from Arabidopsis were used as queries to search the corresponding genome databases using Blast+ pairwise sequence alignment software^[51]. The parameters used in the current study were as follows: E-value ≤ 1e-5. The putative GSK3 proteins in the selected plant species were further validated via SMART (http://smart.embl-heldelberg.de/), Pfam (http://pfam.xfam.org/) and CDD (https://www.ncbi.nlm.nlh.gov/cdd) databases.

Phylogenetic analysis GSK3s

Multiple sequence alignment of GSK3 proteins from the selected species was performed using ClustalW^[52]. The phylogenetic tree was constructed by MEGA11 software^[53] under the JTT+G model based on 1,000 bootstrap replications.

Conserved motif analysis

Gene structure and conserved motif analysis of GSK3s in celery and Arabidopsis were performed. The structure of the GSK3 genes was predicted using Gene Structure Displayer Server 2.0 (GSDS). The General Feature Format Version 3 (GFF3) file of the GSK3 genes was uploaded to the GSDS website in order to distinguish the positions of untranslated regions (UTRs), exons, and introns in each GSK3 gene. The online MEME tool (https://meme-suite.org/meme/) was employed to perform the conserved motif analysis of GSK3 proteins. The parameters used in the current study were as below: the maximum number of motifs was 15, and the width of motifs varied from six to 200 residues. CFVisual software was used to visualize the results.

Analysis of collinearity

Whole genome protein sequences of the selected species were aligned using the Blastp program with E-value ≤ 1e-5 as the parameters. The MCScanX software^[54] was used to search the collinear blocks and genes based on the alignment results with Num of Blast Hits ≤ 20, and the special GFF files were processed. The paralogous and orthologous of GSK3 genes among three Apiaceae species were identified from the results of collinearity with E-value ≤ 2E-120 and visualized by Circos software. The gene duplication type of GSK3 genes was evaluated using a sub-program incorporated in the MCScanX software.

Cis-acting element analysis

The sequence of 2,000 bp upstream from the initiation codon in celery GSK3 genes was obtained through BLAST searching in the celery genome, and then was submitted to the PlantCARE website (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to evaluate the cis-acting elements. Venn diagrams were drawn to explore the significance of the corresponding elements.

Gene expression pattern analysis

The transcriptome sequencing data of celery was obtained from the Genome Sequence Archive database (GSA, https://ngdc.cncb.ac.cn/gsa) under the accession number CRA001996, which includes the expression data from the root, petiole, and leaf three specific tissues in celery. The expression of the AgSK genes were determined by Fragments Per Kilobase of transcript per Million fragments (FPKM) values and then the log2 transformed FPKM values were used to do the clustering and to generate the heatmaps with Multi-Experiment View.

Plant growth, abiotic stresses treatments and quantitative real-time PCR

Two-week-old celery seedlings were immersed in 1/2 Ms liquid medium containing 200 mM NaCl or 20% polyethylene glycol (PEG) under 22 °C for salt or drought stresses, respectively. Cold treatment was performed under 4 °C conditions. Samples were harvested at 0, 1, 6, and 12 h for all the treatments, and the total RNA was further extracted for each sample. Then qRT-PCR assays were conducted to assess the expression patterns of the AgSKs according to a standard protocol^[55]. A Real-Time PCR machine ABI-7900HT (Applied Biosystems, CA, US) and the SYBR Super Real PreMix Plus (TIANGEN Bio, Inc.) were employed in this study. Three biological replicates were performed for each timepoint in each treatment. The oligos used in this study are listed in Supplemental Table S10, and the AgTUB-B gene was used as the internal reference. For the germination assay, the sterilized seeds of AgSK22 transgenic plants and wild type were planted on 1/2 MS solid media supplemented with 0 or 150 mM NaCl, then the germination rate was calculated on a daily basis. Statistical analysis was conducted with the SPSS software. The various data sets were analyzed by one-way ANOVA (analysis of variance) with Duncan post hoc multiple comparisons.

Yeast two hybrid assay

AgBZR1-pGBKT7 and AgSK-pGADT7 vectors were co-transformed into the AH109 yeast cells. The colonies grown on SD/-Leu-Trp medium were screened on SD/-Leu-Trp-His medium to determine the interaction between AgBZR1 and AgSKs.

Subcellular localization analysis of AgSK proteins

Subcellular localization analysis was performed in Nicotiana benthamiana leaf epidermal cells according to a published protocol^[56].