Identification of NPF Family Genes in Brassica rapa Reveal Their Potential Functions in Pollen Development and Response to Low Nitrate Stress

Nitrate Transporter 1/Peptide Transporter Family (NPF) genes encode membrane transporters involved in the transport of diverse substrates. However, little is known about the diversity and functions of NPFs in Brassica rapa. In this study, 85 NPFs were identified in B. rapa (BrNPFs) which comprised eight subfamilies. Gene structure and conserved motif analysis suggested that BrNFPs were conserved throughout the genus. Stress and hormone-responsive cis-acting elements and transcription factor binding sites were identified in BrNPF promoters. Syntenic analysis suggested that tandem duplication contributed to the expansion of BrNPFs in B. rapa. Transcriptomic profiling analysis indicated that BrNPF2.6, BrNPF2.15, BrNPF7.6, and BrNPF8.9 were expressed in fertile floral buds, suggesting important roles in pollen development. Thirty-nine BrNPFs were responsive to low nitrate availability in shoots or roots. BrNPF2.10, BrNPF2.19, BrNPF2.3, BrNPF5.12, BrNPF5.16, BrNPF5.8, and BrNPF6.3 were only up-regulated in roots under low nitrate conditions, indicating that they play positive roles in nitrate absorption. Furthermore, many genes were identified in contrasting genotypes that responded to vernalization and clubroot disease. Our results increase understanding of BrNPFs as candidate genes for genetic improvement studies of B. rapa to promote low nitrate availability tolerance and for generating sterile male lines based on gene editing methods.

During plant development and responses to environmental stresses, nutrients and other substrates are transported according to altered metabolic pathways involved in the synthesis, storage, mobilization, and conversion, thereby requiring transporters with capacities to transport a wide diversity of chemical substrates [10,11]. NPFs are one of the largest transporter groups in plants and can transport a wide range of substrates. For example, AtNPF2.10/GTR1 (glucosinolate transporters 1) and AtNPF2.11/GTR2 are essential for the  (Table S1). The numbers of Br/Bol/BniNPF proteins were 1.6, 2.0, and 1.8 times that in Arabidopsis, respectively, owing to the gene expansion of Brassica during their evolution. The gene IDs, genome locations, coding sequence lengths, protein lengths, and other characteristics of the identified Br/Bol/BniNPF genes are listed in Table S1. The subcellular localization of these 292 proteins were predicted using the WoLF PSORT program, revealing that most of the NPF proteins were located on the plasma membrane (269), chloroplast (11), or in vacuoles (8) (Table S1). The lengths of BrNPFs ranged from 98 to 1070 amino acids, with an average length of 534 amino acids, and MWs ranging from 8.77 to 383.66 kDa, with an average weight of 62.08 kDa. Likewise, a wide range of PIs was also observed (4. 73-9.74) which may be due to the considerable range in protein lengths (Table S1).
The conserved motifs of the BrNPFs were analyzed with the MEME program, yielding the identification of 10 conserved motifs among 85 BrNPFs from B. rapa ( Figure 2 and Table S2). The 10 motifs of typical BrNPF proteins followed the order of Motif9-Motif4-Motif2-Motif10-Motif6-Motif5-Motif3-Motif7-Motif1-Motif8. Apart from motifs 7, 8, and 9, the other motifs contained the core PTR2 domain conserved sequence (Table S3). No significantly conserved motifs were identified for the 7, 8, and 9 domains based on BLAST searches with NCBI and Pfam.
In order to investigate the structural diversity of BrNPFs, exon-intron organizations were analyzed. The number of BrNPF introns ranged from 0 to 14. The most common exonintron organization comprised four exons separated by three introns, which were present in 41 BrNPFs ( Figure 2 and Table S2). Except for BrNPF2.8, BrNPF2.18, and BrNPF7.4, most BrNPFs contained more than one intron, indicating the possible existence of alternative splicing during expression. Eight types of gene structures were identified in the NPF5 subfamily, implying they exhibited diverse functions. BrNPF members from the same subgroups exhibited similar gene structures for the other subfamilies, indicating the potential for conserved functions ( Figure 2).

