Genome-Wide Identification and Expression Profiling of Cytokinin Oxidase/Dehydrogenase (CKX) Genes Reveal Likely Roles in Pod Development and Stress Responses in Oilseed Rape (Brassica napus L.)

Cytokinin oxidase/dehydrogenases (CKXs) play a critical role in the irreversible degradation of cytokinins, thereby regulating plant growth and development. Brassica napus is one of the most widely cultivated oilseed crops worldwide. With the completion of whole-genome sequencing of B. napus, genome-wide identification and expression analysis of the BnCKX gene family has become technically feasible. In this study, we identified 23 BnCKX genes and analyzed their phylogenetic relationships, gene structures, conserved motifs, protein subcellular localizations, and other properties. We also analyzed the expression of the 23 BnCKX genes in the B. napus cultivar Zhong Shuang 11 (‘ZS11’) by quantitative reverse-transcription polymerase chain reaction (qRT-PCR), revealing their diverse expression patterns. We selected four BnCKX genes based on the results of RNA-sequencing and qRT-PCR and compared their expression in cultivated varieties with extremely long versus short siliques. The expression levels of BnCKX5-1, 5-2, 6-1, and 7-1 significantly differed between the two lines and changed during pod development, suggesting they might play roles in determining silique length and in pod development. Finally, we investigated the effects of treatment with the synthetic cytokinin 6-benzylaminopurine (6-BA) and the auxin indole-3-acetic acid (IAA) on the expression of the four selected BnCKX genes. Our results suggest that regulating BnCKX expression is a promising way to enhance the harvest index and stress resistance in plants.


Introduction
Cytokinins promote cell division, regulate root and shoot differentiation synergistically with auxin [1][2][3][4][5], release apical dominance, and delay leaf senescence [6][7][8]. Cytokinins also function in responses to various biotic and abiotic stresses [9,10]. Most naturally occurring cytokinins are adenine derivatives whose N 6 -side chain is substituted by an isoprenoid or aromatic side-chain. Cytokinin biosynthetic and metabolic enzymes play critical roles in maintaining cytokinin concentrations in individual cells or tissues at the proper level for normal functioning [11]. Isopentenyl transferases (IPT) catalyze a step in cytokinin biosynthesis, and cytokinin oxidase/dehydrogenases (CKXs) catalyze the irreversible degradation of cytokinin [12][13][14]. CKXs act on specific substrates such as isopentenyladenine and isopentenyladenine riboside [15]. These enzymes catalyze the irreversible degradation of cytokinins and their derivatives by removing their N 6 -substituted isoprene chains [16]. In this catalytic reaction, CKX enzymes display a dual catalytic mode in conjunction with molecular oxygen or other specific substances, which act as electron acceptors [17]. CKXs exhibit a relatively high degree of thermal stability, and the optimal pH for their activity depends on their natural substrates [18]. CKXs can function as flavin enzymes by covalently binding to flavin adenine dinucleotide (FAD) via a histidine residue, and they contain both FAD-and CK-binding domains [19]. Most CKX proteins undergo post-translational glycosylation, which might influence their molecular weights and localization [15]. CKX proteins are encoded by a small gene family [20]. CKX family genes and some putative members have been identified or functionally expressed in heterologous hosts, including various plant species such as Arabidopsis thaliana [5], Oryza sativa [21], Nicotiana tabacum [22], Zea mays [23,24], Triticum aestivum [25,26], Hordeum vulgare [25,26], Brassica rapa [27], Gossypium hirsutum [28], Setaria italica [29], Fragaria vesca [30], Dendrobium orchid [31], Pisum sativum [32] and Glycine max [33]. A. thaliana contains seven CKX genes (AtCKX1-AtCKX7) [5] with differing biochemical characteristics, subcellular localizations and expression patterns. Detailed analysis of transgenic Arabidopsis plants expressing each of these genes revealed that all of the transformants exhibited reduced cytokinin levels and developmental changes in roots and shoots, suggesting that the spatial and temporal expression patterns of CKXs might be related to the diverse functions of cytokinins [5,34]. Furthermore, overexpressing AtCKX1 and AtCKX2 in transgenic Centaurium erythraea led to decreased morphogenetic potential but had no effect on biomass production compared to the untransformed control [35]. The reduced expression of OsCKX2 in rice greatly increased cytokinin accumulation in the inflorescence meristem and increased the number of reproductive organs [21]. Similarly, silencing of HvCKX1 in barley led to lower CKX activity, thereby increasing seed yields and root weight [36].
Oilseed rape (Brassica napus L.), one of the most widely cultivated oilseed crops worldwide, was derived from natural hybridization between B. rapa and Brassica oleracea. A major goal in breeding oilseed rape is to enhance its yield and harvest index. The number of siliques per plant, number of seeds per pod and thousand-seed weight affect yield in this crop, and there is a significant positive correlation between silique length and seed weight [37,38]. The silique pericarp functions as a source organ to provide photosynthates to the developing seed; therefore, silique length may affect seed filling [39].
Numerous studies have shown that CKX genes strongly contribute to yield and stress responses [40]. In the current study, to identify genes that might help determine silique length and function in pod development in B. napus, we identified 23 CKX genes in B. napus and constructed a phylogenetic tree to evaluate the phylogenetic relationships between BnCKX proteins and those of other species. We also analyzed the properties and expression patterns of BnCKXs and identified differentially expressed genes (DEGs) between plants with long versus short siliques. Finally, we investigated the responses of selected BnCKX genes to hormone treatment.

