Carbohydrate metabolism and fertility related genes high expression levels promote heterosis in autotetraploid rice harboring double neutral genes

Background Autotetraploid rice hybrids have great potential to increase the production, but hybrid sterility is a major hindrance in the utilization of hybrid vigor in polyploid rice, which is mainly caused by pollen abortion. Our previous study showed that double pollen fertility neutral genes, Sa-n and Sb-n, can overcome hybrid sterility in autotetraploid rice. Here, we used an autotetraploid rice line harboring double neutral genes to develop hybrids by crossing with auto- and neo-tetraploid rice, and evaluated heterosis and its underlying molecular mechanism. Results All autotetraploid rice hybrids, which harbored double pollen fertility neutral genes, Sa-n and Sb-n, displayed high seed setting and significant positive heterosis for yield and yield-related traits. Cytological observations revealed normal chromosome behaviors and higher frequency of bivalents in the hybrid than parents during meiosis. Transcriptome analysis revealed significantly higher expressions of important saccharides metabolism and starch synthase related genes, such as OsBEIIb and OsSSIIIa, in the grains of hybrid than parents. Furthermore, many meiosis-related and specific genes, including DPW and CYP703A3, displayed up-regulation in the hybrid compared to a parent with low seed setting. Many non-additive genes were detected in the hybrid, and GO term of carbohydrate metabolic process was significantly enriched in all the transcriptome tissues except flag leaf (three days after flowering). Moreover, many differentially expressed genes (DEGs) were identified in the yield-related quantitative trait loci (QTLs) regions as possible candidate genes. Conclusion Our results revealed that increase in the number of bivalents improved the seed setting of hybrid harboring double pollen fertility neutral genes. Many important genes, including meiosis-related and meiosis-specific genes and saccharides metabolism and starch synthase related genes, exhibited heterosis specific expression patterns in polyploid rice during different development stages. The functional analysis of important genes will provide valuable information for molecular mechanisms of heterosis in polyploid rice. Electronic supplementary material The online version of this article (10.1186/s12284-019-0294-x) contains supplementary material, which is available to authorized users.


Background
Heterosis, or hybrid vigor, is a complex biological phenomenon, which is improved or superior phenotypic performance of hybrid in comparison to one or both parents, such as enhanced grain yield, stress tolerance and biomass production. Heterosis has been extensively applied to increase the rice yield in the world (Cheng et al. 2007). However, productivity of rice has been stagnant in the past few years. Polyploid species play an important role in breeding programs, such as cotton (Flagel et al. 2008), wheat (Goncharov et al. 2007), and rapeseed (Albertin et al. 2006). Rice polyploidization is an effective method to increase the size of rice genome and improve the wide adaptability (Cai et al. 2001;Guo et al. 2017).
The polyploid rice hybrids showed stronger biological advantage and yield potential compared with diploid rice hybrids, and has attracted the attention of many rice researchers (Shahid et al. 2011;Wu et al. 2013;Guo et al. 2017). However, autotetraploid rice has many unfavorable traits, especially low seed setting, which limits its commercial utilization (Shahid et al. 2010(Shahid et al. , 2013aChen et al. 2018a). Recently, our research team found that polyploidy enhanced pollen sterility loci interactions and increased chromosomal abnormalities in autotetraploid hybrid rice (Wu et al. 2015), and also revealed that intersubspecific diploid and autotetraploid hybrid rice sterility could be overcome by double neutral genes (Shahid et al. 2013b;Wu et al. 2017). After years of efforts, our research group have developed few neo-tetraploid rice lines with high seed setting (> 80%) (Guo et al. 2016(Guo et al. , 2017, and two new neo-tetraploid rice have been registered for the PVP (Protection for new varieties of plant) in China . Moreover, neo-tetraploid rice could overcome the sterility and produce high heterosis in autotetraploid hybrid rice (Guo et al. 2017). Two photoperiod-and thermo-sensitive genic male sterile lines (PS006 and PS012) of polyploid rice showed stronger hybrid vigor and great potential for improving rice quality and productivity ).
Next generation high-throughput sequencing, such as RNA sequencing (RNA-seq) and microarray technology, is widely used to investigate gene expression and function. The RNA-seq enabled us to understand differentially expressed genes associated with abiotic stresses and pollen development in diploid rice (Jin et al. 2013;Hu et al. 2016;Fu et al. 2017). The RNA-seq has also been applied to detect differentially expressed genes between diploid and autotetraploid rice during pollen development Chen et al. 2018a;Li et al. 2018). Furthermore, RNA-seq has been widely used to investigate heterosis in various plants, such as wheat , maize (Ma et al. 2018), rapeseed (Shen et al. 2017), tobacco (Tian et al. 2018) and rice (Wei et al. 2009). The differentially expressed genes were found to be closely associated with heterosis in super rice LYP9 and its parents through RNA-seq (Wei et al. 2009). Later, many differentially expressed genes closely related to root heterosis at tillering and heading stage were detected in super hybrid XY9308 and its parents through RNA-seq (Zhai et al. 2013). Chen et al. (2018b) compared the transcriptomes between super hybrid Wufengyou T025 (WFYT025) and its parents during young panicle development, and suggested that carotenoid biosynthesis and plant hormone signal transduction were enriched in differentially expressed genes, and these genes were related to the grain number heterosis. A number of genes associated with leaf, anthers and ovary heterosis in neo-tetraploid rice hybrids were identified by RNA-seq, which were related to photosynthesis and metabolic process and transport (Guo et al. 2017).
Our previous study indicated that saccharide abnormal distribution and down-regulation of saccharide transport genes may cause pollen sterility and lead to low seed setting in autotetraploid rice (Chen et al. 2018a). In this study, to increase yield of autotetraploid rice, we used an autotetraploid rice line harboring double neutral genes for pollen fertility at Sa and Sb loci, which could overcome F 1 sterility when it crossed with low fertility autotetraploid rice Chen et al. 2018a). We thus primarily aimed to evaluate heterosis mechanism of neo-tetraploid and autotetraploid rice harboring double neutral genes, and to observe the role of chromosome configuration and behavior in heterosis and fertility. In addition, we detected differentially expressed genes between parents and hybrid in nine tissues at three development stages using RNA-seq, which would provide insights into the molecular mechanism underlying heterosis in autotetraploid and neo-tetraploid rice, and provide new germplasm for polyploid rice breeding.

