Divergent Response Strategies of CsABF Facing Abiotic Stress in Tea Plant: Perspectives From Drought-Tolerance Studies

In plants, the bZIP family plays vital roles in various biological processes, including seed maturation, flower development, light signal transduction, pathogen defense, and various stress responses. Tea, as a popular beverage, is widely cultivated and has withstood a degree of environmental adversity. Currently, knowledge of the bZIP gene family in tea plants remains very limited. In this study, a total of 76 CsbZIP genes in tea plant were identified for the whole genome. Phylogenetic analysis with Arabidopsis counterparts revealed that CsbZIP proteins clustered into 13 subgroups, among which 13 ABFs related to the ABA signaling transduction pathway were further identified by conserved motif alignment and named CsABF1-13, these belonged to the A and S subgroups of CsbZIP and had close evolutionary relationships, possessing uniform or similar motif compositions. Transcriptome analysis revealed the expression profiles of CsABF genes in different tissues (bud, young leaf, mature leaf, old leaf, stem, root, flower, and fruit) and under diverse environmental stresses (drought, salt, chilling, and MeJA). Several CsABF genes with relatively low tissue expression, including CsABF1, CsABF5, CsABF9, and CsABF10, showed strong expression induction in stress response. Thirteen CsABF genes, were examined by qRT-PCR in two tea plant cultivars, drought-tolerant “Taicha 12” and drought-sensitive “Fuyun 6”, under exogenous ABA and drought stress. Furthermore, CsABF2, CsABF8, and CsABF11, were screened out as key transcription factors regulating drought tolerance of tea cultivars. Subsequently, some potential target genes regulated by CsABFs were screened by co-expression network and enrichment analysis. This study update CsbZIP gene family and provides a global survey of the ABF gene family in tea plant. The resolution of the molecular mechanism of drought resistance in different varieties could be helpful for improving stress resistance in tea plant via genetic engineering.


