Genome-wide identification of bZIP transcription factors and their responses to abiotic stress in celery

Abstract Celery (Apium graveolens L.) is one of the most important vegetables in the Apiaceae family, rich in nutrients and widely grown around the world. bZIP transcription factors family plays an important role in the transcription regulation of plant growth and development, as well as adaptation to the external environment. In this paper, 62 bZIP family transcription factors were screened and identified based on the whole genome sequence of celery. The bZIP proteins of celery and Arabidopsis thaliana were divided into 10 subfamilies according to the phylogenetic tree. Phylogenetic and evolutionary analysis showed that the number of bZIP family members gradually expanded from lower plants to higher plants during the long evolution process. Based on the homology of celery and A. thaliana bZIP genes, the interaction network between celery bZIP transcription factors and other proteins in the genome was constructed, and the correlation data of protein interaction were also obtained. The expression profiles of 12 selected AgbZIP genes were detected and analyzed under abiotic stress treatments and different tissues using RT-qPCR. The results showed that AgbZIP can respond to high temperature, low temperature, drought, and high salt stress.


Introduction
The environment that has a negative influence on plant growth and development is called stress. Several environmental factors adversely affect plant growth and development and final yield performance of a crop. Soil salinization, low and high temperatures, and drought are among the major environmental constraints to crop productivity worldwide [1]. The bZIP transcription factors is one of the largest regulatory families of transcription factors and plays important roles in physiological processes in plants, such as plant biology, abiotic stress response, growth, and development [2]. At present, a large number of bZIP transcription factors have been identified in eukaryotes. Many members of the bZIP family have been identified in plants, including 75 bZIP genes in Arabidopsis (A. thaliana) [3], 125 in maize (Zea mays) [4], 121 in banana (Musa paradisiaca) [5], 55 in grapevine (Vitis vinifera) [6], and 116 in apple (Malus domestica) [7].
TFs are a group of key regulatory proteins that play important roles in controlling the expression of signal response genes [8]. The DNA-binding domain of plant transcription factors contains amino acid residues that are in contact with the DNA bases of the cis-acting elements, which determine the specificity of the protein [9,10]. Depending on the specific domain, transcription factors can be divided into several families, such as MYB, MADS, WRKY, bZIP, etc. Members in the same family are highly conserved in certain regions of the amino acid sequence [11]. bZIP family usually includes a leucine zipper dimerisation domain and an alkaline domain that binds to DNA [12]. It is named according to its common bZIP conserved domain, which contains 60-80 amino acids [13,14]. The bZIP transcription factors form a homo-or heterodimer and bind to a specific gene promoter through their basic region, thereby regulating the expression of related downstream genes [14].
Celery (Apium graveolens L.) is one of the most important leaf vegetables in Apiaceae family. Celery is cultivated worldwide and utilized in food and cosmetic industries because it is an excellent source of vitamins, phenolic compounds, volatile oils, and other nutrients [15]. In this study, 62 bZIP family transcription factors were identified in celery, and then the structural domain and phylogenetic relationship were analyzed in detail. The expression patterns of some selected AgbZIP genes were detected under abiotic stress treatments. This study provides a reference to further investigate the function of celery bZIP transcription factors and the improvement of celery breeding.

Plant materials and growth conditions
The celery cultivar 'Jinnan shiqin' was used as the experimental material. This variety was selected from Jinnan district, Tianjin city, China, with developed root system, large leaf area, low crude fibre content, and strong biotic/abiotic resistance [16]. Plants were grown in an artificial climate chamber at State Key Laboratory of Crop Genetics and Germplasm Enhancement, Nanjing Agricultural University. Twomonth-old seedlings were subjected to high temperature (38 C), low temperature (4 C), salt (0.2 molÁL À1 NaCl), and drought (200 gÁL À1 polyethylene glycol (PEG)) conditions. Plant samples were collected at 0, 1, 2, 4, 8, and 24 h after treatment. The whole treated plants were harvested. To study the expression patterns of AgbZIP genes in different tissues, samples of roots, petioles and leaf blades were collected from growing plants under normal conditions. The samples were frozen in liquid nitrogen, and stored at À80 C for further study.

