Search
Submit Manuscript
5.1CiteScore
3.2Impact Factor
Advanced Search
MAP  >  Epigenetics Insights
2024 Volume 17
Article Contents
Abstract
Introduction
Materials and methods
Plant materials
Protein extraction
Trypsin digestion
Affinity enrichment of ubiquitinated peptides
LC-MS/MS and data analysis
Bioinformatics methods
Results
Overview of the ubiquitomic analysis
Properties of ubiquitinated proteins and motifs
A temperature-regulated rice endosperm ubiquitome
Functions for differentially expressed ubiquitinated proteins
Discussion
Ubiquitination of heat shock proteins and transcription factors under high-temperature stress
Sucrose and starch metabolism are modified by ubiquitination under high-temperature stress
Ubiquitination affects sucrose hydrolysis
Ubiquitination of ADP-glucose pyrophosphorylase (AGPase)
Ubiquitination regulates amylose biosynthesis
Ubiquitination regulates amylopectin biosynthesis
Ubiquitination affects ribosomal systems under high-temperature stress
Author contributions
Data availability
Acknowledgments
Conflict of interest
Supplementary information
Rights and permissions
References
About this article
Next Previous
ARTICLE   Open Access    

Genetic and epigenetic variations underlying flavonoid divergence in Beihua and Sijihua honeysuckles

  • 1.

    Horticulture Research Institute, Sichuan Academy of Agricultural Sciences, Chengdu 610066, Sichuan, China

  • 2.

    State Key Laboratory for Conservation and Utilization of Subtropical Agro-Bioresources, College of Agriculture, Guangxi University, Nanning 530004, China

  • 3.

    Yazhouwan National Laboratory, Sanya 572024, Hainan, China

  • # Authors contributed equally: Xianyun Yu, Hang Yu

