Transcriptomic and Systematic Analyses Reveal TabZIP Family Members Involved in the Response to Different Concentrations of Auxin in Wheat Root

gene structure and conserved motifs among the different related to functional of TabZIP genes. TabZIPs perceive Transcriptomic studies of wheat roots were carried out to compare their physiological and transcriptional responses to different concentrations of auxin. Our results showed that different auxin concentrations had different effects on roots’ genes of wheat. The inhibitory effect of 1 µM auxin on roots manifested as an enrichment of the glycolysis/gluconeogenesis and glutathione metabolism pathways. Enrichment of the phenylpropanoid biosynthesis and plant-pathogen interaction pathways in response to 50 µM auxin indicated a stress response. Plant hormone signal transduction, starch and sucrose metabolism and photosynthesis were enriched in all the proles. Furthermore, TabZIP TFs in respond to different concentrations of auxin was also detected. The gene structure and phylogenetic classication analysis of TabZIPs further conrmed that the complexity of the regulatory machinery involved in different concentrations of auxin and cross-regulation of multicomponent signalling pathways. Determining which of these TabZIPs are the primary responses of different cincentrations auxin is necessary before the molecular basis can be fully described. Such genes would represent candidates for interfering in wheat root regulation.

The bZIP family is one of the largest and most diverse TF families, and its members are involved in seed germination, light signalling, hormone and sugar signalling [16,17]. bZIP proteins are characterized by a conserved bZIP domain, which is 60 to 80 amino acids in length and has two distinct structures [18,19]. The rst structure is a highly conserved basic DNA-binding region comprising 16 amino acid residues, and the most striking characteristic is the presence of an invariant N-x7-R/K motif. The other structure is a more diversi ed leucine zipper dimerization region. It has been thoroughly studied in a variety of biological systems. Humans have 56 bZIP proteins [20], whereas there are a total of 77 AtbZIP genes [21] and 89 OsbZIP genes [18]. The role of some bZIP genes in the regulation of speci c auxin response has been preliminarily studied. VIP1 is a bZIP TF in Arabidopsis which changes root cap structures and local auxin responses, and decreases the root vertical growth index (VGI) [22]. In addition, the HY5 gene, which contain bZIP domain, is involved not only in the light response but also in gravitropic and secondary thickening, in which auxin plays an important role [23,24]. OsbZIP46 exhibits an auxin response, suggesting a possible role of OsbZIP46 in mediating cross-talk between abscisic acid (ABA) and indoleacetic acid (IAA) in rice [25]. Although studies on the role of bZIP TFs in the regulation of auxin are limited, it would be of great importance to understand the various roles of bZIP family members and their structures for further investigation of the molecular mechanisms of auxin involved in root regulation.
Transcriptomic analysis is an effective method for global expression pro ling of genes involved in the response of auxin in plants, and compared with microarrays, RNA-seq can provide more information at a more affordable cost and is emerging as a powerful tool for transcriptomic analyses [26]. There have been more comprehensive studies on the time effect and tissue speci city of single concentrations of auxin. However, in our work, we focused on the effects of different auxin concentrations (1 µM and 50 µM) on transcriptomic differences in wheat roots in order to reveal the similarities and differences of molecular networks of wheat roots in response to different concentrations of auxin. Furthermore, we made a comprehensive bioinformatics analysis of the TabZIP genes from the wheat genome. Phylogenetic and conserved motif analyses were performed to reveal the differences among subgroups of TabZIP genes. Additionally, TabZIP genes involved in auxin response were found from the transcriptome data, as well as the expression differences of these auxin-responsive TabZIP genes under different concentrations were compared and analyzed. Subsequently, we discuss the relationship of these auxin-responsive TabZIP genes in gene structures and functional characteristics to reveal their potential role in the auxin regulatory network. The results of this study provided the basis for further study of TabZIP genes and provided some promising candidate genes for transgenic technology to improve the root system of wheat.