Results
Heterosis evaluation of hybrids generated by crossing of autotetraploid with neo-tetraploid rice Analysis of the agronomic traits of five hybrids, which were developed by crossing autotetraploid rice line (T449) with five neo-tetraploid rice lines, showed significant improvement in important yield-related traits, including number of filled grains per plant, yield per plant and seed setting. Evaluation of heterosis indicated that the values for mid-parent heterosis (MPH) were positive for all the traits except grain length and total grains per plant, and the highest MPH was found for grain yield per plant (170.89%). The high-parent heterosis (HPH) values were positive for the filled grains and grain yield per plant, and the highest HPH was detected for filled grains per plant (71.10%) (Additional file 1: Table S1).
Then, we further selected a hybrid (T449 × H1) for transcriptome analysis to analyze the heterosis mechanism in detail. The hybrid had high seed setting, pollen and embryo sac fertility ( Fig. 1; Table 1 and Additional file 2: Table S2), although the maternal line, T449, had low pollen fertility and seed setting. The hybrid displayed significant positive MPH for all the traits, and the values for MPH were very high for filled grains per plant, grain yield per plant and seed setting (Table 1). Meanwhile, the hybrid also showed positive HPH values for most of the traits except effective number of panicles per plant, 1000-grain weight, and grain length and width (Table 1).
Chromosome configuration at diakinesis and metaphase I in F 1 hybrid compared to its parents The tetravalent was the most type of chromosome configuration at diakinesis and metaphase I in the parents.
The numbers of bivalent chromosomes were higher in H1 (high seed setting) than T449 (low seed setting), and univalent type of chromosomes were higher in T449 than both H1 and F 1 (Table 2, Fig. 2). However, bivalent chromosomes were the most frequent in F 1 (Table 2), and 9.07 and 12.57 bivalents were found in each pollen mother cell (PMC) at diakinesis and metaphase I, respectively. There were significant differences in the numbers of chromosome configurations between F 1 hybrid and parents, and these results indicated that the increase of bivalent type may improve pollen fertility and seed setting. The chromosome configurations of tetravalent types were divided into five types, including ring shape (Fig. 3a1-a4), chain shape ( Fig. 3b1-b4), frying Fig. 1 Comparisons of morphological characteristics between F 1 hybrid and its parents. a Plant appearance of F 1 , T449 and H1. Pollens of T449 (b), F 1 (c), and H1 (d), blue arrows indicate normal pollens, red arrows indicate abnormal pollens. Bar = 100 μm. Embryo sacs of T449 (e), F 1 (f) and H1 (g), white arrows indicate antipodal cells and red arrows indicate two polar nuclei. Bar = 100 μm pan shape ( Fig. 3c1-c4), "X" shape ( Fig. 3d1-d4) and "OK" shape ( Fig. 3e1-e4). The frequency of ring shape was the most frequent chromosome configuration of tetravalent types at diakinesis and metaphase I (Table 2).

