Genome-Wide Identification and Characterization of FBA Gene Family in Polyploid Crop Brassica napus

Fructose-1,6-bisphosphate aldolase (FBA) is a versatile metabolic enzyme involved in multiple important processes of glycolysis, gluconeogenesis, and Calvin cycle. Despite its significance in plant biology, the identity of this gene family in oil crops is lacking. Here, we performed genome-wide identification and characterization of FBAs in an allotetraploid species, oilseed rape Brassica napus. Twenty-two BnaFBA genes were identified and divided into two groups based on integrative analyses of functional domains, phylogenetic relationships, and gene structures. Twelve and ten B. napus FBAs (BnaFBAs) were predicted to be localized in the chloroplast and cytoplasm, respectively. Notably, synteny analysis revealed that Brassica-specific triplication contributed to the expansion of the BnaFBA gene family during the evolution of B. napus. Various cis-acting regulatory elements pertinent to abiotic and biotic stresses, as well as phytohormone responses, were detected. Intriguingly, each of the BnaFBA genes exhibited distinct sequence polymorphisms. Among them, six contained signatures of selection, likely having experienced breeding selection during adaptation and domestication. Importantly, BnaFBAs showed diverse expression patterns at different developmental stages and were preferentially highly expressed in photosynthetic tissues. Our data thus provided the foundation for further elucidating the functional roles of individual BnaFBA and also potential targets for engineering to improve photosynthetic productivity in B. napus.


Introduction
Fructose-1,6-bisphosphate aldolase (FBA, EC4.1.2.13 or aldolase) is an essential metabolism enzyme in the glycolytic pathway [1]. FBA catalyzes the reversible aldol cleavage of fructose-1,6-bisphosphate (FBP) into dihydroxyacetone phosphate (DHAP) and glyceraldehyde-3-phosphate (G3P), two important intermediates for oil biosynthesis [2]. DHAP could be further converted to diacylglycerol (DAG) by multiple enzymatic reactions, and DAG is a key substrate of diacylglycerol acyltransferase (DGAT) for the synthesis of triacylglycerols (TAGs). Meanwhile, G3P could be converted to malonyl CoA (coenzyme A) that is then used to produce fatty acids. Thus, FBA is not only one of the key regulatory enzymes in the glycolysis pathway but also may control the flux of carbohydrates and, therefore, play an important role in the oil yield of oilseeds [3,4].

