Transcriptome Analysis Reveals Key Seed-Development Genes in Common Buckwheat (Fagopyrum esculentum)

Seed development is an essential and complex process, which is involved in seed size change and various nutrients accumulation, and determines crop yield and quality. Common buckwheat (Fagopyrum esculentum Moench) is a widely cultivated minor crop with excellent economic and nutritional value in temperate zones. However, little is known about the molecular mechanisms of seed development in common buckwheat (Fagopyrum esculentum). In this study, we performed RNA-Seq to investigate the transcriptional dynamics and identify the key genes involved in common buckwheat seed development at three different developmental stages. A total of 4619 differentially expressed genes (DEGs) were identified. Based on the results of Gene Ontology (GO) and KEGG analysis of DEGs, many key genes involved in the seed development, including the Ca2+ signal transduction pathway, the hormone signal transduction pathways, transcription factors (TFs), and starch biosynthesis-related genes, were identified. More importantly, 18 DEGs were identified as the key candidate genes for seed size through homologous query using the known seed size-related genes from different seed plants. Furthermore, 15 DEGs from these identified as the key genes of seed development were selected to confirm the validity of the data by using quantitative real-time PCR (qRT-PCR), and the results show high consistency with the RNA-Seq results. Taken together, our results revealed the underlying molecular mechanisms of common buckwheat seed development and could provide valuable information for further studies, especially for common buckwheat seed improvement.


