Genome-Wide Survey and Expression Analysis of the KT/HAK/KUP Family in Brassica napus and Its Potential Roles in the Response to K+ Deficiency

The KT/HAK/KUP (HAK) family is the largest potassium (K+) transporter family in plants, which plays key roles in K+ uptake and homeostasis, stress resistance, and root and embryo development. However, the HAK family has not yet been characterized in Brassica napus. In this study, 40 putative B. napus HAK genes (BnaHAKs) are identified and divided into four groups (Groups I–III and V) on the basis of phylogenetic analysis. Gene structure analysis revealed 10 conserved intron insertion sites across different groups. Collinearity analysis demonstrated that both allopolyploidization and small-scale duplication events contributed to the large expansion of BnaHAKs. Transcription factor (TF)-binding network construction, cis-element analysis, and microRNA prediction revealed that the expression of BnaHAKs is regulated by multiple factors. Analysis of RNA-sequencing data further revealed extensive expression profiles of the BnaHAKs in groups II, III, and V, with limited expression in group I. Compared with group I, most of the BnaHAKs in groups II, III, and V were more upregulated by hormone induction based on RNA-sequencing data. Reverse transcription-quantitative polymerase reaction analysis revealed that the expression of eight BnaHAKs of groups I and V was markedly upregulated under K+-deficiency treatment. Collectively, our results provide valuable information and key candidate genes for further functional studies of BnaHAKs.


Introduction
Potassium (K + ) is the most abundant nutrient element in plants, accounting for 10% of the dry weight [1], and plays essential roles in protein synthesis [2], phloem transport [3], osmoregulation [4], and plant development [5]. K + can also enhance the resistance to biotic and abiotic stress for improving the plant growth status. Therefore, plants often exhibit signs of chlorosis and necrosis under a state of K + deficiency, resulting in a functional decline in leaf photosynthesis, and consequently compromising crop yield and quality [6,7]. The soil K + concentration is generally considered to vary between 0.1 and 1 mM [8]; however, a large portion of cultivated land worldwide remains K + -deficient [9,10]. Consequently, the proper use of K + fertilizer and increasing utilization of K + will contribute to sustainable crop production and high-quality crops. Brassica napus is an important oil crop that requires

Identification of HAK Family Members in the B. napus Genome
To identify the HAK-encoding genes in the B. napus genome, we performed a preliminary BLASTP search using the Arabidopsis HAK proteins (AtKUPs) [31] as queries. After removing the redundant and severely deleted/missing sequences, and verifying the sequences using the Pfam tool, we finally obtained a total of 40 putative BnaHAKs with complete functional domains. To distinguish the candidates, these genes were labeled BnaHAK01 to BnaHAK40 according to their chromosome locations. Using the same method, we also identified 21 and 23 non-redundant HAKs in B. rapa (BrHAKs) and B. oleracea (BoHAKs), respectively (Table S1).
Physicochemical property analysis showed that the length of BnaHAK proteins (BnaHAKs) ranged from 402 to 864 amino acids. The molecular weight of the BnaHAKs ranged from 44.40 to 171.55 kDa, and the isoelectric points were concentrated at 5. 22-9.44. Subcellular localization analysis showed that all 40 BnaHAKs are present in the vacuoles, and approximately half of them were also localized on the cell membrane (Table 1), which is consistent with their role in K + uptake and homeostasis.
Overall, we identified a relatively larger number of HAKs in the B. napus genome compared with that reported for many other species [38]. This might be attributed to the fact that B. napus (AACC. n = 19) is a "young" species, which was produced by the recent hybridization between B. rapa (AA. n = 10) and B. oleracea (CC. n = 9).
On the basis of the topology and bootstrap values of the NJ tree, HAK proteins from these four Brassicaceae species were classified into four groups: groups I-III and V ( Figure 1). This confirmed that group IV was lost in the Brassicaceae lineage during evolution [38,39]. Group I consisted of two BnaHAKs (5%), one BrHAK, one BoHAK, and one AtKUP; group II was the largest, including 18 BnaHAKs (45%), nine BrHAKs, 12 BoHAKs, and six AtKUPs; group III had 14 BnaHAKs (35%), eight BrHAKs, seven BoHAKs, and three AtKUPs; and group V contained six BnaHAKs (15%), three BrHAKs, three BoHAKs, and three AtKUPs ( Figure 1). In the phylogenetic tree, the number of BnaHAKs was double that of AtKUPs in groups I and V; however, it was three times that of AtKUPs in group II and was 4.5 times the number in group III. This suggested that the amplification rate in distinct groups varies, which might reflect the specific functional needs of each group during evolution.
On the basis of the topology and bootstrap values of the NJ tree, HAK proteins from these four Brassicaceae species were classified into four groups: groups I-III and V ( Figure 1). This confirmed that group IV was lost in the Brassicaceae lineage during evolution [38,39]. Group I consisted of two BnaHAKs (5%), one BrHAK, one BoHAK, and one AtKUP; group II was the largest, including 18 BnaHAKs (45%), nine BrHAKs, 12 BoHAKs, and six AtKUPs; group III had 14 BnaHAKs (35%), eight BrHAKs, seven BoHAKs, and three AtKUPs; and group V contained six BnaHAKs (15%), three BrHAKs, three BoHAKs, and three AtKUPs (Figure 1). In the phylogenetic tree, the number of BnaHAKs was double that of AtKUPs in groups I and V; however, it was three times that of AtKUPs in group II and was 4.5 times the number in group III. This suggested that the amplification rate in distinct groups varies, which might reflect the specific functional needs of each group during evolution.

