Search
Submit Manuscript
1.6CiteScore
Advanced Search
MAP  >  Grass Research
2023 Volume 3
Article Contents
Abstract
Introduction
Types of functional rice
High resistant starch rice
Colored rice
Giant embryo rice
Low glutelin rice
Golden rice
Fe and Zn fortified rice
Breeding technology of functional rice
Traditional polymerization breeding
Transgenic polymerization breeding
Molecular marker-assisted selection polymerization breeding
Cultivation and regulation technology of functional rice
Effect of light on the accumulation of active components in rice
Effect of temperature on the accumulation of active components in rice
Effect of water management on the accumulation of active components in rice
Effect of fertilizer application on the accumulation of active components in rice
Effect of harvesting time on the accumulation of active components in rice
Prevention and control of pests and diseases for functional rice
Conclusions
Acknowledgments
Conflict of interest
Supplementary information
Rights and permissions
References
About this article
Next Previous
ARTICLE   Open Access    

Integration analysis of full-length transcriptomics and metabolomics provides new insights into the mechanism of sex differentiation in buffalograss (Buchloe dactyloides)

  • 1.

    Institute of Grassland, Flowers, and Ecology, Beijing Academy of Agriculture and Forestry Sciences, Beijing 100097, China

  • 2.

    School of Grassland Science, Beijing Forestry University, Beijing 100083, China

