The “Bcajh97-01A/B” is a GMS A/B line (sister line) of B. rapa, with the progenies segregating into sterile and fertile types in a 1:1 ratio. The sterile “Bcajh97-01A” has been demonstrated as a mutant of male meiosis cytokinesis without mature pollen, which is the only difference detected between “Bcajh97-01A” and the fertile “Bcajh97-01B” (Huang et al., 2008a, 2009). This characteristic of “Bcajh97-01A/B” GMS line is ideal to study gene expression during pollen development.

“Bcajh97-01A/B” was planted in experimental green-house of Zhejiang University. At the flowering stage, siliques at 48 h after pollination (HAP), roots, stems, leaves, and inflorescences were collected from 15 individual plants of the fertile line for reverse transcription PCR. At the same time, according to the longitudinal diameters, the flower buds were divided into five stages corresponding to five distinct pollen developmental stages (Huang et al., 2008a). The pollen development in flower buds was confirmed by microscopic examination.

Pistils in the fertile plants were collected at 1, 3, and 10 HAP. One HAP pistils correspond to the period of pollen germination on the stigma, 3 HAP pistils represent pollen tube growth phase, and 10 HAP pistils correspond to the period when the pollen tube reaches the ovule and fertilization occurs (Jiang et al., 2013). The materials were frozen in liquid nitrogen immediately after harvested and stored at −75°C for subsequent analysis.

For phytohormone treatments, one of the following solutions of 100 μM GA, 100 μM MeJA, and 100 μM ABA was sprayed on all leaves (Liu et al., 2014; Duan et al, 2015). Control plants were sprayed with distilled water. The concentration of phytohormone treatments chosen was based on the studies by Suzuki et al. (2002), Liu and Mehdy (2007), and van Hengel et al. (2004). Seedlings in three containers were used in each treatment, and each container had 12 seedlings. The treatments were repeated three times. Control and treated leaves were harvested at 0, 4, and 12 h after treatment (HAT). The samples were frozen in liquid nitrogen and stored at −80°C until use.

Total RNA was extracted from tissues mentioned above using Trizol reagent (Invitrogen, Carlsbad, CA). The RNA integrity was analyzed on agarose gel. SuperScript IV (Invitrogen) and oligo (dT) primer was used to synthesize the cDNA. The expression of cDNA was detected using BiometraT Professional Thermo cycler and Taq (Dingguo, Shanghai). Reverse transcription polymerase chain reaction (RT-PCR) was employed to detect the expression of BrAGPs in different tissues of B. rapa, while quantitative real-time PCR (qRT-PCR) was carried out to test the expression variation of BrAGPs in treatment with endogenous hormones. RT-PCR and qRT-PCR were both carried out in triplicates using gene specific primers (Supplementary data 1: Tables S1, S2). BrUBC10 was used as the reference gene. The results of qRT-PCR were calculated using the 2^−ΔΔCt method (Livak and Schmittgen, 2001) and normalized to the corresponding distilled water treatment. They were further gene-wise normalized, mean-centered and clustered hierarchically using the centroid linkage clustering method in Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/index.html).