BrNPF Chromosomal Location and Gene Duplication Analysis
Eighty-five BrNPF genes were present on the ten chromosomes of Brassica rapa and were non-uniformly distributed ( Figure A1). Chromosome A09 harbored the largest number of BrNPFs (18 members), followed by chromosomes A07, A02, and A06, which contained 14, 13, and 11 BrNPFs, respectively. Chromosome A04 carried the smallest number of BrNPFs (two). BLAST and MCScanX analysis indicated that BrNPF gene duplication events were present in the B. rapa genome. Briefly, 22 tandem duplicated genes (25.9%) were identified that belonged to ten clusters (Figures 3 and A1). Among the tandem duplicated genes, two, two, and three clusters were located on chromosomes A02, A07, and A09, respectively. The other three clusters were located on chromosomes A03, A05, and A06, respectively (Figures 3 and A1). These results suggest that tandem duplication is related to NPF expansion in Brassica genomes.  (53). The neighbor-joining phylogenetic tree was generated in MEGA6 with full-length NPF protein sequences, and branch support was evaluated with 1000 bootstrap replicates.

BrNPF Chromosomal Location and Gene Duplication Analysis
Eighty-five BrNPF genes were present on the ten chromosomes of Brassica rapa and were non-uniformly distributed ( Figure A1). Chromosome A09 harbored the largest number of BrNPFs (18 members), followed by chromosomes A07, A02, and A06, which contained 14, 13, and 11 BrNPFs, respectively. Chromosome A04 carried the smallest number of BrNPFs (two). BLAST and MCScanX analysis indicated that BrNPF gene duplication events were present in the B. rapa genome. Briefly, 22 tandem duplicated genes (25.9%) and B. oleracea, collinear gene pairs were identified using the MCScanX software package. A total of 26, 18, and 13 gene pairs were identified between B. rapa and B. oleracea, B. rapa and B. nigra, and B. rapa and Arabidopsis, respectively ( Figure 3 and Table S3). All syntenic NPF genes in B. rapa were located on chromosomes A01, A02, A07, and A09 ( Figure 3). Six segmental duplication events were also identified within the B. oleracea genome ( Figure  3), indicating that the greater numbers of NPFs in the B. oleracea genome could be due to segmental gene duplication events in its genome. To evaluate the collinear relationships of all NFP genes in Arabidopsis, B. rapa, B. nigra, and B. oleracea, collinear gene pairs were identified using the MCScanX software package. A total of 26, 18, and 13 gene pairs were identified between B. rapa and B. oleracea, B. rapa and B. nigra, and B. rapa and Arabidopsis, respectively ( Figure 3 and Table S3). All syntenic NPF genes in B. rapa were located on chromosomes A01, A02, A07, and A09 ( Figure 3). Six segmental duplication events were also identified within the B. oleracea genome (Figure 3), indicating that the greater numbers of NPFs in the B. oleracea genome could be due to segmental gene duplication events in its genome.