Gene Structure and Intron Pattern of BnaHAKs
Gene structural diversity provides an important clue into the function and evolutionary history of multi-gene families. Thus, we analyzed the gene structures of the identified BnaHAKs using the online software GSDS 2.0.
The intron number varied widely among the 40 BnaHAKs, ranging from 4 in BnaHAK37 to 11 in BnaHAK18 ( Figure 2B). The exon-intron structures of groups I, III, and V were relatively conserved compared with those of group II, with only a few exceptions that might have been due to the low genome sequence quality. The two BnaHAKs in group I contained 10 and 8 introns, respectively; the members of group III generally possessed 7-9 introns, except for BnaHAK16, BnaHAK36, and BnaHAK37, which contained 4 or 5 introns; group V contained 8-9 introns; and 18 BnaHAKs in group II contained 5-11 introns ( Figure 2B). Notably, although the intron number was variable, the intron insertion sites and phase of most introns (10 introns) were conserved across the different groups ( Figure 2C). Among them, intron 4 was completely conserved in all four groups; introns 3 and 6-9 were highly conserved in the four groups; and introns 1 and 10 were only conserved in groups I, and III and V respectively ( Figure 2C). The insertion site of intron 2 was slightly variable, resulting in independent conserved patterns in the four groups, with three typical types found in group II ( Figure 2C). Similarly, intron 5 contained two insertion types in group II. Moreover, the intron insertion patterns of AtKUPs, BrHAKs, and BoHAKs further validated our results in BnaHAKs, indicating their conservation during evolution ( Figure S1). Therefore, we speculated that these conserved intron insertion sites might exist in HAKs of plant species. In addition, these introns were commonly inserted outside of the TM domains, except for intron 6, which was highly conserved in the same TM domain across HAKs. In addition, there were some unconserved intron insertion sites identified that were generally concentrated on the C-terminal ( Figure S2).
Overall, these results revealed that the intron number of BnaHAKs differs among different groups or even within a group, but that most of the intron insertion sites are highly conserved.