Introduction
Seed or grain development is one of the most crucial processes in plant development. It not only determines the successful racial continuation of seed plants, but also affects the final seed yield and seed quality [1,2]. Generally, seed development can be broadly divided into two phases: Embryogenesis and maturation [3,4]. The embryogenesis phase mainly involves cell division and expansion, while maturation is usually an accumulation process of seed storage materials [3]. Seed development is predominantly regulated by genetic factors and is also greatly influenced by physiological pathways and environmental cues. At the physiological level, phytohormones-including auxins, cytokinins (CKs), gibberellins (GAs), brassinolides (BRs), abscisic acid (ABA), and ethylene (ET)-contribute to seed development [2]. At the molecular level, in the past two decades, many genes involved in seed development have been identified and their molecular pathways deciphered in the model plants Arabidopsis and rice through analysis of mutants, transcriptomes, and quantitative trait locus (QTLs) [2, [5][6][7][8][9]. Furthermore, in the last few years, investigations of the molecular mechanisms of seed development have also been undertaken in some crop plants such as maize [10,11], wheat [12,13], soybean [14][15][16], barley [4], Tartary buckwheat [17,18], chickpea [19], and Brassica napus [20,21] by transcriptome analysis. These studies have identified large numbers of related regulatory and functional genes, and provided a better understanding of the molecular mechanisms of seed development.
Common buckwheat (Fagopyrum esculentum) is an annual minor crop which belongs to the eudicot family Polygonaceae, mainly cultivated in temperate regions of Asia, Europe, and North America [22]. As a major part of the human diet, common buckwheat seeds contain high levels of starch, protein, dietary fibers, vitamins, minerals, and flavonoids, which are beneficial to human health [23][24][25]. Furthermore, the seed of common buckwheat is larger than that of another cultivated buckwheat, tartary buckwheat (Fagopyrum tataricum). Therefore, it is of considerable interest to identify genes and dissect the molecular mechanisms of seed development in common buckwheat, as this will contribute to common buckwheat and tartary buckwheat seed improvement, including seed formation, morphogenesis, size, and nutrient accumulation. To date, several studies have investigated the transcriptional profile of the developing seeds of common buckwheat though RNA-Seq [26][27][28]. However, these studies only focused on one developmental time point of common buckwheat seed, and the differentially expressed genes that might be involved in seed development have not been identified.
In the present study, we carried out a transcriptomic analysis of common buckwheat seed at three developmental stages (pre-filling stage (S1), filling stage (S2), and initial maturity stage (S3)) with the purpose of investigating its transcriptional dynamics and identifying candidate genes involved in seed development. In addition, some vital seed-development candidate genes were verified by quantitative real-time PCR (qRT-PCR). The results obtained in this study will enhance our understanding of the molecular mechanisms of common buckwheat seed development and provide candidate genes for future common buckwheat seed improvement.

Overview of Sequencing Data Analysis
To gain insight into the transcriptional dynamics of seed development in common buckwheat, nine libraries from three samples (three biological replicates for each sample) were constructed and subjected to high-throughput RNA-Seq. An overview of the sequencing data is shown in Table 1. A total of 248.53 million raw reads, with an average of 27.61 million raw reads for each library, were yielded. After removing the adaptor sequences and low-quality sequence reads, 242.86 million clean reads were obtained, with an average of 26.98 million clean reads for each library. The Q30 percentages for all libraries exceeded 93%, and the range of GC (Guanine and Cytosine) content for each library was 47.65 to 50.58%. A total of 63.32 to 69.46% of clean reads were mapped to the reference genome for each library. The Pearson's rank correlation between any two of the three replicates for each sample was >0.95, with 0.979, 0.985, and 0.972 for S1, 0.952, 0.964, and 0.968 for S2, and 0.983, 0.993, and 0.998 for S3 (Table S1), indicating the high reliability and reproducibility of the data. Item S1-1 S1-2 S1-3 S2-1  S2-2  S2-3  S3-1  S3-2  S3-3   Raw Reads  25,672,949 22,929,589 24,865,374 25,424,457 27,829,981 30,527,439 25,646,270 37,234,140 28,939,774  Clean Reads  25,494,672 22,747,870 24,566,718 25,250,033 27,625,287 30,212,779 21,406,245 36,901, Notes: S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, and -1, -2, and -3 represent the three replicates per sample. GC Content (%) represents the percentage of Guanine and Cytosine in clean reads. Q_30 (%) represents the percentage of nucleotides with quality value ≥30.

Identification and Annotation of Differentially Expressed Genes (DEGs)
To identify the DEGs during the common buckwheat seed development, the gene expression profile clustering was subjected to first analysis. All 54,582 genes were assigned to seven different clusters by the K-means method, which led to the identification of genes with similar expression during all stages of seed development, and indicated that a large number of genes were differentially expressed between different samples ( Figure 1).  Item S1-1 S1-2 S1-3 S2-1  S2-2  S2-3  S3-1  S3-2  S3-3  Raw  Reads  25,672,949  22,929,589  24,865,374  25,424,457  27,829,981  30,527,439  25,646,270  37,234,140  28,939,774   Clean  Reads  25,494,672  22,747,870  24,566,718  25,250,033  27,625,287  30,212,779  21,406,245  36,901, Notes: S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, and -1, -2, and -3 represent the three replicates per sample. GC Content (%) represents the percentage of Guanine and Cytosine in clean reads. Q_30 (%) represents the percentage of nucleotides with quality value ≥30.

Identification and Annotation of Differentially Expressed Genes (DEGs)
To identify the DEGs during the common buckwheat seed development, the gene expression profile clustering was subjected to first analysis. All 54,582 genes were assigned to seven different clusters by the K-means method, which led to the identification of genes with similar expression during all stages of seed development, and indicated that a large number of genes were differentially expressed between different samples ( Figure 1).

Figure 1.
Clustered gene expression profiles from developing tartary buckwheat seeds. The clusters were defined on the basis of genes' temporal expression profiles in R using K-means. The Y-axis represents the standardized FPKM value of genes, and the X-axis represents the different sample. S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, respectively. Then, gene expression levels among the three developmental stages were subjected to a pairwise comparison. A total of 4619 DEGs were identified, of which 745 (492 up-and 253 down-regulated), 2449 (1506 up-and 943 down-regulated), and 3680 (2162 up-and 1518 down-regulated) were found between S1 and S2, S2 and S3, and S1 and S3, respectively ( Figure 2B). Of these DEGs, 304, 576, and 1572 DEGs specifically existed in the S1 vs. S2, S2 vs. S3, and S1 vs. S3 comparisons, respectively; 147, 382, and 1814 DEGs existed both in S1 vs. S2 and S2 vs. S3, S1 vs. S2 and S1 vs. S3, and S2 vs. S3 and S1 vs. S3, respectively; and 88 DEGs were contained in all three comparisons (Figure 2A).

Figure 1.
Clustered gene expression profiles from developing tartary buckwheat seeds. The clusters were defined on the basis of genes' temporal expression profiles in R using K-means. The Y-axis represents the standardized FPKM value of genes, and the X-axis represents the different sample. S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, respectively.
Then, gene expression levels among the three developmental stages were subjected to a pairwise comparison. A total of 4619 DEGs were identified, of which 745 (492 up-and 253 down-regulated), 2449 (1506 up-and 943 down-regulated), and 3680 (2162 up-and 1518 down-regulated) were found between S1 and S2, S2 and S3, and S1 and S3, respectively ( Figure 2B). Of these DEGs, 304, 576, and 1572 DEGs specifically existed in the S1 vs. S2, S2 vs. S3, and S1 vs. S3 comparisons, respectively; 147, 382, and 1814 DEGs existed both in S1 vs. S2 and S2 vs. S3, S1 vs. S2 and S1 vs. S3, and S2 vs. S3 and S1 vs. S3, respectively; and 88 DEGs were contained in all three comparisons (Figure 2A).  To further explore the function of these DEGs, Gene Ontology (GO) enrichment analysis was carried out. As a result, 395 (S1 vs. S2), 1244 (S2 vs. S3), and 1878 (S1 vs. S3) DEGs were divided into biological process, cellular component, and molecular function categories ( Figure 3 and Table S2). As shown in Figure 3, the biological process category contained 20 functional groups, in which 'metabolic process', 'cellular process', and 'single-organism process' were the top three represented terms. In the cellular component category, 16 functional groups were enriched, and the most abundant terms were 'cell part', 'cell', and 'organelle'. For the molecular function category, 16 functional groups were identified, and the terms with greatest abundance were 'catalytic activity', 'binding', and 'transporter activity'. To further explore the function of these DEGs, Gene Ontology (GO) enrichment analysis was carried out. As a result, 395 (S1 vs. S2), 1244 (S2 vs. S3), and 1878 (S1 vs. S3) DEGs were divided into biological process, cellular component, and molecular function categories ( Figure 3 and Table S2). As shown in Figure 3, the biological process category contained 20 functional groups, in which 'metabolic process', 'cellular process', and 'single-organism process' were the top three represented terms. In the cellular component category, 16 functional groups were enriched, and the most abundant terms were 'cell part', 'cell', and 'organelle'. For the molecular function category, 16 functional groups were identified, and the terms with greatest abundance were 'catalytic activity', 'binding', and 'transporter activity'. . GO enrichment of identified DEGs in developing common buckwheat seeds. S3. S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, respectively.

DEGs Involved in Hormone Signal Transduction Pathways
Hormone signals play an important role in plant seed development. Based on the KEGG enrichment analysis results, 30 DEGs were assigned to "plant hormone signal transduction" (ko04075), including the auxin, CK, GA, BR, ABA, ET, JA, and SA signaling pathways (Table 3). In the auxin signaling pathway, one gene encoding auxin-responsive protein (IAA) was up-regulated at S2, five IAA genes were up-regulated at S3, and one ortholog of AFR encoding for auxin response factor was up-regulated both at S2 and S3 (Table 3). Moreover, one ortholog of AUX1 encoding for auxin influx carrier and two IAA were down-regulated at S3. In the cytokinin pathway, one ARR (two-component response regulator) and one AHP (histidine-containing phosphotransfer protein) were increased at S2 and S3, respectively (Table 3). In the GA pathway, one PIF3 (phytochrome-interacting factor 3) showed up-regulated trends during seed development, and one ortholog of bHLH127 was increased only at S3 (Table 3). In the ABA pathway, there were two up-and one down-regulated transcripts at S2, which belonged to PP2C (protein phosphatase 2C) and ABF (encoding for ABA responsive element binding factor), respectively. In addition, one PYR/PYL (abscisic acid receptor), one PP2C, and one SnRK2 (serine/threonine protein kinase) were down-regulated at S3 (Table 3). In the ET pathway, two orthologs of EIN3 (Ethylene-Insensitive Protein 3) were up-regulated at S2 and S3, respectively. Furthermore, two ETR (Ethylene Receptor), one EIN4 (Ethylene-Insensitive Protein 4), and one CTR1 (serine/threonine-protein kinase) were down-regulated at S3 (Table 3). In the BR pathways, one BRI1 (Brassinosteroid Insensitive 1) was down-regulated at S3, and other BRI1 displayed up-regulated at S2 and down-regulated at S3 (Table 3). In addition, one DEG was involved in both the JA and SA pathways, which was up-regulated (JAZ) at S3 and down-regulated (TGA) at S2 (Table 3), respectively.

DEGs Involved in TFs
TFs are the key regulatory proteins in seed development. Therefore, the expression dynamics of TF genes in the development of common buckwheat were investigated. According to the RNA-seq data, a total of 2453 TFs were identified as expressed in at least one development phase (Table S4) regulated) TFs had significantly differential expression in S1 vs. S2, S2 vs. S3, and S1 vs. S3, respectively (Table S3). Among these, the TF families with the largest numbers were AP2 (9), bHLH (11), and AP2 (22) for these three comparison groups, respectively (Table S5).