The searching criteria for AGPs referred to the guidelines that are widely used (Schultz et al., 2002, 2004; Tan et al., 2003; Ma and Zhao, 2010). A Perl script (amino acid bias) was written to calculate the PAST (Pro, Ala, Ser, Thr) percentage for the candidate proteins (Supplementary data 1: Supplement Methods). All annotated proteins were downloaded from the Brassica database (BRAD, http://brassicadb.org, V1.2). Perl is a high-level, general-purpose, interpreted, dynamic programming language for UNIX, Windows, and OSX operating systems (http://www.activeperl.com/). For this study, the Perl language compiler ActivePerl-5.16.2.3010812913 was used.

The amino acid bias program was compiled with two different output files of the lists of proteins. The first output file named “AGPs” includes that all the proteins above 50% PAST threshold and >90 amino acid residues in length. Another output file named “AG-peptide” includes all proteins above 35% PAST and between 50 and 90 amino acid residues in length. The N-terminal and C-terminal signals of AG-peptide account for more than 50% of the deduced proteins and the length of AG-peptide is shorter than classical AGPs, so we reduced the PAST level and shortened the searching length of the sequences. Classical AGPs and lys-rich AGPs were extracted from file “AGPs,” in which lys-rich AGPs had a lys-rich domain of approximately 16 amino acid residues that was flanked on both sides by AGP glycomodules, and classical AGPs covered the other genes without the domains like fascilin-like domain and early nodulin-like domain. Pfam (http://pfam.xfam.org/) and SMART (http://smart.embl-heidelberg.de/) were used for analyzing the conserved domains. Moreover, orthologous genes of all the annotated classical AGPs, AG peptides and lys-rich AGPs in A. thaliana were screened in the proteome of B. rapa.

The signal peptide of proteins was predicted by the SignalP 4.1 Server (http://www.cbs.dtu.dk/services/SignalP/) (Petersen et al., 2011). In addition, proteins with predicted signal peptide cleavage site in the first 50 amino acids were considered as the presence of N-terminal signal sequences (Nielsen et al., 1997). To determine the presence of a C-terminal GPI anchor additional signal, all the proteins in the “AGPs” and “AG-peptide” files were subjected to the big-PI Plant Predictor (http://mendel.imp.ac.at/gpi/plant_server.html) and PSORT Prediction (http://psort.hgc.jp/form.html) (Nakai and Horton, 1999; Eisenhaber and Eisenhaber, 2007). Proteins having both signal peptide and GPI anchor were treated as “high confidence” AGPs.

The analysis of synteny using information of 24 conserved collinear blocks of ancestral karyotype (AK) in B. rapa and A. thaliana was carried out according to the previous study (Cheng et al., 2012, 2013). Chromosomal locations of the syntenic genes of A. thaliana AGPs (AtAGPs) and BrAGPs were gathered from The Arabidopsis Information Resource (TAIR, http://arabidopsis.org) and BRAD, respectively. The syntenic relationships between or inside the genomes were illustrated by Circos 5.05 (Krzywinski et al., 2009).

A set of 458 core eukaryotic genes and 458 random genes in A. thaliana (required for retained proportion analysis), were downloaded from CEGMA (http://korflab.ucdavis.edu/Datasets/cegma) (Parra et al., 2007) and were used to search for syntenic genes in the BRAD database. The orthologous genes of several AtAGPs in A. lyrata, Capsella rubella, Camelina sativa, B. napus, B. oleracea, Thellungiella halophile, Thellungiella salsuginea, Sisymbrium irio, Schrenkiella parvula, and Aethionema arabicum were detected through syntenic analysis in BRAD.

Multiple sequence alignment was performed by ClustalW (http://www.ch.embnet.org/software/ClustalW.html) in standard settings with some modification. Phylogenetic trees were constructed with both neighbor-joining (NJ) and maximum likelihood (ML) method with 1,000 bootstraps using MEGA6 software (Tamura et al., 2013). Poisson model were used in Model/Method and Uniform rates were chosen in Rates among sites. Treatment of Gaps/Missing Data was selected as Pairwise deletion. The clans consisting only one orthologous gene were treated as orphan genes, while other clans were named as clades in sequence. Motif analysis was performed in Multiple Em for Motif Elicitation (MEME) (http://meme-suite.org/tools/meme) (Bailey et al., 2009) for the identification of the short conserved motifs in the subfamilies of classical AGPs and AG peptides.

The molecular evolutionary rates between orthologous genes were estimated by calculating the ratio of nonsynonymous substitution rate (Ka) and synonymous substitution rate (Ks) between orthologous gene pairs. The protein sequences were aligned by ClustalW and the resulting alignment was used as a reference to align the nucleotide sequences. The Ka/Ks value was then determined using the YN86 module in PAML 4.7 (Yang, 2007). Significant differences were calculated by one-way ANOVA tests in SPSS 19 (IBM, USA).

Initially, 76 proteins longer than 90 amino acids were retrieved in the B. rapa genome, with biased amino acid compositions of at least 50% PAST. Similarly, 115 potential AG peptides were identified by searching for proteins between 50 and 90 amino acids in length with biased amino acid compositions of at least 35% PAST (Table 1).

Signal peptides are important for glycosylation of AGPs in the endoplasmic reticulum, while GPI anchor is required for anchoring AGPs to the plasma membrane. Therefore, these 191 proteins were further screened for the presence of a signal peptide and a GPI anchor. Fifty eight proteins longer than 90 amino acids and 24 proteins shorter than 90 amino acids were identified as they contain both the signal peptide and the GPI anchor (Table 1). By eliminating the proteins contained other domains such as fascilin-like domain and early nodulin-like domain, 22 classical AGPs, 24 AG peptides, and three lys-rich AGPs were identified as “high confidence” AGPs (Tables 1, 2). Additionally, the orthologous genes of all the annotated classical AGPs, AG peptides and lys-rich AGPs in A. thaliana were screened in the proteome of B. rapa. With this approach, 15 additional AGPs have been revealed as AGP candidates, including 11 classical AGPs, and four AG peptides (Table 2).

All AGPs in B. rapa were named after their orthologous proteins in A. thaliana. Bra038390 (81 amino acids in length) was identified as the ortholog of a classical AGP—AtAGP10. Hence, for convenience, Bra038390 was regarded as a classical AGP but not as an AG peptide. Moreover, sequence homology analysis did not result in any orthologous proteins in the annotated AtAGPs for three BrAGPs, Bra028633, Bra035584, and Bra036401, so they were designated as BrAGP59, BrAGP60, and BrAGP61, respectively. Overall, apart from the chimeric proteins, 64 AGPs were detected in B. rapa, including 33 classical AGPs, 28 AG peptides and three lys-rich AGPs (Table 2).

Most of the BrAGPs were distributed on chromosomes, with the exception of BrAGP9.2, BrAGP11.2, and BrAGP12.2, which were located on scaffolds. On average, the chromosomes contained four to nine BrAGPs, with chromosome A06 harboring only three and chromosome A03 having up to 11. The 24 conserved collinear blocks in ancestral karyotype were used to identify the syntenic relationship between BrAGPs and their corresponding orthologous genes in A. thaliana. The result illustrated that all BrAGPs are located in the same blocks as their A. thaliana orthologous genes (Figure 1), suggesting that the expansion of AGP gene family in B. rapa was depended on the WGT event.

No orthologous genes of AtAGP5, AtAGP7, AtAGP19, AtAGP41, AtAGP51, and AtAGP56 were identified in B. rapa. To investigate whether these AtAGPs were newly acquired in A. thaliana or were lost in B. rapa, their orthologous genes in other sequenced Brassica species were analyzed. Brassica species can be classified into three lineages by spectrum analysis, and most sequenced species were concentrated on Lineage I and Lineage II. Specifically, A. thaliana, A. lyrata, C. rubella, and C. sativa belonged to Lineage I, while B. rapa, B. napus, B. oleracea, T. halophile, T. salsuginea, S. irio, and S. parvula were in Lineage II (Koch and German, 2013). Analysis of these orthologous genes illustrated that AtAGP5, AtAGP7, AtAGP19, AtAGP41, and AtAGP56 exist not only in the species belonging to Lineage I, but also in the species belonging to Lineage II as well as in the primitive Brassica species, A. arabicum. This finding indicated that these genes may be lost in B. rapa (Supplementary data 1: Figure S1). The orthologs of AtAGP51 could only be detected in species in Lineage I, indicating that this gene arose after the separation of A. thaliana and B. rapa. Based on this result, the retained proportion of AGPs in B. rapa was calculated. The AGP gene family had a retained proportion of 50%, which was much higher than that of randomly selected genes (45%) and similar to that of core eukaryotic genes (52%) (Figure 2A). Particularly, classical AGPs, AG peptides, and lys-rich AGPs had a proportion of 48, 58, and 33%, respectively (Figure 2B).

As illustrated in the whole genome level, BrAGPs also displayed differentially-retained proportions among the three subgenomes, namely least fractionated subgenome (LF), medium fractionated subgenome (MF1) and most fractionated subgenome (MF2) of B. rapa (Figure 2C). The AGP gene family contained more genes of LF subgenome (63%) than those of random selected genes (59%) and core eukaryotic genes (59%), while it had similar proportion of core eukaryotic genes in the MF1 and MF2 subgenomes for AGP genes and core eukaryotic genes, respectively. This indicated that the highly retained proportion of AGPs is mainly attributed to the gene reservation in LF.

Full length protein sequences of all identified BrAGPs and their orthologs in A. thaliana were used to construct phylogenetic trees. The topologies of NJ and ML trees were mainly consistent, thus only the NJ trees are presented here (Figure 3). The 55 classical AGPs in B. rapa and A. thaliana were classified into six clades and three orphan genes with high bootstrap support (Figure 3). Motif analysis showed that the AGPs in Clade I mainly have two types of structural features. AGP2, AGP3, AGP4, and AGP7 contained motif 1, motif 4 and motif 3, while AGP1, AGP5, and AGP10 included motif 1 and motif 4. Clade II was constituted of AGP54 and AGP57 with motif 1 and motif 2. Clade III had three members, AGP6, AGP11, and AGP50, in which AGP6 and AGP11 had motif 1 and motif 4, similar to some of the members in Clade I, while AGP50 only had motif 1. AGP9, AGP58, and AGP62 belonged to Clade IV with AGP58 having only motif 4, while AGP9 and AGP62 having motif 1 and motif 4. Clade V was composed of AGP25, AGP26, AGP27, which only hold motif 1. Tandem gene pair, AGP52 and AGP53 contained motif 4. Phylogenetic analysis indicated that the tandem repeat event occurred after the separation of B. rapa and A. thaliana, as the AGPs belonging to B. rapa had a closer resemblance to each other than to the orthologs in A. thaliana.

Phylogenetic analysis of AG peptides led to the division of three clades and five orphan genes, and almost all AG peptides had motif 6 and motif 7 (Figure 4). Four members of Clade I, AGP12, AGP13, AGP14, and AGP21, had motif 10 in addition to motif 6 and motif 7. AGP16, AGP20, AGP22, and AGP41 in Clade II all contained each motif 6, motif 7 and motif 8. Clade III included AGP23 and AGP43, which all had motif 8.

Estimation of molecular evolution rates forms an important basis for our understanding of the evolution processes in B. rapa. YN86 module was used for estimating the molecular evolutionary rates of the orthologous gene pairs among AGPs (Supplementary data 1: Table S3). All three types of AGPs encoding genes were evolved at a Ka/Ks value lower than 1, indicating that they all have evolved through purifying selection (Figure 5A). Classical AGPs showed an average Ka/Ks value of 0.30, AG peptides exhibited and Lys-rich AGPs were evolved at an average Ka/Ks value of 0.20 and 0.38, respectively. Statistical analysis illustrated that classical AGPs have evolved faster than AG peptides, but there is no significant differences existed between classical AGPs and lys-rich AGPs. Furthermore, the average molecular evolutionary rates of each clade were compared. In classical AGPs, Clade VI had a significantly higher Ka/Ks value than the other five clades (Figure 5B). In AG peptides, Clade I, Clade II, and Clade III shared a similar molecular evolutionary rate, while the orphan genes evolved at a faster rate of evolution (Figure 5C).

Half of the classical BrAGPs showed expression in all of the tissues detected, including root, stem, leaf, inflorescence and silique. Seven AGPs were expressed in more than one but not all tissues. BrAGP6, BrAGP11.1, BrAGP11.2, and BrAGP54.1 expression was restricted in inflorescence (Figure 6A). Eleven of the 24 AG peptides encoding genes were expressed in all tissues with BrAGP12.2, BrAGP13.2, BrAGP15.1, BrAGP15.2, BrAGP22.1, and BrAGP22.2 displaying similar expression levels among different tissues. Half of the genes encoding AG peptides were expressed in different tissues, of which BrAGP15.3 was inflorescence-specific and BrAGP14 was expressed exclusively in the stem. The expression of BrAGP44 was not detected in any of the tissues (Figure 6A).

Among the three lys-rich AGPs encoding genes, BrAGP18.1 and BrAGP18.2 were expressed in all tissues, while BrAGP17 was only detectable in stem, inflorescence and silique (Figure 6A). Interestingly, most paralogous gene pairs displayed different expression patterns. Specifically, from the classical AGPs group, the gene pairs of BrAGP1.1 and BrAGP1.2, BrAGP2.1 and BrAGP2.2, BrAGP3.1 and BrAGP3.2, BrAGP4.1 and BrAGP4.2/BrAGP4.3 had differential expression patterns. Likewise, the paralogous gene pairs of BrAGP12.1 and BrAGP12.2, BrAGP13.1 and BrAGP13.2, BrAGP15.1, BrAGP15.2 and BrAGP15.3, BrAGP16.1 and BrAGP16.2, BrAGP21.1, BrAGP21.2 and BrAGP21.3, BrAGP40.1 and BrAGP40.2, BrAGP43.1 and BrAGP43.2, encoding for AG peptides, showed differential expression patterns.

For the five inflorescence-specific genes, including four classical AGPs and one AG peptide encoding genes, their expression were further analyzed among different developmental stages of pollen in subdivided flower buds from the “Bcajh97-01A/B”GMS line (Figure 6B). These five BrAGPs were expressed in the fertile but not in sterile flower buds. BrAGP6 and BrAGP11.2 were only expressed at Stage I which corresponds to the mature pollen stage, while BrAGP11.1, BrAGP15.3, and BrAGP54.1 were detected at Stage III to V, which means that they are expressed throughout the uninucleate pollen stage to the mature pollen stage. Since some genes expressed in mature pollen usually continue to be expressed in germinating pollen and pollen tubes, the expression levels of the five BrAGPs were investigated in pistils at three different time points after pollination in a fertile line. BrAGP11.1 and BrAGP54.1 were found to be expressed in pistils during the whole fertilization processes; however the other three genes were not detected in the pistils (Figure 6C).

Phytohormones serve as key signals for the plants to respond to environmental stimuli (Wania et al., 2016). In order to investigate whether BrAGPs were involved in phytohormone responses, we chose 15 classical AGPs, 13 AG peptides and two lys-rich AGPs, which according to our qRT-PCR displayed expression in the leaf. Three vital phytohormones, ABA, GA, and MeJA, involved in abiotic and biotic stress resistance, were used in this study to investigate the leaf BrAGPs expression in response to phytohormone treatments.

Following ABA treatment, 22 BrAGPs were up-regulated, including 13 classical AGPs, eight AG peptides and one lys-rich AGPs encoding genes (Figure 7A). Expression of 11 classical AGPs, four AG peptides, and a lys-rich AGPs encoding genes were markedly increased at 4 HAT and then recovered at 12 HAT; whereas, the expression levels of two classical AGPs and three AG peptides encoding genes were increased at 12 HAT. Eight BrAGPs, including two classical AGPs, five AG peptides, and one lys-rich AGPs encoding genes were down-regulated upon ABA treatment. Their expression levels decreased mainly at 4 HAT and didn't recover at 12 HAT, while BrAGP22.3 decreased at 12 HAT. Except BrAGP2.1 and BrAGP2.2 that shared similar responses, most paralogous gene pairs displayed different expression patterns after ABA treatment.

Up-regulation of 12 classical AGPs, seven AG peptides, and one lys-rich AGPs genes was observed after GA treatment (Figure 7B). The expression levels of three classical AGPs encoding genes were increased at 4 HAT and subsequently at 12 HAT. Three classical AGPs and three AG peptides encoding genes were also up-regulated at 4 HAT but their transcript levels dropped back at 12 HAT. There was no effect on the expression of genes encoding for six classical AGPs, four AG peptides, and one lys-rich AGP, until 12 HAT. In contrast, the expression of encoding genes for three classical AGPs, six AG peptides, and one lys-rich AGP was down-regulated. The expression level of BrAGP2.2 decreased at 4 HAT and recovered at 12 HAT, while BrAGP12.1 was down-regulated at 12 HAT. The remaining genes were down-regulated at 4 HAT and sustained low expression levels at 12 HAT. Although the expression of most paralogous gene pairs was different, the response pattern between BrAGP40.1 and BrAGP40.2, and between BrAGP22.2 and BrAGP22.3, was similar (Figure 7B).

Nineteen BrAGPs were up-regulated in MeJA treatment. Two classical AGPs and two AG peptides encoding genes were induced at 12 HAT, while 10 classical AGPs, four AG peptides and one lys-rich AGPs encoding genes, were up-regulated at 4 HAT and recovered at 12 HAT (Figure 7C). In addition, 11 AGPs were down-regulated, following MeJA treatment. One classical AGP and two AG peptides were down-regulated at 4 HAT and recovered at 12 HAT, while the expression levels of the encoding genes for two classical AGPs, five AG peptides, and one lys-rich AGP decreased at 4 HAT and didn't recover at 12 HAT. BrAGP1.1 and BrAGP1.2 showed the same response to ABA, but the rest of the paralogous gene pairs displayed different expression patterns.

Interestingly, the ABA-, GA-, and MeJA-induced expression of BrAGPs were largely overlapped. Seventeen BrAGPs were up-regulated by all three phytohormones, including genes encoding 10 classical AGPs, six AG peptides, and one lys-rich AGPs. In addition, BrAGP2.2 and BrAGP26 were up-regulated by both ABA and MeJA, and BrAGP13.1 was up-regulated both by ABA and GA (Figure 7D). Five AG peptides encoding genes and one lys-rich AGP encoding gene were down-regulated by all three phytohormones. However, BrAGP3.2 and BrAGP60 were solely down-regulated by ABA and MeJA, the expression levels of BrAGP10.1 and BrAGP12.2 decreased only in response to GA and MeJA, respectively (Figure 7E).

Comparison of the gene expression induction kinetics of different hormones revealed that the overlap between the response to ABA and MeJA is much bigger than that between ABA and GA, or MeJA and GA (Figure 7). In general, most of the transcriptional response of BrAGPs to ABA and MeJA treatments was rapid and transient with expression changes detected at 4 HAT, while these of GA were relatively slow (expression changes were not detectable until 12 HAT).

AGPs are characterized by the repetitive nature of their protein backbones with less than 40% similarity between them (Schultz et al., 2002). The low similarity makes it difficult to be identified using the basic local alignment search tool (BLAST) or profile hidden Markov model (HMM) methods. “Hyp contiguity hypothesis” has provided a guidance to identify AGPs by Perl scripts (Kieliszewski and Lamport, 1994). However, it may lead to variable results if different programs and standards are employed (Schultz et al., 2002; Showalter et al., 2010). Here, based on the threshold values set in Showalter et al. (2010), a whole genome wide screening of classical AGPs, AG peptides, and lys-rich AGPs was performed in B. rapa, resulting in the identification of 64 BrAGPs including 49 with high confidence.

B. rapa contains more AGPs than A. thaliana, especially classical AGPs and AG peptides that is mainly attributed to the WGT event. However, we could not identify orthologs for three BrAGPs in A. thaliana, and orthologs for six AtAGPs in B. rapa. Ortholog comparison with closely related Brassica species revealed that five AtAGPs are lost in B. rapa, while one AtAGP and three BrAGPs were independently evolved after the separation of A. thaliana and B. rapa. These results indicate that besides the WGT event, other processes, for instance, dosage effect, might have also contributed to the evolution of the AGP family. In this study we have shown that classical AGPs, AG peptides and lys-rich AGPs are highly represented in the genome of B. rapa. Similar results were also reported by two previous studies where 33 FLAs and 52 ENODLs were identified in B. rapa. The corresponding retained proportions for these two types of chimeric AGPs were described to be 52 and 67%, respectively (Jun and Xiaoming, 2012; Li et al., 2013).

According to previous study, the over retained genes or gene families were always dosage sensitive or usually undergone subfunctionalization and/or neofunctionalization (Kim et al., 2014). In our investigation of the expression of BrAGPs, we revealed that 4/5 “high confidence” classical AGP and 6/7 AG peptide paralogous gene pairs had distinct expression patterns in different tissues. Such results indicated that the high retained proportions of classical AGPs and AG peptides are probably due to subfunctionalization and/or neofunctionalization, but not due to the dosage effect, which may yield increased expression of a gene (Edger and Pires, 2009). However, the result of molecular evolutionary rate analysis indicated that these AGPs were all under purifying selection, illustrating that neofunctionalization may have hardly occurred in these proteins. We thus suggest that subfunctionalization is the most likely mechanism that facilitated the retention of duplicated genes in classical AGPs and AG peptides. Unfortunately, in the previous studies on FLAs and ENODLs, the paralogous gene pairs in B. rapa and A. thaliana were treated as one gene for their gene expression pattern (Jun and Xiaoming, 2012; Li et al., 2013). Therefore, we were unable to verify whether the high retained proportion of the chimeric AGPs are also depended on subfunctionalization. Relevant research work should be focused on this issue in the future.

Comprehensive studies on gene expression can provide useful information for predicting gene function. In this study, we first investigated the expression of BrAGPs in five different tissues and found tissue-specific and non-tissue specific genes. To identify candidate AGPs that may be involved in pollen development, those inflorescence-specific BrAGPs including four classical AGPs encoding gene BrAGP6, BrAGP11.1, BrAGP11.2, BrAGP54.1, and one AG peptide encoding gene BrAGP15.3, were further investigated for their in male fertile/sterile flower buds at five different developmental stages. All these genes were specifically expressed in fertile flower buds at the late developmental stages, indicating a role in pollen formation or pollen function. Some genes expressed in mature pollen usually continued to be expressed in germinating pollen and pollen tubes, showing that they are specifically involved in the regulation of the subsequent pollination and fertilization processes (Pina et al., 2005; Becker and Feijo, 2007; Wang et al., 2008). We also detected the expression of these genes in pistils at three different time points after pollination. BrAGP11.1 and BrAGP54.1 were shown to be expressed in pistils during the whole fertilization processes, which indicated that they might also function in subsequent sexual reproductive events. A previous study has demonstrated that AtAGP6 and AtAGP11 are the only two pollen-specific classical AGPs and are essential for pollen and pollen tube function in A. thaliana (Pereira et al., 2006; Levitin et al., 2008). Recent studies using qRT-PCR and in situ hybridization in our laboratory also revealed that the homologous gene BrAGP6 (formerly named as BcMF18) was specifically expressed in pollen (unpublished data), and BrAGP11.1 (or BcMF8) was expressed in the developing pollen and the pollen tube (Huang et al., 2008b; Lin et al., 2014). The functional disruption of BrAGP6 by antisense RNA technology caused collapse of the pollen due to the absence of cellular content and nucleus, and failure of intine layer formation due to lack of cellulose deposition (unpublished data). The inhibition of BrAGP11.1 resulted in slipper-shaped and bilaterally sunken pollen with abnormal intine development and aperture formation, an arrest of pollen germination and unstable pollen tube formation (Lin et al., 2014). Although the functions of BrAGP15.3 and BrAGP54.1 are not fully understood, their expression patterns are similar to the known pollen-related AGP genes. Therefore, we suggest that these two genes may play important roles in pollen development and/or fertilization. Our results also revealed that only a small number of BrAGPs was inflorescence-specific. This is in agreement with a previous expression analysis of the FLAs and ENODLs in B. rapa, which also showed that only two BrENODL genes and five FLA genes were specifically expressed in inflorescence (Jun and Xiaoming, 2012; Li et al., 2013). Similarly, in A. thaliana, only a subset of AGP genes is expressed in pollen grain and pollen tubes (Pereira et al., 2006). Therefore, we speculate that, despite of the diversity of AGPs, only a small number of AGPs are endowed with a specific and defined role in male sexual reproduction in B. rapa as those in A. thaliana.

AGPs are also implicated in signal recognition and transduction (Suzuki et al., 2002; van Hengel and Roberts, 2003), although the exact mechanisms are still elusive. The mRNA level of AtAGP31 decreases at 8 h after ABA, MeJA, and wounding treatment (Liu and Mehdy, 2007), while AtAGP30 specifically responds to ABA treatment (van Hengel and Roberts, 2003). Moreover, a mutation of this gene showed suppression of the ABA-induced delay in germination and alteration of the expression of some ABA-regulated genes (van Hengel and Roberts, 2003). In rice, two classical AGPs, OsAGP1, and OsAGP15 were significantly up-regulated under the ABA treatment, indicating that they may response to ABA (Ma and Zhao, 2010). Moreover, AGPs are also involved in response to GA, for example, AGPs are suggested to play important roles in GA-induced α-amylase production in barley aleurone cells (Suzuki et al., 2002). The ATH1 microarray analysis results (GEO accession number: GSE39384) also revealed that many AGP-encoding genes are regulated by ABA, GA, and MeJA in A. thaliana (Supplementary data 2: Table S4). A similar trend was observed for some homologous AGPs between B. rapa and A. thaliana. For example, AGP1 and AGP2 were up-regulated in A. thaliana after 3 h treated by ABA, while their orthologs in B. rapa were also up-regulated at 4 HAT by ABA (Figure 7). The conserved response to the phytohormone treatments in different plant species indicates that AGPs may be important functional genes for phytohormone signaling pathways. However, many homologous BrAGP genes show diverse change trends under phytohormone treatments with their A. thaliana homologs. This result demonstrated that the differentiation of the AGP genes might be intended to acclimatize the plant to specific environmental changes. In our study, 17 BrAGPs were up-regulated and six were down-regulated by all three types of phytohormones, ABA, GA, and MeJA treatments. These genes comprise the majority (approximately 76%) of AGPs which are expressed in the leaf of B. rapa. These findings indicate that quite a few AGPs might be directly associated or indirectly interact with the hormone signaling network.

The molecular mechanism by which AGPs are involved in phytohormone signaling remains to be elucidated. The relative rapid expression of AGPs under exogenous phytohormone treatments suggested that they may function in maintaining the organism's homeostasis rather than solely being structural components. In plants, sugar levels reflect the plant's physiological conditions and sugar signaling is closely associated with multiple phytohormone signaling cascades and abiotic stress signaling (Smeekens, 2000; Gibson, 2005; Rolland et al., 2006). Interestingly, sugar can induce ABA in the sugar-ABA signaling cascade (Laby et al., 2000). Changes in the expression of each BrAGP gene by exogenous ABA treatment indicates that the AGPs may participate in the signaling cascade by the dynamic release and uptake of glycosyl radicals. It was clear that there is complex cross-talk among plant hormones, which are important in plant growth regulation and defense responses (Kohli et al., 2013; Wania et al., 2016). In this study, we have shown that there is a large overlap of the BrAGPs induced or suppressed by ABA, GA, and MeJA treatments. It indicates that these BrAGPs might participate in different signaling pathways or in overlapping processes that are controlled by phytohormonal interactions.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.