INTRODUCTION
Tea plants [Camellia sinensis (L.) O. Kuntze] originated in the southwest region of China and have been cultivated and utilized for thousands of years. Today, tea is widely grown worldwide (Reygaert, 2014). The tea plant is sensitive to environmental conditions and various adverse environmental factors, such as extreme temperature, drought, and salinity, can affect the growth of tea in the wild (Wang W. et al., 2016). Due to global climate change, the frequency of extreme weather such as high temperature and drought, has increased in recent years, causing serious harm to tea production. The study of the environmental tolerance of tea plants is helpful for improving the yield and quality of tea leaves. The potential stress resistance of plants is usually determined by the expression of stress-induced genes regulated by specific transcription factors (TFs) (Yamaguchi-Shinozaki and Shinozaki, 2006). TFs are widespread in plants and regulate the expression of target genes by binding directly or indirectly with specific cis-regulatory elements.
The basic (region) leucine zippers (bZIPs) are a class of TFs that are widely present in the life cycle of plants and play an important regulatory role, participating in many important biological processes, such as seed maturation, flower development, light signaling, pathogen defense, and response to various biotic and abiotic stresses . With respect to the number of TFs in the whole genome, the size of the bZIP-family varies considerably among species, with 78 in Arabidopsis thaliana (Droege-Laser et al., 2018), 91 in Oryza sativa (Liu et al., 2018), 247 in Brassica napus (Zhou et al., 2017), and 191 in Triticum aestivum (Agarwal et al., 2019).
Abscisic acid responsive element-binding factors (ABFs) are the most prominent subgroup in the A-bZIP family and have been found to act at the core of ABA signaling, which implements adaptive responses to counteract water deficits in vegetative tissues in response to abiotic stress, such as drought, salinity or extreme temperature (Droege-Laser et al., 2018). ABFs are regarded as a transit station for transmitting ABA signals to downstream functional genes. In short, members of the clade A serine/threonine protein phosphatases 2Cs (PP2Cs) have been found to negatively control ABA signaling by inactivating the SnRK2 (SNF1-related kinases 2) via dephosphorylation. ABA is captured by PYR/PYL/RCAR (pyrabactin resistance proteins/PYR-like proteins/regulatory component of the ABA receptor) co-receptors and forms ternary complexes with PP2C, thereby preventing their interaction with SnRK2 kinases (Fujii et al., 2009;Umezawa et al., 2010). SnRK2s directly phosphorylate five conserved serine/threonine kinase phosphorylation sites (RXXS/T) in ABFs to strongly enhance their transactivation properties (Takashi et al., 2006). Furthermore, ABFs regulate the expression of the ABA signaling pathway and other downstream stress-responsive target genes via direct binding to ABA responsive element (ABRE: ACGTGG/TC) cis-elements (Choi et al., 2000).
In recent years, several ABFs involved in plant response to drought stress have been identified. MEABL5, a transcription factor belonging to the A-clade of the bZIP family in Manihot esculenta Crantz, was found to specifically bind to the ABRE ciselement in the promoter of the major cell wall invertase gene MeWINV3 in vitro and in vivo to positively regulate drought tolerance (Liu et al., 2019). In Brassica oleracea (Zhou et al., 2013), BolABI5 was induced by ABA or drought stress and combined with ABRE cis-acting element in the promoter to resist drought stress. Overexpression of AtABF3 in Medicago sativa can reduce the transpiration rate and accumulation of reactive oxygen species (ROS) and improve drought tolerance (Wang Z. et al., 2016). In Ipomoea batatas, overexpression of IbABF4 can promote the root elongation of Arabidopsis seedlings under drought stress . FtbZIP5 found in Fagopyrum tataricum can enhance the drought and salt resistance in transgenic Arabidopsis, mainly by regulating the accumulation of ROS and inducing the plant antioxidant system. In addition, FtbZIP83 can respond to drought stress by up-regulating the transcription abundance of downstream ABA-induced genes, and overexpression of transgenic Arabidopsis was found to improve drought tolerance by reducing oxidative damage Li Q. et al., 2020).
We used the genomic data to identify the 76 bZIPs and then analyzed the expression level of these genes using RNAseq publicly available, 15 more than those reported by Hou et al. (2018). A phylogenetic tree was constructed with the Arabidopsis bZIP gene family as a reference. Then, based on the analysis of the evolutionary relationship and conserved motifs, 13 ABF genes related to the ABA signal transduction pathway were screened for transcriptome analysis. The expression patterns of CsABFs in eight different tea tissues and organs under different environmental stresses (drought, salt, chilling, and MeJA) were analyzed by transcriptome data analysis. The expression differences of CsABFs between drought-sensitive "Fuyun 6" and drought-tolerant "Taicha 12" were further compared in order to screen key TFs regulating the drought resistance of C. sinensis and to identify potential regulatory candidate genes using co-expression network analysis. This study provides a new concept of the structure and function of the ABF gene in tea plants, valuable resources for improving the stress resistance of tea plants, and a basis for improving the quality of tea.