Differentially expressed genes in F 1 hybrid and parents
In order to investigate transcriptome changes in F 1 and its parents, the transcriptome profiles of F 1 and parents were analyzed in nine tissues, including anthers (P1) and flag leaves (L1) at meiosis stage, and flag leaves (L2), leaf sheath (Z2), anther (P2) and embryo sac (E2) at pre-flowering stage, and flag leaves (L3), leaf sheath (Z3) and grain (P3) at three days after flowering (Additional file 4: Figure S2). In total, more than 3.7 billion clean reads were detected in different samples. We aligned clean reads against the Nipponbare reference genome (MSU 7.0), and 92.52% to 96.32% annotated transcripts of the reference genome was obtained in our materials (Additional file 6: Table S4). The correlation for the gene expression level from three biological replicates of each line was more than 0.8 (Additional file 7: Table S5), and principal component analysis (PCA) indicated that three biological replicates were clustered together, and flag leaves and leaf sheath were also clustered together (Additional file 8: Figure S3). The correlations between F 1 and its parents were investigated in different samples by hierarchical cluster analysis using Cluster 3.0 software. The results demonstrated that F 1 and its parents always assembled into primary groups at  the same tissue, and the transcriptome profiles of F 1 were similar to H1, and these results were consist with the morphological and cytological observations (Additional file 9: Figure S4). A total of 12 DEGs were randomly selected for qRT-PCR validation. We compared the qRT-PCR results, and the expression trends were consistent with RNA-seq data, and a correlation coefficient was R 2 = 0.8806 (Additional file 10: Figure  S5), which demonstrated that RNA-seq data is reliable.
Differentially expressed anther specific genes were associated with meiosis stage-specific genes in F 1 compared to its parents T449 is an autotetraploid rice line with low fertility and H1 is a neo-tetraploid rice line with high fertility, so we specifically investigated meiosis-specific genes associated with high fertility. We compared expression levels of genes between F 1 vs T449 and H1 vs T449, and 381 genes were found to be commonly up-regulated between both comparisons in anthers during meiosis (Additional file 16: Table S8). GO analysis of these 381 genes showed that six biological process categories and nine molecular function categories were significantly enriched (Additional file 17: Fig. 4 Gene ontology (GO) enrichment heat map for DEG FPU in 9 tissues (GO terms were selected based on their appearance at least in three tissues or more). L1 and P1 represent flag leaves and anthers at meiosis stage, respectively. L2, P2, E2 and Z2 represent flag leaves, anther, embryo sac and leaf sheath at pre-flowering stage, respectively. L3, P3 and Z3 represent flag leaves, grain and leaf sheath at three days after flowering, respectively Table S9). We compared 381 up-regulated genes with microarray data of wild type rice anther meiosis stage-specific expression, and meiosis-related genes (Fujita et al. 2010;Tang et al. 2010;Deveshwar et al. 2011;Yant et al. 2013;Luo et al. 2014;Wright et al. 2015), and identified 4 meiosis-related genes including LOC_Os09g32020 (OsDFR), LOC_Os01g68870 (MSP1) and LOC_Os12g24420 and LOC_Os10g06770 (CDKG1), and 26 meiosis-specific expressed genes (Additional file 18: Table S10). Four pollen-related genes were identified from further analysis of 26 meiosis-specific expressed genes, including LOC_ Os03g07140 (DPW), LOC_Os04g24530 (OsACOS12), LOC_Os06g40550 (OsABCG15; PDA1) and LOC_ Os08g03682 (CYP703A3).
Among 381 genes, 47 genes exhibited down-regulation in T449 compared to its diploid counterpart (E249) (Chen et al. 2018a) (Fig. 6a). Interestingly, two meiosis-related and 19 meiosis-specific genes were identified in these genes (Additional file 19: Table S11), including LOC_Os09g32020 (OsDFR) and LOC_ Os12g24420 (CDKG1), and analysis of 19 meiosis-specific genes revealed four pollen-related genes, including LOC_Os03g07140 (DPW), LOC_Os04g24530 (OsA-COS12), LOC_Os06g40550 (OsABCG15; PDA1) and LOC_Os08g03682 (CYP703A3). We performed the predicted protein-protein interactions of 47 genes using STRING v10, and the results showed that 17 genes constituted genetic sub-networks, including one meiosis-related and six meiosis-stage specific genes. The meiosis-related gene (OsDFR, LOC_Os09g32020) interacted with two meiosis-specific gene, including LOC_Os08g03682 (CYP703A3), encodes cytochrome P450 hydroxylase and LOC_Os04g24530 (OsACOS12), encodes acyl-CoA synthetase 12, which interacted with six differently expressed genes, including three cytochrome P450, two AMP-binding enzyme and a protein binding protein (Fig. 6b). Human Diseases Endocrine and metabolic diseases 3 2 3 2 1 1 1 0 1 * and ** denote significant enrichment of DEG FPU among functional categories with P < 0.05 and P < 0.01, respectively. L1 and P1 represent flag leaves and anthers at meiosis stage, respectively. L2, P2, E2 and Z2 represent flag leaves, anther, embryo sac and leaf sheath at pre-flowering stage, respectively. L3, P3 and Z3 represent flag leaves, grain and leaf sheath at three days after flowering, respectively Co-expression network analysis of differentially expressed genes in different tissues by weighted gene coexpression network analysis (WGCNA) WGCNA, which is a systems biology tool, was used to understand the relationships and networks in a set of genes. In this study, WGCNA was constructed using RNA-seq data, and 28 WGCNA modules were identified (Fig. 7a, Additional file 20: Table S12). The gene numbers in these modules were ranged from 30 (MEwhite module) to 9294 (MEgrey module). Interestingly, the turquoise and grey modules consist of 67.72% of the genes in the network analysis. We found that some modules showed correlation with the different tissues in F 1 and parents (Fig. 7b), for example, MEbrown module in leaf, MEblue module in mature anther, MEred module in mature embryo sac, and MEturquoise and MEgrey module in three tissues (anther, embryo sac and grain), which indicated that these modules may play putatively important roles in tetraploid rice leaf and reproductive organs. Furthermore, a total of 1335 genes were involved in the MEbrown module, and GO enrichment analysis showed significant enriched terms that were related to photosynthesis light harvesting, light reaction, carboxylic acid metabolic process and chlorophyll metabolic process. These results indicated that brown module genes may play an important role in the photosynthesis in tetraploid rice. In total, 9291 genes were involved in the MEturquoise module, and GO analysis revealed significant terms associated with DNA repair, carbohydrate metabolic process and transport, which indicated that MEturquoise module plays an important role in the fertility and yield of tetraploid rice (Additional file 21: Table S13).
Association between DNA sequence variations and differentially expressed genes in F 1 hybrid compared to parents A total of 100,395,344 and 150,936,032 clean reads were obtained in T449 and H1 by using genome re-sequencing, respectively. Approximately 98.16% (T449) and 96.11% (H1) of clean reads were mapped onto the Nipponbare reference genome, and the reads coverage depths were 35× and 49× in T449 and H1, respectively (Additional file 22: Table S14). A total of 912,892 SNPs and 195,976 InDels were detected between T449 and H1 by using the two filter conditions (coverage ≥10 and ≤ 100, and removal of heterozygous SNPs and InDels). We found that about 5% of SNPs and 6% of InDels were detected in intergenic regions, and 60% SNPs and 65% InDels were identified in up or down regulatory regions, which might be related to the differentially expressed genes (Additional file 23:  Table S17). These results were nearly consistent with the GO enrichment of DEG FPU .  Expression patterns of saccharides metabolism and starch synthase related genes in the hybrid compared to parents The GO and KEGG analyses of DEG FPU showed that there were significant differences for carbohydrate metabolic process in nine tissues between F 1 and its parents, and DEG FPU were involved in sucrose synthase, cell wall invertase, 6-phosphofructokinase, and hexokinase. Many saccharide metabolic genes were up-regulated in the F 1 compared to its parents in L1, P1 and P3 (Additional file 26: Table  S18). Interestingly, the saccharide transporters were up-regulated in the F 1 compared to its parents in L1, P1 and P3, and these results were consistent with saccharide metabolic genes (Additional file 26: Table S18). The two saccharide transporter genes (LOC_Os02g10800 and LOC_Os03g07480) were also up-regulated in F 1 compared to both parents in P3 (Additional file 26: Table S18). In addition, the invertase (OsINV3 and OsINV4), sucrose synthase (OsSUS3 and OsSUS4), hexokinase gene (OsHXK6), starch branching enzyme (OsBEIIb) and two starch synthase genes (OsSSIIIa and wx) displayed higher levels of expressions in F 1 than parents in the grains (three days after flowering) (Fig. 8a). Moreover, the promoter regions of OsSUS3, OsBEIIb and OsSSIIIa also exhibited differences between maternal and paternal rice lines by re-sequencing (Fig. 8b).

