Identification of Differentially Expressed Genes Involved in the Molecular Mechanism of Pericarp Elongation and Differences in Sucrose and Starch Accumulation between Vegetable and Grain Pea (Pisum sativum L.)

Pea (Pisum sativum L.), as a major source of plant protein, is becoming one of the major cultivated crop species worldwide. In pea, the pericarp is an important determinant of the morphological characteristics and seed yield. To investigate the molecular mechanism of pericarp elongation as well as sucrose and starch accumulation in the pods of different pea cultivars, we performed transcriptomic analysis of the pericarp of two types of pea cultivar (vegetable pea and grain pea) using RNA-seq. A total of 239.44 Gb of clean sequence data were generated, and were aligned to the reference genome of Pisum sativum L. In the two samples, 1935 differentially expressed genes (DEGs) were identified. Among these DEGs, three antioxidant enzyme superoxide dismutase (SOD) were detected to have higher expression levels in the grain pea pericarps at the pod-elongating stages. Otherwise, five peroxidase (POD)-encoding genes were detected to have lower expression levels in the vegetative pericarps at the development stage of pea pod growth. Furthermore, genes related to starch and sucrose metabolism in the pea pod, such as SUS, INV, FBA, TPI, ADPase, SBE, SSS, and GBSS, were found to be differentially expressed. The RNA-seq data were validated through real-time quantitative RT-PCR of 13 randomly selected genes. Our findings provide the gene expression profile of, as well as differential expression information on, the two pea cultivars, which will lay the foundation for further studies on pod development and nutrition accumulation in the pea and provide valuable information for pea cultivar improvement.


Introduction
Pea (Pisum sativum L.), belonging to the legume family, is an annual or biennial crop that originated in the southern part of Europe, the Mediterranean coast, and West Asia [1]. Pea seeds can not only be used as food but also as fodder, and their tender seeds, pods, and seedlings are consumed as edible vegetables. Pea is a widely cultivated crop species and a major source of plant protein for both animal and human consumption [2]. Moreover, the seeds are also used in treating wrinkled skin, diabetes, 1.5 cm in width for the grain pea and vegetable pea, respectively. During pea growth, the size of the pericarp increased, and the largest size was reached at 28 DAP, at which the length and width of the WDZY-04 pea pericarp were 7.9 and 1.54 cm, and those of WDZY-14 were 8.23 and 1.56 cm, respectively. As for the grain pea, WDZY-04, the pericarp size changed most significantly from 7 to 28 DAP, and the length of the pericarp at 28 DAP was almost 1.27 times that at 7 DAP. This result means that, compared with vegetable pea, the pericarp of grain pea elongated more quickly from 7 to 28 DAP ( Figure 1D). Meanwhile, during the growth development of the two pea cultivars, the trend of the pericarp weight was increased from 7 to 21 DAP and then decreased. The pericarp weight of WDZY-14 at each stage was heavier than that of WDZY-04 except at 35 DAP.
The pericarp sucrose contents of the two cultivars differed at each development stage. WDZY-14 accumulated more sucrose than WDZY-04 during growth development ( Figure 1B). However, the two cultivars had a similar trend at each development stage, which is that the content of sucrose increased from 7 to 14 DAP, where it reached its peak, followed by a decreasing trend until 28 DAP, whereby it began to increase. Meanwhile, the starch content of the two cultivars presented a single peak pattern during their growth development ( Figure 1C). The starch content of the cultivar WDZY-04 was higher than that of WDZY-14. width of the WDZY-04 pea pericarp were 7.9 and 1.54 cm, and those of WDZY-14 were 8.23 and 1.56 cm, respectively. As for the grain pea, WDZY-04, the pericarp size changed most significantly from 7 to 28 DAP, and the length of the pericarp at 28 DAP was almost 1.27 times that at 7 DAP. This result means that, compared with vegetable pea, the pericarp of grain pea elongated more quickly from 7 to 28 DAP ( Figure 1D). Meanwhile, during the growth development of the two pea cultivars, the trend of the pericarp weight was increased from 7 to 21 DAP and then decreased. The pericarp weight of WDZY-14 at each stage was heavier than that of WDZY-04 except at 35 DAP.