More Information
  • The study of sexual and evolutionary differences has long been imperative in the field of biology. Unlike animals, dioecious angiosperms account for only about 6% of the total. Buffalograss (Buchloe dactyloides) plays a vital role in environmental restoration, creating economic benefits and promoting the high-quality development of the grassland and turf industries. Its natural populations contain differing ratios of male, female, and monoecious plants. The value of buffalograss for studying the sex differentiation mechanism in plants cannot be ignored. However, few studies have investigated transcript annotation and complete mRNA structure in B. dactyloides, and the pathways of species-specific factors in sex differentiation remain unknown. We integrated the full-length transcriptome, second-generation transcriptome, and metabolome to specify candidate factors influencing sex differentiation. We identified 110,870 full-length transcripts and obtained 100,362 (90.52%) transcript and annotation information. Then we identified 49,448 differentially expressed genes (DEGs) and 3,070 differentially accumulated metabolites (DAMs) in female, male, and monoecious leaf samples. The co-enrichment analysis indicated that sexual differentiation was regulated by glutathione metabolism, photosynthesis, plant hormone biosynthesis, catabolism, and signaling. The identification of DEGs and DAMs that participate in glutathione metabolism, photosynthesis, abscisic acid (ABA), cytokinin (CTK), and gibberellin (GA) biosynthesis, catabolism, and signaling has helped illuminate the roles of plant hormones in the sex differentiation of B. dactyloides. The full-length transcriptomic data will facilitate additional studies on functional genes. Integration of transcriptomic and metabolomic data advances knowledge of the molecular mechanism of sex differentiation and provides information for B. dactyloides breeding programs.

  • Introduction

    With the development of the world economy, people's lifestyles have changed dramatically, and long-term high-intensity work has put many people's bodies in a sub-healthy state. The increasing incidence of various chronic diseases has not only put enormous pressure on society's healthcare systems but also caused endless suffering to people[1]. Therefore, people's demands on the functionality and safety of food are increasing, and it has become the consensus of people that 'not just eating enough, but more importantly eating well'.

    Rice is the staple food for more than half of the world's population and the main economic source for a large number of rural people[2]. However, due to the rising cost of rice cultivation, farmers are gaining less and less economic benefits from growing rice, which seriously undermines their incentive to grow rice and poses a serious threat to world food security. Increasing the added value of rice not only helps to increase farmers' income but also helps to ensure world food security. The presence of a large number of functional ingredients in rice makes it possible to increase the added value of rice, and functional rice has therefore been widely noticed.

    Functional rice refers to rice containing certain specific components that play a regulatory and balancing role in human physiological functions in addition to the nutrients necessary for human growth and development in the endosperm, embryo, and rice bran. They can increase human physiological defense mechanisms, prevent certain diseases, help recovery, delay aging, and boost physical strength and energy levels[3]. Rice is a staple food for more than half of the world's population[4], and its functional components have a great potential to be exploited for human welfare. Using functional rice as a carrier to address health problems and realize 'medicine-food homology' is an excellent motivation for promoting functional rice. The current typical functional rice is introduced in this paper. It also summarizes the breeding and cultivation technologies of functional rice.

    Types of functional rice

    High resistant starch rice

    Rice has a high glycemic index. Its long-term consumption leads to obesity, diabetes, and colon disease in many people[5]. However, the consumption of rice rich in resistant starch (RS) can greatly reduce the risk of these diseases[6]. Therefore, breeding rice varieties with high RS content has attracted considerable attention from breeders in various countries. However, the variability of RS content between different rice varieties is low, and there are few germplasm resources available for selection, thus making it challenging to breed rice varieties with high RS content using traditional breeding methods. Combining traditional and modern molecular breeding techniques can greatly improve the successful production of high RS rice breeds. Nishi et al.[7] selected a high RS rice variety EM10 by treating fertilized egg cells of Kinmaze with N-methyl-N-nitrosourea. However, its yield was very low, and it was not suitable for commercial production. Wada et al.[8] crossed 'Fukei 2032' and 'EM129' as parents and selected Chikushi-kona 85, a high RS rice variety with a higher yield than EM10. Miura et al.[9] bred ultra-high RS BeI-BEIIB double mutant rice by crossing the Abe I and Abe IIB mutant strains, and the content of RS in the endosperm reached 35.1%. Wei et al.[10] found that the simultaneous inhibition of starch branching enzyme (SBE) genes SBEIIb and SBEI in Teqing by antisense RNA could increase the RS content in rice to 14.9%. Zhu et al.[11] used RNAi technology to inhibit the expression of SBEI and SBEII genes in rice, which increased the content of RS in rice endosperm from 0 to 14.6 %. Zhou et al.[6] found that rice RS formation is mainly controlled by soluble starch synthase (SSIIA). However, its regulation is dependent on the granule-bound starch synthase Waxy (Wx), and SSIIA deficiency combined with high expression of Wxa facilitates the substantial accumulation of RS in the rice. The results of Tsuiki et al.[12] showed that BEIB deficiency was the main reason for the increased accumulation of RS in rice. Itoh et al.[13] developed new mutant rice lines with significantly higher levels of RS in rice by introducing genes encoding starch synthase and granule-bound starch synthase in the rice into the BEIB-deficient mutant line be2b.

    Colored rice

    The accumulation of anthocyanins/proanthocyanidins in the seed coat of the rice grain gives brown rice a distinct color[14]. Most common rice varieties lack anthocyanins in the seed coat, and so far, no rice variety with colored endosperm in its natural state has been identified. However, Zhu et al.[15] bred rice with purple endosperm using transgenic technology. Red rice contains only proanthocyanidins, while black and purple rice contain anthocyanidins and proanthocyanidins[16]. Red seed coat of rice was found to be controlled by the complementary effects of two central effect genes Rc and Rd. The loss of function of the Rc gene prevented the synthesis of proanthocyanidins, while the Rd gene could enhance the effect of the Rc gene in promoting proanthocyanidins synthesis[17]. Purple seed coat color is controlled by two dominant complementary genes Pb and Pp. Pb determines the presence or absence of seed coat color, and Pp determines the depth of seed coat color[18]. In addition, phycocyanin synthesis is also regulated by transcription factors such as MYB, bHLH, HY5, and WD40[14], but the exact regulatory mechanism is not clear. Colored rice is rich in bioactive components, such as flavonoids, phenolic acids, vitamin E (VE), glutelin, phytosterols, and phytic acid (PA). It also contains large amounts of micronutrients such as Ca, Fe, Zn, and Se[19], and has a much higher nutritional and health value than ordinary white rice. In addition, Zhu et al.[20] successfully developed rice with enriched astaxanthin in the endosperm by introducing the genes sZmPSY1, sPaCrtI, sCrBKT, and sHpBHY. This achievement has laid a solid foundation for the further development of functional rice industry.

    Giant embryo rice

    Giant embryo rice refers to rice varieties whose embryo volume is more than twice that of ordinary rice[21]. Rice embryo contains more nutrients than the endosperm; therefore, the nutritional value of giant embryo rice greatly exceeds that of ordinary rice. Studies have found that the levels of γ-aminobutyric acid (GABA), essential amino acids, VE, γ-oryzanol, phenols, and trace elements in giant embryo rice are considerably higher than that in ordinary rice[21]. Satoh & Omura[22] used the chemical mutagen N-methyl-N-nitrosourea to treat the fertilized egg cells of the rice variety Kinmaze to obtain a 'giant embryo' mutant. The mutants’ embryo occupied 1/4–1/3 of the rice grain volume and was 3–4 times larger than normal rice embryo[23]. Its GABA content increased dramatically after the rice was soaked in water. Maeda et al.[24] crossed the giant embryo mutant EM40 of Kinmaze with the high-yielding variety Akenohoshi to produce the giant embryo rice variety 'Haiminori'. The embryo size of 'Haiminori' is 3–4 times that of ordinary rice, and the GABA content of its brown rice is 3–4 times higher than that of 'Nipponbare' and 'Koshihikari' after soaking for four hours in water. A few genes that can regulate the size of rice embryos have been identified, and GE is the first identified rice giant embryo gene[25]. Nagasawa et al.[26] found that the loss of GE gene function resulted in enlarged embryos and smaller endosperm in rice. Lee et al.[27] found that the inhibition of LE gene expression by RNAi technology could lead to embryo enlargement in rice, but the regulatory mechanism remains to be investigated.

    Low glutelin rice

    Protein is the second most crucial nutrient in rice, accounting for 7–10% of the grain weight, and glutenin accounts for 60%–80% of the total protein content in rice grains[28]. Compared to other proteins, glutenin is more easily digested and absorbed by the body[29]. Therefore, higher glutenin content in rice can improve its nutritional value. However, people with renal disease (a common complication of diabetes) have impaired protein metabolism, and consumption of rice with lower glutelin content can help reduce their protein intake and metabolic burden[30]. Japanese breeders treated Nihonmasari with the chemical mutagen ethyleneimine and selected the low-glutelin rice mutant NM67[31]. Iida et al.[31] developed a new rice variety LGC-1 (Low glutelin content-1) with a glutelin content of less than 4% by backcrossing the NM67 mutant with the original variety 'Nihonmasari'. According to Miyahara[32], the low glutelin trait in LGC-1 is controlled by a single dominant gene Lgc-1 located on chromosome 2. Subsequently, Nishimura et al.[33] produced two rice varieties, 'LGC Katsu' and 'LGC Jun' with lower glutelin content by crossing LGC1 with a mutant line Koshikari (γ-ray induction) lacking 26 kDa globulin (another easily digestible protein).

    Golden rice

    Vitamin A (VA) is one of the essential nutrients for the human body[34]. However, rice, a staple food, lacks VA, leading to a VA deficiency in many people. β-carotene is a precursor for VA synthesis and can be effectively converted into VA in the human body[35]. Therefore, breeding rice varieties rich in β-carotene has attracted the attention of breeders in various countries. Ye et al.[36] simultaneously transferred phytoene synthase (psy), phytoene desaturase (crt I), and lycopene β-cyclase (lcy) genes into rice using the Agrobacterium-mediated method and produced the first generation of golden rice with a β-carotene content of 1.6 µg·g−1 in the endosperm. However, due to the low content of β-carotene in rice, it is difficult to meet the human body's demand for VA. To increase β-carotene content in rice, Paine et al.[37] introduced the phytoene synthase (psy) gene from maize and the phytoene desaturase (crt I) gene from Erwinia into rice. They obtained the second generation of golden rice with 37 µg g−1 of β-carotene in the endosperm, with nearly 23-fold increase in β-carotene content compared to the first generation of golden rice.

    Fe and Zn fortified rice

    Fe and Zn are essential trace elements for human beings. The contents of Fe and Zn in common rice are about 2 μg·g−1 and 16 μg·g−1, respectively[38], which are far from meeting human needs. In 2004, to alleviate micronutrient deficiencies among underprivileged people in developing countries, the Consultative Group on International Agricultural Research launched the HarvestPlus international collaborative program for improving Fe, Zn, and β-carotene levels in staple crops, with breeding targets of 13 μg·g−1 and 28 μg·g−1 for Fe and Zn in rice, respectively. Masuda et al.[39] found that expression of the nicotianamine synthase (NAS) gene HvNAS in rice resulted in a 3-fold increase in Fe and a 2-fold increase in Zn content in polished rice. Trijatmiko et al.[38] overexpressed rice OsNAS2 gene and soybean ferritin gene SferH-1 in rice, and the Fe and Zn content in polished rice of rice variety NASFer-274 reached 15 μg·g−1 and 45.7 μg·g−1, respectively. In addition, it has been found that increasing Fe intake alone does not eliminate Fe deficiency but also decreases the amount of Fe absorption inhibitors in the diet or increases the amount of Fe absorption enhancers[40]. The negatively charged phosphate in PA strongly binds metal cations, thus reducing the bioavailability of Fe and Zn in rice[41], while the sulfhydryl group in cysteine binds Fe, thereby increasing the absorption of non-heme Fe by the body[42]. To improve the bioavailability of Fe and Zn, Lucca et al.[40] introduced a heat-tolerant phytase (phyA) gene from Aspergillus fumigatus into rice and overexpressed the cysteine-rich protein gene (rgMT), which increased the content of phytase and cysteine residues in rice by 130-fold and 7-fold, respectively[40].

    Breeding technology of functional rice

    The functional quality of rice is highly dependent on germplasm resources. Current functional rice breeding mainly adopts transgenic and mutagenic technologies, and the cultivated rice varieties are mainly enriched with only one functional substance and cannot meet the urgent demand by consumers for rice enriched with multiple active components. The diversity of rice active components determines the complexity of multifunctional rice breeding. In order to cultivate multifunctional rice, it is necessary to strengthen the application of different breeding technologies. Gene polymerization breeding is a crop breeding technology that can polymerize multiple superior traits that have emerged in recent years, mainly including traditional polymerization breeding, transgenic polymerization breeding, and molecular marker-assisted selection polymerization breeding.

    Traditional polymerization breeding

    The transfer of beneficial genes in different species during traditional polymeric breeding is largely limited by interspecific reproductive isolation, and it is challenging to utilize beneficial genes between different species effectively. Gene transfer through sexual crosses does not allow accurate manipulation and selection of a gene and is susceptible to undesirable gene linkage, and in the process of breed selection, multiple backcrosses are required[43]. Thus, the period of selecting target plants is long, the breeding cost is high, and the human resources and material resources are costly[44]. Besides, it is often difficult to continue the breakthrough after a few generations of backcrossing due to linkage drag. Thus, there are significant limitations in aggregating genes by traditional breeding methods[45].

    Transgenic polymerization breeding

    Transgenic technology is an effective means of gene polymerization breeding. Multi-gene transformation makes it possible to assemble multiple beneficial genes in transgenic rice breeding rapidly and can greatly reduce the time and workload of breeding[46]. The traditional multi-gene transformation uses a single gene transformation and hybridization polymerization method[47], in which the vector construction and transformation process is relatively simple. However, it is time-consuming, laborious, and requires extensive hybridization and screening efforts. Multi-gene-based vector transformation methods can be divided into two major categories: multi-vector co-transformation and multi-gene single vector transformation[47]. Multi-vector co-transformation is the simultaneous transfer of multiple target genes into the same recipient plant through different vectors. The efficiency of multi-vector co-transformation is uncertain, and the increase in the number of transforming vectors will increase the difficulty of genetic screening, resulting in a reduced probability of obtaining multi-gene co-transformed plants. Multi-gene single vector transformation constructs multiple genes into the T-DNA region of a vector and then transfers them into the same recipient plant as a single event. This method eliminates the tedious hybridization and backcrossing process and solves the challenges of low co-transformation frequency and complex integration patterns. It can also avoid gene loss caused by multi-gene separation and recombination in future generations[47]. The transgenic method can break through the limitations of conventional breeding, disrupt reproductive isolation, transfer beneficial genes from entirely unrelated crops to rice, and shorten the cycle of polymerizing target genes significantly. However, there are concerns that when genes are manipulated, unforeseen side effects may occur, and, therefore, there are ongoing concerns about the safety of transgenic crops[48]. Marker-free transgenic technology through which selective marker genes in transgenic plants can be removed has been developed. This improves the safety of transgenic crops, is beneficial to multiple operations of the same transgenic crop, and improves the acceptance by people[49].

    Molecular marker-assisted selection polymerization breeding

    Molecular marker-assisted selection is one of the most widely used rice breeding techniques at present. It uses the close linkage between molecular markers and target genes to select multiple genes directly and aggregates genes from different sources into one variety. This has multiple advantages, including a focused purpose, high accuracy, short breeding cycle, no interference from environmental conditions, and applicability to complex traits[50]. However, few genes have been targeted for the main effect of important agronomic traits in rice, and they are mainly focused on the regulation of rice plant type and the prevention and control of pests and diseases, and very few genes related to the synthesis of active components, which can be used for molecular marker-assisted selection are very limited. Furthermore, the current technical requirements and costs for analyzing and identifying DNA molecular markers are high, and the identification efficiency is low. This greatly limits the popularization and application of functional rice polymerization breeding. Therefore, to better apply molecular marker-assisted selection technology to breed rice varieties rich in multiple active components, it is necessary to construct a richer molecular marker linkage map to enhance the localization of genes related to functional substance synthesis in rice[51]. Additionally, it is important to explore new molecular marker technologies to improve efficiency while reducing cost.

    It is worth noting that the effects of gene aggregation are not simply additive. There are cumulative additive effects, greater than cumulative epistatic effects, and less than cumulative epistatic effects among the polymerization genes, and the effects are often smaller than the individual effect. Only with a clearer understanding of the interaction between different QTLs or genes can functional rice pyramiding breeding be carried out reasonably and efficiently. Except for RS and Se, other active components of rice mainly exist in the rice bran layer, and the content of active components in the endosperm, the main edible part, is extremely low. Therefore, cultivating rice varieties with endosperm-enriched active components have broad development prospects. In addition, because crops with high quality are more susceptible to pests and diseases[52], the improvement of rice resistance to pests and diseases should be considered during the polymerization breeding of functional rice.

    Cultivation and regulation technology of functional rice

    The biosynthesis of active components in rice is influenced by rice varieties but also depends on cultivation management practices and their growth environment.

    Effect of light on the accumulation of active components in rice

    Environmental conditions have a greater effect on protein content than genetic forces[53]. Both light intensity and light duration affect the synthesis and accumulation of active components in rice. Low light intensity in the early stage of rice growth is not conducive to the accumulation of glutelin in rice grains but favors the accumulation of amylose, while the opposite is true in the late stage of rice growth[54]. Low light intensity during the grain-filling period reduces the accumulation of total flavonoids in rice[55] and decreases Fe ions' movement in the transpiration stream and thereby the transport of Fe ions to rice grains[56]. An appropriate increase in light intensity is beneficial to the accumulation of flavonoids, anthocyanins, and Fe in rice, but the photostability of anthocyanins is poor, and too much light will cause oxidative degradation of anthocyanins[57]. Therefore, functional rice is best cultivated as mid-late rice, which would be conducive to accumulating active components in rice.

    Effect of temperature on the accumulation of active components in rice

    The temperature has a great influence on the synthesis of active components in rice. An appropriate increase in the temperature is beneficial to the accumulation of γ-oryzanol[58] and flavonoids[59] in rice. A high temperature during the grain-filling period leads to an increase in glutelin content in rice[60], but an increase in temperature decreases the total phenolic content[61]. The results regarding the effect of temperature on the content of PA in rice were inconsistent. Su et al.[62] showed that high temperatures during the filling period would increase the PA content, while Goufo & Trindade[61] reported that the increase in temperature would reduce the PA content. This may be due to the different growth periods and durations of temperature stress on rice in the two studies. The synthesis of anthocyanins/proanthocyanidins in colored rice requires a suitable temperature. Within a certain range, lower temperatures favor the accumulation of anthocyanins/proanthocyanidins in rice[63]. Higher temperatures will lead to degradation, and the thermal stability of proanthocyanidins being higher than that of anthocyanins[64]. In addition, cold or heat stress facilitates GABA accumulation in rice grains[65]. Therefore, in actual production, colored rice and low-glutelin rice are best planted as late rice, and the planting time of other functional rice should be determined according to the response of its enriched active components to temperature changes.

    Effect of water management on the accumulation of active components in rice

    Moderate water stress can significantly increase the content of glutelin[66] and GABA[67] in rice grains and promote the rapid transfer of assimilation into the grains, shorten the grain filling period, and reduce the RS content[68]. Drought stress can also induce the expression of the phytoene synthase (psy) gene and increase the carotenoid content in rice[69]. Soil moisture is an important medium in Zn diffusion to plant roots. In soil with low moisture content, rice roots have low available Zn, which is not conducive to enriching rice grains with Zn[70]. Results from studies on the effect of soil water content on Se accumulation in rice grains have been inconsistent. Li et al.[71] concluded that flooded cultivation could significantly increase the Se content in rice grains compared to dry cultivation. However, the results of Zhou et al.[72] showed that the selenium content in rice grains under aerobic and dry-wet alternative irrigation was 2.44 and 1.84 times higher than that under flood irrigation, respectively. This may be due to the forms of selenium contained in the soil and the degree of drought stress to the rice that differed between experiments[73]. In addition, it has been found that too much or too little water impacts the expression of genes related to anthocyanin synthesis in rice, which affects the accumulation of anthocyanins in rice[74]. Therefore, it is recommended to establish different irrigation systems for different functional rice during cultivation.

    Effect of fertilizer application on the accumulation of active components in rice

    Both the amount and method of nitrogen application affect the accumulation of glutelin. Numerous studies have shown that both increased and delayed application of nitrogen fertilizer can increase the accumulation of lysine-rich glutelin to improve the nutritional quality of rice (Table 1). However, this improvement is not beneficial for kidney disease patients who cannot consume high glutelin rice. Nitrogen stress can down-regulate the expression of ANDs genes related to the anthocyanins biosynthesis pathway in grains, resulting in a decrease in anthocyanins synthesis[55]. Increased nitrogen fertilizer application can also increase the Fe, Zn, and Se content in rice[75,76]. However, some studies have found that increased nitrogen fertilizer application has no significant effect on the Fe content of rice[77], while other studies have shown that increased nitrogen fertilizer application will reduce the Fe content of rice[78]. This may be influenced by soil pH and the form of the applied nitrogen fertilizer. The lower the soil pH, the more favorable the reduction of Fe3+ to Fe2+, thus promoting the uptake of Fe by rice. Otherwise, the application of ammonium fertilizer can improve the availability of soil Fe and promote the absorption and utilization of Fe by rice. In contrast, nitrate fertilizer can inhibit the reduction of Fe3+ and reduce the absorption of Fe by rice[79].

    Table 1.  Effect of nitrogen fertilizer application on glutelin content of rice.
    SampleN level
    (kg ha−1)    		Application timeGlutelin content
    (g 100 g−1)    		References
    Rough rice05.67[66]
    270Pre-transplanting : mid tillering : panicle initiation : spikelet differentiation = 2:1:1:16.92
    300Pre-transplanting : mid tillering : panicle initiation : spikelet differentiation = 5:2:2:16.88
    Brown rice05.35[83]
    90Pre-transplanting : after transplanting = 4:16.01
    Pre-transplanting : after transplanting = 1:16.60
    180Pre-transplanting : after transplanting = 4:16.53
    Pre-transplanting : after transplanting = 1:17.29
    270Pre-transplanting : after transplanting = 4:17.00
    Pre-transplanting : after transplanting = 1:17.66
    Rough rice05.59[84]
    187.5Pre-transplanting : after transplanting = 4:16.47
    Pre-transplanting : after transplanting = 1:16.64
    300Pre-transplanting : after transplanting = 4:17.02
    Pre-transplanting : after transplanting = 1:17.14
    Polished rice03.88[85]
    90Pre-transplanting : tillering : booting = 2:2:14.21
    180Pre-transplanting : tillering : booting = 2:2:14.43
    270Pre-transplanting : tillering : booting = 2:2:16.42
    360Pre-transplanting : tillering : booting = 2:2:14.87
    Brown rice09.05[86]
    120Flowering22.14
     | Show Table
    DownLoad: CSV

    Appropriate application of phosphorus fertilizer is beneficial in promoting the translocation of Fe and Zn from leaves to rice grains, thus increasing the content in rice grains[80]. However, the excessive application of phosphate fertilizer will reduce the availability of Fe and Zn in soil, resulting in less uptake by the roots and a lower content in the rice grains[81]. The content of PA in rice increased with a higher phosphorus fertilizer application rate[80]. Increasing the phosphorus fertilizer application rate would increase the values of [PA]/[Fe] and [PA]/[Zn] and reduce the effectiveness of Fe and Zn in rice[80]. Currently, there are few studies on the effect of potassium fertilization on the synthesis of active components in rice. Available studies report that increased application of nitrogen fertilizer can increase the Zn content in rice[82]. Therefore, the research in this area needs to be strengthened.

    Because the iron in soil mainly exists in the insoluble form Fe3+, the application of iron fertilizer has little effect on rice biofortification[87]. There are different opinions about the effect of Zn fertilizer application methods. Phattarakul et al.[88] believed that foliar spraying of Zn fertilizer could significantly improve the Zn content in rice grains. Jiang et al.[89] concluded that most of the Zn accumulated in rice grains were absorbed by the roots rather than from the reactivation of Zn in leaves. In contrast, Yuan et al.[90] suggested that soil application of Zn fertilizer had no significant effect on Zn content in rice grains. The different results may be affected by the form of zinc fertilizer applied and the soil conditions in the experimental sites. Studies have found that compared with the application of ZnEDTA and ZnO, zinc fertilizer in the form of ZnSO4 is most effective for increasing rice's Zn[70]. In addition, the application of zinc fertilizer reduces the concentration of PA in rice grains[70].

    The form of selenium fertilizer and the method and time of application will affect the accumulation of Se in rice grains. Regarding selenium, rice is a non-hyperaccumulative plant. A moderate application of selenium fertilizer can improve rice yield. However, the excessive application can be toxic to rice, and the difference between beneficial and harmful supply levels is slight[91]. Selenite is readily adsorbed by iron oxide or hydroxide in soil, and its effectiveness in the soil is much lower than selenite[92]. In addition, selenate can migrate to the roots and transfer to rice shoots through high-affinity sulfate transporters. In contrast, selenite is mainly assimilated into organic selenium in the roots and transferred to the shoots in smaller amounts[93]. Therefore, the biological effectiveness of Se is higher in selenate-applied soil than in selenite application[94] (Table 2). Zhang et al.[95] found that the concentration of Se in rice with soil application of 100 g Se ha-1 was only 76.8 μg·kg-1, while the concentration of Se in rice with foliar spray of 75 g Se ha-1 was as high as 410 μg·kg-1[73]. However, the level of organic selenium was lower in rough rice with foliar application of selenium fertilizer compared to soil application[96], while the bioavailability of organic selenium in humans was higher than inorganic selenium[97]. Deng et al.[73] found that the concentrations of total selenium and organic selenium in brown rice with selenium fertilizer applied at the full heading stage were 2-fold higher than those in brown rice with selenium fertilizer applied at the late tillering stage (Table 2). Although the application of exogenous selenium fertilizer can rapidly and effectively increase the Se content of rice (Table 2), it can easily lead to excessive Se content in rice and soil, which can have adverse effects on humans and the environment. Therefore, breeding Se-rich rice varieties is a safer and more reliable way to produce Se-rich rice. In summary, functional rice production should include the moderate application of nitrogen and phosphorus fertilizer and higher levels of potassium fertilizer, with consideration to the use of trace element fertilizers.

    Table 2.  Effect of selenium fertilizer application on the selenium content of rice.
    SampleSe level (g Se ha−1)Selenium fertilizer formsApplication methodSe content (μg·g−1)References
    Rough rice00.002[98]
    18SeleniteFoliar spray at full heading0.411
    Polished rice00.071[99]
    20SeleniteFoliar spray at full heading0.471
    20SelenateFoliar spray at full heading0.640
    Rough rice75SeleniteFoliar spray at late tillering0.440[73]
    75SeleniteFoliar spray at full heading1.290
    75SelenateFoliar spray at late tillering0.780
    75SelenateFoliar spray at full heading2.710
    Polished rice00.027[100]
    15SeleniteFoliar spray at full heading0.435
    45SeleniteFoliar spray at full heading0.890
    60SeleniteFoliar spray at full heading1.275
     | Show Table
    DownLoad: CSV

    Effect of harvesting time on the accumulation of active components in rice

    The content of many active components in rough rice is constantly changing during the development of rice. It was found that the content of total flavonoids in brown rice increased continuously from flowering stage to dough stage and then decreased gradually[101]. The γ-oryzanol content in rice decreased by 13% from milk stage to dough stage, and then gradually increased to 60% higher than milk stage at full maturity[101]. The results of Shao et al.[102] showed that the anthocyanin content in rice reached its highest level at two weeks after flowering and then gradually decreased. At full ripeness, and the anthocyanins content in brown rice was only about 50% of the maximum level. The content of total phenolics in rice decreased with maturity from one week after flowering to the fully ripe stage, and the loss of total phenolics reached more than 47% by the fully ripe stage. In contrast, the content of total phenolics in black rice increased with maturity[102]. Moreover, RS content in rough rice decreases during rice maturation[68]. Therefore, the production process of functional rice should be timely and early harvested to obtain higher economic value.

    Prevention and control of pests and diseases for functional rice

    Pests and diseases seriously impact the yield and quality of rice[103]. At present, the two most effective methods to control pests and diseases are the use of chemical pesticides and the planting of pest and disease-resistant rice varieties. The use of chemical pesticides has greatly reduced the yield loss of rice. However, excessive use of chemical pesticides decreases soil quality, pollutes the environment, reduces soil biodiversity[104], increases pest resistance, and aggravates the adverse effects of pests and diseases on rice production[105]. It also increases residual pesticide levels in rice, reduces rice quality, and poses a severe threat to human health[106].

    Breeding pest and disease-resistant rice varieties are among the safest and effective ways to control rice pests and diseases[107]. In recent years, many pest and disease resistance genes from rice and microorganisms have been cloned[47]. Researchers have used these genes to breed rice varieties resistant to multiple pests and diseases through gene polymerization breeding techniques. Application in production practices delivered good ecological and economic benefits[108].

    Green pest and disease control technologies must consider the synergies between rice and water, fertilizer, and pest and disease management. In this regard, the rice-frog, rice-duck, and other comprehensive rice production models that have been widely used in recent years are the most representative. These rice production models significantly reduced chemical pesticide usage and effectively controlled rice pests and diseases[109]. The nutritional imbalance will reduce the resistance of rice to pests and diseases[110]. Excessive application of nitrogen fertilizer stimulates rice overgrowth, protein synthesis, and the release of hormones, increasing its attractiveness to pests[111]. Increased soluble protein content in rice leaves is more conducive to virus replication and increases the risk of viral infection[112]. Increasing the available phosphorus content in the soil will increase crop damage by pests[113], while insufficient potassium supply will reduce crop resistance to pests and diseases[114]. The application of silica fertilizer can boost the defense against pests and diseases by increasing silicon deposition in rice tissue, inducing the expression of genes associated with rice defense mechanisms[115] and the accumulation of antifungal compounds in rice tissue[116]. The application of silica fertilizer increases the release of rice volatiles, thereby attracting natural enemies of pests and reducing pest damage[117]. Organic farming increases the resistance of rice to pests and diseases[118]. In addition, rice intercropping with different genotypes can reduce pests and diseases through dilution and allelopathy and changing field microclimate[119].

    In conclusion, the prevention and control of rice pests and diseases should be based on chemical and biological control and supplemented by fertilizer management methods such as low nitrogen, less phosphorus, high potassium and more silicon, as well as agronomic measures such as rice-aquaculture integrated cultivation, organic cultivation and intercropping of different rice varieties, etc. The combined use of multiple prevention and control measures can improve the yield and quality of functional rice.

    Conclusions

    Functional rice contains many active components which are beneficial to maintaining human health and have high economic and social value with broad market prospects. However, the current development level of the functional rice industry is low. The development of the functional rice requires extensive use of traditional and modern polymerization breeding techniques to cultivate new functional rice varieties with endosperm that can be enriched with multiple active components and have broad-spectrum resistance to pests and diseases. It is also important to select suitable planting locations and times according to the response characteristics of different functional rice active components to environmental conditions.

    Acknowledgments

    This work is supported by the National Natural Science Foundation of China (Project No. 32060430 and 31971840), and Research Initiation Fund of Hainan University (Project No. KYQD(ZR)19104).

  • The authors declare that they have no conflict of interest.

  • Supplemental Table S1 Primers used for qRT-PCR validation.
    Supplemental Table S2 Prediction of Simple sequence repeats (SSRs) out of transcript datasets.
    Supplemental Table S3 List of targets of lncRNA in B. dactyloides.
    Supplemental Table S4 Summary of AS events in B. dactyloides.
    Supplemental Table S5 Illumaina Sequencing Data Assessment.
    Supplemental Table S6 Differentially expressed genes identified between different sample groups.
    Supplemental Table S7 Differentially accumulated metabolites identified between different sample groups.
    Supplemental Fig. S1 Summary of SMRT sequencing. (A) CCS read length distribution. (B) FLNC sequences read length distribution. (C) Consensus isoforms read length distribution.
    Supplemental Fig. S2 The quality evaluation of transcriptomic and metabolomic data. PCA analysis of all transcripts (A) and metabolites (B), respectively. Heatmap analysis of correlation between samples in all transcripts (C) and metabolites (D), respectively.
    Supplemental Fig. S3 KEGG enrichment of the up-regulated and down-regulated DEGs identified in Male vs Female (A), Monoecious vs Female (B), and Monoecious vs Male (C), respectively.
    Supplemental Fig. S4 Prediction of regulatory proteins. A Number and family of top ten regulatory proteins predicted by SMRT sequencing data. B Venn diagram of differentially expressed regulatory proteins. C Number and family of top ten differentially expressed regulatory proteins. D-F The differentially expressed regulatory proteins predicted in Male vs Female, Monoecious vs Female, and Monoecious vs Male, respectively.
    Supplemental Fig. S5 GO enrichment of sex-specific genes. (A) Genes present only in female plants. (B) Genes absent in female plants. (C) Genes present only in male plants. (D) Genes absent in male plants. (E) Genes present only in monoecious plants. (F) Genes absent in monoecious plants.
    Supplemental Fig. S6 KEGG enrichment of sex-specific genes. (A) Genes present only in female plants. (B) Genes absent in female plants. (C) Genes present only in male plants. (D) Genes absent in male plants. (E) Genes present only in monoecious plants. (F) Genes absent in monoecious plants.
    Supplemental Fig. S7 KEGG enrichment analysis of DAMs in the positive (left) and negative (right) ion mode in the comparison of Male vs Female (A), Monoecious vs Female (B), and Monoecious vs Male (C) sample groups.
  • [1]

    Liao Q, Du R, Gou J, Guo L, Shen H, et al. 2020. The genomic architecture of the sex-determining region and sex-related metabolic variation in Ginkgo biloba. The Plant Journal 104:1399−409

    doi: 10.1111/tpj.15009

    CrossRef   Google Scholar

    [2]

    Massonnet M, Cochetel N, Minio A, Vondras AM, Lin J, et al. 2020. The genetic basis of sex determination in grapes. Nature communications 11:2902

    doi: 10.1038/s41467-020-16700-z

    CrossRef   Google Scholar

    [3]

    Marshall C, Warnke S, Amundsen K. 2022. Simple sequence repeat marker development and diversity analysis in buffalograss. Crop Science 62:1373−82

    doi: 10.1002/csc2.20725

    CrossRef   Google Scholar

    [4]

    Johnson PG, Kenworthy KE, Auld DL, Riordan TP. 2001. Distribution of buffalograss polyploid variation in the southern Great Plains. Crop Science 41:909−13

    doi: 10.2135/cropsci2001.413909x

    CrossRef   Google Scholar

    [5]

    Johnson PG, Riordan TP, Johnson-Cicalese JJ. 2000. Low-mowing tolerance in buffalograss. Crop Science 40:1339−43

    doi: 10.2135/cropsci2000.4051339x

    CrossRef   Google Scholar

    [6]

    Budak H, Shearman RC, Parmaksiz I, Gaussoin RE, Riordan TP, et al. 2004. Molecular characterization of buffalograss germplasm using sequence-related amplified polymorphism markers. Theoretical and Applied Genetics 108:328−34

    doi: 10.1007/s00122-003-1428-4

    CrossRef   Google Scholar

    [7]

    Johnson PG, Riordan TP, Arumuganathan K. 1998. Ploidy level determinations in buffalograss clones and populations. Crop Science 38:478−82

    doi: 10.2135/cropsci1998.0011183X003800020034x

    CrossRef   Google Scholar

    [8]

    Stark R, Grzelak M, Hadfield J. 2019. RNA sequencing: the teenage years. Nature Reviews Genetics 20:631−56

    doi: 10.1038/s41576-019-0150-2

    CrossRef   Google Scholar

    [9]

    Yan H, Sun M, Zhang Z, Jin Y, Zhang A, et al. 2023. Pangenomic analysis identifies structural variation associated with heat tolerance in pearl millet. Nature Genetics 55:507−18

    doi: 10.1038/s41588-023-01302-4

    CrossRef   Google Scholar

    [10]

    Guan J, Yin S, Yue Y, Liu L, Guo Y, et al. 2022. Single-molecule long-read sequencing analysis improves genome annotation and sheds new light on the transcripts and splice isoforms of Zoysia japonica. BMC Plant Biology 22:263

    doi: 10.1186/s12870-022-03640-7

    CrossRef   Google Scholar

    [11]

    Teng K, Teng W, Wen H, Yue Y, Guo W, et al. 2019. PacBio single-molecule long-read sequencing shed new light on the complexity of the Carex breviculmis transcriptome. BMC Genomics 20:789

    doi: 10.1186/s12864-019-6163-6

    CrossRef   Google Scholar

    [12]

    Liu L, Teng K, Fan X, Han C, Zhang H, et al. 2022. Combination analysis of single-molecule long-read and Illumina sequencing provides insights into the anthocyanin accumulation mechanism in an ornamental grass, Pennisetum setaceum cv. Rubrum. Plant Molecular Biology 109:159−75

    doi: 10.1007/s11103-022-01264-x

    CrossRef   Google Scholar

    [13]

    Wenger AM, Peluso P, Rowell WJ, Chang PC, Hall RJ, et al. 2019. Accurate circular consensus long-read sequencing improves variant detection and assembly of a human genome. Nature Biotechnology 37:1155−62

    doi: 10.1038/s41587-019-0217-9

    CrossRef   Google Scholar

    [14]

    Jozefczuk S, Klie S, Catchpole G, Szymanski J, Cuadros-Inostroza A, et al. 2010. Metabolomic and transcriptomic stress response of Escherichia coli. Molecular Systems Biology 6:364

    doi: 10.1038/msb.2010.18

    CrossRef   Google Scholar

    [15]

    Bradley PH, Brauer MJ, Rabinowitz JD, Troyanskaya OG. 2009. Coordinated concentration changes of transcripts and metabolites in Saccharomyces cerevisiae. PLoS Computational Biology 5:e1000270

    doi: 10.1371/journal.pcbi.1000270

    CrossRef   Google Scholar

    [16]

    Sun M, Yan H, Zhang A, Jin Y, Lin C, et al. 2023. Milletdb: a multi-omics database to accelerate the research of functional genomics and molecular breeding of millets. Plant Biotechnology Journal 21:2348−57

    doi: 10.1111/pbi.14136

    CrossRef   Google Scholar

    [17]

    Sharon D, Tilgner H, Grubert F, Snyder M. 2013. A single-molecule long-read survey of the human transcriptome. Nature Biotechnology 31:1009−14

    doi: 10.1038/nbt.2705

    CrossRef   Google Scholar

    [18]

    Li W, Godzik A. 2006. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22:1658−59

    doi: 10.1093/bioinformatics/btl158

    CrossRef   Google Scholar

    [19]

    Hackl T, Hedrich R, Schultz J, Förster F. 2014. proovread: large-scale high-accuracy PacBio correction through iterative short read consensus. Bioinformatics 30:3004−11

    doi: 10.1093/bioinformatics/btu392

    CrossRef   Google Scholar

    [20]

    Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, et al. 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nature Protocols 8:1494−512

    doi: 10.1038/nprot.2013.084

    CrossRef   Google Scholar

    [21]

    Liu X, Mei W, Soltis PS, Soltis DE, Barbazuk WB. 2017. Detecting alternatively spliced transcript isoforms from single-molecule long-read sequences without a reference genome. Molecular Ecology Resources 17:1243−56

    doi: 10.1111/1755-0998.12670

    CrossRef   Google Scholar

    [22]

    Beier S, Thiel T, Münch T, Scholz U, Mascher M. 2017. MISA-web: a web server for microsatellite prediction. Bioinformatics 33:2583−85

    doi: 10.1093/bioinformatics/btx198

    CrossRef   Google Scholar

    [23]

    Zheng Y, Jiao C, Sun H, Rosli HG, Pombo MA, et al. 2016. iTAK: a program for genome-wide prediction and classification of plant transcription factors, transcriptional regulators, and protein kinases. Molecular Plant 9:1667−70

    doi: 10.1016/j.molp.2016.09.014

    CrossRef   Google Scholar

    [24]

    Kong L, Zhang Y, Ye Z, Liu X, Zhao S, et al. 2007. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Research 35:W345−W349

    doi: 10.1093/nar/gkm391

    CrossRef   Google Scholar

    [25]

    Sun L, Luo H, Bu D, Zhao G, Yu K, et al. 2013. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Research 41:e166

    doi: 10.1093/nar/gkt646

    CrossRef   Google Scholar

    [26]

    Wang L, Park HJ, Dasari S, Wang S, Kocher JP, et al. 2013. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Research 41:e74

    doi: 10.1093/nar/gkt006

    CrossRef   Google Scholar

    [27]

    Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, et al. 2016. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Research 44:D279−D285

    doi: 10.1093/nar/gkv1344

    CrossRef   Google Scholar

    [28]

    Li J, Ma W, Zeng P, Wang J, Geng B, et al. 2015. LncTar: a tool for predicting the RNA targets of long noncoding RNAs. Briefings in Bioinformatics 16:806−12

    doi: 10.1093/bib/bbu048

    CrossRef   Google Scholar

    [29]

    Boeckmann B, Bairoch A, Apweiler R, Blatter MC, Estreicher A, et al. 2003. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucleic Acids Research 31:365−70

    doi: 10.1093/nar/gkg095

    CrossRef   Google Scholar

    [30]

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, et al. 2000. Gene ontology: tool for the unification of biology. Nature Genetics 25:25−29

    doi: 10.1038/75556

    CrossRef   Google Scholar

    [31]

    Tatusov RL, Galperin MY, Natale DA, Koonin EV. 2000. The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Research 28:33−36

    doi: 10.1093/nar/28.1.33

    CrossRef   Google Scholar

    [32]

    Koonin EV, Fedorova ND, Jackson JD, Jacobs AR, Krylov DM, et al. 2004. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biology 5:R7

    doi: 10.1186/gb-2004-5-2-r7

    CrossRef   Google Scholar

    [33]

    Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. 2012. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Research 40:D109−D114

    doi: 10.1093/nar/gkr988

    CrossRef   Google Scholar

    [34]

    Huerta-Cepas J, Szklarczyk D, Forslund K, Cook H, Heller D, et al. 2016. eggNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences. Nucleic Acids Research 44:D286−D93

    doi: 10.1093/nar/gkv1248

    CrossRef   Google Scholar

    [35]

    Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, et al. 2013. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29:15−21

    doi: 10.1093/bioinformatics/bts635

    CrossRef   Google Scholar

    [36]

    Bray NL, Pimentel H, Melsted P, Pachter L. 2016. Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology 34:525−27

    doi: 10.1038/nbt.3519

    CrossRef   Google Scholar

    [37]

    Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550

    doi: 10.1186/s13059-014-0550-8

    CrossRef   Google Scholar

    [38]

    Wang J, Zhang T, Shen X, Liu J, Zhao D, et al. 2016. Serum metabolomics for early diagnosis of esophageal squamous cell carcinoma by UHPLC-QTOF/MS. Metabolomics 12:116

    doi: 10.1007/s11306-016-1050-5

    CrossRef   Google Scholar

    [39]

    Thévenot EA, Roux A, Xu Y, Ezan E, Junot C. 2015. Analysis of the human adult urinary metabolome variations with age, body mass index, and gender by implementing a comprehensive workflow for univariate and OPLS statistical analyses. Journal of Proteome Research 14:3322−35

    doi: 10.1021/acs.jproteome.5b00354

    CrossRef   Google Scholar

    [40]

    Yu G, Wang L, Han Y, He Q. 2012. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology 16:284−87

    doi: 10.1089/omi.2011.0118

    CrossRef   Google Scholar

    [41]

    Yu G, Xie Z, Lei S, Li H, Xu B, et al. 2022. The NAC factor LpNAL delays leaf senescence by repressing two chlorophyll catabolic genes in perennial ryegrass. Plant Physiology 189:595−610

    doi: 10.1093/plphys/kiac070

    CrossRef   Google Scholar

    [42]

    Shearman RC, Riordan TP, Johnson PG. 2004. Buffalograss. Warm‐Season (C4) Grasses 45:1003−26

    doi: 10.2134/agronmonogr45.c31

    CrossRef   Google Scholar

    [43]

    Lee JT. 2012. Epigenetic regulation by long noncoding RNAs. Science 338:1435−39

    doi: 10.1126/science.1231776

    CrossRef   Google Scholar

    [44]

    Di C, Yuan J, Wu Y, Li J, Lin H, et al. 2014. Characterization of stress-responsive lncRNAs in Arabidopsis thaliana by integrating expression, epigenetic and structural features. The Plant Journal 80:848−61

    doi: 10.1111/tpj.12679

    CrossRef   Google Scholar

    [45]

    Reddy ASN, Marquez Y, Kalyna M, Barta A. 2013. Complexity of the alternative splicing landscape in plants. The Plant Cell 25:3657−83

    doi: 10.1105/tpc.113.117523

    CrossRef   Google Scholar

    [46]

    Wang M, Wang P, Liang F, Ye Z, Li J, et al. 2018. A global survey of alternative splicing in allopolyploid cotton: landscape, complexity and regulation. New Phytologist 217:163−78

    doi: 10.1111/nph.14762

    CrossRef   Google Scholar

    [47]

    Zhou F, Xu D, Xiong S, Chen C, Liu C, et al. 2023. Transcriptomic and metabolomic profiling reveal the mechanism underlying the inhibition of wound healing by ascorbic acid in fresh-cut potato. Food Chemistry 410:135444

    doi: 10.1016/j.foodchem.2023.135444

    CrossRef   Google Scholar

    [48]

    Xu D, Lin H, Tang Y, Huang L, Xu J, et al. 2021. Integration of full-length transcriptomics and targeted metabolomics to identify benzylisoquinoline alkaloid biosynthetic genes in Corydalis yanhusuo. Horticulture Research 8:16

    doi: 10.1038/s41438-020-00450-6

    CrossRef   Google Scholar

    [49]

    Schmitt B, Vicenzi M, Garrel C, Denis FM. 2015. Effects of N-acetylcysteine, oral glutathione (GSH) and a novel sublingual form of GSH on oxidative stress markers: a comparative crossover study. Redox Biology 6:198−205

    doi: 10.1016/j.redox.2015.07.012

    CrossRef   Google Scholar

    [50]

    Knoop H, Gründel M, Zilliges Y, Lehmann R, Hoffmann S, et al. 2013. Flux balance analysis of cyanobacterial metabolism: the metabolic network of Synechocystis sp. PCC 6803. PLoS Computational Biology 9:e1003081

    doi: 10.1371/journal.pcbi.1003081

    CrossRef   Google Scholar

    [51]

    Urano K, Maruyama K, Jikumaru Y, Kamiya Y, Yamaguchi-Shinozaki K, et al. 2017. Analysis of plant hormone profiles in response to moderate dehydration stress. The Plant Journal 90:17−36

    doi: 10.1111/tpj.13460

    CrossRef   Google Scholar

    [52]

    Sun TP. 2008. Gibberellin metabolism, perception and signaling pathways in Arabidopsis. The Arabidopsis Book 6:e0103

    doi: 10.1199/tab.0103

    CrossRef   Google Scholar

    [53]

    Kieber JJ, Schaller GE. 2014. Cytokinins. The Arabidopsis Book 12:e0168

    doi: 10.1199/tab.0168

    CrossRef   Google Scholar

  • Cite this article

    Guan J, Yue Y, Yin S, Teng W, Zhang H, et al. 2023. Integration analysis of full-length transcriptomics and metabolomics provides new insights into the mechanism of sex differentiation in buffalograss (Buchloe dactyloides). Grass Research 3:24 doi: 10.48130/GR-2023-0024
    Guan J, Yue Y, Yin S, Teng W, Zhang H, et al. 2023. Integration analysis of full-length transcriptomics and metabolomics provides new insights into the mechanism of sex differentiation in buffalograss (Buchloe dactyloides). Grass Research 3:24 doi: 10.48130/GR-2023-0024
    shu