Cis-Elements in Promoters of BrNPFs
In order to identify the cis-regulatory elements in BrNPF promoters, cis-elements were analyzed using the PlantCARE platform (https://bioinformatics.psb.ugent.be/webtools/ plantcare/html/ accessed on 16 August 2022). A total of 904 phytohormone-responsive elements, 1605 environmental-responsive elements, 205 plant growth, and developmentrelated elements, and 979 transcriptional factor binding sites were predicted within the promoters of the 85 BrNPFs (Table S4). Among these, light-responsive, MYB transcriptionbinding, MYC transcription-binding, and ABA-responsive cis elements were the four most prevalent (Table S4). Thus, most BrNPFs respond to diverse environmental stresses and are regulated by various transcriptional factors (TFs).

Tissue Expression of BrNPFs Reveals Their Potential Functions during Pollen Development
In order to investigate the expression of BrNPFs, their expression patterns were compared using RNA-sequencing data from 59 different organ or tissues samples, including callus, roots, stems, stem leaves, flowers, siliques, head leaves (24 samples), developmental stages of floral buds (10 samples), pistils (four samples), unfertilized ovules, embryos (seven samples), and seed coats (seven samples) (Table S5). BrNPFs were differentially expressed among groups of the 59 tissue samples ( Figure 4). The BrNPF1.1, BrNPF1.2, BrNPF1.3, and BrNPF6.5 genes were generally predominantly expressed in all tissues except embryos and seed coats. BrNPF2.5 and BrNPF2.21 were mostly expressed in stem leaves and opened flowers. BrNPF7.3 and BrNPF7.5 were predominantly expressed in roots, while BrNFP7.1 was predominantly expressed in callus tissues. BrNPF2.5 and BrNPF2.21 were mostly expressed in all developmental floral buds and late development seed coats. Further, BrNPF2.6, BrNPF2.15, BrNPF7.6, and BrNPF8.9 were predominantly expressed in fertile buds from the uninucleate to binucleate pollen stages. The primary difference between fertile and sterile floral buds is the presence of pollen grains [40]. Accordingly, BrNPF2.6, BrNPF2.15, BrNPF7.6, and BrNPF8.9 were predominantly expressed per normal pollen development in B. rapa. Semi-quantitative RT-PCR was conducted to confirm the expression patterns of these four genes based on the different developmental stages of floral buds from male genetic sterility (GMS) lines ( Figure 5A,B). The expression patterns of these four genes were similar between RNA-Seq and RT-PCR datasets. Briefly, BrNPF2.6, BrNFP7.6, and BrNPF8.9 were specifically expressed in the floral buds from F2 (floral buds containing tread stage pollen) to F3 (floral buds after the tetrad stage, but before containing mature pollen) stages ( Figure 5B). In addition, BrNPF2.15 was only specifically expressed in F3 floral buds ( Figure 5B).

Expression of BrNPFs during B. rapa Growth during Vernalization and P. brassicae
NPFs have been previously suggested to be involved in plant responses during vernalization [21]. In order to identify BrNPFs that might be responsive to vernalization in B. rapa, RNA-Seq data for BrNPFs during vernalization were re-calculated, as previously described [51]. Briefly, BrNPF2. 5 Figure 6A).

Expression of BrNPFs during B. rapa Growth during Vernalization and P. brassicae
NPFs have been previously suggested to be involved in plant responses during vernalization [21]. In order to identify BrNPFs that might be responsive to vernalization in B. rapa, RNA-Seq data for BrNPFs during vernalization were re-calculated, as previously described [51]. Briefly, BrNPF2.5, BrNPF2.21, BrNPF3.3, BrNPF6.1, BrNPF6.3, BrNPF7.1, and BrNPF7.2 were up-regulated during vernalization in both genotypes, while BrNPF7.3 and BrNPF7.5 were only up-regulated in the late bolting type (Table S7 and Figure 6A). BrNPF2.5 and BrNFP2.21 are homologs of AtNPF2.13 that are responsible for source-to-sink remobilization of nitrate [15], indicating that they might play important roles during nitrate remobilization during vernalization. BrNPF3.3 is a homolog of AtNPF3.1 that enables the transport of GA [28], indicating that BrNPF3.3 might play important role in GA transport during the developmental phase transition.
The development of clubroot disease can be influenced by nitrogen fertilization [52], and NPFs are one of the main transporters of nitrogen, indicating that NPFs might be involved in a response mechanism to P. brassicae infection. RNA-seq data suggested that many NPF members were expressed after inoculation with P. brassicae (Table S7 and Figure 6B). Briefly, BrNPF2.23 was only responsive to P. brassicae in the susceptible line, while BrNPF5.3 was induced in the susceptible line after infection by P. brassicae ( Figure 6B). BrNPF2.21 and BrNPF7.4 expressions were both up-regulated in the resistant line compared to the susceptible line ( Figure 6B). BrNPF2.24 exhibited down regulation in both the susceptible and resistant lines. Thus, NPF genes may participate in clubroot disease responses via the transport of their specific substrates.

Expression of BrNPF Responses to Low Nitrate Conditions
Most NPF genes in Arabidopsis are related to nitrate transport, indicating that nitrate is the primary substrate of NPF proteins [53]. To assess the potential functions of BrNPFs in nitrate uptake and use, B. rapa seedlings (accession Chiifu-401-42) were hydroponically cultured in Hoagland's nutrient solution [21] and treated with low nitrate conditions. After treatment, plant heights, leaf areas, fresh weights, and nitrogen contents were significantly lower, while root lengths increased compared to normal growth conditions ( Figure 7A-F). The root and shoot components of seedlings were then separately sampled for RNA-seq analysis. Considering the criteria of TPM ≥ 1 and fold change values > 2.0, a total of 39 BrNPFs were identified as responsive to low nitrate conditions either in the shoots or roots (Table S8 and Figure 7G). Among these, ten BrNPFs (BrNPF2. 12 Figure 7G). Three BrNPFs (BrNPF6.6, BrNPF6.7, and BrNPF7.3) were up-regulated in shoots, but downregulated in roots, indicating that they may exhibit contrasting roles in roots and shoots during low nitrate conditions.

