Genome-Wide Analysis of Purple Acid Phosphatase Genes in Brassica rapa and Their Association with Pollen Development and Phosphorus Deprivation Stress

Genes Their Abstract: PAPs (purple acid phosphatases) belong to the metallo-phosphoesterase superfamily and play important roles in developmental processes, phosphorus foraging, and recycling. However, the speciﬁc functions of BrPAPs in Brassica rapa are poorly understood. In this study, 39 BrPAPs were identiﬁed and divided into three major clades and nine subgroups. In 8 of the 39 BrPAPs, some invariant amino acid residues were lost or shifted. Based on an expression proﬁling analysis, BrPAP11 , 14 , 20 , 24 , 29 , and 34 were speciﬁcally expressed in fertile ﬂoral buds, indicating their critical roles during pollen development. A total of 21 BrPAPs responded to Pi deprivation in either shoots or roots. Of these, BrPAP4 , 5 , 19 , and 21 were upregulated in roots under Pi depravation conditions, while BrPAP12 was upregulated in the roots in normal conditions. BrPAP28 was upregulated in shoots under Pi depravation conditions, indicating its function shifted compared with its Arabidopsis homolog, AtPAP26 . The present work contributes to further investigation of BrPAPs as candidate genes for genetic improvement studies of low phosphorus tolerance as well as for creating male sterile lines based on gene editing methods in Brassica rapa .


Introduction
Purple acid phosphatases (PAPs) belong to the metallo-phosphoesterase superfamily and are a type of acid phosphatase (APase) comprised of a binuclear metal center binding Fe (III)-M (II) complex (where M = Fe, Zn, or Mn) at the active site [1,2]. Due to the existence of a tyrosine residue ligated to a ferric iron, this group of acid phosphatases has a purple color [3]. They have been isolated from most eukaryotic organisms and some bacteria [3,4]. APases can hydrolyze a wide range of anhydrides and phosphate esters [5]. Therefore, improvement of APase activity is an important avenue for improving the efficiency of phosphorus utilization of plants in a low Pi (inorganic phosphate) environment. This enables plants to utilize extracellular and intracellular organic phosphorus (Po) [5][6][7].