Plant Materials and Abiotic Stress Treatment
Two-year-old tea cutting seedlings (Fuyun 6# and Taicha 12#) were obtained from Qianhe Tea Seedling Co., Ltd. in Quanzhou, Fujian Province. Seedlings with a height of 30-35 cm, complete root systems, and consistent growth were screened for domestication culture. Transition culture was performed with clean water and 25% Shigeki Konish nutrient solution every 7 days. Finally, the seedlings were transferred to the phytotron for hydroponic culture, and the conditions were as follows: illumination duration of 16 h light and 8 h dark alternately, room temperature of 23 • C ± 1, air humidity of 45-50%, ventilation at intervals of 12 h, and the nutrient solution was replaced every 3 days.
After the tea cutting seedlings grew new fibrous roots, 0.25 mM ABA was added to the culture solution. The second and third leaves were sampled after treatment for 0, 0.5, 6, and 24 h and then immediately placed in liquid nitrogen. In addition, 25% PEG-6000 was added to the culture solution to perform the drought treatment, the second and third leaves were sampled after treatment for 0, 1, and 3 days. The control group was cultured in normal nutrient solution. Total samples were stored at −80 • C for qRT-PCR detection. Three biological replicates of each sample were taken for subsequent testing.

Chlorophyll Fluorescence Detection
The leaves of Fuyun 6# and Taicha 12# treated with drought for 0, 1, 3, and 5 days were measured by a portable pulse-modulated chlorophyll fluorescence analyzer (PAM 2500, Zealquest Scientific Technology Co., Ltd., Shanghai, China). Five sample leaves were tested in each treatment, and the results were averaged. First, the leaf was clamped by a dark adapted leaf clip. After 30 min of dark adaptation, the leaf was measured with IMA-min/B and blue at 450 nm. Then the initial fluorescence (F o ), maximum fluorescence (F m ), variable fluorescence (F v ), and photosystem II (PSII) maximum light energy conversion rate (F v /F m ) were obtained. The calculation formula used was as follows:

Total RNA Extraction and qRT-PCR Analysis
Total RNA was extracted by the RNA-prep pure plant kit (polysaccharide and polyphenolic rich, TIANGEN, Shanghai, China), and reverse transcription was performed using the Goldenstar TM RT6 cDNA Synthesis Kit Ver.2 (TSINGKE, Beijing, China). The CsGAPDH gene was applied as an internal reference for real-time quantitative PCR (qRT-PCR) analysis. Suitable primers were designed for all selected genes by Primer plus 3.0 1 , and only primers with a single melting curve and the amplification efficiency between 95-105% were selected for subsequent qRT-PCR analysis. All screened genespecific primers are listed in Supplementary Table 2. The 2 × TsingKe Master qPCRmixKit (TsingKe, Beijing, China) was used for qRT-PCR on the QuanStudio 6-Flex fluorescence quantitative PCR platform (Thermo Fisher Scientific, Singapore). Every sample was tested in three biological replicates and three technical replicates, which showed good consistency, and the relative expression level was calculated using the 2 − CT method.

Identification and Biological Analysis of CsbZIP and CsABF
To screen the bZIP genes in the Camellia sinensis genome, the reported AtbZIP sequences were downloaded from the Arabidopsis gene annotation database TAIR 2 . Based on 78 AtbZIP protein sequences, a blast search was performed in the C. sinensis genome annotation database of TPIA 3 (Xia et al., 2019) to obtain homologous CsbZIP coding sequences, genomic sequences and peptide sequences. Seventy-six annotated bZIP protein sequences were obtained. The integrity of the bZIP domain of each protein sequence was verified by SMART 4 and PFAM 5 and confirmed as a member of the bZIP gene family. Through the KEGG pathway on the TPIA website, 13 ABF transcription factors related to ABA signal transduction pathway were identified and named CsABF1-13 (Supplementary Table 1).
To investigate the phylogenetic relationship between the CsABF and CsbZIP gene families in tea plant, 78 bZIP in Arabidopsis were applied to construct the phylogenetic tree with 76 bZIP in C. sinensis. The phylogenetic tree was generated and displayed by the Mega 7.0 program with the following parameters: neighbor-joining (NJ) method, P-distance model, pairwise deletion, and 1000 bootstrap replicates. Multiple sequence alignments of CsABF were performed by DNAMAN. The Multiple Em for Motif Elicitation (MEME) suite (version 4.10.1) 6 was used to identify the conserved motifs in all CsABFs. The parameters were as follows: maximum number of motifs: 6; number of repeats: any; maximum number of motifs: 50; optimal motif width: 6-200 amino acid residues.

Transcriptome Data Processing and Expression Pattern Analysis
The transcriptome data of the tissue-specific expression of CsABFs were downloaded from the TPIA website, and the eight organs included apical bud, young leaf, mature leaf, old leaf, stem, root, flower, and fruit in appropriate seasons. The transcriptome data of CsABFs in response to drought and salinity stress (Zhang et al., 2017), chilling stress , and MeJA treatment (Shi et al., 2017) were obtained from NCBI. The data acquisition and homogenization processing methods refer to Li M. et al. (2020) Table 3). Multiple Experiment Viewer 4.9 software was used for heatmap clustering analysis.

Construction of Co-expression Network and Prediction of Functional Genes
Co-expression analysis was applied in the tea plant drought stress transcriptome to construct the network related to CsABFs, which were then intersected with the genes screened from TeaCon website 7 to narrow down the target range. Co-expressed genes were screened with the corresponding Pearson correlation coefficients (PCC-value) of more than 0.7 and then were used to construct the co-expressed network of CsABFs with Cytoscape 3.6.1 (Shannon et al., 2003). Some genes co-expressed by key transcription factors were selected for functional enrichment analysis in order to identify potential regulatory target genes (Supplementary Table 4). GO enrichment analysis was carried out by R language 3.6.3 (Bell Laboratories).