DEGs Involved in Seed Size
Seed size is a key factor to determine seed yield in crops. To identify the genes involved in the seed size of common buckwheat, 95 known seed size-related genes from different plant species were used to homologously search the DEGs (Table S6). As a result, 17 DEGs were identified as the orthologs of 15 known seed size-related genes (AtANT, AtAP2, AtDA1, AtDET2, AtDWF4, AtEOD1, AtIKU2, AtTTG2, OsSLG, OsGRF4, OsGIF1, OsSRS5, OsWRKY53, OsCYP78A13, and OsCYP724B1) (Table S7). Of these, two genes (DWF4 and CYP724B1) and one gene (GIF1) showed up-regulation only at S2 or at S3, respectively; five genes (DET2, GRF4, WRKY53, IKU2, and TTG2) were up-regulated during seed development stages; four genes (SSR5, SSR5, IKU2, and SLG), with high expression at S1 and S2, were down-regulated at S3; two genes (CYP78A13 and ANT) were up-regulated at S1 compared to S2, and followed by down-regulation at S3; however, three genes (AP2, DA1, and EOD1) displayed the opposite expression pattern with CYP78A13 and ANT, which were down-regulated at S1 compared to S2, and followed by up-regulation at S3 ( Figure 5 and Table S7). Table 3. DEGs involved in hormone signaling pathways. S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, respectively. FDR stands for false discovery rate. Log2FC stand for Log2 fold change.  Figure 5. Heat map diagrams of relative expression levels of the differentially expressed orthologs of seed size-related genes in the different development stages of common buckwheat. The heat map was constructed using MeV 4.9.0 based on the Log2(FPKM) value of genes. S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, respectively.