Chromosomal Distribution and Duplication of BnaHAKs
To understand the duplication relationship of the HAK gene family in B. napus, we analyzed the chromosomal location and collinearity relationship of the BnaHAKs.
Our results showed that 38 of the 40 BnaHAKs (the detailed locations of BnaHAK39 and BnaHAK40 are currently unknown) were distributed unevenly across 16 of the 19 B. napus chromosomes, except for A n 06, A n 10, and C n 05 ( Figure 3). The chromosomes A n 01, A n 03, and C n 01 possessed the most BnaHAKs (4 genes), whereas chromosomes A n 04, A n 09, C n 08, and C n 09 only had one BnaHAK each. Moreover, the number of BnaHAKs between the A n (18 genes) and C n (12 genes) subgenomes also appeared to be uneven ( Figure 3). Collinear relationship analysis showed that 35 of the 40 BnaHAKs had a collinear relationship with B. napus, B. rapa, and/or B. oleracea homologs (Table S2). Among these 35 BnaHAKs, 13 (37%) were derived from whole-genome duplication (WGD), with nine genes originating from B. rapa and four from B. oleracea. The remaining BnaHAKs (22) originated from small-scale duplication events in the B. napus genome, including 12 (34%), due to segmental exchange (SE) events, six (17%) from homologous exchange (HE) events, and four (12%) from segmental duplication (SD) events. No tandem duplication (TD) events were noted in the BnaHAKs. Interestingly, all members in group V were derived from small-scale duplications, displaying a special evolutionary law. Moreover, most of these small-scale duplicated genes were also dated from B. rapa, suggesting a bias demand and/or retention for the HAKs from B. rapa in B. napus. The non-synonymous (Ka) to synonymous (Ks) substitution rate ratios of the collinear gene pairs were all less than 1 (Table S3), indicating that the duplicated BnaHAKs mainly underwent purification selection during evolution. the BnaHAKs. Interestingly, all members in group V were derived from small-scale duplications, displaying a special evolutionary law. Moreover, most of these small-scale duplicated genes were also dated from B. rapa, suggesting a bias demand and/or retention for the HAKs from B. rapa in B. napus. The non-synonymous (Ka) to synonymous (Ks) substitution rate ratios of the collinear gene pairs were all less than 1 (Table S3), indicating that the duplicated BnaHAKs mainly underwent purification selection during evolution.  BnaHAKs were mapped onto 16 chromosomes, except for BnaHAK39 and BnaHAK40 whose chromosomal locations are currently unclear. The C nn chromosomes represent chromosomal fragments mapped to the C n subgenome, but the detailed locations are still unclear. A black frame represents the genes that originated from a whole-genome duplication event; the grey lines and grey background represent the genes that originated from segmental exchange duplication; the black lines represent the genes that originated from homologous exchange; and the dashed lines represent the genes that originated from segmental duplication.
In summary, these results demonstrated that both a WGD event and small-scale duplication events (SE, HE, and SD) were the major driving forces for the large HAK expansion in B. napus; in particular, the BnaHAKs derived from B. rapa tended to be retained in the B. napus genome.

Regulation Mechanism in the Promoter Regions of BnaHAKs
We predicted the potential transcriptional regulators of the BnaHAKs in the PlantTFDB database using the promoter sequences (-1500 bp), and constructed a regulation network based on the results ( Figure 4). Overall, we identified 269 potential TFs belonging to 26 TF families that might bind to the promoter regions of 37 BnaHAKs (no TF was predicted on the promoters of BnaHAK09, BnaHAK13, and BnaHAK19) ( Table S4). The most abundant TFs belonged to the ERF (46 genes), Dof (45 genes), and WRKY (37 genes) families ( Figure 4A). In addition, the Dof, MIKC_MADS, and B3 families regulated the most BnaHAKs, 19, 10, and 9, respectively ( Figure 4B), suggesting their important roles in regulating the expression of BnaHAKs. The remaining TF families (including MYB, G2-like, and Trihelix) might be rare, and appear to regulate the expression of only a few BnaHAKs. In most cases, these TFs might bind to multiple potential BnaHAKs in different groups, with a few TFs only targeting a specific BnaHAK ( Figure 4B). For instance, the ZF-HD TF might only bind to the promoter of BnaHAK15, SRS TF might only bind to BnaHAK40, and HD-ZIP TF might only bind to BnaHAK04 ( Figure 4B).
families regulated the most BnaHAKs, 19, 10, and 9, respectively ( Figure 4B), suggesting their important roles in regulating the expression of BnaHAKs. The remaining TF families (including MYB, G2-like, and Trihelix) might be rare, and appear to regulate the expression of only a few BnaHAKs. In most cases, these TFs might bind to multiple potential BnaHAKs in different groups, with a few TFs only targeting a specific BnaHAK ( Figure 4B). For instance, the ZF-HD TF might only bind to the promoter of BnaHAK15, SRS TF might only bind to BnaHAK40, and HD-ZIP TF might only bind to BnaHAK04 ( Figure 4B). To further explore the regulatory mechanism of BnaHAKs, we also analyzed the potential ciselements in the promoter regions of BnaHAKs using PlantCARE online software. Accordingly, we identified a total of 61 types comprising 2717 cis-elements in the 40 BnaHAK promoters (Table S5). In addition to the most common basic core and light-responsive elements, we identified nine types of hormone-responsive cis-elements, including TGA-Element, ABRE, and P-Box, suggesting that the expression of BnaHAKs might be regulated by diverse hormones. Moreover, seven other types of cis-elements (e.g., LTRs, AREs, TC-rich repeats) were identified that are involved in biotic and abiotic responses such as low-temperature, anaerobic induction, defense and stress, and wound responses ( Figure 4C).
These results showed that most BnaHAKs might be regulated by multiple TFs, implying a complex regulatory network for BnaHAKs in K + uptake and transport in B. napus.