High heterosis and the frequency of bivalents in F 1 hybrid harboring double neutral genes
In the current study, the paternal line was high fertility neo-tetraploid rice and maternal line was autotetraploid rice harboring double neutral genes. The hybrid displayed stronger heterosis for yield and yield-related traits, such as filled grains per plant, total grains per plant, grain yield per plant and seed setting. These results were consistent with the previous studies, where autotetraploid hybrids exhibited high heterosis for filled grains per panicle, grain yield per plant and seed setting (Shahid et al. 2012;Guo et al. 2017). In addition, the pollen and embryo sac fertilities were investigated to evaluate the fertility of hybrid and its parents. Our results showed that the embryo sac fertility of hybrid and its parents was higher than 89%, and pollen fertility of hybrid and paternal line was high, while pollen fertility of maternal line was low. These results revealed that pollen fertility has a greater impact on seed setting than embryo sac fertility in autotetraploid rice hybrid and its parents.
It is well known that meiosis process has a great effect on plant reproductive development, and chromosome behavior and configuration play an important role in the plant meiosis and directly correlated with pollen fertility. The way of quadrivalent separation depends on chromosome configuration. Many quadrivalent Fig. 6 Gene expression levels of 47 genes and predicted protein-protein interaction network. a, The distribution of 47 genes exhibited up-regulation in F 1 hybrid and H1 compared to T449 and down-regulation in T449 compared to E249 during meiosis stage, b, Predicted protein-protein interaction network of differently expressed genes (black), meiosis-specific (blue) and meiosis-related (red) genes. T represents T449, F represents hybrid, and H represents H1 chromosomes were found in Triticum monococcum during diakinesis and metaphase I (Kim and Kuspira 1993), while chain, ring and frying pan shapes were three main types of chromosome configurations in autotetraploid rice (Luan et al. 2007;He et al. 2011a). Ring, chain, frying pan, "X" and "OK" shapes were observed in the present study, and we have drawn these shape models based on the observation. The ring shape quadrivalent was found more frequently compared to other quadrivalent types in the present study. Our results were in agreement with the previous studies, who also observed ring shape quadrivalent in autotetraploid rice (Luan et al. 2007;He et al. 2011b). He et al. (2011b) revealed that high frequency of bivalent was related to high pollen fertility and seed setting in autotetraploid rice hybrid. We also observed higher frequency of bivalents in the autotetraploid hybrid than parents at diakinesis and metaphase I. Chromosome behavior had a direct relationship with pollen fertility and seed setting in autotetraploid rice hybrids and interspecific hybrid between Brachiaria ruziziensis and B. brizantha (Adamowski et al. 2008;He et al. 2011b;Guo et al. 2016). Here, the frequency of normal chromosome behavior was significantly higher in the hybrid with high fertility than maternal line with double neutral genes. It is worth to mention that the frequency of abnormal chromosome behavior is much higher during anaphase II Fig. 7 WGCNA based gene expression matrix between hybrid and parents. a Hierarchical cluster tree showing co-expression modules identified by WGCNA. Each leaf in the tree represents one gene. The major tree branches constitute 28 modules labeled with different colors. b Modulesample relationship. Each row corresponds to a module. Each column corresponds to a sample. The color of each cell at the row-column intersection indicates the correlation coefficient between the module and the sample than other stages in this study. Here, many types of abnormal chromosome behaviors, including asynchronous meiocytes, abnormal spindles and straggling chromosomes were observed at anaphase II (Additional file 5: Table S3). Of these abnormalities, asynchronous meiocytes was the highest, and the frequencies of asynchronous meiocytes were 66.9%, 65.3% and 56.0% in T449, F 1 and H1, respectively. We inferred that most of asynchronous meiocytes could develop normal tetrad according to their morphology, which might be the major reason for lower frequency of the normal cells at anaphase II than other phases.
The expression patterns of meiosis and meiosis-related genes promote high fertility in F 1 hybrid A number of meiosis-related and meiosis-specific genes were detected in rice (Fujita et al. 2010;Tang et al. 2010;Deveshwar et al. 2011;Yant et al. 2013;Luo et al. 2014;Wright et al. 2015). A total of 55 meiosis-related or meiosis-stage-specific genes were found to be down-regulated, which increased pollen sterility loci interactions in autotetraploid rice hybrids (Wu et al. 2015). In the present study, the meiotic stages were determined by the floret length and according to the observation by 4′, 6-diamidino-2-phenylindole (DAPI) staining (He et al. 2011a). The DAPI staining could clearly distinguish between meiotic stages and other pollen development stages. For the proper understanding of meiosis-related and meiosis-specific genes between hybrid and parents, we dissected anthers from floret for RNA-seq analysis. A total of four meiosis-related and 26 meiosis-stage specific genes were identified, which were found to be up-regulated in hybrid and paternal line compared to maternal line. Interestingly, of the four meiosis-related and 26 meiosis-specific genes, two meiosis-related and 19 meiosis-specific genes were also found to be down-regulated in autotetraploid rice compared to diploid progenitors (Chen et al. 2018a). For example, DPW gene encodes a fatty acyl ACP reductase, and was found to be essential for anther cuticle, pollen wall and pollen sporopollenin biosynthesis (Shi et al. 2011), and maternal and paternal lines showed upstream and intron variations in DPW gene. OsACOS12, which is an acyl-CoA synthetase, is essential for sporopollenin synthesis in rice (Li et al. 2016a;Yang et al. 2017), and maternal and paternal lines displayed downstream, upstream, intron and non-synonymous variations in OsACOS12. PDA1 encodes an ABC transporter (OsABCG15) and required for the transport of lipidic precursors for anther cuticle and pollen exine development (Zhao et al. 2015). CYP703A3, cytochrome P450 hydroxylase, is involved in the tapetum degeneration retardation, a known pollen exine formation (Yang et al. 2014). The meiosis-related gene (LOC_Os12g24420) encoded cyclin-dependent kinase, which is homolog to CDGK1 in Arabidopsis, and CDKG1 protein kinase is crucial for synapsis and recombination in Arabidopsis during meiosis (Zheng et al. 2014), and we observed changes in the downstream and intron region of CDKG1 in maternal and paternal lines. These results suggested that the expression profiles of important meiosis-related or meiosis-specific genes have a significant effect on the fertility of polyploidy hybrid rice.