Data Source and Statistical Analysis
One-way analysis of variance (ANOVA) was performed on the experimental data using SPSS 22 (IBM, United States), and multiple comparisons were performed using Duncan test to analyze the significance of the difference between the two treatments at the 0.05 level.

Identification and Phylogenetic Analysis of CsbZIP Family Members in Tea Plant
According to the reported AtbZIP genes as queries, 76 bZIP genes were identified from the latest assembled version database of the tea plant information archive (TPIA 8 ). The integrity of all the bZIP domains was determined using the online tool SMART (see Text Footnote 4). The relevant information is listed in Supplementary Table 1.
The bZIP transcription factor gene family is huge and complex, but all members contain the typical leucine zipper region N-X7-R/K-X9-L-X6-L-X6-L (Jakoby et al., 2002). Droege-Laser et al. (2018) divided 78 members of Arabidopsis bZIP gene family into the 13 subfamilies of A, B, C, D, E, F, G, H, I, J, K, M, and S according to the basic domain and other conserved domains. Subgroup A is mainly involved in ABA and stress signal transduction, subgroup C is mainly involved in seed development and pathogen defense, subgroup D is involved in disease defense and physiological growth, and subgroup H plays a key role in photosynthesis; subfamily S is the largest bZIP subfamily in Arabidopsis and is involved in sugar signaling and stress pathways (Jakoby et al., 2002). To elucidate the phylogenetic relationships of the bZIP family in tea plant, an unrooted tree was established based on the aligned 76 CsbZIPs with 78 AtbZIP protein sequences (Figure 1). The results indicated that most of the obtained subgroups were consistent with previous phylogenetic analyses (Hou et al., 2018). As shown in Figure 1, the phylogenetic tree clustered all of the CsbZIP members into 13 subgroups similarly to AtbZIPs. In addition, thirteen ABF genes related to ABA signal transduction were identified and named CsABF1-13. Noteworthy, CsABF1-13 were mainly clustered in the A subgroup (specifically CsABF1/2/3/5/6/7/8/11/12/13), meanwhile CsABF10 and CsABF4/9 clustered in subgroup K and S, respectively (blue diamond). Most of the genes in the A and S subgroups of the ABA-activated signaling pathway are ABF (ABRE binding factor) or AREB (ABA response element binding protein) proteins, which can bind to different ABRE-containing promoters (Choi et al., 2000;Uno et al., 2000). Under all stress conditions, the differentially expressed genes were enriched in the A and S subgroups, indicating that these two subgroup members might act as key factors in stress responses .