Comprehensive Analysis of miRNAs Targeting BnaHAKs
Several studies have shown that miRNAs play a significant role in targeting gene expression. Therefore, we predicted the potential miRNAs of BnaHAKs using the psRNATarget online software. The results indicated that 38 of the 40 BnaHAKs might be the targets of 145 miRNAs (no miRNA was predicted to target BnaHAK15 and BnaHAK27) (Table S6), revealing the important roles of these miRNAs in BnaHAKs expression. The candidate miRNAs belonged to 79 families. Among them, miR2673 appears to regulate the greatest number of BnaHAKs (11), with the others regulating 1-7 BnaHAKs. Moreover, in groups II, III, and V, an individual BnaHAK gene could be the target of multiple miRNAs. For example, BnaHAK05 in group II, BnaHAK16 in group III, and BnaHAK04 in group V might be the targets of eight, seven, and four miRNAs, respectively. Conversely, the two members in group I, BnaHAK14 and BnaHAK32, might be targeted only by a single miRNA, miR8001 and miR8024, respectively. In addition, the complementarity of the miRNAs and their targeted BnaHAKs in the four groups were all mainly between 70% and 85% (Table S5), indicating that the sequences of miRNAs had low rates of differentiation and mutation across distinct groups during evolution. These findings provide valuable information for further elucidating the function of miRNAs in regulating BnaHAKs expression.

Characteristics of the Predicted BnaHAK Proteins
We next explored the dimensional structure characteristics of the BnaHAKs, using the Protter online tool and SWISSMODEL online server to predict the secondary and tertiary structures, respectively. All 40 BnaHAKs had highly conserved secondary structures, which were commonly equipped with a conserved N-terminal TM domain and a C-terminal cytoplasmic domain, with a few exceptions that might be attributed to the deletion or lack of protein sequences ( Figure 5A). Coincident with the secondary structures, the tertiary structures of BnaHAKs were also highly conserved. Significantly, there were two conserved types of the C-terminal cytoplasmic domain in groups II, III, and V identified (type "a" and type "b"), whereas the members of group I had only the type "a" domain ( Figure 5B), suggesting the common origin of type "a" during evolution. In most cases, each TM domain contained 12 typical TM helices with an obvious topologically inverted repeat. Consistent with a previous study [35], these helices traversed the intra-and extracellular sections of the membranes, and the final helix extended into the cytoplasm and connected to the cytoplasmic domain ( Figure 5B,C). Moreover, each TM domain generally had three consecutive K + -binding sites that were highly conserved in groups I and V, but 1-3 K + -binding site(s) were lost in groups II and III due to one or more key amino acid residue substitution(s) in these regions ( Figure 5D). Considering the crucial roles of the K + -binding sites for BnaHAKs in K + absorption, we hypothesized that these substitutions might indicate functional diversity. Furthermore, introns 3, 6, and 9 of BnaHAKs were conservatively distributed in the three K + -binding sites, implying that they might affect the function of these regions. absorption, we hypothesized that these substitutions might indicate functional diversity. Furthermore, introns 3, 6, and 9 of BnaHAKs were conservatively distributed in the three K + -binding sites, implying that they might affect the function of these regions. Overall, our results suggested that the dimensional structures of the BnaHAKs are highly conserved, and substitutions in the K + -binding site(s) might reflect their structural and even functional diversity during evolution.

Spatiotemporal Expression Profiles of BnaHAKs across Different Developmental Stages
On the basis of published RNA-Seq data (BioProject: PRJNA358784), we examined the Overall, our results suggested that the dimensional structures of the BnaHAKs are highly conserved, and substitutions in the K + -binding site(s) might reflect their structural and even functional diversity during evolution.