Dominance, non-additive and yield-related genes/QTLs contribute to heterosis
The DEG FP were divided into five basic groups according to their expression profiles, including over-dominance (HBP), under-dominance (LBP), dominance (CHP and CLP), and mid-parent (BBP) . In our data, the dominance expression was the most prevalent class among DEG FP (46.34-77.71%). Similarly, the dominance expression patterns were found to be the most abundant among DEG FP in wheat and rice hybrids by RNA-seq analysis (Zhai et al. 2013;Liu et al. 2018). These results indicated that dominance expression have a great effect on the performance of hybrids. According to gene expression levels of hybrid and its parents, the gene expression profiles of F 1 could be divided into two types, the first type is called as additive expression, which is contributed by each allele from its parents in a hybrid, and another is non-additive expression that differed from the mid-parent value (MPV) (Wei et al. 2009). In the previous study, whether or not a transcript shows non-additive expression is most likely to be affected by the contributions of cisand trans-acting element of a gene (Zhang et al. 2008). In this study, NAGs only accounted for 2.1-6.5% of the total expressed genes, but 57.3-72.4% of DEG FPU were NAGs in each tissue. These results were consistent with diploid rice (Wei et al. 2009), wheat ) and maize heterosis (Swanson-Wagner et al. 2006), who also detected NAGs in F 1 . The previous studies have shown that NAGs play vital roles in heterosis (Zhang et al. 2008;Liu et al. 2018), and NAGs were associated with circadian rhythm, flowering time, and panicle branching in rice (Li et al. 2016b). Overall less number of NAGs detected in this study, but major portion of DEG FPU was constituted of NAGs. Therefore, we speculated that NAGs play important roles in high F 1 heterosis of polyploid rice.
The potential relationships between differently expressed genes and QTLs have been proposed in many yield-related QTL regions using RNA-seq (Zhai et al. 2013;Chen et al. 2018b). A recent study showed that nine genes (Hd3a, TAC1, Ghd8, Sd-1, NAL1, Hd1, GW6a, IPA1 and DEP1) have a major impact on heterosis in diploid hybrid rice (Huang et al. 2016). In our study, DEP1, which is related to rice panicle and located in SKPNB (spikelet number, AQBK037), was found to be up-regulated in grains (three days after flowering) of F 1 hybrid compared to maternal line and also up-regulated in leaf of hybrid than parents during meiosis stage. GW6a, which involved in grain-weight and located in TSDWT (1000-seed weight, AQEB012), revealed much higher expression in anthers of hybrid than parents during meiosis. The semi-dwarf gene (Sd-1), which is involved in biosynthesis of gibberellin and located in GRYLD (grain yield, AQQ005), exhibited higher levels of expressions in leaf sheath (three day after flowering) of hybrid than parents. Hd3a, which is related to rice flowering and regional adaptation, and located in FGRNB (filled grain number, AQCF008), was up-regulated in leaf (before flowering) of hybrid compared to parents. These results showed that these candidate genes in QTL regions may contribute to heterosis in autotetraploid hybrid rice.

