Diurnal metabolic regulation of isoflavones and soyasaponins in soybean roots (2025)

2.1 Chemicals

Malonyldaidzin and malonylgenistin were purchased from Nagara Science. Soyasaponin Ab, soyasapogenol A, and soyasapogenol B were purchased from Funakoshi. Soyasaponin Bb was purchased from ChromaDex. The other chemicals were purchased from Wako Pure Chemical Industries Ltd. or Nacalai Tesque Inc., unless otherwise stated.

2.2 Plant materials and growth conditions

The soybean seeds (cv. Enrei) used in the study were purchased from Tsurushin Shubyo (Matsumoto, Japan). The growth condition of the soybeans in hydroponic cultures was set up according to the description of Sugiyama etal.(2016). After 7days of the growth in autoclaved vermiculite containing water at 25°C with 16/8-hr photoperiods, the seedlings were rinsed and transferred to a hydroponic culture system where the soybeans were grown in 450-mL plastic containers filled with a mineral nutrient medium consisting of 3.0mM MgSO₄, 6.3mM KNO₃, 0.87mM KCl, 1.4mM KH₂PO₄, 2.5mM Ca(NO₃)₂, 21μM Fe–EDTA, 4.5μM KI, 28μM MnCl₂, 19μM H₃BO₃, 2.3μM ZnSO₄, 0.5μM CuSO₄, and 0.003μM Na₂MoO₄, with pH 6.0. The soybeans were kept in a cultivation room set at 25°C with a 12-hr light/dark cycle. After 2weeks, the plants were transferred to a new medium 6hr before sampling for analysis. The leaf and root tissues and root exudates were sampled at ZT0 (6:00 a.m.), ZT6 (0:00p.m.), ZT12 (6:00p.m.), and ZT18 (0:00 a.m.) for 48hr with five biological replicates.

2.3 RNA extraction and transcriptome sequencing

The total RNA was derived from soybean leaves and roots using RNeasy Plant Mini Kits (Qiagen, CA) according to the manufacturer's instruction. The DNA in each total RNA sample was digested using DNase I (RNase-free DNase sets, Qiagen). The RNA-seq library was prepared through the Lasy-Seq v1.1 protocol (Kamitani etal.,2019; https://sites.google.com/view/lasy-seq/) using 500ng total RNA. The library was sequenced by paired-end 150bp+150bp mode of HiSeqX platform (Illumina).

2.4 Transcriptome data analysis