Spatiotemporal Expression Profiles of BnaHAKs across Different Developmental Stages
On the basis of published RNA-Seq data (BioProject: PRJNA358784), we examined the spatio-temporal expression of 40 BnaHAKs in 50 B. napus tissues from seven organs at five developmental stages (seeding, budding, initial flowering, full-bloom, and silique stages). We excluded four BnaHAKs from the heatmap with no detectable transcript or weak expression levels [fragments per kilobase of exon model per million reads mapped (FPKM) < 1].
All BnaHAKs showed clear spatial and temporal expression characteristics, and the expression patterns of the four groups were quite different. The expression profile of group I was relatively narrow; the two members of this group, BnaHAK14 and BnaHAK32, had detectable transcript levels in the root and/or seed tissues ( Figure 6A). By contrast, most of the genes in groups II and III exhibited broad expression profiles among the 50 tissues investigated ( Figure 6A), indicating their wide-reaching roles in B. napus. The expression patterns of group V were divided into three categories, which were highly expressed in the seed, leaf, and stem tissues ( Figure 6A). In general, the expression patterns of the homologs in each group were generally different, such as BnaHAK14 and BnaHAK32 in group I, implying their functional divergence during evolution. However, the expression profiles of sister pair genes (i.e., homologs with collinear relationships) were generally similar, with Pearson correlation coefficients greater than 0.8 (Table S7), including BnaHAK08/BnaHAK26, BnaHAK12/BnaHAK29, and BnaHAK04/BnaHAK22 ( Figure 6A), suggesting possible functional redundancy of these duplication pairs.
To explore the hormone-induced expression patterns of the 40 BnaHAKs, we examined their transcript levels in the seedling roots of the B. napus ZS11 ecotype under five exogenous hormone treatments (IAA, ACC, ABA, GA 3 , and 6-BA). Four BnaHAKs were excluded from the heatmap as they had no or weak expression levels (FPKM < 1). The members of group I were less sensitive to hormone treatments, as only BnaHAK32 was induced by IAA, ABA, and GA 3 treatments to different degrees ( Figure 6B). In contrast, the expression of the majority of members in groups II, III, and V was clearly upregulated under the five hormone treatments ( Figure 6B). In groups II and V, the expression of more than half of the genes was upregulated by IAA, ACC, and GA 3 treatments ( Figure 6B). In group III, the expression of most of the genes was upregulated by IAA, ABA, and GA 3 induction ( Figure 6B). These results demonstrated that most of the BnaHAKs in these four groups were sensitive to exogenous IAA, ACC, and GA 3 , and a few genes could also be induced by ABA and 6-BA. In addition, the sensitivity of each group to the five hormone treatments was generally divergent.

Expression Levels of BnaHAKs under Low-K + Conditions
The expression levels of HAKs are generally regulated by the concentration of K + . Among the four groups, group I is mainly involved in high-affinity K + uptake, whereas most of the genes in groups II and III encode low-affinity transporters. Several genes in group V have also been related to K + uptake under low-K + (-K + ) conditions [20,31]. Thus, to gain insight into the potential roles of BnaHAKs in response to -K + stress in B. napus seedling roots, eight BnaHAKs in group I (BnaHAK14/BnaHAK32) and V (BnaHAK07/BnaHAK25, BnaHAK04/BnaHAK22, and BnaHAK01/BnaHAK19) were selected for reverse transcription-quantitative polymerase chain reaction (RT-qPCR) assays.
As shown in Figure 7, the eight genes were all upregulated to varying degrees at different times under the -K + condition. The expression patterns of these proteins could be divided into two main types: members of group I showed an upregulation trend at the six time points with significantly higher levels at the early stages, whereas members of group V exhibited a polarized pattern with a transient upregulation pattern under K + starvation. This pattern indicated that group I might play an important role in K + uptake under -K + conditions. Notably, the expression profiles of the sister pairs of BnaHAK14/BnaHAK32 and BnaHAK04/BnaHAK22 were very similar, whereas those of BnaHAK07/BnaHAK25 and BnaHAK01/BnaHAK19 were quite different. This indicated the existence of functional differentiation of the duplicated genes during evolution in B. napus. In summary, the expression of group I and V genes was sensitive to the -K + environment, which provides a foundation for further exploration of their roles in K + uptake in the roots.