Identification, Properties, and Genomic Distribution of BnaFBA Genes
To identify each member of the FBA gene family in B. napus, the glycolytic domain (PF00274) and fructose-bisphosphate aldolase class-II domain (PF01116) from the Pfam database (http://pfam.xfam. org/) were employed as queries to search against B. napus var. Darmor-bzh protein dataset. The peptides of putative B. napus FBAs (BnaFBAs) with the best hit of A. thaliana FBAs (AtFBAs) and TaFBAs were further used to predict functional domains by using the Pfam and SMART databases to confirm the presence of the fructose-bisphosphate aldolase domains. As indicated in Table 1, we identified 22 BnaFBA genes in B. napus. Meanwhile, we found 14 BraFBAs and fourteen BolFBAs in B. rapa var. Z1 and B. oleracea var. HDEM, respectively (Table S1).
Based on the functional domains contained, the 22 BnaFBAs could be divided into two classes-class-I and class-II. All the BnaFBAs of class-I harbored one glycolytic domain (PF00274) across the proteins, and BnaFBA9a and BnaFBA9b of class-II had one fructose-bisphosphate aldolase class-II domain (PF01116). Except for the glycolytic domain (PF00274), some members of class-I FBA, such as BnaFBA5b, as the longest FBA of class-I with 870 aa, also harbored PRK (phosphoribulokinase kinase) (PF00485) and UPRTase (uracil phosphoribosyltransferase) (PF14681) domains in the N terminal. Likewise, in B. oleracea, BolFBA5 contained three functional domains of PF00274, PF00485, and PF14681, and BolFBA6b had two tandem copies of glycolytic domains (PF00274). Interestingly, BraFBA6a contained both glycolytic domain (PF00274) and RNA polymerase II-binding domain (PF04818), suggesting that it might be involved in RNA processing in the nucleus (Table S1).
To attain a general view of the distribution of BnaFBA genes on the genome of B. napus, the 22 BnaFBAs were mapped on the corresponding chromosomes according to their physical positions. Of the 22 BnaFBA genes, 19 were evenly distributed on the 15 B. napus chromosomes, while the three other BnaFBAs were assigned to the random chromosomes (two on the An chromosomes and one on the Cn chromosome). The number of BnaFBAs per chromosome ranged from one to two. Among them, eleven BnaFBA genes were distributed over the seven A subgenome chromosomes of B. napus, including A01, A02, A04, A06, A07, A08, and A09, as well as Ann_random. Equally, the other 11 BnaFBAs were distributed over eight C subgenome chromosomes of B. napus, including C01, C02, C03, C04, C05, C06, C07, and C08, as well as Cnn_random. Most of the BnaFBAs were not localized in the terminal regions of the chromosomes, where the gene density was relatively high in B. napus (Figure 1).

Phylogenetic and Structure Analysis of BnaFBAs
To explore the molecular evolution of the FBA gene family in B. napus, a total of 79 FBA genes from B. napus, B. rapa, B. oleracea, A. thaliana, and wheat, were used to construct an unrooted phylogenetic tree. According to the phylogenetic relationships of these FBA genes, they could be divided into two independent classes, consistent with the classification by functional domains they contained ( Figure 2A). Furthermore, the FBA genes of class-I could be further classified into four subclasses, namely class-Ia, class-Ib, class-Ic, and class-Id. Therefore, in B. napus, there were seven BnaFBA genes in the class-Ia group, two in class-Ib, nine in class-Ic, two in class-Id, and two in class-II. Consistent with its polyploid origin, except for the FBA4 gene, the genome of B. napus maintained homologs of each FBA gene derived from the diploid parents, B. rapa and B. oleracea (Figure 2A).
Based on the gene information of the genome available in the GENOSCOPE database, we performed gene structure analysis by comparing the coding sequence (CDS) of BnaFBAs, BraFBAs, and BolFBAs. In B. napus, A subgenome homologs and C subgenome homologs of FBA genes that were in the adjacent branches of the phylogenetic tree exhibited the same gene structure ( Figure 2B). Within each class of FBA genes, the features of exons, such as order, length, and number, were largely conserved except for FBA6 genes ( Figure 2B). Besides, the organization of the introns of BnaFBA genes was highly variable. The length of introns varied extensively in different members of the BnaFBA gene family, ranging from 30 to 4152 bp, with the number ranging from 2 to 41. Compared to the class-I BnaFBA genes, the class-II genes were much longer and had much more exons and introns ( Figure 2B). To further explore the higher-order structure of the BnaFBA proteins, the three-dimensional (3D) structural models for BnaFBA1a, BnaFBA8a, and BnaFBA9a were generated using SWISS-MODEL. Based on the experimental structure of class-I rabbit muscle aldolase, the SWISS-MODEL analysis results revealed that BnaFBA1a and BnaFBA8a of the class-I group could form tetramers structures, and interfaces A and B were observed in the BnaFBA1a and BnaFBA8a tetramers. Different from BnaFBA1a and BnaFBA8a, BnaFBA9a belonged to class-II group and form dimers based on its relatively high similarity to class-II aldolases of Thermus aquaticus ( Figure 3A-C). Furthermore, the catalytic residues of D74-K147-K186-R188-E226-E228-K269-S301-R331 and D30-K103-K142-R144-E183-E185-K225-S266-R298 were observed in BnaFBA1a and BnaFBA8a, respectively. Additionally, similar to the class-II FBAs of Thermus aquaticus, BnaFBA9a contained active sites of H1181-E1232-H1277-H1309 that also serve as the divalent metal cation binding sites ( Figure 3E,F). Multiple sequence alignment results indicated that class-I BnaFBAs, BraFBAs, BolFBAs, and AtFBAs had high conserved catalytic residues ( Figure S1). The active sites among class-II BnaFBAs, BraFBAs, and BolFBAs were the same as FBA isozymes in the Thermus aquaticus ( Figure S2).

Synteny and Gene Duplication of BnaFBA Genes
A. thaliana is the most prominent model system for plant molecular biology and genetics research, whose structural genes have been identified and functionally characterized. Thus, we traced the orthologous gene pairs between A. thaliana and Brassica species to investigate the evolutionary history by syntenic gene analysis. A total of 29 orthologous gene pairs were identified between A. thaliana and B. napus, 20 between A. thaliana and B. rapa, and 19 between A. thaliana and B. oleracea. In addition, we also obtained 29, 15, and 15 paralogous gene pairs within B. napus, B. rapa, and B. oleracea, respectively ( Figure 4A). The previous study revealed that crucifer (Brassicaceae) lineage experienced two whole-genome duplications (named α and β) and one triplication event (γ) shared by most dicots [37]. Moreover, Brassica species experienced an extra whole-genome triplication (WGT) event compared with A. thaliana [38]. As WGT of the Brassica ancestor, FBA genes in the A. thaliana genome might have triplicated orthologous copies in B. rapa and B. oleracea. Consequently, some FBA genes (i.e., FBA2 and FBA8) existed triple the number of those in A. thaliana, while the other genes (i.e., FBA1, FBA3, FBA5, and FBA6) had double or equal the number ( Figure 4B). The FBA genes of B. napus were inherited from its diploid ancestors; thus, most of the BnaFBA genes were double the number of those in B. rapa and B. oleracea (i.e., FBA2, FBA3, FBA5, FBA8, and FBA9). However, both FBA1 and FBA6 genes lost one copy in the C subgenome of B. napus, while FBA4 lost all the copies in B. napus compared to its two ancestors ( Figure 4B). Gene duplication analysis with syntenic and phylogenomic approaches using tool DupGen_finder in B. napus showed that all BnaFBA genes had corresponding duplicate genes. In B. napus, a total of 42 gene pairs were derived from whole-genome duplication (WGD), with one gene pair of BnaFBA2e-BnaFBA2f being derived from transposed duplication (Table S2).

Cis-Acting Elements in the Putative Promoter Regions of BnaFBA Genes
As important molecular switches, cis-acting elements in the promoter region could provide useful information to understand the function and regulation of the genes during plant development and responses to various stresses. The 1.5 kb genomic DNA sequences identified from upstream of the BnaFBA genes were extracted and deployed in cis-acting regulatory elements analysis with PlantCARE. Various cis-acting regulatory elements existed within the promoter regions of BnaFBA genes ( Table 2 and Table S3). For example, BnaFBA genes contained multiple phytohormone responsive elements, such as ABRE (abscisic acid-responsive element), AuxRE (auxin-responsive element), ERE (ethylene-responsive element), GARE (gibberellin-responsive element), MeJARE (MeJA-responsive element), and SARE (salicylic acid-responsive element). This suggested that the expression of BnaFBAs might be induced by different phytohormones. Besides, the cis-acting regulatory elements involved in stress-responsive elements, such as ARE (anoxic-responsive element), DRE (damage-responsive element), DIRE (drought-responsive element), DSRE (drought-and stress-responsive element), HSRE (heat stress-responsive element), LTRE (low-temperature-responsive element), and WRE (wound-responsive element), were also found within the promoters of BnaFBA genes, suggesting that expression levels of BnaFBAs might be also regulated by various environmental factors like drought, heat, and low-temperature. Globally, three phytohormone-related elements (i.e., ABRE, ERE, and MeJARE) and two stress-responsive elements (ARE and HSRE) were detected with high frequency in the promoter regions of BnaFBA genes. Notably, each BnaFBA had multi-copy LREs (light-responsive elements) ranging from 8 to 26, implying that BnaFBAs might play roles in light responses.

Natural Variations of BnaFBA Family Genes in B. Napus
Critical sequence polymorphism across the gene and its flanking regions may reflect the evolutionary process of species adapting to different environments. The public resequencing datasets of 991 B. napus germplasm accessions covering three main ecotype groups, namely winter, semi-winter, and spring types, were collected for variation analysis of the BnaFBA family genes [39]. The polymorphism sites of CDS in BnaFBA family genes ranged from two (BnaFBA2c) to 130 (BnaFBA6) ( Table 3). The π and θ w of nucleotide diversity parameters extended from 0.00057 (BnaFBA3b) to 0.04158 (BnaFBA6) and from 0.00054 (BnaFBA3b) to 0.02255 (BnaFBA6), respectively. Some members of BnaFBA gene family were conserved, such as BnaFBA3b (θ w = 0.00054), BnaFBA2c (θ w = 0.00092), and BnaFBA2d (θ w = 0.00098), while others had high polymorphism, such as BnaFBA9b , and BnaFBA6 (θ w = 0.0225). Besides, the BnaFBA6 variation ratio reached a peak of 12.60 among the 22 BnaFBA genes, whereas BnaFBA2c had the lowest variation ratio of 0.17, with only two polymorphic sites (Table 3). Generally, due to the longer length of genes, BnaFBA9a and BnaFBA9b of class-II had much more variations than FBA genes of class-I; however, the variation ratio and nucleotide diversity of the coding regions of BnaFBA9a and BnaFBA9b showed no difference with the BnaFBA genes of class-I except for BnaFBA6 (Table 3). In addition to CDS regions, a total of 1029 and 814 variations were also identified in the upstream/downstream 1.5 kb regions and intronic regions of BnaFBA genes, respectively (Table S4). Notably, BnaFBA8a, BnaFBA8b, and BnaFBA6 each harbored one stop-gain mutation that led to premature stop codons, which indicated that these genes exhibited loss-of-function ( Figure 5, Table S5). Only a few of the Indels had been detected in BnaFBAs. For example, BnaFBA1c had three non-frameshift Indels, and BnaFBA3b contained only one non-frameshift Indel, with BnaFBA6 harboring two frameshift deletion/insertion variations ( Figure 5). To study the population selection pressure, we conducted neutral testing using Tajima's D. Tajima's D values of all the BnaFBAs were positive, with significant Tajima's D values (p < 0.01) being observed in BnaFBA1a, BnaFBA1e, BnaFBA2d, and BnaFBA2e (Table 3). Particularly, the Tajima's D values of BnaFBA1b and BnaFBA8a reached extremely significant levels (p < 0.001 and p < 0.0001). Based on the variation of BnaFBA genes, we performed principal component analysis (PCA), in which no significant difference in the polymorphism of BnaFBA genes was seen between different ecotype groups of B. napus ( Figure S3).

Expression Patterns of BnaFBA Genes
The members of the FBA gene family are playing important roles in diverse significant physiological and biochemical processes in plants [21,25]. Besides, FBAs also participate in response to various environmental stimuli [30][31][32]. To identify the function of BnaFBA genes under various conditions, the expression levels of BnaFBAs were evaluated during growth and development, as well as in response to biotic stresses and phytohormones in B. napus.
To explore the expression patterns of the BnaFBAs during growth in B. napus, we analyzed their expression levels in twelve various tissues (leaf, root, stem, cauline leaf, pistil, stamen, petal, flower bud, axillary bud, silique wall, embryo, and seed) at different developmental stages. As a result, at the trefoil stage and the flowering stage, BnaFBA genes showed opposite expression patterns between leaves and roots. For example, BnaFBA1a/b/c, BnaFBA2a/b/c/d/f, and BnaFBA5a/b displayed relatively high expression level in leaves than that in roots. On the contrary, BnaFBA3a/b and BnaFBA8a/b/c/d/e/f exhibited relatively high expression levels in roots than that in leaves. At the flowering stage, the expression levels of BnaFBA1a/b/c and BnaFBA2a/b/c/d/f were higher in pistil than that in stamen and petal tissues. At the pod stage, BnaFBA1a/b/c and BnaFBA2a/b/c/d/f were highly expressed in the silique wall, with their expression levels being increased during the development of pod. BnaFBA3a, BnaFBA3b, BnaFBA8c, and BnaFBA8d were highly and constantly expressed in all tissues, particularly showing the highest expression levels in the embryo at 25 days after pollination (DAP). In addition, the expression of BnaFBA2a/b/c/d/e/f increased during the development of seeds and reached high levels at 4 weeks after pollination (WAP), then decreased in the following mature stages of seed development ( Figure 6A).
BnaFBA9a and BnaFBA9b of class-II showed constitutive expression at relatively low levels in all the tissues examined. BnaFBA6 of class-I exhibited unexpressed in most tissues. Taken together, these results suggested that BnaFBA genes displayed diverse spatiotemporal expression patterns during the growth and development of different tissues in B. napus.  Figure 6B). Under drought stress, the expression levels of BnaFBA1a/c, BnaFBA2a/b/c/d/e/f, and BnaFBA5a/b decreased in the leaves. Similarly, the expression levels of BnaFBA1b, BnaFBA2b/d/e/f also decreased in roots compared to the control ( Figure 6B). We further examined the expression of BnaFBAs under exogenous SLs treatments in B. napus. Notably, the members of the BnaFBA2 group showed diverse expression patterns after SLs treatments in B. napus. At 12 h after SLs treatments, the expression levels of BnaFBA2a/b/c/f were significantly down-regulated, while the expression of BnaFBA2d was slightly decreased relative to the control. At 1 day after SLs treatment, BnaFBA2a/b/c/d expression levels were significantly up-regulated compared to the control, whereas BnaFBA2e/f showed no difference. At 4 days after SLs treatment, all members of the BnaFBA2 group exhibited up-regulated expression compared to the control. Then, at 7 days after SLs treatments, the expression levels of BnaFBA2a/b/c/d turned to be down-regulated, while BnaFBA2e/f expression levels were still up-regulated compared to the control ( Figure 6B).