Plant Growth and Pi Deprivation Treatment
Brassica rapa L. ssp. pekinensis 'Chiifu' seeds were germinated in Petri dishes at 23 ± 1 • C in darkness. The germinated seeds were sown into plastic pots (7 cm × 7 cm) filled with vermiculite and grown for 20 d at 23 ± 2 • C with a light intensity of 6000-7000 Lux on a 16:8 h (L:D) photoperiod. During this growth phase, the plants were watered every two days with Hoagland solution (in mmol/L: KNO 3 6, Ca(NO) 3 [29]. The pH of the solution was adjusted to 6.0 using HCl and NaOH. For Pi deprivation treatment, the KH 2 PO 4 in the Hoagland solution was replaced with KCl. After 20 d of cultivation, 15 to 20 plants were selected randomly. Five plants were used for root length, root weight, fresh weight, leaf area, and phosphorus content measurements. The roots and shoots of the remaining plants were separated and sampled. After sampling, the roots and shoots were frozen immediately in liquid nitrogen and stored at −80 • C until use. All Pi deprivation treatments were carried out with three independent biological replicates.
For the fertile and sterile plants of B. rapa, seeds were germinated in Petri dishes at 23 ± 1 • C in darkness, then the germinated seeds were transferred to a 4 • C growth chamber in darkness for 30 d to induce vernalization. After vernalization, the seeds were sown into pots (15 × 15 × 18 cm) containing potting soil and grown in a greenhouse at 23 ± 1 • C with a light intensity of 6000-7000 Lux and a 16:8 (L:D) photoperiod. The floral buds were then collected from five plants using three biological replicates, as previously described [30]. Root and shoot tissues were collected from three-week-old seedlings without vernalization. Stems and stem leaf tissues were sampled from these plants seven days after bolting. All tissues were stored at −80 • C until use.

Leaf Area and Phosphorus Content
For leaf area determination, the outermost leaf of Chiifu seedlings was dissected, and the leaf area was determined with a Yaxin-1241 leaf meter (Beijing Yaxinliyi, China) following manufacturer's instructions. Total phosphorus concentration was determined as described previously [31]. In brief, 0.2 g of leaves was ground and extracted using a laboratory-scale high-performance microwave digestion system (Milestone Ethos UP, Italy, Sorisole) with 8 mL HNO 3 -H 2 O 2 (3:1) reagent. Afterward, the mixtures were reacted with molybdenum blue reagent and measured 20 min later at 700 nm by using a spectrophotometer (UV-1200, Instrument Co., Ltd., Shanghai, China).

Identification of BrPAP Genes in Brassica rapa
To identify the PAP genes from B. rapa, all putative protein sequences of B. rapa (version 3.0) were downloaded from BRAD (http://brassicadb.agridata.cn/brad/ accessed on 28 January 2021) [32]. All of the protein sequences of the 29 AtPAPs (A. thaliana purple acid phosphatase) were retrieved from TAIR (https://www.arabidopsis.org/index.jsp; Araport11 accessed on 28 January 2021) based on a previous publication [2] and used as queries to search against all putative protein sequences of B. rapa using BLASTP. A total of 43 putative BrPAP proteins were obtained based on two criteria: an E-value below 10 −5 and sequence identity above 20%. Next, the amino acid sequences of the 43 putative BrPAP proteins were analyzed to determine the presence of conserved sequence motifs (DXG/GDXXY/GNH(E/D)/VX2H/GHXH) defined previously for BrPAPs [2], and BraA02g012470.3C, BraA04g002700.3C, BraA06g026740.3C, and BraA10g003360.3C were removed. As a result, 39 proteins containing these conserved sequence motifs were identified. To find all potential PAPs in B. rapa, each of the 39 BrPAPs were used as queries for BLAST searching at NCBI and Phytozome 13 (https://phytozome-next.jgi.doe.gov/phytomine/ begin.do accessed on 28 January 2021). However, no additional predicted BrPAPs were found. All BrPAPs were named BrPAP1 to BrPAP39 based on their genomic locations.
For co-expression analysis, RNA-seq data from male sterility lines of B. rapa were reassembled [40,41,43,44], and the genes with no differences between fertile and sterile floral buds were removed. Pearson's correlation coefficient (PPC) values were calculated in Excel using BrPAP candidates as queries. The threshold values for PPCs were set between 0.60 and −0.60. GO (Gene Ontology) enrichment analysis was carried out using agriGO (http://systemsbiology.cau.edu.cn/agriGOv2/ accessed on 28 January 2021) [49] and the clusterProfiler package in R [50].

RNA Extraction, Semi-RT-PCR, and qRT-PCR
Total RNA was isolated using the RNAiso Plus Reagent (Takara Biomedical Technology Co., Ltd., Beijing, China), based on manufacturer protocols. The PrimeScript™ RT reagent Kit with gDNA Eraser (Takara Biomedical Technology Co., Ltd., Beijing, China) was used to synthesize 1st strand cDNA with 1 µg of total RNA. After synthesis, this cDNA was diluted to 10 ng/µL for PCR analysis. Semi-RT-PCR was carried out using 20 ng/µL cDNA with the following program: denaturation at 94 • C for 5 min, followed by 28 cycles of 94 • C for 30 s, 55 • C for 30 s, and 72 • C for 60 s. For qRT-PCR analysis, 10 ng/µL of cDNA was used and performed with a 30 s pre-denaturation at 95 • C, followed by 40 cycles of 95 • C for 5 s, 55 • C for 30 s, and 72 • C for 60 s. All of the primers used are listed in Table S1. The B. rapa actin gene was used as an internal control. The products of semi-RT-PCR were separated on 1.5% agarose gels and stained with GeneGreen (Tiangen Technology Co., Ltd., Beijing, China). The qRT-PCR results were analyzed using the 2 −∆Ct method with three biological replicates. After removing the genes with Ct value above 35 and less than two-fold change, the results are represented by a heatmap using R scripts with row normalization.

Brassica rapa PAP Genes
After BLASTP search and conserved motif determination, 39 PAP genes were identified from the B. rapa genome and designated as BrPAP1 to BrPAP39 according to their positions on chromosomes ( Figure 1). Among them, 31 BrPAP proteins contained the seven conserved metal-binding residues (D, D, Y, N, H, H, and H) in five consensus blocks (Table 1). Another three BrPAPs (BrPAP14, BrPAP17, and BrPAP21) lacked the first conserved block. The fourth block was lost in BrPAP9. For BrPAP7, one conserved residue in the fourth block was changed from H to S and the fifth block was missing. In the proteins from BrPAP22, BrPAP32, and BrPAP39, one conserved residue in the second block changed from Y to F, which has been reported in Arabidopsis [2], soybean [51], rice [11], maize [12], chickpea [13], and tea [15]. sizes of the genomic DNA sequences for BrPAPs varied from 1196 bp (BrPAP16) to 6766 bp (BrPAP23), with an average length of 2443 bp. Molecular weight of identified BrPAPs ranged from 34.8 kDa (BrPAP21) to 75.79 kDa (BrPAP20). Based on predictions from SignalP, signal peptides were found in 34 BrPAPs except for BrPAP2, BrPAP9, BrPAP21, BrPAP25, and BrPAP39. N-glycosylation site prediction analysis indicated that BrPAPs exhibited different numbers of glycosylation sites, from 0 to 6 ( Table 2).

Gene ID
Gene Name

PAP Conserved Motifs GDXG GDXXY GNH(D/E) VXXH GHXH
General information for BrPAPs, including gene size, molecular weight, Arabidopsis homologs, and potential N-linked glycosylation sites, were also analyzed ( Table 2). The sizes of the genomic DNA sequences for BrPAPs varied from 1196 bp (BrPAP16) to 6766 bp (BrPAP23), with an average length of 2443 bp. Molecular weight of identified BrPAPs ranged from 34.8 kDa (BrPAP21) to 75.79 kDa (BrPAP20). Based on predictions from SignalP, signal peptides were found in 34 BrPAPs except for BrPAP2, BrPAP9, BrPAP21, BrPAP25, and BrPAP39. N-glycosylation site prediction analysis indicated that BrPAPs exhibited different numbers of glycosylation sites, from 0 to 6 ( Table 2).

Phylogenetic and Gene Structure Analyses of BrPAPs
A phylogenetic analysis was carried out using amino acid sequences from BrPAPs and AtPAPs with the neighbor-joining method. All PAP proteins were classified into three major clades (I to III), and a further classification of these three major clades produced nine subgroups ( Figure 2). This analysis was similar to the findings of a previous study in Arabidopsis [2], with one exception. AtPAP10, 12, and 26 were classified into the same group in Arabidopsis [2], whereas AtPAP26 was separated from AtPAP10 and 12 and classified into a new subgroup (I a-3) in this study. Previously, a single gene was identified as a potential ortholog of AtPAP26 in maize and rice, indicating its conserved functions in the Pi deprivation response, which was also observed in B. rapa and was grouped into I a-3 ( Figure 2). Based on our phylogenetic analysis, one AtPAP was related to one to three BrPAP homologs, which may have been due to whole-genome polyploidization and a gene loss event of B. rapa after separation from Arabidopsis, although this was not the case for AtPAP10 (Table 2 and Figure 2). For AtPAP10, five homologs (BrPAP2, BrPAP8, BrPAP9, BrPAP10, and BrPAP33) were identified in B. rapa, which indicated the expansion of AtPAP10 in B. rapa over evolutionary time.
During the evolution of a multigene family, gene structure usually diversifies. We therefore wanted to study BrPAPs related to their functional diversification and evolution. The intron size and number of BrPAPs were highly variable, showing 11 distinct exonintron organization patterns. Seven exons separated by six introns was the most common organizational pattern and was present in 10 BrPAPs (Table 2 and Figure 3). Except for BrPAP6, most BrPAPs contained more than one intron, indicating the possible existence of alternative splicing during gene expression. For Ia-1, IIb, and IIIb, genes from the same subgroups showed similar gene structures, indicating that they might have conserved functions. The diversity of gene structures of genes from IIIa implied their diverse functions ( Figure 3). classified into a new subgroup (I a-3) in this study. Previously, a single gene was identified as a potential ortholog of AtPAP26 in maize and rice, indicating its conserved functions in the Pi deprivation response, which was also observed in B. rapa and was grouped into I a-3 ( Figure 2). Based on our phylogenetic analysis, one AtPAP was related to one to three BrPAP homologs, which may have been due to whole-genome polyploidization and a gene loss event of B. rapa after separation from Arabidopsis, although this was not the case for AtPAP10 (Table 2 and Figure 2). For AtPAP10, five homologs (BrPAP2, BrPAP8, BrPAP9, BrPAP10, and BrPAP33) were identified in B. rapa, which indicated the expansion of AtPAP10 in B. rapa over evolutionary time.  During the evolution of a multigene family, gene structure usually diversifies. We therefore wanted to study BrPAPs related to their functional diversification and evolution. The intron size and number of BrPAPs were highly variable, showing 11 distinct exonintron organization patterns. Seven exons separated by six introns was the most common organizational pattern and was present in 10 BrPAPs (Table 2 and

Potential Functions of BrPAPs during Pollen Development
To study the functional divergence of BrPAPs, the spatial and temporal expression patterns of BrPAPs in B. rapa were analyzed. RNA-seq data derived from shoots (vernalized and non-vernalized), leaves, stems, roots, flowers, differential developmental stages of seeds, and male sterility lines were explored [22,39,[41][42][43][44]. Based on this RNA-

Potential Functions of BrPAPs during Pollen Development
To study the functional divergence of BrPAPs, the spatial and temporal expression patterns of BrPAPs in B. rapa were analyzed. RNA-seq data derived from shoots (vernalized and non-vernalized), leaves, stems, roots, flowers, differential developmental stages of seeds, and male sterility lines were explored [22,39,[41][42][43][44]. Based on this RNA-seq data, BrPAP transcripts accumulated in all tissues examined (Figure 4), as previously reported for AtPAPs under standard growth conditions [4]. This indicated that BrPAPs might play diversified functions during development. Interestingly, BrPAP11, 24, 29, and 34 from subgroup Ia-1 showed specific expression in fertile buds in B. rapa (Figure 4). In addition to these four genes, BrPAP14 and 20 also showed high expression in fertile buds in B. rapa. This indicated that these genes specifically expressed in fertile buds might be related to pollen development. To confirm the expression patterns of these six genes, different developmental stages of floral buds from genetic male sterility (GMS) lines were collected based on a previous study ( Figure 5A) [30]. According to our semi-RT-PCR results, BrPAP11 and 24 were first detected in F2 floral buds (with tetrad pollen), and these BrPAPs were highly expressed in F3 floral buds (with pollen grain after tetrad formation but before maturation) ( Figure 5B). BrPAP14 and 29 were highly expressed in both F2 and F3 floral buds. BrPAP20 and BrPAP34 were only expressed in F4 (with pollen grain after tetrad formation but before maturation) and F3 (with mature pollen grains) floral buds, respectively ( Figure 5B). To narrow down and predict the functions of BrPAP11, 14, 20, 24, 29, and 34 during pollen development, we carried out co-expression analysis using previously published male-sterility RNA-seq data in B. rapa [22,41,43,44]. When the threshold value for PPC (Pearson's correlation coefficient) was set between 0.60 and −0.60, 1511, 979, 2355, 1593, 1886, and 1667 genes were found to be co-expressed with BrPAP11, 14, 20, 24, 29, and 34, respectively (Table S2 and Figure 6A). GO enrichment analysis was carried out to highlight the biological processes influenced by targeted BrPAPs ( Figure 6B). The pollen development process was overrepresented by the co-expressed genes of BrPAP11, 20, 24, and 34. The pollination process was only highlighted by the co-expressed genes of BrPAP20, 24, 29, and 34. Lipid metabolic and fatty acid biosynthetic processes were enriched by co-expressed genes of BrPAP14. Carbohydrate transport, pollen sperm cell differentiation, pectin catabolic, polysaccharide catabolic, cell wall organization, pollen germination, and carbohydrate metabolic processes were only overrepresented by the co-expressed genes of BrPAP20. results, BrPAP11 and 24 were first detected in F2 floral buds (with tetrad pollen), and these BrPAPs were highly expressed in F3 floral buds (with pollen grain after tetrad formation but before maturation) ( Figure 5B). BrPAP14 and 29 were highly expressed in both F2 and F3 floral buds. BrPAP20 and BrPAP34 were only expressed in F4 (with pollen grain after tetrad formation but before maturation) and F3 (with mature pollen grains) floral buds, respectively ( Figure 5B). To narrow down and predict the functions of BrPAP11, 14,20,24,29, and 34 during pollen development, we carried out co-expression analysis using previously published male-sterility RNA-seq data in B. rapa [22,41,43,44]. When the threshold value for PPC (Pearson's correlation coefficient) was set between 0.60 and −0.60, 1511, 979, 2355, 1593, 1886, and 1667 genes were found to be co-expressed with BrPAP11, 14,20,24,29, and 34, respectively (Table S2 and Figure 6A). GO enrichment analysis was carried out to highlight the biological processes influenced by targeted BrPAPs ( Figure  6B). The pollen development process was overrepresented by the co-expressed genes of BrPAP11, 20, 24, and 34. The pollination process was only highlighted by the co-expressed genes of BrPAP20, 24, 29, and 34. Lipid metabolic and fatty acid biosynthetic processes were enriched by co-expressed genes of BrPAP14. Carbohydrate transport, pollen sperm cell differentiation, pectin catabolic, polysaccharide catabolic, cell wall organization, pollen germination, and carbohydrate metabolic processes were only overrepresented by the co-expressed genes of BrPAP20.  F3, after the tetrad stage but before containing mature pollen. F4, containing mature pollen. S1-S4, floral buds from sterile plants. S1, before the tetrad stage. S2, at the tetrad stage. S3 and S4, after the tetrad stage with aborted pollen grains.  F3, after the tetrad stage but before containing mature pollen. F4, containing mature pollen. S1-S4, floral buds from sterile plants. S1, before the tetrad stage. S2, at the tetrad stage. S3 and S4, after the tetrad stage with aborted pollen grains.

Expression of BrPAP Genes in Response to Pi Deprivation in Brassica rapa
PAPs are considered to play an essential role in the processes of phosphorus foraging and recycling [5,14]. To identify the potential functions of BrPAPs related to phosphorus acquisition and utilization, germinated B. rapa seeds were sown in pots containing vermiculite and watered with Hoagland solution [28]. For Pi deprivation treatments, the KH 2 PO 4 in Hoagland solution was replaced with KCl. After treatments, the expression patterns of BrPAPs and physiological characteristics were measured. With Pi deprivation treatment, the plant height, fresh weight, leaf area, and Pi-content were all significantly decreased, while the root-shoot ratio of biomass increased compared with the normal growth condition (Figure 7). The shoot and root parts of seedlings were sampled separately for qRT-PCR analysis. Using a Ct value below 35 and fold change above two as thresholds, 21 BrPAPs responded to Pi deprivation either in shoots or roots (Figure 8). Four genes were mainly upregulated in roots under Pi deprivation ( Figure 8A). Three genes were upregulated both in shoots and roots under Pi deprivation ( Figure 8B). One gene, BrPAP12, was preferentially upregulated by Pi deprivation ( Figure 8C). Ten genes were induced by Pi deprivation in shoots ( Figure 8D). Three genes were consistently upregulated in shoots ( Figure 8E).

Expression of BrPAP Genes in Response to Pi Deprivation in Brassica rapa
PAPs are considered to play an essential role in the processes of phosphorus foraging and recycling [5,14]. To identify the potential functions of BrPAPs related to phosphorus acquisition and utilization, germinated B. rapa seeds were sown in pots containing vermiculite and watered with Hoagland solution [28]. For Pi deprivation treatments, the KH2PO4 in Hoagland solution was replaced with KCl. After treatments, the expression patterns of BrPAPs and physiological characteristics were measured. With Pi deprivation treatment, the plant height, fresh weight, leaf area, and Pi-content were all significantly decreased, while the root-shoot ratio of biomass increased compared with the normal growth condition (Figure 7). The shoot and root parts of seedlings were sampled separately for qRT-PCR analysis. Using a Ct value below 35 and fold change above two as thresholds, 21 BrPAPs responded to Pi deprivation either in shoots or roots (Figure 8). Four genes were mainly upregulated in roots under Pi deprivation ( Figure 8A). Three genes were upregulated both in shoots and roots under Pi deprivation ( Figure 8B). One gene, BrPAP12, was preferentially upregulated by Pi deprivation ( Figure 8C). Ten genes were induced by Pi deprivation in shoots ( Figure 8D). Three genes were consistently upregulated in shoots ( Figure 8E).

Identification and Analysis of BrPAPs
Based on conserved sequences and motifs, genome-wide analysis of PAPs has been performed on many plants, including Arabidopsis [2], soybean [51], maize [12], tea [15], rice [11], Brassica napus [10], and ten vegetable species [9]. In this study, 39 BrPAPs were identified and classified into three different clades with nine subgroups based on a phylogenetic analysis using PAPs from Arabidopsis and B. rapa ( Figure 2). Previously, 37 BrPAPs were identified in B. rapa [9], which could be explained by the fact that a higher genome version was used [32], leading to significantly higher accuracy of identification of BrPAPs [9]. The classification results in our study are mostly consistent with previous results from Arabidopsis and rice [2,11]. AtPAP26 from Ia-2 clustered into Ia-3 with BrPAP28 ( Figure 2) [2]. Due to the whole-genome duplication events and gene loss events that have occurred in B. rapa over evolutionary time, one Arabidopsis gene usually has one to three B. rapa homologs [32], and the results from this study also demonstrated this observation except for AtPAP10. Based on our phylogenetic analysis, AtPAP10 corresponded to five BrPAPs (BrPAP2, BrPAP8, BrPAP9, BrPAP10, and BrPAP33), indicating the functional expansion of AtPAP10 (Figure 2). Consistent with studies in maize and rice, only one single gene, BrPAP28, was identified to be the ortholog of AtPAP26 in subgroup Ia-3, indicating its highly conserved function during evolution ( Figure 2) [12].
Based on the results from Arabidopsis, rice, maize and other plants, some PAPs lack the invariant amino acid residues involved in the coordination of the di-metal nuclear

Identification and Analysis of BrPAPs
Based on conserved sequences and motifs, genome-wide analysis of PAPs has been performed on many plants, including Arabidopsis [2], soybean [51], maize [12], tea [15], rice [11], Brassica napus [10], and ten vegetable species [9]. In this study, 39 BrPAPs were identified and classified into three different clades with nine subgroups based on a phylogenetic analysis using PAPs from Arabidopsis and B. rapa ( Figure 2). Previously, 37 BrPAPs were identified in B. rapa [9], which could be explained by the fact that a higher genome version was used [32], leading to significantly higher accuracy of identification of BrPAPs [9]. The classification results in our study are mostly consistent with previous results from Arabidopsis and rice [2,11]. AtPAP26 from Ia-2 clustered into Ia-3 with BrPAP28 ( Figure 2) [2]. Due to the whole-genome duplication events and gene loss events that have occurred in B. rapa over evolutionary time, one Arabidopsis gene usually has one to three B. rapa homologs [32], and the results from this study also demonstrated this observation except for AtPAP10. Based on our phylogenetic analysis, AtPAP10 corresponded to five BrPAPs (BrPAP2, BrPAP8, BrPAP9, BrPAP10, and BrPAP33), indicating the functional expansion of AtPAP10 (Figure 2). Consistent with studies in maize and rice, only one single gene, BrPAP28, was identified to be the ortholog of AtPAP26 in subgroup Ia-3, indicating its highly conserved function during evolution (Figure 2) [12].
Based on the results from Arabidopsis, rice, maize and other plants, some PAPs lack the invariant amino acid residues involved in the coordination of the di-metal nuclear center of known PAPs [2,11,12]. In B. rapa, 31 of 39 BrPAPs contained seven invariant amino acid residues ( Table 1). The first block was absent in BrPAP14, BrPAP17, and BrPAP21. BrPAP9 did not possess the fourth and fifth blocks, and the fifth block was also lost in BrPAP7.
The third invariant amino acid residue Y shifted to T or F in BrPAP17, BrPAP21, BrPAP22, BrPAP32, and BrPAP39. In BrPAP7, the fourth and fifth conserved residues N and H were replaced by E and S, respectively. Compared with Arabidopsis, all BrPAPs had Arabidopsis homologs except for BrPAP9, and there was no conserved block absence in Arabidopsis (Tables 1 and 2) [2]. This indicates that these atypical BrPAPs might not be biochemically active PAPs and their true functions need further study.

BrPAPs and Pollen Development
Previous studies have suggested that PAPs play roles in the processes of carbon metabolism [25], cell wall synthesis [26], ROS metabolism, and response to salt stress [27].
Little is known about their functions related to pollen development, even for AtPAP6, 11,14,19,23,24, and 25, which are predominantly expressed in flowers [4]. In this study, BrPAP11, 14, 20, 24, 29, and 34 were specifically expressed in fertile floral buds in B. rapa (Figures 4 and 5). Combining the differences between fertile and sterile floral buds in the current and previous studies ( Figure 5A) [30], in addition to our semi-RT-PCR results, we concluded that these BrPAPs might function and play different roles during the process of pollen development ( Figure 5B). In support of this hypothesis, BrPAP11 was initially expressed in F2 floral buds and highly expressed in F3 floral buds, and its co-expressed genes were involved in the biological processes of organic substance transport, transmembrane transport, gametophyte development, and pollen development. This suggested that BrPAP11 might be associated with nutrition transport during pollen development ( Figure 5B). BrPAP14 was specifically expressed in F2 and F3 floral buds, and the biological processes of lipid metabolic, nitrogen compound transport, fatty acid biosynthetic, and pollen development were enriched in its co-expressed genes ( Figures 5B and 6B). This implies that BrPAP14 might be involved in secondary metabolism during pollen development. The expression levels of BrPAP20 were high in F4 floral buds (containing mature pollen grains), and carbohydrate metabolic, polysaccharide catabolic, and pectin catabolic were enriched in its co-expressed genes ( Figures 5B and 6B), suggesting that this gene might have functions during pollen maturation. The normal processing of cell wall degradation and reorganization are required for the development and functions of pollen grains [30,41]. In tobacco, NtPAP12 were found binding to the cell wall and enhanced the activities of cellulose and callose synthases [26]. The highlighted cell wall organization process by BrPAP20 and its co-expressed genes indicate that the functions of BrPAP20 might have similar functions during pollen development to NtPAP12 ( Figure 6B). Even BrPAP24 had similar expression patterns to those of BrPAP11, and pollination and pollen tube growth were only enriched by BrPAP24 and its co-expressed genes. The processes of transmembrane transport, gametophyte development, and pollen development were enriched by the co-expressed genes of BrPAP11 and BrPAP24. This observation suggests that the functions of BrAPA11 and BrPAP24 partially overlap. A similar phenomenon was observed for BrPAP14 and BrPAP29. For BrPAP34, the biological processes of pollen development, pollination, pollen tube growth, and pollen tube development were revealed by its co-expressed genes based on GO enrichment analysis, highlighting its possible roles during pollen germination and pollen development. Previously, AtPAP5 and 11 were highly expressed in anthers and involved in pollen tube growth processes [28]. In this study, BrPAP11, BrPAP24, BrPAP34, AtPAP11, and AtPAP5 were classified into subgroup Ia-1, indicating the importance of subgroup Ia-1 during pollen tube growth. Phytic acid is one of the most important components of pollen, AtPAP15 and AtPAP23 were expressed in pollen or floral buds with phytase activity in Arabidopsis [4,18], indicating the specifically expressed BrPAP11, 14,20,24,29, and 34 in fertile floral buds might have similar functions. Collectively, these analyses suggest that BrPAP11, 14, 20, 24, 29, and 34 have multiple functions during pollen development with some overlap. These findings will assist future molecular breeding for GMS in B. rapa.

Responses of Brassica rapa and BrPAPs under Pi Deprivation Conditions
Phosphorus is one of most important essential nutrients for plant growth and is involved in many metabolic processes such as respiration and photosynthesis [12]. However, the content of inorganic phosphorus (Pi) is lower than that required for plant growth in most soils, and most phosphorus is present in the Po (organic phosphorus) form or is fixed by iron, calcium, and other elements. This means that plants cannot directly utilize it [52]. APases can enable plants to use intracellular and extracellular organic phosphorus to survive in low Pi conditions [5]. Among the plant APases, PAPs are one of the most important classes and they can be induced by low Pi conditions [53]. In Arabidopsis, several AtPAPs can be induced by Pi deprivation such as AtPAP10, 11, 12, 17, and 26 [2,17,54]. In rice, the overexpression of OsPAP10a and OsPAP10c can improve the utilization of external organic phosphorus, and both of these can be induced by Pi starvation [21,22,24]. In the current study, 21 of the 39 BrPAPs were found to be responsive to Pi deprivation and these can be divided into five groups (Figure 8). The first group showed a high expression in roots under Pi deprivation, including BrPAP4, 5, 19, and 21 ( Figure 8A). BrPAP1, 23, and 33 were clustered into a second group and their expression patterns were higher in both shoots and roots under Pi deprivation compared with normal growth conditions ( Figure 8B). The third group contained only one gene, BrPAP12, which was upregulated in roots in normal conditions and shoots with Pi deprivation ( Figure 8C). Most of the Pi deprivationresponsive genes were clustered into the fourth group, which included ten BrPAPs (BrPAP3, 8, 11, 13, 14, 16, 22, 24, 36, and 38); these genes were predominantly upregulated in shoots under Pi deprivation ( Figure 8D). BrPAP17, 28 and 30 were classified into a fifth group and were upregulated in shoots in both normal and Pi deprivation conditions ( Figure 8D). In previous studies, AtPAP10 was specifically induced by Pi limitation on the root surface and was considered to be one of the major Pi starvation-induced APases in Arabidopsis [17,55]. Among BrPAPs, five homologs of BrPAP10 were found and only two of them, BrPAP8 and BrPAP33, were Pi deprivation-responsive genes (Figures 2 and 8), indicating the functional expansion of the AtPAP10 group during evolution. The expression levels of AtPAP12 increased in low phosphate conditions in Arabidopsis, and its homologs, BrPAP12 and BrPAP13, were also induced in shoots under Pi-deprivation (Figures 2 and 8) [2]. AtPAP26 functions as a dual-targeted PAP and is the predominant intracellular APase as well as a major secreted APase in Pi limited conditions [19]. BrPAP28, the homolog of AtPAP26, was only upregulated in shoots in Pi deprivation, indicating that its function had shifted compared with AtPAP26 (Figures 2 and 8E). Future studies exploring BrPAP responses to Pi deprivation will provide insights into the functions of different BrPAPs and novel gene sources for engineering crops with increased tolerance to Pi deprivation.

Conclusions
We systematically identified and analyzed the BrPAP genes from the B. rapa genome at the genome-wide level. A total of 39 BrPAPs were identified and divided into three major clades with nine subgroups. Expression profiles revealed that all BrPAPs were expressed in all examined tissues, thus indicating their various functions under standard conditions. Interestingly, BrPAP11, 14, 20, 24, 29, and 34 were specifically expressed in fertile floral buds, indicating their critical roles during pollen development. Their potential functions were further highlighted by co-expression analysis. Twenty-one BrPAPs responded to Pi deprivation in either shoots or roots. Further studies are required to determine the functions of BrPAPs. The information presented here provides a foundation for future investigations to gain a deeper understanding of the roles of each BrPAP as well as the molecular mechanisms underlying pollen development and low Pi tolerance.