A total of 21 FBA genes from Triticum aestivum L. (TaFBA) were identified by BLAST analysis against the URGI (https://urgi.versailles.inra.fr) and Ensembl Plants (http://plants.ensembl.org) databases using the gene sequences of previously annotated FBA genes in wheat (GenBank No. FJ403591.1, AJ420778.1, FJ625793.3, KR139996.1, KR139997.1) and Arabidopsis (GenBank No: AT2G21330, At4G38970, AT2G01140, AT5G03690, AT4G26530, AT2G36460, AT4G26520, AT3G52930, AT1G18270) with E-value < 1e⁻¹⁰. Using these 21 TaFBA genes as queries, another 31 FBA genes from T. monococcum (6 genes), T. urartu (7 genes), A. speltoides (5 genes), A. sharonensis (6 genes), and A. tauschii (7 genes) were further identified from the URGI database. Using the genome sequences of these FBA genes downloaded from the URGI database as queries, the FBA cDNA sequences of T. aestivum, T. urartu, A. sharonensis, and A. tauschii were obtained from the WheatExp database (http://wheat.pw.usda.gov/WheatExp/) and NCBI's TSA and SRA databases. The FBA cDNA sequences of T. monococcum and A. speltoides were predicted by homology analysis.

We used URGI, Ensembl Plants databases, wheat_2014_90K SNP gene chip, and Chinese Spring_Deletion gene chip to predict the chromosomal localization of the TaFBA genes. To further confirm our results, gDNA PCR amplification of 21 nulli-tetrasomic lines (NA-TD, NB-TD, and ND-TB) of T. aestivum L.var. Chinese Spring was performed. The gene-specific primers used in PCR analysis are listed in Table S2. Nulli-tetrasomic Chinese Spring genotypes contain four copies of one chromosome pair (tetrasomic) to compensate for the lack of a homoeologous chromosome pair (nullisomic) (Sears et al., 1966). For instance, Nulli1A-Tetra1D (N1AT1D) lacks chromosome 1A (nullisomic-1A), but contains two chromosome pairs of chromosome 1D (tetrasomic-1D), resulting in a line containing no copies of chromosome 1A, two doses of chromosome 1B, and four doses of chromosome 1D. PCR analysis of genomic DNA from nulli-tetrasomic lines allows one to determine which chromosome gene is located on based on failure to amplify specific nullisomic lines (Ohnishi et al., 2008). PCR amplification was performed as follows: 94°C for 5 min; 94°C for 30 s, annealing at 60°C for 30 s, and extension at 72°C for 1 min, 35 cycles; and 72°C for 10 min for the final extension. PCR reactions were performed using the proofreading enzymes LA Taq or PrimeSTAR® HS DNA Polymerase (TaKaRa, Japan) depending on the application.

21 TaFBA genes (including 1.5-kb promoter regions) were cloned by genomic PCR of T. aestivum L. var. Chinese Spring. PCR protocols were the same as that used in chromosomal localization described in 2.2. The promoter sequences of two TaFBA genes (TaFBA6 and TaFBA12) were obtained by high-efficiency thermal asymmetric interlaced PCR (hiTAIL-PCR) (Liu and Chen, 2007) due to the incomplete sequencing of the wheat genome. The primers used in this study are listed in Table S2. The PCR products were purified using an AxyPrep DNA gel extraction kit (Axygen, USA) and ligated into the pMD19-T vector (TaKaRa, Japan). The recombinant plasmid DNA clones were sequenced by Invitrogen (USA). The generated TaFBA gene sequences were then submitted to GenBank of NCBI (Accession No. were listed in Table 1).

Multiple sequences alignment of 21 TaFBA genes and 31 FBA genes from wheat relatives were performed by MUSCLE (http://www.ebi.ac.uk/Tools/msa/muscle/). A rooted UPGMA phylogenetic tree was constructed with MEGA 6 under default parameters and bootstrapping with 1,000 replicates (Tamura et al., 2011). Genome sequences, cDNA sequences and protein primary sequences of FBAs were aligned by DNAMAN 6.0 software with pairwise method and were clustered by MeV 4.9 using the average linkage hierarchical clustering (HCL) method.

The exon-intron structure of TaFBA genes was determined by aligning the cDNA sequences to the corresponding genomic sequences using Gene Structure Display Server (GSDS, http://gsds.cbi.pku.edu.cn). The transposable elements of TaFBA genes were predicted on TREP (http://botserv2.uzh.ch/kelldata/trep-db/blast/blastTREP.html) and LTR_FINDER (http://tlife.fudan.edu.cn/ltr_finder/). Putative cis-acting regulatory DNA elements in the promoter sequences (1.5-kb upstream of the 5′-UTR) of TaFBA genes were annotated with PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/), and displayed in Argo Genome Browser v1.0.31 (http://www.broadinstitute.org/annotation/argo/). The exon-intron structure of the TaFBA genes was determined by aligning the cDNA sequences to the corresponding genomic sequences using the Gene Structure Display Server (GSDS, http://gsds.cbi.pku.edu.cn). The transposable elements of the TaFBA genes were predicted by using TREP (http://botserv2.uzh.ch/kelldata/trep-db/blast/blastTREP.html) and LTR_FINDER (http://tlife.fudan.edu.cn/ltr_finder/). Putative cis-acting regulatory DNA elements in the promoter region (1.5-kb upstream of the 5'-UTR) of the TaFBA genes were annotated with PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) and displayed in the Argo Genome Browser v1.0.31 (http://www.broadinstitute.org/annotation/argo/).

The physicochemical characteristics of the protein sequence of each FBA such as the number of amino acids (aa), molecular weight (MW), isoelectric point (pI), and grand average of hydropathy (GRAVY) value were calculated with the ProtParam tool in ExPASy (http://web.expasy.org/protparam/). The subcellular location of each FBA was predicted by using WoLF PSORT (Horton et al., 2007). Functional domains and conserved motifs were annotated with the Pfam database (Finn et al., 2016) and the MEME website (http://meme-suite.org/tools/meme).

A total of 245 FBA proteins of archaea, bacteria, and eukaryotes were identified by PHMMER searches on the Ensembl Plants database and key word searches of “fructose-1,6-bisphosphate aldolase” in the UniprotKB/Swiss-Prot databases. HMMER 3.0 was employed to perform an HMM search (Zhang and Wood, 2003), with the family-specific Glycolytic (PF00274) and F_bP_aldolase domain (PF01116).

Amino acid sequence alignment of the FBAs was performed by using the Clustal Omega method (Sievers et al., 2011), and a rooted maximum likelihood (ML) phylogenetic tree was constructed using MEGA 6. The phylogenetic tree files were edited by Evolview (http://www.evolgenius.info/evolview/). We used the I-TASSER program (Zhang, 2008) and Swiss-PdbViewer v4.1.0 (Schwede et al., 2003) to remodel and analyze the tertiary structure, and SWISS-MODEL to construct the quaternary structure of the TaFBAs.

To analyze the expression profiles of the TaFBA genes, the microarray data of the TaFBA genes were downloaded from PLEXdb (Dash et al., 2012) and the RNA-seq data from wheat transcriptome database WheatExp. Combined with the analysis using Genevestigator (Hruz et al., 2008), the expression profiles of different tissues and developmental stages of the TaFBA genes were generated by using OriginPro 9.1.

Seeds of Chinese Spring were sterilized with 75% alcohol for 30 s and 15% sodium hypochlorite for 5 min, rinsed 5 times, placed on moistened filter paper in Petri dishes, and germinated in the dark for 2 days. Then, the wheat seedlings were cultivated in Hoagland liquid culture medium at 23°C and 16 h light/8 h dark cycles. Treatment of the wheat seedlings was performed as described Zeng et al. (2016), with minor modifications. After 10 days of growth, salinity (200 mM NaCl), drought (15 % PEG 6000), ABA (50 μM ABA), heat (40°C), and cold (4°C) treatments were conducted. The shoots and roots of the seedlings were collected at 0, 1, 2, 6, 12, 24, and 48 h after treatment.

We investigated the influence of light on the expression of the TaFBA genes as described by Zhao et al. (2014). The etiolated wheat seedlings were cultured under continuous light, whereas the normal-cultured seedlings were placed in the dark after being exposure to light for 24 h. Wheat leaves were collected at 1, 2, 4, and 8 h after light or dark treatment, immediately frozen in liquid nitrogen, and kept at −80°C until RNA isolation. Next, total RNA were isolated by using the RNAiso plus reagent (TaKaRa, Japan) and treated with gDNA Eraser (TaKaRa, Japan) for 2 min at 42°C to degrade any residual genomic DNA. Then, first-strand cDNAs were synthesized from the total RNA by using a PrimeScript™ RT reagent kit (TaKaRa, Japan) according to the manufacturer's instructions. Finally, qRT- PCR was performed in optical 96-well plates (VIOX, UK) with a CFX96 Touch Real-Time PCR Detection System (BIO-RAD, USA) by using the SYBR Green method. The wheat ADP ribosylation factor (ADP-RF) gene was used as reference gene (Paolacci et al., 2009) to normalize the expression data. The thermal cycle protocol was set up as follows: 95°C for 30 s; 95°C for 5 s, and 60°C for 30 s for 40 cycles. The gene specific primers used for qRT-PCR were listed in Table S8. The quantitative analysis was accomplished with the 2^−ΔΔCT method. PCR was performed based on three technical replicates for each of the biological duplicates.

To investigate the TaFBA gene families, BLAST analysis with FBA genes of wheat and Arabidopsis thaliana (GenBank Accession Nos. are shown in 2.1) was performed using the IWGSC database. A total of 21 TaFBA genes (named as TaFBA1~21) were identified and cloned (Data sheet DOC1, TXT2). The cDNA sequences (Data Sheet TXT1) were obtained from NCBI and the WheatExp Database.

We aligned the cDNA and predicted protein sequences of the TaFBA genes by using DNAMAN with the pairwise method (Table S1). Phylogenetic and sequence analyses revealed that the TaFBA genes could be classified into two subfamilies, classes I and II. There were 18 TaFBA genes in class I, with TaFBA1~3, TaFBA4~6, TaFBA7~9, and TaFBA10~16 showing high identity (>90%) respectively. Meanwhile, three class II TaFBA genes (TaFBA19~21) were highly homologous and showed lower identity (< 30%) with those of class I.

The molecular characteristics of the TaFBA genes and their subcellular location are listed in Table 1. The predicted MW of class I TaFBA proteins ranged from 25.93 to 55.39 kDa, with TaFBA17 apparently being the smallest polypeptide, whereas class II TaFBAs were about 148 kDa, bigger than class I. WoLF PSORT was used to predict the subcellular location of 21 putative TaFBA proteins. Similar to Arabidopsis, rice, and tomato, 18 class I TaFBAs were located in two compartments: TaFBA1~9 were in the chloroplast/plastid and were thus designated as CpFBAs, whereas TaFBA10~18 were localized to the cytoplasm, and hereby named as cFBA. Three class II proteins TaFBA19~21 might also be located in the chloroplast/plastid, and thus described as CpFBAs.

Based on the 0.6× Chinese Spring genome draft sequence (Mayer et al., 2014), we predicted the chromosomal location of the TaFBA genes (Table 1). In addition, the TaFBA genes were queried on the mapped loci using EST-derived probes in GrainGenes-SQL databases. Therefore, the TaFBA1, 2, 3, 4, 7, 8, 11, 19, and 20 genes were mapped onto different gene loci through a wheat_2014_90K SNP gene chip (Figure S1A), and the TaFBA10, 12, 14, and 15 genes were mapped using a Chinese Spring_Deletion gene chip (Figure S1B).

To confirm the chromosomal location and identify the primer specificity of TaFBA genes, 21 Chinese Spring nulli-tetrasomic lines were used to amplify genomic DNA (The gene -specific primers used are listed in Table S2). The PCR results (Figure 1) confirmed the location of each TaFBA gene copy, as indicated by the failed amplification in the specific nulli-tetrasomic line. These findings validated that the genes were located within the same locus as predicted in Table 1 and Figure S1.

A total of 21 TaFBA genes were localized to 12 chromosomes of 4 homoeologous groups: TaFBA1~3 genes in the short arms and TaFBA10~16 in the long arms of homoeologous group 3; the TaFBA5/6 genes in the short arm, and TaFBA4 in the long arm of homoeologous group 4; the TaFBA7~9 and TaFBA17/18 genes in the short arms of homoeologous group 5; and the TaFBA19~21 genes located in the short arms of homoeologous group 7.

For the investigation of the molecular evolution of the TaFBA gene family, 31 FBA genes from other closely related wheat species were identified: six TmFBA genes from T. monococcum (A^mA^m), seven TuFBA genes from T. urartu (A^uA^u), five AspFBA genes from A. speltoides (SS), six AshFBA genes from A. sharonensis (S^shS^sh), and seven AtaFBA genes from A. tauschii (DD). The gene IDs and sequences are listed in Table S3, and Data Sheets TXT4, TXT 5. Then, an unrooted phylogenetic tree was constructed using MUSCLE in MEGA6 (Figure 2), and pairwise alignments were conducted using the hierarchical clustering method by MeV (Figure S2).

Owing to the transcript annotation of the wheat genome, the analysis and comparison of the structural features of the FBAs in different subfamilies were also conducted (Figure 2). The organization (number, order, and length) of the exons were mostly conserved within different FBA groups, whereas the introns and UTRs showed variable lengths and distribution. Based on the exon numbers and intron lengths, a total of 52 FBA genes from wheat and other closely related species could be divided into six independent subgroups: class Ia ~ class Ie and class II. 21 TaFBA genes belonged to different subgroups: TaFBA7/8/9 genes in class Ia, TaFBA4/5/6 genes in class Ib, TaFBA1/2/3 genes in class Ic, TaFBA10/12/15/17/18 genes in class Id, TaFBA11/13/14/16 genes in class Ie, and TaFBA19/20/21 genes in class II. Each subgroup contained five wheat relative FBA genes, except for six in class Ie. Class II FBA genes existed in both wheat and its closely related species.

Various patterns of exon-intron architecture were identified. For instance, class Ia genes and most cytoplasmic TaFBA genes had one intron, and class Ib and class Ic genes had five introns. On the other hand, class II TaFBA genes contained 41 introns because of the integration of F_bP_aldolase domain and other functional domains. Generally, a higher number of introns could cause alternative splicing and different splice variants, and regulate gene expression at the posttranscriptional level. Whether there are different transcripts of class II genes is thus worth exploring further.

In addition, the annotation of genome sequences suggested that different kinds of transposable elements (TEs) were distributed unevenly in the promoter and intron regions of the TaFBA genomes (Figure S3, Table S4) such as long terminal repeats (LTRs), short interspersed nuclear elements (SINEs), terminal inverted repeats (TIRs), and other TEs. Several TEs were detected in the TaFBA19 and TaFBA20 intron regions, which might influence RNA splicing.

Using HMMER 3.0, a Uniprot KB/Swiss-Prot database search was performed using the query “fructose 1,6-bisphosphate aldolase,” and 211 sequences from different species were downloaded for FBA protein evolutionary analysis. Approximately 23 of these sequences were removed by pairwise comparison and the absence of complete adolase-type TIM domain. Together with plant FBA protein sequences identified by PHMMER from EnsemblPlants database, a total of 245 FBA proteins were obtained, which consisted of 19 archaeal FBAs, 97 bacterial FBAs, and 129 eukaryotic FBAs (listed in Table S5, Data Sheets TXT3, TXT6). Multiple primary sequences of FBA protein were aligned by using Clustal Omega, and a rooted ML phylogenetic tree was constructed by using MEGA6 with default parameters (Figure 3).

Evolutionary analysis revealed that the FBAs in various species were clearly classified into two groups, namely, classes I and II (Figure 3A). The class I FBAs in plants were classified into two subgroups, class I cFBAs, and class I CpFBA, occurring in plastids and the cytoplasm, respectively. The class I cFBAs in plants were more closely related to the FBAs in animals than class I CpFBAs. There were two subgroups (types A and B) in class II FBAs. Most FBAs in fungi belonged to type A, and 10 FBAs in plants (including TaFBA19~21) belonged to type B. Both plant cFBAs and CpFBAs were more distantly related to the FBAs in bacteria and fungi.

To investigate the conserved domain and residues of TaFBA proteins, we aligned the wheat FBA protein sequences and conserved sites by using Pfam and MEME. A total of 8 conserved protein domains and 30 motifs were defined (Figure 3). Most of the class I FBAs possessed a glycolytic protein domain, and the other class I FBAs possessed a DeoC (deoxyribose-phosphate aldolase) domain. Type A of class II FBAs possessed one F_bP_aldolase domain or FBPase_3 (fructose-1,6-bisphosphatase) domain. All of the type B of class II FBAs possessed an FBPase_3 domain, but the class II FBAs of plants possessed the other 4 domains, including NAD_binding_2 (NAD-binding domain of 6-phosphogluconate dehydrogenase), NAD_binding_11 (NAD-binding of NADP-dependent 3-hydroxyisobutyrate dehydrogenase), DUF1537 (putative sugar-binding N-terminal domain), and DUF1357_C domain (putative nucleotide-binding of sugar-metabolizing enzyme), which suggest that class II FBAs are NADP-dependent and probably participate in sugar sensing and signaling in plants. Several identical regions of plant FBAs were fairly conserved such as “EGTLLKPNMVTPG,” “GARFAKWR,” “LSGGQSEEEA,” “GILAADES,” “LVPIYEPE,” “RALQ,” “ILFEET,” and “RAKANS” for class I FBAs (Figure 3B). Among these regions, the “LSGGQSEEEA” (TIM _phosphate_binding superfamily), “LVPIYEPE,” and “ILFEET” matched well with that of animal/fungi/bacterial FBAs, suggesting that FBAs were highly conserved during evolution.

The quaternary structure of three wheat FBA proteins (TaFBA4, TaFBA10, and TaFBA20) without N-terminal transit peptides were predicted by using SWISS-MODEL and I-TASSER (Figure 4). TaFBA4 and TaFBA10 belong to class I subfamily and could form tetramers based on their relatively high similarity to class I rabbit muscle aldolase (Blom and Sygusch, 1997) and human muscle aldolase (Dalby et al., 1999) respectively. TaFBA20 is a class II protein that is predicted to form dimers based on its similarity to class II aldolases of Thermus aquaticus (Izard and Sygusch, 2004).

Similar to other tetrameric FBAs, interfaces A and B were observed in the TaFBA4 and TaFBA10 tetramers (Figures 4D,E). In contrast to interface A, which comprised helix-packing interactions, interface B exhibited loop-loop interactions that were formed by buried hydrogen bonds. In TaFBA4, each interface B had six hydrogen bonds that were involved in the carboxylate “O” atoms of Asp126, interacting with the backbone amides of three consecutive residues (G124-L125-D126) (Sherawat et al., 2008) (Figures 4A–C). Different from TaFBA4, interface B of TaFBA10 was related to Asp124, which interacts with G122-F123-D124.

The catalytic residues in the active region of class I FBAs are highly conserved (Maurady et al., 2002). TaFBA4 contained nine catalytic residues, D32-K105-K144-R146-E184-E186-K226-S268-R298. The catalytic residues in TaFBA10 included D30-K103-K142-R144- E183-E185-K225-S266-R298. Although, class I TaFBAs had the same catalytic residue sequences, different conformations in the active sites would affect the spatial structures (Figures 4G,H). Additionally, class II TaFBAs were the same as the FBAs in Thermus aquaticus, which had active sites including imidazole rings of H81, H178, and H210, and the carboxylate of E132, and the sites also acted as divalent metal ion-binding sites (Figure 4I).

The amino acid sequence alignment of the TIM barrel domain (Figure S4) indicated that the active sites of TaFBA1~TaFBA16 (class I) were homologous to those of rabbit and human FBA isozymes, and the active sites among TaFBA19~21 (class II) corresponded to those of FBA isozymes in Thermus aquaticus. A few substitution mutations and deletion sites were observed in TaFBA17 and TaFBA18, and the catalytic active site sequence of TaFBA17 was D-K-T-P-K-S-R, and that of TaFBA18 was D-K-K-R-K-S-R.

Cis-elements are important molecular switches that are involved in the regulation of gene transcription during plant growth, development, and responses to biotic and abiotic stresses. The 1.5-kb promoter sequences upstream of the 5′-UTR of 19 TaFBA genes (excluded TaFBA19/21) were identified from wheat genome sequences and used in plant cis-acting regulatory elements analysis in the PlantCARE database (Figure 5, Table 2).

Various putative environment stimulus-responsive cis-elements were identified (Table S6). Because TaFBA genes mostly possess multiple phytohormone and abiotic stress-responsive elements, we considered that the expression pattern of TaFBA genes could be regulated by various environmental factors. LREs were highly enriched in 18 TaFBA gene promoters. These findings indicate that the TaFBA genes might be involved in light responses.

The expression profiles of different tissues and developmental stages were downloaded from the WheatExp database and PLEXdb database respectively, and the corresponding probe ID of the TaFBA genes are shown in Table S7. Combined with the results of Genevestigator, we analyzed the spatiotemporal expression profiles of the members of the TaFBA gene family (Figure 6).

The expression of the TaFBA genes showed tissue specificity (Figure 6A). TaFBA4-9 (class I CpFBA genes) showed similar expression patterns; these were upregulated in leaves and downregulated in the roots, particularly TaFBA 4 and 5. On the other hand, most class I cFBA genes were apparently upregulated in the roots. TaFBA12 showed significantly higher expression levels in all tissues, particularly in the roots. TaFBA19/20/21 (class II CpFBA genes) were slightly and constantly expressed in all tissues.

The wheat TaFBA genes were expressed at different developmental stages (Figure 6B). The expression level of TaFBA10/12/18 was higher than that of the other genes, whereas that of TaFBA19/20/21 was lower. The expression profiles of TaFBA4/5/6 and TaFBA7/8/9 varied with different developmental stages. Based on the results of differential expression and phylogenetic analyses, the TaFBA genes were classified into eight gene clusters: TaFBA1/2/3, TaFBA4/5/6, TaFBA7/8/9, TaFBA14/15/17, TaFBA10/12/18, TaFBA11, TaFBA13/16, and TaFBA19/20/21.

Most TaFBA gene promoters were highly enriched with LREs (Figure 5 and Table 2). Thus, the expression patterns of the TaFBA genes under light and dark treatment were investigated by using qRT-PCR (Figure 7). TaFBA1-9 (CpFBA) and TaFBA10/12/18 (cFBA) were downregulated under short-term light treatment (Figure 7A), whereas the other TaFBAs were upregulated. Meanwhile, all TaFBA genes were immediately upregulated within 1 h of transferring from light to dark, especially TaFBA4-9 (Figure 7B). TaFBA14/15/17, TaFBA11, TaFBA13/16, and TaFBA19/20/21 were induced by light, whereas TaFBA4/5/6, TaFBA7/8/9, and TaFBA10/12/18 were induced by short-term exposure to the dark. These results demonstrate that the LRE cis-elements of the TaFBA genes were not directly related to light responses. Therefore, FBA isozymes might play different roles in various biological processes involved in light signal transduction.

The expression profiles of TaFBA gene clusters under different conditions (salt, drought, ABA, heat, and low temperature) were detected using qRT-PCR. The results revealed the complicated expression patterns of the TaFBA gene family in response to various abiotic stresses (Figure 8). The gene specific primers used for qRT-PCR are listed in Table S8.

In shoots (Figure 8A), the class I cytosol TaFBA genes displayed active responses to multiple abiotic stresses. The expression of TaFBA11 and TaFBA13/16 was significantly induced by high temperature. TaFBA14/15/17 and TaFBA10/12/18 were expressed during exposure to drought and ABA. Under salt stress, most TaFBA genes were downregulated except for TaFBA10/12/18. On the other hand, all class I CpFBA genes were apparently sensitive to multiple stresses, as indicated by a repression in expression. On the other hand, class II CpFBA genes TaFBA19/20/21 were induced after short-term exposure to ABA, heat, and low temperature, implying the role of class II TaFBA genes in stress responses.

The expression pattern of TaFBA genes in roots (Figure 8B) was obviously different from shoots. ABA might lead to the similar change as drought to each gene cluster. Both cFBA and CpFBA were induced by drought, ABA, heat and cold, but only cFBA was upregulated by salt. Among cFBA genes, TaFBA10/12/18 were up-regulated in multiple conditions except heat, TaFBA11 and TaFBA13/16 were up-regulated in different stresses except low temperature. In CpFBA genes, TaFBA1/2/3 could be induced by drought and ABA, TaFBA4/5/6 could be induced by drought and heat, and TaFBA19/20/21 could be induced by all stresses, especially by heat and cold.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.