Discussion
Given their key roles in plant K + uptake, homeostasis, translocation, and stress resistance, the HAK family has been identified in many plant species such as Arabidopsis [31], rice [14], wheat [32], and Zea mays (maize) [40]. Previously, the plant HAK family was generally classified into four groups (I-IV), such as in rice and maize [14,40]. However, recent studies further divided HAKs into five groups, with the newly defined Group V separated from the old group III [38,39]. In this study, we first systematically identified the HAK family in an important oil crop, B. napus, and identified 40 members, which constitutes the second largest HAK family identified in plants to date [19,38]. In summary, the expression of group I and V genes was sensitive to the -K + environment, which provides a foundation for further exploration of their roles in K + uptake in the roots.

Discussion
Given their key roles in plant K + uptake, homeostasis, translocation, and stress resistance, the HAK family has been identified in many plant species such as Arabidopsis [31], rice [14], wheat [32], and Zea mays (maize) [40]. Previously, the plant HAK family was generally classified into four groups (I-IV), such as in rice and maize [14,40]. However, recent studies further divided HAKs into five groups, with the newly defined Group V separated from the old group III [38,39]. In this study, we first systematically identified the HAK family in an important oil crop, B. napus, and identified 40 members, which constitutes the second largest HAK family identified in plants to date [19,38]. Based on our phylogenetic, intron insertion pattern, and protein characteristic analyses, we confirmed the existence of group V, which supports this new classification of the plant HAK family.
In previous studies, SD and TD events were reported as the major contributors to the large HAK gene expansion in many plant genomes [40][41][42][43]. However, in this study, we found that not only small-scale duplication events (including SE, 34%; HE, 17%; and SD, 12%) but also genome-wide duplication (hybridization between B. rapa and B. oleracea) were the major contributors to the rapid gene expansion of the BnaHAK family. This same scenario was also recently demonstrated in the sugarcane HAK family [33]. Given that Brassicaceae species have experienced a common whole-genome triplication (WGT) event, the 13 AtKUPs should be expanded to~40 and~80 genes in the B. rapa/B. oleracea and B. napus genomes, respectively. However, we found that 48% (19) of the BrHAKs and 43% (17) of the BoHAKs were lost after the WGT event, while most of the duplicated BnaHAKs were retained after the WGD event (hybridization). This might be attributed to the fact that the WGD event between B. rapa and B. oleracea occurred only about~7500 years ago [44], which might be too short for the loss of duplicates.
Since the first HAK gene was cloned from barley, many homologs have been identified and functionally analyzed in various plant species (Table S8). The HAK genes in groups I and V were proven to be mainly involved in high-affinity K + uptake and translocation. For example, AtHAK5 (I) [16], HvHAK1 (I) [15], LeHAK5 (I) [17], and AtKUP7 (V) [19] can respond to K + deficiency, maintaining K + uptake and transport. Similarly, we found that the expression of BnaHAKs in groups I and V was upregulated under the -K + condition, suggesting a similar function of these genes in K + uptake and transport. Interestingly, according to their expression patterns, BnaHAK14 and BnaHAK32 seem to have complementary functions of K + uptake with each other in the roots; namely, BnaHAK14 was highly expressed and absorbed K + under control conditions, whereas the expression of BnaHAK32 was greatly upregulated to maintain K + absorption in the -K + environment.
Moreover, due to more extensive study of group I, other functions of group I members have also been identified. For instance, AtHAK5 [45], SIHAK5 [46], and OsHAK1 [47] in group I can mediate high-affinity Cs + uptake because of the similarity between Cs + and K + . Some genes can also enhance plant stress resistance; for example, OsHAK1 [24] and HvHAK1 [48] are involved in drought stress, and AtHAK5 [49] and OsHAK16 [25] play a role in the salt stress response. In contrast, most members of group II are generally engaged in low-affinity K + absorption and complement K + channels [22,23]; the genes in group III have been shown to participate in the maintenance of K + /Na + homeostasis, such as AtKUP11, OsHAK11, and OsHAK12 [21]; and the genes in group IV are involved in Na + uptake, such as PpHAK13 from Physcomitrella patens [50], PhaHAK5 from Phragmites australis [51], and ZmHAK4 from maize [26]. Remarkably, more genes in groups II, III, and IV have been reported to be associated with various plant developmental processes. For instance, AtKUP4 (II) is crucial for both gravitropic responses and root hair formation by mediating auxin efflux [30]; AtKUP2, AtKUP6, AtKUP8 [4], and CnHAK1 [23] in group II are involved in K + -dependent cell expansion and mediate K + efflux in roots, which negatively regulate plant growth and cell size; AtKUP4 and OsHAK10 in group II play a role in seed maturation and reproductive processes, respectively [14,52]; AtKUP9 in group III can maintain root growth and meristem activity by regulating K + and auxin homeostasis under low-K + stress [27]; the cotton GhKT1 in group III plays a role in fiber elongation [53]; and LjKUP in group IV from Lotus japonicus is involved in nodulation development [54]. As mentioned above, group IV is absent in Brassicaceae [38]. On the basis of our findings, we suspected that this may be due to its functional redundancy with groups II and III. In summary, the five groups of the HAK family commonly play a role in plant K + uptake and transport, as well as in the stress response, growth, and development.
This study provides a useful basis for future studies on the functions of the BnaHAKs and can contribute to the long-term goal of improving the K + utilization efficiency of oil crops.