Identification of CKXs in B. napus, B. rapa, and B. oleracea
The nucleotide and protein sequences of the AtCKXs were acquired from The Arabidopsis Information Resource database (https://www.arabidopsis.org/index.jsp). BnCKXs, BrCKXs and BoCKXs were identified through BLASTP analysis [41] of AtCKXs against the Brassica database (BRAD, http://brassicadb.org/brad/index.php) [42]. The CKX enzymes contain both a FAD-binding and a cytokinin-binding domain. Therefore, InterProScan (http://www.ebi.ac.uk/interpro/) was used to exclude any genes that did not contain these domains [43].

Multiple Sequence Alignment and Phylogenetic Analysis
Multiple sequence alignment of the protein sequences from the four species mentioned above was performed with MEGA7.0 [44] using default parameters. The conserved blocks of all predicted sequences were identified using the Gblocks program [45], and substitution saturation tests were performed with DAMBE program. To increase the accuracy of the predicted phylogenetic relationships of CKX family members of the four species, phylogenetic trees were constructed using three methods. A Neighbor-Joining tree (NJ tree) was generated via bootstrap analysis with 1000 replicates using MEGA7.0. A Maximum Likelihood tree (ML tree) was generated using the online program ATGC phyML with default parameters (http://www.atgc-montpellier.fr/phyml/) [46]. A Bayesian Inference tree (BI tree) was generated using the MrBayes program and ModelGenerator software [47]. Finally, the optimum topology was chosen using Tree-puzzle-5.3 and CONSEL software [48].

Chromosomal Locations, Gene Structures and Protein Profiles of BnCKXs
Chromosomal position information for the BnCKXs was obtained from the B. napus Genome Browser, and the genes were mapped onto chromosomal linkage groups with MapChart [49]. To investigate the exon/intron organization of the BnCKXs, the gene structures were analyzed using the Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/) [50]. Protein motifs were identified using MEME (http://meme-suite.org/tools/meme) [51]: the maximum number of motifs was set to 20 and default settings were used for other parameters. Annotations of the identified motifs were obtained from InterProScan. The isoelectric points (pI) and molecular weights (Mw) of the BnCKXs were calculated using the ExPASy proteomics server database (http://www.expasy.org/tools/) [52]. The subcellular localizations of the BnCKX proteins were determined using WoLF PSORT (http://www.genscript.com/ tools/wolf-psort) [53].

RNA-seq Analysis
To identify DEGs and to analyze the tissue-specific expression of the BnCKXs, publically available B. napus RNA-sequencing (RNA-Seq) data were downloaded and used to construct a heatmap with R-Studio. Data for the 23 BnCKX genes were extracted from the RNA-Seq database (accession number SRP072900), including expression data from eight different tissues collected from four B. napus materials with an extremely high or low harvest index grown in two different environments (Chongqing and Yunnan province, China) [57].

RNA Isolation and qRT-PCR Analysis
Total RNA was extracted with an RNeasy Extraction kit (Invitrogen, Carlsbad, CA, USA). First-strand cDNA was synthesized using M-MLV transcriptase (TaKaRa Biotechnology, Dalian, China) after DNase I treatment according to the manufacturer's instructions. Quantitative reverse-transcription polymerase chain reaction (qRT-PCR) was then conducted to measure the expression levels of the 23 BnCKXs in various tissues of ZS11. After analyzing the RNA-Seq data and the expression patterns of all BnCKXs in ZS11, qRT-PCR was also conducted to measure the expression of selected genes (BnCKX5-1, BnCKX 5-2, BnCKX 6-1 and BnCKX7-1) in seeds and silique pericarps at different stages in B. napus plants with extremely long versus short siliques. The qRT-PCR was performed as described by Ma [58]. Gene-specific primer pairs for the 23 BnCKXs were designed using Primer Premier 6.0 (Table S2) [59]. Melting curve analysis and agarose gel electrophoresis were performed to confirm the specificity of the product. Three biological replicates with three technical replicates were performed for each reaction. The BnActin7 gene served as an endogenous reference gene. The calculation of the gene relative expression levels followed the 2 −∆∆CT method described by Livik and Schmittgen. [60] 3. Results

Identification and Phylogenetic Analysis of BnCKXs
Using A. thaliana AtCKX protein sequences as queries, seven AtCKXs, 23 BnCKXs, 12 BoCKXs and 12 BrCKXs were identified via BLASTP analysis. Among the three phylogenetic trees, the BI tree had better topology than NJ, and ML trees and was therefore used to analyze the evolutionary relationships of CKX proteins. The 54 CKX proteins were grouped into seven distinct clades, with AtCKX1-7 and the corresponding B. rapa, B. oleracea, and B. napus subfamily members in each clade, hereafter referred to as groups I-VII, respectively. (Figure 1). Basic information about the BnCKXs is shown in Table 1. The encoded proteins vary from 336 (BnCKX2-2) to 768 (BnCKX2-1) amino acid (aa) residues in length, and the corresponding molecular weights of the BnCKXs range from 38.17 kDa (BnCKX2-2) to 87.02 kDa (BnCKX2-1). The BnCKXs also exhibit a wide range of pI values from 4.92 (BnCKX7-2) to 9.2 (BnCKX1-4). The pI values of proteins in groups I, II, VI and VII (except BnCKX7-3, pI value is 7.93) are less than 6.0, while those in group III and IV are greater than 7.0.

Chromosomal Localization
Chromosome mapping revealed that the BnCKXs are unevenly distributed throughout the chromosomes. Chromosome A02, A03, A10, C04 and C09 each contain two BnCKX genes, whereas chromosomes A04, A05, A09, C03, C06, C07, C08, A07-random, A09-random, A10-random, Ann-random, C07-random and Cnn-random contain only one. BnCKX1-3 and BnCKX1-4 (Group III) are located on chromosome C04, while BnCKX7-2 and BnCKX7-3 (Group VI) are located on chromosome A10 (Figure 2).  Figure 3 shows the gene structures of the BnCKXs. Overall, the closer the relationship between genes, the more similar their structures. Among the BnCKXs, the gene pairs BnCKX1-2 and BnCKX1-3, BnCKX4-1 and BnCKX4-2, BnCKX7-2 and BnCKX7-4, and BnCKX3-2 and BnCKX3-3 showed the highest similarity in terms of exon/intron number, distribution and length. However, BnCKX7-3 and BnCKX2-1 are strikingly different from other members of the same group in terms of exon numbers. BnCKX7-3 and BnCKX2-1 contained 18 and 11 exons, respectively. The MEME tool identified 20 conserved motifs in the BnCKXs (Figure 4). Motif 2 is present in all family members, whereas motifs 1, 6, and 7 are present in all family members except BnCKX2-2. Motifs 3, 4, 5, 8, 9 and 10 are present in all family members except BnCKX2-2. Motifs 13 and 20 are present only in BnCKX5-3 and BnCKX5-4, and motif 17 is found only in BnCKX5-7. Overall, except for BnCKX2-1, BnCKX2-2, and BnCKX7-3, the BnCKX proteins contain highly similar motifs in terms of type, order and number. We annotated the conserved motifs using the InterProScan program. Motifs 2 and 3 are associated with cytokinin dehydrogenase and the vanillyl-alcohol oxidase domain, whereas motifs 1, 2, 4, and 6 are annotated as FAD/cytokinin binding domains that can bind both FAD and cytokinins as substrates. The other motifs could not be annotated.  Figure 3 shows the gene structures of the BnCKXs. Overall, the closer the relationship between genes, the more similar their structures. Among the BnCKXs, the gene pairs BnCKX1-2 and BnCKX1-3, BnCKX4-1 and BnCKX4-2, BnCKX7-2 and BnCKX7-4, and BnCKX3-2 and BnCKX3-3 showed the highest similarity in terms of exon/intron number, distribution and length. However, BnCKX7-3 and BnCKX2-1 are strikingly different from other members of the same group in terms of exon numbers. BnCKX7-3 and BnCKX2-1 contained 18 and 11 exons, respectively. The MEME tool identified 20 conserved motifs in the BnCKXs (Figure 4). Motif 2 is present in all family members, whereas motifs 1, 6, and 7 are present in all family members except BnCKX2-2. Motifs 3, 4, 5, 8, 9 and 10 are present in all family members except BnCKX2-2. Motifs 13 and 20 are present only in BnCKX5-3 and BnCKX5-4, and motif 17 is found only in BnCKX5-7. Overall, except for BnCKX2-1, BnCKX2-2, and BnCKX7-3, the BnCKX proteins contain highly similar motifs in terms of type, order and number. We annotated the conserved motifs using the InterProScan program. Motifs 2 and 3 are associated with cytokinin dehydrogenase and the vanillyl-alcohol oxidase domain, whereas motifs 1, 2, 4, and 6 are annotated as FAD/cytokinin binding domains that can bind both FAD and cytokinins as substrates. The other motifs could not be annotated.  Subcellular localization analysis showed that seven of the 23 BnCKX proteins were predicted to be located in the vacuole, six in the mitochondria, five in the chloroplast, two in the nucleus, and two   Subcellular localization analysis showed that seven of the 23 BnCKX proteins were predicted to be located in the vacuole, six in the mitochondria, five in the chloroplast, two in the nucleus, and two Subcellular localization analysis showed that seven of the 23 BnCKX proteins were predicted to be located in the vacuole, six in the mitochondria, five in the chloroplast, two in the nucleus, and two in the cytosol. The remaining proteins might be excreted into the extracellular space. BnCKX1-5 and 1-6 might be located in the mitochondria, and the BnCKX3s are localized to the vacuole.

Cis-Acting Elements in the BnCKXs
Cis-acting elements in gene promoters, which serve as the binding sites of transcription factors, are important for regulating the initiation of gene transcription and for transcription efficiency. We identified 82 types of cis-acting elements in the promoters of the BnCKXs (Table S3). TATA-box and CAAT-box are core promoter elements, while light-responsive elements are the most varied cis-acting elements (25 types) and are present at the highest frequency (168 instances). Hormone-responsive cis-acting elements also account for a large proportion of elements, including auxin-responsive elements (TGA-element, AuxRR-core and TGA-box), MeJA-responsive elements (CGTCA-motif and TGACG-motif), gibberellin-responsive elements (GARE-motif, P-box and TATC-box), ethylene-responsive elements (ERE), salicylic acid-responsive elements (TCA-element) and abscisic acid-responsive elements (ABRE).
We also identified some environmental stress-responsive element in these genes. For instance, we identified: box-W1 and box E, which are fungal elicitor-responsive elements; ELI-box3, an elicitor-responsive element; WUN-motif, a wounding-responsive element; box S, a wounding-and pathogen-responsive element; LTR, involved in low-temperature responses; HSE, which functions in heat stress responses, MBS, a drought-inducible element and TC-rich repeats, which function in defense and stress responses.
Notably, some cis-acting elements in BnCKX might function in tissue-specific expression, such as the following: CAT-box and CCGTCC-box, which function in meristem-specific activation; Skn-1_motif and GCN4_motif, which are involved in endosperm expression; RY-element, involved in seed-specific regulation; as 1, involved in root-specific expression and HD-Zip 1 and HD-Zip 2 control palisade mesophyll cell differentiation and leaf morphology, respectively.

RNA-Seq Analysis
As shown in the heatmaps ( Figure 5), BnCKX1, BnCKX2, BnCKX4 exhibited almost no expression in any tissue. Other BnCKX family members were expressed at high or low levels in different tissues. For example, BnCKX5-1 and 5-2 were highly expressed in seeds and silique pericarps, while BnCKX6-1 and 6-2 were mainly expressed in leaves, stems and silique pericarps but were not expressed in seeds or buds. BnCKX7-1 was highly expressed in stems and leaves. Members in the same subfamily also showed diverse expression patterns. For example, BnCKX3-2 and 3-3 were not expressed in any tissues, while BnCKX3-1 and 3-4 were expressed specifically in buds and seeds. The expression patterns of BnCKX7-2 and 7-4 were very similar; neither was expressed in seeds or stems. By contrast, BnCKX7-3 was expressed in all tissues examined at similar levels. In addition, eight DEGs were found in plants with different yields or harvest indices. Notably, the expression levels of BnCKX5-1, 6-1, 6-2, and 7-1 in siliques were significantly different in the siliques of plants with different harvest indices.

Expression of selected BnCKX genes in different plant materials and in response to IAA and 6-BA treatment
Based on the RNA-Seq and ZS11 qRT-PCR data, we selected four genes (BnCKX5-1, 5-2, 6-1, and 7-1) that might function in silique development and analyzed their expression patterns in 'ZS4' and 'NY12' (the average silique length can be seen in Table S1) via qRT-PCR. Variance analysis revealed that the expression levels of these four genes showed significant differences or extremely significant (p value < 0.05 or 0.01) between material 'ZS4' and 'NY12' (Figure 7 and Table S4), suggesting that BnCKX5-1, 5-2, 6-1, and 7-1 might help determine the length and development of siliques. Notably, BnCKX5-1 and 7-1 were more highly expressed in 14S, 21S, 30S in 'ZS4', whereas BnCKX5-2 was more highly expressed in 7S, 21S, 30S in 'NY12'. Moreover, BnCKX6-1 was mainly expressed in the silique pericarp, with extremely significant differences in expression between these two materials. BnCKX6-1 was strongly expressed in 14P and 21P in 'ZS4' but in 7P and 30P in 'NY12'. Perhaps lower expression levels of BnCKX6-1 during the early stages of silique development and higher expression levels in the middle-stage of this process lead to greater final silique lengths.

Expression of Selected BnCKX Genes in Different Plant Materials and in Response to IAA and 6-BA Treatment
Based on the RNA-Seq and ZS11 qRT-PCR data, we selected four genes (BnCKX5-1, 5-2, 6-1, and 7-1) that might function in silique development and analyzed their expression patterns in 'ZS4' and 'NY12' (the average silique length can be seen in Table S1) via qRT-PCR. Variance analysis revealed that the expression levels of these four genes showed significant differences or extremely significant (p value < 0.05 or 0.01) between material 'ZS4' and 'NY12' (Figure 7 and Table S4), suggesting that BnCKX5-1, 5-2, 6-1, and 7-1 might help determine the length and development of siliques. Notably, BnCKX5-1 and 7-1 were more highly expressed in 14S, 21S, 30S in 'ZS4', whereas BnCKX5-2 was more highly expressed in 7S, 21S, 30S in 'NY12'. Moreover, BnCKX6-1 was mainly expressed in the silique pericarp, with extremely significant differences in expression between these two materials. BnCKX6-1 was strongly expressed in 14P and 21P in 'ZS4' but in 7P and 30P in 'NY12'. Perhaps lower expression levels of BnCKX6-1 during the early stages of silique development and higher expression levels in the middle-stage of this process lead to greater final silique lengths. To gain more insight into the roles of B. napus CKX genes in stress responses, we used qRT-PCR to examine the expression of four DEGs in response to treatment with the synthetic cytokinin 6-BA or the auxin IAA ( Figure 8). BnCKX5-1, 5-2, and 6-1 were moderately induced by 6-BA. The expression patterns of BnCKX5-1 and 5-2 were similar, both were slightly induced at the beginning of 6-BA treatment and reached a peak at 3 h, followed by a decrease, and another peak at 24 h. Perhaps this phenomenon is related to the complex dynamic balance that shifts after application of the exogenous hormone. By contrast, BnCKX7-1 was downregulated under this treatment. Under IAA treatment, BnCKX5-1 and 5-2 were significantly induced at 3 h, while BnCKX7-1 was slightly downregulated, followed by upregulation, reaching a peak at 12 h, and was again downregulated at 24 h. Compared with 6-BA treatment, IAA more strongly induced BnCKX6-1 expression, with a >10fold increase in expression at 3 h.  To gain more insight into the roles of B. napus CKX genes in stress responses, we used qRT-PCR to examine the expression of four DEGs in response to treatment with the synthetic cytokinin 6-BA or the auxin IAA ( Figure 8). BnCKX5-1, 5-2, and 6-1 were moderately induced by 6-BA. The expression patterns of BnCKX5-1 and 5-2 were similar, both were slightly induced at the beginning of 6-BA treatment and reached a peak at 3 h, followed by a decrease, and another peak at 24 h. Perhaps this phenomenon is related to the complex dynamic balance that shifts after application of the exogenous hormone. By contrast, BnCKX7-1 was downregulated under this treatment. Under IAA treatment, BnCKX5-1 and 5-2 were significantly induced at 3 h, while BnCKX7-1 was slightly downregulated, followed by upregulation, reaching a peak at 12 h, and was again downregulated at 24 h. Compared with 6-BA treatment, IAA more strongly induced BnCKX6-1 expression, with a >10-fold increase in expression at 3 h. To gain more insight into the roles of B. napus CKX genes in stress responses, we used qRT-PCR to examine the expression of four DEGs in response to treatment with the synthetic cytokinin 6-BA or the auxin IAA ( Figure 8). BnCKX5-1, 5-2, and 6-1 were moderately induced by 6-BA. The expression patterns of BnCKX5-1 and 5-2 were similar, both were slightly induced at the beginning of 6-BA treatment and reached a peak at 3 h, followed by a decrease, and another peak at 24 h. Perhaps this phenomenon is related to the complex dynamic balance that shifts after application of the exogenous hormone. By contrast, BnCKX7-1 was downregulated under this treatment. Under IAA treatment, BnCKX5-1 and 5-2 were significantly induced at 3 h, while BnCKX7-1 was slightly downregulated, followed by upregulation, reaching a peak at 12 h, and was again downregulated at 24 h. Compared with 6-BA treatment, IAA more strongly induced BnCKX6-1 expression, with a >10fold increase in expression at 3 h.

Discussion
The B. napus (AACC, 2n=38) genome comprises the B. rapa and B. oleracea genomes, which descended from a common ancestor [61,62]. Genome-wide analysis and expression profiling of the CKX gene family have been performed in Arabidopsis and the diploid B. rapa, revealing 7 and 12 CKX genes in Arabidopsis and B. rapa, respectively [5,27]. By contrast, we identified only 23 BnCKXs from the tetraploid B. napus in the current study, and it contains only 2-6 copies homologs of each AtCKX genes due to gene loss and genome shrinkage. By comparing the number of BnCKX and BoCKX genes identified in each subfamily in the present work with those identified by Liu et al in B. rapa, we found that the number of BrCKX and BoCKX genes identified in each subfamily is equal and that the number of BnCKXs is about the sum of BrCKXs and BoCKXs. In this study, only 3 BnCKX2 genes were identified, whereas B. rapa contains 2 BrCKXs and B. oleracea contains 2 BoCKXs. These results suggested that there may be gene elimination in BnCKX2s. By constructing a phylogenetic tree and predicting the gene structures, protein profiles, and conserved motifs of the BnCKXs, we found that genes belonging to the same subfamily exhibit similar features. The structures and conserved motifs of these genes might be closely related to their expression patterns. For example, BnCKX7-1, 7-2, and 7-4 share similar gene structures and conserved motifs, as well as expression levels and tissue-specific expression patterns. However, the gene structure and conserved motifs of BnCKX7-3 are quite different from those of the other CKX members, as is its expression pattern, pointing to a functional division in this gene family.
To better understand the functions and interactions of the 23 BnCKX proteins, we predicted their subcellular localizations. The BnCKX proteins are localized in different cellular compartments, including the chloroplast, cytosol, mitochondria, vacuole, and nucleus. In addition, BnCKX3 proteins contain vacuole-targeting peptides, and AtCKX3-GFP is localized to the vacuole [5]. By contrast, BnCKX6 and BnCKX1-2, -4, -5 and -6 are localized to the mitochondria, indicating that BnCKX6 likely interacts with BnCKX1 in the mitochondria.
Analysis of cis-elements in the promoters of the BnCKXs suggested that they play a role in regulating biotic and abiotic stress responses, as well as plant growth. Indeed, CKX1 expression is induced by cytokinins, abscisic acid, and abiotic stress in maize [63], and CKXs are thought to function in the crosstalk between gibberellic acid and cytokinin signaling, thereby regulating plant growth and development [64,65]. In the current study, we predicted that six types of hormone-responsive cis-acting elements are present in the promoters of BnCKX genes, including auxin-responsive, jasmonate-responsive, gibberellin-responsive, ethylene-responsive, and salicylic acid-responsive elements. However, we did not detect cytokinin-responsive elements, perhaps because only a few such elements have thus far been identified in plants. CKX-overexpressing plants exhibit increased drought and salinity tolerance, as well as enhanced heat stress tolerance [66,67]. In the current study, we identified MBS elements in all BnCKX genes, as well as HSE and TCA-elements in 15 and 10 BnCKX genes, respectively. These results suggest that most BnCKX genes are associated with plant growth and stress responses.
Based on the RNA-Seq data (Table S5), we constructed two heat maps, revealing the diverse gene expression patterns of the BnCKXs. In general, the expression patterns of BnCKXs in two environments were highly similar, with expression of only half of the 23 BnCKX genes detected. However, the expression patterns of members of the same subfamily sometimes differed, such as for the BnCKX3 and BnCKX7 subfamilies, suggesting some functional diversification among gene family members. In addition, most BnCKX genes were highly expressed in reproductive organs such as buds, flowers or siliques. In Arabidopsis, ckx3 and ckx5 double mutants form large inflorescences and floral meristems and contain supernumerary ovules, thus leading to an increase in seed set per silique [68]. There were also some discrepancies between the expression patterns of RNA-seq data and quantitative polymerase chain reaction (qPCR) analysis, especially the expression level of BnCKX1-2, 1-3, 2-4 and 3-4. One possible reason is that the total RNA for RNA-Seq was pooled from four types B.napus lines with extremely high and low harvest index, while the samples for qRT-PCR analysis were collected from 'ZS11'. On the other hand, the seeds and silique pericarps for RNA-Seq were harvested at 20 days after flowering, and the mature leaves, buds and stems were harvested at flowering stage 63-65. These discrepancies also indicated the complex expression profiles of BnCKX genes in different materials and plant development stages. In the current study, we found that the expression level of BnCKX5 markedly differed not only at different stages of pod development but also in different plant materials. Thus, we hypothesize that BnCKX5 might play an essential role in modulating cytokinin levels, thus regulating silique development. Notably, BnCKX7-1 was downregulated following exogenous 6-BA treatment, and Arabidopsis overexpressing AtCKX7 showed increased CKX enzymatic activity, ultimately leading to cytokinin deficiency [69]. Based on these findings, we hypothesize that cytokinin deficiency is the ultimate result of CKX7-1 expression and that the process of cytokinin degradation by CKX enzymes is a dynamic and complex process.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/9/3/168/s1, Table S1: Silique length of extremely materials at several stages of pod development, Table S2: Primers used for quantitative reverse-transcription polymerase chain reaction (qRT-PCR) analysis of BnCKX genes, Table S3: Stress-and hormone-responsive elements in the promoter regions of BnCKX genes, Table S4: The expression levels differences of selected four genes by variance analysis, Table S5: RNA-Seq data of 23 BnCKX genes.