Identification and Analysis of BrNPFs
NFPs are one of the largest groups of transporter family genes in plants. Genome-wide identification of NFP family genes has been conducted in many plants based on sequence and motif conservation, including in Arabidopsis [15], rice [54], apple [17], sugarcane [18], Brassica napus [19][20][21], hexaploid wheat (Triticum aestivum L.) [22], spinach (Spinacia oleracea L.) [23] and tea plants (Camellia sinensis) [24]. In this study, 85 NFPs were identified in B. rapa, in addition to 110 members in B. oleracea and 97 in B. nigra (Table S1). Ninety-five NPFs were previously identified in B. rapa ("Chiifu-401", version 1.5) [1]. The differences in identification could be explained by the use of a higher quality genome version ("Chiifu-401", version 3.0) [55], leading to greater accuracy in BrNPF identification. Similar results were also observed for B. napus, in which 169, 193, and 199 NPFs were identified based on different versions of B. napus genomic data [19][20][21]. Brassica species have undergone an extra genomic duplication event compared to Arabidopsis [56]. Thus, one Arabidopsis gene should theoretically have one to three orthologs in Brassica genomes. However, the expansion of NPFs in B. rapa, B. oleracea, and B. nigra represent 1.6-, 2.0, and 1.8-fold increases relative to the Arabidopsis genome, respectively. Thus, duplicated genes may have been lost during Brassica evolution. Consistently, AtNPF2.1, AtNPF2.2, and AtNPF2.5 did not have any orthologs in the three species of Brassica analyzed here, while six, six, and four orthologs of AtNPF5.2 were identified in the B. rapa, B. oleracea, and B. nigra genomes, respectively ( Figure 1 and Table S1).
Segmental duplication and tandem duplication are two major mechanisms of gene family duplication in plants [57]. In this study, 22 tandem duplicated genes (25.9%) were identified in the B. rapa genome, while no segmental duplication events were identified ( Figure S1). Thus, tandem duplication may be a primary driving force in the expansion of BrNPFs. In contrast, six segmental duplication events were observed within the genome of B. oleracea, which may explain the greater number of NPFs in the B. oleracea genome compared to the B. rapa and B. nigra genomes (Figure 3).

BrNPF Functions
NPF genes represent one of the largest transporter gene families in plants and participate in the transport of diverse substrates across membranes over short or long distances [11,19]. Cis-elements play key roles in controlling gene expression during plant development and are responsive to stress [58]. Numerous cis-elements have been previously observed in the promoter regions of NPFs from B. napus, tea plant, apple, and spinach plants, in addition to others [14,19,20,23,52]. In the present study, many phytohormone-responsive, environmental-responsive, plant growth and development-related, and transcriptional factor binding cis-elements were identified in the promoters of the BrNPFs (Table S4). The existence of abundant and diverse cis-elements in NPF gene promoters could be related to their multiple functions during plant development and responses to environmental stresses.
Gene expression patterns can provide clues for predicting gene functions. Consequently, expression profiles for BrNPFs from 59 diverse tissues were analyzed, in addition to expressional profiles responsive to vernalization and P. brassicae infection between contrasting genotypes and under low nitrate availability stress (Figures 4, 6 and 7). Some BrNPFs exhibited tissue-specific expression, while others exhibited differential expression due to vernalization, P. brassicae infection, and low nitrate availability. For example, the homolog of AtNPF3.1, BrNPF3.3, exhibited induction by vernalization ( Figure 6A). Nitrate was previously reported to delay flowering time via the GA signaling pathway [59]. Further, AtNPF3.1 in Arabidopsis has been reported to be involved in GA transport [60], with GA increasing during B. rapa vernalization [61]. The up-regulation of BrNPF3.3 due to vernalization suggests that it might contribute to GA presence during vernalization. Further, nitrogen fertilization has been shown to affect the susceptibility of B. napus [52]. Homologs of AtNPF2.13 and BrNPF2.21 exhibited up-regulation in the resistant genotype ( Figure 6B). AtNPF2.13 is involved in the remobilization of nitrate from sources to sinks [62], indicating that BrNPF2.21 might respond to clubroot disease through the redistribution of nitrate. These results collectively provide new insights for the future functional prediction and characterizations of BrNPFs.