Saccharides metabolism and starch synthase related genes play an important role in heterosis
Carbohydrate metabolism plays an essential role in the plant growth and development. RNA-seq analysis showed that processes of carbohydrate metabolism are related to heterosis in rice (Wei et al. 2009;Zhai et al. 2013). In addition, our research group reported that abnormal distribution of saccharides and saccharides-related genes may influence pollen fertility and cause decrease in the yield of autotetraploid rice (Chen et al. 2018a). Here, the DEG FPU was significantly enriched in the carbohydrate metabolism process between hybrid and its parents in nine tissues sampled across different development stages. The grain filling stage is an important part of growth and development in rice, and a large quantity of carbohydrates are synthesized and transported into the grain at this stage, particularly the embryo starts exponential growth at three days after flowering that shows dual rhythmicity (Itoh et al. 2005). The stage of three days after flowering is important stage for grain development, so we focused on saccharides metabolism in the grains three days after flowering. The carbohydrates supply from the leaves to pollen or grain involves sucrose transport and degradation, monosaccharides formation and transport, and starch generation (Ruan 2012). There are two types of enzymes that catalyze the sucrose degradation in plants, one is sucrose synthase (SUS) and the other is invertase (Ruan et al. 2010). Here, the invertase (OsINV3 and OsINV4) and sucrose synthase (OsSUS3 and OsSUS4) were found to be up-regulated in the grains of hybrid compared to parents. After sucrose degradation, the resulting hexoses undergo phosphorylation by hexokinase for starch synthesis. Subsequently, hexokinase plays important role in hexose signaling and sensing (Cho et al. 2009;Kim et al. 2016). The hexokinase gene (OsHXK6) was up-regulated in grains (three after days flowering) of the hybrid compared to its parents. The starch synthase, starch debranching and starch branching enzyme have a great influence on the starch generation and metabolism (Zeeman et al. 2010). The starch branching enzyme (OsBEIIb) and two starch synthase (OsSSIIIa and wx) genes were up-regulated in the grains of hybrid compared to parents. In addition, we detected differences between maternal and paternal rice lines in the promoter regions of OsSUS3, OsBEIIb and OsSSIIIa by re-sequencing. Consequently, the genetic effects of OsSUS3, OsBEIIb and OsSSIIIa may cause allelic heterozygosity in promoter regions of hybrid.
Sucrose and monosaccharide transporters are important proteins for the translocation of saccharides from source to sink organs (Ruan et al. 2010). OsSUT1 primarily play a role in the transport pathway (Scofield et al. 2007). In our study, sucrose transporter (OsSUT1) displayed much higher expression patterns in the grain (three after days flowering) of hybrid than parents. OsBT1, which encodes putatively ADP-glucose transporter and localizes in the amyloplast envelope membrane, plays a crucial role in starch synthesis (Cakir et al. 2016;Li et al. 2017a). We found that OsBT1 was up-regulated in the grain (three after days flowering) of hybrid compared to parents. Transcriptome profiling showed high expression levels of saccharides metabolism and starch synthase related genes in the hybrid, which might be an indication of enhanced source for sink tissues and resulted in high yield of hybrid.
In our previous studies, we found that the double neutral genes can overcome the hybrid sterility in autotetraploid rice , and detected specific differentially expressed genes associated with fertility and heterosis in neo-tetraploid rice by RNA-seq (Guo et al. 2017). Here, an autotetraploid rice line (T449), harboring Sa-n and Sb-n double neutral genes for pollen sterility loci, was used to generate the hybrids by crossing with neo-tetraploid rice, and investigated the heterosis and fertility by cytological and RNA-seq methods. We further want to understand the role of double neutral genes in heterosis and fertility of neo-tetraploid rice. Therefore, we observed chromosome behavior and gene expression patterns during important growth stages. The results showed that seed setting of F 1 hybrid improved with the increase in number of bivalents, and many important genes, including meiosis-related and meiosis-specific genes and saccharides metabolism and starch synthase related genes, exhibited heterosis specific expression patterns in polyploid rice during different development stages.