Identification and Phylogenetic Analysis of the HAK Family in B. napus
The AtKUPs were obtained from the TAIR database (http://www.arabidopsis.org). To identify the HAK-encoding genes in the B. napus genome, we performed a BLASTP [55] search in the GENOSCOPE database (http://www.genoscope.cns.fr/brassicanapus/), using the known AtKUPs as queries with a low-stringency criterion (cutoff p < 0.1). After deleting the redundant sequences, the remaining sequences were examined by the Pfam tool (http://pfam.xfam.org/search/sequence) to ensure that they contained the typical domains of the HAK family. We excluded the severely missing sequences from our dataset, and finally obtained 40 BnaHAKs. The DNA and coding sequences of the candidates were obtained from the GENOSCOPE database. Using the same method, we also identified BrHAKs and BoHAKs from the Phytozome v12.1 database (http://www.Phytozome.net) [56]. The biochemical properties of the candidates were predicted using the ExPASy tool (http://www.expasy.org/tools/) [57], and the subcellular localization was investigated using Plant-mPLoc (http://www.csbio.sjtu.edu.cn/ bioinf/plant-multi/).
To explore the evolutionary relationship of the HAK family in B. napus, B. oleracea, B. rapa, and Arabidopsis, we performed multiple sequence alignment of the obtained protein sequences of BnaHAKs, BrHAKs, BoHAKs, and AtKUPs using the MAFFT version 7 tool with default parameters (https://mafft.cbrc.jp/alignment/server/). Subsequently, a phylogenetic tree was built using MEGA v5.2 with the NJ method based on the multiple sequence alignment. The parameters used in the phylogenetic analyses were as follows: Poisson correction, bootstrap with 1000 replicates, and pairwise deletion. Finally, FigTree v1.4.0 (http://tree.bio.ed.ac.uk/software/figtree/) was applied to view and edit the tree file.

Gene Structure Analysis of BnaHAKs
The gene structures of BnaHAKs and AtKUPs were analyzed using Gene Structure Display Server (GSDS) 2.0 using with the DNA and coding sequences (http://gsds.cbi.pku.edu.cn/) [58]. The TM regions of BnaHAKs and AtKUPs were predicted by SMART (http://smart.embl-heidelberg.de/). By comparing the DNA and coding sequences of each BnaHAK using MEGA v5.2, the intron insertion sites in the corresponding protein sequences were manually located, and the intron insertion information for the AtKUPs, BrHAKs, and BoHAKs was acquired from Phytozome v12.1 (https: //phytozome.jgi.doe.gov/pz/portal.html).

Chromosomal Location and Collinearity Analysis of BnaHAKs
We acquired information on the chromosome locations of candidate BnaHAKs from the GENOSCOPE database. Mapchart v2.2 software was applied to draw the chromosome map of BnaHAKs. The cross-genome collinearity relationship of BnaHAKs, BrHAKs, BoHAKs, and AtKUPs was calculated and identified using the CoGe online tool (https://genomevolution.org/coge/). The duplication events of BnaHAKs were defined based on the collinearity relationship. The nucleotide substitution rate (Ka/Ks) of BnaHAKs was calculated using KaKs_Calculator2.0 software.

Structure Prediction Analysis of BnaHAK Proteins
To explore the tertiary structure features, the protein sequences of BnaHAKs were submitted to SWISS-MODEL, a protein-modeling server (https://swissmodel.expasy.org/interactive). To validate the secondary structural information, further analysis was performed by submitting the protein sequences of BnaHAKs to the Protter v1.0 online tool (http://wlab.ethz.ch/protter/) [60].

Gene Expression Analysis
We used an RNA-Seq dataset in the National Center of Biotechnology Information (NCBI) (BioProject: PRJNA358784) to detect the temporal and spatial expression patterns of BnaHAKs across all developmental stages of the B. napus cultivar 'Zhongshuang 11 (ZS11). Similarly, an RNA-Seq dataset of B. napus ZS11 seedling roots under five hormone inductions (IAA, GA3, 6-BA, ABA, and ACC) was obtained from NCBI (BioProject ID PRJNA608211) to explore the hormone-responsive expression patterns of BnaHAKs. The R package v3.5.3 [61] was used to analyze and draw a heatmap based on log2-transformed data. Genes with FPKM < 1 might be pseudogenes or only expressed under specific conditions; thus, they were excluded from the heatmap. The Pearson correlation coefficient was obtained by calculating the expression levels of homologous genes in different tissues/organs.

Plant Materials and Growth Condition
Seeds of ZS11 were obtained from the College of Agriculture and Biotechnology, Southwest University. The seeds were germinated in individual plastic pots filled with vermiculite, grown in an artificial climatic chamber at 25 • C with a 16:8 h photoperiod (day:night), and watered with 1/2-strength Hoagland solution every four days. Then, seedlings at the four-leaf stage were changed from soil culture to hydroponic culture with 1/2-strength Hoagland solution. The solutions were changed every three days. The seedlings at the five-leaf stage were used for the -K + treatment. In this treatment, three independent repeated trials were performed, each with three plants harvested at each sampling. The normal nutrient solution comprised 1.25 mM KNO 3 , 1.5 mM Ca(NO 3 ) 2 , 0.75 mM MgSO 4 , 0.5 mM KH 2 PO 4 , 75 µM FeEDTA, 50 µM H 3 BO 3 , 10 µM MnCl 2 , 2 µM ZnSO 4 , 1.5 µM CuSO 4 , and 0.075 µM (NH 4 ) 6 Mo 7 O 24 (control, CK), whereas for the -K + treatment, 1.25 mM KNO 3 and 0.5 mM KH 2 PO 4 were replaced by 0.5 mM phosphoric acid [31]. The pH was adjusted to 5.8 with Tris. The root tissues were harvested at 0.5, 6, 12, 24, 48, and 72 h after the treatments, immediately frozen in liquid nitrogen, and then stored at −80 • C for RNA isolation.

RT-qPCR Analysis of BnaHAKs under Low-K + Conditions
The EASYspin total RNA Extraction kit (Biomed, Beijing, China) was used to extract the total RNA from each sample. The concentration and quality of the total RNA were tested using gel electrophoresis and a NanoDrop 2000 spectrophotometer to confirm that the A260/280 ratio remained at 1.8-2.1, and that the A260/230 ratio exceeded 2.0. The RNA sample was treated with DNase I (Promega, Beijing, China), and was then used for cDNA synthesis by reverse transcription in a 20-µL reaction system according to the manufacturer's instructions of the M-MuLV RT kit (Takara Biotechnology, Beijing, China). The primers used in this experiment were designed using Primer Premier 5 software and are listed in Table S9. BnaActin7 (GenBank accession no. AF024716) and BnaUBI (GenBank accession no. NC027770) served as double reference genes. The SYBR-Green PrimeScript RT-PCR Kit (Takara Biotechnology, Beijing, China) was used for real-time PCR analysis using the CFX Connect™ Real-Time System (Bio-Rad, Chongqing, China), and every reaction system consisted of three technical replicates. The thermocycling parameters included initial denaturation at 95 • C for 5 min, followed by 45 cycles of denaturation at 95 • C for 15 s and annealing at 60 • C for 15 s (the annealing temperature of BnaHAK32 and BnaHAK19 was 55 • C). Finally, we obtained the data (mean ± standard deviation) of all three independent repeated trials and calculated the relative expression of BnaHAKs using the 2 (−∆∆Ct) method. Error bars represent standard errors from three independent repeated trials. Differences in expression levels in BnaHAKs according to K + treatments were assessed by one-way analysis of variance (* p < 0.05; ** p < 0.01) using Excel 2010.

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