Sequence database searches
The sequence of Arabidopsis bZIP transcription factors were acquired from Arabidopsis information resource (http://www.arabidopsis.org/). Celery data were downloaded from CeleryDB, which is a genomic database of celery [17].

Motif recognition and phylogenetic analysis
Motifs of the selected genes were analyzed using MEME (Version 4.9.1) [18]. Multiple alignments of bZIP protein sequences were performed by Clustal X 1.83 software [19] and the phylogenetic tree was constructed with MEGA 5.0 [20]. The protein interaction network was built with STRING software [21].

RNA isolation and RT-qPCR analysis
Total RNA was extracted using the RNA kit (Beijing Tiangen, China) on the basis of Manufacturer's instructions. RNA from each of the treated samples was reverse transcribed into cDNA, using the Prime Script RT reagent kit (TaKaRa, Dalian, China). Primers for selected AgbZIP genes were designed using Primer Premier 5.0 Software. Celery internal reference gene was selected as the internal control [22]. We provide the sequences of all the primers in supplementary material Table S1. Primers were synthesized by Nanjing Genscript Corporation (Nanjing China). The RT-qPCR assay was performed using the SYBR Premix Ex Taq on a CFX96 Real-time PCR system (Bio-Rad). The reaction conditions were as follows: 95 C for 30 s followed by 40 cycles of 95 C for 5 s and 60 C for 30 s for annealing and extension. The experiments were performed with three independent biological replicates, and the relative expression level was calculated using the 2 ÀDDCt method [23].

Phylogenetic analysis of the bZIP transcription factors family in celery
In order to analyze the bZIP transcription factors family in celery, we constructed a phylogenetic tree using the amino acid sequences of 68 proteins from A. thaliana and 62 proteins from celery. According to the classification of A. thaliana bZIP family [3], the phylogenetic tree clustered all bZIP members into 10 groups (A-H) ( Figure 1). Based on this phylogenetic tree, we analyzed the composition of each subfamily (group) of the bZIP in celery. This included Group A, the largest group in the family (23%) with 15 members (Figure 2).

Motif analysis of the celery bZIP transcription factors family
The conservative motifs of 62 bZIP transcription factors in celery were analyzed by MEME software. The symbols of these motifs and the composition of each transcription factor in the bZIP family are shown in Figures 3 and 4, respectively. Each bZIP contains motif 1. Members of Groups H, C, B, and E all contain motif 1 and motif 3. The module composition of each family is similar. For example, motif 2 only appears in the D family. The D family members include the motifs 1, 2, 3, 5, 7, and 8. It also proves the adjacent relationship between the same families in the evolutionary relationship.

Evolution of bZIP transcription factors among different plant species
To analyze the evolution of the bZIP transcription factors family in plants, we searched 11 plants bZIP   During evolution, plant genomes must be altered to adapt to changes in the environment and become more complex. The development and growth of plants depend on the proper regulation of several genes [24]. The study found that bZIP transcription factors play a crucial role in these regulatory processes [25]. As more and more plants are genome-wide sequenced, the bZIP transcription factors family has been identified across the genome-wide range of many plants. In this study, we selected a total of 11 plants for the comparative analysis of bZIP transcription factors family. The number of family members varies greatly from plant to plant, with only 19 members in the alga Chlamydomonas reinhaurdtii. The number of bZIP family members was significantly higher than that in the algae, suggesting that the large-scale expansion of the bZIP family occurred after the evolution of higher plants.

Interaction network analysis of AgbZIP and related proteins
To study the relationship between AgbZIP and related proteins, all AgbZIP proteins were investigated in an Arabidopsis association model in STRING software to identify the functional interactions ( Figure 6). A total of 62 AgbZIP proteins are correlated with each other.
As shown in Figure 5, the combined score between bZIP60 and bZIP17 protein was the highest at 0.985. It shows that the correlation between Agr23923, Agr23022, and Agr29095 is very strong, indicating that these proteins can work together to regulate some biological processes. Agr23022 and Agr29095 belong to Group B together. bZIP17 and bZIP60 play an important role in the resistance of higher plants to high salt stress [26]. bZIP17 is an ER-membrane anchored transcription factors that responds to salt stress. In response to salt treatment, bZIP17 moves to Golgi apparatus and is processed by S1P and S2P protease to produce a soluble form of transcription factors, which is then transferred to the nucleus [27]. bZIP60 is a transcription factor involved in UPR signalling and appears to be involved in plant responses to salt stress [28]. This suggests that the AgbZIP transcriptions may be involved in the defence of celery against abiotic stress.
Transcription factors play an important role in plant growth and development and in response to the external environment through their interaction with DNA and other proteins. In view of that interaction between bZIP transcription factors and other proteins in Arabidopsis, we constructed the network map of the interactions between the celery bZIP transcription factors and the other proteins of the celery genome. The results showed that the total 62 celery bZIP transcription factors and 30 celery proteins constitute complex network regulation, which illustrates that the bZIP transcription factors of celery may participate in multiple biological processes. These results provide a feasible basis for studying the relationship between celery bZIP transcription factors and protein, and for further specific studies on how celery bZIP transcription factors interacts with other genes, so as to understand the molecular mechanism of celery bZIP transcription factors in the process of plant growth and development and its regulation under abiotic stress.

RT-qPCR analysis of selected AgbZIP genes under abiotic stress
Environmental stress can adversely affect plant growth and productivity and trigger a series of pathological, physiological, biochemical, and molecular changes  [29]. Related studies have shown that the bZIP transcription factors in A. thaliana can regulate many downstream salt tolerance and drought tolerance related genes [30]. bZIP plays an important role in plant response to salt damage, drought, cold damage, mechanical damage, and osmotic stress [31]. In order to clarify the expression pattern of AgbZIP genes under different abiotic stresses, the expression profiles of AgbZIP genes (Agr35790, Agr40452, Agr31176, Agr18552, Agr10824, Agr05322, Agr00336, Agr09027, Agr24495, Agr24735, Agr24967, and Agr38129) were detected under abiotic stress treatments, including heat, cold, drought, salinity (Figure 7).
High-temperature treatment: After high-temperature treatment, 12 AgbZIP genes showed different expression profiles. Agr18552, Agr31176, Agr09027, and Agr38129 all showed a significant increase in relative expression levels at 2 h and 24 h after high-temperature treatment, and all reached peak relative expression levels after 24 h of treatment. Low-temperature treatment: Under the condition of cold treatment, the relative expression levels of Agr05322, Agr35790, and Agr10824 decreased slightly compared with control after 1 h and 2 h treatment, but began to rise slightly after 4 h of treatment, and the relative expression level increased significantly after 8 h and 24 h, and reached the peak at 24 h. The relative expression levels of Agr40452 and Agr31176 also decreased slightly after 1 h and 2 h of treatment, and the relative expression level increased significantly and reached the peak after 4 h of treatment, but decreased again after 8 h and 24 h of treatment. Agr00336, Agr09027, and Agr24967 were upregulated significantly compared with control after 1 h and 2 h of treatment. The relative expression levels of Agr09027 and Agr24967 reached the peak after 1 h of treatment, but began to decrease after 4 h of treatment.
Drought treatment: After drought treatment, the relative expression levels of Agr31176, Agr05322, Agr35790, and Agr10824 genes were decreased compared with control. Agr18552, Agr09027, and Agr38129 were upregulated significantly compared with control after 1 h of treatment, reached the peak at 1 h, but decreased again after 2 h of treatment.
Salt treatment: After salt treatment, 12 AgbZIP genes showed different expression profiles. The relative expression levels of Agr05322 and Agr10824 decreased in the early stage of treatment, increased significantly after treatment of 8 h and 24 h, and reached the peak after 8 h of treatment. The relative expression level of Agr35790 and Agr24495 in each period under salt treatment did not change significantly compared with control, and remained at a stable level. Agr18552 and Agr24967 were upregulated compared with control at each period under salt treatment.
In this experiment, the results of RT-qPCR indicate that the AgbZIP genes can be up-regulated or downregulated in response to various abiotic stresses, and their relative expression levels were different in roots, petioles, and leaf blades. It indicates that AgbZIP can respond to high temperature, low temperature, drought, and high salt stress. In view of the important function of bZIP transcription factors in plant stress response, gene engineering can be used to overexpress bZIP or regulate the expression level of specific bZIP in crops. It can improve the multi-stress resistance and comprehensive quality of crops for cultivation of new varieties of good crops.

Gene expression analysis of celery roots, petioles, and leaf blades
As natural properties of each gene, tissue-specific expression is closely associated with the growth and development of the particular tissue [32]. In order to clarify the expression pattern of AgbZIP genes in celery tissues (roots, petioles, and leaf blades), the expression profiles of AgbZIP genes (Agr35790, Agr40452, Agr31176, Agr18552, Agr10824, Agr05322, Agr00336, Agr09027, Agr24495, Agr24735, Agr24967, and Agr38129) in celery roots, petioles, and leaf blades were detected ( Figure 8).
In the roots, petioles, and leaf blades of celery, 12 AgbZIP genes showed different expression profiles. Agr35790, Agr40452, Agr10824, Agr24495, Agr00336, and Agr24735 were expressed significantly higher in the petioles compared with the roots and leaf blades. The relative expression levels of Agr24967 and Agr38129 were significantly higher in the leaf blades compared with the roots and petioles. The relative expression levels of Agr31176, Agr05322, Agr18552, and Agr09027 were significantly higher in the root compared with the petioles and leaf blades.

cis-regulatory elements analysis of selected AgbZIP genes
The cis-regulatory elements in 7 AgbZIP genes (Agr18552, Agr10824, Agr00336, Agr09027, Agr24495, Agr24735, and Agr38129) promoters were analyzed using the online software PlantCARE based on the celery genome database, and the 9 most common elements are represented in Figure 9. These elements included a cis-acting regulatory element essential for the anaerobic induction (ARE), a cis-acting element involved in low-temperature responsiveness (LTR), a cis-acting regulatory element involved in seed-specific regulation (ABRE), a motif named CAAT-box associated with promoter and enhancer regions of common cisacting elements, and five light responsive elements (Box 4, G-box, G-Box, GA-motif, and MRE). The frequency of MRE, MYB, and MYC components was higher in the promoter of seven genes. Most AgbZIP genes contained more than one cis-element in their promoter regions.
The bZIP protein is one of the largest transcription factor families in eukaryotic cells [33]. It is important in plant response to abiotic stress [34]. At present, the study of bZIP transcription factors has made rapid progress, and bZIP has been found in A. thaliana [35], Capsicum annuum [36], Brassica rapa pekinensis [37], Glycine max [38], Zea mays [39], Helianthus annuus [40], and Solanum lycopersicum [41]. In the family of transcription factors, members of the same subfamily have the same or similar domains and may function similarly in binding to cis-elements and interacting with other proteins [42].

Conclusions
In this study, 62 AgbZIP transcription factors were identified from the celery genome database. According to the classification of A. thaliana bZIP family, the phylogenetic tree clustered all bZIP members into 10 groups. The physicochemical properties of AgbZIP were predicted by bioinformatics analysis. The results showed that the AgbZIP of the same subgroup had similar physicochemical properties, and the number of bZIP family transcription factors increased gradually from lower plants to higher plants. RT-qPCR results showed that AgbZIP genes could respond to various abiotic stresses. Our study can provide evidence for studying the regulatory mechanism of bZIP transcription factors in celery.

Disclosure statement
No potential conflict of interest was reported by the authors.

Compliance with ethical standards
This article does not contain any studies with human participants or animals performed by any of the authors.

Notes on contributor
XAS and YQQ initiated and designed the research; YQQ, FK, XZS, DAQ and LJX performed the experiments; YQQ, FK and XZS analyzed the data; XAS contributed reagents/ materials/analysis tools; YQQ wrote the paper; XAS and XZS revised the paper. All authors read and approved the final manuscript.