Figures(9)

Download PDF

Article Metrics

Article views(3265) PDF downloads(250)

Other Articles By Authors

ARTICLE   Open Access    

Integration analysis of full-length transcriptomics and metabolomics provides new insights into the mechanism of sex differentiation in buffalograss (Buchloe dactyloides)

  • 1.

    Institute of Grassland, Flowers, and Ecology, Beijing Academy of Agriculture and Forestry Sciences, Beijing 100097, China

  • 2.

    School of Grassland Science, Beijing Forestry University, Beijing 100083, China

  • Received: 25 September 2023
  • Accepted: 02 November 2023
  • Published online: 05 December 2023
Grass Research  3,  Article number: 24  (2023)  |  Cite this article

Abstract: The study of sexual and evolutionary differences has long been imperative in the field of biology. Unlike animals, dioecious angiosperms account for only about 6% of the total. Buffalograss (Buchloe dactyloides) plays a vital role in environmental restoration, creating economic benefits and promoting the high-quality development of the grassland and turf industries. Its natural populations contain differing ratios of male, female, and monoecious plants. The value of buffalograss for studying the sex differentiation mechanism in plants cannot be ignored. However, few studies have investigated transcript annotation and complete mRNA structure in B. dactyloides, and the pathways of species-specific factors in sex differentiation remain unknown. We integrated the full-length transcriptome, second-generation transcriptome, and metabolome to specify candidate factors influencing sex differentiation. We identified 110,870 full-length transcripts and obtained 100,362 (90.52%) transcript and annotation information. Then we identified 49,448 differentially expressed genes (DEGs) and 3,070 differentially accumulated metabolites (DAMs) in female, male, and monoecious leaf samples. The co-enrichment analysis indicated that sexual differentiation was regulated by glutathione metabolism, photosynthesis, plant hormone biosynthesis, catabolism, and signaling. The identification of DEGs and DAMs that participate in glutathione metabolism, photosynthesis, abscisic acid (ABA), cytokinin (CTK), and gibberellin (GA) biosynthesis, catabolism, and signaling has helped illuminate the roles of plant hormones in the sex differentiation of B. dactyloides. The full-length transcriptomic data will facilitate additional studies on functional genes. Integration of transcriptomic and metabolomic data advances knowledge of the molecular mechanism of sex differentiation and provides information for B. dactyloides breeding programs.

    HTML

    Introduction

    • Unlike animals, dioecious angiosperms account for only about 6% of the total.[1]. In many cases, dioecious plants provide better tools for understanding the genetics of sex determination than dioecious animals[2]. Buffalograss (Buchloe dactyloides (Nutt.) Engelm.) is a perennial, stoloniferous, low herbaceous, and low-input ecological grass[3]. Buffalograss is a dioecious plant, and a few are hermaphroditic. Buffalograss comprises polyploid polymorphisms (x = n = 10), including diploid, tetraploid, pentaploid, and hexaploid individuals[4]. The turf quality of the female plant in B. dactyloides is better than that of the male as when the male plant enters the flowering stage, there is a distinct yellow color due to the flower shaft being higher than the bush, and the female plant does not have such a problem[5]. It was reported that there are no sex-determining heteromorphic chromosomes, and sex determination is linked to gene expression and gibberellin-induced morphological alterations in buffalograss[4]. Studies on B. dactyloides have focused on variety breeding and improvement[6,7], physiological differentiation in response to biotic and abiotic stress, and phenotypic identification. However, the transcriptional information of B. dactyloides, especially in complete mRNA sequences, alternative splicing (AS) events, and long non-encoding RNAs (lncRNAs), is still largely unknown. Moreover, the large number of unannotated transcripts limits the research on the mechanisms of sex differentiation of B. dactyloides.

      Plants have complex transcriptional processes; transcriptome research is critical to understanding plant life processes[8]. The RNA-seq technology, based on Illumina sequencing platforms, has significantly higher sequencing depth and a reduced error rate, resulting in substantially higher genome coverage[9], allowing us to better comprehend critical gene expression levels and regulatory mechanisms[10]. However, RNA-seq data cannot systematically and accurately collect or assemble complete transcripts or identify transcripts of isoform, homologous, superfamily, or allele expression, making it insufficient to understand the transcriptome and underlying biological processes. Instead of interrupting RNA fragments, full-length transcriptomes employing PacBio single-molecule real-time (SMRT) sequencing platforms can directly reverse full-length cDNA. Ultra-long reads (median 10 kb) from the platform contain a single complete transcript sequence[11,12]. Assembly is not required for post-analysis, and what is measured is[13]. Moreover, the response to metabolites is more specific than the general response seen at the transcript level[14]. At present, comprehensive investigations of alterations in metabolite levels remain rare. It is particularly true for the joint and parallel analysis of sex differentiation at two levels of genome information processing: the transcriptomic and the metabolomic. However, the integrative analysis of metabolites and transcripts helps identify genes and small molecules involved in critical biological processes[15,16]. We completed the full-length transcriptome sequencing analysis using the transcriptome analysis methods. We explored the complete mRNA structure features of B. dactyloides and completed the identified transcript function annotation. We further investigated the underlying transcriptional mechanism of B. dactyloides sex differentiation by integrating Illumina RNA-seq data and metabolomic data. This work will effectively provide helpful information about the transcriptome of B. dactyloides. We expect it will offer essential resources for further research on sex differentiation, functional gene identification, and variety breeding studies in B. dactyloides.

    Materials and methods

      Plant materials and RNA sample preparation

    • Buchloe dactyloides (cv. Texoka) were vegetatively propagated in an intelligent greenhouse at the Beijing Academy of Agricultural and Forestry Sciences (Beijing, China). In April 2022, we harvested female, male, and monoecious plants' roots, stolon, flowers, and leaves as materials. All samples were frozen in liquid nitrogen and stored at −80 °C for further investigations. To identify as many transcripts as possible, we extracted the total RNA from all organ materials of buffalograss for SMRT sequencing. However, in practical application, we completed sex identification before blooms. We selected leaves of female, male, and monoecious plants for Illumina and metabolite sequencing. First, we extracted total RNA using the plant RNA kit (Omega, GA, USA). Second, the purity, concentration, and integrity of total RNA were assessed using ND-1000 spectrophotometer (NanoDrop Technologies, DE, USA) and 2100 bioanalyzer (Agilent Technologies, CA, USA).

      Full-length library construction and PacBio SMRT sequencing

    • We mixed RNA from female, male, and monoecious plants' roots, stolons, flowers, and leaves into one sample. Following the qualification of the sample, a full-length library was constructed with no size selection. After the full-length library passed quality control, we performed full-length transcriptome sequencing based on the PB PacBio Sequel II platform. Then full-length transcriptome sequencing data were analyzed from raw reads in three steps[17]. First, we extracted and polished circular consensus (CCS) sequences using full passes ≥ three and sequence accuracy > 0.9. Next, we classified the full-length transcripts by detecting whether 5' and 3' cDNA primers and poly(A) tails in CCS sequences. Finally, we classified full-length non-chimeric (FLNC) transcripts to recognize consensus isoforms, high-quality isoforms (accuracy > 0.99), and low-quality isoforms using the SMRTLink software. We used CD-HIT v4.6.1 software[18] to merge highly similar sequences and remove redundancy. The low-quality consensus isoforms were polished based on the corresponding Illumina sequencing data using proovread software[19].

      Full-length transcriptome structure analysis

    • We conducted the structure analysis with non-redundant FLNC transcripts identified by the full-length transcriptome, including CDS prediction, AS events, and SSR identification and analysis. We used TransDecoder (v5.0.0) software[20] to predict the reliable CDS and their corresponding protein sequences of the transcript. The AS events prediction method referred to the report following the three alignments[21]. We screened transcripts with a length > 1,000 bp and completed SSR analysis utilizing MISA software[22].

      Identification of regulatory proteins, lncRNAs and lncRNA target

    • We used iTAK software[23] with default parameters to predict regulatory proteins, including transcription factors (TFs), transcriptional regulators (TRs), and protein kinases (PKs) in B. dactyloides. We filtered out the putative protein-coding RNAs before lncRNAs prediction. We used CPC[24], CNCI[25], CPAT[26], and Pfam[27], to screen > 200 bp transcripts in length that contain two or more exons. The intersection of the four sets was applied in the following lncRNA analysis. In addition, RNA targets of lncRNAs were predicted using LncTar software[28].

      Function annotation of transcripts

    • To acquire function annotation information, we compared the non-redundant transcripts to non-redundant protein sequences (NR), Swissprot[29], Gene Ontology (GO)[30], Clusters of Orthologous proteins Groups (COG)[31], Eukaryotic Ortholog Groups (KOG)[32], Pfam, Kyoto Encyclopedia of Genes and Genomes (KEGG)[33], and eggNOG[34] databases utilizing DIAMOND software. Moreover, the protein sequences of transcripts were blasted (p-value < 1E-5) against the TrEMBL database based on the Swissprot database to obtain the annotation information.

      Illumina cDNA library construction and sequencing

    • Nine cDNA libraries (three biological replicates for female, male, and monoecious plants) were constructed for Illumina sequencing. After verifying the quality of nine libraries, we used the Illumina NovaSeq6000 platform (San Diego, CA, USA) to perform Illumina sequencing. The specific sequencing data quality controls were: cut the sequencing adapter and the primer sequence in reads; filtrate the low-quality reads to ensure the read quality. Finally, we obtained clean data that is sufficient for subsequent accurate analysis. We further calculated GC content, Q30, and sequence duplication levels of clean data.

      Quantification of transcript expression and analysis of differential expression transcripts

    • We used non-redundant transcript sequences identified from PacBio SMRT as a reference and applied NGS reads to them to acquire transcript position information using STAR[35] software. The method of quantification of transcript expression levels was as follows: first, we used Kallisto[36] software to compare NGS reads with transcripts generated from PacBio SMRT and count directly; then, we calculated transcript expression levels using fragments per kilobase of transcript per million fragments mapped (FPKM). To evaluate the reproducibility of biological replicates, we set the Pearson Correlation Coefficient (R2) as the correlation index of samples.

      Then we used DESeq2[37] software to identify differential expressions with replicates. We detected DEGs using the threshold: fold change (FC) of not less than two and a false discovery rate (FDR) of less than 0.01. We adjusted the p-value using the Benjamini-Hochberg method. Volcano Plot assessed the differences in transcript expression between two samples and their corresponding significance. Moreover, GO and KEGG functional annotation and enrichment analysis on DEGs we performed.

      qRT-PCR analysis

    • DEGs involved in abscisic acid (ABA) metabolic and signaling pathways (BdABA1, BdAOG, and BdSAPK8), gibberellin (GA) metabolic and signaling pathways (BdKAO, BdGID1, and BdDWARF8), cytokinin (CTK) related pathways gene (BdORR9), photosynthesis-related pathways (BdPsbP which encodes PSII reaction center subunit and PetC which encodes cytochrome b6/f complex subunits), and one bHLH transcription factor bHLH137 related to sex differentiation were subjected to qRT-PCR analysis. We reverse-transcribe RNA to cDNA using PrimeScriptTM Reverse Transcriptase (TaKaRa, Dalian, China). The primers used in the study were designed using Primer Premier 5 and listed in Supplemental Table S1. qRT-PCR analysis was performed on a Bio-Rad CFX Connect™ RealTime System using SYBR® Premix Ex Taq™ II (TaKaRa, Dalian, China). Three biological replicates were used for all gene expression analyses. We calculated the relative expression level of genes using the 2−ΔΔCᴛ method.

      Metabolites extraction and LC-MS/MS analysis

    • Three experimental groups (male, female, and monoecious plants) contained six samples per group harvested for sequencing and separate analyses. We performed metabolomics experiments using the Waters UPLC Acquity Ⅰ-Class PLUS and UPLC Xevo G2-XS QTOF. Then we analyzed samples using the Waters Acquity UPLC HSS T3 column (1.8 um, 2.1 mm × 100 mm) in positive and negative codes[38]. The mobile phases were 0.1% formic acid aqueous solution (A) and 0.1% formic acid acetonitrile (B), and the injection volume was 1.0 μL. We acquired primary and secondary mass spectrometry data in MSe mode using MassLynx V4.2 software (Waters). Dual-channel data acquisition, including low and high collision energies, can be simultaneously accomplished during each data acquisition cycle. The low and high collision energies were 2 V and 10–40 V, respectively, and the scanning frequency was 0.2 s for a mass spectrum. The parameters of the ESI ion source were according to a previously described report[38]. We processed raw data for peak extraction and alignment and completed identification based on the online METLIN and Biomark's self-built databases using Progenesis QI software. In addition, we identified the theoretical fragment; the mass deviation of the parent ion and fragment ion is less than 100 ppm and 50 ppm, respectively[38].

      Quantification of metabolites and analysis of differential metabolites

    • We used the R package (ropls) to perform orthogonal projections to latent structures-discriminant analysis (OPLS-DA) modeling between groups, and we completed 200 permutation tests to ensure the reliability of the model[39]. We determined the model's variable importance in projection (VIP) value using multiple cross-validations. We compared the difference multiples based on the grouping and used the T-test to calculate the difference significance p-value of metabolites. We further identified DAMs using the threshold: p-value less than 0.05 and VIP value greater than one. The hypergeometric distribution test was used to calculate the DAMs of KEGG pathway enrichment significance[40].

      Data preprocessing and annotation, and correlation analysis between DEGs and DAMs

    • After normalizing the original peak area information with the total peak area, we evaluate the overall quality of the data. We evaluated the overall data quality using principal component analysis (PCA) and Spearman correlation analysis. We searched the classification and pathway information of identified compounds utilizing the KEGG[33] database.

      To calculate the correlation between all DEGs and DAMs, we preprocessed the data using the z-value transformation method and screened according to the correlation coefficient (CC) and p-value of correlation. The screening threshold was |CC| greater than 0.80 and CCP less than 0.05. We mapped DEGs and DAMs simultaneously to related pathways using MapMan software[41].

    Results

      SMRT sequencing analysis and function annotation of transcripts

    • We accomplished the PacBio SMRT sequencing of B. dactyloides and obtained 45.10 Gb of clean data. The PacBio SMRT reads generated in this study were submitted to the BioProject database of the National Center for Biotechnology Information (accession numbers PRJNA991319). The cDNA size of the full-length library was 1–6 kb. First, 675,685 CCS sequences were acquired from raw reads. The average read length of CCS was 1,691 bp, and the mean number of passes was 31. Next, 83.26% (562,568) of CCS sequences were identified as FLNC sequences. Moreover, FLNC sequences were clustered using the SMRTLink software, and 179,367 consensus isoforms were identified. To improve the accuracy of the sequences, we used proovread software to correct the low-quality consensus isoforms corresponding to Illumina sequencing data. 179,248 polished high-quality isoforms were identified, and the percentage of polished high-quality isoforms was 99.93%. The mean read length of consensus isoform was 1,667 bp. The read length distributions of the CCS sequence, FLNC sequence, and consensus isoform are shown in Supplemental Fig. S1. Finally, we obtained 110,870 non-redundant transcripts generated using CD-HIT software.

      Using DIAMOND software, we compared 110,870 non-redundant transcript sequences with all databases and obtained 100,362 (90.52%) transcript annotation information (Fig. 1a). We compared the transcript sequences with the NR database. We identified 99,534 homologous transcripts (Fig. 1b). NR homologous species result showed that the most abundant transcripts (49,236) were distributed in Eragrostis curvula, accounting for 49.47%. GO annotation and enrichment analysis divided 82,896 transcripts into three categories (Fig. 1c). First, the most abundant transcripts were involved in the 'biological process' term. Transcripts involved in the 'biological process' were associated with the cellular process (41,643), metabolic process (38,118), and biological regulation (13,703). Then, transcripts involved in the 'cellular component' were related to cellular anatomical entity (45,140) and intracellular (24,922). Finally, for 'molecular function', transcripts were mainly involved in binding (42,921) and catalytic activity (38,620). We identified 136 KEGG pathways involved in 69,000 transcripts using the KEGG database. We listed 13 KEGG pathways involving more than 1,000 transcripts in Fig. 1d. The results showed that plant-pathogen interaction (3,058), plant hormone signal transduction (2,349), carbon metabolism (2,096), biosynthesis of amino acids (1,537), starch and sucrose metabolism (1,535), and MAPK signaling pathway-plant (1,501) were the six abundant pathways for transcripts.

      Figure 1. 

      Function annotation of non-redundant transcripts. (a) Function annotation of transcripts in all databases. (b) Nr homologous species distribution diagram of transcripts. (c) Distribution of GO terms for all annotated transcripts. Transcripts (82,896) were classified into three groups: 'biological process', cellular component', and 'molecular function'. (d) Thirteen KEGG pathways involving more than 1,000 transcripts.

      Full-length transcriptome structure analysis

    • A total of 101,909 ORFs and 56,478 carried complete ORFs. One hundred to 200 amino acid transcripts were the most abundant, accounting for 23.92%. 94.98% of the transcript length was more than 800 bp. The length distribution of amino acid sequences in the coding region of ORFs is shown in Fig. 2a. We analyzed the over 1,000 bp transcripts to predict SSR using MISA software. Six types of SSRs, totaling 32,942 SSRs, were identified (Supplemental Table S2). Most of these SSRs were mono- and tri-nucleotide repeats, accounting for 47.1% and 31.1%, respectively. The least abundant SSRs were hexa- and penta-nucleotide repeats, accounting for 0.2% and 0.3%, respectively.

      Figure 2. 

      Structure analysis of full-length transcripts. (a) Length distribution of amino acid sequences in the coding region of ORFs. (b) Venn diagram of lncRNAs predicted by CPC, CPAT, CNCI, and Pfam methods. (c) GO and (d) KEGG enrichment of the transcripts possessed AS events.

      For lncRNA prediction, we used the method of filtering step by step: first, the intersection of CPAT and CPC prediction, then CNCI prediction, and finally, Pfam to predict CNCI results. We predicted 7,682 lncRNAs; the Venn diagram is shown in Fig. 2b. We predicted the targets of 6,474 lncRNAs using LncTar (Supplemental Table S3), accounting for 84.27%. Using BLAST software, we identified 4,091 AS events involved in 6,188 transcripts (Supplemental Table S4). To discover the critical biological processes and functions of transcripts in B. dactyloides, we completed GO and KEGG annotation and enrichment analysis of the transcripts presented in AS events. A total of 2,557 transcripts were classified into three groups: 'biological processes', 'cellular components', and 'molecular function'. We satisfied the top 20 enriched pathways (Fig. 2c), and the top five most enriched pathways were 'pyruvate, phosphate dikinase activity', 'pyruvate metabolic process, regulation of proteolysis', 'ubiquitin-like protein ligase binding', and 'regulation of photoperiodism, flowering'. Using the KEGG database, we identified 5,833 transcript annotations, including 128 KEGG pathways. KEGG pathway enrichment analysis showed that the three most enriched pathways were 'glycosaminoglycan degradation', 'Glycosphingolipid biosynthesis – ganglio series', and 'Spliceosome' (Fig. 2d).

      Acquisition and mining of transcriptome data

    • The Illumina NGS reads generated in this study were submitted to the BioProject database of the National Center for Biotechnology Information (accession number PRJNA991464). Leaf samples of female, male, and monoecious plants were determined for transcriptome analysis to clarify the sex differentiation mechanism in B. dactyloides (Fig. 3a−c). Using the Illumina NovaSeq6000 platform, we obtained 191.62 million clean reads and 53.44 GB of clean data from nine samples. The GC content was more than 53.36%, and the Q30 value was more than 94.05% (Supplemental Table S5). We identified 110,870 transcript expression profiles. We evaluated the quality of transcriptomic data using PCA and Spearman rank correlation. 110,870 transcripts were divided into three groups, with good repeatability within each group based on PCA results (Supplemental Fig. S2a). The correlation diagram showed that the correlation values in the female group are > 0.918; these values are > 0.977 in the male group and > 0.975 in the monoecious group (Supplemental Fig. S2b). The above results indicated that the sequencing data was sufficient for DEG analysis. We identified 49,448 DEGs using DESeq2 software (Fig. 3d). There were 31,174 DEGs (15,670 up-regulated and 15,504 down-regulated), 30,789 DEGs (15,992 up-regulated and 14,797 down-regulated), and 31,882 DEGs (16,265 up-regulated and 15,617 down-regulated) identified in the groups 'male vs female', 'monoecious vs female', and 'monoecious vs male', respectively (Supplemental Table S6).

      Figure 3. 

      (a)–(c) Differences in the appearance and the number of (d) DEGs and (e) DAMs between female, male, and monoecious sample groups in B. dactyloides. Red boxes indicate the flower of B. dactyloides. (f) Association analysis between RPKM ratio and relative expression levels of the qRT-PCR. Data represent the mean ± SD (n = 3). Bars represent RNA-seq data, and lines represent qRT-PCR data. The lower-case letters above the bars indicate a significant difference between the sample groups at p < 0.05.

      We measured the expression levels of ten DEGs by qRT-PCR, including three ABA-related pathway genes (ABA1, AOG, and SnRK2), three GA-related pathways genes (KAO, GID1, and DELLA), two photosynthesis-related genes (PsbP and PetC), one CTK-related pathways gene (ARR), and one bHLH transcription factor. The qRT-PCR results validated the RNA-seq data (Fig. 3f).

      We performed KEGG enrichment analysis of up-and down-regulated DEGs to reveal further the sex differentiation mechanism in B. dactyloides (Supplemental Fig. S3). For the group 'male vs female', 'synthesis and degradation of ketone bodies', 'photosynthesis – antenna proteins', and 'circadian rhythm – plant' were activated, while 'biosynthesis of various secondary metabolites – part 2', 'circadian rhythm – plant', and spliceosome were inhibited. For the group 'monoecious vs female', 'photosynthesis – antenna proteins', 'porphyrin and chlorophyll metabolism', 'circadian rhythm – plant' were activated, while diterpenoid biosynthesis, glutathione metabolism, arginine and proline metabolism were inhibited. For the group 'monoecious vs male', photosynthesis – antenna proteins, porphyrin and chlorophyll metabolism, and photosynthesis were activated, while isoflavonoid biosynthesis, monoterpenoid biosynthesis, and linoleic acid metabolism were inhibited.

      Regulatory proteins dynamics during B. dactyloides sex differentiation

    • We predicted 9,916 regulatory proteins from 212 families using the iTAK prediction tool, including 4,331 TFs, 1,469 TRs, and 4,116 PKs. The top ten most numerous regulatory protein families are shown in Supplemental Fig. S4a. The top five regulatory protein families were RLK-Pelle_DLSV (520), bHLH (388), Others (309), MYB-related (300), and bZIP (252). Further, we analyzed the overlapping differentially expressed regulatory proteins between male vs female, monoecious vs female, and monoecious vs male. The result showed 688 differentially expressed regulatory proteins from 129 families in three groups (Supplemental Fig. S4b). The top ten most numerous regulatory protein families are shown in Supplemental Fig. S4c. RLK-Pelle_DLSV (62), others (31), NAC (26), bHLH (25), and MYB-related (19) were the top five regulatory protein families.

      Furthermore, we concentrated on the expression patterns of regulatory protein of sex differentiation in B. dactyloides. For the group 'male vs female', 3,224 differentially expressed regulatory proteins were from 177 families (Supplemental Fig. S4d), including 1,719 up-regulated regulatory proteins from 166 families and 1,505 down-regulated regulatory proteins from 162 families. For the group 'monoecious vs female', there were 3,473 differentially expressed regulatory proteins from 179 families (Supplemental Fig. S4e), including 1,647 up-regulated regulatory proteins from 165 families and 1,826 down-regulated regulatory proteins from 158 families. For the abundant RLK-Pelle_DLSV regulatory protein family, there were 198 down-regulated regulatory proteins, whereas 80 up-regulated regulatory proteins. The number of down-regulated regulatory proteins in the WRKY family was more than seven times that of the up-regulated regulatory proteins. There were 3,515 differentially expressed regulatory proteins from 180 families for the group 'monoecious vs male' (Supplemental Fig. S4f), including 1,575 up-regulated regulatory proteins from 161 families and 1,940 down-regulated regulatory proteins from 160 families. For the abundant RLK-Pelle_DLSV regulatory protein family, there were 193 down-regulated regulatory proteins, whereas 57 up-regulated regulatory proteins. The number of down-regulated regulatory proteins in the WRKY family was more than ten times that of up-regulated regulatory proteins.

      Sex-specific genes analysis

    • We found variations in gene expression among the sexes. Female plants had 588 genes relative to male and monoecious plants that deleted 2,020 genes. In addition, 1,961 genes were absent from female and monoecious plants vs male plants, and 526 genes were present only in male plants. Monoecious plants had 344 genes relative to female and male plants that deleted 1,629 genes. To identify the gene function related to plant sexual differentiation, we performed GO and KEGG enrichment of sex-specific genes (Supplemental Figs S5 & S6). Genes present only in female plants play significant roles in ADP binding (47) and defense response (41). Genes absent in female plants mainly participate in ADP binding (87) and defense response (76). Genes present only in male plants mainly participate in ADP binding (20). Genes absent in male plants are primarily involved in ADP binding (90), defense response (75), polysaccharide binding (18), rDNA binding (6), and 5S rDNA binding (6). Genes present only in monoecious plants mainly participate in ADP binding (19) and defense response (17). Genes absent in monoecious plants are the main participants in ADP binding (43), defense response (23), modification-dependent protein catabolic process (16), and dioxygenase activity (11).

      In addition, for genes present only in female plants, the top three enriched pathways were plant-pathogen interaction, fatty acid elongation, and pentose phosphate. The top three enriched pathways for genes absent in female plants were phenylalanine, tyrosine and tryptophan biosynthesis, fructose and mannose metabolism, and tyrosine metabolism. Sesquiterpenoid and triterpenoid biosynthesis, ribosome, and isoquinoline alkaloid biosynthesis were the top three enriched pathways for genes present only in female plants. The top three enriched pathways for genes absent in male plants were glutathione metabolism, alpha-Linolenic acid metabolism, and peroxisome. Citrate cycle (TCA cycle), Riboflavin metabolism, and alpha-Linolenic acid metabolism were the top three enriched pathways for genes present only in monoecious plants. Alpha-Linolenic acid metabolism, fatty acid biosynthesis, and carbon fixation in photosynthetic organisms were the top three enriched pathways for genes absent in monoecious plants.

      Acquisition and mining of metabolomic data

    • The metabolomics data generated in this study were submitted to MetaboLights database (accession number MTBLS8147). We scanned 10,011 peaks and 1,966 metabolites in positive code, 11,130 peaks, and 1,787 metabolites through qualitative and quantitative metabolome analysis of 18 samples using the LC-QTOF platform (Fig. 3e). We evaluated the quality of metabolomic data using PCA and Spearman rank correlation. A total of 3,753 metabolites were divided into three groups, with good repeatability within each group based on PCA results (Supplemental Fig. S2c). In addition, the correlation diagram showed that the correlation values in the female group are > 0.945; these values are > 0.922 in the male group and > 0.953 in the monoecious group (Supplemental Fig. S2d). We scanned 485 and 516 (positive and negative) up-regulated DAMs, as well as 369 and 339 down-regulated DAMs in the male vs female comparison. We observed 727 and 670 (positive and negative, the same below) up-regulated DAMs and 558 and 488 down-regulated DAMs in the monoecious vs female comparison. In the monoecious vs male comparison, we identified 691 and 626 up-regulated DAMs and 559 and 515 down-regulated DAMs (Supplemental Table S7). We performed KEGG enrichment pathway analysis based on identified DAMs. The top 20 enriched metabolic pathways in different groups are shown in Supplemental Fig. S7. For the group male vs female, the top four enriched pathways of DAMs (positive) were 'one carbon pool by folate', 'glutathione metabolism', 'pyrimidine metabolism', and 'anthocyanin biosynthesis'. The top two enriched pathways (negative) were flavone and flavonol biosynthesis and monobactam biosynthesis (Supplemental Fig. S7a). For the group 'monoecious vs female, the top four pathways (positive) were terpenoid backbone biosynthesis, glycerophospholipid metabolism, arachidonic acid metabolism, and carotenoid biosynthesis, while arachidonic acid metabolism was most enriched pathways (negative) (Supplemental Fig. S7b). For the group 'monoecious vs male', 'arginine and proline metabolism' was the top enriched pathways (positive), whereas 'phenylalanine, tyrosine and tryptophan biosynthesis', 'aflatoxin biosynthesis', and 'tyrosine metabolism' were the top three pathways (negative) (Supplemental Fig. S7c).

      To further clarify the metabolic changes in response to sex differentiation, we completed KEGG pathway-based analysis with metabolites in female, male, and monoecious plants combined with positive and negative codes. We focused on the pathways involved in at least five DAMs and calculated differential abundance scores. Then we captured the tendency for DAMs in a pathway to increase or decrease than another group (Fig. 4). Among the 46 pathways identified in the group 'male vs female', three were up (> 0.5 differential abundance score, red, the same below), including bisphenol degradation, benzoate degradation, and protein digestion and absorption. During the 60 pathways found in the group 'monoecious vs female', three were up, including protein digestion and absorption, betalain biosynthesis, and central carbon metabolism in cancer. During the 71 pathways in the group 'monoecious vs male', five were up, including 'phosphotransferase system (PTS)', 'protein digestion and absorption', 'phenylalanine, tyrosine, and tryptophan biosynthesis', 'valine, leucine, and isoleucine degradation', and 'biosynthesis of various alkaloids'.

      Figure 4. 

      Pathway-based analysis of metabolic changes on different groups. The differential abundance score captures the average, gross changes for all metabolites in a pathway. A score of 1 indicates all measured metabolites in the pathway increase, and –1 indicates all measured metabolites in a pathway decrease.

      Correlation analysis of transcriptomic and metabolomic data

    • To further identify metabolic pathways associated with sex differentiation, we completed co-enrichment analyses of DEGs and DAMs (Fig. 5). Six significantly enriched pathways (p < 0.01) in the sample group (male vs female) included biosynthesis of amino acids and various secondary metabolites – part 2, carbon metabolism, photosynthesis, riboflavin metabolism, and Vitamin B6 metabolism (Fig. 5a). Seven significantly enriched pathways (p < 0.01) in sample group (monoecious vs female) included alpha-Linolenic acid metabolism, biosynthesis of amino acids and various secondary metabolites – part 2, carbon, glutathione, glyoxylate, dicarboxylate, porphyrin, and chlorophyll metabolism (Fig. 5b). Eleven significantly enriched pathways (p < 0.01) in sample group (monoecious vs female) included alpha-Linolenic acid metabolism, butanoate metabolism, glutathione metabolism, isoflavonoid biosynthesis, monoterpenoid biosynthesis, porphyrin, chlorophyll, propanoate, and riboflavin metabolism, sesquiterpenoid and triterpenoid biosynthesis, tyrosine metabolism, and Vitamin B6 metabolism (Fig. 5c).

      Figure 5. 

      Co-enrichment p-value histogram of DEGs and DAMs. (a) Male vs female. (b) Monoecious vs female. (c) Monoecious vs male.

      Since these pathways are associated with glutathione metabolism, photosynthesis, and plant hormone biosynthesis and signal transduction, we mapped DEGs and DAMs simultaneously to related pathways using MapMan software better to understand the sex differentiation mechanisms in B. dactyloides.

      Glutathione metabolism and photosynthetic capacity changes in response to sex differentiation

    • We draw the networks of DEGs and DAMs implicated in the glutathione metabolism pathway using MapMan software (Fig. 6a) to evaluate the impact of glutathione metabolism among sexes. The results showed that 154 genes and two metabolites were differentially expressed among the three groups. The glutamate content in monoecious plants was lower relative to male and female plants. Oxidized glutathione (GSSG) levels in female plants were lower than in male plants. GCL (encoding for glutamate-cysteine ligase catalytic subunit) and DHAR (encoding glutathione dehydrogenase/transferase) were less expressed in male plants relative to female and monoecious plants. In addition, the expression of GSS (encoding glutathione synthase) and GSR (encoding glutathione reductase) was lower in female plants than in male and monoecious plants. Many GSTs encoding glutathione S-transferase showed the highest in monoecious plants, followed by female plants, and the lowest in male plants. The results showed that the DEGs and DAMs related to glutathione metabolism were involved in response to sexual differentiation.

      Figure 6. 

      Comparison of DEGs and DAMs involved in (a) glutathione metabolism and (b) photosynthesis pathways in response to sex differentiation. The red and blue text indicated DAM and DEG, respectively. The red color means up-regulation, while the blue color means down-regulation.

      The DEGs and DAMs involved in photosynthesis are shown in Fig. 6b. 78.34% of DEGs in photosynthesis showed higher expression levels in male plants relative to monoecious plants, including PsbA, PsbB, PsbP, PsbR, PsbW, Psb27 proteins in photosystem Ⅱ (PSⅡ), PetH and PetJ responsible for photosynthetic electron transport, PsaB, PsaE, PsaG, and PsaO proteins in photosystem Ⅰ (PSⅠ), and the gamma and delta complexes of F-type ATPase. 73.38% DEGs in photosynthesis showed higher expression levels in female plants relative to monoecious plants, including PsbC, PsbB, PsbK, PsbO, PsbP, PsbS, PsbW, and Psb28 proteins in PSⅡ, PetA encoding cytochrome b6/f complex (Cyt b6f) subunits, PsaD, PsaG, and PsaO proteins in PSⅠ, and the delta and b sub-complexes of F-type ATPase. In addition, the expression level of photosynthesis-related genes varied in female and male plants. The PsbA, PetJ, PsaB, and delta content were higher, whereas PsbC, PsbK, PetA, and PsaA were less in female plants than in male plants. Phosphoric acid (Pi), ATP, and NADP were related to carbon fixation in photosynthetic organisms and varied among sexes. The Pi levels showed the highest in male plants, then female plants, and the lowest in monoecious plants. In female plants, the content of ATP was higher in monoecious plants, and the content of NADP was higher in male plants.

      Plant hormones in response to sex differentiation

    • We examined the changes in hormone metabolites in female, male, and monoecious plants and found that abscisic acid (ABA), cytokinin (CTK), and gibberellin (GA) were differentially expressed. DEGs and DAMs involved in ABA, CTK, and GA metabolism and signaling pathways were mapped to related pathways using MapMan software to better understand the mechanism of plant hormones in response to sex differentiation.

      We analyzed the levels of DAMs and DEGs encoding enzymes that participate in ABA metabolism in female, male, and monoecious plants (Fig. 7). Antheraxanthin content and the expression of DEGs encoding zeaxanthin epoxidase (ABA1/ZEP) were highest in female plants, then male plants, and lowest in monoecious plants. The neoxanthin levels were lower than in female and male plants. Compared with male plants, three out of five DEGs encoding 9-cis-epoxycarotenoid dioxygenases (NCEDs) were down-regulated in female plants. Compared with monoecious plants, four out of five DEGs encoding NCEDs were up-regulated in female plants; this value was eight out of nine in male plants. The number of overexpression DEGs encoding short-chain alcohol dehydrogenase (ABA2) and abscisic aldehyde oxidase (AAO3) was highest in monoecious plants, followed by female plants, and lowest in male plants. Compared with monoecious plants, the abscisic alcohol, ABA, and ABA glucosyl ester (ABA-GE) levels were lower than in female and male plants. 8' -hydroxy ABA was slightly down-regulated in male plants than in monoecious plants. Then we focus on the ABA catabolism enzymes abscisate beta-glucosyltransferases (AOGs) and Cytochrome P450 (CYP707A). Compared with male plants, two out of three DEGs encoding AOGs were up-regulated in female plants. Compared with monoecious plants, two DEGs encoding AOGs were up-regulated in female plants, and one DEGs encoding AOGs were up-regulated in male plants. Compared with male plants, three out of five DEGs encoding CYP707A were down-regulated in female plants. Compared with monoecious plants, three out of four DEGs encoding AOGs were up-regulated in female or male plants.

      Figure 7. 

      Comparison of DEGs and DAMs involved in ABA biosynthetic and signaling pathways in response to sex differentiation. The red and blue text indicated DAMs and DEGs, respectively. The red color means up-regulation, while the blue color means down-regulation.

      The expression patterns of ABA signaling genes were further analyzed. ABA acceptor protein PYR/PYL, Type 2C protein phosphatase (PP2C), SNF1-related protein kinase 2 (SnRK2), and transcriptional factor ABF were all differentially expressed. 71.25% and 55.22% of DEGs showed higher expression levels in male and female plants relative to monoecious plants, respectively. Furthermore, 71.61% of DEGs showed lower expression levels in female plants than in male plants.

      We analyzed the GA biosynthesis, catabolism, and signal transduction pathways involving 275 DEGs and seven DAMs (Fig. 8). The expression of KS is less in the monoecious plant than in the female plant. Of the genes in the GA biosynthesis pathway, KAO and GA20ox are the most noteworthy. The essential genes, KAO and GA20ox, are more active in the male plant than in the female or monoecious plant (adjusted p-value < 0.01). PB_transcript_123979, encoding GA3ox, was expressed only in the monoecious plant. GA catabolism genes GA2ox also exist with specific expression. PB_transcript_514 was expressed only in the female plant. In contrast, PB_transcript_85857 was expressed only in male and monoecious plants. Compared with monoecious plants, about 80% of transcripts encoding GID1 (GA receptor) proteins were down-regulated in female plants, and this value was 75% in male plants. For GA signal repressors, DELLAs, there was little difference in the number of up-down and down-down DEGs among the three groups. Compared with phytochrome interacting transcription factors (PIFs) interacting with DELLA in monoecious plants, there were twice as many up-regulated DEGs as down-regulated DEGs.

      Figure 8. 

      Comparison of DEGs and DAMs involved in GA biosynthetic and signaling pathways in response to sex differentiation. The red and blue text indicated DAM and DEG, respectively. The red color means up-regulation, while the blue color means down-regulation. GGDP, geranylgeranyl diphosphate; CPS, ent-copalyl diphosphate synthase; KS, ent-kaurene synthase; KO, ent-kaurene oxidase; KAO, ent-kaurene acid oxidase; GA20ox, GA 20-oxidase; GA3ox, GA 3-oxidase; GA2ox, GA 2-oxidase; GID1, GA INSENSITIVE DWARF1.

      We examined the expression of transcripts encoding enzymes and metabolites involved in CTK metabolism in female, male, and monoecious plants (Fig. 9). The expression of three critical enzymes involved in CTK biosynthesis, adenosine phosphate-isopentenyltransferases (IPTs), was significantly reduced in female and male plants compared with monoecious plants. The substrate synthesizing iPRPs, ATP, exhibited a higher level in male plants than in monoecious plants. Similarly, one of the iPRPs, the iPRMP in male plants, was higher than that in female and monoecious plants. The primary forms of CTK in plants are isopentenyladenine (iP), trans-zeatin (tZ), and dihydrozeatin (DZ). We found that iP exhibited a lower content in female plants, and tZ exhibited a lower level in male plants. Compared with monoecious plants, six out of seven DEGs encoding CTK oxidase (CKX) were up-regulated in female plants, while this value was six out of nine in male plants. UDP-glucose, conjugated to active CTKs, either tZ or DZ, showed a high content in male plants. The metabolites of CTK metabolism, UDP, trans-Zeatin-O-glucoside (tZOG), and O-Xylosylzeatin (OXZ), exhibited a higher level in monoecious plants than in male and female plants. The content of UDP in male plants was higher than in female plants. Likewise, the expression of DEG-encoding enzymes (UGT73C and UGT85A1) conjugated to CTKs in glucose to inactivate them was highest in monoecious plants, then male plants, and lowest in female plants. Transcripts encoding cytokinin response 1 (CRE1), histidine-containing phosphotransfer proteins (AHPs), and type-B response regulators (B-ARRs) were differentially expressed in female, male, and monoecious plants. Then we focused on the negative regulator of CTK signaling, the type-A response regulator (A-ARRs). Compared with male plants, 64.7% of DEGs encoding A-ARRs were up-regulated in female plants. Compared with monoecious plants, 57.9% of DEGs encoding A-ARRs were up-regulated in female plants, whereas this value was 33.3% in male plants.

      Figure 9. 

      Comparison of DEGs and DAMs involved in CTK biosynthetic and signaling pathways in response to sex differentiation. The red and blue text indicate DAM and DEG, respectively. The red color means up-regulation, while the blue color means down-regulation.

    Discussion

    • Buffalograss is a perennial warm-season grass species of Gramineae used commonly as pasture, ecological, or lawn grass[42]. Buffalograss plays a vital role in environmental restoration, creating economic benefits and promoting the high-quality development of the grassland and turf industries[3]. Buffalograss is a dioecious plant, and a few are monoecious. The visual and functional quality of the female plant in B. dactyloides is better than the male plant because when the male plant enters the flowering stage, there is a distinct yellow color due to the flower shaft being higher than the bush, and the female plant does not have such a problem[5]. There are both dioecious and monoecious plants in buffalograss, providing great convenience for studying sex differentiation. The identification of sex differentiation candidate genes (e.g., male sterility) is an essential step to proclaim the practical value of B. dactyloides and may provide necessary information for sex differentiation-related pathway elucidation.

      The missing genome sequence and full-length cDNAs of B. dactyloides limit our research. This study first provided the full-length transcriptome of B. dactyloides and completed 90.52% of transcript function annotation using PacBio SMRT technology. Recent studies have shown that long non-coding RNAs (lncRNAs) play critical functions in various biological processes[43,44]. lncRNAs could interact with proteins, RNAs, and DNA to perform complex and diverse functions[28]. Therefore, identifying lncRNAs and their targets helps understand the biological processes and mechanisms of plants. We predicted 7,682 lncRNAs and target genes for 6474 lncRNAs. AS allows the plant to enhance gene coding potential, regulate gene expression via different mechanisms, and play essential functions in plant development[45,46]. We identified 4,091 AS events involved in 6,188 transcripts. Interestingly, 32 transcripts in AS events participate in the regulation of photoperiodism, flowering with a p-value of 1.08 × 10−8.

      The correlation analysis of transcriptomic and metabolomic data better explains transcription regulation mechanisms in metabolic pathways[47,48]. To characterize sex differentiation in the buffalograss transcriptome and its regulation, we measured and analyzed the transcriptional and metabolomic responses of B. dactyloides to sex differentiation. The most important pathways responding to sex differentiation were glutathione metabolism, photosynthesis, and plant hormone metabolism in B. dactyloides.

      The study demonstrated that glutathione metabolism was the most significantly enriched pathway in response to sex differentiation. One hundred and fifty four genes and two metabolites (glutamate and GSSG) related to glutathione metabolism had altered contents. GCL, GSS, GSR, DHAR, and GST participate in glutathione metabolism to generate glutamate, cysteine, and glycine, which are involved in several critical cellular functions, such as scavenging H2O2 and protecting cells from oxidative stress[47,49]. Based on our data, we can infer that glutathione metabolism was in response to sexual differentiation.

      Notably, most DEGs were upregulated in male and female plants relative to monoecious plants, including those encoding PSⅡ proteins (such as PsbB and PsbO), cytochrome b6/F complex (such as PetA and PetC), electron transport proteins (such as PetE and PetJ), PSⅠ proteins (such as PsbB and PsbO), and ATP production proteins (gamma and delta). These genes govern light harvesting, electron transport, reduction-oxidation reactions, and ATP production[41,50]. The increased ATP and Pi levels aligned with this result. The above results supported the positive functions of maintaining the photosynthetic activity of leaves in male and female plants.

      Three hormones (ABA, GA, and CTKs) varied in metabolite content among the sexes. NCEDs are a crucial enzyme in ABA biosynthesis[51]. We found the disruption of multiple NCEDs and the lowest abscisic alcohol and ABA levels in monoecious plants relative to female and male plants. Interestingly, we found the expression profile between NCEDs and ABA2/AAO3 in female, male, and monoecious plants was opposite, which suggests NCEDs, ABA2, and AAO3 were candidate genes involved in ABA homeostasis. ABA-GE, a stored form of ABA[51], consistently performs with abscisic alcohol and ABA in response to sexual differentiation. The lowest ABA content and disruption of multiple AOGs and CYP707A were found in monoecious plants compared to female and male plants, which suggests that reduced ABA leads to reduced ABA catabolism. In addition, we found that ABA signal transduction was more active in female and male plants than in monoecious plants. The above results showed that ABA biosynthesis, catabolism, and signaling synergize in response to sexual differentiation.

      The active GAs (mainly GA4) are essential for stamen development and control of pistil growth in Arabidopsis[52]. The GA biosynthesis essential genes, KAO and GA20ox, are more active in male plants than female plants (adjusted p-value < 0.01). Furthermore, GA4 signal transduction was most active in female plants, followed by male plants, and inhibited in male plants. The above results showed that GA biosynthesis is inconsistent with GA signaling in response to sexual differentiation.

      Compared with male and female plants, we found the overexpression of all three IPTs, the disruption of multiple CKXs, and the highest tZ and iP levels in monoecious plants. The results were in line with the previous report that tZ and iP levels are regulated and fine-tuned by IPTs and CKXs[53]. In addition, tZ deactivation was associated with the up-regulation of CTK modification enzymes, UGT73C and UGT85A1, and the results corresponded with the previous report[51]. The cytokinin receptor CRE1 and critical regulators (AHPs and B-ARRs) function in CTK perception and signaling[53]. Surprisingly, there was little difference in the number of DEGs encoding the proteins CRE1, AHPs, and B-ARRs in female, male, and monoecious plants, which is likely due to compensatory changes in the expression of CTK signaling genes to maintain the appropriate level of CTK function[53]. A-ARRs are considered CTK signaling negative regulators[51,53]. Compared with female and monoecious plants, we found disruption of multiple A-ARRs and the lowest tZ levels in male plants. The results suggest that reduced tZ leads to reduced CTK signaling[51]. The results above suggest that active CTKs and genes encoding IPTs, CKXs, UGT73C, UGT85A1, and A-ARRs involved in CTK metabolism and signaling play critical roles in response to sex differentiation.

    Conclusions

    • Few studies have investigated transcript annotation and complete mRNA structure in B. dactyloides, and the pathways of species-specific factors in sex differentiation remain unclear. We performed full-length transcriptome, second-generation transcriptome, and metabolome analysis to specify candidate factors influencing sex differentiation. We first provided the full-length transcriptome of B. dactyloides, identified 110,870 full-length transcripts, and obtained 90.52% of transcript annotation information. Then we identified 49,448 DEGs and 3,070 DAMs in female, male, and monoecious leaf samples. The co-enrichment analysis indicated that sexual differentiation was regulated by glutathione metabolism, photosynthesis, plant hormone biosynthesis, catabolism, and signaling. The identification of DEGs and DAMs that participate in glutathione metabolism, photosynthesis, ABA, CTK, and GA biosynthesis, catabolism, and signaling has helped illuminate the roles of plant hormones in the sex differentiation of B. dactyloides. The full-length transcriptomic data will allow further functional gene studies. Integrating transcriptomic and metabolomic data advances knowledge of the molecular mechanism of sex differentiation and provides information on breeding programs for B. dactyloides.

    Author contributions

    • The authors confirm contribution to the paper as follows: methodology, investigation, writing-original draft preparation, writing-reviewing and editing: Guan J; investigation, supervision, validation: Muye Liu; investigation: Yue Y; investigation, formal analysis, visualization: Yin S; investigation, validation: Teng W; investigation, validation: Zhang H, Wen H; resources, investigation, validation: Wu J; conceptualization, resources, writing-reviewing and editing, funding acquisition: Teng K; conceptualization, supervision, writing - review & editing: Fan X. All authors reviewed the results and approved the final version of the manuscript.

    Data availability

    • The PacBio SMRT reads and the Illumina NGS reads generated in this study were submitted to the BioProject database of the National Center for Biotechnology Information (accession numbers PRJNA991319, https://dataview.ncbi.nlm.nih.gov/object/PRJNA991319?reviewer=d9n2nbtdcpp1l5q6keuit0ih0b; PRJNA991464, https://dataview.ncbi.nlm.nih.gov/object/PRJNA991464?reviewer=7lj1iamddr3jk1nhtvf6fikp7i). The metabolomics data generated in this study were submitted to MetaboLights database (accession number MTBLS8147, https://www.ebi.ac.uk/metabolights/editor/MTBLS8147/descriptors).

      Acknowledgments

      • This work was funded by the Scientific Funds of Beijing Academy of Agriculture and Forestry Sciences (Grant Number: KJCX20210431) and Key R&D projects of Yunnan Province (Grant Number: 202203AP140153).

      Conflict of interest

      • The authors declare that they have no conflict of interest.

      Supplementary information
      • Supplemental Table S1 Primers used for qRT-PCR validation.
      • Supplemental Table S2 Prediction of Simple sequence repeats (SSRs) out of transcript datasets.
      • Supplemental Table S3 List of targets of lncRNA in B. dactyloides.
      • Supplemental Table S4 Summary of AS events in B. dactyloides.
      • Supplemental Table S5 Illumaina Sequencing Data Assessment.
      • Supplemental Table S6 Differentially expressed genes identified between different sample groups.
      • Supplemental Table S7 Differentially accumulated metabolites identified between different sample groups.
      • Supplemental Fig. S1 Summary of SMRT sequencing. (A) CCS read length distribution. (B) FLNC sequences read length distribution. (C) Consensus isoforms read length distribution.
      • Supplemental Fig. S2 The quality evaluation of transcriptomic and metabolomic data. PCA analysis of all transcripts (A) and metabolites (B), respectively. Heatmap analysis of correlation between samples in all transcripts (C) and metabolites (D), respectively.
      • Supplemental Fig. S3 KEGG enrichment of the up-regulated and down-regulated DEGs identified in Male vs Female (A), Monoecious vs Female (B), and Monoecious vs Male (C), respectively.
      • Supplemental Fig. S4 Prediction of regulatory proteins. A Number and family of top ten regulatory proteins predicted by SMRT sequencing data. B Venn diagram of differentially expressed regulatory proteins. C Number and family of top ten differentially expressed regulatory proteins. D-F The differentially expressed regulatory proteins predicted in Male vs Female, Monoecious vs Female, and Monoecious vs Male, respectively.
      • Supplemental Fig. S5
      • Supplemental Fig. S5 GO enrichment of sex-specific genes. (A) Genes present only in female plants. (B) Genes absent in female plants. (C) Genes present only in male plants. (D) Genes absent in male plants. (E) Genes present only in monoecious plants. (F) Genes absent in monoecious plants.
      • Supplemental Fig. S6 KEGG enrichment of sex-specific genes. (A) Genes present only in female plants. (B) Genes absent in female plants. (C) Genes present only in male plants. (D) Genes absent in male plants. (E) Genes present only in monoecious plants. (F) Genes absent in monoecious plants.
      • Supplemental Fig. S7 KEGG enrichment analysis of DAMs in the positive (left) and negative (right) ion mode in the comparison of Male vs Female (A), Monoecious vs Female (B), and Monoecious vs Male (C) sample groups.
      Rights and permissions
      • Copyright: © 2023 by the author(s). Published by Maximum Academic Press, Fayetteville, GA. This article is an open access article distributed under Creative Commons Attribution License (CC BY 4.0), visit https://creativecommons.org/licenses/by/4.0/.
    Figure (9)  References (53)
  • About this article
    Cite this article
    Guan J, Yue Y, Yin S, Teng W, Zhang H, et al. 2023. Integration analysis of full-length transcriptomics and metabolomics provides new insights into the mechanism of sex differentiation in buffalograss (Buchloe dactyloides). Grass Research 3:24 doi: 10.48130/GR-2023-0024
    Guan J, Yue Y, Yin S, Teng W, Zhang H, et al. 2023. Integration analysis of full-length transcriptomics and metabolomics provides new insights into the mechanism of sex differentiation in buffalograss (Buchloe dactyloides). Grass Research 3:24 doi: 10.48130/GR-2023-0024
    shu

Catalog

  • HTML

  • Introduction
  • Materials and methods

  • Plant materials and RNA sample preparation

  • Full-length library construction and PacBio SMRT sequencing

  • Full-length transcriptome structure analysis

  • Identification of regulatory proteins, lncRNAs and lncRNA target

  • Function annotation of transcripts

  • Illumina cDNA library construction and sequencing

  • Quantification of transcript expression and analysis of differential expression transcripts

  • qRT-PCR analysis

  • Metabolites extraction and LC-MS/MS analysis

  • Quantification of metabolites and analysis of differential metabolites

  • Data preprocessing and annotation, and correlation analysis between DEGs and DAMs

  • Results

  • SMRT sequencing analysis and function annotation of transcripts

  • Full-length transcriptome structure analysis

  • Acquisition and mining of transcriptome data

  • Regulatory proteins dynamics during B. dactyloides sex differentiation

  • Sex-specific genes analysis

  • Acquisition and mining of metabolomic data

  • Correlation analysis of transcriptomic and metabolomic data

  • Glutathione metabolism and photosynthetic capacity changes in response to sex differentiation

  • Plant hormones in response to sex differentiation

  • Discussion
  • Conclusions
  • Author contributions
  • Data availability
  • Acknowledgments
  • Conflict of interest
  • Supplementary information
  • Rights and permissions
  • About this article
logo

Contents

Policies

Services

Partnerships

© 2025 Maximum Academic Press, except Open Access articles. MAP is a member of ALPSP, COPE, CrossRef, cOAlition S, OASPA and STM.

Privacy Policy  |  Site Map  |  Terms and Conditions

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return