Phenotypic responses to different concentrations of auxin in wheat roots
The seedlings with the same growth rate (Fig. 1a) were selected for treatments with different concentrations of auxin (1 µM and 50 µM). Fig. 1b shows the changes in the appearance of wheat leaves and roots affter 1week treatments. Compared with control, the leaves of wheat treated with 1 µM auxin didn't show signi cant difference, however, the length of its' seminal root was much shorter. This indicated that 1 µM auxin signi cantly inhibited root growth but had little effect on the growth of leaves. In response to 50 µM, signi cant changes had taken place in both leaf and root tissues. The leaves and the length of the seminal roots were shorter than both 1 µM and contorl. The data in Fig. 1d show there are signi cant differences in the root length of wheat samples under the different treatments. The effects of these different concentrations of auxin are re ected by the wheat phenotype, and signi cant genetic changes took place at the same time.
Overview of libraries data sets by RNA-seq As shown in Table 1, a total of 100.78G nucleotides, equivalent to 686,382,894 raw reads and 678,304,326 clean reads were generated from the libraries, an average of 75 million reads (11.2 Gb) per sample. The size of data greatly exceeded the transcript size of the wheat genome, therefore, it was assumed to cover most transcribed genes within the wheat genome. The percentage of clean reads with a Q30 was greater than 89% for all nine of the samples, and the GC content ranged from 53.11 to 53.71%. After aligning the nine sets of clean reads to the reference genome (TAGCv1) individually, the percentage of reads mapping to the reference genome ranged from 88.31 to 92.36%, with uniquely matches were from 81.46 to 85.25%. In summary, the sequencing quality of the nine samples was high and suitable for follow-up analysis.
Overview of the data and identi cation of differentially expressed genes (DEGs) The results showed that 9783 genes were classi ed as DEGs. The expression of most of these genes was down-regulated in the two treatments (Fig. 1e). There were 308 and 326 genes whose expression was upregulated and downregulated, respectively, were identi ed under 1uM auxin treatment (Additional le 1: TableS1). There were 6118 genes whose expression was downregulated and 2333 genes whose expression was upregulated in response to 50 µM (Additional le 2: Table S2). As shown in Fig. 1e, there were several overlapping genes among the DEGs in the two treatments. We performed a Gene Ontology (GO) enrichment analysis to investigate the distribution of DEGs in terms of their involvement in biological process (BPs), cellular components (CCs) and molecular functions (MFs). Twenty-three GO terms (MF, 6 terms; CC, 6 terms;  BP, 11 terms) were signi cantly enriched in the 1 µM auxin treatment compared with the control, and catalytic  activity, cell part and metabolic process were the most signi cantly enriched terms in the MF, CC, and BP categories, respectively (Fig. 2a). In addition, as shown in Fig. 2b, the number of GO terms in the 50 µM auxin treatment increased to 34 (MF, 11 terms; CC, 9 terms; BP, 14 terms). The most signi cantly enriched terms in the three categories were the same between the two treatments.
To investigate the major pathways of the DEGs, we aligned all the DEGs to Kyoto Encylcopedia of Genes and Genomes (KEGG) pathways (Additional le 3: Table S3 and Table S4). In response to 1 µM auxin, the DEGs were enriched in 42 pathways, in which glycolysis/gluconeogenesis (6 DEGs) and glutathione metabolism (4 DEGs) were the top two pathways containing the greatest number of DEGs (Fig 3a. Q-value < 0.05). In response to 50 µM auxin, thirty-six KEGG pathways were identi ed (Fig. 3b) and the top three pathways were phenylpropanoid biosynthesis (56 DEGs), plant-pathogen interaction (45 DEGs) and photosynthesis (39 DEGs). As shown in Fig. 3a, b, different enriched pathways were identi ed in the two different treatments, with the exceptions of the insulin signalling and MAPK signalling pathway-plant pathways (Q-value < 0.05).
The starch and sucrose metabolism pathway was the second most enriched pathway, and 5 DEGs that accounted for 7.69% of 65 DEGs in pro le 6 were annotated to this pathway.
The wheat genome contains 176 TabZIP transcription factors.
A total of 468 sequences (including splice variants) were identi ed on the basis of the functional annotation (PFAM domains) of the TGAC wheat genome. This dataset was simpli ed by keeping only one splice variant from each genomic locus for further analyses ( Additional le 5: Table S5). A total of 185 sequences were identi ed by the SMART, PFAM and CDD programs. To facilitate the analysis and understanding of evolutionary trees, we used a consistent naming pattern for all TabZIP genes, taking into account their subfamily association, phylogenetic relationships and subgenomic location. Each gene name starts with an abbreviation for the species name, e.g., Triticum aestivum (Ta). The gene names include A, B or D, indicating the chromosomes on which they are located, as shown in Table S5 (Additional le 5). AtbZIP TFs have been classi ed into 10 groups based on the sequence similarity of basic regions and certain conserved motifs [16]. Similarly, OsbZIP TFs have been classi ed into 11 groups (I-XI) on the basis of their amino acid sequences and DNA-binding sites within the bZIP basic domain [18]. Domain analysis revealed that all but ve of the potential TabZIP TFs had a typical bZIP domain with an invariant N-x7-R/K motif within the basic region and a heptad repeat of Leu or other bulky hydrophobic amino acids (Ile, Val, Phe, or Met) positioned exactly nine amino acids upstream of R/K toward the C-terminus. Of the remaining ve, in the case of TabZIP4A_73, TabZIP4D_91, TabZIP5A_95 and TabZIP5D_119, the conserved Arg/Lys in the basic region was replaced with an Ile residue. In the case of TabZIP3A_51, the conserved Arg/Lys in the basic region was replaced with a Cys residue [27]. When generating the unrooted neighbour-joining phylogenetic tree via MEGA 7, the short sequences of TabZIP1A_1, TabZIP2B_29, TabZIP2B_30, and TabZIP3B_62 were arti cially eliminated, and a total of 176 TabZIP genes were ultimately determined.
Gene structure, motif composition, and phylogenetic classi cation analysis of TabZIP genes To analyse the evolutionary relationships among the bZIP family members in wheat, Arabidopsis and rice, 196 amino acid sequences were used to construct a phylogenetic tree. Recently, studies have reported that the bZIP TFs from Arabidopsis and rice are divided into 10 and 11 major clades, respectively, each of which contains more than 12 members (multi-gene clades), 4 to 10 members (medium-gene clades), or 1 to 3 members (oligo-gene clades) [16,28]. bZIP domains of 10 different AtbZIP genes and 10 different OsbZIP genes (Additional le 5: Table S6) from each of the subgroups were selected as representatives for further comparison [16,18]. All sequences information is shown in Table S6. Based on the phylogenetic relationships and gene structure information obtained here, TabZIP family can be divided into 9 groups (1-9, Fig. 5). Gene phylogenetic classi cation plays an important role in the evolution of multiple gene families [29].
Gene structural diversity is an vital component of gene family evolution and further supports phylogenetic groupings [30]. The intron/exon organization and number of introns are typical imprints of gene evolution within gene families [31]. The percentage of intronless TabZIP family members obtained in the present study was 18.75%, which was similar to that in rice (19.1% [18]); TabZIPU_185 in group 3 had the highest number of introns (9). Most TabZIP genes contained less than 5 introns (Fig. 7). Generally, closely clustered TabZIPs in the phylogenetic tree exhibited similar exon-intron structures, including those for exon length and intron numbers. The multi-gene clades included groups 1, 5, 7, 8, and 9. Most intronless TabZIPs clustered into group 8, except for TabZIP7D_174, whose motif position obviously differed from that of the other genes in the same group. The structure of genes in group 8 is more compact than that of genes in group 5, and the group 5 genes have one intron. Oligo-gene clades include groups 3 and 6, which contain seven to nine introns.
To further investigate the structural features of TabZIPs, the 196 bZIP protein sequences were analysed using the Multiple Expectation maximization for Motif Elicitation (MEME) tool, which revealed 10 conserved motifs ( Fig. 7). The features of these protein motifs are listed in Additional le 5: Table S7. Motif 2 was present within the TabZIPs in most groups, indicating evolutionary conservation. Compared with group 9, whose motifs are more numerous and complex, group 1 is characterized by the presence of motif 1. To further elucidate the phylogenetic mechanisms of TabZIP family, we constructed comparative syntenic maps of wheat together with Arabidopsis and rice (Fig. 6a). A total of 14 and 130 bZIP genes showed syntenic relationships with those in Arabidopsis and rice respectively, suggesting that TabZIPs were more closely allied with OsbZIPs than AtbZIPs and that the large-scale expansion of these families may not have occurred before the monocot-dicot split [32]. Importantly, related studies of OsbZIPs can also provide a theoretical basis for TabZIP family.