BrNPFs and Pollen Development
AtNPF2.8 was previously suggested by co-expression analysis to be involved in the transport of flavonol-3-O-sophoroside from tapetum cells to pollen walls in Arabidopsis [13]. Two orthologs of AtNPF2.8 were identified in B. rapa, including BrNFP2.6 and BrNPF2.15 ( Figure 1 and Table S1). Both orthologs exhibited expression only in fertile floral buds of MS lines ( Figure 5B). Further, BrNPF2.6 was highly expressed in floral buds containing pollen grains from tetrads prior to the mature stage (from the F2 to F3 stages). In addition, BrNPF2.15 only exhibited expression in F3-stage floral buds that contain pollen grains after the tetrad to before the mature stages ( Figure 5B). Additionally, gene co-expression profiles and associated GO enrichment biological processes differed between BrNFP2.6 and BrNPF2.15. For example, pollen development was represented by genes co-expressed with BrNFP2.6, while pollen tube development was suggested by genes co-expressed with BrNPF2.15 ( Figure 5C,D). Thus, the function of AtNPF2.8 might have expanded or diversified over Brassica evolution. Flavonol diglucosides are essential for maintaining pollen fertility and increasing pollen tolerance to environmental stresses [63]. Thus, expanded AtNPF2.8 genes suggest that the regulatory network related to flavonol diglucoside metabolism during pollen development might be more complex in B. rapa compared to Arabidopsis.
BrNFP7.6 and BrNFP8.9 were also specifically expressed in fertile floral buds (Figures 4 and 5B) and were the only identified orthologs of AtNPF7.1 and AtNPF8.2, respectively ( Figure 1). AtNPF8.2 is also referred to as AtPRT5 and facilitates peptide transport into germinating pollen and possibly into maturing pollen, ovules, and seeds [36]. Pollen development and pollen tube development were suggested by genes co-expressed with BrNFP8.9 ( Figure 5C-E). Based on the phylogenetic and Semi-quantitative RT-PCR analysis results, the functions of AtNPF8.2 and BrNFP8.9 may be highly conserved during pollen development and pollen tube growth. The processes of pollination, pollen tube growth, lipid oxidation, transport, and response to nutrient levels were also represented by genes co-expressed with BrNPF7.6, indicating that they might also play important roles in pollen development. Taken together, the fertile floral buds exhibited genes that were specifically expressed in these tissues, and these may provide new breeding targets for creating sterile male lines in B. rapa.

Responses of BrNPFs during Low Nitrate Stress
Nitrate was previously identified as the main substrate of NPF genes [22,53]. Here 45.9% of BrNPFs (39 of 85) exhibited differential expression due to low nitrate availability stress ( Figure 7G). Among these, ten were induced in both shoots and roots, eleven were specifically induced in shoots, and seven were only up-regulated in roots, indicating that they might play positive roles in nitrate absorption, uptake, homeostasis, and redistribution under low nitrate availability conditions (Table S8 and Figure 7G). AtNPF2.3 in Arabidopsis functions as a root stele transporter and contributes to nitrate translocation to shoots during salt stress [12]. In this study, a homolog of AtNPF2.3, BrNPF2.3, exhibited specific induction in roots, indicating that it might serve a similar function (Figures 1 and 7G). AtNPF6.2 plays key role in regulating leaf nitrate homeostasis [32]. The expression of its homolog, BrNPF6.5, was repressed in both shoots and roots, indicating that low nitrate availability may lead to decreased NPF transporter activity (Figures 1 and 7G). AtNPF2.13 is reportedly involved in the transport of nitrate and GA [15]. A homolog of AtNPF2.13, BrNPF2.21, was up-regulated in both shoots and roots, indicating that BrNPF2.21 might function in response to low nitrate availability by coupling to hormone signaling. AtNPF6.3 can repress lateral root growth under low nitrate availability by promoting basipetal auxin transport out of roots [3]. In this study, a homolog of AtNPF6.3, BrNPF6.6, was down-regulated in roots but up-regulated in shoots (Figures 1 and 7G). Taken together, the expression of BrNPFs under low nitrate availability conditions suggests that there is a crosstalk between low nitrate stress responses and phytohormone signaling pathways, consistent with results from previous studies [59].