qRT-PCR Analysis of Selected BnaFBA Genes under Abiotic Stresses
To further validate the functional roles of BnaFBAs under abiotic stresses, four BnaFBA genes from different clusters were selected for the examination of their expression levels under three stress conditions using quantitative real-time PCR (qRT-PCR) in B. napus. These genes included BnaFBA2a, BnaFBA5b, and BnaFBA8a of class-I, and BnaFBA9a of class-II. The qRT-PCR analysis was carried out using rapeseed plants exposed to salt, heat, and drought stresses. At 3 days after NaCl treatment, the expression level of BnaFBA8a was significantly up-regulated by approximately 11-fold compared with the control, while BnaFBA2a, BnaFBA5b, and BnaFBA9a were significantly down-regulated, particularly, the expression levels of BnaFBA5b and BnaFBA9a were largely reduced (Figure 7). In addition, at three days after heat and drought stress treatments, the expression level of BnaFBA5b was significantly reduced nearly by half, whereas the expression of BnaFBA8a was slightly increased compared to the control (Figure 7).

Co-Regulatory Networks of BnaFBA Genes
Based on the publically available RNA-seq datasets collected from different tissues, biotic and abiotic stresses, and phytohormone treatments in B. napus, we calculated the Pearson correlation coefficients (PCCs) of the expression levels of BnaFBA genes and constructed co-regulatory networks. Positive correlations were observed between members of class-Ia and class-Ib, such as BnaFBA3a/b, BnaFBA8a/b/c/d/e/f, and BnaFBA6. Likewise, the BnaFBA genes of class-Ic and class-Id also showed positive correlations between each other. Particularly, members of class-Id, including BnaFBA1a/b/c and BnaFBA2a/b/c/d, exhibited strong positive correlations. BnaFBA2b and BnaFBA2f showed significant negative correlation with BnaFBA3a/b and BnaFBA8a/b/e/f. However, they displayed a significant positive correlation with other BnaFBA genes of class-Ic and class-Id ( Figure 8A). All the significant PCCs (p-value ≤ 0.05 and |PCC| > 0.5) of BnaFBAs were extracted and used to construct co-regulatory networks delineated by Cytoscape