Gene Structure and Expression Profile of the CsABF Gene Family
The function of ABA in plants depends on signal perception and transduction. Existing research has shown that when plants suffer stress due to adversity, especially drought stress, the body increases ABA synthesis levels, senses, and transmits signals 8 http://tpdb.shengxin.ren/ through the PYL-PP2C-SnRK2 complex. ABF transcription factors are activated to regulate the expression of downstream genes, thus enabling plants to adapt to harsh environments (Fahad et al., 2017). In order to better explore the response strategies of tea plant to drought stress, we focused on the ABFs, which are closely related to ABA signal transduction in the bZIP family. The diversity of plant protein sequences and structures generated in biological evolution is the possible mechanism for the formation of polygene families and functional diversity. The characteristics of plant diversity lead to the efficient utilization of natural resources or adaptation to adverse environments (Shang et al., 2013). Therefore, we performed multiple sequence alignment for 13 full-length CsABF protein sequences using DNAMAN and predicted the conserved motifs and their phylogenetic relationships with MEME in order to further understand the structural features of the ABF proteins in tea plant.
Phylogenetic analyses of CsABFs were performed (Figure 2A), and corresponding genetic structure analyses by the GSDS website were used to reveal the exon and intron structures of CsABFs ( Figure 2C). It's worth noting that the CsABF3, CsABF4, and CsABF11 genes did not have any introns in their genomic sequences. Other CsABF genes contain anywhere from one to multiple introns each. These differences indicate that CsABFs are complex in evolution, and have different expression regulation modes.
Six conserved motifs were identified in CsABFs with widths ranging from 6 to 50, and the locus distribution of each sequence was zero or one. The height of each letter represents the conservatism of a specific amino acid in each motif ( Figure 2B). Similar to other plants, Motif1 is typical in bZIP domains detected by PFAM. The basic domain is close to the N-terminal of the leucine zipper domain and acts as a nuclear localization signal by binding the fixed N-X7-R/K structure with a specific DNA sequence (Landschulz et al., 1988). Most of the CsABF protein sequences share common motifs, but the conserved domains of CsABF4, CsABF9, and CsABF10 are obviously different from others ( Figure 2B), which may be the reason why they do not belong to the A subfamily of CsbZIP, which also means that they have genetic functional diversity. This result can provide a reference for further studies on functional differentiation among bZIP subfamilies.
Based on TPIA transcriptome data, the differential expression profile of CsABF genes in eight tissues was analyzed. The relative expression of all genes in a single tissue is shown in Figure 2D. The expression levels of CsABF12 were higher in old leaf, those of CsABF2, CsABF4, and CsABF8 were highly expressed in root, and CsABF2, CsABF8, and CsABF12 were highly expressed in fruit, which suggests that the CsABFs have functional diversification. There are also some gene family members with very low expression that may be regulated by different developmental processes and stresses. When plants suffer from drought stress, the root is the first to feel a water deficit, and the root stem cell niche, meristem and vascular system will coordinate their response to drought. Taking into account that CsABF2 and CsABF8 are both expressed in roots and fruits, it may be worthwhile to investigate their roles in the regulation of drought resistance in tea plants.

Expression Pattern of CsABFs Under Different Abiotic Stress
In order to analyze the different response strategies of CsABFs to different abiotic stresses, transcriptome data of tea plant under different stress treatments were downloaded from NCBI, including drought, salinity, chilling, and MeJA treatment. After homogenizing the RPKM values of corresponding unigenes, a heatmap was drawn for further analysis. From Figure 3, it can be seen that divergent response strategies of CsABFs to face abiotic stresses. The expression levels of CsABF2, CsABF8, and CsABF10 were significantly up-regulated from 12 to 24 h after MeJA treatment and then decreased to normal levels. CsABF1, CsABF5, and CsABF9 were significantly down-regulated within 12 h after treatment, and the downregulation of CsABF5 expression persisted for a long time, indicating the negative regulation of MeJA signaling. Under drought treatment, the expression levels of CsABF2, CsABF8, CsABF9, and CsABF13 were significantly up-regulated, while the expression levels of CsABF4, CsABF7, and CsABF10 were decreased. Under cold stress treatment, the expression levels (D) Tissue-specific expression. Different colors represent different organs, and the circle size represents the relative expression level. YL, young leaves (1st and 2nd leaves); ML, mature leaves (3rd and 4th leaves); OL, old leaves. The relative expression data were log2 transformed and diagram was drawn using the Multiple Array Viewer software. of CsABF6 and CsABF7 were significantly up-regulated during full acclimation (CA1) and de-acclimation (CA3), while CsABF1 and CsABF5 showed no significant changes in CA1, but were increased in CA3. The expression of CsABF10 showed significant down-regulation in both CA1 and CA3. Under salt treatment, the expression levels of CsABF1, CsABF6 and CsABF9 were increased, and CsABF6 reached the highest expression level at 48 h. The expression level of CsABF7 obviously decreased. The results showed that there was no redundancy in the stress responses at transcriptional level of the CsABF family members, and no single gene could deal with more than two stresses. Based on our analysis, we consider that each gene show a specific answer to the stresses imposed and they might regulate downstream genes to deal with a complex environment.