Plant Growth and Low Nitrate Treatments
Uniform B. rapa seeds (accession Chiifu-401-42) were germinated in Petri dishes at 23 ± 1 • C in the dark for two days, followed by hydroponic cultivation of germinated seeds in Hoagland's nutrient solution for one week [64]. To establish low nitrate treatments, KNO 3 and Ca(NO 3 ) 2 in Hoagland's solutions were replaced with KCl and CaCl 2 , respectively. The final concentration of NO 3 − in the treatments was 0.1 mM. During cultivation, growth conditions within growth chambers were set as previously described [65]. After treatments, the shoots and roots in the low nitrogen and control treatments were individually harvested and immediately frozen in liquid nitrogen, followed by storage at −80 • C.
In order to collect materials from male sterile (MS) lines, seeds of MS lines from our previous study were germinated in Petri dishes at 23 ± 1 • C in the dark [40]. Vernalization was then induced with germinated seeds at 4 • C in the dark for 30 days. After vernalization, seeds were sown into pots (15 × 15 × 18 cm) containing potting soil and transferred to a greenhouse, followed by growth at 23 ± 1 • C with a light intensity of 6000-7000 Lux under a long day photoperiod (light/dark, 16 h/8 h). After flowering, floral buds were collected from MS line plants using three biological replicates and with previously reported criteria [40]. Root and shoot tissues were collected from three-week-old seedlings without vernalization. Stem and leaf tissues were sampled from plants one week after bolting. The siliques were collected two weeks after pollination. After sampling, tissues were immediately frozen in liquid nitrogen and stored at −80 • C until further analysis.

Identification of BrNPFs in Three Prototypical Diploid Species of Brassica
To identify NPF proteins from three prototypical diploid species of Brassica, all putative protein sequences encoded by the B. rapa ("Chiifu-401", version 3.0) [55], B. oleracea ("JZS," version 2.0) [66] and B. nigra ("Ni100", version 2.0) [67] genomes were downloaded from the Brassicaceae Database (BRAD; www.brassicadb.cn accessed on 20 April 2022). The previously identified 53 AtNPFs were used as queries to search against Brassica protein sequences (using an E-value < 10 −5 and identity > 20%). Search hits without a protondependent transport 2 (PTR2) domain (PF00854) were excluded based on HMM analysis with an E-value cut-off of 10 −5 . In order to identify all potential NPFs among the three Brassica, the search baits were used as BLAST queries for searching against the Phytozome 13 and NCBI databases with an E-value < 10 −5 and identity > 20%. No additional predicted NPFs were identified at this stage. All NPFs were identified according to previously reported rules [1].