Discussion
The Calvin cycle is the initial pathway of photosynthetic carbon fixation, which plays a requisite role in the growth and development of plants [40]. Extensive efforts for seeking a breakthrough in the regulation of this cycle have been made to substantially enhance photosynthetic CO 2 capacity and plant productivity. The carboxylation capacity of ribulose bisphosphate carboxylase oxygenase (Rubisco) and the regenerative capacity of ribulose diphosphate (RuBP) are uncovered to be essential for maintaining high photosynthetic CO 2 fixation capacity. Previous studies revealed that three non-regulated enzymes, including fructose−1,6-bisphosphate aldolase (FBA or aldolase), sedoheptulose 1,7-bisphosphatase (SBPase), and transketolase (TK), had significantly higher control coefficient on photosynthesis than the other Calvin cycle enzymes, which indicated that they could limit photosynthetic rate and exert significant control over photosynthetic carbon flux other than Rubisco [41]. Particularly, FBA could catalyze the reversible conversion of DHAP and FBP and also catalyze DHAP and erythrose 4-phosphate (E4P) to sedoheptulose 1,7-bisphosphate (SBP). Thus, FBA is not only one of the key regulatory enzymes in the glycolysis pathway but also may lie in a vital strategic position to determine the carbon partitioning in the Calvin cycle, which made FBA probably an important candidate target of engineering to boost photosynthetic carbon CO 2 fixation capacity [42].
B. napus is the second-largest source of vegetable oil crops and is cultivated around the world [33,34]. B. napus seeds contain oil, carbohydrates, and proteins as major storage reserves. In most seeds, glycolysis in plastids supplies carbon for fatty acid (FA) synthesis [43]. FBA, as a key enzyme in the glycolytic metabolism, provides precursors for amino acid and fatty acid synthesis [26]. To our best knowledge, there is not any report in the literature describing FBA family genes and their function in B. napus so far. In this study, we performed a comprehensive identification and characterization of the FBA gene family in B. napus. Compared to the number of FBA genes identified in other plant species, such as eight in A. thaliana [17], seven in rice [22], eight in tomato [19], and 21 in wheat [11], B. napus had more FBA genes, including 22 BnaFBA genes distributed in fifteen B. napus chromosomes and two random chromosomes. The BnaFBAs can be classified into two classes based on the functional domains contained, class-I and class-II. Enzyme kinetics analysis of aldolase 1 (class-I FBAs) and aldolase 2 (class-II FBAs) in Escherichia coli revealed that aldolase 1 and 2 hydrolyze fructose 1,6-bisphosphate by the aldol cleavage reaction [12]. Class-I BnaFBA proteins could form tetramer structures with high conserved catalytic residues of D-K-K-R-E-E-K-S-R that were homologous to those of rabbit FBA isozymes [44]. Class-II BnaFBA could form dimers with the active sites of H-E-H-H that corresponded to those of FBA isozymes in the Thermus aquaticus [45,46]. Eleven BnaFBAs of class-I and BnaFBA9a of class-II were predicted to be localized in chloroplast, while nine class-I BnaFBAs and BnaFBA9b of class-II were in the cytoplasm. This is consistent with the subcellular localization of their homologs in A. thaliana [17] and wheat [11]. Although AT1G18270 in Arabidopsis is the homolog of BnaFBA9a and BnaFBA9b, and it is annotated as ketose-bisphosphate aldolase class-II family protein; however, it contains fructose-bisphosphate aldolase class-II (pfam01116) domain and zinc-binding site and is further assigned as fructose-bisphosphate aldolase (NCBI Reference Sequence: NP_001117303.1). More accurately, AT1G18270 is fructose-bisphosphate aldolase class-II family protein. This inaccuracy might result from the homologous protein annotation in bacteria, such as Variovorax paradoxus B4. Now, the ketose-bisphosphate aldolase, class-II protein, was further assigned as fructose-1,6-bisphosphate aldolase, class-II in Variovorax paradoxus B4 (GenBank: AGU50177.1). Therefore, BnaFBA9a and BnaFBA9b were classified as class-II FBA genes in Brassica napus.
In contrast to the model plant A. thaliana, except for the paleo-polyploidization of alpha (α)-beta (β)-gamma (γ) WGD events, the Brassica species, such as B. rapa and B. oleracea, experienced an extra whole-genome triplication (WGT) event at approximately 15.9 Mya [47][48][49]. B. napus is a relatively new species of Brassica genus with a short history of post-Neolithic speciation (~7500 years) and domestication (~700 years) and a recent allotetraploid formed from hybridization between B. rapa and B. oleracea [35]. Therefore, B. napus has experienced five genome duplication events (3 × 2 × 2 × 3 × 2) at times during the evolution process. Gene duplication that resulted from whole-genome duplication (WGD), tandem duplication (TD), proximal duplication, transposed duplication, or dispersed duplication is one of the primary driving forces to the evolution of morphological and physiological diversity in plants [50]. Our results showed that Brassica species had an extended FBA gene family, containing extra copies of FBA1, FBA2, FBA6, and FBA8. Most gene duplications of the FBA gene family resulted from WGD (Table S2). Considering the collinear correlations and subgenome, transposed duplication and triplication events in the FBA family had contributed to this expansion prior to the speciation of Brassica species. WGD is often followed by intensive gene loss to adapt to continuously changing environments. Compared to A. thaliana and its two ancestors, B. napus species lost some copies of FBA genes (i.e., FBA1, 3, 4, 5, and 6). The remaining duplicated or triplicated FBA genes identified might be conducive to the adaptation of B. napus to various adverse environments during speciation and domestication.
Genome resequencing provides an effective way to identify a large number of variations, which lay a foundation for further identification and functional validation of candidate genes contributed to important traits in crop plants, such as rice, tomato, and soybean [51][52][53]. Critical sequence polymorphisms across the gene and its flanking regions may reflect the evolutionary trends and breeding selection effects on the genes. Based on the resequencing of a worldwide collection of 991 B. napus germplasm accessions that released recently, we explored the pattern of genetic polymorphisms of the BnaFBA family genes. The results showed that the 22 BnaFBAs diversified in sequence polymorphism, with the polymorphism sites ranging from two to 130 (Table 3). The BnaFBA genes showed different levels of polymorphism. For example, BnaFBA6, BnaFBA1c, and BnaFBA8d with the highest levels of genetic variation polymorphism might have experienced weak selection pressure during the evolution process, whereas BnaFBA2c, BnaFBA2d, and BnaFBA3b were highly conserved after strong selection. Particularly, BnaFBA6 harbored the highest levels of genetic variation, including two missense polymorphism sites with a frequency of 12.12% and 5.72%, respectively (Table S5). Moreover, a one-stop gained site with a frequency of 1.02% was also found in the CDS regions of the BnaFBA6 gene. These results suggested that the BnaFBA6 gene might be under the pseudogenization process. Furthermore, BnaFBA genes might have experienced balancing selection or population shrinkage. For example, six FBA family genes, including BnaFBA8b, BnaFBA1b, BnaFBA1c, BnaFBA2e, BnaFBA1a, and BnaFBA2d, harbored some selected signals and probably underwent breeding selection in the process of natural selection and domestication in B. napus.
Fructose-1,6-bisphosphate aldolase (FBA) is a non-regulated enzyme in the Calvin cycle, whose activity is not regulated by effectors or posttranslational modification but by expressional regulation or protein degradation [42]. Recent studies using transgenic plants with reduced enzyme levels have revealed that FBAs play important roles in regulating carbon flux through the Calvin cycle. Elevated plastid aldolase activity accelerated RuBP regeneration and resulted in increased photosynthetic capacity, growth rate, and biomass yield in tobacco and cyanobacterium [27,54]. Increased activity of FBA in Anabaena sp. strain PCC 7120 stimulated photosynthetic yield [55]. Generally, the expression levels of BnaFBA family genes were higher in the overground tissues than that in the underground. In B. napus, some members of the FBA family, such as BnaFBA5a, BnaFBA5b, BnaFBA2b, BnaFBA2c, and BnaFBA2d, showed higher expression levels in mature leaves and the top of stems than that in young leaves and the base of stems. Notably, BnaFBA genes were preferentially highly expressed in the photosynthetic tissues and stages, particularly leaves and silique wall in B. napus. For example, during the development of silique, BnaFBA5a and BnaFBA5b exhibited the highest expression levels at 10 days after DAP, while BnaFBA1a, BnaFBA2a, and BnaFBA2b had highest expression levels at 15 days after DAP in the silique wall. Besides, BnaFBA1a, BnaFBA2a, BnaFBA2b, BnaFBA2c, BnaFBA2d, BnaFBA2e, and BnaFBA2f showed the highest expression levels at four weeks after DAP in the seed tissues when seed fill begins with rapid embryo growth, oil biosynthesis, and protein accumulation [56]. In addition, BnaFBA3a, BnaFBA3b, BnaFBA8c, and BnaFBA8d also showed the highest expression levels at 25 days in the embryo tissues. Thus, the members of BnaFBAs that have a strong correlation with photosynthetic events are potential targets for engineering to improve photosynthetic capacity in the future, which would be further investigated in our next study.
Previous studies revealed FBA family genes also involved in plant defense and response to various biotic and abiotic stresses, such as cold and heat stress [17], salt stress [21], drought stress [21], water-deficit stress [30], stress with Rhizoctonia solani Kuhn [31], and high light acclimation stress [32]. Here, we found that the members of the BnaFBA family showed diverse expression patterns in response to Sclerotinia sclerotiorum infection and drought stress in B. napus. For example, BnaFBA1a/b/c, BnaFBA2a/b/c/d/e/f, BnaFBA5a/b, and BnaFBA9b were down-regulated, while BnaFBA3a/b and BnaFBA8a/b/c/d/e/f were up-regulated compared to the control in the leaves at 24 h after Sclerotinia sclerotiorum infection. Under drought stress, the expression levels of 10 BnaFBA genes (BnaFBA1a/c, BnaFBA2a/b/c/d/e/f, and BnaFBA5a/b) decreased in the leaves and five BnaFBA genes (BnaFBA1b, BnaFBA2b/d/e/f ) lowered in roots. Besides, we applied salt, heat, and drought stresses on B. napus seedlings to explore the expression changes of BnaFBAs. The expression of the BnaFBAs selected was induced by the three stresses. The expression levels of BnaFBA2a, BnaFBA5b, and BnaFBA9a were down-regulated, while BnaFBA8a was up-regulated compared to the control after treatments. In addition, FBA family genes were reported to play roles in response to phytohormones, such as ABA [17] and GA [29]. Various cis-acting regulatory elements related to phytohormone responsive elements were observed in the promoter regions of BnaFBA genes in B. napus, including ABRE, AuxRE, ERE, GARE, MeJARE, and SARE, indicating that the expression levels of BnaFBA family genes might be affected by ABA, IAA (indolylacetic acid), ethylene, GA, MeJA, and SA. Furthermore, our results showed that BnaFBAs, such as BnaFBA2a/b/c/d/e/f, could be also induced by strigolactones (SLs), a new class of plant hormones playing functional roles in the development of root.
In wheat, some FBA genes showed close correlations in expression patterns and could be classified into different clusters, such as TaFBA1/2/3, TaFBA14/15/17, and TaFBA19/20/21 [11]. Similarly, based on the co-regulatory networks of BnaFBAs, the 22 BnaFBA family genes could be divided into three gene clusters. Strong positive correlations between the expression levels were observed among the members within each gene cluster of BnaFBAs. Meanwhile, only a few significant negative correlations appeared between minority members of cluster 1 and cluster 2. Cluster 3, consisting of BnaFBA9a/b of class-II, was independent among the three clusters. The results suggested that the BnaFBA genes might have functional redundancy within gene clusters and functional diversification among gene clusters in B. napus.
In summary, we performed a genome-wide identification of FBA family genes in B. napus, as well as B. rapa and B. oleracea, and further investigated their gene structures and phylogenetic relationships. The cis-acting regulatory elements in the promoter regions, natural variations, and expression patterns of BnaFBAs in different tissues or under treatments were analyzed. Our findings provided useful information regarding FBA family genes in B. napus. Remarkably, we found that some members of the FBA family that underwent breeding selection and were highly expressed in leaves and silique wall probably had positive roles in the promotion of photosynthetic capacity in B. napus. These BnaFBA genes might be utilized in the development and selection of high-yield B. napus cultivars.