The raw-reads data were quality controlled by removing low-quality bases using Trimmomatic (Bolger etal.,2014) with default parameters. The trimmed reads were aligned to the soybean genome (Glycine_max_v2.1 assembly) (Schmutz etal.,2010) using STAR v2.7.0f (Dobin etal.,2013) based on Ensembl Plants release 43 (Monaco etal.,2014) gene annotations. Gene expression levels were estimated as transcripts per million (TPM) (Wagner etal.,2012) using RSEM v1.3.1 (Li & Dewey, 2011) with default parameters. The transcriptome data set supporting the results of this study is publicly available at the DNA Data Bank of Japan (https://www.ddbj.nig.ac.jp) (DRA010744). Principal component analysis (PCA) was performed using the whole transcriptome data set (30,362 genes; TPM>1) after removing low-expression genes. Diurnally rhythmic genes were detected using the JTK_CYCLE algorithm (Hughes etal.,2010) with a MetaCycle package (Wu etal.,2016) in R environment, with 24-hr periodicity and false discovery rate (FDR) < 0.01. Rhythmic genes were classified into eight patterns based on the phase obtained by JTK_CYCLE. Gene Ontology (GO) enrichment analysis of gene set was performed using the SoyBase GO Term Enrichment Tool (http://www.soybase.org) in accordance with a diurnal pattern of interests. Gene IDs for isoflavone and soyasaponin biosynthesis, ABC and MATE transporters, and ICHG were collected based on the procedures described in previous literature (Ahmad etal.,2017; Krishnamurthy etal.,2019; Liu etal.,2016; Mishra etal.,2019; Sundaramoorthy etal.,2019; Yoo etal.,2013). The networks were constructed using the network visualization software Cytoscape (v. 3.7.2; Shannon et al., 2003).

2.5 Preparation of root extracts and exudates

The preparation of root extracts and exudates was performed according to the procedure previously described by Sugiyama etal.(2016). The medium containing root exudates was filtered through Omnipore membrane filters (Millipore). The medium was passed through a Sep-Pak C18 Plus short cartridge (Waters), which was eluted with 2ml of MeOH. The eluant was dried under nitrogen and reconstituted in 50µl of MeOH for LC-MS/MS analysis.

2.6 LC-MS/MS analysis

The samples were separated using an ACQUITY UPLC BEH C18 Column (2.1×50mm, 1.7µm, Waters) on an LC system (ACQUITY H-Class System, Wnters). The LC mobile phase consisted of (C) water containing 0.1% (v/v) formic acid and (D) acetonitrile. The gradient program was isocratic at 10% D, Initial; linear 1at 0%–85% D, 0–15min; isocratic at 100% D, 15–16min; and isocratic at 100% D, 16–20.5min. The injection volume of each sample was 5µl, and the flow rate was 0.2mlmin. The isolated samples were detected using a tandem quadrupole MS (Xevo TQ-S, Waters) in the Multiple Reaction Monitoring (MRM) mode. The MRM conditions for the respective compounds are listed in TableS1.

2.7 Statistical analysis

Statistical differences were calculated using the Tukey's HSD test at p< 0.05 implemented in R (v. 3.6.1; R Core Team, 2019). The outliers (defined as >1.5*IQR) in the boxplots were excluded from the calculated data.

3.1 Time-dependent transcriptome analysis

The RNA-seq analysis yielded over 1.5 billion reads from 80 samples, that is, eight time points, two tissues, and five replicates. We conducted PCA and found that the leaf and root transcriptomes were clearly separated (FigureS1). The circadian clock-related transcription factors in soybeans include Late elongated hypocotyl/circadian clock associated (LCL), TOC1, PRR, GIGANTEA (GI), LUX ARRHYTHMO (LUX), and Early flowering 4 (ELF4) (Cheng etal.,2019; Li etal.,2013; Liew etal.,2017; Wang etal.,2020). The expression of these key circadian clock-related genes showed diurnal changes both in the leaves and roots and displayed a similar expression tendency as observed in previous studies (Cheng etal.,2019; Li etal.,2013; Liew etal.,2017; Locke etal.,2018; Marcolino-Gomes etal.,2014; Wang etal.,2020) (Figure2).

To provide an overview of the diurnal changes of gene expressions in soybean, rhythmic genes with 24-hr periodicity were examined by the JTK_CYCLE algorithm (FDR<0.01) (Hughes etal.,2010). Rhythmic genes were sorted into eight phases according to the peak time of the gene expression. A total of 6,777 and 1,240 genes were detected as diurnally oscillating genes among soybean genome (55,897 genes) in leaves and roots, respectively (Figure3, FigureS2, TableS2 and Table S3). To understand the characteristic biological processes according to time of day, we carried out GO enrichment analysis for the diurnally rhythmic genes (TableS2 and Table S3). In the leaves, the GO terms related to light response were represented at ZT0 and the GO terms associated with photosynthesis were frequently represented from ZT3 to ZT9, which consistent with the GO enrichment analysis on diurnal transcripts in rice leaves (Xu etal.,2011) (Figure4a, TableS4). During nighttime, the GO terms involved in protein synthesis, respiration, and transport were enriched in chronological order (Figure4a, TableS4). In roots, GO terms related to transport were highlighted at ZT0 and the GO terms associated with circadian rhythm were enriched at ZT9 (Figure4b, TableS5). The GO terms linked to secondary metabolism were not enriched at any time of day, although Xu etal.(2011) represented that secondary metabolic process was induced from ZT4 (10:00a.m.) to ZT6.

3.2 Diurnal variation of biosynthetic gene expression for isoflavone and soyasaponin

To analyze the diurnal metabolic changes of isoflavones and soyasaponins in soybean, we investigated the fluctuations of gene expression for all the reported genes involved in isoflavone and soyasaponin biosynthesis. No diurnal variation pattern common to the isoflavone or soyasaponin biosynthetic genes was observed in the leaves (Figures5a and 6a). In contrast, the expression profiles in the roots exhibited clear diurnal variation patterns for the genes involved in isoflavone and soyasaponin biosynthesis (Figures5b and 6b, TableS6 and Table S7). Isoflavone biosynthetic genes, such as chalcone reductase (CHR), chalcone synthase (CHS), isoflavone synthase (IFS), and 2-hydroxyisoflavanone dehydratase (HIDH) showed high expression at ZT6 and low expression at night (Figure5b, TableS6). This variation pattern is inconsistent with that of the flavonoid biosynthetic genes in Arabidopsis (Harmer etal.,2000), which were highly expressed at night under constant light conditions after 7days of culture in a 12-hr light/dark condition. The expression levels of genes encoding β-Amyrin synthase (BAS), cytochrome P450 (CYP) for triterpenes, and UDP-glucuronosyltransferase (UGT) involved in soyasaponin biosynthesis (Krishnamurthy etal.,2019; Sundaramoorthy etal.,2019) increased from ZT18 (0:00a.m.) to ZT0 (6:00a.m.) and decreased at ZT6, which was an inverse pattern to that of isoflavone biosynthesis (Figure6b, TableS7).

The transcription factors of the MYB family play crucial roles in the regulation of isoflavone biosynthesis. GmMYB29, GmMYB102, GmMYB133, GmMYB176, GmMYB280, and GmMYB502 are positive regulators of isoflavone biosynthetic genes (Anguraj Vadivel etal.,2019; Bian etal.,2018; Chu etal.,2017; Sarkar etal.,2019; Yi etal.,2010), while GmMYB29B, GmMYB39, and GmMYB100 negatively regulate those genes (Jahan etal.,2020; Liu etal.,2013; Yan etal.,2015). Among these MYBs, GmMYB176 showed the highest expression level in the roots (Figure7, TableS8). GmMYB176 exhibited a remarkable diurnal variation, with a higher expression level at ZT0 and ZT6 than at ZT12 (Figure7, TableS8), indicating the roles of this MYB member in the induction of isoflavone biosynthetic genes during the daytime.

3.3 Co-expression network of isoflavone and soyasaponin metabolism-related genes

The transporters engaged in the accumulation or secretion of both isoflavones and soyasaponins in soybeans have not been reported to date. In this study, co-expression network analysis was performed using Spearman's correlation coefficient threshold value of 0.8 to explore the candidate genes responsible for isoflavone and soyasaponin transport. Isoflavone biosynthetic genes displayed a separated co-expression network from those of soyasaponin biosynthesis (Figure8). The isoflavone biosynthetic gene cluster included three genes coding for ABC transporters (Glyma.10G019000, Glyma.13G361900, and Glyma.15G011900), and a MATE-type transporter gene Glyma.01G026200 (Figure8a, Table1). Moreover, the cluster of soyasaponin biosynthetic genes only contained an ABC transporter gene, Glyma.15G148500 (Figure8b, Table1). When co-expression network analysis was performed using a Spearman's correlation coefficient threshold value of 0.7, the isoflavone biosynthetic gene cluster contained five additional ABC transporter genes, namely Glyma.01G008200, Glyma.03G101000, Glyma.07G233900, Glyma.13G043800, and Glyma.20G242000, whereas the soyasaponin biosynthetic gene cluster included four more ABC transporter genes, namely Glyma.04G069800, Glyma.08G101500, Glyma.19G021500, and Glyma.19G184300 (FigureS2, TableS9). The correlation analysis using publicly open transcriptome data has recently become available in SoyCSN for soybean (Wang etal.,2019). The co-expression of the ABC transporter genes (Glyma.10G019000, Glyma.13G043800, Glyma.13G361900, and Glyma.15G011900) with CHS7 (Glyma.01G228700) and IFS1 (Glyma.07G202300), and the co-expression of two other ABC transporter genes (Glyma.15G148500 and Glyma.19G021500) with BAS1 (Glyma.07G001300) and CYP93E1 (Glyma.08G350800) were also obtained from SoyCSN (Wang etal.,2019). These findings suggest that the ABC transporter genes are co-expressed with the biosynthetic genes of isoflavone or soyasaponin not only in diurnal variations but also in the symbiosis or other tissues. It is thus expected that those transporters mediate in the vacuolar accumulation and secretion into the rhizosphere of isoflavones or soyasaponins.

The coordinated gene expression for biosynthesis and transport has been reported in several studies. For instance, MtMATE1 of Medicago truncatula was reported as an epicatechin 3'-O-glucoside transporter identified from the co-expression with its biosynthetic genes (Zhao & Dixon,2009). A glucosinolate transporter 1 (GTR1) in Arabidopsis was identified as jasmonoyl-isoleucine and gibberellin transporter based on its co-expression with jasmonate biosynthetic genes (Saito etal.,2015). Furthermore, Catharanthus roseus CrNPF2.9 was identified by co-expression analysis of the biosynthetic genes of monoterpene indole alkaloids (Payne etal.,2017), and saffron (Crocus sativus) CsABCC4a was revealed as a transporter for crocins by narrowing down the candidate transporter genes that were highly expressed in pistils (Demurtas etal.,2019). These reports prompted us to genetically and biochemically investigate candidate transporter genes for isoflavone and soyasaponin transport further.

ICHG is grouped within the soyasaponin biosynthetic gene cluster, although it hydrolyzes isoflavone glucosides and not soyasaponins, presumably (Suzuki etal.,2006). Nevertheless, the expression level of ICHG in the roots from ZT18 to ZT0 was high, whereas it was low at ZT6 (Figure9), which was consistent with the pattern of the soyasaponin biosynthetic genes (Table1, Figures6 and 8). This could be a coincidence because ICHG is the most downstream gene in the secretion of isoflavones and could be expressed later than those for isoflavone biosynthesis, while there remains a possibility that ICHG hydrolyzes the glucoside linkage of soyasaponins in the apoplast.

3.4 Diurnal variation of isoflavones and soyasaponins in roots and root exudates

The gene expression in isoflavone and soyasaponin biosynthesis displayed diurnal variations in the roots. We analyzed the contents of these specialized metabolites in both the roots and root exudates using LC-MS/MS. The daidzein content was lowest at ZT0 and increased from ZT6 to ZT18 (Figure10a), followed by the highest expression of isoflavone biosynthetic genes at ZT6 (Figure5b). This temporal difference between gene expression and metabolite accumulation was consistent with the flavonoid biosynthesis in Arabidopsis (Nakabayashi etal.,2017). The contents of its glycosides, daidzin, and malonyldaidzin, differed among time points, but did not have a variation pattern. The contents of genistein and its glucosides were about 5- to 10-fold lower than those of daidzein and its glucosides (Figure10a). Genistein and genistin contents exhibited similar trends with that of daidzein, while malonylgenistin did not show apparent diurnal variations (Figure10a). In contrast to the accumulation of isoflavone glucosides in the roots, the aglycone daidzein was the predominant isoflavone form found in root exudates as observed previously (Sugiyama etal.,2016). The amount of daidzein and genistein in root exudates did not show clear diurnal variations (Figure10b). In contrast, the contents of glucosides, such as daidzin, malonyldaidzin, genistin, and malonylgenistin, showed a tendency to increase during ZT6 to ZT12 (Figure10b), which is in concordance with the decreased ICHG expression at ZT6 (Figure9). These findings implicate the possible involvement of two modes in the secretion of isoflavones to the rhizosphere, that is, one for the secretion of aglycones and another for the secretion of glucosides, because daidzein content constantly remained at its highest among isoflavones even at ZT6 when ICHG expression was suppressed.

We also determined the contents of soyasaponins using LC-MS/MS. We focused on two major soyasaponins in root exudates of soybean, soyasaponin Ab and soyasaponin Bb, and their aglycones, soyasapogenol A and soyasapogenol B (Tsuno etal.,2018), due to the limited availability of authentic samples for soyasaponins. The root content of soyasaponin Bb was low at ZT0 but was significantly increased at ZT6 and maintained until ZT18 (Figure11a). The root content of soyasaponin Ab did not exhibit a clear diurnal variation pattern. In root exudates, contents of soyasaponin Ab and Bb slightly varied but showed no apparent diurnal pattern (Figure11b). The amount of soyasaponin aglycones in the roots was much lower than that of their glucosides. Soyasapogenol A was not detected. Soyasapogenol B was present in small quantities, and it was high at ZT6 and low at night (Figure11a). The content of soyasapogenol B in the root exudates was relatively high from ZT0 to ZT6, which would reflect the induction of soyasaponin biosynthetic genes during the night and the slight increase of soyasapogenol B in roots at ZT6 (Figure11b).