Phylogenetic and Bioinformatic Analysis of BrNPFs
Phylogenetic analysis was conducted with the NPF protein sequences from Arabidopsis, B. oleracea, B. nigra, and B. rapa after alignment with the MUSCLE program, implemented in MEGA6 with default parameters [68]. An unrooted phylogenetic reconstruction was then constructed using MEGA6 with neighbor-joining methods and analysis parameters, including pairwise taxa deletion, 1000 bootstrap replicates, and the use of the Jones Taylor Thornton (JTT) amino acid substitution model [69]. The chromosomal positions of each NPF from the three genomes were identified among those from the BRAD (www.brassicadb.cn accessed on 16 August 2022) and visualized with a custom Python script. The isoelectric point (PI) and molecular weights (MWs) of the NPFs were analyzed using the ProtParam tool (Expasy, the Swiss Bioinformatics Resource Portal, https://web.expasy.org/protparam/ accessed on 16 August 2022) [70]. Subcellular localization predictions of NPFs were conducted using the WoLF PSORT software package (https://wolfpsort.hgc.jp/ accessed on 16 August 2022) with default settings. Conserved motifs in the BrNPFs were identified using the MEME software program (Suite 5.1.1, http://meme-suite.org/ accessed on 16 August 2022) [71]. BrNPFs gene structures were drawn using the Gene Structure Display Server (GSDS; version 2.0, http://gsds.cbi.pku.edu.cn/ accessed on 16 August 2022) [72].

RNA Extraction, Leaf Area, and Nitrate Content
Total RNA was isolated from 100 mg of homogenized leaves using the RNAiso Plus Reagent (Takara Biomedical Technology Co., Ltd., Beijing, China) according to the manufacturer's instructions. The outermost leaves of Chiifu seedlings were dissected for leaf area determination using a Yaxin-1241 leaf meter (Beijing Yaxinliyi, Beijing, China), following the manufacturer's instructions. The nitrogen concentrations of oven-dried shoots were measured using the Kjeldahl method with a JK9830 Kjeldahl Auto Analyzer (ELITE-Lab Instrument Co., Ltd., Jinan, China) [73] and are expressed as concentrations of per hundred dry matter (g/100 g).

RNA-Sequencing and Assembly
RNA samples from low nitrate treatments were sent to Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China) for RNA-Seq analysis. The RNA libraries were constructed and sequenced on the Illumina platform. Sequencing and analyses were conducted following standard protocols at Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China). Filtered clean reads were then mapped to the reference genome ("Chiifu-401", version 3.0) using the HISAT2 software program [74], and transcripts per million (TPM) values were calculated using DESeq [75]. Raw sequencing data were deposited in the China National Center for Bioinformation (CNCB) under the project ID PRJCA012597.

Expression of BrNPFs within RNA-Seq Data
To analyze the expression of BrNPFs among various tissues, publicly available RNAseq data from 59 different organs or tissues of B. rapa were retrieved from the NCBI database via the Bioproject accessions PRJNA185152, PRJNA778186, PRJNA641876, and PRJNA473318 [76][77][78][79]. Data were included from tissues comprising callus, stems, stem leaves, opened flowers, siliques, different developmental stages of heading leaves, floral buds, pistils, embryos, and seed coats. Gene expression levels were re-calculated using the transcripts per million (TPM) metric. Heatmaps for expression profiles of BrNPFs were generated using the TBtools software program (version 1.0987663) [80].
The expression of BrNPFs after vernalization was recalculated for a previous study (BioProject PRJNA615255) using the TPM metric [51]. To achieve vernalization, two inbred B. rapa accessions, including a late bolting type (JWW) and an early bolting type (XBJ), were investigated. Prior to vernalization, both inbred lines were grown at 25 ± 2 • C for 32 days under natural light conditions. Both inbred lines were then transferred to a growth chamber at 4 • C with 150 µmol m −2 s −1 light intensity under long daylight conditions (16/8 h, day/night) for vernalization, followed by a collection of the third fully expanded leaves from the center for subsequent analyses. JWW leaves were collected after 0, 25, 30, 35, 45, and 50 days following treatment. XBJ leaves were collected 0, 10, 15, 25, 40, and 50 days after treatment.
The TPM values of BrNPF genes were recalculated after infection with Plasmodiophora brassicae, as described in previous studies (Bio-Project PRJNA743585) [81]. In order to initiate P. brassicae infection, 20-day-old healthy plants of resistant (BrT24) and susceptible (Y510-9) B. rapa genotypes were inoculated with 20 mL of a P. brassicae (race 4) solution.
For the control group, 20 mL of sterile water was used for inoculation. The root samples for each genotype were then collected at 0, 3, 9, and 20 d after inoculation, based on the four-time points of disease development [81].

Conclusions
Here, a total of 85, 110, and 97 NFP proteins were identified in the genomes of B. rapa, B. oleracea, and B. nigra, respectively. The gene structures, chromosomal locations, conserved motifs, cis-elements, evolutionary relationships, gene duplications, and expression patterns of the BrNFPs were systematically analyzed. These results provide new targets for future studies to elucidate the molecular mechanisms underlying BrNPF functions in pollen development, nitrate utilization, responses to vernalization, and P. brassicae infection response in B. rapa, especially for BrNPF2.6, BrNPF2.15, BrNPF7.6, and BrNPF8. 9 showing potential for generating sterile male lines based on gene editing methods in B. rapa and, possibly, other crops.