RNA-seq of the Pea Pericarp Transcriptome
The total RNA was extracted from the pea pericarp of WDZY-14 and WDZY-04 at five growth developmental stages. Then, cDNA libraries were prepared and sequenced using an Illumina HiSeq 4000 platform. We subsequently obtained transcriptomic data from 30 samples that contained the two cultivars of pea pericarp at five stages and their biological replicates. As a result, 239.4 Gb of clean data were produced, and 1641,965,412 clean reads were generated after removing adaptor sequences, ambiguous reads, and low-quality reads ( Table 1). The quality of base calling was mostly above Q30, with >96% of the reads having a quality score above Q30. The GC content ranged from 43.05% to 44.16%. The distribution of the base qualities and base percentage composition of the reads of each sample are shown in Figures S1 and S2. Moreover, a coverage analysis of the genes and an assessment of the sequencing randomness were conducted in our study to evaluate the quality of the sequencing. Each clean read for each sample was evenly distributed in the gene body ( Figure S3). Meanwhile, Pearson's correlation analysis was performed to evaluate the reproducibility of the biological replicates. As shown in Figure S4, the correlations between samples among the same biological replicates were good, with a value ranging from 0.91 to 1.000. However, the correlations between samples of the biological replicates varied significantly, with a value ranging from 0.04 to 0.88. Meanwhile, qRT-PCR analysis was used to validate the quality of the RNA-seq data. A total of 13 genes were selected. As expected, most of these candidate genes had similar expression tendencies (Figure 2A). While the exact change did not exactly match that of the others, the expression trends of all genes from the qRT-PCR and RNA-seq analyses were largely consistent (Pearson's correlation coefficient, R 2 = 0.85) ( Figure 2B). The strong correlation between the RNA-seq and qRT-PCR data indicates the reliability of the transcriptomic profiling data.
All of the clean reads were then mapped to the reference genome of Pisum sativum Linn (downloaded from NCBI: https://www.ncbi.nlm.nih.gov/nuccore/PUCA000000000.1/), and more than 85% of the clean reads perfectly matched the reference genome (Table 1). Moreover, a coverage analysis of the genes and an assessment of the sequencing randomness were conducted in our study to evaluate the quality of the sequencing. Each clean read for each sample was evenly distributed in the gene body ( Figure S3). Meanwhile, Pearson's correlation analysis was performed to evaluate the reproducibility of the biological replicates. As shown in Figure S4, the correlations between samples among the same biological replicates were good, with a value ranging from 0.91 to 1.000. However, the correlations between samples of the biological replicates varied significantly, with a value ranging from 0.04 to 0.88. Meanwhile, qRT-PCR analysis was used to validate the quality of the RNA-seq data. A total of 13 genes were selected. As expected, most of these candidate genes had similar expression tendencies (Figure 2A). While the exact change did not exactly match that of the others, the expression trends of all genes from the qRT-PCR and RNA-seq analyses were largely consistent (Pearson's correlation coefficient, R 2 = 0.85) ( Figure 2B). The strong correlation between the RNA-seq and qRT-PCR data indicates the reliability of the transcriptomic profiling data.
(A) All of the clean reads were then mapped to the reference genome of Pisum sativum Linn (downloaded from NCBI: https://www.ncbi.nlm.nih.gov/nuccore/PUCA000000000.1/), and more than 85% of the clean reads perfectly matched the reference genome (Table 1).

Analysis of the Expression Level and Differentially Expressed Genes (DEGs)
To quantify the expression levels of the transcripts, the Bowtie 2 program was used with RSEM [26,27]. Then, the numbers of mapped reads and the FPKM (fragments per kb per million reads) values were obtained for the following analysis. The statistical results on the expression level (FPKM) of the transcripts for each sample are shown in Table S1.
The genes that were differentially expressed in the vegetable pea pod and grain pea pod at the five growth developmental stages were compared using |log2(fold change)| ≥ 1 and FDR (false discovery rate), with ≤0.05 as a significant cutoff criterion. After data filtering, we detected 6842, 6287, 8767, 8101, and 8417 differentially expressed genes (DEGs) in the two cultivars at the five stages (Table S2). Among these DEGs, 2523, 3036, 3562, 3965, and 3766 were log2 fold change ≥1 (more expressed in the vegetable pea pericarp), and 4319, 3215, 5205, 4136, and 4651 were log2 fold change ≤ −1 (less expressed in the vegetable pea pericarp) (Table S2). A total of 1935 DEGs were common in both cultivars and were putatively considered to be associated with the phenotypic trait differences in this species ( Figure 3A). Of the DEGs in the vegetable pea pericarp, 730 out of a total of 1935 were more expressed and 1158 were less expressed in the vegetable pea pericarp.

Analysis of the Expression Level and Differentially Expressed Genes (DEGs)
To quantify the expression levels of the transcripts, the Bowtie 2 program was used with RSEM [26,27]. Then, the numbers of mapped reads and the FPKM (fragments per kb per million reads) values were obtained for the following analysis. The statistical results on the expression level (FPKM) of the transcripts for each sample are shown in Table S1.
The genes that were differentially expressed in the vegetable pea pod and grain pea pod at the five growth developmental stages were compared using |log2(fold change)| ≥ 1 and FDR (false discovery rate), with ≤0.05 as a significant cutoff criterion. After data filtering, we detected 6842, 6287, 8767, 8101, and 8417 differentially expressed genes (DEGs) in the two cultivars at the five stages (Table S2). Among these DEGs, 2523, 3036, 3562, 3965, and 3766 were log2 fold change ≥1 (more expressed in the vegetable pea pericarp), and 4319, 3215, 5205, 4136, and 4651 were log2 fold change ≤ −1 (less expressed in the vegetable pea pericarp) (Table S2). A total of 1935 DEGs were common in both cultivars and were putatively considered to be associated with the phenotypic trait differences in this species ( Figure 3A). Of the DEGs in the vegetable pea pericarp, 730 out of a total of 1935 were more expressed and 1158 were less expressed in the vegetable pea pericarp.
In order to understand the expression pattern of DEGs in the five stages for the two cultivars, we conducted hierarchical clustering for the DEGs for the five compared groups ( Figure 4). The differentially expressed genes in the five groups were mainly classified into high-expression genes (red) and low-expression genes (green). We grouped genes with a similar expression pattern into a set and used six, five, four, four, and three model profiles to summarize the expression pattern of this subcluster of genes. The DEGs in the five groups fluctuated obviously (more expressed or lower expressed). The expression pattern of the sub-cluster genes in each group is shown in Figure S5. In order to understand the expression pattern of DEGs in the five stages for the two cultivars, we conducted hierarchical clustering for the DEGs for the five compared groups ( Figure 4). The differentially expressed genes in the five groups were mainly classified into high-expression genes (red) and low-expression genes (green). We grouped genes with a similar expression pattern into a set and used six, five, four, four, and three model profiles to summarize the expression pattern of this subcluster of genes. The DEGs in the five groups fluctuated obviously (more expressed or lower expressed). The expression pattern of the sub-cluster genes in each group is shown in Figure  S5.

GO and KEGG Classification
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed the biological processes, cellular components, molecular functions, and metabolic pathways associated with the identified transcripts from the vegetable pea pericarp and grain pea pericarp. All of the DEGs were divided into three GO categories: Biological process, cellular component, and molecular function ( Figure 3B). In the GO annotation, 45, 44, 46, 45, and 45 terms were categorized for the DEGs of developmental stages I-V, respectively. In the biological process, 'metabolic process', 'cellular process', and 'single-organism process' were the most highly In order to understand the expression pattern of DEGs in the five stages for the two cultivars, we conducted hierarchical clustering for the DEGs for the five compared groups (Figure 4). The differentially expressed genes in the five groups were mainly classified into high-expression genes (red) and low-expression genes (green). We grouped genes with a similar expression pattern into a set and used six, five, four, four, and three model profiles to summarize the expression pattern of this subcluster of genes. The DEGs in the five groups fluctuated obviously (more expressed or lower expressed). The expression pattern of the sub-cluster genes in each group is shown in Figure  S5.

GO and KEGG Classification
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed the biological processes, cellular components, molecular functions, and metabolic pathways associated with the identified transcripts from the vegetable pea pericarp and grain pea pericarp. All of the DEGs were divided into three GO categories: Biological process, cellular component, and molecular function ( Figure 3B). In the GO annotation, 45, 44, 46, 45, and 45 terms were categorized for the DEGs of developmental stages I-V, respectively. In the biological process, 'metabolic process', 'cellular process', and 'single-organism process' were the most highly

GO and KEGG Classification
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed the biological processes, cellular components, molecular functions, and metabolic pathways associated with the identified transcripts from the vegetable pea pericarp and grain pea pericarp. All of the DEGs were divided into three GO categories: Biological process, cellular component, and molecular function ( Figure 3B). In the GO annotation, 45, 44, 46, 45, and 45 terms were categorized for the DEGs of developmental stages I-V, respectively. In the biological process, 'metabolic process', 'cellular process', and 'single-organism process' were the most highly represented terms, and 'membrane part', 'membrane', 'organelle', 'binding', 'catalytic activity', and 'transporter activity' were the most enriched in the cellular component and molecular function of the level-two GO term ( Figure 3B).
To investigate the DEG-related pathways, we conducted KEGG annotation of these transcripts. As a result, 19 pathways belonging to 5 categories were retrieved for each compared group. Among them, 'carbohydrate metabolism', 'biosynthesis of other secondary metabolites', 'signal transduction', 'energy metabolism', 'folding, sorting and degradation', 'amino acid metabolism', and 'translation' were the top three annotation terms in each group ( Figure S6). Similarly, we also analyzed the total DEGs (higher expressed and lower expressed in the vegetable pea pericarp, respectively), and the top three level-two GO terms of more expressed DEGs were 'biosynthesis of other secondary metabolites', 'carbohydrate metabolism', 'signal transduction', 'folding, sorting and degradation', 'amino acid metabolism', 'lipid metabolism', and 'translation'. The category with the highest number of lower expressed DEGs in the five groups was 'carbohydrate metabolism', followed by 'energy metabolism', 'amino acid metabolism', 'signal transduction', 'translation', and 'folding, sorting, and degradation' ( Figure S6).

DEGs Related to Pod Elongation
Reactive oxygen species (ROS) play an important role in such plant functions as cell wall loosening and elongation [28]. The antioxidant enzymes superoxide dismutase (SOD) and peroxidase (POD) were found to be associated with pod growth through the regulation of ROS generation and transformation [12]. In our study, we compared the difference in growth between two cultivars of pericarp and identified the DEGs associated with pericarp elongation. As a result, three superoxide dismutase (SOD)-encoding genes were differentially expressed ( Table 2). The gene XLOC_013249 (NR database_id: CAD42655.1), XLOC_030516 (NR database_id: XP_013461079.1, and XLOC_038705 (NR database_id: XP_004508271.1) were more highly expressed in the grain pea pericarp than in vegetable pea pericarp at stages III and V ( Table 2). The expression levels of these three genes were not significantly different at stages I and II (supplementary dataset 1).  For peroxidase (POD), we identified 12 POD-encoding genes that included 8 DEGs in vegetable pea and grain pea. Five POD-encoding genes, XLOC_016611 (NR database_id: BAD97435.1), XLOC_006267 (NR database_id: BAD97438.1), XLOC_018196 (NR database_id: XP_010104370.1), XLOC_007148, and XLOC_006821 (NR database_id: BAD97436.1), were more highly expressed in grain pea at pod development stages I-V (Table 2). However, the three genes, XLOC_034947 (NR database_id: BAD97436.1), XLOC_037479 (NR database_id: BAD97439.1), and XLOC_031402 (NR database_id: CAA09881.1), exhibited the opposite trend at pod development stage IV, but these three genes were also more highly expressed in grain pea at stages I and III (Table 2).
In order to further investigate the relationship between SOD or POD and pod elongation, Pearson's correlation coefficient analysis was conducted using data of the gene expression level and plant growth. We selected the coefficients related to the phenotypes that met the requirements of statistical significance (multiple testing corrections and an adjusted p-value < 0.05). One SOD-encoding gene showed a significant positive correlation with pod elongation at WDZY-04, and three POD-encoding genes showed a significant negative correlation with pod elongation at WDZY-04 and WDZY-14, respectively (Supplementary dataset 2). To show the correlation clearly, we selected data of one SOD-encoding gene (gene_id: XLOC_013249) and one POD-encoding gene (gene_id: XLOC_007148) to perform the scatter diagram ( Figure 5A). SOD-encoding gene showed a significant positive correlation with pod elongation at WDZY-04, and three POD-encoding genes showed a significant negative correlation with pod elongation at WDZY-04 and WDZY-14, respectively (Supplementary dataset 2). To show the correlation clearly, we selected data of one SOD-encoding gene (gene_id: XLOC_013249) and one POD-encoding gene (gene_id: XLOC_007148) to perform the scatter diagram ( Figure 5A).

DEGs Related to Pod Sucrose Metabolism
Sucrose synthase (SUS) and invertase (INV) play an important role in sucrose metabolism. All of the sugar precursors that participate in sucrose metabolism must be decomposed to hexoses, such as glucose and fructose, or the hexoses (e.g., UDP-glucose) must be ramified by SUS or INV [29]. A total of three encoding sucrose synthase genes were deemed to be the significant DEGs in the two compared groups. The three genes, XLOC_034373 (NR databse_id: O24301.1), XLOC_025267 (NR databse_id: CAC32462.1), and XLOC_002091 (NR databse_id: XP_003591492.2), were more highly expressed in the vegetable pea pericarp at development stages I and V, respectively. The expression level of these three genes in the vegetable pea pericarp were 1.81-, 3.12-, and 1.66-fold higher than that of the grain pea pericarp (Table S3, Figure 5B, Supplementary Dataset 1).

DEGs Related to Pod Sucrose Metabolism
Sucrose synthase (SUS) and invertase (INV) play an important role in sucrose metabolism. All of the sugar precursors that participate in sucrose metabolism must be decomposed to hexoses, such as glucose and fructose, or the hexoses (e.g., UDP-glucose) must be ramified by SUS or INV [29]. A total of three encoding sucrose synthase genes were deemed to be the significant DEGs in the two compared groups. The three genes, XLOC_034373 (NR databse_id: O24301.1), XLOC_025267 (NR databse_id: CAC32462.1), and XLOC_002091 (NR databse_id: XP_003591492.2), were more highly expressed in the vegetable pea pericarp at development stages I and V, respectively. The expression level of these three genes in the vegetable pea pericarp were 1.81-, 3.12-, and 1.66-fold higher than that of the grain pea pericarp (Table S3, Figure 5B, Supplementary Dataset 1).

DEGs Related to Pod Sucrose Metabolism
Sucrose synthase (SUS) and invertase (INV) play an important role in sucrose metabolism. All of the sugar precursors that participate in sucrose metabolism must be decomposed to hexoses, such as glucose and fructose, or the hexoses (e.g., UDP-glucose) must be ramified by SUS or INV [29]. A total of three encoding sucrose synthase genes were deemed to be the significant DEGs in the two compared groups. The three genes, XLOC_034373 (NR databse_id: O24301.1), XLOC_025267 (NR databse_id: CAC32462.1), and XLOC_002091 (NR databse_id: XP_003591492.2), were more highly expressed in the vegetable pea pericarp at development stages I and V, respectively. The expression level of these three genes in the vegetable pea pericarp were 1.81-, 3.12-, and 1.66-fold higher than that of the grain pea pericarp (Table S3, Figure 5B, Supplementary Dataset 1).
In addition, genes of the soluble sugar metabolism pathway were analyzed. These include fructose-bisphosphate aldolase (FBA)-and triosephosphate isomerase (TPI)-encoding genes, which participate in fructose and mannose metabolism (ko00051), and the uridylyl transferase-encoding gene, which participates in galactose metabolism (ko00052). Our results showed that five fructose-bisphosphate aldolase-encoding genes, XLOC_028238 (NR database_id: P46257.1), XLOC_023185 (NR database_id: P46256.1), XLOC_003731(NR database_id: Q43088.1), XLOC_027726, and XLOC_033268 (NR database_id: XP_003607065.1), were differentially expressed in the grain pea pericarp and vegetable pea pericarp at various development stages. Among them, one gene, P46257.1 (Transcript_id: XLOC_028238), was more expressed in the vegetable pea pericarp from development stage I to stage V. With plant growth, the expression level of XLOC_028238 decreased form stage I to II, and then increased steadily and reached the second highest at stage III, and then decreased (Supplementary Dataset 1). Similarly, a varying trend of gene XLOC_028238 expression occurred in grain pea, although its FPKM value was lower than that in the vegetable pea pericarp (Table S3, Figure 5B). The other four genes were highly expressed in vegetable pea pericarp at all stages, except for XLOC_027726, which was repressed in stage V. Moreover, one triosephosphate isomerase-encoding gene, XLOC_011807 (NR database_id: XP_013465404.1), was identified, and it was found that this gene was more expressed in the vegetable pea pod at stages I, II, and V but was slightly higher expressed at stages III and IV than the grain pea pod. A total of six uridylyl transferase-encoding genes were detected, but these genes exhibited no great difference in their expression levels in the two compared cultivars.
Three starch branching enzyme-encoding genes, XLOC_032511, XLOC_039401, and XLOC_024258 (NR_id: Q41058.1), were all highly expressed in the grain pea pericarp and had a low expression in the vegetable pea pericarp from stage I to V. The expression levels of these genes in the grain pea pericarp were 2.80-to 5.14-fold higher than those in the vegetable pea pericarp (Supplementary Dataset 1). In addition, the GBSS-and SSS-encoding genes were also more highly expressed in the grain pea pericarp at stages III and IV, respectively (Table S3, Figure 5B, Supplementary Dataset 1).

Discussion
The pea pericarp is an important determinant of the morphological characteristics and seed yield. Mature pericarp usually consists of three parts: The exocarp, mesocarp, and endocarp. As the homologous organ of the leaves, the pod pericarp has a stronger photosynthetic capacity than the leaf at the late stage of seed development, and it can continuously input nutrient products into seeds [32]. The pod has been demonstrated to have a complete functional photosynthetic system, and its contribution to seed yield cannot be ignored [33]. Currently, most studies on the pod have focused on soybeans, beans, and other legumes [34,35]. By contrast, reports on the green pea pericarp are rare as it is generally considered to be a biological waste material. Thus, the pea pericarp remains largely uncharacterized. In our research, we selected two domestic pea cultivars, the vegetable pea cultivar WDZY-14, and the grain pea cultivar WDZY-04, which have significant phenotypic characteristics, to investigate the molecular mechanism of different phenotypic traits. The whole transcriptome data of the two cultivars of pea pericarp at five developmental stages were obtained by high-throughput sequencing technology. About 7.98 Gb of clean reads for each sample were filtered, and approximately 87.78% of the sequences were successfully mapped to the pea reference genome ( Table 1). The analysis of the differentially expressed genes reveals that 1935 DEGs co-existed in the five developmental stages. GO and KEGG annotation revealed that these DEGs were mainly involved in the 'metabolic process', 'cellular process', 'single-organism process', 'membrane part', 'membrane', 'organelle', 'binding', 'catalytic activity', 'transporter activity', 'carbohydrate metabolism', 'biosynthesis of other secondary metabolites', 'signal transduction', 'energy metabolism', 'folding, sorting, and degradation', 'amino acid metabolism', and 'translation'.
Reactive oxygen species (ROS) can impact metabolism and plant growth by interacting with proteins, carbohydrates, and other components in the cell and play an important role in cell wall loosening and elongation [36]. Research into the pericarp elongation of Pisum sativum suggested that high levels of O 2 − and ·OH may have an impact on cell wall loosening and cell elongation, and that superoxide dismutase (SOD) and peroxidase (POD) were associated with pericarp growth through the regulation of ROS generation and transformation [28]. A high SOD activity and low-level POD can increase pod wall thickness by regulating ROS [12]. In our study, we identified three SOD-encoding genes that were differentially expressed in the two cultivars. More SOD-encoding genes were expressed in the grain pea pericarp at development stages III and IV. According to our measurements, the size of the grain pea WDZY-04 pericarp experienced the greatest change from stage I to stage IV, especially in the early developmental stages. Thus, we suspect that more SOD-encoding genes in the pod of grain pea contributed to the faster pod elongation. In addition, five POD-encoding genes were lower expressed in the vegetable pea pericarp at the developmental stage. A previous study indicated a positive relation between pod wall thickness and SOD activity and a negative relation between pod wall thickness and POD activity, and these two enzymes are synergic and responsible for pod growth through the regulation of ROS generation [12]. Similar to the findings of Liu et al. [12], our Pearson's correlation coefficient analysis also showed a significant positive correlation between the expression level of one SOD-encoding gene and pod elongation at WDZY-04 and a significant negative correlation between the expression level of three POD-encoding genes at WDZY-04 and WDZY-14. Thus, we think that the more POD-encoding genes found in the grain pea pod is one of the reasons for the vegetable pea pod being longer than the grain pea pod. These results further suggest that POD may have contributed to pod elongation.
Vegetable pea typically has a higher content of sugar than grain pea, so it is sweeter than grain pea [25]. This remarkable phenotype characteristic has led to the division of the food-oriented pea into two types (vegetable pea and grain pea). Pea pericarp consists of wrapped seeds that are directly linked to each other, so it plays an important role in both pea seed yield and quality. The reason for peas being sweeter is closely related to sucrose metabolism in plant cells. The distribution of sucrose to storage organs, such as seeds, fruits, and tubers, is one of the most important factors determining crop yield and quality, so sucrose metabolism also plays a key role in plant growth and development [37,38]. After unloading the phloem into sinks, sucrose must be degraded into hexoses or their derivates to become metabolically available [39]. Sucrose synthase (SUS) and invertase (INV) play an important role in this process. Sucrose synthase is a glycosyltransferase, and can convert sucrose into UDP-glucose and fructose in the presence of UDP [40,41]. Invertase is a hydrolytic enzyme that hydrolyzes sucrose into glucose and fructose [42]. In our study, we identified three sucrose synthase-encoding genes that were more highly expressed in the vegetable pea pericarp than in grain pea, and hypothesized that these three encoding genes are involved in sugar synthesis in pea pericarp as one of the molecular mechanisms that contribute to the sweet taste of peas. Meanwhile, we also identified four invertase-encoding genes that were more highly expressed in the vegetable pea pericarp, and suggested that these genes may be one of the factors responsible for the higher sweetness of the vegetable pea pericarp. Furthermore, the other genes related to sucrose metabolism were also analyzed. Uridylyl transferase, fructose-bisphosphate aldolase, and triosephosphate isomerase were all more highly expressed in the vegetable pea pericarp than in the grain pea cultivar. These results suggested that these genes play a role in the sweet phenotype. At the same, the identified expression of these genes, associated with sucrose metabolism in the pea pericarp, was consistent with previous studies on pea seeds [25]. We also hypothesized that the traits of the pea pericarp could also affect the sweet traits of pea seeds.
Similarly, considering the difference in starch content between the two peas, we identified the DEGs related to starch metabolism in the pea pericarp to investigate the starch synthesis mechanism in the pea pericarp and the relation between pea and the pea pericarp. Starch is the main carbohydrate that accumulates in mature seeds. Genes involved in starch biosynthesis have been reported in a previous study, and include SUS, GBSS, SS, BE, ADPase, and DBE [31]. However, pea pericarp and pea starch biosynthesis remain largely unknown. We characterized the expression levels of ADPase, GBSS, the starch branching enzyme phosphoglucomutase, and soluble starch synthase-encoding genes, and the results revealed that most of these genes related to starch synthesis were highly expressed in the grain pea pericarp, especially at the later stage of plant growth. At the later stage of pea seed development, the grain pea seed accumulated a lot of starch and its starch content was higher than that of the vegetable pea. Interestingly, we observed that the starch-related genes in the grain pea pericarp were more highly expressed at a later development stage compared to vegetable pea, and we also presumed that the mechanism of pea pericarp starch metabolism may affect pea seed starch accumulation.
In short, our study determined the phenotype characteristics of two pea pericarps, and a transcriptomic analysis was conducted. The aim was to research the molecular mechanisms involved in different traits of the vegetable pea pericarp and grain pea pericarp. Numerous DEGs were identified as being involved in ROS generation and sucrose and starch biosynthesis in the pea pericarp. Taken together, these results may facilitate the understanding of the molecular mechanisms involved in the sweetness and growth differences between the two types of peas, as well as aid in the construction of a genetic map for pea.

Plant Material
The domesticated pea (Pisum sativum L.) samples, including the vegetable pea cultivar "WDZY-14" and the grain pea cultivar "WDZY-04", were planted in an experimental farm at Northwest A&F University, Yangling (E108 • 4 /N34 • 16 /511 m), Shaanxi, China, in the 2018 growing season. Three biological replicates were set up for each cultivar. The pea is a self-pollination species, and the first day of the flower being fully open is defined as the first day of fertilization (0 DAP) [43]. After pollination on 17 April (the flower fully open, 0 DAP), we collected pea pericarp at five developmental stages according to the pea growth of the two cultivars for further study: Stage I (7 days after pollination, 7 DAP), stage II (14 days after pollination, 14 DAP), stage III (21 days after pollination, 21 DAP), stage IV (28 days after pollination, 28 DAP), and stage V (35 days after pollination, 35 DAP). Each pericarp was considered one biological replicate, and we collected three pericarps in three plants at five developmental stages. Thus, 30 samples were collected in total (two cultivars at five developmental stages, with three replicates, e.g., the number, WDZY-04-I#1, was the cultivar, WDZY-04, from the first replicate at stage I). Figure 1A shows the different developmental stages of the pericarp in the two pea cultivars. After harvest, the tissues were immediately frozen in liquid nitrogen and stored at −80 • C until required. For RNA extraction, we divided the pericarp and seed of the pea of each cultivar.

Measurement of the Phenotypic Traits of the Two Cultivars
An estimation of the sugars and starches was carried out on the whole pericarps of the two cultivars harvested at each developmental stage according to a previous methodology [44]. A total of 0.1 g (dry weight) of the samples were ground into fine powder and placed into a 10-mL centrifuge tube, and then homogenized in 4.0 mL of ethanol (80%). After extraction in a water bath (80 • C) for 40 min and centrifugation, the supernatant was collected, and 4.0 mL of 80% ethanol were added to the precipitate, for another extraction. After decolorization in an 80 • C water bath for 30 min, the supernatants were prepared to a constant volume of 25 mL. After filtration, the filtrate was taken for analysis. The soluble total sugar content and starch content were determined by the anthrone method, and the sucrose content was determined by the resorcinol method.

Illumina Sequencing and Mapping
The total pea pericarp RNA was extracted using an RNAprep Pure Plant kit (Bio TeKe, Beijing, China) following the manufacturer's instructions. The purity and concentration of each RNA sample was determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The equal amounts of RNA from each sample were pooled for cDNA library construction. Stranded cDNA libraries were constructed using a NEBNext Ultra Directional RNA Library Prep Kit (cat#E7420, NEB, UK) according to the manufacturer's protocols. Briefly, the mRNA was fragmented into 250 to 450 bp followed by first strand cDNA synthesizing. Then, dUTP was added as a marker during the syntheses of the second strand cDNA. Finally, the double strand cDNA was digested with Uracil -DNA Glycocasylase (UDG) before PCR reaction. In this way, only the first strands of cDNA were kept and sequenced. Transcriptome sequencing was carried out on the Hiseq 4000 (Illumina, San Diego, CA, USA) platform using a paired-end run (2 × 150bp).
After removing the adaptor sequences and low-quality sequences (Q < 30), the clean reads were aligned to the reference genome sequences of the Pisum sativum Linn genome (downloaded from NCBI: https://www.ncbi.nlm.nih.gov/nuccore/PUCA000000000.1/) using Tophat2 software with default parameters [45]. All of the raw data, generated by sequencing, were deposited in NCBI SRA under the following accession: PRJNA548433.

Expression Level Analysis and Gene Annotation
A fragments per kilobase per million (FPKM) analysis, which simultaneously considers the sequencing depth and length of a count, was used to measure the gene expression levels [46]. Genes with an expression level of at least 1 FPKM in at least one sample were retained after removing genes with low expression levels. The differentially expressed genes (DEGs) were identified by deseq2 in R software, with an FDR (false discovery rate) ≤ 0.05 and |log2 (fold change)| ≥1 between two samples.
Gene ontology (GO) enrichment analysis and normalized gene expression data were used to identify the function of and relationships between differentially expressed genes. The identified DEGs were subjected to GO (Gene Ontology: http://geneontology.org/) and KEGG pathway enrichment analysis using phyper in R software. Moreover, the NR (NCBI nonredundant protein sequences (https: //www.ncbi.nlm.nih.gov/refseq/about/nonredundantproteins/) and GO and KO (KEGG ORTHOLOGY, https://www.kegg.jp/kegg/ko.html)) annotation of transcripts were carried out using BLAST (cutoff E-value < 1 × 10 −5 ).

qRT-PCR Validation
A total of 13 genes were selected to validate the RNA-seq data. The HiScript ® II Q RT SuperMix for qPCR (+g DNA wiper) Kit (Nanjing, China) was used according to the manufacturer's instructions to generate the first cDNA after extracting the total RNA from six samples subjected to RNA-seq. The gene-specific primers were produced by Primer 5.0 software (Premier, Canada), and the primer sequence is listed in Table S4. The gene-specific primer was synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). ChamQTM SYBR ® Color qPCR Master Mix (10 µL; Vazyme, Nanjing, China) was mixed with gene-specific primers, sterilized water, and the synthesized cDNA, with 20 µL as the total reaction volume. The reaction were performed on an qTOWER 2.2 (Analytik Jena AG, Jena, Germany). The two-step quantitative RT-PCR program began at 95 • C for 30 s, followed by 40 cycles of 95 • C for 10 s and 60 • C for 30 s. Each reaction was carried out with three biological replicates and three technical replicates. Tubulin was used as the internal reference gene. The data were analyzed using the 2 −∆∆Ct method to obtain relative mRNA expression data.

Conclusions
We investigated the molecular mechanism responsible for the phenotypic traits of the pericarp of two pea cultivars through a transcript analysis of five developmental stages. We analyzed the DEGs of the two cultivars and determined their relative expression levels. Then, hierarchical clustering was used to analyze the expression pattern of the DEGs. As a result, we identified 1935 DEGs common to the five developmental stages of pea. Moreover, we identified important genes related to pod elongation and starch and sucrose synthase that may influence the seed quality. Our research will provide a basis for further studies on starch biosynthesis in pea and a reference for heredity breeding.

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