Differential expression of TabZIP genes in different auxin concentration treatments
The fragments per kilobase of transcript per million mapped reads (FPKM) values of the 22 TabZIP genes are listed in Additional le 5: Table S8, and a hierarchical clustering analysis and heatmap were generated to display the expression patterns of the TabZIP genes (Fig. 8a). There were only four TabZIP genes (TabZIP6B_141, TabZIP6D_147, TabZIP7A_158, and TabZIP7D_169) changed in response to 1 µM auxin, and the expression of them was downregulated. In addition, each has 3 introns. TabZIP7A_158 in group 9 didn't respond to 50 µM auxin, but the expression of TabZIP6B_141, TabZIP6D_147 and TabZIP7D_169, which contain a common motif 4 in group 4, remained downregulated in response to 50 µM auxin. The expression of the TabZIP6B_141 and TabZIP6D_147 parallel pair was downregulated to a great extent under 50 uM than under 1 µM, while TabZIP7D_169 was downregulated to a similar extent under the two concentrations. 18 of the 22 TabZIP genes did not respond to 1 µM auxin; their expression changed only in response to 50 µM auxin, with the expression of 4 and 14 upregulated and downregulated, respectively. With the exceptions of TabZIP1D_15, whose expression was upregulated, and TabZIP1A_5, whose expression was downregulated, all of the 18 TabZIP genes contain motif 2. The genes whose expression was upregulated were concentrated in groups 1, 5, and 6, but there were also genes whose expression was downregulated in group 5. The expression of TabZIP7B_161 and TabZIP7B_167 decreased the most. Moreover, the genes whose expression was downregulated showed a more complicated motif composition than did the genes whose expression was upregulated, and their intron number was more varied.
Validation of differential transcription using quantitative real-time PCR (qRT-PCR) To further verify the ndings from Illumina RNA-Seq results, we have performed 9 DEGs for their key roles in response to different concentrations of auxin by qRT-PCR. According to the response results of TabZIPs to different concentrations of auxin and the grouping of TabZIPs, we rst selected 15 representative TabZIPs to enhance the reliability of the results by qPCR, but the ampli cation effect of 6 of them was not good, and nally the results of 9 TabZIPs were obtained (Fig. 8b). Although TabZIP6B_141, TabZIP6D_147, TabZIP5B_118, TabZIP6A_134 and TabZIP7A_158 were down-regulated at different auxin concentrations, TabZIP6B_141 and TabZIP6D_147 were down-regualated much more under 1uM treatment than 50 uM treatment, TabZIP5B_118 was down-regualated much more under 50uM treatment than 1 uM treatment.
TabZIP5A_103, TabZIP5D_130 and TabZIPU_182 only responded speci cally under 50uM. Most of qPCR outcomes correlated closely with the transcriptomic results calculated from the RNA-Seq output indicated that these genes are indeed involved in the auxin response.

