ORIGINAL RESEARCH article
Sec. Systems Biology Archive
Volume 11 - 2020 | https://doi.org/10.3389/fgene.2020.589326
Temporospatial Flavonoids Metabolism Variation in Ginkgo biloba Leaves
- 1Co-Innovation Centre for Sustainable Forestry in Southern China, Nanjing Forestry University, Nanjing, China
- 2College of Forestry, Nanjing Forestry University, Nanjing, China
- 3Department of Forest & Conservation Sciences, Faculty of Forestry, The University of British Columbia, Vancouver, BC, Canada
Ginkgo (Ginkgo biloba L.) is a high-value medicinal tree species characterized by its flavonoids beneficial effects that are abundant in leaves. We performed a temporospatial comprehensive transcriptome and metabolome dynamics analyses of clonally propagated Ginkgo plants at four developmental stages (time: May to August) across three different environments (space) to unravel leaves flavonoids biosynthesis variation. Principal component analysis revealed clear gene expression separation across samples from different environments and leaf-developmental stages. We found that flavonoid-related metabolism was more active in the early stage of leaf development, and the content of total flavonoid glycosides and the expression of some genes in flavonoid biosynthesis pathway peaked in May. We also constructed a co-expression regulation network and identified eight GbMYBs and combining with other TF genes (3 GbERFs, 1 GbbHLH, and 1 GbTrihelix) positively regulated the expression of multiple structural genes in the flavonoid biosynthesis pathway. We found that part of these GbTFs (Gb_11316, Gb_32143, and Gb_00128) expressions was negatively correlated with mean minimum temperature and mean relative humidity, while positively correlated with sunshine duration. This study increased our understanding of the molecular mechanisms of flavonoids biosynthesis in Ginkgo leaves and provided insight into the proper production and management of Ginkgo commercial plantations.
Ginkgo (Ginkgo biloba L.) leaves contain a variety of medicinal compounds, which have been used in healthcare and food industries. Flavonoids are the major bioactive ingredients in Ginkgo leaves, including flavonols, flavones, and anthocyanins (Meng et al., 2019). These molecules have been reported to have beneficial effects in preventing metabolic syndrome at different levels such as early stage Alzheimer’s and cardiovascular diseases (Tian et al., 2017; Gruenwald et al., 2020). Flavonoids also act as growth regulators controlling single organ and whole plant development (Agati et al., 2012). Therefore, it is essential to understand the molecular mechanisms of flavonoids accumulation during Ginkgo leaves development to ultimately improve the production and management of Ginkgo plantations.
Recently, considerable efforts have been dedicated to improving Ginkgo leaves flavonoids for commercial production. Studies showed that several agronomic measures could increase the flavonoids content, such as alternative partial root-zone irrigation (Wang et al., 2016), fertilization (Guo et al., 2016), and foliar fertilization (Wu et al., 2020). Treatments with salicylic acid, UV-B, and NaCl, all have shown a positive effect on increasing Ginkgo leaves flavonoids content (Ni et al., 2017, 2018; Zhao et al., 2020). More importantly, additional efforts have been directed at the molecular level to achieve the same objective. For example, transcriptome libraries have been constructed for various Ginkgo tissues (Ye et al., 2019) and leaves with different flavonoid contents (Wu et al., 2018), for improving the understanding of flavonoid biosynthesis. Another strategy for improving Ginkgo leaves metabolites yield is through genetic engineering; however, detailed information on gene expression profiling and transcriptional dynamics that regulate flavonoids accumulation is scarce.
Leaves undergo a series of developmental and physiological changes during their lifespans, involving complex, but highly regulated molecular processes to maximize fitness in a given ecological setting (Leopold, 1961; Fenner, 1998). It was found that leaves from the same Ginkgo tree could exhibit differences in flavonoids content at different developmental stages (young vs. mature leaves) (Guo et al., 2020). Additionally, Ginkgo leaves from plants growing at different elevations (different environments), and the same growing period also displayed substantial differences in their flavonoids accumulation (Zou et al., 2019). However, the current understanding of flavonoids accumulation regulation mechanism, which varies according to the development stage and geographical distribution, is limited. Transcriptomes can provide information regarding gene expression and regulation at specific developmental stages or under specific physiological conditions (Sato et al., 2011). Furthermore, integration of different-omics data, such as metabolome, will help elucidate the complex mechanism controlling flavonoid biosynthesis (Weckwerth, 2008).
Here, we conducted comprehensive temporospatial transcriptome and metabolome dynamics analyses of clonally propagated Ginkgo plants at four developmental stages (May to August) across three different environments (test-sites) to unravel leaves flavonoids biosynthesis spatial-temporal variation. The study-specific objectives are to: (1) quantify the transcriptional responses to spatial (environmental cues) and temporal (development stages) conditions; (2) explore the association between flavonoids accumulation and expression of flavonoid related structural genes; and (3) elucidate the regulatory network involved in gene expression associated with flavonoids biosynthesis. The broader aim of this work is intended to improve our understanding of the transcriptional dynamics that regulate flavonoids accumulation at the molecular level and provide insightful information for enhancing flavonoids content of Ginkgo leaves for the proper production and management of Ginkgo commercial plantations.
Materials and Methods
Plant Materials and Sample Collection
Generally, the optimum age of flavonoids production in Ginkgo leaf-harvest plantations is trees under 5-year-old (Zou et al., 2019), thus older trees are considered suboptimum. Therefore, in the present study we utilized leaf samples collected from 2-year-old clonally propagated (grafted) Ginkgo trees. Trees are spatially replicated over three test sites (i.e., different environments) (Table 1). These sites are: (1) Yi Ning (YN), located in northwestern China (lat.: 43.41°N, long.: 81.11°E), characterized by a typical mid-temperate continental semi-arid climate; (2) Pi Zhou (PZ), located in central China (lat.: 34.21°N, long.: 117.58°E) characterized by a warm temperate monsoon climate; and (3) Qu Jing (QJ), located in southern China (lat.: 25.52°N, long.: 103.58°E), characterized by a subtropical plateau monsoon climate. In each site, the experiment is planted as a complete randomized block design with three blocks (replicates), each harboring 20 Ginkgo clones.
Table 1. Geographical distribution and climate factors [mean annual temperature (MAT), mean annual precipitation (MAP), and mean annual sunshine duration (MASD)] of the studied three test sites [Yi Ning (YN), Pi Zhou (PZ), and Qu Jing (QJ)].
Samples were conducted between leaf expansion (May, after majority of leaves expansion) and leaf “commercial” ripening (August, before autumnal senescence). During this biological window, Ginkgo leaves are at their substance’s peak activity and are easy to harvest and store (Ellnain-Wojtaszek et al., 2002). Leaves were collected on a clear day in the middle of each month (May to August) to represent four temporal leaf developmental stages. A single clone was randomly selected across the three blocks (i.e., 3 biological replications) and the collected leaf samples provided the material for the metabolomics and transcriptomics analyses. Each sampled tree was represented by three crown positions (top, middle, and bottom), each provided a single complete and healthy leaf. In total, 36 samples (4 development stages × 3 environments × 3 biological replicates) were used for the metabolomics and transcriptomics analyses. Collected leaves were immediately preserved in liquid nitrogen, freeze-dried, and kept at −80°C until further use. To measure the temporal variation in total flavonoid glycosides (TFG) content, a monthly time series sampling was conducted (at mid-month between May and August) on the PZ site and nine leaves were randomly collected from the 20 Ginkgo clones planted in each block. Additionally, to measure the spatial changes of TFG content across environments, nine leaves were randomly collected from the same 20 Ginkgo clones from the three blocks and the same sampling scheme was conducted across the three sites (sampling was conducted in mid-August). Leaves were oven-dried (70°C, 48 h), crushed, sieved through a 100-mesh sieve, and vacuum packed. All experiments were performed with three biological replicates.
Total Flavonoid Glycosides Measurement
Ginkgo leaves flavonoids were extracted following the Pharmacopoeia of the People’s Republic of China (PPRC) procedures (Commision, 2010), and flavonoid glycosides content were determined by high-performance liquid chromatography (HPLC). In brief, approximately 0.5 g of oven-dried leaf powder per sample was immersed in petroleum ether and refluxed at 70°C for 2 h to remove impurities. Samples were then soaked in methanol and each sample’s extract was evaporated on a rotary evaporator after refluxed at 70°C for 4 h. Subsequently, the pellet was washed with 25 mL of a 25% methanol-HCl (4:1, v/v) mixture and the eluent was collected and refluxed for 30 min. After cooling to room temperature, the eluent was brought to 50 mL with methanol, then used for determination by HPLC. HPLC (Waters 1525, United States) conditions were set as follows: the mobile phase was methanol and 0.4% H3PO4 solution (56:44, v/v) at 1.0 mL min–1; the column temperature was 30°C; the detection was performed at 360 nm. Quercetin, kaempferol, and isorhamnetin were selected as standard substances following the supplier’s specifications (Yuanye Biological Co., Shanghai, China).
Total flavonoid glycosides content = (quercetin + kaempferol + isorhamnetin) × 2.51 (Commision, 2010). Means and standard errors for each sample were calculated. Differences among samples were determined using one-way ANOVA and significant differences were detected (defined as P < 0.05) using the least significant difference (LSD) test.
The supernatant extraction for each sample was performed as previously described (Guo et al., 2020). In summary, about 50 mg freeze-dried sample was put into an EP tube after grinding. After the addition of 1 mL of extract solvent (acetonitrile-methanol-water, 2:2:1, containing 0.1 mg L–1 lidocaine as an internal standard), the samples were swirled for 30 s, homogenized at 45 Hz for 4 min, and sonicated for 5 min in an ice water bath. The homogenate and sonicate circle was repeated three times, followed by incubation at −20°C for 1 h and centrifugation at 12,000 rpm and 4°C for 15 min. The resulting supernatants were transferred to LC-MS vials and stored at −80°C for later use. LC-MS/MS analyses were performed using an UHPLC system (1290, Agilent Technologies) with a UPLC HSS T3 column coupled to Q Exactive (Orbitrap MS, Thermo). The mobile phase A was 0.1% formic acid in water for positive, and 5 mmol/L ammonium acetate in water for negative, and the mobile phase B was acetonitrile. The elution gradient was set as follows: 0 min, 1% B; 1 min, 1% B; 8 min, 99% B; 10 min, 99% B; 10.1 min, 1% B; 12 min, 1% B [see Supplementary Figure S1 for a UHPLC chromatogram of standards and samples from Yi Ning (YN) site at different sampling stages]. MS raw data files were converted to the mzML format using ProteoWizard, and processed by R package XCMS. OSI-SMMS (version 1.0, Dalian Chem Data Solution Information Technology Co. Ltd.) was used for peak annotation after data processing with an in-house MS/MS database. The metabolites were mapped to the Kyoto Encyclopedia of Genes and Genomics (KEGG) metabolic pathways to identify the substances in the related pathways of flavonoid biosynthesis (ko 00941- ko 00944).
Total RNA extraction, library preparation, and sequencing for each sample (36 libraries: four developmental stages in three different environments with three biological replicas) were performed as previously described (Guo et al., 2020). Total RNA was extracted from the freeze-dried samples using Trizol reagent kit (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s protocol. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and checked using RNase free agarose gel electrophoresis. Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA were synthesized by DNA polymerase I, RNase H, dNTP and buffer. Then the cDNA fragments were purified with QiaQuick PCR extraction kit (Qiagen, Venlo, Netherlands), end repaired, poly (A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeq2500. The Ginkgo Illumina raw sequencing data were submitted to the NCBI BioProject database under project number PRJNA657336.
An index of the reference genome was built, and paired-end clean reads were mapped to the Ginkgo’s reference genome1 using Hisat2. The mapped outputs were processed via StringTie software to obtain FPKM (fragment per kilobase of transcript per million mapped reads) for all the Ginkgo genes in each sample. Based on gene expression, principal component analysis (PCA) and hierarchical clustering analysis were performed with R packages, gmodels and pheatmap2, which were also used to reveal the relationship among samples. The FPKM data were directly used to estimate the differential expression of genes (DEGs) between samples. FDR < 0.05 and | log2FC| > 1 were used as thresholds to identify significant DEGs. The Short Time-series Expression Miner (STEM) software was used to obtain the temporal expression profile of DEGs. Subsequently, DEGs in enriched clustered profiles were used for KEGG pathway enrichment analysis (Q value ≤ 0.05) to assess metabolic pathways and related gene functions.
Weighted gene co-expression network analysis (WGCNA) was performed in the R environment. After filtering with the R package DCGL, a total of 23,182 genes (FPKM > 0) were reserved for subsequent analysis. The adjacency matrix between different genes was constructed with a threshold power of 10. A dynamic tree cut procedure (merge cut height = 0.70, min module size = 50) in R package WGCNA was used to identify similar modules in the hierarchical tree. The expression profile of module genes in each sample was displayed by module eigengene, which was defined as the first principal component of a given module. The Pearson correlations between the eigengenes of each module and the abundance of flavonoids were plotted by R package ggplot2. Subsequently, we identified the encoding transcription factor (TF) genes and the structural genes in the biosynthesis pathways of related flavonoids from the target modules. The gene regulatory networks were generated by Cytoscape software (Version 3.7.1).
The promoter region of TF genes was analyzed for presence of cis-acting regulatory elements by PlantCARE3 and visualized by TBtools software. Additionally, to explore the regulatory effect of environmental factors on the expression of TF genes, Pearson’s product-moment correlation analysis was conducted between TF genes expression and environmental factors during development (daily meteorological data for each area from May to August 20194).
Quantitative Real-Time PCR (qRT-PCR) Analysis
Ten genes involved in flavonoid biosynthesis were randomly chosen for validation by qRT-PCR. According to the manufacturer’s instructions, cDNA was obtained using MonScript RTIII All-in-One Mix with dsDNase kits (Monad, China) and qRT-PCR analysis was carried out using an Applied BiosystemsTM 7500 Real-Time PCR Systems (Monad, China). Primers used were designed in Primer Premier 5 (United States), and the primer sequences are provided in Supplementary Table S2. Glyceraldehyde 3- phosphate dehydrogenase (GAPDH, GenBank Accession No. L26924) gene was used as an internal standard. The relative transcript abundance was calculated using the 2–ΔΔC_T method (Livak and Schmittgen, 2001). August samples from YN, PZ, and QJ sites were used for qRT-PCR analysis. As designed, each sample included three biological replicates and three independent technical repetitions.
All statistical analyses were conducted in R environment (R Core Team, 2019). Differences among TFGs content were determined using one-way analysis of variance (ANOVA) and significant differences were calculated using the least significant difference (LSD) test (defined as P < 0.05). The complex relationships between gene expression profiles were intuitively displayed by a PCA plot and a cluster heat map. Relationships between expression of structural genes and abundance of flavonoids were evaluated using the Pearson’s product-moment correlation analysis (P < 0.05, significant correlation).
Changes in Total Flavonoid Glycosides Content
Temporarily, TFG content showed a declining trend with sampling time during the leaf development process. PZ site’s TFG content time-course analysis showed the highest value occurred in May, followed by a drastic drop in June (P < 0.05) and a significant recovery in July (Figure 1A). Compared to May samples, the TFG content of June, July, and August samples decreased by 66.40, 15.98, and 21.50%, respectively (Figure 1A). Spatially, apparent differences in TFG content were observed across the three growing environments. Compared to August’s PZ and QJ samples, the TFG content of YN was larger by 52.67 and 140.45%, respectively (Figure 1B).
Figure 1. Differences in total flavonoid glycosides (TFG) content among Ginkgo leaf samples. (A) Temporal differences in TFG content among developmental stages in PZ site. (B) Spatial differences in TFG content in August-samples across the three environments. Measurements were performed in triplicates. Error bars indicate standard deviations, and different capital and small letters represent a significant difference (P < 0.05) between developmental stages and environments, respectively.
Changes in Gene Expression Profile
Through transcriptional dynamics analysis, we identified approximately 2.0 billion clean reads from the 36 cDNA libraries that were mapped to the Ginkgo genome. The mapping rates of each library ranged from 91.14 to 95.39% (Supplementary Table S1). Spatially, PCA analysis results showed clear separation on the PC biplot, accounting for 61.1% of total gene expression variance in the data set (Figure 2A). Samples were spatially separated along the PC1 axis with YN, PZ, and QJ positioned on the left, middle, and right, respectively (Figure 2A). Temporally, within each site, the four developmental stage samples tended to follow the same left-to-right trend along PC1, while this trend did not persist for PC2 (Figure 2A). In the hierarchical clustering analysis, we did not detect any evidence of clustering among samples at either the different developmental stages (temporal) or at any given environment (spatial) (Figure 2B). Interestingly, PZ samples of the earlier stage (May) exhibited a closer correlation with the YN samples, whereas the PZ samples of a later stage (August) tended to correlate with the QJ samples, suggesting that major transcriptional program differences existed among development stages within each environment.
Figure 2. Temporospatial differences in gene expression profile. (A) Principal component analysis (PCA) plot showing clustering of leaf transcriptomes at four developmental stages (temporal) under three different environments (spatial). Red squares, green triangles, and blue circles represent samples from YN, PZ, and QJ sites, respectively. The change in color from light to dark represents the four leaf development stages (May to August). (B) Pearson product-moment correlation coefficients and clusters of the RNA-seq data from leaf samples. The redder the rectangle, the stronger the correlation between the samples, whereas the bluer the rectangle, the weaker the correlation.
Differential Gene Expression During Leaf Development
At each leaf developmental stage, we identified different expression genes (DEGs) among samples from different environments (Figure 3A). We found more DEGs differences existed between QJ and YN samples (number of stage-specific genes varied from 644 to 3,318), while fewer DEGs differences between PZ and QJ samples (number of stage-specific genes varied from 74 to 1,097). The variable number of DEGs differences suggested that each stage of Ginkgo clones development had an independent strategy in response to their respective different environmental conditions.
Figure 3. Temporospatial expression pattern of differential expressed genes (DEGs). (A) Venn diagrams and column charts showing DEGs between samples from different growth environments (spatial) at four development stages (temporal). The gray circle/rectangle represents the difference between samples from YN and QJ; the blue one represents the difference between samples from PZ and QJ; the orange one represents the differences between the PZ and YN samples. (B) Profile blocks with a colored background are significant clusters of the P ≤ 0.05, and the same color represents that the profiles are the same cluster. (C) An enriched KEGG map shows significant pathways among the genes of eight profiles. The red rectangle represents a significant enrichment pathway (P < 0.05).
To analyze the temporal expression pattern of DEGs, the 24,958 DEGs were further clustered by Short Time-series Expression Miner (STEM) software. There were 8 identifiable statistically significance (P < 0.05) temporal expression patterns, which were divided into 5 clusters containing a total of 1,908 DEGs (Figure 3B). The expression of DEGs contained in the 0 profile was gradually down-regulated during leaf development, while the temporal expression pattern of the 19 profile showed an opposite pattern. The KEGG pathway enrichment analysis of 1,908 DEGs revealed that 33 pathways were significantly enriched, including a large number of secondary metabolites, carbohydrate, and lipid metabolic pathways (Figure 3C). In particular, the phenylalanine and flavonoid biosynthesis pathways (ko 00940 and ko 00941) were enriched in several profiles. Therefore, these results suggested that the expression of some genes in the flavonoids-related biosynthesis pathways varied as a function of environmental factors (spatial) and developmental stages (temporal). The reliability of the RNA-seq results and the differentially expression analysis was further verified by qRT-PCR (Supplementary Figure S2).
Identification and Screening of Gene Co-expression Modules
Twelve modules were identified in a dendrogram comprising 105 – 3,908 genes, and each module harbored genes encoding the number of transcription factors (TFs) varying from 7 to 210 (Figures 4A,B). In most modules, TF-encoding genes accounted for more than 5% of the total genes, indicating that the transcriptional activity was strictly regulated. Also, the eigengene of each module was associated with the abundance of 17 flavonoids revealed by Pearson product-moment correlation coefficient analysis. Remarkably, three modules (Black, Blue, and Brown) exhibited a strong correlation (r > | 0.7|, P < 0.05) between gene expression and flavonoids accumulation (Figure 4C).
Figure 4. Identification and screening of gene co-expression modules. (A) Dendrogram showing modules identified by the weighted gene co-expression network analysis (WGCNA) and clustering dendrogram of expressed genes; (B) The number of genes and transcription factors (TFs) contained in each module; (C) Correlation coefficient between the abundance of flavonoids and module eigengenes presented with a color scale with red and blue representing positive and negative correlations, respectively.
To better understand the function of genes in these three modules, we assigned the genes to KEGG terms. The top 20 enriched pathways in each module were revealed by bubble maps. The genes from Black and Blue modules were significantly enriched in pathways related to translation, folding, sorting and degradation, signal translation, amino acid metabolism, and energy metabolism (Supplementary Figures S3A,B). In these pathways, some unigenes encoding glutathione S-transferase (GST), vacuolar sorting receptors (VSR), multi-antimicrobial extrusion protein (MATE) were found, which were thought to be involved in the transportation of flavonoids from cytosolic biosynthesis to their vacuolar accumulation (Petrussa et al., 2013). Notably, genes from the Brown module were significantly enriched in pathways related to secondary metabolites biosyntheses, such as flavonoid biosynthesis and phenylpropanoid biosynthesis (Supplementary Figure S3C). Interestingly, a large number of genes encoding key enzyme (flavonoids-related structural genes) had been identified in these pathways, including genes encoding phenylalanine ammonia-lyase (PAL), cinnamate 4-hydroxylase (4CH), chalcone synthase (CHS), chalcone isomerase (CHI), flavonol synthase (FLS), dihydroflavonol 4-reductase (DFR), anthocyanin synthase (ANS), anthocyanidin reductase (ANR), and UDP-glycosyltransferase (UGT).
Construction of Flavonoid-Related Gene Regulation Network
After screening the target modules, we constructed the biosynthesis pathways of six flavonoids (flavone, isoflavone, flavanonol, anthocyanins, flavonol, and flavonol glycoside) and identified the structural genes involved in these pathways from the Brown module (Figure 5A). The developmental stage specificity of the 12 flavonoids accumulation and the 15 structural genes expression was visualized (Figure 5B). We found that three flavones (luteolin, apigenin, and vitexin), one flavanonol (taxifolin), and two flavonol glycosides (astragalin and rutin) had the highest accumulation in May. In contrast, two isoflavones (daidzin and genistein), two anthocyanins (cyanidin and delphinidin), and one flavonol glycosides (isoquercitrin) had the lowest accumulation in leaves at the same developmental stage. Additionally, the accumulation of one flavonol (quercetin) was the highest in July (see Supplementary Table S3 for quantitative values of the 12 identified flavonoids). The 15 structural genes identified in the Brown module had similar developmental expression patterns, with high expression in May and low expression in August. Correlation analysis of transcriptome and metabolome indicated that some structural genes were significantly correlated with specific flavonoids (r > | 0.6|, P < 0.05). For example, quercetin content was positively correlated with the expression of a gene encoding FLS enzyme, while cyanidin content was negatively correlated with the expression of some genes encoding ANS, DFR, and UGT enzyme. Thus, these structural genes may play crucial roles in the accumulation of some flavonoids.
Figure 5. Temporospatial expression patterns of flavonoid-related structural genes in different leaf developmental stages (temporal) and association analysis with the accumulation of several kinds of flavonoids. (A) The biosynthesis pathways of several kinds of flavonoids. (B) The accumulation patterns of 12 flavonoids and the expression patterns of 15 flavonoid-related structural genes in different growth stages and the association analysis between them. The change in color of circle/rectangle from red to blue represents a gradual decrease in the abundance of flavonoids/expression of structural genes. The red line represents the positive correlation between the expression of structural genes and the abundance of flavonoids, and the blue line represents the negative correlation between them (r > | 0.6|, P < 0.05).
In Brown module, a total of 13 genes belonging to four transcription factor families were identified, including those encoding MYB (8 genes), ERF (3 genes), bHLH (1 gene), and Trihelix (1 gene), which may be involved in the regulation of flavonoids accumulation. To explore the regulatory effect of TFs on flavonoid biosynthesis, a co-expression regulation sub-network was established among TF genes and flavonoid-related structural genes according to the correlation analysis (Figure 6A). We observed that GbMYB (Gb_40628) had the highest connectivity and was closely associated with 10 structural genes. Additionally, we observed that a structural gene was regulated by multiple TFs simultaneously, such as GbCHI (Gb_21115) was positively correlated with five GbMYBs (Gb_11316, Gb_32143, Gb_33428, Gb_39081, and Gb_40628), three GbERFs (Gb_00128, Gb_26438, and Gb_37188), and one GbTrihelix (Gb_02053). Therefore, we suggested that these TFs participated in the regulation of gene expression in the flavonoid biosynthesis pathway.
Figure 6. Environmental factors may trigger the activation of TFs to control the expression of flavonoid-related structural genes. (A) The co-expression network between TFs (blue circles) and flavonoid-related structural genes (red circles). The circle size is positively correlated with the connectivity of genes in the regulatory network. (B) The correlation analysis between environmental factors and the expression of TF genes. The red and blue square indicate positive and negative correlations between specific rows and columns, respectively. Climate factors include sunshine duration (SD), relative humidity (Hum), air pressure (Press), average temperature (Tave), minimum temperature (Tmin), and maximum temperature (Tmax). (C) Cis-regulatory elements analysis of TF genes.
Further, the correlation between genes encoding TFs and climate factors was analyzed (Figure 6B). We found that several GbMYBs (Gb_11316, Gb_26833, Gb_32143, Gb_33428, and Gb_40628), GbERFs (Gb_26438 and Gb_37188), GbbHLH (Gb_17233), and GbTrihelix (Gb_02053) were significantly and negatively (r > | 0.6|, P < 0.05) correlated with mean minimum temperature (Tmin); three GbMYBs (Gb_11316, Gb_32143, and Gb_33428), one GbERF (Gb_00128), and one GbTrihelix (Gb_02053) had significant negative correlations with mean relative humidity (Hum); and two GbMYBs (Gb_11316 and Gb_32143) and one GbERF (Gb_00128) had significant positive correlations with sunshine duration (SD). These results suggested that the sunny environment was favorable to the expression of genes encoding TFs, while conversely the cold and humid environment was unfavorable to their expression. The promoter analysis also supported our hypothesis, as multiple light responsiveness (G-Box, Box 4, AE-box, I-box, L-box, Gap-box, Box II, and G-box) and low-temperature responsiveness (LTR) elements were found in the promoter regions of TF genes (Figure 6C).
Flavonoids represent one of the main classes of secondary metabolites that play an important role in plant defense against environmental stresses (e.g., temperature, precipitation, and light) (Akula and Ravishankar, 2011). Additionally, flavonoids extracted are also beneficial compounds for human health as cardioprotective, antihypertensive, and antioxidants (Fuchs et al., 2016). While studying the metabolic process of flavonoids in Ginkgo has been the subject of intense investigations (Wu et al., 2018; Meng et al., 2019; Guo et al., 2020), limited information is available at the genomic level. In the present study, we investigated the temporospatial (four leaves developmental stages and three contrasting test sites) transcriptome and metabolome dynamics biological processes to increase our understanding of Ginkgo’s flavonoids regulatory networks and to provide additional information of the molecular mechanisms of flavonoids accumulation during its leaves development. We used clonally propagated plant material, so the observed differences are attributable to either time (four leaves developmental stages) or space (three contrasting test sites).
Temporospatial Gene Expression Profiles Variation
We used the RNA-seq approach to detect changes in the studied samples gene expression profiles and the principal component analysis, PCA-plot, clearly showed a spatial separation among samples growing at different environments, suggesting that gene expression at the transcriptome level is strongly influenced/modified by environmental conditions (Figure 2A). Observations confirming previously held view that transcriptional regulatory cascades may be key components of differential resilience shown by plants to changing environments (D’esposito et al., 2017). Furthermore, the four developmental stages also were separated (PCA-plot: Figure 2A), suggesting that the observed temporal differences provided a reasonable description of transcription activities during leaf development stages (Garg et al., 2017). Interestingly and more strikingly, the transcriptome from different spatial (environmental) samples showed similar stage-specific expression patterns that were gradually separated in the same direction (PCA-1) during leaf development. These findings reflected the dynamic nature and flexibility of gene expression in response to internal (genetic) and external (environmental) cues at the transcription level during leaf development (Bar and Ori, 2014).
Flavonoids Metabolism Is Temporospatially Influenced
The observed differentially expressed genes (DEGs) among spatially different samples (environments) were identified and exhibited a collinear pattern with the environmental differences between sites. As the environmental differences between sites increase, this was accompanied by a concomitant increase in the differences of DFGs of their respective samples. For example, the difference between DEGs from YN and QJ sites was in line with the observed differences between these two sites environments (Figure 3A). Ginkgo may have developed a genetic control system as a survival strategy in response to different environments (Cho et al., 2018; Wang et al., 2020). Further, the eight temporal expression patterns (DEG sets) identified by STEM analysis (Figure 3B), contained genes significantly affected by environmental (spatial) and developmental (temporal) processes. The KEGG pathway enrichment analysis (Figure 3C) indicated that the expression of genes from Profile 0 presented a down-regulated trend with leaf development, and these genes were significantly enriched in both phenylpropanoid and flavonoid biosynthesis pathways (ko 00940 and ko 00941). These results indicated that a set of genes related to flavonoid biosyntheses, such as structural genes and TF genes, performed the stage-specific (temporally sequenced) function under external environmental stimuli. It has been reported that there is a rhythm, a time-distribution character, to the biosynthesis and metabolism of flavonoids in Ginkgo leaves (Cheng et al., 2012; Sati et al., 2013). Our results indicated that flavonoid-related metabolism was more active at the transcriptional level in the early stage of leaf development, consistent with previous studies (Yan et al., 2019; Zhu et al., 2020). These findings were also confirmed by the result of HPLC analysis; where we found that the content of TFG peaked in the early stage (Figure 1A). We also found that samples from YN and QJ sites had the greatest difference in TFG content (Figure 1B).
A Regulated Transcriptional Network for Flavonoids Biosynthesis
Previous studies have focused on the identification of flavonoid-related structural genes in Ginkgo leaves, such as genes encoding PAL, C4H, 4CL, CHS, CHI, and F3H in early flavonoid biosynthesis pathway, and genes encoding DFR, ANS, and ANR in downstream steps of the pathway (Li et al., 2018; Wu et al., 2018). In the present study, an intensive association network was observed between the expression of 15 selected structural genes and abundance of flavonoids (Figure 5B), suggesting that these structural genes (GbANR, GbANS, GbC4H, GbCHI, GbCHS, GbDFR, GbFLS, GbPAL, and GbUGT) may play crucial roles in the accumulation of specific flavonoids. More specifically, we identified one gene (Gb_41796) encoding PAL, which is an upstream key enzyme and rate limiting of the flavonoids biosynthesis pathway (Wang et al., 2014), whose expression was positively correlated with the abundance of three flavones (luteolin, apigenin, and vitexin), but negatively correlated with the two isoflavones (daidzin and genistein). Additionally, we found two genes encoding UGT (Gb_24439 and Gb_26148) whose expression was significantly and positively correlated with the abundance of flavonoid glycosides (astragalin), consistent with and supporting previous studies (Cui et al., 2016; Zhu et al., 2020). The accumulation of anthocyanins has been reported to be positively correlated with the expression of GbDFRs (Ni et al., 2020), while we found the gene encoding DFR (Gb_09086) was negatively associated with cyanidin accumulation.
Flavonoids biosynthesis is mainly regulated by transcription factors at the transcription level (Xu et al., 2014; Cao et al., 2020). We constructed a co-expression regulation network among TF genes and flavonoid-related structural genes to explore their regulatory relationship (Figure 6A). We discovered eight GbMYBs that positively regulated the expression of multiple structural genes in the flavonoid biosynthesis pathway. MYB TFs represent one of the largest families of a transcription factor in plants, involving in the regulation of different biological processes (Dubos et al., 2010). The large number of GbMYBs in the Ginkgo genome indicated that each of them may involve unique functions. Meng et al. (2019) found that the GbMYB5 was involved in the positive regulation of flavonoid biosynthesis, while Xu et al. (2014) suggested that the GbMYBF2 was responsible for repressing flavonoid biosynthesis. MYB and bHLH can act individually or in concert with other TFs to regulate a series of structural genes involved in flavonoid metabolism (Terrier et al., 2009; Carletti et al., 2013). Similarly, this co-expression network showed positive correlations between GbERF, GbbHLH, and GbTrihelix and certain structural genes associated with flavonoids.
TFs are considered as the major regulators of gene expression in response to environmental changes. MYB, ERF, and bHLH have been shown to play important roles in regulating environmental stress responses (Nakashima et al., 2009; Agarwal and Jha, 2010). In this study, we found that some GbTFs expression was negatively correlated with mean minimum temperature but positively correlated with sunshine duration (Figure 6B). Meanwhile, we also identified abundant light responsiveness elements and LTRs elements in GbTFs promoter regions (Figure 6C). It has been confirmed in previous studies, anthocyanins accumulation in Pinus contorta seedlings grown under short sunlight was significantly lower than those growing in the long sunlight area; long-term light irradiation (16 h) on leaves of Ipomoea batatas generated a dramatic increase in flavonoids content (Camm et al., 1993; Carvalho et al., 2010). As the amount of sunlight increases, there is a concomitant rise in temperature, and the composition of flavonoids in Ribes nigrum has been found to be positively correlated with temperature (Zheng et al., 2012). Our findings further support that proper control of gene expression by TFs was essential for the flavonoids biosynthesis, which played an important role in response to environmental changes (López-Maury et al., 2008).
Our investigation of the temporospatial transcriptome and metabolome dynamics biological processes provided new insights into the biosynthesis of flavonoids in Ginkgo leaves. We indicated that flavonoids content varied greatly at different developmental stages (temporally) and in different growth environments (spatially). Therefore, the careful selection of planting region(s) and optimization of leaf harvesting time are expected to substantially improve the benefits of Ginkgo utilization as a non-timber forest product. We constructed a co-expression regulation network and identified 13 TF genes having crucial roles in controlling the transcriptomic regulation of flavonoids by activating the expression of multiple structural genes. These results provide candidate genes for future enhancement of flavonoids production by genetic strategies in Ginkgo. Furthermore, the large amount of data resources generated will serve as the foundation for a system biology approach to study the dynamics of leaf development and flavonoids accumulation in other plants.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA657336, SUB7912734.
GW and YE-K conceived the study. YG collected the field samples, analyzed the data, and drafted the manuscript. GW, YE-K, TW, and F-FF modified the manuscript. All the authors have approved the manuscript.
This research was funded by the National Natural Science Foundation of China (grant no. 31971689), the National Key Research and Development Program of China (2017YFD0600700), and the Doctorate Fellowship Foundation of Nanjing Forestry University.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We are grateful to B. Ratcliffe for critical review and suggestions.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.589326/full#supplementary-material
Supplementary Figure 1 | UHPLC chromatogram of standards and samples at different sampling stages. (A) Standards; (B–E) samples collected from the Yi Ning (YN) site in May, June, July, and August, respectively.
Supplementary Figure 2 | The qRT-PCR analysis results of key genes involved in the flavonoid biosynthesis. The FPKM (A) and relative expression levels (B) were presented in bar plot with the error bar by calculating the mean and stand deviation (SD) of three independent replicates.
Supplementary Figure 3 | Bubble maps show the top 20 significantly enriched KEGG pathways among the genes in three modules highly associated with flavonoid accumulation.
Supplementary Table 1 | Summary of the sequencing quality of 36 RNA libraries of Ginkgo leaves.
Supplementary Table 2 | Primers used for the qRT-PCR assay.
Supplementary Table 3 | Summary of the quantitative values (peak area) of the identified 12 flavonoids.
- ^ http://gigadb.org/dataset/100613
- ^ https://cran.r-project.org/package=pheatmap
- ^ http://bioinformatics.psb.ugent.be/webtools/plantcare/html/
- ^ http://data.cma.cn/
Camm, E., Mccallum, J., Leaf, E., and Koupai-Abyazani, M. (1993). Cold-induced purpling of Pinus contorta seedlings depends on previous daylength treatment. Plant Cell Environ. 16, 761–764. doi: 10.1111/j.1365-3040.1993.tb00497.x
Carletti, G., Lucini, L., Busconi, M., Marocco, A., and Bernardi, J. (2013). Insight into the role of anthocyanin biosynthesis-related genes in Medicago truncatula mutants impaired in pigmentation in leaves. Plant Physiol. Biochem. 70, 123–132. doi: 10.1016/j.plaphy.2013.05.030
Carvalho, I. S., Cavaco, T., Carvalho, L. M., and Duque, P. (2010). Effect of photoperiod on flavonoid pathway activity in sweet potato (Ipomoea batatas (L.) Lam.) leaves. Food Chem. 118, 384–390. doi: 10.1016/j.foodchem.2009.05.005
Cheng, S., Feng, X., Linling, L., Cheng, H., and Zhang, W. (2012). Seasonal pattern of flavonoid content and related enzyme activities in leaves of Ginkgo biloba L. Notulae Botanicae Horti Agrobotanici Cluj-Napoca 40, 98–106. doi: 10.15835/nbha4017262
Cho, S. M., Lee, H., Jo, H., Lee, H., Kang, Y., Park, H., et al. (2018). Comparative transcriptome analysis of field-and chamber-grown samples of Colobanthus quitensis (Kunth) Bartl, an Antarctic flowering plant. Sci. Rep. 8, 1–14. doi: 10.1038/s41598-018-29335-4
Cui, L., Yao, S., Dai, X., Yin, Q., Liu, Y., Jiang, X., et al. (2016). Identification of UDP-glycosyltransferases involved in the biosynthesis of astringent taste compounds in tea (Camellia sinensis). J. Exp. Bot. 67, 2285–2297. doi: 10.1093/jxb/erw053
D’esposito, D., Ferriello, F., Dal Molin, A., Diretto, G., Sacco, A., Minio, A., et al. (2017). Unraveling the complexity of transcriptomic, metabolomic and quality environmental response of tomato fruit. BMC Plant Biology 17:66. doi: 10.1186/s12870-017-1008-4
Ellnain-Wojtaszek, M., Kruczyński, Z., and Kasprzak, J. (2002). Variations in the free radical scavenging activity of Ginkgo biloba L. leaves in the period of complete development of green leaves to fall of yellow ones. Food Chem. 79, 79–84. doi: 10.1016/S0308-8146(02)00181-4
Fuchs, D., Nyakayiru, J., Draijer, R., Mulder, T. P., Hopman, M. T., Eijsvogels, T. M., et al. (2016). Impact of flavonoid-rich black tea and beetroot juice on postprandial peripheral vascular resistance and glucose homeostasis in obese, insulin-resistant men: a randomized controlled trial. Nutr. Metab. 13:34. doi: 10.1186/s12986-016-0094-x
Garg, R., Singh, V. K., Rajkumar, M. S., Kumar, V., and Jain, M. (2017). Global transcriptome and coexpression network analyses reveal cultivar-specific molecular signatures associated with seed development and seed size/weight determination in chickpea. Plant J. 91, 1088–1107. doi: 10.1111/tpj.13621
Gruenwald, J., Eckert, A., and Kressig, R. W. (2020). The effects of standardized Ginkgo biloba extracts (GBE) on subjective cognitive decline (SCD) in middle-aged adults: a review. Adv. Aging Res. 9, 45–65. doi: 10.4236/aar.2020.93005
Guo, Y., Gao, C., Wang, M., Fu, F., El-Kassaby, Y. A., Wang, T., et al. (2020). Metabolome and transcriptome analyses reveal flavonoids biosynthesis differences in Ginkgo biloba associated with environmental conditions. Ind. Crops Prod. 158:112963. doi: 10.1016/j.indcrop.2020.112963
Li, W., Yang, S., Lu, Z., He, Z., Ye, Y., Zhao, B., et al. (2018). Cytological, physiological, and transcriptomic analyses of golden leaf coloration in Ginkgo biloba L. Hortic. Res. 5, 1–14. doi: 10.1038/s41438-018-0015-4
López-Maury, L., Marguerat, S., and Bähler, J. (2008). Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation. Nat. Rev. Genet. 9, 583–593. doi: 10.1038/nrg2398
Meng, J., Wang, B., He, G., Wang, Y., Tang, X., Wang, S., et al. (2019). Metabolomics integrated with transcriptomics reveals redirection of the phenylpropanoids Metabolic flux in Ginkgo biloba. J. Agric. Food Chem. 67, 3284–3291. doi: 10.1021/acs.jafc.8b06355
Nakashima, K., Ito, Y., and Yamaguchi-Shinozaki, K. (2009). Transcriptional regulatory networks in response to abiotic stresses in Arabidopsis and grasses. Plant Physiol. 149, 88–95. doi: 10.1104/pp.108.129791
Ni, J., Dong, L., Jiang, Z., Yang, X., Sun, Z., Li, J., et al. (2018). Salicylic acid-induced flavonoid accumulation in Ginkgo biloba leaves is dependent on red and far-red light. Ind. Crops Product. 118, 102–110. doi: 10.1016/j.indcrop.2018.03.044
Ni, J., Hao, J., Jiang, Z., Zhan, X., Dong, L., Yang, X., et al. (2017). NaCl induces flavonoid biosynthesis through a putative novel pathway in post-harvest Ginkgo leaves. Front. Plant Sci. 8:920. doi: 10.3389/fpls.2017.00920
Ni, J., Ruan, R., Wang, L., Jiang, Z., Gu, X., Chen, L., et al. (2020). Functional and correlation analyses of dihydroflavonol-4-reductase genes indicate their roles in regulating anthocyanin changes in Ginkgo biloba. Ind. Crops Product. 152:112546. doi: 10.1016/j.indcrop.2020.112546
Petrussa, E., Braidot, E., Zancani, M., Peresson, C., Bertolini, A., Patui, S., et al. (2013). Plant flavonoids-biosynthesis, transport and involvement in stress responses. Int. J Mol. Sci. 14, 14950–14973. doi: 10.3390/ijms140714950
Sati, P., Pandey, A., Rawat, S., and Rani, A. (2013). Phytochemicals and antioxidants in leaf extracts of Ginkgo biloba with reference to location, seasonal variation and solvent system. J. Pharm. Res. 7, 804–809. doi: 10.1016/j.jopr.2013.09.001
Sato, Y., Antonio, B., Namiki, N., Motoyama, R., Sugimoto, K., Takehisa, H., et al. (2011). Field transcriptome revealed critical developmental and physiological transitions involved in the expression of growth potential in japonicarice. BMC Plant Biol. 11:10. doi: 10.1186/1471-2229-11-10
Terrier, N., Torregrosa, L., Ageorges, A., Vialet, S., Verries, C., Cheynier, V., et al. (2009). Ectopic expression of VvMybPA2 promotes proanthocyanidin biosynthesis in grapevine and suggests additional targets in the pathway. Plant Physiol. 149, 1028–1041. doi: 10.1104/pp.108.131862
Tian, J., Liu, Y., and Chen, K. (2017). Ginkgo biloba extract in vascular protection: molecular mechanisms and clinical applications. Curr. Vasc. Pharmacol. 15, 532–548. doi: 10.2174/1570161115666170713095545
Wang, G., Cao, F., Li, C., Guo, X., Wang, J., and Lu, Z. (2014). Temperature has more effects than soil moisture on biosynthesis of flavonoids in Ginkgo (Ginkgo biloba L.) leaves. New For. 45, 797–812. doi: 10.1007/s11056-014-9437-5
Wang, L., Cui, J., Jin, B., Zhao, J., Xu, H., Lu, Z., et al. (2020). Multifeature analyses of vascular cambial cells reveal longevity mechanisms in old Ginkgo biloba trees. Proc. Natl. Acad. Sci. U.S.A. 117, 2201–2210. doi: 10.1073/pnas.1916548117
Wang, L., Shi, H., Wu, J., and Cao, F. (2016). Alternative partial root-zone irrigation enhances leaf flavonoid accumulation and water use efficiency of Ginkgo biloba. New For. 47, 377–391. doi: 10.1007/s11056-015-9521-5
Weckwerth, W. (2008). Integration of metabolomics and proteomics in molecular plant physiology–coping with the complexity by data-dimensionality reduction. Physiol. Plant. 132, 176–189. doi: 10.1111/j.1399-3054.2007.01011.x
Wu, D., Feng, J., Lai, M., Ouyang, J., Liao, D., Yu, W., et al. (2020). Combined application of bud and leaf growth fertilizer improves leaf flavonoids yield of Ginkgo biloba. Ind. Crops Product. 150:112379. doi: 10.1016/j.indcrop.2020.112379
Wu, Y., Guo, J., Zhou, Q., Xin, Y., Wang, G., and Xu, L.-A. (2018). De novo transcriptome analysis revealed genes involved in flavonoid biosynthesis, transport and regulation in Ginkgo biloba. Ind. Crops Product. 124, 226–235. doi: 10.1016/j.indcrop.2018.07.060
Xu, F., Ning, Y., Zhang, W., Liao, Y., Li, L., Cheng, H., et al. (2014). An R2R3-MYB transcription factor as a negative regulator of the flavonoid biosynthesis pathway in Ginkgo biloba. Funct. Integr. Genom. 14, 177–189. doi: 10.1007/s10142-013-0352-1
Yan, J., Yu, L., He, L., Zhu, L., Xu, S., Wan, Y., et al. (2019). Comparative transcriptome analysis of celery leaf Blades identified an R2R3-MYB transcription factor that regulates apigenin metabolism. J. Agric. Food Chem. 67, 5265–5277. doi: 10.1021/acs.jafc.9b01052
Ye, J., Cheng, S., Zhou, X., Chen, Z., Kim, S. U., Tan, J., et al. (2019). A global survey of full-length transcriptome of Ginkgo biloba reveals transcript variants involved in flavonoid biosynthesis. Ind. Crops Product. 139:111547. doi: 10.1016/j.indcrop.2019.111547
Zheng, J., Yang, B., Ruusunen, V., Laaksonen, O., Tahvonen, R., Hellsten, J., et al. (2012). Compositional differences of phenolic compounds between black currant (Ribes nigrum L.) cultivars and their response to latitude and weather conditions. J. Agric. Food Chem. 60, 6581–6593. doi: 10.1021/jf3012739
Zhu, J., Xu, Q., Zhao, S., Xia, X., Yan, X., An, Y., et al. (2020). Comprehensive co-expression analysis provides novel insights into temporal variation of flavonoids in fresh leaves of the tea plant (Camellia sinensis). Plant Sci. 290:110306. doi: 10.1016/j.plantsci.2019.110306
Zou, K., Liu, X., Zhang, D., Yang, Q., Fu, S., Meng, D., et al. (2019). Flavonoid biosynthesis is likely more susceptible to elevation and tree age than other branch pathways involved in phenylpropanoid biosynthesis in Ginkgo leaves. Front. Plant Sci. 10:983. doi: 10.3389/fpls.2019.00983
Keywords: Ginkgo biloba, flavonoids biosynthesis, leaf development, transcriptome dynamics, temporospatial variation
Citation: Guo Y, Wang T, Fu F-F, El-Kassaby YA and Wang G (2020) Temporospatial Flavonoids Metabolism Variation in Ginkgo biloba Leaves. Front. Genet. 11:589326. doi: 10.3389/fgene.2020.589326
Received: 30 July 2020; Accepted: 05 November 2020;
Published: 27 November 2020.
Edited by:Gorji Marzban, University of Natural Resources and Life Sciences, Vienna, Austria
Copyright © 2020 Guo, Wang, Fu, El-Kassaby and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.