Conclusions
In this study, we observed the chromosome behavior and configuration in hybrid and its parents, and found higher frequency of bivalent and normal chromosome behavior in hybrid than parents, which promoted high fertility (heterosis) in the hybrid harboring double neutral genes. Furthermore, we systematically investigated the global transcriptome of hybrid and its parents by RNA-seq. We obtained a large number of DEG FPU , and detected substantial candidate genes, including meiosis-related and meiosis-specific genes, saccharides metabolism and starch synthase related genes, which were up-regulated in hybrid having improved fertility and yield. Our results provided new resource for polyploid rice breeding and exploring of these candidate genes will provide valuable information for revealing molecular mechanisms of heterosis in polyploidy rice.

Rice material
An autotetraploid rice line, DN18-4x (T449), harboring Sa-n and Sb-n double neutral genes for pollen sterility loci, was used to generate the hybrids by crossing with five neo-tetraploid rice lines, including Huaduo 1 (H1), Huaduo 2 (H2), Huaduo 3 (H3), Huaduo 4 (H4) and Huaduo 8 (H8). All the materials were planted at the experimental farm of South China Agricultural University (SCAU) under natural conditions, and management practices followed the recommendations for the area.

Investigation of agronomic traits and data analysis
Agronomic traits, including plant height, effective number of panicles per plant, grain length and width, 1000-grain weight, filled grains per plant, total grains per plant, grain yield per plant and seed setting, were investigated. The standard for investigating these agronomic traits was according to the protocols of People's Republic of China for the registration of a new plant variety DUS (Distinctness, Uniformity and Stability) test guidelines of rice (Guidelines for the DUS test in China, 2012) Guo et al. 2017). The mid-parent heterosis (MPH) and high-parent heterosis (HPH) were estimated by the following formula: MPH = (F 1 − MP)/MP × 100%, and HPH = (F 1 − HP)/HP × 100%, where F 1 related to the performance of hybrid, HP was defined as the performance of the best parent, and MP was defined as an average performance of two parents .

Cytological observation
The chromosome configuration and behaviors were observed according to Wu et al. (2014). The inflorescences of F 1 and its parents lines were collected from the shoots of rice plants with 0 to 4 cm between their flag leaf cushion and the second to last leaf cushion, and fixed in Carnoy's solution (ethanol: acetic acid = 3:1) for 24 h, and the samples were stored in 70% ethanol at 4°C after washing two times. The anther was removed from the floret and placed in a small drop of 1% acetocarmine on a glass slide. The glass slide was covered with a slide cover after 2-3 min, and observed under a microscope (Motic BA200).
The pollen fertility was observed according to Shahid et al. (2013b). The normal and abnormal pollens were observed by staining with 1% I 2 -KI under a microscope (Motic BA200). The whole mount eosin B confocal laser scanning microscopy (WE-CLSM) was used to investigate the embryo sac fertility in F 1 and its parents according to Li et al. (2017b) with minor modifications. The ovary was dissected from the floret, and was hydrated in 70%, 50%, 30%, 10% ethanol and distilled water for 30 min each. Then, the samples were dehydrated in 10%, 30%, 50%, 70%, 90% and 100% ethanol for 30 min after eosin B staining for 12 h. Finally, the samples were shifted into a mixture solution (ethanol and methyl salicylate = 1:1) for 2 h, and then keep in pure methyl salicylate and observed under a laser scanning confocal microscope (Leica SPE).