Effect of Exogenous ABA and Drought Stress on CsABFs Expression in Different Tea Cultivars
According to earlier study, the drought-sensitive variety "Fuyun 6" and drought-tolerant variety "Taicha 12" were screened from many tea varieties (Zhang et al., 2019). In this experiment, they were selected for further studies on the changes of ABF gene family expression patterns during drought and ABA stress treatments. By phenotypic observation and chlorophyll fluorescence monitoring, it can be seen that after 3 days of drought treatment, the leaves of "Fu Yun 6" showed an obvious dehydration phenotype and gradually dried up, while the dehydration of leaves of "Taicha 12" was only observed on the fifth day of drought stress ( Figure 4A). F m represents the maximum fluorescence value that can be absorbed by plants , and its change is affected by environment . In addition, F v /F m is also an important index to study the effects of various environmental stresses on photosynthesis (Li et al., 2014). As shown in Figure 4B, Fm of both tea varieties showed the same slight decreasing trend under drought stress. However, from the third day of drought treatment, F v /F m dropped sharply in "Fuyun 6" but had no significant change in "Taicha 12" indicating that "Taicha 12" had significantly higher drought resistance than "Fuyun 6". This difference was also evident in the leaf fluorescence diagram. The results indicate that the drought resistance of the two tea varieties does exist, providing powerful test materials for the subsequent investigation of molecular mechanisms of drought regulation.
In order to analyze the molecular mechanism of different tea varieties in response to drought stress, and to explore the key factors involved in the regulation of the CsABF genes family, drought-sensitive variety "Fuyun 6" and drought-tolerant variety "Taicha 12" were selected as experimental materials to undergo exogenous ABA and drought stress treatment. The expression patterns of CsABFs in different tea cultivars were detected by qRT-PCR. As can be seen from Figure 5A, all family members except CsABF1, CsABF4, and CsABF12 showed a significant upregulation trend after ABA treatment. The difference between the two tea cultivars is that CsABF7, CsABF8, CsABF9, CsABF10, and CsABF11 induced significantly higher up-regulated expression in "Taicha 12" than in "Fuyun 6", while CsABF3 showed a much higher expression level in "Fuyun 6" than in "Taicha 12". These results suggested that some CsABF gene family members exhibited completely different expression patterns in two different tea cultivars, implying that they may be potentially key resistance factors by regulating downstream genes and generating different resistance responses. The response model of CsABF family members to drought stress was very similar to the results after exogenous ABA treatment, and the expression changes of some genes showed a certain lag effect considering that the signaling conduction induced by exogenous ABA was significantly faster than the endogenous signaling conduction induced by drought stress.  Figure 5B, the expression levels of CsABF2, CsABF8, CsABF9, CsABF10, and CsABF11 reached the highest point within 6-24 h in "Taicha 12" after exogenous ABA treatment and then decreased. They did not change significantly at first in drought-treated "Taicha 12" until their expression levels reached the highest point 72 h after treatment, indicating that they both responded to ABA signal transduction caused by drought. However, they changed less in "Fuyun 6". It is thus speculated that when "Taicha 12" is subjected to drought stress, the expression of these genes will be activated to initiate downstream functional genes to exercise power and then promote the enhancement of drought tolerance in tea plant, which is the key for "Taicha 12" to become a drought-tolerant variety. Therefore, the key transcription factors from these genes for tea plants to cope with drought were screened for followup studies.

Co-expression Network
The gene co-expression network can be divided into a condition-dependent network and a condition-independent network according to the experimental treatment (Proost and Mutwil, 2016). Condition-dependent gene networks are useful for revealing the relationship between genes under a specific condition such as biotic or abiotic stress. Weighted gene coexpression networks have proven to be a valuable tool for revealing unknown gene relationships in condition-dependent models (Zhang et al., 2020).
The co-expression genes of CsABF2, CsABF8, CsABF9, and CsABF11 were screened out and network visualization was constructed as shown in Figure 6A. According to the same PCC-value, we screened 32, 19, 30, and 10 functional genes that were highly correlated with CsABF2, CsABF8, CsABF9, and CsABF11, respectively. The results of the co-expression network showed that although these four ABF transcription factors should stimulate the drought tolerance of "Taicha12" varieties, they were independent of each other, and there were no functional genes that acted synergistically. This indicates that tea plants have a complex molecular regulation mechanism to cope with drought stress. Among these co-expression networks, some of the genes are potentially stress-responsive genes regulated by ABF transcription factors, which need to be characterized by further studies. As shown in Figure 6B, among these candidate genes, CsABF2 may be involved in regulating transport and transduction processes, including intracellular protein transport, oligopeptide transport, and chloride transport, and may also be involved in brassinosteroid FIGURE 6 | (A) Co-expression network analysis of key TFs. The core loci are CsABF2, CsABF8, CsABF9, and CsABF11. (B) GO enrichment analysis table of candidate genes regulated by TF CsABF2/8/9/11. GO analysis will return a P-value for each GO with a differential gene, and a small P-value means that the differential gene is enriched in the GO. mediated signal regulation; CsABF8 may be involved in the regulation of energy metabolism processes such as the cation transport and the polysaccharide catabolic process; CsABF9 may be involved in the amino-glycan catabolic process; CsABF11 may be involved in regulating DNA-templated transcription. Meanwhile, most of the downstream target genes screened were hypothetical proteins that had not yet been named or verified, which posed a challenge for the subsequent verification work.