More Information
  • Received: 30 July 2024
    Revised: 02 September 2024
    Accepted: 24 September 2024
    Published online: 11 October 2024
    Epigenetics Insights  17,  Article number: e002 (2024)  |  Cite this article
  • Flavonoids are important antibacterial and antiviral active substances which are the most crucial medicinal components of honeysuckle. However, the content of medicinally active substances in different honeysuckle cultivars is significantly different. Genetic variations and epigenetics play essential roles in plant evolution and trait improvement. Here, we performed multi-omics sequencing of two honeysuckle cultivars (Beihua and Sijihua) at different stages (WB and GF). The results revealed 9,909,981 SNPs in the genomes of the two cultivars, and 12,688 high-impact SNPs were found to regulate genes involved in important biological pathways, such as plant stress resistance. Furthermore, it was found that the majority of differentially methylated cytosines (DMCs, 81%) between Beihua and Sijihua were associated with SNPs. SNP-related DMCs were associated with 76% of the genes, among which 3,325 DEGs (e.g., LjPAL, LjCHI, and LjFLS) were significantly enriched in the flavonoid biosynthesis pathway. The presence of a large number of SNP-related DMCs in the flanking and gene regions of these genes may have led to the overexpression of the genes in Beihua, which increased the accumulation of flavonoids in Beihua. In summary, the present study provides theoretical and technical support for improving the genetic and epigenetic traits of honeysuckle.

  • Introduction

    Rice (Oryza sativa L.) is a world staple crop that feeds over half of the world's population[1,2]. In recent years, high-temperature events are becoming more frequent and intensive as a result of global warming, which can severely affect rice grain yield and quality[3,4]. During flowering and grain filling stages, high-temperature stress can result in a significant reduction in seed setting rate and influence amylose content, starch fine structure, functional properties and chalkiness degree of rice[5−7]. Transcriptome and proteome analysis in rice endosperm have also been used to demonstrate the differences in high-temperature environments at gene and protein expression levels[8−10]. In addition, as an important post-translational modification, protein phosphorylation has proven to be involved in the regulation of starch metabolism in response to high-temperature stress[11]. However, little is known about whether protein ubiquitination regulates seed development under high-temperature stress.

    Ubiquitination is another form of post-translational modification that plays key roles in diverse cellular processes[12]. Several reports have described the functions of ubiquitination in rice defense responses based on ubiquitome analysis. Liu et al.[13] investigated relationships between ubiquitination and salt-stress responses in rice seedlings using a gel-based shotgun proteomic analysis and revealed the potential important role of protein ubiquitination in salt tolerance in rice. Xie et al.[14] identified 861 peptides with ubiquitinated lysines in 464 proteins in rice leaf cells by combining highly sensitive immune affinity purification and high resolution liquid chromatography-tandem mass spectrometry (LC-MS/MS). These ubiquitinated proteins regulated a wide range of processes, including response to stress stimuli. A later study revealed the relationships between ubiquitination status and activation of rice defense responses, and generated an in-depth quantitative proteomic catalog of ubiquitination in rice seedlings after chitin and flg22 treatments, providing useful information for understanding how the ubiquitination system regulates the defense responses upon pathogen attack[15]. Although many studies have shown that ubiquitination plays improtant roles in the heat response of plant[16,17], there has been little systematic discussion on the ubiquitome of rice endosperm in the context of global climate change.

    In this study, we examine the high-temperature induced ubiquitination change in two rice varieties with different starch qualities, through a label-free quantitative ubiquitome analysis. This study provides a comprehensive view of the function of ubiquitination in high-temperature response of rice developing seed, which will shed new light on the improvement of rice grain quality under heat stress.

    Materials and methods

    Plant materials

    Two indica rice varieties with different starch quality, 9311 and Guangluai4 (GLA4), were used as materials. 9311 is a heat-sensitive variety, which displays low amylose content with good starch quality; while GLA4 is known to be the parental variety of HT54, an indica breeding line with heat tolerance, and thus GLA4 is possibly heat tolerant, which shows high amylose content with poor starch quality[18,19]. Rice growth conditions, sample treatment and collection were conducted as previously described[11].

    Protein extraction

    Husk, pericarp and embryo were detached from immature rice grains on ice[20]. Rice endosperm was then ground with liquid nitrogen, and the cell powder was sonicated 10 rounds of 10 s sonication and 15 s off-sonication on ice in lysis buffer (6 M Guanidine hydrochloride, pH 7.8−8.0, 0.1% Protease Inhibitor Cocktail) using a high intensity ultrasonic processor. Following lysis, the suspension was centrifuged at 14,000 g for 40 min at 4 °C to remove the debris. The supernatant was collected, and the protein concentration was estimated using BCA assay (Pierce BCA Protein assay kit, Thermo Fisher Scientific, Waltham, MA, USA) before further analysis.

    Trypsin digestion

    The protein mixture was reduced by DTT with the final concentration of 10 mM at 37 °C for 1 h, alkylated by iodoacetamide with a final concentration of 50 mM at room temperature in the dark for 0.5 h, and digested by trypsin (1:50) at 37 °C for 16 h. Then the sample was diluted by adding trifluoroacetic acid (TFA) to the final concentration of 0.1%. The enzymatic peptides were desalted on a Sep-Pak C18 cartridge (Waters, Milford, MA, USA), concentrated by lyophilization and reconstituted in precooled IAP buffer (50 mM MOPS-NaOH PH 7.2, 10 mM Na2HPO4, 50 mM NaCl) for further analysis.

    Affinity enrichment of ubiquitinated peptides

    The peptides solution was incubated with prewashed K-ε-GG antibody beads (PTMScan Ubiquitin Remnant Motif (K-ε-GG) Kit), and gently shaken at 4 °C for 1.5 h. The suspension was centrifuged at 2,000 g for 30 s, and the supernatant was removed. The Anti-K-ε-GG antibody beads were washed with IAP Buffer three times and with ddH2O three times. The peptides were eluted from the beads with 0.15% trifluoroacetic acid (TFA). Finally, the eluted fractions were combined and desalted with C18 Stage Tips.

    LC-MS/MS and data analysis

    LC-MS/MS analysis were performed using the methods of Pang et al.[11]. Raw mass spectrometric data were analyzed with MaxQuant software (version 1.3.0.5) and were compared with the indica rice protein sequence database (Oryza sativa subsp. indica-ASM465v1). Parameters were set according to Pang et al.[11]. All measurements were obtained from three separate biological replicates.

    Bioinformatics methods

    Quantification of the modified peptides was performed using the label-free quantification (LFQ) algorithm[11]. Differentially ubiquitinated sites (proteins) in response to high-temperature were identified by Student's t-test (p < 0.05, log2(fold-change) > 1) with at least two valid values in any condition or the ubiquitination sites that exhibited valid values in one condition (at least two of three replicates) and none in the other.

    Subcellular localization was performed using CELLO database (http://cello.life.nctu.edu.tw). Gene Ontology (GO) annotation proteome was derived from the AgriGO (http://bioinfo.cau.edu.cn/agriGO/). The differential metabolic profiles were visualized with MapMan software (version 3.6.0RC1). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation was performed by using KEGG Automatic Annotation Server (KAAS) software. A p-value of < 0.05 was used as the threshold of significant enrichment. SWISS-MODEL was used to generate the tertiary structure of GBSSI (SWISS-MODEL, http://swissmodel.expasy.org/). The figures were annotated with Adobe Illustrator (Adobe Systems, San Jose, CA, USA).

    Results

    Overview of the ubiquitomic analysis

    To elucidate how high-temperature stress influences rice developing endosperm at the ubiquitination level, a label-free analysis was performed to quantify ubiquitome from two indica rice varieties under normal (9311-C and GLA4-C) and high-temperature conditions (9311-H and GLA4-H). The distribution of mass error of all the identified peptides was near zero and most of them (71.5%) were between -1 and 1 ppm, suggesting that the mass accuracy of the MS data fits the requirement (Fig. 1a). Meanwhile, the length of most peptides distributed between 8 and 42, ensuring that sample preparation reached standard conditions (Fig. 1b).

    Figure 1.  Characteristics of the ubiquitinated proteome of rice endosperm and QC validation of MS data. (a) Mass error distribution of all identified ubiquitinated peptides. (b) Peptide length distribution. (c) Frequency distribution of ubiquitinated proteins according to the number of ubiquitination sites identified.
    DownLoad: Full-Size Img PowerPoint

    In all endosperm samples, a total of 437 ubiquitinated peptides were identified from 246 ubiquitinated proteins, covering 488 quantifiable ubiquitinated sites (Supplemental Table S1). Among the ubiquitinated proteins, 60.6% had only one ubiquitinated lysine site, and 18.7%, 8.1%, 5.3%, or 7.3% had two, three, four, or five and more ubiquitinated sites, respectively. In addition, four proteins (1.6%, BGIOSGA004052, BGIOSGA006533, BGIOSGA006780, BGIOSGA022241) were ubiquitinated at 10 or more lysine sites (Fig. 1c, Supplemental Table S1). The proteins BGIOSGA008317 had the most ubiquitination sites with the number of 16. It was noted that besides ubiquitin, two related ubiquitin-like proteins NEDD8 and ISG15 also contain C-terminal di-Gly motifs generated by trypsin cleavage, and the modifications of these three proteins cannot be distinguished by MS[21]. Here, the di-Gly-modified proteome therefore represents a composite of proteins modified by these three proteins. However, the sites from NEDD8 or ISG15 modifications were limited because they mediate only a few reactions in cells[21].

    Properties of ubiquitinated proteins and motifs

    To better understand the lysine ubiquitome changes in rice endosperm induced by high-temperature, we performed a Gene Ontology (GO) functional annotation analysis on all identified ubiquitinated proteins (Fig. 2a). In the biological process GO category, 'metabolic process' and 'cellular process' were mainly enriched, accounting for 75.1% and 74.1% of ubiquitinated proteins, respectively. In addition, 34.6% proteins were associated with 'response to stimulu', emphasizing the regulatory role of ubiquitination modification in response to high-temperature stress. From the cellular component perspective, ubiquitinated proteins were mainly associated with 'cellular anatomical entity' (99.4%), 'intracellular' (84.4%) and 'protein-containing complex' (29.9%). The molecular function category result suggested that these proteins were largely involved in 'binding' (62.7%), 'catalytic activity' (43.4%) and 'structural molecule activity' (16.5%). Furthermore, subcellular location annotation information indicated that 34.7%−39.4% proteins were located in the cytoplasm, and other were mostly located in the nucleus (23.5%−27.7%), plasma membrane (9.4%−11.4%), and chloroplast (9.6%−12.8%) (Fig. 2b). It is noteworthy that the ubiquitinated proteins located in the cytoplasm were decreased in high-temperature environments in both varieties.

    Figure 2.  Analysis of ubiquitinated proteins and motifs. (a) Gene ontology (GO) functional characterization of ubiquitinated proteins. (b) Subcellular localization of ubiquitinated proteins. From the inside out, the ring represents 9311-C, 9311-H, GLA4-C and GLA4-H, respectively. (c) Motif enrichment analysis of ubiquitinated proteins.
    DownLoad: Full-Size Img PowerPoint

    The following two significantly enriched motifs from all of the identified ubiquitinated sites were identified using MoDL analysis: [A/S]xKub and Kubxx[E/Q/R/V]x[E/G/L/P/Q/R/Y], which covered 84 and 100 sequences, respectively (Fig. 2c). Further analysis showed that the conserved alanine (A) and glutamic acid (E) were included in upstream and downstream of the ubiquitinated lysine sites in rice endosperm. A similar phenomenon also occurred in rice leaf[14], wheat leaf[22], and petunia[23], indicating that alanine (A) and glutamic acid (E) were likely to be the specific amino acids in conserved ubiquitination motifs in plants. Additionally, serine (S) was enriched at the position -2 (upstream) of the ubiquitinated lysine, while various amino acids such as arginine (R), glutamic acid (E), glutamine (Q), valine (V) were found at positions +3 and +5 (downstream).

    A temperature-regulated rice endosperm ubiquitome

    To detect possible changes in rice endosperm ubiquitome attributable to high-temperature stress, we then performed LFQ analysis on all quantifiable ubiquitination sites within our dataset. As shown in Fig. 3a, more ubiquitinated proteins, peptides and sites were detected in the treatment groups (9311-H and GLA4-H), suggesting that exposure to high-temperature stress may increase the ubiquitination events in rice endosperm. Only 282 common ubiquitinated sites in 158 proteins were quantifiable for all sample groups due to reversible ubiquitination induced by high-temperature (Fig. 3b). Principal component analysis (PCA) showed that three repeats of each sample clustered together, and four groups were clearly separated (Fig. 3c). Furthermore, the differentially expression profiles of ubiquitination sites (proteins) in 9311 and GLA4 under high-temperature stress were depicted to further understand the possible changes (Fig. 3d). Where LFQ values were missing, the data were filtered to identify those ubiquitination sites with a consistent presence/absence expression pattern. These analyses yielded 89 ubiquitination sites that were only present in 9311-H and six that were only present in 9311-C (Fig. 3d, Supplemental Table S2). Similarly, 51 differentially expressed ubiquitination sites were present in GLA4-H and 13 ubiquitination sites only occurred in GLA4-C (Fig. 3d & Supplemental Table S3). Beyond that, a total of 113 and 50 significantly changed ubiquitination sites (p < 0.05, log2(fold-change) >1) were screened out in 9311 and GLA4, respectively (Fig. 3d, Supplemental Tables S4 & S5). For subsequent comparative analysis, the ubiquitination expression profiles with consistent presence/absence and ubiquitination sites with significant differences in statistical testing were combined and named as 9311-Up, 9311-Down, GLA4-Up, and GLA4-Down, respectively (Fig. 3d). The number of significantly up-regulated ubiquitination sites was far greater than down-regulated ubiquitination sites in both 9311 and GLA4 varieties. These findings indicated that high temperature not only induced the occurrence of ubiquitination sites, but also significantly upregulated the intensity of ubiquitination. Beyond that, the magnitude of the up-regulation in 9311 was higher than that in GLA4 (Fig. 3b & d), indicating that the ubiquitination modification of heat-sensitive varieties was more active than heat-resistant varieties in response to high-temperature stress.

    Figure 3.  A temperature regulated rice endosperm ubiquitome. (a) The number of ubiquitinated proteins, peptides and sites detected in four group samples. (b) Venn diagram of ubiquitination sites (proteins) detected in four group samples. (c) PCA based on ubiquitination intensity across all four sample groups with three biological repetitions. (d) Differentially expression profiles of ubiquitination sites (proteins) in 9311 and GLA4 under high-temperature stress. The expression profiles of selected ubiquitination sites (p < 0.05, log2(fold-change) >1) were normalized using the Z-score and presented in a heatmap.
    DownLoad: Full-Size Img PowerPoint

    To further investigate the ubiquitination regulatory pattern under high temperature stress in two varieties, four groups with significantly regulated sites were analyzed. There were 37 ubiquitination sites showed the same regulatory trend in 9311 and GLA4, accounting for 17.8% and 32.5% of the total differentially expressed sites in 9311 and GLA4, respectively. Among them, 36 ubiquitination sites were upregulated and one site was downregulated (Fig. 4a). In addition, 159 upregulated ubiquitination sites and three downregulated sites were only present in 9311, while 53 upregulated sites and 15 downregulated sites were only present in GLA4. Moreover, nine ubiquitination sites showed opposite regulatory trends in 9311 and GLA4. A similar regulatory trend of ubiquitination proteins is shown in Fig. 4b. It is noted that some proteins had both upregulated and downregulated ubiquitination sites (Supplemental Tables S6 & S7), indicating that significant differences in ubiquitination were, to some extent, independent of protein abundance.

    Figure 4.  Comparison of differentially ubiquitinated sites and proteins in 9311 and GLA4 under high-temperature stress.
    DownLoad: Full-Size Img PowerPoint

    Functions for differentially expressed ubiquitinated proteins

    To understand the function of ubiquitination in response to the high-temperature stress of rice endosperm, we conducted GO enrichment-based clustering analysis of the differentially ubiquitinated proteins in 9311 and GLA4 at high temperature, respectively (Fig. 5). In the biological process category of 9311, proteins were relatively enriched in the carbohydrate metabolic process, polysaccharide metabolic process, starch biosynthetic process, cellular macromolecule localization, protein localization, intracellular transport, and phosphorylation (Fig. 5). For the molecular function analysis, we found that the proteins related to kinase activity, nucleotidyltransferase activity, phosphotransferase activity, and nutrient reservoir activity were enriched (Fig. 5). The two principal cellular components were intrinsic component of membrane and integral component of membrane (Fig. 5). There was no significantly enriched GO term in the GLA4 group due to the dataset containing relatively few proteins, and thus, further enrichment analysis was conducted on the proteins that were common to both varieties. The results showed that proteins were over-represented in carbon metabolism, including starch biosynthesis and metabolism, glucan biosynthesis and metabolism, and polysaccharide biosynthesis and metabolism (Fig. 5), indicating the importance of carbohydrate synthesis and metabolism in the ubiquitination regulatory network.

    Figure 5.  Enrichment analysis of differentially expressed ubiquitinated proteins based on Gene Ontology (GO) terms.
    DownLoad: Full-Size Img PowerPoint

    To identify pathways which were differentially ubiquitinated under high-temperature stress, the KEGG pathway-based clustering analysis was conducted. The results showed that the differentially ubiquitinated proteins in both 9311 and GLA4 were mostly abundant in the pathways of carbohydrate metabolism, starch and sucrose metabolism, folding, sorting and degradation, translation, ribosome, and protein processing in endoplasmic reticulum (Fig. 6a). In the 9311 group, the pathways of carbohydrate metabolism, starch and sucrose metabolism, glycosyltransferases, glycolysis, and energy metabolism were enriched in the differentially ubiquitinated proteins (Fig. 6b); while there was no significantly enriched KEGG pathway in the GLA4 group. We further found the proteins that were common to both varieties were only significantly enriched in the starch and sucrose metabolism pathways (p = 0.04). The ubiquitination proteins involved in the starch and sucrose metabolism mainly include: sucrose hydrolysis (SUS, FK, UGPase), and starch synthesis (AGPase, GBSSI, BEI, BEIIb, PUL, PHO1), which are discussed below.

    Figure 6.  KEGG classification and enrichment analysis of differentially ubiquitinated proteins. (a) Number of differentially ubiquitinated proteins based on KEGG classification in 9311 and GLA4. (b) KEGG enrichment analysis of differentially ubiquitinated proteins in 9311.
    DownLoad: Full-Size Img PowerPoint

    Discussion

    Although many reports have described specific examples of ubiquitination in rice defense responses[13,15,16], our knowledge on global changes in the developing endosperm ubiquitome under high-temperature stress is still lacking. In this study, a label-free quantitative proteomic analysis of ubiquitination was applied to examine the high-temperature induced ubiquitination change of two indica rice varieties (9311 and GLA4) with distinct starch quality. We identified many new lysine modification sites on proteins involved in various pathways, highlighting the complexity of the ubiquitination-mediated regulatory system in high-temperature stress responses in rice.

    Ubiquitination of heat shock proteins and transcription factors under high-temperature stress

    Heat shock proteins accumulate under various stresses and play important roles in plant defenses against abiotic stresses[24,25]. Research has shown that a number of heat shock proteins were prominent in the rice ubiquitome network, of which OsHSP71.1, and OsHSP82A showed increased ubiquitination levels under chitin and flg22 treatment[15]. Here, seven lysine residues on five heat shock proteins possessed ubiquitination modification in rice endosperm. Three sites (BGIOSGA011420-K78, BGIOSGA026764-K99, BGIOSGA029594-K106) showed significant up-regulation in 9311 under high-temperature stress, while in GLA4, the ubiquitination level of BGIOSGA011420-K78 was down-regulated. This differential ubiquitination of heat-tolerant and heat-sensitive varieties provided a basis for studying the regulation of post-translational modification of heat shock proteins under high-temperature stress, despite the regulatory role of those heat shock proteins being still unclear.

    Transcription factors (TFs) play an essential role in the regulation of gene expression. A total of three transcription factors were identified in the ubiquitination dataset, of which two were NAC family members. As one of the largest plant-specific TF families, NAC is involved in the responses to abiotic and biotic stresses[26]. The ubiquitination modification of K173 in BGIOSGA018048 was specifically expressed in the 9311-H group, which may affect the stress resistance level in high-temperature environments. In addition, two sites K148 and K149 of ERF TF family member BGIOSGA024035 was downregulated in GLA4 under high-temperature stress. This differential ubiquitination was likely to affect the expression of related genes regulated by this transcription factor.

    Sucrose and starch metabolism are modified by ubiquitination under high-temperature stress

    Ubiquitination affects sucrose hydrolysis

    The results of GO and KEGG enrichment analysis of differential ubiquitination proteins indicated that the sucrose and starch metabolic pathway was largely affected by ubiquitination regulation under high-temperature stress (Figs 5 & 6). The ubiquitination sites involved in sucrose and starch metabolism are listed in Table 1. To assess how high-temperature stress affects the crucial pathway, the significantly differential ubiquitination sites in 9311 and GLA4 were displayed in the heatmap of specific proteins (Fig. 7).

    Table 1.  Ubiquitination sites related to sucrose and starch metabolism in rice endosperm.
    Gene nameAnnotationProtein entryModification site(s)
    SUS1Sucrose synthase 1BGIOSGA010570K172, K177
    SUS2Sucrose synthase 2BGIOSGA021739K160, K165, K176, K804
    SUS3Sucrose synthase 3BGIOSGA026140K172, K177, K541, K544, K588
    FKFructokinaseBGIOSGA027875K143
    UGPaseUDP-glucose pyrophosphorylaseBGIOSGA031231K27, K150, K303, K306
    AGPS1ADP-glucose pyrophosphorylase small subunit 1BGIOSGA030039K94, K464, K484
    AGPS2ADP-glucose pyrophosphorylase small subunit 2BGIOSGA027135K106, K132, K385, K403, K406, K476, K496
    AGPL2ADP-glucose pyrophosphorylase large subunit 2BGIOSGA004052K41, K78, K134, K191, K227, K254, K316, K338, K394, K396, K463, K508, K513
    AGPL3ADP-glucose pyrophosphorylase large subunit 3BGIOSGA017490K509
    GBSSIGranule bound starch synthase IBGIOSGA022241K130, K173, K177, K181, K192, K258, K371, K381, K385, K399, K462, K517, K530, K549, K571, K575
    BEIStarch branching enzyme IBGIOSGA020506K103, K108, K122
    BEIIbStarch branching enzyme IIbBGIOSGA006344K134
    PULStarch debranching enzyme:PullulanaseBGIOSGA015875K230, K330, K431, K736, K884
    PHO1Plastidial phosphorylaseBGIOSGA009780K277, K445, K941
     | Show Table
    DownLoad: CSV
    Figure 7.  Sucrose and starch pathway at the ubiquitination levels in rice endosperm under high-temperature stress.
    DownLoad: Full-Size Img PowerPoint

    In cereal endosperm, sucrose is the substrate for the biosynthesis of starch. The formation of glucose 1-phosphate (G1P, used in starch synthesis, see below) from sucrose requires a series of enzymes[27]. Here, we found that sucrose synthase 1 (SUS1), SUS2, SUS3, fructokinase (FK), and UDP-glucose pyrophosphorylase (UGPase) were ubiquitinated among all sample groups (Table 1, Fig. 7). In the ubiquitome of seedling and leaf in japonica rice, ubiquitination sites have been found in SUS1, SUS2, UGPase, and FK, which were related to sucrose hydrolysis[14,15]. SUS catalyzed the process of cleaving sucrose into UDP-glucose (UDPG) and fructose. Two ubiquitination sites, K172 and K177, were identified in SUS1 in rice endosperm, which were also found in rice leaves[14]. A total of four ubiquitination sites were identified in SUS2, two of which were also reported in rice seedling and leaf, indicating the conservation of the lysine residues in different rice tissues. It was noted that all four ubiquitination sites in SUS2 were upregulated in high-temperature environments, although the regulated sites of 9311 and GLA4 were different. In 9311, the ubiquitination levels of K160, K174, and K804 were increased, while GLA4 was only upregulated in K176. The ubiquitination sites K541, K544, and K588 in SUS3 were screened from developing rice seeds for the first time. In addition, SUS3 had two completely overlapping sites K172 and K177 with SUS1, and it was difficult to determine which enzymes the two sites belonged to. The ubiquitination levels of SUS3-K541 and SUS3-K544 in 9311 significantly increased in high-temperature environments, while there was no significant difference in the ubiquitination level of SUS3 in GLA4. Overall, the ubiquitination sites of SUS in rice endosperm were located in the functional domain except for SUS2-K804, reflecting the importance of ubiquitination regulation in SUS.

    UGPase catalyses the conversion of glucose 1-phosphate and UTP into UDPG[28]. Research has shown that the mutation of UGPase gene lead to chalky endosperm[29]. As shown in Table 1, four ubiquitination sites K27, K150, K303, and K306 were identified in rice endosperm, which were completely inconsistent with the seven ubiquitination sites in rice seedlings and two in leaves[14,15], reflecting the tissue specificity. We speculated that the UGPase with different modification sites may play different regulatory roles in metabolic pathways in different tissues. Under high-temperature stress, the ubiquitination level of UGPase-K27 was 8.1-fold up-regulated. Liao et al.[10] demonstrated that the expression of UDPase was down-regulated in both heat-tolerant and heat-sensitive rice lines under high temperature conditions, which could reasonably explain the significant up-regulation of UGPase-K27 ubiquitination level. The ubiquitination site K143 of FK was also reported in seedling tissues[15].

    Ubiquitination of ADP-glucose pyrophosphorylase (AGPase)

    The AGPase reaction represents the first committed step of starch biosynthesis[27]. A total of 22 lysine ubiquitination sites were identified in four AGPase subunits (AGPL2, AGPL3, AGPS1, AGPS2). AGPL2 had 13 ubiquitination sites, of which six were located in NTP_transferase domain, including K254, K338, K191, K134, K227, and K316. High-temperature stress resulted in an increase in the ubiquitination level of K254 in both 9311 and GLA4, and significant upregulation of K508 and K513 in 9311, as well as K191, K227, and K316 in GLA4. In contrast, AGPL2-K394 were significantly downregulated in GLA4. AGPL3 contained one ubiquitination site K509, and the modification level of AGPL3-K509 was up-regulated in high-temperature environments in 9311. AGPS1 had one specific ubiquitination site K464 and another two sites K94 and K484 that completely overlapped with AGPS2-K106 and AGPS2-K496, respectively. The modification levels of K464 and K484 significantly increased in high-temperature environments in 9311, and K94 was significantly up-regulated in both varieties. There were seven ubiquitination sites in AGPS2 in rice endosperm, which were different with the sites found in rice leaves[14]. In addition to the two sites that overlapped with AGPS1, AGPS2 had another two ubiquitination sites (K406 and K132) that upregulated in high-temperature environments.

    Ubiquitination regulates amylose biosynthesis

    Amylose content is one of the key determinants that strongly influence rice grain quality[30]. The biosynthesis of amylose requires the catalytic effect of granule-bound starch synthase I (GBSSI)[30,31]. Here, a total of 16 ubiquitination sites were identified in GBSSI (Table 1, Fig. 8a). Among these ubiquitination sites, six lysine residues (K130, K173, K177, K181, K192, K258) were located in glycosyltransferase 5 (GT5) domain, and three sites (K399, K462, K517) were located in GT1 domain (Fig. 6), indicating the important role of ubiquitination regulation of GBSSI. Under high-temperature stress, the ubiquitination levels of six sites (K130, K177, K399, K381, K385, K549) increased in two indica rice varieties, while one sites (K258) showed downregulation in 9311 (Fig. 8a). Numerous studies had described that the amylose content was reduced under high-temperature stress in rice[5,7], which might be due to the degradation of GBSSI proteins caused by the increased significantly up-regulated ubiquitination sites. These ubiquitination sites identified in rice GBSSI with significant differences under high-temperature stress were expected to become a new breakthrough point for the improvement of starch quality.

    Figure 8.  Structure of GBSSI. (a) Domain structure of GBSSI and ubiquitination sites with significant differences in response to high-temperature stress. (b) 3D model of GBSSI and the relationship between ubiquitination sites K462 and ADP, SO4 (salt bridge or hydrogen bond).
    DownLoad: Full-Size Img PowerPoint

    To further determine the regulatory role of the ubiquitination sites in GBSSI, SWISS-MODEL was used to predict 3D structural model. As shown in Fig. 8b, GBSSI had three SO4 (sulfate ions) and one ADP ligand. These ligands interact with GBSSI through hydrogen bonds and salt bridges. Three sites, K447, R458, and K462, were associated with SO4 through salt bridges, while G100, N265, Q412, K462, and Q493 interact with the hydrogen bonds of ADP in GBSSI[32,33]. Based on this finding, it can be reasonably inferred that the K462 site with ubiquitination modification located in the GT1 domain played an important role in the interaction between GBSSI, SO4, and ADP. An in-depth investigation was necessary to gain a more comprehensive understanding of the regulatory function of ubiquitination modification at GBSSI-K462, although there was no significant difference in the ubiquitination level under high-temperature stress.

    Ubiquitination regulates amylopectin biosynthesis

    Amylopectin, the major component of starch, is synthesized by the coordinated action of multiple enzymes including soluble starch synthase (SSs), starch branching enzyme (BEs), starch debranching enzyme (DBEs), and phosphorylases (PHOs or Phos) with ADPG as a substrate. In this study, ubiquitination sites were detected in BEs, DBEs, and Phos.

    BEs, covering two isoforms, BEI and BEII, are responsible for catalyzing the formation of α-1,6-glucosidic linkages of amylopectin[34]. There were three ubiquitination sites (K103, K108, and K122) identified in BEI (Fig. 9a). K122 was the first amino acid in the carbohydrate-binding module 48 (CBM48) domain. Sequence alignment analyses of BEs from eight plants revealed K122 was conserved among all plants' BEI (Fig. 9a), suggesting a high probability of the functional effects of ubiquitination modification at this site. In high-temperature environments, ubiquitination levels of K108 and K122 were significantly up-regulated in 9311, while no significantly regulated ubiquitination sites of BEI were observed in GLA4. Only one ubiquitination site, K134, was found in BEIIb (Fig. 9a). The ubiquitination levels showed a slightly upward trend with no significant differences in high-temperature environments in both varieties. These changes could be one of the reasons for increased gelatinization temperature and relative crystallinity of rice starch in response to high-temperature[5].

    Figure 9.  Domain structure of (a) BEs, (b) PUL and (c) Pho1 as well as their ubiquitination sites with significant differences in response to heat stress. Residues in red indicate the ubiquitination site. Non-ubiquitinated residues are shown in dark grey.
    DownLoad: Full-Size Img PowerPoint

    DBEs consists of isoamylase (ISA) and pullulanase (PUL) with catalytic function for hydrolyzing α-1,6-glucosic linkages[35]. In the present study, we found that only PUL was ubiquitinated in rice endosperm (Fig. 9b). Among five ubiquitination sites (K230, K330, K432, K736, and K884) identified in PUL, K230 was located in the PULN2 domain, while K330 was in the CBM48 domain. Under high-temperature stress, K330 showed completely opposite regulatory trends in two cultivars. In addition, the ubiquitination level of K884, located in the DUF3372 domain, was significantly up-regulated in 9311. Previous study has reported that the expression of PUL was significantly up-regulated in 9311 under high-temperature stress, while GLA4 showed down-regulation in PUL abundance[11]. Consequently, there might be two possible functions of these ubiquitination sites. One possibility is that ubiquitination sites were unrelated to protein degradation; instead, they regulated the biosynthesis of amylopectin by affecting other functions of the protein. Secondly, ubiquitination sites were associated with protein degradation, and the levels of ubiquitination modification were based on protein abundance, resulting in a completely consistent regulation of ubiquitination modification and protein abundance under high-temperature stress.

    PHOs, including two types, Pho1/PHO1 and Pho2/PHO2, are responsible for the transfer of glucosyl units from Glc-1-P to the non-reducing end of a-1,4-linked glucan chains[36]. Pho1 is a temperature-dependent enzyme and considered crucial not only during the maturation of amylopectin but also in the initiation process of starch synthesis[37,38]. The three ubiquitination sites (K277, K445, K941) identified in Pho1 were located in two phosphorylase domains. We found that two sites, Pho1-K277 and Pho1-K445, were only ubiquitinated in high-temperature environments in 9311 and GLA4, respectively. Pang et al.[11] has demonstrated that the protein abundance of Pho1 decreased under high-temperature stress, especially in GLA4. Satoh et al.[38] reported that the functional activity of Pho1 was weakened under conditions of high temperature and its function might be complement by one or more other factors. Hence, these ubiquitination modifications that specifically occurred in high-temperature environments might be related to the degradation of Pho1 proteins.

    Ubiquitination affects ribosomal systems under high-temperature stress

    As a factory for protein synthesis in cells, the ribosome is an extremely crucial structure in the cell[39]. It has been proven that multiple ribosomal subunits were abundantly ubiquitinated in Arabidopsis and wheat[22]. In the present study, 57 ubiquitination sites involving 33 ribosome subunits were identified in 40S and 60S ribosome complexes in rice. Under high-temperature stress, the ubiquitination levels of some sites were significantly upregulated or downregulated, implying that ubiquitination of ribosomal proteins is likely to be an important regulatory mechanism in high-temperature response in rice endosperm. The results of GO and KEGG enrichment analysis indicated that the ribosome system was one of the most active systems for ubiquitination regulation under high-temperature stress. We speculated that the ubiquitin-proteasome system might be involved in the removal of subunits or entire ribosomes that were improperly folded in high-temperature environments. As shown in Fig. 10, the S10e, L18Ae, S27Ae, L9e, S3e, S28e, S20e, and S2e subunits were significantly up-regulated in 9311, while L13e subunits showed a completely opposite regulatory trend at the ubiquitination sites K81 and K88. In GLA4, the ubiquitination levels of S10e, S27Ae, L10Ae, L9e, S3e, S2e, and L4e showed a significant increase, while the ubiquitination level of L17e was significantly down-regulated under high-temperature stress. A total of seven ubiquitination sites involving S10e, S27Ae, L9e, S3e, and S2e subunits were jointly up-regulated in both two varieties. These sites might be related to the degradation of improperly folded ribosome subunits under high-temperature stress, while other ubiquitination sites with variety specificity might be associated with ribosomal function.

    Figure 10.  Ribosome system at the ubiquitination levels in rice endosperm under high-temperature stress. Grey shadings represent ubiquitinated proteins with no significant differences under heat stress. Red and blue shadings indicate up-regulated and down-regulated ubiquitinated proteins, respectively. Orange shading displays a combination of up- and down-regulated ubiquitinated sites in the same ubiquitinated protein.
    DownLoad: Full-Size Img PowerPoint

    In conclusion, this study provides the first comprehensive view of the ubiquitome in rice developing endosperm, and demonstrated that ubiquitination has diverse functions in the high-temperature response of rice endosperm by modulating various cellular processes, especially the sucrose and starch metabolism. Comparative analysis of the temperature-induced ubiquitination status revealed some similarities and more interesting differences between 9311 and GLA4. These differences might be the reason for the different qualities formation of the two indica rice varieties, which could provide potential genetic resources for the improvement of the heat resistance in rice. Considering the diversity of ubiquitination modification, it is worthwhile to further validate and explore the function and regulatory mechanism of the key targets and key pathways. The findings provide valuable insights into the role of ubiquitination in response to high-temperature stress and lay a foundation for further functional analysis of lysine ubiquitination in rice.

    Author contributions

    The authors confirm contribution to the paper as follows: study conception and design: Bao J, Pang Y; data collection: Pang Y; analysis and interpretation of results: Pang Y; draft manuscript preparation: Ying Y; Revised manuscript preparation: Ying Y, Pang Y, Bao J. All authors reviewed the results and approved the final version of the manuscript.

    Data availability

    The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

    Acknowledgments

    This work was financially supported by the AgroST Project (NK2022050102) and Zhejiang Provincial Natural Science Foundation (LZ21C130003).

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

  • Supplemental Table S1 Summary of Sijihua and Beihua RNA-seq reads mapping.
    Supplemental Table S2 Summary of whole genome resequencing mapping.
    Supplemental Table S3 Classification of SNPs based on genome location.
    Supplemental Table S4 Number of genes affected by SNP mutation.
    Supplemental Table S5 Summary of Sijihua and Beihua BS-seq reads mapping.
    Supplemental Table S6 Mining and analyzing flavonoid pathway genes of honeysuckle.
    Supplemental Fig. S1 The significance analysis between Beihua WB and GF, as well as Sijihua WB and GF.
    Supplemental Fig. S2 DNA methylation landscape of Beihua and Sijihua.
    Supplemental Fig. S3 The correlation between SNPs and DNA methylation.
    Supplemental Fig. S4 DNA methylation patterns of protein-coding genes and TEs in Sijihua and Beihua in the WB.
    Supplemental Fig. S5 DNA methylation patterns of protein-coding genes and TEs in Sijihua and Beihua in the GF.
    Supplemental Fig. S6 SNP lead to changes in the type and level of DNA methylation and characterization of the diversity of DMC sites.
    Supplemental Fig. S7 CG≥TG-type mutations in relation to CG methylation and gene expression during the GF stage.
    Supplemental Fig. S8 GO enrichment analysis of CG≥TG-related genes.
    Supplemental Fig. S9 The pattern of differentially expressed genes related to SNP-associated DMCs.
    Supplemental Fig. S10 The comparison of DNA methylation levels of 35 key genes of the flavonoid pathway in two stages of Sijihua and Beihua.
    Supplemental Fig. S11 Heatmaps showing the expression levels of 35 key genes, the number of SNPs in the gene body of the genes, the number of SNP-associated DMCs, and the number of non-SNP-associated DMCs in Beihua and Sijihua during the WB and GF stages (* represents genes that showed significant differences in both stages).
    Supplemental Fig. S12 Genome browser showing DNA methylation and expression levels of LjFLS (EVM0027194) and representative screenshots of two types of DMCs in the promoter region (Chr01:50581502-50581602).
  • [1]

    Saito K, Yonekura-Sakakibara K, Nakabayashi R, Higashi Y, Yamazaki M, et al. 2013. The flavonoid biosynthetic pathway in Arabidopsis: structural and genetic diversity. Plant Physiology and Biochemmistry 72:21−34

    doi: 10.1016/j.plaphy.2013.02.001

    CrossRef   Google Scholar

    [2]

    Andersen JR, Zein I, Wenzel G, Darnhofer B, Eder J, et al. 2008. Characterization of phenylpropanoid pathway genes within European maize (Zea mays L.) inbreds. BMC Plant Biology 8:2

    doi: 10.1186/1471-2229-8-2

    CrossRef   Google Scholar

    [3]

    Hoang VL, Innes DJ, Shaw PN, Monteith GR, Gidley MJ, et al. 2015. Sequence diversity and differential expression of major phenylpropanoid-flavonoid biosynthetic genes among three mango varieties. BMC Genomics 16:561

    doi: 10.1186/s12864-015-1784-x

    CrossRef   Google Scholar

    [4]

    Kato M, Miura A, Bender J, Jacobsen SE, Kakutani T. 2003. Role of CG and non-CG methylation in immobilization of transposons in Arabidopsis. Current Biology 13:421−26

    doi: 10.1016/S0960-9822(03)00106-4

    CrossRef   Google Scholar

    [5]

    Adato A, Mandel T, Mintz-Oron S, Venger I, Levy D, et al. 2009. Fruit-surface flavonoid accumulation in tomato is controlled by a SlMYB12-regulated transcriptional network. PLoS Genetics 5:e1000777

    doi: 10.1371/journal.pgen.1000777

    CrossRef   Google Scholar

    [6]

    Huang J, Zhang C, Zhao X, Fei Z, Wan K, et al. 2016. The Jujube Genome Provides Insights into Genome Evolution and the Domestication of Sweetness/Acidity Taste in Fruit Trees. PLoS Genetics 12:e1006433

    doi: 10.1371/journal.pgen.1006433

    CrossRef   Google Scholar

    [7]

    Huang X, Kurata N, Wei X, Wang ZX, Wang A, et al. 2012. A map of rice genome variation reveals the origin of cultivated rice. Nature 490:497−501

    doi: 10.1038/nature11532

    CrossRef   Google Scholar

    [8]

    Hufford MB, Xu X, van Heerwaarden J, Pyhäjärvi T, Chia JM, et al. 2012. Comparative population genomics of maize domestication and improvement. Nature Genetics 44:808−11

    doi: 10.1038/ng.2309

    CrossRef   Google Scholar

    [9]

    Lei Y, Yang L, Duan S, Ning S, Li D, et al. 2022. Whole-genome resequencing reveals the origin of tea in Lincang. Frontiers in Plant Science 13:984422

    doi: 10.3389/fpls.2022.984422

    CrossRef   Google Scholar

    [10]

    Ren G, Zhang X, Li Y, Ridout K, Serrano-Serrano ML, et al. 2021. Large-scale whole-genome resequencing unravels the domestication history of Cannabis sativa. Science Advances 7:eabg2286

    doi: 10.1126/sciadv.abg2286

    CrossRef   Google Scholar

    [11]

    Zhao H, Sun S, Ding Y, Wang Y, Yue X, et al. 2021. Analysis of 427 genomes reveals moso bamboo population structure and genetic basis of property traits. Nature Communications 12:5466

    doi: 10.1038/s41467-021-25795-x

    CrossRef   Google Scholar

    [12]

    Kaeppler SM, Kaeppler HF, Rhee Y. 2000. Epigenetic aspects of somaclonal variation in plants. Plant Molecular Biology 43:179−88

    doi: 10.1023/A:1006423110134

    CrossRef   Google Scholar

    [13]

    Iwasaki M, Paszkowski J. 2014. Epigenetic memory in plants. The EMBO Journal 33:1987−98

    doi: 10.15252/embj.201488883

    CrossRef   Google Scholar

    [14]

    Jones MJ, Goodman SJ, Kobor MS. 2015. DNA methylation and healthy human aging. Aging Cell 14:924−32

    doi: 10.1111/acel.12349

    CrossRef   Google Scholar

    [15]

    Kulis M, Esteller M. 2010. DNA methylation and cancer. Advances in Genetics 70:27−56

    doi: 10.1016/B978-0-12-380866-0.60002-2

    CrossRef   Google Scholar

    [16]

    Weber M, Schübeler D. 2007. Genomic patterns of DNA methylation: targets and function of an epigenetic mark. Current Opinion in Cell Biology 19:273−80

    doi: 10.1016/j.ceb.2007.04.011

    CrossRef   Google Scholar

    [17]

    Lin L, Wang S, Zhang J, Song X, Zhang D, et al. 2022. Integrative analysis of transcriptome and metabolome reveals the effect of DNA methylation of chalcone isomerase gene in promoter region on Lithocarpus polystachyus Rehd flavonoids. Synthetic and Systems Biotechnology 7:928−40

    doi: 10.1016/j.synbio.2022.05.003

    CrossRef   Google Scholar

    [18]

    Strygina K, Khlestkina E. 2022. Flavonoid biosynthesis genes in Triticum aestivum L.: methylation patterns in cis-regulatory regions of the duplicated CHI and F3H genes. Biomolecules 12:689

    doi: 10.3390/biom12050689

    CrossRef   Google Scholar

    [19]

    Jia H, Jia H, Lu S, Zhang Z, Su Z, et al. 2022. DNA and histone methylation regulates different types of fruit ripening by transcriptome and proteome analyses. Journal of Agricultural and Food Chemistry 70:3541−56

    doi: 10.1021/acs.jafc.1c06391

    CrossRef   Google Scholar

    [20]

    An YQC, Goettel W, Han Q, Bartels A, Liu Z, et al. 2017. Dynamic changes of genome-wide DNA methylation during soybean seed development. Scientific Reports 7:12263

    doi: 10.1038/s41598-017-12510-4

    CrossRef   Google Scholar

    [21]

    Huang H, Liu R, Niu Q, Tang K, Zhang B, et al. 2019. Global increase in DNA methylation during orange fruit development and ripening. Proceedings of the National Academy of Sciences of the United States of America 116:1430−36

    doi: 10.1073/pnas.1815441116

    CrossRef   Google Scholar

    [22]

    Povilus RA, Friedman WE. 2022. Transcriptomes across fertilization and seed development in the water lily Nymphaea thermarum (Nymphaeales): evidence for epigenetic patterning during reproduction. Plant Reproduction 35:161−78

    doi: 10.1007/s00497-022-00438-3

    CrossRef   Google Scholar

    [23]

    Ossowski S, Schneeberger K, Lucas-Lledó JI, Warthmann N, Clark RM, et al. 2010. The rate and molecular spectrum of spontaneous mutations in Arabidopsis thaliana. Science 327:92−94

    doi: 10.1126/science.1180677

    CrossRef   Google Scholar

    [24]

    Kiefer C, Willing EM, Jiao WB, Sun H, Piednoël M, et al. 2019. Interspecies association mapping links reduced CG to TG substitution rates to the loss of gene-body methylation. Nature Plants 5:846−55

    doi: 10.1038/s41477-019-0486-9

    CrossRef   Google Scholar

    [25]

    Xue Y, Shi Y, Qi Y, Yu H, Zou C, et al. 2022. Epigenetic and Genetic Contribution for Expression Bias of Homologous Alleles in Polyploid Sugarcane. Agronomy 12:2852

    doi: 10.3390/agronomy12112852

    CrossRef   Google Scholar

    [26]

    Zhong Z, Feng S, Mansfeld BN, Ke Y, Qi W, et al. 2023. Haplotype-resolved DNA methylome of African cassava genome. Plant Biotech nology Journal 21:247−49

    doi: 10.1111/pbi.13955

    CrossRef   Google Scholar

    [27]

    Becker C, Hagmann J, Müller J, Koenig D, Stegle O, et al. 2011. Spontaneous epigenetic variation in the Arabidopsis thaliana methylome. Nature 480:245−49

    doi: 10.1038/nature10555

    CrossRef   Google Scholar

    [28]

    He L, Xu X, Li Y, Li C, Zhu Y, et al. 2013. Transcriptome analysis of buds and leaves using 454 pyrosequencing to discover genes associated with the biosynthesis of active ingredients in Lonicera japonica Thunb. PLoS One 8:e62922

    doi: 10.1371/journal.pone.0062922

    CrossRef   Google Scholar

    [29]

    Liu T, Yang J, Liu S, Zhao Y, Zhou J, et al. 2020. Regulation of chlorogenic acid, flavonoid, and iridoid biosynthesis by histone H3K4 and H3K9 methylation in Lonicera japonica. Molecular Biology Reports 47:9301−11

    doi: 10.1007/s11033-020-05990-7

    CrossRef   Google Scholar

    [30]

    Huang W, Xiong L, Zhang L, Zhang F, Han X, et al. 2022. Study on content variation of flavonoids in different germplasm during development of Lonicerae Japonicae Flos. Chinese Traditional and Herbal Drugs 53:3156−64

    doi: 10.7501/j.issn.0253-2670.2022.10.026

    CrossRef   Google Scholar

    [31]

    Yu H, Cui N, Guo K, Xu W, Wang H. 2023. Epigenetic changes in the regulation of carotenoid metabolism during honeysuckle flower development. Horticultural Plant Journal 9:577−88

    doi: 10.1016/j.hpj.2022.11.003

    CrossRef   Google Scholar

    [32]

    Yu H, Guo K, Lai K, Shah MA, Xu Z, et al. 2022. Chromosome-scale genome assembly of an important medicinal plant honeysuckle. Scientific Data 9:226

    doi: 10.1038/s41597-022-01385-4

    CrossRef   Google Scholar

    [33]

    Xanthopoulou A, Manioudaki M, Bazakos C, Kissoudis C, Farsakoglou AM, et al. 2020. Whole genome re-sequencing of sweet cherry (Prunus avium L.) yields insights into genomic diversity of a fruit species. Horticulture Research 7:60

    doi: 10.1038/s41438-020-0281-9

    CrossRef   Google Scholar

    [34]

    Xu Q, Wu L, Luo Z, Zhang M, Lai J, et al. 2022. DNA demethylation affects imprinted gene expression in maize endosperm. Genome Biology 23:77

    doi: 10.1186/s13059-022-02641-x

    CrossRef   Google Scholar

    [35]

    Wang ZH, Zhang D, Bai Y, Zhang YH, Liu Y, et al. 2013. Genomewide variation in an introgression line of rice-Zizania revealed by whole-genome re-sequencing. PLoS One 8:e74479

    doi: 10.1371/journal.pone.0074479

    CrossRef   Google Scholar

    [36]

    Wang H, Beyene G, Zhai J, Feng S, Fahlgren N, et al. 2015. CG gene body DNA methylation changes and evolution of duplicated genes in cassava. Proceedings of the National Academy of Sciences of the United States of America 112:13729−34

    doi: 10.1073/pnas.1519067112

    CrossRef   Google Scholar

    [37]

    Gent JI, Ellis NA, Guo L, Harkess AE, Yao Y, et al. 2013. CHH islands: de novo DNA methylation in near-gene chromatin regulation in maize. Genome Research 23:628−37

    doi: 10.1101/gr.146985.112

    CrossRef   Google Scholar

    [38]

    Schmitz RJ, He Y, Valdés-López O, Khan SM, Joshi T, et al. 2013. Epigenome-wide inheritance of cytosine methylation variants in a recombinant inbred population. Genome Research 23:1663−74

    doi: 10.1101/gr.152538.112

    CrossRef   Google Scholar

    [39]

    Schmitz RJ, Schultz MD, Lewsey MG, O'Malley RC, Urich MA, et al. 2011. Transgenerational epigenetic instability is a source of novel methylation variants. Science 334:369−73

    doi: 10.1126/science.1212959

    CrossRef   Google Scholar

    [40]

    Selvaraj S, Krishnaswamy S, Devashya V, Sethuraman S, Krishnan UM. 2014. Flavonoid–metal ion complexes: a novel class of therapeutic agents. Medicinal Research Reviews 34:677−702

    doi: 10.1002/med.21301

    CrossRef   Google Scholar

    [41]

    Zhong S, Fei Z, Chen YR, Zheng Y, Huang M, et al. 2013. Single-base resolution methylomes of tomato fruit development reveal epigenome modifications associated with ripening. Nature Biotechnology 31:154−59

    doi: 10.1038/nbt.2462

    CrossRef   Google Scholar

    [42]

    Zheng X, Wang T, Cheng T, Zhao L, Zheng X, et al. 2022. Genomic variation reveals demographic history and biological adaptation of the ancient relictual, lotus (Nelumbo Adans). Horticulture Research 9:uhac029

    doi: 10.1093/hr/uhac029

    CrossRef   Google Scholar

    [43]

    Zilberman D, Coleman-Derr D, Ballinger T, Henikoff S. 2008. Histone H2A.Z and DNA methylation are mutually antagonistic chromatin marks. Nature 456:125−29

    doi: 10.1038/nature07324

    CrossRef   Google Scholar

    [44]

    Bewick AJ, Ji L, Niederhuth CE, Willing EM, Hofmeister BT, et al. 2016. On the origin and evolutionary consequences of gene body DNA methylation. Proceedings of the National Academy of Sciences of the United States of America 113:9111−16

    doi: 10.1073/pnas.1604666113

    CrossRef   Google Scholar

    [45]

    Kim KD, El Baidouri M, Abernathy B, Iwata-Otsubo A, Chavarro C, et al. 2015. A Comparative Epigenomic Analysis of Polyploidy-Derived Genes in Soybean and Common Bean. Plant Physiology 168:1433−47

    doi: 10.1104/pp.15.00408

    CrossRef   Google Scholar

    [46]

    Wang ZL, Wang S, Kuang Y, Hu ZM, Qiao X, Ye M. 2018. A comprehensive review on phytochemistry, pharmacology, and flavonoid biosynthesis of Scutellaria baicalensis. Pharmaceutical Biology 56:465−84

    doi: 10.1080/13880209.2018.1492620

    CrossRef   Google Scholar

    [47]

    Ji H, Shin Y, Lee C, Oh H, Yoon IS, et al. 2021. Genomic Variation in Korean japonica Rice Varieties. Genes 12:1749

    doi: 10.3390/genes12111749

    CrossRef   Google Scholar

    [48]

    Li R, Maioli A, Lanteri S, Moglia A, Bai Y, et al. 2023. Genomic analysis highlights putative defective susceptibility genes in tomato germplasm. Plants 12:2289

    doi: 10.3390/plants12122289

    CrossRef   Google Scholar

    [49]

    Skarzyńska A, Pawełkowicz M, Pląder W. 2021. Influence of transgenesis on genome variability in cucumber lines with a thaumatin II gene. Physiology and Molecular Biology of Plants 27:985−96

    doi: 10.1007/s12298-021-00990-8

    CrossRef   Google Scholar

    [50]

    Cui Y, Ge Q, Zhao P, Chen W, Sang X, et al. 2021. Rapid mining of candidate genes for verticillium wilt resistance in cotton based on BSA-Seq analysis. Frontiers in Plant Science 12:703011

    doi: 10.3389/fpls.2021.703011

    CrossRef   Google Scholar

    [51]

    Mas-Gómez J, Cantín CM, Moreno MÁ, Martínez-García PJ. 2022. Genetic diversity and genome-wide association study of morphological and quality traits in peach using two Spanish peach germplasm collections. Frontiers in Plant Science 13:854770

    doi: 10.3389/fpls.2022.854770

    CrossRef   Google Scholar

    [52]

    Eichten SR, Stuart T, Srivastava A, Lister R, Borevitz JO. 2016. DNA methylation profiles of diverse Brachypodium distachyon align with underlying genetic diversity. Genome Research 26:1520−31

    doi: 10.1101/gr.205468.116

    CrossRef   Google Scholar

    [53]

    Hu W, Ji C, Shi H, Liang Z, Ding Z, et al. 2021. Allele-defined genome reveals biallelic differentiation during cassava evolution. Molecular Plant 14:851−54

    doi: 10.1016/j.molp.2021.04.009

    CrossRef   Google Scholar

    [54]

    Yin YC, Zhang XD, Gao ZQ, Hu T, Liu Y. 2019. The research progress of chalcone isomerase (CHI) in plants. Molecular Biotechnology 61:32−52

    doi: 10.1007/s12033-018-0130-3

    CrossRef   Google Scholar

    [55]

    Jiang W, Yin Q, Wu R, Zheng G, Liu J, et al. 2015. Role of a chalcone isomerase-like protein in flavonoid biosynthesis in Arabidopsis thaliana. Journal of Experimental Botany 66:7165−79

    doi: 10.1093/jxb/erv413

    CrossRef   Google Scholar

    [56]

    Wang M, Zhang Y, Zhu C, Yao X, Zheng Z, et al. 2021. EkFLS overexpression promotes flavonoid accumulation and abiotic stress tolerance in plant. Plant Physiology 172:1966−82

    doi: 10.1111/ppl.13407

    CrossRef   Google Scholar

    [57]

    Jia D, Li Z, Dang Q, Shang L, Shen J, et al. 2020. Anthocyanin biosynthesis and methylation of the MdMYB10 promoter are associated with the red blushed-skin mutant in the red striped-skin "Changfu 2" apple. Journal of Agricultural and Food Chemistry 68:4292−304

    doi: 10.1021/acs.jafc.9b07098

    CrossRef   Google Scholar

    [58]

    Muir SR, Collins GJ, Robinson S, Hughes S, Bovy A, et al. 2001. Overexpression of petunia chalcone isomerase in tomato results in fruit containing increased levels of flavonols. Nature Biotechnology 19:470−74

    doi: 10.1038/88150

    CrossRef   Google Scholar

    [59]

    Yuan Y, Zuo J, Zhang H, Li R, Yu M, et al. 2022. Integration of Transcriptome and Metabolome Provides New Insights to Flavonoids Biosynthesis in Dendrobium huoshanense. Frontiers in Plant Science 13:850090

    doi: 10.3389/fpls.2022.850090

    CrossRef   Google Scholar

    [60]

    Schilbert HM, Schöne M, Baier T, Busche M, Viehöver P, et al. 2021. Characterization of the Brassica napus Flavonol Synthase Gene Family Reveals Bifunctional Flavonol Synthases. Frontiers in Plant Science 12:733762

    doi: 10.3389/fpls.2021.733762

    CrossRef   Google Scholar

  • Cite this article

    Yu X, Yu H, Lu Y, Zhang C, Wang H. 2024. Genetic and epigenetic variations underlying flavonoid divergence in Beihua and Sijihua honeysuckles. Epigenetics Insights 17: e002 doi: 10.48130/epi-0024-0002
    Yu X, Yu H, Lu Y, Zhang C, Wang H. 2024. Genetic and epigenetic variations underlying flavonoid divergence in Beihua and Sijihua honeysuckles. Epigenetics Insights 17: e002 doi: 10.48130/epi-0024-0002
    shu

Figures(5)

Download PDF

Article Metrics

Article views(1164) PDF downloads(207)

Other Articles By Authors

ARTICLE   Open Access    

Genetic and epigenetic variations underlying flavonoid divergence in Beihua and Sijihua honeysuckles

  • 1.

    Horticulture Research Institute, Sichuan Academy of Agricultural Sciences, Chengdu 610066, Sichuan, China

  • 2.

    State Key Laboratory for Conservation and Utilization of Subtropical Agro-Bioresources, College of Agriculture, Guangxi University, Nanning 530004, China

  • 3.

    Yazhouwan National Laboratory, Sanya 572024, Hainan, China

  • Received: 30 July 2024
  • Revised: 02 September 2024
  • Accepted: 24 September 2024
  • Published online: 11 October 2024
Epigenetics Insights  17,  Article number: e002  (2024)  |  Cite this article

Abstract: Flavonoids are important antibacterial and antiviral active substances which are the most crucial medicinal components of honeysuckle. However, the content of medicinally active substances in different honeysuckle cultivars is significantly different. Genetic variations and epigenetics play essential roles in plant evolution and trait improvement. Here, we performed multi-omics sequencing of two honeysuckle cultivars (Beihua and Sijihua) at different stages (WB and GF). The results revealed 9,909,981 SNPs in the genomes of the two cultivars, and 12,688 high-impact SNPs were found to regulate genes involved in important biological pathways, such as plant stress resistance. Furthermore, it was found that the majority of differentially methylated cytosines (DMCs, 81%) between Beihua and Sijihua were associated with SNPs. SNP-related DMCs were associated with 76% of the genes, among which 3,325 DEGs (e.g., LjPAL, LjCHI, and LjFLS) were significantly enriched in the flavonoid biosynthesis pathway. The presence of a large number of SNP-related DMCs in the flanking and gene regions of these genes may have led to the overexpression of the genes in Beihua, which increased the accumulation of flavonoids in Beihua. In summary, the present study provides theoretical and technical support for improving the genetic and epigenetic traits of honeysuckle.

    HTML

    Introduction

    • Genetic variations and epigenetics play important roles in regulating flavonoid accumulation in plants. Genetic variation refers to DNA sequence differences present in the genomes of different individuals, impacting coding sequences, promoter regions, and regulatory elements. These differences can directly or indirectly influence gene expression levels and regulatory patterns, thereby impacting the plant's phenotype. The most common genetic variation is single nucleotide polymorphism (SNP), resulting from changes in a single nucleotide in the DNA sequence. They can be classified as deletion, insertion, or substitution mutations. Research has indicated that SNPs in coding regions can alter amino acid sequences, impacting protein function. Numerous studies on the diversity of gene sequences related to flavonoid biosynthesis have been conducted in plants such as Arabidopsis thaliana[1], maize[2], mango[3], barley[4], and tomato[5]. The results revealed a significant correlation between SNP content in genes associated with flavonoid biosynthesis and the accumulation of flavonoid compounds. Furthermore, genetic variations lead to phenotypic differences among cultivars of important crops like jujube[6], camellia[7], hemp[8], bamboo[9], rice[10], and maize[11]. Plant phenotypic changes are influenced not only by genetic variations but also by epigenetic factors[12]. Epigenetics refers to heritable changes in gene expression patterns without altering the DNA sequence, encompassing DNA methylation, histone modifications, and chromatin remodeling[13]. DNA methylation is a conserved epigenetic marker across eukaryotes, and it is involved in many important biological processes, such as human aging[14], cancer development[15], genome integrity[16], gene imprinting, and plant stress responses. Studies have shown that DNA methylation impacts the accumulation of flavonoid compounds in plants. For instance, in Lithocarpus polystachyus Rehd[17], Triticum aestivum L.[18] and tomatoes[19], DNA methylation occurring in promoter and coding regions regulates gene expression associated with flavonoid biosynthesis, impacting flavonoid compound accumulation. DNA methylation plays a pivotal role in crucial biological processes in seeds[20], embryos[21], and flowers[22]. Consequently, both genetic and epigenetic variations significantly contribute to flavonoid compound accumulation in plants.

      However, sequence variations and epigenetic modifications are not independent phenomena, they interact intricately and culminate in phenotypic transformation. Methylated cytosines are more prone to induce a higher mutation rate at the sequence level, leading to changes in methylation levels and consequently influencing gene expression. Previous studies have shown that in Arabidopsis[23], lotus[24], and sugarcane[25], methylated cytosines are more susceptible to mutation compared to unmethylated cytosines, resulting in CG=>TG type variations and causing changes in methylation that impact plant gene expression. However, SNPs can also lead to changes in the sequence type or methylation levels of DNA. When cytosine is mutated to other types of bases (A, T, and G), the methylation pattern or level can also change, resulting in a difference in methylation levels at these sites[26]. Similarly, when other bases (A, T, and G) mutate to cytosine, the methylation type and methylation level were acquired, resulting in differences between the two sites. In addition, non-cytosine mutations also result in changes in methylation patterns and methylation levels, for example, CHG=>CHH and CHH=>CG. Studies in Arabidopsis[27] and cassava[26] have revealed that SNPs can change both methylation type and level. This indicated a close relationship between genetic variations and epigenetic during plant growth and development. However, current research on the co-regulatory mechanisms of genetic and epigenetic variations in plant phenotypes remains relatively limited.

      Honeysuckle (Lonicera japonica Thunb.), belongs to the genus Lonicera of the family Caprifoliaceae is named for its flower development process (Fig. 1a), in which the color changes from silver-white to golden-yellow. Honeysuckle is rich in medicinally active compounds, including luteoloside, chlorogenic acid, flavonoids, and sesquiterpenes. Pharmacopeia records indicate that the flowering stage, known as the white bud (WB) stage, has the highest content of medicinally active compounds[28]. Among these, flavonoids play a crucial role in antiviral activity[29] and have been used to treat various viruses. Sijihua and Beihua No. 1 (Beihua) are the two most common honeysuckle cultivars in China and show significant phenotypic differences. Sijihua has a short WB stage, lasting only 2−3 d, whereas Beihua can have a WB stage lasting up to 20 d, allowing for a longer storage period of high-content medicinal active compounds. Additionally, the flowers exhibit conspicuous color variations at distinct developmental stages. Research has demonstrated that Beihua consistently displays higher levels of total flavonoids and other essential medicinal compounds compared to Sijihua[30]. A recent study revealed that DNA methylation in Sijihua affects the expression of carotenoid-related genes, leading to variations in flower color during different stages of development[31]. Honeysuckle, a crucial antiviral medicinal plant source, exhibits significant phenotypic differences among cultivars, including variations in metabolite accumulation and other traits. Growth and development, particularly flowering, are regulated by DNA methylation. This makes the honeysuckle an ideal material for studying the co-regulatory mechanisms of genetic variations and epigenetics in flavonoid compound synthesis.

      The release of the honeysuckle genome has laid the foundation for studying the genomic structure of different cultivars[32]. Whole-genome resequencing technology to sequence honeysuckle will provide a comprehensive and high-resolution view of genetic variations among different cultivars. We aimed to elucidate the molecular mechanisms underlying phenotypic differences and accumulation of active compounds in different honeysuckle cultivars from the perspectives of genetic and epigenetic variations. In this study, multi-omics approaches were employed, including whole-genome resequencing, whole-genome bisulfite sequencing, and transcriptome sequencing. This comprehensive investigation focused on Beihua No. 1 (Beihua) and Sijihua during two key stages (WB and GF), characterized by differences in phenotypes and metabolite accumulation. Using whole-genome resequencing, a large number of genetic variations were identified a pseudo-reference genome was constructed for Beihua, single-base resolution DNA methylation profiles were compared between Beihua and Sijihua at different flowering stages, the characteristics of DNA methylation changes examined in different contexts. The relationship between SNPs and differentially methylated cytosines (DMCs) in honeysuckle genomes during two developmental stages were explored. Furthermore, the impact of SNPs on key enzyme-encoding genes involved in the flavonoid biosynthesis pathwaywere analyzed, thereby revealing crucial factors contributing to the differences in flavonoid biosynthesis and accumulation between the two cultivars. This study reveals the molecular mechanisms by which genetic variations and epigenetic co-regulation to phenotypic differences and variations in metabolite accumulation in the two honeysuckle cultivars. This study provides new insights into honeysuckle breeding and cultivation techniques and serves as a reference for exploring regulatory mechanisms of growth and development through the application of genetic variations and epigenetics in other species.

    Materials and methods

      The source of Beihua's materials and Sijihua's sequencing data

    • In this study, flower tissues in different stages of honeysuckle cultivar 'Beihua No.1' (abbreviated as 'Beihua' throughout the entire article) were collected, including juvenile bud (JB), green bud (GB), white bud (WB), silver flower (SF), and golden flower (GF), from Zhongke Honeysuckle Planting Cooperative in Pingyi County, Shandong, China (35°31'02'' N, 117°36'55'' E). Whole-genome bisulfite sequencing (WGBS) and transcriptome sequencing were performed using flower tissues from two stages of the honeysuckle cultivar Beihua, namely, the WB and GF stages. Young leaf tissue from Beihua was used for the resequencing.

      Three biological WGBS libraries were sequenced using the Hiseq X10 sequencer (Illumina, San Diego, CA, USA) as paired-end 150-bp reads. Three biological libraries were generated using the VHTS Universal V6 RNA-seq Library Prep Kit following the manufacturer's instructions (Illumina). All libraries were sequenced using the Novaseq 6000 platform (Illumina, San Diego, CA, USA), and 150 bp paired-end reads were generated.

      The reference genome for the cultivar 'Sijihua' of honeysuckle was obtained from the NCBI database under accession number PRJNA794868. Additionally, transcriptome data for the bud stage (WB: SRX14408207 − SRX14408209) and the flowering stage (GF: SRX144082013 − SRX144082015) were downloaded from the same database. Whole-genome bisulfite sequencing (WGBS) data for the WB stage (SRR18684863 − SRR18684865) and GF stage (SRR18684866 − SRR18684868) were also retrieved from NCBI. It is worth noting that the collection location and environmental conditions for the cultivation of the Sijihua cultivar of honeysuckle are identical to those of Beihua (Honeysuckle Planting Cooperative in Pingyi County, Shandong, China (35°31'02'' N, 117°36'55'' E)[31].

      Library construction and sequencing

    • For resequencing library construction, DNA isolation was fragmented using Bioruptor (ThermoFisher Scientific, Waltham, MA, USA), resulting in libraries with approximately 300 bp fragment sizes. Quality control of the libraries was performed using the Qubit dsDNA HS Assay Kit (ThermoFisher Scientific) and Agilent 2100 Bioanalyzer System (Agilent Technologies, Santa Clara, CA, USA). High-quality DNA libraries were then sequenced on the BGISEQ-500 platform, generating reads of 150 bp in length.

      WGBS libraries were constructed and prepared using the TruSeq DNA L T kit (Illumina) as described previously. Three biological WGBS libraries were sequenced on Illumina Hiseq X10 sequencers as paired-end 150-bp reads.

      Total RNA was extracted from the honeysuckle cultivar 'Beihua' petals at the WB and GF stages using the cetyltrimethylammonium bromide (CTAB) method. Three biological libraries were prepared using the VHTS Universal V6 RNA-seq Library Prep Kit following the manufacturer's instructions from Illumina. Subsequently, all libraries were sequenced using the Novaseq 6000 platform from Illumina, located in San Diego, CA, USA, producing 150 bp paired-end reads.

      Analysis of whole-genome resequencing data and construction of Beihua pseudogenome

    • We used FASQC (www.bioinformatics.babraham.ac.uk/projects/fastqc) for quality control of the generated FASTQ files, after which low-quality reads were filtered out, including those with more than 10% N content and more than 50% of low-quality bases (< 10). The net data obtained after filtering were compared to the Sijihua reference genome using BWA software (v0.7.12). Duplicates generated by PCR amplification were labeled and removed using the Picard package (https://sourceforge.net/projects/picard/). The 'HaplotypeCaller' function of GATK4 (version 4.1.4.1, https://hub.docker.com/r/broadinstitute/gatk/) was employed to generate GVCF files. Raw variant calling sets underwent hard filtering with parameters 'QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0'. The VCF file generated by hard filtering included the chromosome number of SNPs, SNP positions, reference bases, mutant bases, etc. Next, bedtools maskfasta was used to align the SNP information of the VCF file in the reference genome of Sijihua and finally constructed the pseudo-genome of Beihua. The downstream analysis script add_ka_ks.pl from MCScanX was used for calculating Ka/Ks (non-synonymous/synonymous) values for each gene pair.

      SNP annotation was performed using SNPEff software (https://sourceforge.net/projects/snpeff/) based on the constructed Beihua pseudogenome, SNPs were classified into intergenic regions, upstream, downstream, Splice-site-donor, Splice-site-acceptor, UTR-5-prime, UTR-3-prime, and SNPs encoding exons were further divided into synonymous and nonsynonymous mutations.

      Whole genome bisulfite sequencing data analysis

    • Raw bisulfite sequencing reads of Beihua was trimmed by Trimmomatic v0.39 to obtain clean data and aligned to the replaced Beihua pseudogenome by BSMAP v2.90. A 4% mismatch rate was allowed per 150 bp of read length and only uniquely mapped reads were retained for subsequent analysis. Next, reads aligned to unmethylated lambda DNA were used to estimate the conversion rate. Whether cytosine is methylated or not was identified mainly based on the conversion rate and binomial distribution. The methylation level of each cytosine was calculated using the methration.py script of the BSMAP software. The formula for calculating the level of cytosine methylation is #C/(#C+#T), with #C representing the methylated cytosine and #T being the unmethylated cytosine. Similarly, whole genome bisulfite sequencing data of Sijihua were calculated by the same method. DMCs between the two cultivars of Beihua and Sijihua were identified using methylkit default parameters. In addition, the absolute methylation levels of white and golden flowers of both cultivars, under the same contexts of CG, CHG, and CHH, were to vary more than 40%, 20%, and 10%, respectively. In this paper, both the loss of CG, CHG, and CHH (CG-loss, CHG-loss, and CHH-loss) and the gain of CG, CHG, and CHH (CG-gain, CHG-gain, and CHH-gain) were defined as DMCs, and this type of DMCs is relative to the two extremes of variation from the presence to the absence and from the absence to the presence of the two cultivars.

      Transcriptome sequencing data analysis

    • For each biological replicate of the Beihua and Sijihua transcriptome data, we used Trimmomatic v0.39 for quality control and aligned to the Beihua pseudogenome and Sijihua reference genome using HISAT2 default parameters, and we kept only uniquely mapped reads to estimate expression values (Supplemental Table S1). The expression of each gene we quantified using stringTie v2.1.7. Next, DESeq2 v1.32.0 was used to identify differentially expressed genes in white buds and golden flowers of both Beihua and Sijihua cultivars and differentially expressed genes were required to satisfy Log2|FC| > 2 and FDR < 0.05.

      GO enrichment analysis

    • GOATOOLS was used for gene ontology (GO) enrichment of overlapping genes of differentially expressed genes (DEGs) in white buds and golden flowers of both Beihua and Sijihua cultivars. Only GO terms with P-values less than 0.05 were retained for analysis. Heatmap, and GO enrichment plot of DEGs were made using R software version 3.5.

    Results

      Identification of genomic variations and pseudo-genome construction

    • To reveal the genomic variations between Sijihua and Beihua, Beihua was resequenced (Fig. 1a), generating 3.7 Gb of high-quality clean data. Reads were mapped to Sijihua's reference genome using BWA software. Over 91% of read pairs successfully aligned at a coverage depth of 30 × (Supplemental Table S2), facilitating subsequent genetic variation analysis. 9,909,981 SNPs were identified in the Sijihua genome. The majority of these SNPs were situated in intergenic regions (72.42%), followed by intronic regions (8.4%), upstream (8.35%), and downstream (7.34%) regions (Fig. 1b, Supplemental Table S3). SNPs were classified into four groups based on their impact on genes (https://pcingola.github.io/SnpEff/se_inputoutput/). 12,688 high-impact SNPs (0.13%) resulted in stop-gained, stop-loss, splice-site acceptor, splice-site donor, and start-lost mutations. 102,188 low-impact SNPs (1.03%) resulted in synonymous coding, start gained, synonymous stop, or non-synonymous start mutations. Additionally, 158,107 moderate-impact SNPs (1.6%) caused non-synonymous mutations in the coding region. The remaining 9,636,998 SNPs (97.25%) had a modifying impact on intergenic, intronic, upstream, downstream, UTR-3-prime, and UTR-5-prime regions. A total of 38,808 genes were affected by four SNP categories. These included 7,555 genes associated with high-impact SNPs, 24,967 genes associated with moderate-impact SNPs, 23,256 genes associated with low-impact SNPs, and 38,787 genes associated with modifier-type SNPs (Fig. 1c, Supplemental Table S4). High-impact SNPs may be involved in the expression of important genes with key functions[33]. Therefore, a Gene Ontology (GO) functional enrichment analysis of genes associated with high-impact SNPs were performed. Regardless of the white bud (WB) or golden flower (GF) stage, genes related to high-impact SNPs were significantly enriched in response to stimuli, response to viruses, gene silencing regulation, and RNA interference regulation (Fig. 1d).

      Previous studies have reported the construction of pseudo-reference genomes with incomplete genome assembly by using genetic variation substitution in phylogenetic closed cultivars, such as to obtain pseudo-reference genomes for maize Mo17 and rice Matsumae, researchers substituted SNP information for maize B73[34], and rice RZ35[35]. Since the genomes of Sijihua and Beihua are very similar (SNP ratio of 1.12%), a pseudo-reference genome was created for Beihua by replacing the SNP sites with the reference genome of Sijihua.

      Figure 1. 

      Identification and classification of SNPs. (a) Different development stages of Beihua and Sijihua. Five stages are juvenile bud (JB), green bud (GB), white bud (WB), silver flower (SF), and golden flower (GF), respectively. (b) Classification of SNPs based on genome location. (c) Number of genes affected by SNPs (SNPs were divided into four impact types: high, moderate, low, and modifier). (d) GO enrichment analysis of DEGs related to high-impact SNPs.

      Comparison of DNA methylation between Sijihua and Beihua

    • To investigate DNA methylation distinctions between the two honeysuckle cultivars, DNA methylation was sequenced at two stages of floral development (WB and GF) in Beihua using the WGBS technique (Supplemental Table S5). Compared with DNA methylation in Sijihua[31], no significant differences in global CG DNA methylation levels were found between the two cultivars at different developmental stages (T-test, WB: p-value = 0.17, GF: p-value = 0.62). we observed Beihua had significantly higher global DNA methylation levels than Sijihua in the CHG and CHH sequence contexts (T-test, p-value < 0.05). This phenomenon was consistently observed in both WB and GF stages (Fig. 2a). The methylation of WB and GF in Beihua and Sijihua were also analyzed. It was found that in Beihua, only the CHH methylation level showed significant differences, while in Sijihua, methylation levels in all three contexts exhibited significant differences (Supplemental Fig. S1). Furthermore, the average DNA methylation levels in 500 Kb windows across the genome were calculated and observed varying degrees of bias in all three sequence contexts. Specifically, the CG methylation distribution deviated to the right (Supplemental Fig. S2a & e), indicating that CG methylation levels were lower in Beihua than in Sijihua. The differential distribution of CHG and CHH methylation showed a significant left deviation (Supplemental Fig. S2b, f, c, & h), indicating that the methylation levels of CHG and CHH were significantly higher in Beihua than in Sijihua. These results were consistent with global DNA methylation levels (Fig. 2a).

      The number of methylated cytosines were further compared between Sijihua and Beihua, revealing that Sijihua had 138,186,587 methylated cytosines, whereas Beihua had 113,582,475. During the WB stage, Sijihua had a higher number of methylated cytosines in the CG (27,577,348), CHG (21,837,724), and CHH (88,771,515) contexts compared to Beihua, which had CG (22,491,550), CHG (17,537,206), and CHH (17,537,206) methylated cytosines. A similar trend was observed during the golden flower stage, indicating that Beihua's DNA structure is less susceptible to methylation than Sijihua (Fig. 2b). When calculating the methylation levels of each cytosine, both CG and CHG methylation in the two honeysuckle cultivars showed frequent occurrences of two extremes (0 and 1), representing the completely methylated and unmethylated states, respectively. Both showed a bimodal distribution pattern. In Sijihua, 72.9% and 36.6% of the CG and CHG methylation sites, respectively, occurred in a completely methylated state, whereas in Beihua, 69% and 33.7% of the CG and CHG methylation sites, respectively, occurred in a completely methylated state. The proportion of completely methylated CG and CHG in Sijihua was higher than that in Beihua, this trend remained consistent at the GF stage (Fig. 2c & d). Furthermore, we found that the higher percentage of completely methylated CG in Sijihua, approximately 80.7% (18,525,847 / 22,956,440) of CG methylation, could be maintained from WB to GF. In contrast, in Beihua, only 76.6% (13,962,522 / 18,227,835) of CG methylation was maintained from WB to GF. This suggests that during the process of maintaining these methylation types, CG methylation in Sijihua may be more stable and faithfully replicated during DNA replication than that in Beihua[36].

      To investigate methylation distribution in two honeysuckle genomes, gene density, transposable elements (TE) density, SNP density, and average methylation level were calculated in the three methylation contexts of Sijihua and Beihua (Fig. 2e). The TE-enriched regions of both Sijihua and Beihua showed higher DNA methylation levels and lower gene density at both WB and GF stages (Fig. 2e, Supplemental Fig. S2d). This phenomenon in many flowering plants, such as Arabidopsis[37], soybean[38], maize[39], rapeseed[40], and tomato[41]. SNPs showed an uneven distribution across the genome, with a high content occurring in regions with elevated CG and CHG methylation levels, while being scarce in CHH-methylated regions (Fig. 2e). To explore the relationship between DNA methylation and SNP density, the correlation between SNP density and methylation in three sequence contexts were analyzed. Significantly, we observed that SNP density was positively correlated with CG and CHG methylation levels, and negatively correlated with CHH methylation (Supplemental Fig. S3). This suggests that highly methylated cytosines may influence the stability of DNA sequences, especially the CG and CHG methylation.

      To examine the methylation patterns in the coding regions of the two honeysuckles, the methylation levels within the gene coding regions and their 2 Kb flanking regions for both Sijihua and Beihua were calculated. CG methylation levels in the coding regions were higher in Sijihua than in Beihua, whereas CHG and CHH methylation levels were lower in Sijihua than in Beihua. This pattern was consistent across both WB and GF (Supplemental Figs S4a−c, S5a-c). Based on previous reports indicating the abundance of transposable element insertions in the coding regions of honeysuckles[31], the genes were classified into two types: TE-related genes and TE-unrelated genes. Subsequently, the DNA methylation patterns of these two gene types were calculated. The methylation pattern of the coding regions of TE-related genes was the same as the genome-wide methylation pattern in both honeysuckle cultivars during the two stages, whereas TE-unrelated genes showed significantly reduced methylation levels. This suggests that TE insertions play a crucial role in maintaining the DNA methylation levels in both Sijihua and Beihua (Supplemental Figs. S4d−f, S5d−f). In addition, the DNA methylation of the TE region and the 2 Kb flanking region were calculated. The CG methylation levels of TE in Beihua was higher than that in Sijihua at two stages. However, the CHG methylation level in the TE regions did not differ significantly between these two cultivars at the WB stage, but showed marked differences at the GF stage, indicating dynamic methylation levels of TEs during honeysuckle flower development. Surprisingly, CHH methylation levels of TE regions in Beihua were significantly higher than in Sijihua at both stages (Supplemental Figs S4g−i, S5g−i). Overall, the methylation levels of Sijihua and Beihua differed significantly and varied with developmental stages.

      Figure 2. 

      DNA methylation landscape of Beihua and Sijihua. (a) Comparison of whole-genome DNA methylation levels in Beihua and Sijihua during the WB and the GF stages. Three biological replicates were calculated as dots over the bar graph (T-test, * p-value < 0.05, NS represents non-significant). (b) The relative proportions of methyl-cytosine in the contexts of CG, CHG, and CHH in Beihua and Sijihua during the WB and GF stages. (c)−(d) Distribution of methyl-cytosine methylation levels in all three sequence contexts in Beihua and Sijihua during the WB(C) and GF(D) stages were compared. (e) Circos plot showing gene density, TE density, SNP density, and DNA methylation levels for all three contexts of the WB stage.

      Differentially methylated cytosines are mainly associated with SNPs

    • By comparing the DNA methylation patterns between Sijihua and Beihua, significant differences in methylation levels of all three sequence contexts were observed. To gain a deeper understanding of the specific sites displaying differential methylation, the study focused on comparing methylation levels between Sijihua and Beihua. During the WB stage, the analysis identified 1,602,266 (19%) DMCs, including 641,008 CG DMCs, 589,493 CHG DMCs, and 371,765 CHH DMCs. Notably, hyper-methylation (Sijihua > Beihua) was more prevalent than hypomethylation in three methylation contexts. This trend persisted in the GF stage (Fig. 3a). However, it is important to consider that genetic variation (SNPs) can also contribute to the variation in DNA methylation sites by influencing the sequence context or altering methylation levels[26]. To further understand how genetic variations (SNPs) impact DNA methylation sites in honeysuckles, 6,679,481 (81%) SNP-associated DMCs were identified at both WB and GF stages. It was observed that 2,539,019 DMCs resulted from cytosine loss (36%), 1,808,116 DMCs were attributed to cytosine gain (26%), and 2,332,346 DMCs were ascribed to other types of mutations (38%). For example, an illustrative snapshot from the Genome Browser displayed various SNP-associated DMCs, including CHH-gain, CG-gain, CHG-loss, CG-loss, and CHH-loss (Fig. 3b). Genetic variations can lead to changes in both methylation type and level. Moreover, it was found that the number of SNP-associated DMCs was much greater than the number of non-SNP-associated DMCs (DMCs caused by different methylation levels in the same context) when comparing differentially methylated cytosines within the two honeysuckles (Fig. 3a). These findings underscore the potential substantial contribution of genetic variations in driving disparities in DNA methylation levels between the two honeysuckles.

      Regardless of the number of cytosine methylation type mutations or non-cytosine mutations, it was found that the CHH type was the most affected in terms of cytosine loss, gain, and other types of mutations, followed by the CHG and CG types (Fig. 3c & d). Among other types of mutations, the significance of guanine (G) loss and gain was particularly evident, playing a pivotal role in the dynamic transition between distinct methylation types. Transitions from CG to CHH and from CHG to CHH are the most prevalent, with a predominant proportion of guanine (G) mutations mutating to adenine (A) or thymine (T). Following this trend, mutations from CHH to CG and CHH to CHG emerged as the next most frequent occurrences, primarily involving the conversion of adenine (A) to guanine (G) (Fig. 3e). It is worth noting that these mutations not only led to changes in methylation type but also influenced the methylation levels themselves.

      Specially, when CG sites underwent mutations, transforming into CHG or CHH sites, they exhibited a higher likelihood of becoming hypomethylated. Conversely, CHH sites transitioning to CHG or CG sites were more prone to hypermethylation (Supplemental Fig. S6a). Extending the examination to the diversity of the 500 bp flanking sequences surrounding the SNP-associated DMC sites, it was observed that downstream nucleotide diversity was higher than that upstream of the DMC sites. Additionally, nucleotide diversity at the DMC sites markedly exceeded that within the adjacent flanking sequences. Intriguingly, this phenomenon was not observed at non-DMC sites. The higher nucleotide diversity at DMC sites underscores their susceptibility to frequent natural selection, implying the potential for intricate interplay between SNPs and DNA methylation dynamics (Supplemental Fig. S6b).

      Figure 3. 

      Identification of differentially methylated cytosines in the honeysuckle genome. (a) The proportion of two DMC types in the honeysuckle genome, including non-SNP-associated DMCs and SNP-associated DMCs. (b) Representative screenshots of DMCs. Methylation site losses, gains, and mC context changes were indicated on the underside of the track (red: Sijihua; blue: Beihua; Light blue shading represents non-SNP-associated DMCs; Light green represents SNP-associated DMCs). (c) The proportion of methylation types of the SNP-associated DMCs. (d) Sankey diagram showing the number of changes in DNA methylation types in three contexts between Sijihua and Beihua. (e) Sunburst chart showing the number of methylation type changes caused by nucleotide mutation.

      Mutations in CG=>TG affect the methylation level of genes

    • To further investigate the effect of cytosine loss on DNA methylation levels, the base substitution rates were calculated for homologous genes between Beihua and Sijihua. The rate of C=>T base substitutions was higher than that of other mutation types (Fig. 4a), consistent with studies in Arabidopsis[23] and lotus[42]. Within the honeysuckle genome, methylated cytosines exhibited a greater propensity to mutate into C=>T compared to unmethylated cytosines (1,155,283 for mC=>T, 636,736 for non-mC=>T, binomial test, p-value = 1.727e-10), resulting in a higher frequency of CG to TG mutations[24]. However, CG=>TG mutations were not evenly distributed throughout the genome. The majority of these mutations occurred in intergenic regions (272,395), followed by introns (19,132) and exons (44,634) (Fig. 4b).

      To validate whether the mutation from CG methylation type to TG in honeysuckle impacted gene methylation levels, we assessed CG methylation levels within the gene body and the 2 Kb flanking regions of both Sijihua and Beihua genes. The analysis revealed significant differences in the gene body regions of both Sijihua and Beihua at both stages (Mann-Whitney test, p-value < 0.001) (Fig. 4c, Supplemental Fig. S7a). When genes unaffected by CG=>TG mutations were considered, as expected, no significant differences (Mann-Whitney test, p-value > 0.05) were observed in the gene-body regions at two stages of Sijihua and Beihua (Fig. 4d, Supplemental Fig. S7b). Subsequently, the study was focused specifically on genes with CG=>TG mutations and calculated the CG methylation levels in their gene body and 2 Kb flanking regions[43]. Intriguingly, significant differences in the gene body methylation levels between Sijihua and Beihua were observed (Mann-Whitney test, p-value < 0.001) (Fig. 4e, Supplemental Fig. S7c). This indicates that genetic variations leading to CG=>TG mutations directly affect CG methylation levels in the gene body. Furthermore, it was observed that genes affected by CG=>TG mutations not only showed reduced CG methylation levels in the gene body but also displayed a significant reduction in gene expression levels in Beihua and Sijihua (Fig. 4f−h, Supplemental Fig. S7d−f). This observation is consistent with previous findings in other plants such as Arabidopsis[44], and Brassica napus[45], where higher CG methylation within the gene body promotes gene transcription.

      In summary, CG=>TG mutations not only altered the methylation levels, but also affected gene expression levels. To understand the functions of these genes affected by CG=>TG mutations, GO enrichment analysis was performed and identified their predominant enrichment in transcriptional regulation, stress response, flavonoid biosynthesis pathways, responses to viruses, biosynthesis of flavonoids, flavonols, and flavonoids (Supplemental Fig. S8). Collectively, these findings underscore the critical role played by the interactions between epigenetic and genetic variations in regulating gene functions in honeysuckle.

      Figure 4. 

      Relationship of CG=>TG type mutations with CG methylation and gene expression. (a) The substitution rate of nucleotides in homologous genes between Beihua and Sijihua. (b) The distribution of CG=>TG mutation sites on the honeysuckle genome. Comparison of CG methylation levels in the body region and the flanking 2 Kb region of (c) all genes, (d) CG=>TG unrelated genes, and (e) CG=>TG-related genes during the WB stage in Sijihua and Beihua. The expression levels of (f) all genes, (g) CG=>TG unrelated genes, and (h) CG=>TG-related genes during the WB stage in Sijihua and Beihua were compared (Mann-Whitney test, *** p-value < 0.001, NS represents non-significant).

      SNPs and DMCs co-regulate genes involved in the flavonoid pathway

    • To investigate deeper into the impact of SNP-associated DMCs, genes were identified that had overlaps with SNP-associated DMCs within their gene body and 2 Kb flanking regions, designating them as genes related to SNP-associated DMCs. It was found that SNP-associated DMCs were linked to 76% (29,899/39,320) of the genes within the honeysuckle genome. Among these genes, 5,324 and 7,067 genes showed differential expression (DEGs, log2|FC| > 2, FDR < 0.05) in the WB and GF stages, respectively, with an overlap of 3,325 genes displaying differential expression in both stages (Fig. 5a). Notably, the majority of genes (3,158/3,325) that exhibited differential expression at both stages showed consistent expression trends. Specifically, in both the WB and GF stages, the expression levels of these genes in Beihua were either higher than those in Sijihua (C3) or lower than those in Sijihua (C4) (Supplemental Fig. S9a). Interestingly, GO enrichment analysis of these genes with consistent expression trends revealed significant enrichment in processes related to flavonoid biosynthesis, flavonol biosynthesis, cellular glucan metabolism, stress response pathways, and floral whorl development (Supplemental Fig. S9b). Overall, genes affected by both SNPs and DMCs play important biological functions in various critical pathways, indicating their substantial regulatory significance.

      The above analysis has highlighted those genes affected by both SNPs and DMCs were significantly enriched in the flavonoid biosynthesis pathway. Flavonoids represent a vital group of secondary metabolites in plants renowned for their natural medicinal properties. They play crucial roles in plant growth, development, and defense against biotic and abiotic stresses. The synthesis of flavonoids originates from phenylalanine through the phenylpropanoid pathway, involving several key enzymes (Fig. 5d)[46]. In the investigation of the effects of SNPs and DMCs on flavonoid biosynthesis in honeysuckle, 35 homologous genes encoding key enzymes in the flavonoid biosynthesis pathway in the honeysuckle genome were identified (Fig. 5d, Supplemental Table S6). Despite the presence of a high number of SNPs in both the gene body and promoter regions of these 35 homologous genes, their Ka/Ks ratios are considerably lower than 1. This suggests that these genes have undergone purifying selection during evolution and possess relatively conserved structures and sequences. Consequently, genetic variation may not be the primary driving factor behind the differences in flavonoid biosynthesis between the two honeysuckles (Fig. 5b & e).

      Furthermore, by comparing the methylation levels of these 35 genes in Sijihua and Beihua at different stages, significant differences were observed in both the gene body and promoter regions. Additionally, some genes showed significant changes in expression. For example, LjPAL (EVM0007831), LjCHS (EVM0026111), LjCHI (EVM0013981), and LjFLS (EVM0027194) showed significantly higher expression levels in Beihua than in Sijihua at both the WB and GF stages (Fig. 5c, Supplemental Fig. S10). Interestingly, it was found that both the gene body and promoter regions of these 35 genes contained SNP-associated DMCs, as well as non-SNP-associated DMCs. These DMCs appear to play crucial roles in the regulation of gene expression (Supplemental Fig. S11). For example, a genome browser showed that LjCHI (EVM0013981) and LjFLS (EVM0027194) have a large number of SNP-associated DMCs and non-SNP-associated DMCs in the promoter region. Additionally, the expression levels of these genes were significantly higher in Beihua than in Sijihua (Fig. 5f & g, Supplemental Fig. S12). Therefore, the promoter regions of these key enzyme genes contain numerous nucleotide mutations that lead to changes in promoter region methylation, directly affecting gene expression. This may be a crucial factor contributing to the higher flavonoid levels in Beihua compared to Sijihua.

      Figure 5. 

      The relationship between genes related to SNP-associated DMCs and the flavonoid pathway. (a) Venn diagram showing the number of overlapping differentially expressed genes (DEGs) associated with SNP-associated DMCs between WB and GF. (b) Boxplot showing Ka/Ks ratios of 35 key genes and genomes in the flavonoid pathway of honeysuckle (Mann-Whitney test. *** p-value < 0.001). (c) Boxplot showing the methylation levels of CG, CHG, and CHH in the promoter regions of 35 key genes in the flavonoid pathway during the WB in Beihua and Sijihua (Mann-Whitney test. *** p-value < 0.001). (d) The pathway of flavonoid synthesis in plants. (e) Heatmaps showing the expression levels of 35 key genes of the flavonoid pathway, the number of SNPs in the upstream region of the genes, the number of SNP-associated with DMCs, and the number of non-SNP-associated with DMCs for Beihua and Sijihua at the WB and GF (* represents DEGs in the WB and GF stages). (f) Genome browser showing DNA methylation level and expression level of LjCHI (EVM0013981). (g) Genome browser showing representative screenshots of two types of DMCs in the promoter region (Chr01: 82302800 − 82302800) of LjCHI (EVM0013981).

    Discussion

    • In this study, whole-genome resequencing, whole-genome bisulfite sequencing, and transcriptome sequencing techniques were employed to sequence different honeysuckle cultivars. From both the genetic and epigenetic perspectives, the factors underlying the phenotypic differences observed between the two honeysuckles were investigated. Previous studies have revealed that high-impact SNPs in the genomes of various plants, such as rice[47], tomato[48], cucumber[49], cotton[50], and peach[51], have significant effects on important agronomic traits. These high-impact SNPs have been observed to play crucial roles in plant resistance. For example, in tomato, resistance genes carrying high-impact SNPs exhibit reduced susceptibility to Oidium neolycopersici infection, resulting in enhanced resistance. In cotton, the GhDRP gene with high-impact SNPs showed decreased expression upon infection with V. dahliae, resulting in leaf yellowing, shedding, and severe cell damage. This indicates a robust lignification response to counteract pathogen attacks. The present study identified a substantial number of genes in the honeysuckle genome affected by high-impact SNPs, significantly enriched in resistance pathways, such as response to fungi, defense response to virus, and response to mechanical stimulus. These findings provide a solid foundation for improving the desirable traits of honeysuckle cultivars.

      Subsequently, the genome-wide single-base resolution DNA methylomes between Beihua and Sijihua were compared. Notably, significant differences in DNA methylation were observed between these two honeysuckles in the three contexts (CG, CHG, and CHH). Some studies have suggested that changes in DNA methylation are influenced by genetic variations, subsequently regulating alterations in gene expression. In the present study, a correlation between SNP density and different methylation types was observed. Specifically, SNP density showed a significant positive correlation with CG and CHG methylation, but a negative correlation with CHH methylation. Previous reports have also indicated a positive correlation between differential methylation regions and the density of SNP variations, suggesting that genetic variations in proximity may influence some of the methylation differences observed between different cultivars[52]. Furthermore, a higher density of SNPs were also observed at sites with differentially methylated cytosines, a phenomenon previously observed in the cassava genome[26].

      In the genome of honeysuckle, there is a pronounced preference for nucleotide mutations (C=>T), primarily driven by cytosine methylation. While DNA methylation plays a critical role in regulating gene expression, it can also have detrimental effects because methylated cytosines are more susceptible to deamination, leading to the conversion of cytosine to thymine. Similar phenomena have been observed in the Arabidopsis[23], canola[24], and sugarcane[25]. The interplay between genetic variation and DNA methylation has a regulatory effect on gene expression[53]. The results revealed that genes related to Beihua affected by CG=>TG had significantly lower methylation levels in the gene body compared to Sijihua, and a significant difference in the gene expression between Beihua and Sijihua.

      Numerous studies have indicated that phenylalanine ammonia-lyase (PAL)[54] and chalcone isomerase (CHI)[55] served as the first and second rate-limiting enzymes in flavonoid synthesis, respectively. Overexpression of these enzymes leads to a significant accumulation of flavonoids. Additionally, the overexpression of flavanol synthase (FLS) results in a substantial accumulation of flavonoids[56]. Therefore, the promoter regions of these key enzyme genes contain numerous nucleotide mutations that lead to changes in promoter region methylation, directly affecting gene expression. This may be a crucial factor contributing to the higher flavonoid levels in Beihua compared to Sijihua. These findings align with prior research on the apple genome[57]. In summary, the interplay between SNP and DMCs in the biosynthesis of flavonoids in honeysuckle regulates the expression changes of genes related to key enzymes, thereby influencing the accumulation of flavonoids.

      A large number of DMCs in the genomes of two honeysuckles were identified, with most of them being SNP-associated DMCs. These DMCs were associated with approximately 80% of the genes and were significantly enriched in several vital pathways, such as the flavonoid biosynthesis pathway, flavonol biosynthesis pathway, anticancer pathway, and cellular glucan metabolic process. In particular, SNP-associated DMCs were found in genes related to key enzymes in the flavonoid synthesis pathway. Some of these genes directly influenced the accumulation of flavonoids and exhibited significant differential expression, including LjPAL (EVM0007831), LjCHI (EVM0013981), and LjFLS (EVM0027194). Previous studies have shown that the expression levels of LjPAL, LjCHI, and LjFLS are positively correlated with flavonoid accumulation in other plants, such as Arabidopsis[55], tomato[58], Dendrobium officinale[59], and Brassica napus[60]. In summary, the interplay between genetic variation and epigenetic regulation exerts a substantial influence on gene expression, leading to noticeable phenotypic differences between the two honeysuckles. These findings provide novel insights into breeding and cultivation techniques for honeysuckles.

    Conclusions

    • In this study, 9,909,981 SNPs were identified in Sijihua, including 12,688 high-impact SNPs that were significantly enriched in the stress resistance pathways. By comparing the DNA methylation patterns of Beihua and Sijihua, significant differences in DNA methylation levels were observed between the two honeysuckles. Thus, a substantial number of DMCs were identified between these two honeysuckles, with 81% of which were SNP-associated DMCs. Furthermore, methylated cytosines are more prone to mutation, resulting in CG=>TG, and altered DNA methylation, further regulating gene expression. SNP-associated DMCs are linked to 76% of protein-coding genes, with 3,325 genes exhibiting differential expression in both stages (WB and GF), and significantly enriched in the biosynthetic pathway of flavonoids. In the flavonoid pathway, important genes affecting flavonoid accumulation such as LjPAL, LjCHI, and LjFLS showed significant differences. The flanking 2 Kb region and body region of these genes produced a large number of SNP-associated DMCs, which is likely to be the reason that SNP-associated DMCs are regulating the overexpression of genes in Beihua, resulting in an increase in the accumulation of flavonoids in Beihua.

    Author contributions

    • The authors confirm contribution to the paper as follows: study conception and design: Wang H; re-sequencing performing and DNA methylation analysis: Yu X; transcriptome performing and biological pathway analysis: Yu X, Yu H; tissues collection and analysis: Lu Y, Zhang C; draft manuscript preparation: Yu X, Wang H. All authors reviewed the results and approved the final version of the manuscript.

    Data availability

    • The raw reads generated in this study have been deposited in the CNCB sequence read archive (SRA) with the accession number PRJCA018541. The reference genome, WGBS, and RNA-seq public data of Sijihua were downloaded from the NCBI database under the project numbers: PRJNA794868, PRJNA824715, and PRJNA813701, respectively.

      Acknowledgments

      • This work was supported by the National Natural Science Foundation of China (32160142), Guangxi Natural Science Foundation (2023GXNSFDA026034), Sugarcane Research Foundation of Guangxi University (2022GZA002), and State Key Laboratory for Conservation and Utilization of Subtropical Agro-bioresources (SKLCUSA-b202302) to Haifeng Wang.

      Conflict of interest

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

      • # Authors contributed equally: Xianyun Yu, Hang Yu

      Supplementary information
      • Supplemental Table S1
      • Supplemental Table S1 Summary of Sijihua and Beihua RNA-seq reads mapping.
      • Supplemental Table S2 Summary of whole genome resequencing mapping.
      • Supplemental Table S3 Classification of SNPs based on genome location.
      • Supplemental Table S4 Number of genes affected by SNP mutation.
      • Supplemental Table S5 Summary of Sijihua and Beihua BS-seq reads mapping.
      • Supplemental Table S6 Mining and analyzing flavonoid pathway genes of honeysuckle.
      • Supplemental Fig. S1 The significance analysis between Beihua WB and GF, as well as Sijihua WB and GF.
      • Supplemental Fig. S2 DNA methylation landscape of Beihua and Sijihua.
      • Supplemental Fig. S3 The correlation between SNPs and DNA methylation.
      • Supplemental Fig. S4 DNA methylation patterns of protein-coding genes and TEs in Sijihua and Beihua in the WB.
      • Supplemental Fig. S5 DNA methylation patterns of protein-coding genes and TEs in Sijihua and Beihua in the GF.
      • Supplemental Fig. S6 SNP lead to changes in the type and level of DNA methylation and characterization of the diversity of DMC sites.
      • Supplemental Fig. S7 CG≥TG-type mutations in relation to CG methylation and gene expression during the GF stage.
      • Supplemental Fig. S8 GO enrichment analysis of CG≥TG-related genes.
      • Supplemental Fig. S9 The pattern of differentially expressed genes related to SNP-associated DMCs.
      • Supplemental Fig. S10 The comparison of DNA methylation levels of 35 key genes of the flavonoid pathway in two stages of Sijihua and Beihua.
      • Supplemental Fig. S11 Heatmaps showing the expression levels of 35 key genes, the number of SNPs in the gene body of the genes, the number of SNP-associated DMCs, and the number of non-SNP-associated DMCs in Beihua and Sijihua during the WB and GF stages (* represents genes that showed significant differences in both stages).
      • Supplemental Fig. S12 Genome browser showing DNA methylation and expression levels of LjFLS (EVM0027194) and representative screenshots of two types of DMCs in the promoter region (Chr01:50581502-50581602).
      Rights and permissions
      • © 2024 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 (5)  References (60)
  • About this article
    Cite this article
    Yu X, Yu H, Lu Y, Zhang C, Wang H. 2024. Genetic and epigenetic variations underlying flavonoid divergence in Beihua and Sijihua honeysuckles. Epigenetics Insights 17: e002 doi: 10.48130/epi-0024-0002
    Yu X, Yu H, Lu Y, Zhang C, Wang H. 2024. Genetic and epigenetic variations underlying flavonoid divergence in Beihua and Sijihua honeysuckles. Epigenetics Insights 17: e002 doi: 10.48130/epi-0024-0002
    shu

Catalog

  • HTML

  • Introduction
  • Materials and methods

  • The source of Beihua's materials and Sijihua's sequencing data

  • Library construction and sequencing

  • Analysis of whole-genome resequencing data and construction of Beihua pseudogenome

  • Whole genome bisulfite sequencing data analysis

  • Transcriptome sequencing data analysis

  • GO enrichment analysis

  • Results

  • Identification of genomic variations and pseudo-genome construction

  • Comparison of DNA methylation between Sijihua and Beihua

  • Differentially methylated cytosines are mainly associated with SNPs

  • Mutations in CG=>TG affect the methylation level of genes

  • SNPs and DMCs co-regulate genes involved in the flavonoid pathway

  • 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