RNA-seq experiments and data analysis
All samples were collected during meiosis, pre-flowering, and three days after flowering. The meiosis stage is a crucial event for the sexual reproduction of eukaryotes to form haploid spores and gametes ). The pre-flowering stage is an important stage for pollen and embryo sac fertility. The carbohydrates are synthesized and transported into the grains in large quantity during grain filling stage (Itoh et al. 2005). The nutrients produced by leaves are transported to other organs through the leaf sheath, so the leaf sheath has an important role in the transportation of energy. Flag leaf is one of the most important photosynthetic organs in rice and has an important impact on crop yield and quality. In addition, pollen and embryo sac have a significant impact on rice fertility and yield (Shahid et al. 2013b;Wu et al. 2015;Li et al. 2017b). Hence, we collected the nine tissues during these development stages from hybrid and its parents in three biological replicates, including anthers (P1) and flag leaves (L1) at meiosis stage, and flag leaves (L2), leaf sheath (Z2), anther (P2) and embryo sac (E2) at pre-flowering stage, and flag leaves (L3), leaf sheath (Z3) and grains (P3) at three days after flowering (Additional file 4: Figure S2).
All tissues of hybrid and its parents were harvested in three biological replicates and immediately kept at − 80°C for RNA extraction. The total RNA was extracted according to the manual instructions of the TRIzol Reagent (Life technologies, California, USA). The library was prepared according to the vendor's recommended protocol. The RNA-seq was performed on the Illumina HiSeq 4000 sequencing platform (LC Sceiences, USA). Using the Illumina paired-end RNA-seq approach, we sequenced the transcriptome that generated millions of paired-end reads. Low quality reads, including reads containing sequencing adaptors, reads containing sequencing primer and nucleotides with quality score lower than 20, were removed. The mapped reads from each sample were assembled using StringTie, and all transcriptome samples were mixed to reconstruct a comprehensive transcriptome by perl scripts. After the generation of transcriptome, the StringTie and Ballgown were used to evaluate the gene expression levels. String-Tie was used to perform expression level for mRNAs by calculating FPKM (fragments per kilobase of transcript per million fragments mapped reads), and false discovery rate (FDR) was used to determine the threshold of the P-value in multiple tests.
The Venny software was used to identify the overlapped differentially expressed genes in different samples (http:// bioinfogp.cnb.csic.es/tools/venny/index.html). Hierarchical analysis was carried out for all genes using Cluster 3.0 software after normalization. Transcription factor analysis was done according to transcription factor data (Jin et al. 2017). Gene Ontology (GO) enrichment analysis was employed for functional categorization by using AgriGO tool (http://systemsbiology.cau.edu.cn/agriGOv2/).

Expression patterns of differentially expressed genes (DEGs)
The expression patterns of DEGs were defined according to Liu et al. (2018). We defined the gene expression in F 1 as EF 1 , and genes expression of both parental lines as ET449 and EH1.We defined the average value of both parental lines as MPV (mid-parental value). If the F 1 was significantly (FDR ≤ 0.05 and fold change ≥2) different from MPV, we defined these genes as non-additive genes (NAGs), if there was non-significant difference between F 1 and MPV, these genes were defined as additive genes. Classification of DEG FP was performed according to the expression of EF1 relative to ET449 and EH1. ">" and "<" represents statistically higher or lower, and "=" represents statistically similar. If EF1 > ET449 > EH1, or EF1 > EH1 > ET449, then expression patterns of these genes were considered as higher than both parents (HBP); if EF1 = ET449 > EH1, or EF1 = EH1 > ET449, then expression patterns of these genes were considered as close to higher parent (CHP); if ET449 > EF1 > EH1, or EH1 > EF1 > ET449, then expression patterns of these genes were considered as between two parents (BBP); if EF1 = ET449 < EH1, or EF1 = EH1 < ET449, then expression patterns of these genes were considered as close to lower parent (CLP); if EF1 < ET449 = EH1, or EF1 < ET449 < EH1, or EF1 < EH1 < ET449, these genes were considered as lower than both parents (LBP).

Real-time qRT-PCR analysis
A total of 12 DEGs were randomly selected for validation of RNA-Seq data by qRT-PCR. The gene-specific primers were designed using Primer Premier 5.0 software, and checked in the NCBI (National Center for Biotechnology Information) for specific primers (Additional file 27: Table S19). Total RNA was taken from sequenced samples, and the first-strand cDNA was synthesized using the Transcriptor cDNA Synth. Kit 1 (Roche) according to the manufacturer's instructions. The qRT-PCR reaction procedure was 30 s at 95°C, with 40 cycles of 95°C denaturation for 10 s and 60°C annealing and extension for 30 s, and performed on the Light-cycler480 system (Roche). The genes relative expression levels were calculated using the 2-ΔΔCt method (Livak & Schmittgen 2001). All qRT-PCR reactions were performed in triplicate.
Mapping DEG FPU to rice QTLs and weighted gene coexpression network analysis (WGCNA) Rice QTL data with physical positions on the MSU Rice Genome Annotation Project Release 6.1 were acquired from Gramene (Youens-Clark et al. 2011). The DEG FPU were mapped onto 1019 yield related QTLs and 26 yield-related traits using gene coordinates from the MSU Rice Genome Annotation Project. The gene co-expression networks were used WGCNA package in R (Langfelder & Horvath 2008). To reduce noise, genes with total FPKM < 5 in 81 samples were removed. The modules were obtained using the automatic network construction with default settings.