Phenotypic and physiological changes
Compared with normal tartary buckwheat cultivar “JQ”, the matured hull of rice-tartary buckwheat “XMQ” showed markedly reduced mechanical strength, which was easily broken and shedding when rubbed by hand (Fig. 1a). Scanning electron microscopy showed that many flaws occurred on the “XMQ” hull that did not exist on the “JQ” (Fig. 1b, c). These results suggested that the lower mechanical strength of the “XMQ” hull contributed to its easy dehulling, and the easy dehulling trait was redefined as “brittle hull” trait. To further investigate whether the lower mechanical strength of the “XMQ” hull was caused by the change of hull compositions, we measured the hull cellulose, hemicellulose, and lignin contents of “XMQ”, “JQ”, and other three normal tartary buckwheat cultivars (“JJQ”, “PQ” and “CQ”). As shown in Fig. 2, when compared with the four normal tartary buckwheat cultivars, the hull cellulose content of “XMQ” was significantly lower than all tested normal tartary buckwheat cultivars, while the hull hemicellulose and lignin contents of “XMQ” were only significantly lower than among two normal tartary buckwheat cultivars. These results suggested that the reduced cellulose, hemicellulose, and lignin contents might contribute to the brittleness and easy dehulling of rice-tartary buckwheat, of which the reduced cellulose content might be the major reason.
Dynamic change of cellulose and hemicellulose contents, and hull thickness in “XMQ” and “JQ”
To investigate the accumulation of cellulose and hemicellulose during hull development, the contents of cellulose and hemicellulose in seed hulls of “XMQ” and “JQ” were determined at four different stages (Fig. 3a). As shown in Fig. 3b, c, the two cultivars displayed a similar trend in the change of cellulose and hemicellulose contents, which were sustained growth during hull development. In addition, the significant difference in cellulose and hemicellulose contents between the two cultivars started to occur at 10 DAP. To further insight into whether the hull ditch thickness changes in the two cultivars during seed development, paraffin section analysis was performed and the ditch thickness of the hull was measured according to paraffin section results (Fig. 4). The hull ditch thickness of “XMQ” showed significantly decreased during seed development, while the hull ditch thickness of “JQ” only decreased at 10 DAP compared with 5 DAP and no obvious difference was observed after 10 DAP (Fig. 4b). In addition, the hull ditch thickness of “XMQ” was significantly lower than that of “JQ” at 15 and 20 DAP, although its hull ditch thickness was obviously higher than that of “JQ” at 5 and 10 DAP (Fig. 4b). Transmission electron microscopy further revealed the cell wall thickness of sclerenchyma cells of “XMQ” hull ditch was significantly thinner than “JQ” (Fig. 5).These results indicated that the cellulose content, hemicellulose content, and hull ditch thickness of “XMQ” and “JQ” were dynamic and different changes, and these hull samples could be used for further transcriptome analyses.
Transcriptome analysis of hull development in tartary buckwheat cultivars
To insight into the transcriptome dynamics during hull development, RNA-Seq analyses of the hull of “XMQ” and “JQ” hull were performed at four different development stages (8 tissues, three independent biological replicates for each tissue, 24 samples in total). A total of 521.39 and 517.72 million clean reads, with an average of 43.43 and 43.09 million clean reads for each sample, were generated for “XMQ” and “JQ”, respectively (Additional file 1: Table S1). 85.84 to 92.74 % of clean reads were mapped to the reference genome for each sample, of which 64.54–73.47 % were uniquely mapped. A total of 29,843 gene loci, including 28,465 known and 1,378 novel gene loci, were generated from all mapped clean reads of 24 samples via Cufflinks and Cuffmerge analyses. PCC between the biological replicates of different tissues changed from 0.90 to 0.99 (except for one replicate of J15 and J20 in “JQ” and one replicate of X15 in “XMQ”, which were filtered for further analyses), indicating the high reliability and reproducibility of the replicates (Additional file 2: Figure S1).
The normalized expression level (FPKM) for each gene was calculated, and genes with FPKM ≥ 0.1 were considered as expressed. A total of 27,955 genes, including 26,589 known and 1,366 novel genes, were expressed in at least one sample. Among them, the number of expressed genes in different tissues changed from 23,575 (X20) to 25,014 (X10) in “XMQ”, and 23,282 (J20) to 24,266 (J5) in “JQ” (Additional file 2: Figure S2a). The proportions of genes presenting very high (FPKM ≥ 100), high (50 ≤ FPKM < 100), moderate (10 ≤ FPKM < 50), low (2 ≤ FPKM < 10), and very low (0.1 ≤ FPKM < 2) expression were relatively similar in all tissues except in X20 tissue (Additional file 2: Figure S2b). In addition, the number of very-low-expression (about 42%) and moderate-expression genes (about 30%) accounted for the two largest proportions in all tissues (Additional file 2: Figure S2b). Taken together, these analyses indicated that we obtained sufficient coverage of the transcriptome of the development hull in these two tartary buckwheat cultivars.
Transcriptome comparison of hull development in tartary buckwheat cultivars revealed the critical stages for hull development difference
To investigate the relationships between the hull development transcriptomes from the “XMQ” and “JQ”, HCA (Additional file 2: Figure S3a) and PCA (Additional file 2: Figure S3b) were performed based on the average FPKM values for the 27,955 genes that were expressed in at least one of the eight tissue samples. As results, a higher correlation of the same development stage between the two cultivars was observed except at 20 DAP (Figure S3A). PCA analysis showed X5 and J5 were grouped together, while clear separations were observed between X10 and J10, X15 and J15, and X20 and J20 (Additional file 2: Figure S3b), suggesting that a higher degree of similarity in transcriptional programs between these two cultivars at 5 DAP and obvious transcriptional differences existed at 10, 15, and 20 DAP. In addition, these analyses also indicated that 10 DAP might be the critical stage to determine the hull difference formation between the two cultivars at the molecular level.
Identification of genes that specifically/preferentially expressed in each stage of hull development in both the cultivars
The genes that are specifically/preferentially expressed in each stage of hull development in both cultivars were identified based on the SS algorithm with SS score ≥0.5. A total of 6,649 and 5,918 specific/preferential genes were identified in all four stages for “JQ” and “XMQ”, respectively. Among these genes, 348 (5.2%) and 302 (5.1%) specific/preferential genes were encoded for transcription factors (TFs). The number of stage-specific/preferential genes ranged from 392 to 4613 for “JQ” and 267 to 3959 for “XMQ” (Fig. 6a). The 5 DAP and 10 DAP had the largest and lowest of number stage-specific/preferential genes, respectively, for the two cultivars (Fig. 6a). Furthermore, a high proportion of stage-specific/preferential genes was common in both cultivars, and cultivar-specific genes were also observed in all four stages of the two cultivars (Fig. 6a). A heatmap of all these stage-specific/preferentially expressed genes in the two tartary buckwheat cultivars was shown in Fig. 6b. These results indicated that each hull development stage had its own independent and common developmental programs for both cultivars and also accurately reflected the accumulation of endogenous mRNAs in the hull development of the two cultivars.
The gene ontology (GO) enrichment analyses were carried out for the stage-specific/preferential genes of the two cultivars at each stage. The results of biological process enrichment showed that the 5 DAP of hull development was marked by cell cycle and cell biosynthetic process (Additional file 2: Figure S4), the 10 DAP was significantly and specifically marked by cell wall and phenylpropanoid metabolic process (Fig. 6c), the 15 DAP was majorly marked by “sulfur metabolic” (Additional file 2: Figure S5), and the 20 DAP major was involved in lipid metabolic, response biotic and abiotic stress and transport (Additional file 2: Figure S6). These results suggested that these stage-specific/preferential genes performed stage-specific functions during hull development of tartary buckwheat and 10 DAP was the key stage in the determination cell wall biogenesis of the tartary buckwheat hull.
DGEs between the two tartary buckwheat cultivars
Genes that had significantly different expression between “XMQ” and “JQ” were identified at each stage of hull development. A total of 9,250 (including 693 TF-encoding genes) and 4,187 (294 TF-encoding genes) genes showed significantly higher and lower expression at different stages of hull development in “XMQ” compared with “JQ”, respectively (Fig. 7). Between the two cultivars, the two largest number of DEGs occurred at 20 DAP (5,884) and 10 DAP (3,916) (Fig. 7), while the fewest occurred at 15 DAP. In addition, some members of most TF families shown significantly different expression in “XMQ”, and the number of different TF families showed an obvious difference, which implied they had diverse functions during hull development (Additional file 2: Figure S7).
To further investigate the biological functions of these DEGs between “XMQ” and “JQ”, GO enrichment analyses were performed. A total of 159, 250, 92 and 161 biological processes were significantly enriched (P <0.05) at 5 DAP, 10 DAP, 15 DAP, and 20 DAP, respectively (Additional file 1: Table S2). Among them, some biological processes were specially/commonly overrepresented at different stages of hull development (Additional file 1: Table S2). Notably, GO terms such as the plant-type primary cell wall biogenesis (GO:0009833), cellulose metabolic process (GO:0030243), cellulose biosynthetic process (GO:0030244), and cellulose catabolic process (GO:0030245) were uniquely enriched at 5 DAP, whereas the hemicellulose metabolic process (GO:0010410) and cell wall modification (GO:0042545) were specially enriched at 10 DAP (Additional file 1: Table S2). In addition, cell wall biogenesis (GO:0042546) and fruit dehiscence (GO:0010047) were significantly enriched at 5 DAP and 10 DAP (Additional file 1: Table S2). These results indicate that 5 DAP and 10 DAP were the key stages that determined the difference in the hull between the two tartary buckwheat cultivars.
Identification and valuation of the key genes involved in the hull difference between the two tartary buckwheat bultivars by gene coexpression analysis
To investigate the gene regulatory network (GRN) during hull development and identify the key genes involved in the hull difference between “XMQ” and “JQ”, 9,549 genes that were differentially expressed in at least one hull development stage between the two cultivars were used to carry out WGCNA. A total of 18 modules (consisting of 34-3317 genes) were identified and labeled with different colors (Fig. 8a). The correlation analyses of these modules with the cellulose and hemicellulose contents revealed that the MEbrown module presented higher correlation with cellulose content (r = 0.8, p = 0.02) (Additional file 2: Figure S8). In this module, genes had specific high expression in the “JQ” hull at 20 DAP. No gene in this module was found to be involved in cellulose and cell wall biosynthesis based on homologous annotation in the Arabidopsis TAIR database, which might be caused by the visible hull difference between the two cultivars having formed at 20 DAP. Considering that (1) easy dehulling is a recessive trait, (2) hard dehulling is the corresponding dominant trait, and (3) 10 DAP was identified as the key stage for hull difference formation between the two cultivars, we assumed that the module with genes that presented specific high expression in the “JQ” at 10 DAP was the key module involved in the different hull formation of the two cultivars. As a result, the MEred module was found to meet these standards, which consisted of 533 genes (Fig. 8b). Based on a homologous annotation in the Arabidopsis TAIR database, 28 TFs were identified in this module, which represented 17 TF families. All 28 TFs were further identified as hub TFs due to the highly connected nodes in this module. These hub TFs included homologs of A. thaliana TFs with known functions of regulating SCW biosynthesis, such as the first-layer NAC regulators (NST1, NST2, and NST3/SND1), the second-layer regulator (MYB46/MYB83), and the third-layer regulators (MYB54, MYB103, C3H14, and C3H15) (Additional file 1: Table S3). The homologs of XND1 and VNI2, two NAC TFs that regulate xylem vessel formation, were also identified in these hub TFs (Fig. 9, Additional file 1: Table S4). In addition, several other hub TFs were found to be homologous with A. thaliana TFs, which were involved in ethylene signaling (EIN3, ERF71, and RAV), gibberellin biosynthesis (ATH1), jasmonic acid signaling (WRKY50), and multiple hormonal responses (MIF2). Notably, 6, 12, 1, and 5 SCW biosynthesis enzyme genes were also identified in this module (Fig. 9, Additional file 1: Table S3), which were homologs of known downstream targets of A. thaliana SCW biosynthesis TFs (NAC and MYB) and catalyzed cellulose biosynthesis, hemicellulose biosynthesis, pectin biosynthesis, and lignin biosynthesis, respectively. Among these enzyme genes, four cellulose biosynthesis genes (CESA4, CESA7, CESA8, and XTH22) and eight hemicellulose biosynthesis genes (IRX8, IRX9, IRX14-L, GXM1, GUX5, TBL3, TBL31 and TBL33) were also identified as hub genes in this module (Fig. 9). By comparison, no pectin and lignin biosynthesis enzyme genes were identified as hub genes in this module (Fig. 9).
The expression of most identified hub TF genes and enzyme genes of SCW biosynthesis was significantly higher in “JQ” than in “XMQ” at 5 DAP or 10 DAP or both 5 DAP and 10 DAP (Additional file 1: Table S4, Additional file 1: Table S5), which accords with the above GO enrichment results and suggested again that the early development stages of the hull were the key in determining the hull difference between the two tartary buckwheat cultivars. The expression of 17 identified SCW-related genes, including 7 regulatory genes and 10 enzyme genes, were further verified to be highly similar (r ≥ 0.78) to those observed in RNA-seq data by RT-qPCR analyses (Fig. 10). This indicated the reliability of the transcriptomic data and the identified genes that caused the hull difference between the two tartary buckwheat cultivars. To further verify that the different expression of these identified SCW biosynthesis genes in the early hull development stages was the reason for the hull difference formation between rice- tartary buckwheat and normal tartary buckwheat, the expression of the first-layer regulators (NST1, NST2, and SND1/NST3) of SCW biosynthesis was tested in the hull of the other three normal tartary buckwheat cultivars. As shown in Fig. 11, all three genes showed the highest expression in normal tartary buckwheat cultivars at 10 DAP, and the significantly different expression (fold change > 2) between “XMQ” and normal tartary buckwheat cultivars occurred both at both 5 DAP and 10 DAP. This suggested that the different expression of these SCW biosynthesis genes at the early hull development stages was the primary reason the hull difference between the rice-tartary buckwheat and normal tartary buckwheat.
‘‘In silico’’ promoter analysis of hub TFs and enzyme genes involved in SCW biosynthesis revealed the presence of binding sites of SCW-related TFs
The SCW-related TF binding cis-elements SNBE (NAC binding)  and SMRE (MYB biding)  were investigated in the promoter sequences of nine SCW-related TFs (NST1, NST2, SND1/SNT3, MYB54, MYB46/MYB83, MYB103 (2), C3H14, and C3H15), four cellulose biosynthesis-related enzyme genes (CESA4, CESA7, CESA8, and XTH22), eight hemicellulose biosynthesis-related enzyme genes (IRX8, IRX9, IRX14-L, GXM1, GUX5, TBL3, TBL31, and TBL33) and other 19 identified hub TFs. As results, the SNBE and SMRE cis-elements presented in the promoter sequences of almost all SCW-related TFs and enzyme genes (Table 1). In addition, a high number of SNBE cis-elements were also found in the promoter sequences of the other 19 hub TFs, while there no or very few SMRE cis-elements displayed in these hub TFs (Table 1). Due to ethylene response EIN3 TF was found in hub TFs and it had been demonstrated playing crucial regulatory role in many developmental processes , so we further investigated the EIN3 binding cis-element in promoter sequences of the above genes. As shown in Table 1, the EIN3 binding cis-element appeared in the promoter sequences of 33 out of the 40 genes, in which the highest number of EIN3 binding cis-elements were found in the first-layer NAC regulators of SCW biosynthesis (NST1, NST2, and NST3/SND1) (15, 15, and 9) and the xylem vessel formation regulator XND1 (11), suggesting that the EIN3 might be the direct up-stream regulator of the first-layer NAC regulators of SCW biosynthesis.