Identification and Property Analysis of BnaFBA Genes
The genome sequence, annotation, and protein datasets of B. napus var. Darmor-bzh, B. rapa var. Z1, and B. oleracea var. HDEM were obtained from the GENOSCOPE database (http://www. genoscope.cns.fr/brassicanapus/ and http://www.genoscope.cns.fr/externe/plants/). The glycolytic domain (PF00274) and fructose-bisphosphate aldolase class-II domain (PF01116) from the Pfam database (http://pfam.xfam.org/) were applied as queries to search against local B. napus protein sequence dataset using HMMER (version 3.2.1, HHMI, Chevy Chase MD, USA) with E-value setting at 1e-5 [57]. Then, for further confirmation of BnaFBA proteins, the sequences of predicted BnaFBA proteins were searched against all the annotated proteins of A. thaliana and wheat (Triticum aestivum L.) using BLASTP (version 2.2.26, Bethesda, MD, USA) with E-value < 1 × 10 -5 . The putative BnaFBAs with best hits of the A. thaliana and wheat FBA proteins remained and were further deployed to determine the fructose-bisphosphate aldolase domains by using the Pfam and SMART databases (http://smart.embl-heidelberg.de/). The molecular weight (MW), isoelectric point (pI), and grand average of hydropathy (GRAVY) of each BnaFBA were calculated using the ProtParam tool (http://web.expasy.org/protparam/).

Cis-Acting Regulatory Elements and Subcellular Localization Analysis
The 1.5 kb promoter sequence upstream from the transcription start site of each BnaFBA gene was extracted from the B. napus genome sequence and used to predict cis-acting regulatory elements by PlantCARE [58]. The subcellular localization of each BnaFBA was predicted by Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant/) [59].