DISCUSSION
With the development of bioinformatics technology, the genome sequencing and assembly results for tea plant have become increasingly complete. Transcriptome sequencing has great potential in the discovery and identification of new genes. The bZIP transcription factor family in tea plant has been identified, with 18 genes (Cao et al., 2015) identified at the beginning, increasing to 61 genes (Hou et al., 2018), and finally to 76 genes in the present study. This promoted our exploration of the regulatory functions of bZIP-TFs. As a perennial woody plant, tea has wide distribution and a rich variety of resources, and it has developed unique adaptive environmental regulation mechanisms during evolution (Shen et al., 2019). Among them, the molecular mechanism of CsABF regulating the drought resistance of tea plant is worthy of further exploration. In this study, a phylogenetic tree comparison of CsbZIP and its corresponding homologous AtbZIP sequences was performed to explore the CsbZIP genes family at the genome-wide level, and the subfamily of the CsABF TFs involved in the ABA signal transduction pathway was identified. The transcriptome data of the CsABF genes in diversified stress treatments were analyzed, as well as the qRT-PCR data of gene expression patterns in drought-sensitive "Fuyun 6" and drought-tolerant "Taicha 12" under drought stress. CsABF2, CsABF8, CsABF9, and CsABF11 were screened out as key TFs in response to drought stress, and their potential downstream regulatory target genes were screened out by the co-expression network.

Co-expression Network Provides Powerful Support for Exploring the New Function Genes
Co-expression network analysis can reveal potential associations between genes and regulatory factors, which undoubtedly helps us to search for potential downstream regulatory genes of CsABF and uncover the association between them to co-create the response network of drought stress. All candidate genes in Figure 6 were further validated by qRT-PCR to confirm that their expression patterns were indeed consistent with the corresponding ABF TFs. Therefore, mining the functions of these genes will help us understand how woody plants differ from herbaceous plants in responding to drought stress. Sequence alignment of the candidate genes screened by the coexpression network revealed that there were 32 highly correlated functional genes with CsABF2, among which TEA000304.1 belonged to the NRT1/PTR gene family, and the function of homologous gene AtNPF4.6 as a hormone transporter was verified in A. thaliana. AtNPF4.6 was initially identified as nitrate transporter NRT1.2, but it was also identified as an ABA transporter (Chiba et al., 2015). AtNPF4.6 has been proven to transport ABA to guard cells, thereby closing stomata to reduce water loss caused by transpiration and improving the drought resistance of Arabidopsis. Three other members of the AtNPF gene family, AtNPF4.1, AtNPF4.2, and AtNPF4.5, have also been shown to transport ABA in yeast (Kanno et al., 2012). Therefore, it is proposed that TEA000304.1 is regulated by CsABF2 and plays a role in the process of drought stress resistance in the "Taicha 12" cultivar.
Among the highly relevant candidate genes by CsABF8, TEA008554.1 was found to have high homology with β-amylase I (BamI) in Arabidopsis. AtBamI can participate in the daily degradation of starch in guard cells to maintain stomatal openings. Under drought stress, AtBamI can also generate a carbon skeleton to support proline biosynthesis through the daily degradation of instantaneous starch, so as to cope with the osmotic pressure imbalance caused by drought and ensure the normal growth of plants. Therefore, we speculate that TEA008554.1 is involved in regulating the physiological process of starch degradation in tea leaves under drought treatment.
Plants under drought stress can produce osmoprotective agents, such as sucrose, proline, LEA protein and antioxidant enzymes (Szabados and Savouré, 2010). Under the protection of osmotic agents, protein denaturation, membrane damage and reactive oxygen damage were reduced in plants, so as to maintain normal cell turgor and improve the drought resistance of plants. This study revealed that many putative proteins were involved in the process of CsABF regulating drought resistance. It was speculated that CsABF2 and CsABF8 were mainly responsible for defense by regulating polysaccharide metabolism and transporter proteins, while CsABF11 was involved in stress response mainly through the regulation of the DNA transcription process. However, the function of these unknown genes remains to be demonstrated in follow-up research.