Discussion
Effects of inhibitory concentrations and stress concentrations of auxin on wheat roots Treatments with exogenous auxin can ablate root elongation, induce lateral root formation, and globally affect the levels of transcripts [33,34]. Our data concerning two auxin concentrations that showed 1uM auxin had an inhibitory effect on wheat roots, while 50 µM auxin caused a stress response. Many auxin-induced physiological or molecular changes that occurred in the two treatments differed, suggesting the importance of the concentration effects of auxin. Here, we analysed the inhibitory concentration (1 µM) and stressconcentration (50 µM) effects of auxin on the biology of and molecular processes in wheat roots.
With respect to inhibitory-concentration effects, we found that relatedly few genes responded to 1 µM auxin, and only 7.57% (48/634) of the genes were signi cantly enriched in KEGG pathways.
Glycolysis/gluconeogenesis (6 DEGs) and glutathione metabolism (4 DEGs) were the top two enriched KEGG pathways (Fig. 2c), indicating that these DEGs are relatively highly sensitive to the treatment with inhibitory concentrations. The expression of most of these 10 DEGs was downregulated by auxin treatment. At inhibitory concentrations, auxin inhibits H + -ATPases and causes alkalinization of cell walls. A relatively high pH of the apoplast decreases the activity of expansins and endotrans glucosylases/hydrolases and increases peroxidase activity, cross-linking hydroxyproline-rich glycoproteins and counteracting cell wall-loosening enzymes [35,36]. Many of these enzymes are associated with plant defence mechanisms in response to biotic or abiotic stress [37]. Moreover, genes involved in glutathione reaction are essential enzymes in plant stress resistance and protect tissues against oxidative stress, toxins, herbicide stress, and other stressors [38].
With respect to stress-concentration effects, we observed that high levels of auxin cause plant stress, which involves inducing the expression of defence-and detoxi cation-related genes to maintain auxin homeostasis.
Phenylpropanoid biosynthesis (56 DEGs), plant-pathogen interaction (45 DEGs) and photosynthesis (39 DEGs) were the top enriched KEGG pathways (Fig. 3b). The expression of most genes enriched in the plantpathogen interaction pathway was strongly upregulated in response to auxin treatment. This pathway in pro le 6 was the most enriched, which implies that there are speci c response mechanisms to very high auxin concentrations. These responses are possibly an adaptation to pathogen infection. In addition, the phenylpropanoid pathway is one of the pathways for lignin biosynthesis in plants, and its stimulation is considered a common feature of several abiotic stress responses [39]. Lignin is the main component of both the cell wall and apoplastic barriers and is the main structure in plant mechanical support and defence systems [40]. Therefore, it could be suggested that high concentrations of auxin may stimulate the phenylpropanoid pathway in wheat roots, leading to a shift in metabolism toward lignin.