Structure and Chromosomal Localization Analysis
The gene structures of the BnaFBA genes were inferred by aligning their coding sequences to the B. napus genomic sequence. Then, a schematic map of the exon-intron structure of each BnaFBA gene was drawn by Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn/) [60]. Three-dimensional (3D) protein models of the BnaFBAs were generated using SWISS-MODEL [61,62]. The physical chromosomal locations of BnaFBAs were obtained from the B. napus genome annotation information. The graphical representation of the BnaFBA genes on chromosomes was plotted using R software with RIdeogram package (https://github.com/TickingClock1992/RIdeogram) [63].

Multiple Alignments and Phylogenetic Analysis
Multiple alignments of FBA protein sequences from B. napus, B. rapa, B. oleracea, A. thaliana, and wheat were performed using MUSCLE (version 3.8, Hinxton, Cambridge, UK) with default parameters [64]. The phylogenetic tree was generated using MEGA7.1 with the neighbor-joining method, and the robustness of each node in the tree was determined using 1000 bootstrap replicates [65].

Expression Analysis of BnaFBA Genes
The publically available RNA-seq dataset of 12 different tissues (sepal, pistil, stamen, ovule, pericarp, blossomy pistil, wilting pistil, root, flower, leaf, silique wall, and stem) collected from different stages of the growth (BioProject ID: PRJNA394926) [36], RNA-seq dataset of leaf and root tissues under drought stress (BioProject ID: PRJNA256233) [74], RNA-seq dataset of seeds across four phases of the development (BioProject ID: PRJNA311067) [56], RNA-seq dataset of stems and leaves after Sclerotinia sclerotiorum infection (BioProject ID: PRJNA321917) [75], and time-series RNA-seq dataset of roots under a synthetic analog of strigolactones (rac-GR24) treatments (BioProject ID: PRJNA484313) [76] were downloaded from the NCBI SRA database, and further used as main sources to perform gene expression profiling of BnaFBA genes in B. napus. The transcriptome reads were mapped to B. napus var. Darmor-bzh reference genome using HISAT2 (version 2.1.0, Baltimore, MD, USA) with the default settings [77]. The read counts per gene were generated by featureCounts [78]. Fragments per kilobase of exon per million fragments mapped (FPKM) was used for the quantification of gene expression. The clustered heatmaps were visualized with expression levels (log2) of BnaFBA genes by R software using the pheatmap function package (https://cran.r-project.org/web/packages/pheatmap/).

Plant Materials and Treatments
Rapeseed seeds were germinated on a filter paper saturated with distilled water in darkness at 22 • C for two days. Then, the seedling plants were transferred to a 4 L hydroponic system containing continuously aerated 1/2 Murashige and Skoog (MS) liquid solution (pH 5.8, without agar and sugar) and grown in an incubator under a photosynthetic flux of 160 µmol photons m −2 s −1 and a humidity of about 50% (16 h light at 25 • C/8 h darkness at 22 • C). The 1/2 MS liquid solution was changed once every two days. After three weeks, the seedlings were transferred to a new 1/2 MS liquid solution (pH 5.8, without agar and sugar) for different stress treatments. For salt, heat, and drought stress treatments, seedlings were exposed to 1/2 MS solution (pH 5.8, without agar and sugar) containing 250-mM NaCl, 40 • C conditions, and 20% (w/v) polyethylene glycol (PEG), respectively. Seedlings exposed to 1/2 MS solution at 22 • C were used as controls. Leaf samples were collected 3 days after each treatment. All collected samples were immediately frozen in liquid nitrogen and stored at −70 • C for further analysis.

RNA Isolation and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) Analysis
Total RNAs were extracted from each sample using an RNA extraction kit (Takara, Dalian, China) following the manufacturer's procedure. Two micrograms of total RNA was used to synthesize the first-strand cDNA using the Prime Script RT reagent Kit (Takara, Dalian, China) according to the manufacturer's protocol. Real-time quantitative PCR was performed using 2 mL of cDNA in a 20 mL reaction volume using an SYBR Green PCR kit (GeneCopoeia Inc., Rockville, MD, USA) with ViiATM 7 Dx platform (ABI, Los Angeles, CA, USA). The qRT-PCR reaction condition was as follows: 95 • C for 5 min, 40 cycles at 95 • C for 30 s, 55 • C for 30 s, and 72 • C for 30 s. Gene-specific primers were designed and listed in Supplementary Materials (Table S6). The relative expression levels of these genes were analyzed by the 2 − Ct method. The BnaTMA7 gene (BnaC05g11560D), which exhibited stable expression in different/same tissues under various experimental conditions, was used as an internal control [79]. All qRT-PCR reactions were assayed in triplicates.

Pearson Correlation Analysis
On the basis of the RNA-seq results, the Pearson correlation coefficients (PCCs) and p-value of the expression levels of BnaFBA gene pairs were calculated using R software with cor and cor.test function packages, respectively. The correlation heatmap was generated by R software with the corrplot function package. The gene co-regulatory networks were constructed by Cytoscape (version 3.1, Seattle, WA, USA) based on the PCCs of BnaFBA gene pairs with a p-value ≤ 0.05 [80].

Conflicts of Interest:
The authors declare no conflict of interest.