DEGs Involved in Starch Biosynthesis
To provide insight into the transcriptional dynamics of starch biosynthesis genes, a total of 25 sucrose synthase (SUS), 3 UDP glucose pyrophosphorylase (UGPase), 14 ADP glucose pyrophosphorylase (AGPase), 8 granule bound starch synthase (GBSS), 24 starch synthase (SS), 4 starch-branching enzyme (SBE), and 3 debranching enzyme (DBE) genes were identified by using the starch biosynthesis-related genes in Arabidopsis to perform a homologous comparative search in the common buckwheat genome database (Table S8). Of these, four SUS, one UGPase, and one SBE genes showed up-regulation only at S2; one SUS, one AGPase, one GBSS, and one SS displayed downregulation only at S3; two SUSs and one SS, and one SBE gene were up-regulated at S2, and followed by down-regulation at S3; one DBE gene was down-regulated with all seed development stages (Table 4 and Figure 6).

Figure 5.
Heat map diagrams of relative expression levels of the differentially expressed orthologs of seed size-related genes in the different development stages of common buckwheat. The heat map was constructed using MeV 4.9.0 based on the Log 2 (FPKM) value of genes. S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, respectively.

DEGs Involved in Starch Biosynthesis
To provide insight into the transcriptional dynamics of starch biosynthesis genes, a total of 25 sucrose synthase (SUS), 3 UDP glucose pyrophosphorylase (UGPase), 14 ADP glucose pyrophosphorylase (AGPase), 8 granule bound starch synthase (GBSS), 24 starch synthase (SS), 4 starch-branching enzyme (SBE), and 3 debranching enzyme (DBE) genes were identified by using the starch biosynthesis-related genes in Arabidopsis to perform a homologous comparative search in the common buckwheat genome database (Table S8). Of these, four SUS, one UGPase, and one SBE genes showed up-regulation only at S2; one SUS, one AGPase, one GBSS, and one SS displayed down-regulation only at S3; two SUSs and one SS, and one SBE gene were up-regulated at S2, and followed by down-regulation at S3; one DBE gene was down-regulated with all seed development stages (Table 4 and Figure 6). Table 4. DEGs involved in starch biosynthesis. S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, respectively. FDR stands for false discovery rate. Log2FC stand for Log2 fold change.  Debranching enzyme. S1, S2, and S3 stand for the seed samples at prefilling stage, filling stage, and initial maturity stage, respectively.