Auxin signals interact with many pathways involved in wheat root growth
The plant root is a dynamic structure which is constantly modi ed in a temporal and spatial manner to respond to environmental stimuli including auxin [41]. Auxin had a relatively strong effect on the expression of genes annotated as participating in ribosome and hormone-dependent signalling pathways (Table 2).
Ribosomal proteins are known translational regulators of the auxin response and provide downstream regulation of lipid-based metabolic processes involved in tissue differentiation and transport [42]. This activity is essential in the developing meristem for the speci city of root stem cells. In addition, evidence for the cross-talk between auxin and light signalling [43] and sugars [44] is growing. An enrichment of genes involved in the response to photosynthesis was observed in response to treatment with auxin. Because of the complexity of the hormone-light interaction effect in root development, additional advances in our understanding of hormonal regulation of the root genome may bene t from the utilization of dark-grown root systems [45]. These investigations may nd direct changes related to belowground development processes by removing the changes in chloroplast-related processes. In Arabidopsis, mutations in two members of a family of 19 GH3related proteins provide further evidence for the proposed link between auxin signalling and photomorphogenesis [43,46]. Moreover, the quantitative level of GH3 transcription has been shown to be regulated by bZIP TFs via G-boxes. Our study implicates several core pathways and TabZIP genes as being involved in auxin stimulation, may revealing the dynamic regulatory process of auxin in the wheat root system.
Evolutionary relationships and gene structure analyses of the TabZIP gene family can provide information for the study of speci c functions bZIP TFs invole in the light response, the stress response and sucrose signalling, pathogen defence, seed maturation and ower development [17]. Within the DNA cis-acting elements preferentially bound by bZIP TFs, the G-box element is ubiquitous in genes induced by ABA, auxin, rosin and salicylic acid [17,47]. Here, we identi ed 176 TabZIP genes, which we assigned to 9 conserved subfamilies (Fig. 5). The gene structures and protein motifs were relatively conserved in the groups, which provided a valuable reference for TabZIPs' analysis and function (Fig. 7). Most of the paralogous pairs had the same exons and motifs; however, some pairs had different numbers of exons but had the same motifs (exons 2 and 5 in TabZIP7D_171 and TabZIPU_179, respectively, and exons 7 and 9 in TabZIPU_185 and TabZIPU_178, respectively). This may be due to the loss or gain of exons during gene evolution [32]. The differences in gene structure and additional conserved motifs among the different groups may be related to the functional diversity of TabZIP genes.