Drought-Tolerant Cultivars Provide a Perspective for Studying Tea Plant Response Strategies to Abiotic Stress
The diversification of plant resources has allowed them to evolve different strategies in response to environmental stress. Many species have been screened for drought-sensitive and drought-tolerant varieties to study the differences in molecular regulation mechanisms. According to current research reports, screening is mainly based on plant phenotypes, physiological data, and changes in metabolite content in vivo, and there is no consistent molecular label that can be used to distinguish drought-sensitive and drought-tolerant varieties. In T. aestivum, the drought-tolerant variety "Luhan 7" can resist drought by responding to light-signaling pathways, while the drought-sensitive "Yangmai 16" variety can improve the yield through pre-drought induction to increase post-flowering drought stress tolerance during the vegetative period (Abid et al., 2016). In Solanum lycopersicum, drought resistance can be improved by reducing chlorophyll concentration and light absorption by chloroplasts (Sánchez-Rodríguez et al., 2012). In Foeniculum vulgare Mill, the abundance of glycoldegradation related proteins in the drought-sensitive type decreased in response to drought stress, while those in the drought-tolerant type increased in response to drought stress (Khodadadi et al., 2017).
The ABA signal transducing pathway is an important pathway of drought stress signal transmission. As shown in Figure 7, there may be regulatory factors highly responsive to ABA signal in different tea varieties, which are involved in regulating downstream functional genes to build stress resistance network to deal with the change of environment. Zhang et al. (2019) identified drought-sensitive and drought-tolerant varieties in tea plant. We conducted research on the molecular mechanism of ABA signaling in response to drought in tea pant in the early stage, but the drought-resistance mechanism of downstream transcription factors remains to be further studied. There are many reports on the involvement of different TFs in drought stress regulation. For example, in T. aestivum, TaCIPK27 partly plays a positive regulatory role to improve the drought resistance of wheat through the ABA-dependent pathway under drought stress . In A. thaliana, AtWRKY46 can independently regulate osmotic stress response and stomatal movement to resist drought stress. In Gossypium hirsutum, GbMYB5 could increase proline content, antioxidant enzyme activity and malondialdehyde (MDA) content to resist drought stress (Chen et al., 2015). In O. sativa, ONAC022 plays a positive role in drought stress tolerance through modulating an ABA-mediated pathway (Hong et al., 2016). OsNAC14 mediates drought tolerance by recruiting factors involved in DNA damage repair and defense responses, thereby improving drought tolerance in rice (Shim et al., 2018). In addition, the transcription factor OsNF-YA10 enhanced drought tolerance in rice . There are many reports on the involvement of different transcription factors in the regulation of drought stress, but whether transcription factors are the key factors for the formation of different drought-resistant varieties of tea remains to be further studied.
In this study, by constructing the potential regulatory network of CsABFs in different varieties (drought-sensitive and droughttolerant), we gained a better understanding of the role of TFs in key varieties of crops and provided a theoretical basis for breeding drought-resistant varieties at the molecular level.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
JL conducted the experiments and drafted the manuscript. JD analyzed the transcriptome data. LT and ML conducted part of the qRT-PCR experiment. XZ provided the selected experimental materials. SZ directed the construction of transcription factor co-expression network. XW designed the experiments. QC directed the experiments and revised the manuscript. All authors have read and approved the final manuscript.