qRT-PCR Validation of RNA-Seq Results
Based on the results mentioned above, 15 genes which had significantly differential expression during seed development were selected by performing qRT-PCR assays to validate the correlation against RNA-Seq. The results showed that 14 genes had high correlation between the qPCR and RNA-Seq data sets (Figure 7), which thus verified the transcriptomic data. Debranching enzyme. S1, S2, and S3 stand for the seed samples at pre-filling stage, filling stage, and initial maturity stage, respectively.

qRT-PCR Validation of RNA-Seq Results
Based on the results mentioned above, 15 genes which had significantly differential expression during seed development were selected by performing qRT-PCR assays to validate the correlation against RNA-Seq. The results showed that 14 genes had high correlation between the qPCR and RNA-Seq data sets (Figure 7), which thus verified the transcriptomic data.

Discussion
The seed development of crop plants affects not only seed fate, but also final seed yield and quality. Thus, an in-depth understanding of the molecular mechanisms of seed development is crucial in improving seed yield and quality of crop plants. To date, many genes that contribute to seed development have been identified and functionally verified, and molecular pathways have been deciphered in the model plants Arabidopsis and rice [4][5][6][7][8][9]. Furthermore, in some other crop plants, transcriptome analysis has been performed to give insight into the transcriptional dynamics of seed development and provided some clues as to the development of seeds [10][11][12][13][14][15][16][17][18][19][20][21]. In this study, we performed a transcriptome analysis during common buckwheat seed development to give insight into its transcriptional dynamics and find candidate genes that might be involved in seed development.
In total, 242.86 million clean reads were obtained from nine sequencing libraries, and 63.32 to 69.46% clean reads for each library were mapped to the reference genome. The map rate was relatively low, which might have been caused by the draft assembly of the common buckwheat genome [22]. The Q30 percentage was over 93% for each library, and the range of GC content was 47.77-50.58% for each library. Furthermore, the reproducibility was very well (r > 0.95) for the three replicates of each sample. These results indicated that the RNA-Seq had higher quality. Through comparison of gene expression patterns between different samples, 4619 DEGs were identified, and the DEGs involved in the Ca 2+ signaling pathway, hormone signaling pathway, TFs, seed size, and starch biosynthesis were the main focus of this study.
The Ca 2+ signal transduction pathways participate in many developmental processes in plants, including pollen tube growth, stem elongation, vascular development, embryogenesis, and seed germination [29,30]. In addition, some studies have also indicated that the Ca 2+ signal transduction pathways play an important role in seed development with a focus on CDPK. For example, several CDPK genes were differential or special expression during seed development in rice, maize, and peanut, some of which expressed dominant in the early stage, while others expressed prominently in the latter stage [29,[31][32][33][34][35][36], suggesting that different CDPK genes might be involved in different development process of seeds. Furthermore, OsCDPK2, OsCDPK11, and OsCDPK23 in rice and RcCDPK1 in castor have been functionally verified to be essential for seed development [32,[34][35][36]. In this study, four CDPK genes were identified as DEGs during common buckwheat seed development. Among them, two were up-regulated at S2 and down-regulated at S3, and one was only decreased at S3, indicating that they might have an important role in the early stages of common buckwheat seed development. In contrast, the other one was down-regulated at S2, suggesting that it might have an opposed function during common buckwheat seed development with the previous three CDPK genes. In addition, we also identified another 18 DEGs involved in the Ca 2+ signaling pathway, including 3 CML, 2 CBL, 2 CCX, 5 CIPK, 4 Ca 2+ -ATPase, and 2 CHX genes, which displayed similar expression patterns with CDPK genes. These results indicate that the Ca 2+ signaling pathway is conserved in regulatory seed development and that these identified DEGs could be key candidate genes for the seed development of common buckwheat.
It is well known that hormone signaling, including auxin, CK, GA, BR, ABA, and ET signaling, plays a crucial role in seed development [8,9]. In tartary buckwheat (Fagopyrum Tararicum), a dozen of these hormone-related DEGs were identified by RNA-seq, which displayed different expression patterns during seed development [17,18]. In this study, based on KEGG pathway enrichment analysis, we also identified a dozen genes expressed differentially in these hormone signaling pathways during common buckwheat seed development, and also found that some of them have the similar expression patterns with those in tartary buckwheat. Both our study and previous studies have suggested that these hormone signaling pathways are conserved for seed development in different seed plants. Notably, we found one JAZ and one TAG gene, which are the key genes in the JA and SA signaling pathways, respectively, were also differentially expressed during common buckwheat seed development, implying that the JA and SA signaling pathways probably also have important roles in common buckwheat seed development.
Seed size is a key factor to determine yield in crops. The identification of genes associated with seed size will support yield improvement. To date, over 90 genes have been identified and functionally validated to regulate seed size in different seed plants [9]. In our study, 17 DEGs were identified as the orthologs of 15 known seed size genes from Arabidopsis and rice. Among these, the expression level of the orthologs of AtAP2, AtDA1, AtEOD1, and OsSLG were down-regulated during seed development. Interestingly, AtAP2 [37,38], AtDA1 [39], AtEOD1 [39,40], and OsSLG [41] have been demonstrated to negatively regulate the seed size in Arabidopsis and rice, respectively. This indicates that these genes from common buckwheat may also negatively regulate the seed size. In contrast, the orthologs of several positively-regulated genes for seed size, such as AtIKU2 [42], OsGRF4 [43,44], OsGIF1 [43,44], OsCYP78A13 [45,46], OsCYP724B1 [47], AtDET2 [48], AtDWF4 [48], OsWRKY53 [49], AtANT [50], and AtTTG2 [51], were increased at express level in at least one development stage, suggesting that they have a conservative function of positively regulating the seed size of common buckwheat. Notably, the expression level of the orthologs of OsCYP78A13 and AtANT were up-regulated at filling stage and down-regulated at initial maturity stage, indicating that they might regulate the accumulation of seed storage materials and ultimately affect seed size. In addition, an ortholog of AtIKU2 and two orthologs of OsSRS5 (positively regulated seed size) were down-regulated at both filling stage and initial maturity stage, suggesting that they might preform function for seed size at the embryogenesis phase, mainly related to cell division and expansion.
The seeds of common buckwheat contain abundant starch, accounting for over 70% of seed dry weight [23][24][25]. However, the mechanism of starch biosynthesis in common buckwheat remains unclear. It is well known that seven class genes, including SUS, UGPase, AGPase, GBSS, SS, SBE, and DBE, catalyze starch biosynthesis [52]. In our data, we identified 25 SUS, 3 UGPase, 14 AGPase, 8 GBSS, 24 SS, 4 SBE, and 3 DBE genes. Of these, seven SUS, one AGPase, one UGPase, one GBSS, two SS, two SBE, and one DBE genes were differentially expressed during the seed development of common buckwheat. This finding suggests that these identified starch biosynthesis DEGs might be the major functional genes for starch biosynthesis in common buckwheat seed, and could accelerate further functional analysis of them as well.
Seed development is predominantly regulated by genetic factors and is also greatly influenced by physiological pathways [1,2]. Usually, seed development is involved in seed size change and various nutrients accumulation, which determine crop yield and quality [2]. In this study, we provided insight into the transcriptional dynamics of seed development in common buckwheat at three different stages, and identified DEGs involved in the development of common buckwheat seed. Of these DEGs, some were prominently expressed in the early stage, and some other were highly expressed in the latter stage, indicating different DEGs play different roles during common buckwheat seed development. More importantly, of these DEGs, some genes were identified as the key candidate genes that might participate in the seed development of common buckwheat, including Ca 2+ signaling, hormone signaling, TFs, seed size, and starch biosynthesis-related genes, and its differential expression patterns during common buckwheat seed development were verified by qRT-PCR. The results are helpful to further dissect the molecular mechanism of buckwheat seed development and provide a basis for buckwheat seed size and quality improvement.

Plant Material and Sample Collection
The common buckwheat cultivar "Chitian No. 1" was used in this study. It was planted in the growth chamber of the Research Center of Buckwheat Industry Technology, Guizhou Normal University (Lat. 26•62 N, 106•72 E, Alt. 1168 m), China, in Spring 2018. All plants were grown under natural daylight and subjected to normal management during the growth periods. Seeds were harvested on the 8th, 14th, and 21st days after full bloom, which corresponded to the pre-filling stage (S1), filling stage (S2), and initial maturity stage (S3). For each stage, approximately 200 mixed seeds were collected from the same five plants, and three independent biological replicates were performed. All samples were immediately frozen in liquid nitrogen and stored at −80 • C for RNA-Seq and qRT-PCR validation analyses.

RNA Extraction, Library Construction, and Sequencing
Total RNA was extracted using the EASYspin Plus Plant RNA Kit (Aidlab, Beijing, China) with three replicates for each sample, according to the manufacturer's instructions. Contaminated DNA in the total RNA samples was removed using DNase I (TAKARA, Dalian, China). The quality and concentration of the total RNA were determined using 1.2% agarose gel electrophoresis and a NanoDrop 2000c spectrophotometer (NanoDrop, Wilmington, DE, USA). Then, the cDNA libraries for Illumina sequencing were constructed using the NEBNext ® Ultra™ II RNA Library Prep Kit for Illumina ® (New England BioLabs, Ipswich, MA, USA). Finally, the constructed cDNA libraries were sequenced using the Illumina XtenPE150 system by Biomarker Technologies Co., Ltd. (Beijing, China).

Identification of DEGs
The transcript abundance of each unigene was calculated and normalized to fragments per kilobase of transcript per killion fragments mapped (FPKM) [62]. Significantly differential gene expression among the three samples was evaluated using the DEseq package [63]. Genes with a threshold of FDR (false discovery rate) value <0.05 and |log2(fold change)| ≥ 1 were assigned as differentially expressed genes (DEGs). The identified DEGs were further subjected to analysis through Gene Ontology (GO) functions, COG (Cluster of Orthologous Groups of proteins), and KEGG pathway analysis.

Validation of DEGs by qRT-PCR
qRT-PCR analysis was performed to verify the RNA-Seq data. These genes that reflected some of the functional categories and different expression levels were selected. The reference cDNA sequences of these genes were obtained from the common buckwheat genome sequence database. The RT-qPCR primers were designed using Primer 5.0 software according to the reference cDNA sequences (Table S9). FeActin7 was used as the internal control. RT-qPCR was conducted using the SYBR premix Ex Taq kit (TaKaRa, Dalian, China) on ABI 7500 Fast Real-time PCR system (ThermoFisher Scientific, Waltham, MA, USA). Each sample was analyzed in triplicate. The relative expression change of each gene was calculated using the 2 −∆∆Ct method [64]. The Pearson correlation coefficients between the fold change among different stages from qRT-PCR and from RNA-Seq were calculated using Origin 8.0.

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