Functional characterization and elucidating the roles of auxin-responsive TabZIP genes in auxin regulatory networks
The distribution of TabZIPs that responded to different concentrations of auxin was not concentrated, indicating that they may play an key role in different mechanisms. Under 1 µM auxin, which exerted an inhibitory effect on the roots, there were only 4 TabZIPs whose expression was down-regulated, while under 50 µM auxin, the number of down-regulated TabZIPs increased, and there were 4 TabZIPs whose expression was up-regulated. By combining gene structures, phylogenetic and synteny analysis, new clues to the biological function of TabZIP genes could be inferred through comparison with those function-known bZIP genes from model plants. Of the 4 genes that responded to 1 µM auxin, 3 genes (TabZIP6B_141, TabZIP6D_147, and TabZIP7D_169) were from group 4 and classi ed with AtbZIP5G_H (or HY5 TFs).
TabZIP6B_141 and TabZIP6D_147 have also been described as homologues of HY5 and they are already sensitive to 1uM auxin (Fig. 8a). The structure of these proteins has a relative long leucine zipper region and contains two Asn residues. Interestingly, HY5 TFs are involved in the light-based regulation and likely play a role in a wide variety of stimulus responses and developmental processes in the hypocotyl and roots [23,48].
The expression of TabZIP5A_103 and TabZIP5D_130, which are classi ed together with AtbZIP4G_G and OsbZIP01g_10 (OSBZ8) in group 6, were upregulated under 50 µM auxin. In Arabidopsis, AtbZIP16 and AtbZIP68 in group G can form homodimers and heterodimers with GBF1, GBF2 and GBF3 to regulate lightresponsive promoters [49][50][51][52]. GBF1/GBF2 is involved in the regulation of plant morphogenesis in response to blue light, which promotes the expansion of cotyledons and the formation of lateral roots [53]. OsbZIP01g_10 was shown to have structural features typical of the GBF-type bZIP proteins [54]. Therefore, TabZIP genes may play an important role in root growth in response to auxin cross-talk with light signals. But we should also consider the in uence of hydroponics on the results, passive exposure of the roots to light may strongly activate photosynthesis in root tissues. TabZIPs in both group 5 and group 6 contain only motif 2, indicating that their genetic relationships are relatively close. However, the members in group 5 are diverse. The expression of TabZIP2D_38 was upregulated in response to 50 µM auxin, while the expression of TabZIPU_182, TabZIP2D_37 and TabZIP7B_161 in group 5 was downregulated. Moreover, the function of ATbZIP3G_B, which is in the same group, is unknown, indicating that there may be other unknown genes in this group, or the function of this protein may be similar to that of group 6 proteins. TabZIP2D_38, TabZIPU_182 and TabZIP7B_161 are in the same subgroup as OsbZIP02g_6. Therefore, there isn't enough evidence on which to speculate their function.
TabZIP genes that are highly homologous to AtbZIP5G_C in Arabidopsis group C and AtbZIP1G_S in groupS were divided into group 7 and group 8. The group S genes tend to heterodimerize within the subgroup or with group C bZIP genes in Arabidopsis, dimerisation is likely to ful l speci c functions in gene regulation [55,56].
Members of these groups have markedly expansive leucine zipper regions with up to nine heptapeptide repeats. AtbZIP11 which was in the C/S bZIP low-energy activated network gained control over the root meristem by directly activating IAA3/SHY2 transcription [57]. Morever, bZIP11 modulate Aux/IAA and GH3 expression by recruiting the histone acetylation (HAT) machinery to target promoters [58]. Additionaly, genes in group 8 contain almost no introns, whereas the expression of individual genes in groups 7 and 8 was downregulated in response to 50 µM auxin. In previous studies, genes with few or no introns were thought to display enhanced levels of expression in plants [59,60]. To respond in a timely manner to various stresses, genes with a compact genetic structure with relatively few introns would be rapidly activated [61]. These genes may be involved in different reaction processes, homologues of group S bZIPs are also transcriptionally activated after stress treatment (e.g. cold, drought, wounding) [62]. In rice, LIP19 and OsOBF which are classi ed with OsbZIP12g_8 interact in response to cold stress [63].
There are also some genes whose expression was downregulated under 50 µM auxin in group 9, and they may complement or replace TabZIP7A_158 in the same group, whose exprsesion was downregulated under 1 µM auxin. The TabZIPs in group 9 were homologous to those in Arabidopsis group A (AtbZIP2G_A) and OsbZIP08g_5 (TRAB1). 7 AtbZIP genes in group A (such as AtbZIP39/ABI5 and AtbZIP36/ABF2/AREB1) and OsbZIP08g_5 (TRAB1) were studied, the available functional information suggests they play roles in the response to ABA or stress signalling, including drought and salinity [54,[64][65][66]. TabZIP6A_134 was described as homologues of ABI5. Interestingly, ABI5 probably acts by recruiting the B1 domain of the ABI3 protein to the LEA promoter. ABI3 transcription may be regulated by auxin, and auxin signalling in lateral root development may have an ABI3-dependent component [67]. Moreover, ABI3 TFs could potentially interact with other ARF or IAA proteins to mediate auxin responses [68,69]. As discussed above, many TabZIPs probably have overlapping functions, Our identi cation of all TabZIP genes is thus a necessary prerequisite for the dissection of individual bZIP protein function. Further investigation of the molecular mechanisms underlying how TabZIPs perceive auxin stimulus-induced signals and how TabZIPs control the expression of target genes will unravel the role TabZIPs play in root development and the auxin signalling pathways in wheat.

Conclusion
Transcriptomic studies of wheat roots were carried out to compare their physiological and transcriptional responses to different concentrations of auxin. Our results showed that different auxin concentrations had different effects on roots' genes of wheat. The inhibitory effect of 1 µM auxin on roots manifested as an enrichment of the glycolysis/gluconeogenesis and glutathione metabolism pathways. Enrichment of the phenylpropanoid biosynthesis and plant-pathogen interaction pathways in response to 50 µM auxin indicated a stress response. Plant hormone signal transduction, starch and sucrose metabolism and photosynthesis were enriched in all the pro les. Furthermore, TabZIP TFs in respond to different concentrations of auxin was also detected. The gene structure and phylogenetic classi cation analysis of TabZIPs further con rmed that the complexity of the regulatory machinery involved in different concentrations of auxin and cross-regulation of multicomponent signalling pathways. Determining which of these TabZIPs are the primary responses of different cincentrations auxin is necessary before the molecular basis can be fully described. Such genes would represent candidates for interfering in wheat root regulation.

Plant Materials, Growth Conditions and Auxin Treatments
The wheat (Triticum aestivum) variety used in this study is "Chinese Spring" which is sequenced and used worldwide, and seeds are obtained from our Lab. The wheat seeds were surface sterilized with 1% NaClO solution for 15 min, rinsed three times with distilled water, kept at 4 ℃ for 2 days, and then transferred in to growth chamber (14/10 h day/night photoperiod under 24/17 ℃ day/night temperatures). After 1 days, 30 healthy seedlings with approximately equal root lengths were transferred to Hoaglands containing different concentrations of auxin (0, 1, 50 uM). After 1 week, images were photographed, and the root elongation was measured by NIH-Image software (https://imagej.nih.gov/ij/index.html). Meanwhile, Roots were collected for transcriptomic and qRT-PCR analysis.

mRNA sequencing by Illumina HiSeq
Total RNA of each sample was extracted using TRIzol Reagent (Invitrogen) / RNeasy Mini Kit (Qiagen) / other kits. Total RNA of each sample was quanti ed and quali ed by Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), NanoDrop (Thermo Fisher Scienti c Inc.) and 1% agrose gel.1 μg total RNA with RIN value above 7 was used for following library preparation. Next generation sequencing library preparations were constructed according to the manufacturer's protocol (NEBNext® Ultra™ RNA Library Prep Kit for Illumina®).
The poly(A) mRNA isolation was performed using NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB) or Ribo-Zero™ rRNA removal Kit (illumina). The mRNA fragmentation and priming was performed using NEBNext First Strand Synthesis Reaction Buffer and NEBNext Random Primers. First strand cDNA was synthesized using ProtoScript II Reverse Transcriptase and the second-strand cDNA was synthesized using Second Strand Synthesis Enzyme Mix. The puri ed double-stranded cDNA by AxyPrep Mag PCR Clean-up (Axygen) was then treated with End Prep Enzyme Mix to repair both ends and add a dA-tailing in one reaction, followed by a T-A ligation to add adaptors to both ends. Size selection of Adaptor-ligated DNA was then performed using AxyPrep Mag PCR Clean-up (Axygen), and fragments of ~360 bp (with the approximate insert size of 300 bp) were recovered. Each sample was then ampli ed by PCR for 11 cycles using P5 and P7 primers, with both primers carrying sequences which can anneal with ow cell to perform bridge PCR and P7 primer carrying a six-base index allowing for multiplexing. The PCR products were cleaned up using AxyPrep Mag PCR Clean-up (Axygen), validated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), and quanti ed by Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA).
Then libraries with different indices were multiplexed and loaded on an Illumina HiSeq instrument according to manufacturer's instructions (Illumina, San Diego, CA, USA). Sequencing was carried out using a 2x150 bp paired-end (PE) con guration; image analysis and base calling were conducted by the HiSeq Control Software (HCS) + OLB + GAPipeline-1.6 (Illumina) on the HiSeq instrument. The sequences were processed and analyzed by GENEWIZ.
De novo assembly of RNA-seq reads and quantifying gene expression In order to remove technical sequences, including adapters, polymerase chain reaction (PCR) primers, or fragments thereof, and quality of bases lower than 20, pass lter data of fastq format were processed by Trimmomatic (v0.30) to be high quality clean data. The obtained quality ltered reads were de novo assembled into contigs by the Trinity Program [70]. Unigenes were de ned after removing redundancy and short contigs from the assembly. Reference genome sequences and gene model annotation les of relative species were downloaded from genome website, such as UCSC, NCBI, ENSEMBL. Secondly, Hisat2 (v2.0.1) was used to index reference genome sequence. Finally, clean data were aligned to reference genome via software Hisat2 (v2.0.1). The number of unique-match reads was calculated and normalized to RPKM (reads per kb per million reads) for gene expression analysis. Differential expression analysis used the DESeq Bioconductor package, a model based on the negative binomial distribution. After adjusted by Benjamini and Hochberg's approach for controlling the false discovery rate, p-value of genes were setted <0.05 to detect differential expressed ones [71]. The differentially expressed genes (DEGs) between control and 1uM, or between control and 50uM, or between 1uM and 50 uM were restricted with FDR (false discovery rate)≤0.001 and the absolute value of log2 Ratio≥1. To examine the expression pro le of DEGs, the expression data υ (from Control, 1uM and 50uM treatment) were normalized to 0, log2 (υ1uM/υControl), log2 (υ50uM/υControl). and then clustered by Short Time-series Expression Miner software (STEM) [72]. GO-Term Finder was used identifying Gene Ontology (GO) terms that annotate a list of enriched genes with a signi cant p-value less than 0.05 and we used scripts in house to enrich signi cant differential expression gene in KEGG pathways.

