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 1-week treatments. Compared with control, the leaves of wheat treated with 1 µM auxin didn’t show significant difference, however, the length of its’ seminal root was much shorter. This indicated that 1 µM auxin significantly inhibited root growth but had little effect on the growth of leaves. In response to 50 µM, significant 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 significant differences in the root length of wheat samples under the different treatments. The effects of these different concentrations of auxin are reflected by the wheat phenotype, and significant 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 identification of differentially expressed genes (DEGs)
The results showed that 9783 genes were classified 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 identified under 1uM auxin treatment (Additional file 1: TableS1). There were 6118 genes whose expression was downregulated and 2333 genes whose expression was upregulated in response to 50 µM (Additional file 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 significantly enriched in the 1 µM auxin treatment compared with the control, and catalytic activity, cell part and metabolic process were the most significantly 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 significantly 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 file 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 identified (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 identified in the two different treatments, with the exceptions of the insulin signalling and MAPK signalling pathway-plant pathways (Q-value < 0.05).
In total, 3279 DEGs were clustered into 8 profiles by Short Time-series Expression Miner software (STEM, Fig. 4a). The general trend diagram of all profiles is shown in Figure S1 (Additional file 4). Of the 3279 DEGs, 2831 were further clustered into 3 profiles (Fig. 4b-d, p-value ≤0.05), which involved two downregulated patterns (profile 1 and profile 0) and an upregulated pattern (profile 6). Profile 1 and profile 0 contained 1743 and 602 DEGs, respectively, while profile 6 contained 486 DEGs. Of the DEGs, 19.18% (629/3279) were annotated. The top 20 KEGG pathways with the highest representation of the DEGs are shown in Table 2. The ribosome (ko03010), starch and sucrose metabolism (ko00500), plant hormonesignal transduction (ko04075), phenylpropanoid biosynthesis (ko00940), photosynthesis (ko00195), and carbon metabolism (ko01200) pathways were significantly enriched. However in profile 6, no DEGs were annotated to ribosome pathway. The starch and sucrose metabolism pathway was the second most enriched pathway, and 5 DEGs that accounted for 7.69% of 65 DEGs in profile 6 were annotated to this pathway.
The wheat genome contains 176 TabZIP transcription factors.
A total of 468 sequences (including splice variants) were identified on the basis of the functional annotation (PFAM domains) of the TGAC wheat genome. This dataset was simplified by keeping only one splice variant from each genomic locus for further analyses ( Additional file 5: Table S5). A total of 185 sequences were identified 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 file 5). AtbZIP TFs have been classified into 10 groups based on the sequence similarity of basic regions and certain conserved motifs. Similarly, OsbZIP TFs have been classified into 11 groups (I–XI) on the basis of their amino acid sequences and DNA-binding sites within the bZIP basic domain. Domain analysis revealed that all but five 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 five, 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. 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 artificially eliminated, and a total of 176 TabZIP genes were ultimately determined.
Gene structure, motif composition, and phylogenetic classification 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 file 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 classification plays an important role in the evolution of multiple gene families.
Gene structural diversity is an vital component of gene family evolution and further supports phylogenetic groupings. The intron/exon organization and number of introns are typical imprints of gene evolution within gene families. The percentage of intronless TabZIP family members obtained in the present study was 18.75%, which was similar to that in rice (19.1%); 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 file 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. Motif 9 of groups 3, 7, and 8 indicates that the proteins have a longer Leu zipper. Most of the closely related members showed similar motifs in terms of alignment and position indicated that they may possess similar biological functions. Overall, the conserved motif compositions and similar gene structures of the TabZIP members within the same group, together with the phylogenetic analysis results, strongly support the reliability of the group classifications.
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. 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 file 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 findings 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 first selected 15 representative TabZIPs to enhance the reliability of the results by qPCR, but the amplification effect of 6 of them was not good, and finally 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 specifically 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.