Sequence retrieval and annotation of TabZIP TFs
The protein sequence datas of wheat were download from the Ensembl Plants (http://plants.ensembl.org/index.html). Functional annotations were ltered for PFAM (http://pfam.xfam.org/) [73] identi ers of the bZIP_1, bZIP_2 and bZIP_C (PF00170, PF07716, and PF12498), respectively. HMMER 3.0 (http://eddylab.org/software/hmmer/Userguide.pdf) was used to search the bZIP genes from wheat genome database. A total of 468 sequences were identi ed. The default parameters were adopted and the cutoff value was set to1.2e-09 to screen 230 sequences. Splice variants were excluded and only the rst variant was kept for further analysis, with seven exceptions. All candidate genes that may contain bZIP domain based on HMMER results were further examined by con rming the existence of the bZIP core sequences using the SMART, PFAM and CDD program. Finally, Each potential gene was then manually examined to ensure the invariant N-x7-R/K motif in the basic region and a heptad repeat of Leu or other bulky hydrophobic amino acids (Ile, Val, Phe, or Met) positioned exactly nine amino acids upstream of R/K toward the C terminus. There are 176 sequences were identi ed.

Phylogenetic Analysis and Classi cation
Presently, the AtbZIP TFs have been classi ed into 10 groups [16] and OsbZIP TFs have been classi ed into 11 groups (I-XI) [18]. Using JGI database website (https://phytozome.jgi.doe.gov/pz/portal.html ) search and download protein sequences of the bZIP TFs of rice and Arabidopsis. AtbZIP and OsbZIP from each of the groups or subgroups, were selected as representatives for the further comparison. Next, the bZIP protein sequences of the characterized bZIP proteins were used to create multiple protein sequence alignments using ClustalW with default parameters. MEGA7 [74] was used to generate the unrooted Neighbor-joining phylogenetic tree with the following parameters: 1000 times bootstrap test, poisson model and pairwise deletions. The classi cations of TabZIP proteins were based on the topology and bootstrap values of the phylogenetic tree.

Gene Structure and Conserved Motif Characterization
The exon-intron organization of TabZIP TFs were determined by comparing predicted coding sequences with their corresponding full-length sequences using the online program Gene Structure Display Server (GSDS: http://gsds.cbi.pku.edu.cn/) [75]. The MEME program (http://meme-suite.org/) for protein sequence analysis was used to identify conserved motifs in the identi ed TabZIP proteins [76]. The optimized parameters were employed as the following: the number of repetitions, any; the maximum number of motifs 10; and the optimum width of each motif, between 6 and 100 residues.
Heatmap of TabZIP TFs and the orthologous bZIP genes obtained from wheat, Arabidopsis and rice.
Expression levels as log2tpm and a heatmap was generated with MEV. Genes were clustered according to their expression using Hierarchical Clustering (metric:Pearson Uncentered, method: complete).Multiple Collinearity Scan toolkit (MCScanX) was adopted to analyze the gene duplication events, with the default parameters. To exhibit the synteny relationship of the orthologous bZIP genes obtained from wheat, Arabidopsis and rice, the syntenic analysis maps were constructed using python written by ourselves.
Sequence of AtbZIP and OsbZIP proteins were obtained based on the description in corresponding literatures and downloaded from the phytozome database.
qRT-PCR veri cation of differential expression Total RNA was isolated from roots of wheat subjected to the various concentration treatments described above. cDNA was synthesized from 1 μg total RNA using PrimeScript® (TaKaRa, Japan) and an oligo (dT) primer, according to the manufacturer's instructions. Primers for quantitative real time PCR (qRT-PCR) were designed by Beacon Designer7.9 software (Premier Biosoft, America) listed in Additional le 5: Table S9 and synthesized by Sangon Biotech (Shanghai) Co., Ltd. SYBR based qRT-PCR reactions (TaKaRa, Japan) were performed on ABI StepOnePlus. Each 20 μL qPCR contained 3 μl diluted cDNA, 1ul of each primer, and 10 μl SYBR qPCR master mix (Vazyme, China), and was following the reaction conditions: 95°C/1 min followed by 40  (https://www.ncbi.nlm.nih.gov/SRA /PRJNA644764/). Other datas generated or analyzed during this study are included in this article and its additional les.

Competing interests
The authors declare that the research was conducted in the absence of any commercial or nancial relationships that could be construed as a potential con ict of interest. Authors' Contributions JW , XT and YF designed the study, performed data analyses and wrote the manuscript. XT carried out the analysis of TabZIP genes, YF performed GO and KEGG pathway analysis, XT and ZJ performed qRT-PCR experiments, JW and ZW revised the manuscript. All authors have read and approved the nal manuscript. The numbers 1-3 after CK, 1uM and 50uM identify the three independent biological replicates for the control and auxin treatments, respectively.