Epigenomic Regulatory Mechanism of Flowering in Apple Demonstrated in Two Varieties with Contrasted Flowering Behaviors

doi:10.21203/rs.3.rs-805966/v1

Download PDF

Research Article

Epigenomic Regulatory Mechanism of Flowering in Apple Demonstrated in Two Varieties with Contrasted Flowering Behaviors

https://doi.org/10.21203/rs.3.rs-805966/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Flowering is the necessary condition and yield basis for woody fruits in their life cycle. Although there has been considerable interest in the regulatory mechanisms underlying floral induction and flowering, the associated epigenetic modifications remain relatively uncharacterized.

Results: We identified the genome-wide of DNA methylation changes and the transcriptional responses in axillary bud of ‘Qinguan’ (QA) and ‘Fuji’ (FA) varieties with contrasted flowering behaviors. The DNA methylations were19.35%, 62.96% and 17.68% for FA, and 19.64%, 62.49% and 17.86% for QA in the CG, CHG and CHH contexts, respectively. Number of hypermethylated or hypomethylated DMRs in different regions were contributed to significantly up/downregulated gene expression. DNA methylation can positively or negatively regulate gene expression based on the CG, CHG and CHH contexts in different regions. Additionally, the huge differences in transcription of MIKC^c-Type MADS-box genes, and multiple flowering genes in multiple flowering pathways (i.e., light, age, GA and sugar) by changing DNA methylation, contributed to contrasted flowering behaviors in both QA and FA. Specifically, the floral meristem identify genes (i.e., FT, LEAFY, AP1 and SOC1) were significantly higher expression in QA than FA, but the floral repressor (i.e., SVP, AGL15, and AGL18) had an opposite result. Significant differences in multiple hormone levels were due to DEGs and their DMRs in their synthesis pathways, leading to both contrasted axillary bud outgrowth and flowering behaviors.

Conclusions: The whole-genome bisulfite sequencing (BS) libraries of QA and FA with diverse flowering capabilities have been constructed for finding whole-genome cytosine methylation profiles. The RNA sequencing of QA and FA and diverse flowering capabilities have been combined together to identify the gene expression patterns and the correlation with their methylation states so that we can better understand the epigenetic regulation mechanisms of floral induction and formation in apple.

Plant Physiology and Morphology

Plant Molecular Biology and Genetics

Apple trees

Epigenomic

DMRs

RNA sequencing

Flowering

DNA methylation

Apple is an economically important fruit tree worldwide, account for more than 50% of the global apple production area (FAO 2019, http://www.fao.org/faostat). However, apple industry still faced some key production problems with poor quality flower buds, and alternate bearing, especially, for the ‘Fuji’ varieties, in which were more than 65% of the total cultivated area in China [1–3]. Additionally, different varieties (i.e., profusely flowering ‘Qinguan’ and weakly flowering ‘Nagafu No. 2’) [4, 5], as well as different types of shoots leaded to diverse flowering capabilities, indicating that flower induction in woody fruit trees was a more complicated regulatory process, which is influenced by multiple factors [2, 3, 6]. Although, the studies on regulation mechanism and network of flowering in woody fruit trees had been reported [5, 7, 8], there is still little known regarding the molecular regulatory mechanisms of floral induction involved in epigenetic modification in apple trees.

Cytosine methylation (DNA methylation) as one of the regulatory mechanisms of epigenetic modification play key roles in the regulation of multiple life processes in plants, such as flower induction and flower development [3], juvenile-to-adult transition [9, 10], germination [11], reproductive development [12], and fruit development and ripen [13–15]. Additionally, DNA methylation can positively or negatively regulate gene expression based on which regions the methylation occurs involved in CG,CHG and CHH contexts in multiple woody plants, such as Apple [3], Pear [16], and Citrus [17, 18].

Now, the model plant, Arabidopsis thaliana, had evolved six major flowering regulation pathways [i.e., photoperiod, vernalization, autonomous, gibberellic acid (GA), thermosensory and aging] to regulate flowering under internal or a variety of environmental signal conditions [8, 19–21]. In this process, some flowering key genes including FT, SOC1, LEAFY, and AP1 as floral integration factor play key roles in regulation of flowering involved in mRNA transcription and epigenetic regulation [22–24]. For example, the expression of FT gene was regulated by epigenetic modification involved in ELF6 with H3K4 demethylases [25], as well as chromatin silencing complex of EMF1c during floral induction [26]. NF-Y transcription factor contributed to regulate flowering by H3K27me3 demethylation of the SOC1 gene [27]. Additionally, GA decreasing by GA-biosynthesis inhibitors contributed to affect global DNA methylation levels in GA pathways to affect floral induction [28].

The MADS-box family genes including MIKC^c, Mα, Mβ, Mγ, and Mδ subclades in plants as transcription factors known to bind to DNA had multiple biological functions, such as floral induction and development, seed development, and fruit ripen [29–31]. The new hypotheses for the genomic origins, biological conservation, and divergence of MADS-box gene family members were proposed associated with multiple functions, such as regulators of floral development (APETALA3 and PI), flowering time (FLC), and root development (AGL17) [31]. Multiple MADS-box genes, including AGL15 and AGL18 as floral repressors along with SVP play a negatively roles in floral induction by directly binding to the FT promoter [32]. SWP73B belong to SWITCH/SUCROSE NONFERMENTING chromatin remodeling complexes can bind to the promoter regions of multiple MADS-box genes, such as AGL24, SOC1, and AP1 to regulate flowering in Arabidopsis thaliana [33]. The SOC1-like genes, AGL42, AGL71 and AGL72 can promote flowering in the shoot apical and axillary meristems through a gibberellins-dependent pathway [22]. The decreasing of the H3K27 trimethylation in FLC gene by downregulation of the CUL4 activity level regulated flowering in the autonomous pathway [34]. The class II TCP transcription factors (TFs) integrated into the FT-FD complex to control floral initiation by directly binding to the promoter of APETALA1 (AP1) gene [35].

In this study, we investigated whether there are any epigenetic changes between axillary bud with profusely flowering ‘Qinguan’ (QA) and no flowering ‘Nagafu No. 2’ (FA) apple varieties, and whether the DNA methylation changes contribute to the expression patterns of genes involved in apple floral induction and development. We constructed whole-genome bisulfite sequencing (BS) libraries of QA and FA with diverse flowering capabilities to find whole-genome cytosine methylation profiles. Additionally, we combined RNA sequencing of QA and FA with diverse flowering capabilities to identify the gene expression patterns, as well as correlation with their methylation states, in which was contributed a better understanding of epigenetic regulation mechanisms of floral induction and formation in apple.

Phenotypes characteristics of ‘Fuji’ and ‘Qinguan’ axillary buds and their flowering ability

‘Fuji’ axillary buds were not flowering in the next year (April 15th, 2019), but ‘Qinguan’ axillary buds can flowering in the same time with ~ 86% flowering rate (Fig. 1A, B). Additionally, there were big differences of morphological characteristics between QA and FA (Fig. 1C), in which the fresh weight, length and width of ‘Fuji’ were significantly larger than ‘Qinguan’ axillary buds (Fig. 1D). Detailed micrographs of longitudinal structural sections of QA and FA in Fig. S1.

Genome-wide examination of cytosine methylation in ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud

About 118.36 and 120.35 million raw reads of two BS libraries for FA and QA were generated by the Illumina HiSeq 2000 platform, respectively (yielding ~ 35.51 and ~ 36.11 data for FA and QA) (Table S1). Additionally, 114,420,604 and 118,150,703 clean reads for FA and QA were produced, respectively, as well as 62,656,722 (54.76%) for FA and 73,170,730 (61.93%) for QA clean reads were uniquely mapped to the Malus x domestica reference genome (https://www.rosaceae.org/species/malus/malus_x_domestica/genome_GDDH13_v1.1) [36] (Table S2). The detail information of gene coverage for FA (76.4%, 63.12% and 20.42%) and QA (77.11%, 63.62% and 18.78%) involved in CG, CHG and CHH contexts, respectively, were indicated in Table S3.

The circus plots showed that the whole chromosome-wide distribution of DNA methylation involved in the density of 5mC in the CG, CHG, and CHH contexts, the density of TEs, the genes and RNAs density of each chromosome for QA and FA (Fig. 2A, B; Fig. S2), indicating that there were significant differences for these indicators among various chromosomal positions. Additionally, there were similar proportion of methylation levels in CG, CHG and CHH contexts for FA (i.e., 19.35%, 62.96% and 17.68% for CG, CHG and CHH, respectively) and QA (19.64%, 62.49% and 17.86% for CG, CHG and CHH, respectively) (Fig. 2C). The detail information of methylated CX percentage for FA and QA involved in CG, CHH and CHG contexts in each chromosome were showed in Fig. 2D, indicating that the percentage of CG context was significantly higher than CHG and CHH for both FA and QA. However, the mean methylation levels involved in CG, CHG and CHH contexts were significantly higher in QA than in FA (Fig. 2E). Bisulfite sequencing data for the mean depth and proportion of covered bases in each chromosome for the QA and FA were in Fig. S3. The methylation densities, and levels and in the all contexts in QA and FA were in Fig. S4.

Distribution of DNA methylation in genic regions

A heat map analysis showed that DNA methylation levels were significantly differences among different gene features including promoter, exon, intron, and repeat regions (Fig. S5A). The DNA methylation levels including CG, CHG and CHH contexts for both FA and QA were highest in the repeat regions, followed by the promoter, exon, and intron regions (Fig. S5A, B). Additionally, the DNA methylation levels in all contexts were higher in the regions 2-kb upstream and downstream of the gene-bodies than in the gene-bodies (Fig. S5C).

Circos plots showed that there were significant differences in CG, CHG and CHH contexts in the whole chromosome-wide distribution of DNA methylation levels between FA and QA (Fig. 3A). DNA methylation levels including the CG, CHG, and CHH contexts were similar between FA and QA in all the different gene regions (i.e., promoter, exon, intron, and other regions) (Fig. 3B), as well as the genebody and the regions 2-kb upstream and downstream of the gene-bodies (Fig. 3B). Specifically, the DNA methylation levels in exon region for CG and CHG contexts were significantly higher in QA than in FA (Fig. 3B), but the DNA methylation levels in promoter region for CHH, and in intron region for CHG and CHH contexts, as well as in repeat region for CHH context had an opposite results (Fig. 3B). Additionally, for CHH context in the genebody and the regions 2-kb upstream and downstream of the gene-bodies regions, the DNA methylation levels in QA were significantly lower than in FA (Fig. 3C), but for CG context in genebody region had an opposite result between QA and FA (Fig. 3C).

Identification of differentially methylated regions (DMRs) between axillary bud with profusely flowering ‘Qinguan’ and no flowering ‘Nagafu No. 2’, and functional annotation analysis

Circos plots showed that the whole chromosome-wide distribution of differentially methylated regions (DMRs) (FDR < 0.05) for both hypermethylated and hypomethylated in CG, CHG and CHH contexts between FA and QA (Fig. 4A; Fig. S6). Venn diagrams present unique or public DMRs in genes (QA_vs_FA_DMRs_genes) and in promoter regions (QA_vs_FA_DMRs_promoter_genes) between FA and QA in the CG, CHG, and CHH contexts (Fig. 4B). Specifically, the QA_vs_FA comparison group comprised 2820, 1426, and 1808 DMR_genes (Fig. 4B) as well as 2853, 2462, and 2547 DMR_promoter_genes in the CG, CHH, and CHG contexts, respectively (Fig. 4B). Additionally, the DMR methylation levels in the CG and CHG contexts of the DMRs were significantly higher in QA than in FA, but the DMRs in CHH contexts had an opposite result between FA and QA (Fig. 4C). A heat map also showed the visual information of the DMRs and their DNA methylation levels in all three contexts (FDR < 0.05) between FA and QA (Fig. 4D). The numbers of the hypermethylated or hypomethylated DMRs in all contexts also had a big difference among different genic regions (i.e., the promoter, TSS, TTS, exon, repeat, and other regions) (Fig. 4E). For example, the numbers of DMRs (i.e., hypermethylated or hypomethylated) in TSS, TTS, utr5, and utr3 regions were significantly less than other regions including promoter, exon, repeat (Fig. 4E). Additionally, the numbers of hypermethylated DMRs in CG and CHG contexts in repeat region were significantly more than hypomethylated DMRs between FA and QA, but for CHH context, there was an opposite result (Fig. 4E).

To know the function of the DMR_genes and DMR_promoter_genes, the most enriched GO terms were identified for both hypermethylated and hypomethylated in QA compared with FA involved in CG, CHG and CHH contexts (Fig. S7, S8). For example, the most enriched GO terms included metabolic process, cellular process, biological regulation for QA_vs_FA_hypermethylated in the CG context (Fig. S7). For the DMR_promoter_genes, metabolic process, cellular process, biological regulation in biological process, cell, cell part, and organelle in cellular component, binding, catalytic activity, transporter activity in molecular function were the most enriched GO terms for QA_vs_FA_hypermethylated in the CG context (Fig. S8). The KEGG pathway also revealed some significantly enriched pathways involved in the DMR_genes and DMR_promoter_genes for the QA_vs_FA comparison group (Fig. S9, S10), including the following: metabolic pathways, biosynthesis of secondary metabolites for hyper/hypo_DMR_genes in the CG context (Fig. S9A, D). Additionally, the information regarding the significantly enriched pathways for the Hyper/Hypo_DMR_promoter_genes in the CG, CHG, and CHH contexts was provided in Fig. S10A–F.

Transcriptome profiling and functional analysis of differentially expressed genes between axillary bud with profusely flowering ‘Qinguan’ and no flowering ‘Nagafu No. 2’ apple

The hierarchical clustering heat maps showed the differentially expressed genes in the QA_vs_FA comparison group including 19,656 (9,543 downregulated and 10,113 upregulated) DEGs (log₂ fold-change > 0; padj < 0.05) (Fig. 5A, B). The GO functional analysis revealed the most enriched terms for the DEGs involved in biological process (i.e., oxoacid metabolic process, organic acid metabolic process, and catabolic process for upregulated genes; as well as signaling, signaling transduction, and single organism for downregulated genes), cellular component (i.e., membrane protein complex, thylakoid, thylakoid part for upregulated genes; as well as cell periphery, cell cortex, plasma membrane part for downregulated genes) and molecular function (i.e., peptidase activity, endopeptidase activity, transferase activity for upregulated genes; as well as calcium ion binding, ubiquitin-protein transferase activity, calmodulin binding for downregulated genes) in QA than in FA (Fig. 5C). Additionally, plant-pathogen interaction, spliceosome, amino sugar and nucleotide sugar metabolism for the upregulated genes, as well as carbon metabolism, biosynthesis of amino acids, protein processing in endoplasmic reticulum for the downregulated genes in QA than in FA were the most enriched KEGG pathways (Fig. 5D). Additionally, the gene expression distribution, as well as pearson correlation among samples of QA and FA were in Fig.S11. Details regarding the RNA-sequencing data for both QA and FA were in Table S4, S5.

Association between DNA methylation and transcript abundance

The gene expression levels were divided into the following four quartiles based on their FPKM values: non-expressed (FPKM < 1), low expression level (1 ≤ FPKM < FPKM_25%), medium expression level (FPKM_25% ≤ FPKM < FPKM_75%), and high expression level (FPKM ≥ FPKM_75%) (Fig. 6A). We found that the lowly expression genes had high methylation levels, as well as highly expression genes had low methylation levels in the gene-body for CG, CHG and CHH contexts, and upstream of the TSS and downstream of the TTS regions for CG and CHG contexts for both QA and FA (Fig. 6A). However, for CHH context, the most highly expression genes had the high methylation levels in upstream of the TSS and downstream of the TTS regions (promoter regions) in both QA and FA (Fig. 6A).

To know the correlation between the gene methylation and expression levels, the methylated genes were divided into the following five quintiles based on the promoter and gene-body methylation levels: group 1 (methyl_level < level_20%), group 2 (level_20% ≤ methyl_level < level_40%), group 3 (level_40% ≤ methyl_level < level_60%), group 4 (level_60% ≤ methyl_level < level_80%), and group 5 (methyl_level ≥ level_80%) (Fig. 6B). We also found that methylation levels had a negative regulatory relationship with gene expression levels involved in the CG, CHG, and CHH contexts in the gene-body region, and CG and CHG contexts in upstream of the TSS and downstream of the TTS regions (promoter regions), but they had a positive regulatory relationship when the methylation in CHH context happened in promoter regions (upstream of the TSS and downstream of the TTS regions) for both QA and FA (Fig. 6B). Additionally, the detail information about the relationship between methylation difference and gene expression (FPKM) involved in CG, CHG and CHH contexts for both hypermethylated and hypomethylated DMRs in different regions (i.e. gene-body and promoter) in QA compared with FA were presented in Fig. 7A, B, C.

Huge differences in transcription of MIKCc-Type MADS-box genes by changing DNA methylation, contributed to contrasted flowering behaviors

We identified 132 MADS-box family genes in apple that belong to five subgroups (i.e., MIKC^c, Mβ, Mα and Mδ) (Fig. 8A; Fig. S12). Among them, these differentially expressed MADS-box family genes were mainly belonged to MIKC^c subgroup between QA and FA, but no differentially expressed MADS-box genes were found in the other two subgroups (Mα and Mδ) (Fig. 8A). Eleven differently expressed MADS-box genes belong to cluster I (i.e., AGL16_MD08G1206400; AGL18_MD00G1106300; SEP1_MD06G1204100/MD09G1009100; SEP4_MD13G1059300) were significantly lower expression in axillary bud with profusely flowering ‘Qinguan’ than ‘Nagafu No. 2’, but other 17 differently expressed MADS-box genes belong to cluster II (i.e., six AGL24, three AGL80, AGL6, AGL79 and AGL24) had an opposite result (Fig. 8B; Additional file 1). Specially, flowering key genes (SOC1_MD07G1123600 and AP1_MD06G1204400), as well as floral meristem identify genes (AG_MD10G1271000 and AP3_MD02G1136500) were belong to cluster II that showed higher expression levels in axillary bud with profusely flowering ‘Qinguan’. However, the key negative regulation gene in floral induction (SVP_MD01G1038600/MD15G1384600) belong to cluster I were significantly lower expression in axillary bud with profusely flowering ‘Qinguan’ than ‘Nagafu No. 2’ (Fig. 8B; Additional file 1). Additionally, the methylation levels of thirteen MADS-box family genes (i.e., AGL21_intron_CHG, AGL21_promoter_CHH and AGL12_intron_CHH and AGL62_exon_CHG) in different gene regions were significantly higher expression in axillary bud with profusely flowering ‘Qinguan’ than ‘Nagafu No. 2’ (Fig. 8C), which were contributed to down-regulation of these genes in axillary bud of ‘Qinguan’ than in ‘Fuji’. However, other eleven DMRs (i.e., AGL30_exon_CG, AGL42_promoter_CG, and AGL42_promoter_CHG) had an opposite methylation pattern between axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’ apple, which may contribute to the regulation of these MADS-box genes expression involved in promoting flowering in ‘Qinguan’ axillary bud (Fig. 8C). Similarly, IGV snapshots presented the representative DMRs that were hypermethylation or demethylation between the axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’, including AGL80_MD00G1098000_CHH, AGL42_MD01G1103300_CG, and AGL42_MD00g1004500_CG in the promoter, and AGL79_MD17G1001300_CHH in the intron, and AGL30_MD13G1070300_CG in the exon (Fig. 8D). Additionally, multiple differently expression of MADS-box genes formed complex gene interaction regulatory networks including floral promoter (i.e., SOC1 and AGL42), floral organ identify (i.e., AG and AP3), and floral repressor (i.e., SVP and AGL18), contribute to regulation of flowering in axillary buds of ‘Qinguan’ (Fig. 9).

Differentially expression genes in multiple flowering pathways between axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’ apple by changing DNA methylation, contributed to contrasted flowering behaviors

We identified these different expression genes (DEGs) and different methylation regions (DMRs) in multiple flowering regulating pathway between axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’ apple (Fig. 10; Additional file 2, 3). Firstly, these genes involved in key floral meristem identify (i.e., LFY, AP1, SOC1 and FT) were significantly higher expression in axillary buds of ‘Qinguan’ than in ‘Nagafu No. 2’ (Fig. 10A). In age pathway, eleven SPL genes that can directly up-regulated LFY and SOC1 showed significantly higher expression in ‘Qinguan’ axillary buds than in ‘Nagafu No. 2’, but there had an opposite result in two AP2-like and one TOE genes involved in negative regulation of SOC1 expression (Fig. 10A). In photoperiodic pathway that centered on the CO gene, seven COL and eleven BBX genes were significantly lower expression in ‘Qinguan’ axillary buds than in ‘Nagafu No. 2’, and the upstream regulator of the CO gene, such as five COPs, three SPAs, PHYA, CRY1 and CRY2 had a similar result (Fig. 10A). And another upstream regulator of CO gene involved in circadian clock pathway including FKF1, two CDF2 and two CDF3, three PRR and TOC1 also had a similar result (Fig. 10A). However, two NF-YA, four NF-YB and four NF-YA that can interact with CO gene were significantly higher expression in axillary buds of ‘Qinguan’ than in ‘Nagafu No. 2’ (Fig. 10A). In vernalization and autonomous pathway, these genes were significantly upregulated in ‘Nagafu No. 2’ compared with ‘Qinguan’ including two MSI1, SWN and VIN3, as well as two FLK and FY that can directly inhibit the expression of FLC gene to promoting flowering (Fig. 10A). Additionally, FAF1-like complex including ELF5, EFS, ATX2, three SDG2, two VIP4 and VIP5, and SWR1/SRCAP-like complex including ARP5, ARP6, SEF and PIE1 that can directly promote the expression of FLC gene to inhibit flowering were significantly higher expression in ‘Nagafu No. 2’ than in ‘Qinguan’, that contributed to promote floral induction in ‘Qinguan’ (Fig. 10A). And these genes including REF and SRP1 involved in histone modification were significantly differentially expression between axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’ (Fig. 10A), suggesting they play key role in regulation of flowering process in ‘Qinguan’. In sugar pathway, two types of genes (i.e., IDD4/5 and SUS4) belong to the upstream of FT and SPL genes were significantly higher expression in axillary buds of ‘Qinguan’ than in ‘Nagafu No. 2’ (Fig. 10A). And five TPS14 genes showed higher expression in QA than in FA, but one TPS2 and TPS4, as well as five TPS21 genes had an opposite result (Fig. 10A), indicating that TPS genes play multifaceted functions on flowering in sugar pathway. In addition, two RGA2 genes in GA pathway were significant differential expression between axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’, contribute to both different flowering ability (Fig. 10A).

We also identified many different methylation regions (DMRs) between axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’ involved in multiple flowering pathways (Fig. 10B). For example, five DMRs(i.e., LFY_MD14G1146700_CG in promoter, and ELF4_MD11G1074300_CHH in exon)belong to cluster I were significantly lower expression in axillary bud with profusely flowering ‘Qinguan’ than ‘Nagafu No. 2’ in flowering pathways (Fig. 10B), but sixteen DMRs (i.e., SPL1_MD16G1137700_CHH, and SLP14_MD14G1060299_CHH in promoter; TFL1_MD12G1023900_CG in promoter and TFL2_MD07g1081800_CG in exon) belong to cluster II were had an opposite result (Fig. 10B). In photoperiodic pathway, twenty-one DMRs belong to cluster II, such as COL4_MD16G1148500_CG, and CRY2_MD16G1217800_CHG, DOF_MD15G126570_CG/CHG in exon, BBX32_MD17G1147000_CG/CHG in exon showed significantly higher methylation levels in QA than in FA (Fig. 10B), but BBX32_MD17G1147000_CHH in promoter, and COP9_MD09G1015400_CG in exon belong to cluster I had an opposite result, suggesting that these DMRs in different regions may contribute to regulate genes expression levels. Similarly, six DMRs (i.e., SPA3_MD11G1247500_CG in exon and SDG2_MD17G1091000_CHG in promoter) belong to cluster I were showed significantly lower levels in QA than in FA, and nineteen DMRs (i.e., ARP4_MD03G1252000_CG in exon, and VIN3_MD08G1199400_CG in promoter ) belong to cluster II had an opposite result (Fig. 2B). In sugar and GA pathways, five TPS21_MD03G1214900_CHG in intron, GRA_MD13G1022100_CG/CHG in exon, and _CHH in promoter, and TPS21_MD03G1207000_CHH in exon belong to cluster I were significantly higher levels in FA than in QA, but these DMRs belong to cluster II (i.e., IDD5_MD16G1009600_CHH in promoter, and GAI_MDMD17G1260700_CHG in exon, and _CG in promoter) had an opposite result (Fig. 10B).

Additionally, IGV snapshots showed these were multiple DMRs in gene different regions between QA and FA that affecting gene expression, such as hypermethylation in QA than in FA (i.e., SPL14_CHH, CRY1_CHH in promoter, PHYB_CG, and TPS21_CHG in exon, as well as BBX32_MD17G1147000 involved in CG, CHH and CHG contexts in exon region) or demethylation in QA than in FA, such as NF-YB_MD00G1101600 involved in CG, CHH and CHG contexts in exon region (Fig. 10C), suggesting that gene expression were regulated by these DMRs in different regions. We also found that these DMRs in different regions involved in CG, CHH and CHG contests were negatively correlated with gene expression levels, such as BBX32_CG _exon, COL6_CG_exon, NF-YB_CG_exon, and COL14_CG_promoter (Fig. 10D), but DMRs in promoter involved in CHH context can promote gene expression, such as SPL14_MD14G1020600_CHH in promoter (Fig. 10D), indicating that DNA methylation can positively or negatively regulate gene expression.

Significant differences in multiple hormone levels were due to DEGs and their DMRs in their synthesis pathways, leading to both contrasted axillary bud outgrowth and flowering behaviors.

There were significant differences in multiple plant hormones between axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’ (Fig. 11A; Additional file 4). The ABA and GA contents were significantly higher in QA than in FA, as well as no significant difference in IAA, but content of CK had an opposite result between them (Fig. 11A), associated with the significantly higher rates of germination and flowering in QA than in FA (Fig. 11A).

In addition, these genes involved in ABA biosynthesis, such as ZEP, three NCED1, two NCED4, two ABA2 were significantly higher in axillary buds of ‘Nagafu No. 2’ than in ‘Qinguan’, a similar result about these key GA biosynthesis genes (i.e., two CPS1 and KS) can be found in our data (Fig. 11A), but an opposite result involved in these genes (i.e., GA2OX1, GA2OX2, and GA20X8) that changing active GAs into inactive GAs (Fig. 11A), contributing to lower GA levels in QA than in FA. We also found that IPT2 gene involved in CK biosynthesis, and YUCC8 involved in IAA biosynthesis were significantly higher expression in QA that in FA, that may lead to higher level of CK in QA than in FA (Fig. 11A). Interesting, these genes involved in SL biosynthesis (i.e., D27, MAX3, MAX4, and MAX1) were no significantly difference between QA and FA (Fig. 11A), but the genes belong to SL signaling including three D14 and two MAX2 were significantly higher expression in FA than in QA (Fig. 11A) involved in regulation of axillary bud outgrowth and formation through targeted regulation of BRC1 as well as it is regulated by auxin. Additionally, multiple PIN genes (i.e., two PIN1, PIN3, and five PIN-like) involved in IAA transport from axillary bud into other organs were significantly higher in QA than in FA (Fig. 11A), that making relative lower IAA level in axillary bud of ‘Qinguan’ to induce axillary bud outgrowth compared with ‘Fuji’.

Additionally, IGV snapshots showed GA2OX8_MD05G1341000 in CG, CHG and CHH in exon region involved in GA biosynthesis were demethylation in QA than in FA (Fig. 11B) associated with relative higher expression of GA2OX8 in QA, indicating that negative correlation between DNA methylation and gene transcription level. Additionally, the PIN1_MD12G1095100 in CHH in promoter region were also significantly higher in FA than in QA, but expression level of PIN1 had an opposite result between FA and QA (Fig. 11B), suggesting that DNA methylation play key roles in gene expression.

Identification of related differentially expressed genes by a qRT-PCR assay

These representative flowering related genes (i.e., FT_MD12G1262000, SOC1_MD07G1123600, AP1_MD06G1204400, SPL4_MD17G1236000, PHYA_MD04G1183900) were significantly higher expression including qRT-PCR and RNA-seq data in QA than in FA, but COL2_MD09G1079200 had an opposite result (Fig. 12A), additionally, a linear relationship between the qRT-PCR and RNA-seq data was detected for these genes between QA and FA, with an R² of 0.8943 (P < 0.01) (Fig. 12B).

Flower bud formation was the necessary condition and yield basis for fruits of woody fruit trees [4, 7]. Numerous studies have reported the regulation mechanism of flowering process in their whole life cycle of plants, including floral induction [37], juvenile-to-adult transition [38], reproduction [12], and floral organ identity process [39] involved in transcriptional and post-transcriptional, and epigenetic regulation in multiple flowering pathways. DNA methylation as epigenetic modifications regulation is critical for floral induction and formation in woody plants, such as Apple [3], Pear [16], and Citrus [18]. Now, genome-wide identification and distribution of DNA methylation and epigenetic variations in all contexts in plants (i.e., A. thaliana [23], rice [40], and woody trees [3] ), have confirmed their key regulatory roles in plant development and floral induction. In the current study, we compared the DNA methylation between axillary bud with profusely flowering ‘Qinguan’ and no flowering ‘Nagafu No. 2’ to elucidate the epigenetic changes associated with floral induction and formation in apple trees (Fig. 2). The proportion of methylation levels for FA with 19.35%, 62.96% and 17.68% for CG, CHG and CHH, as well as for QA with 19.64%, 62.49% and 17.86% for CG, CHG and CHH, in apple axillary bud (Fig. 1A, B; 2C), respectively, with similar results from other studies [3], indicating their correlation with diverse flowering capabilities between QA and FA.

Identification of differentially methylated regions (DMRs) between axillary bud with profusely flowering ‘Qinguan’ and no flowering ‘Nagafu No. 2’, and analysis of their functions (Fig. 4, 5), given us a chance to uncover their key roles in multiple biological processes, such as floral induction and flowering. Similar results were reported for other studies regarding chromatin modifications (DNA methylation) that contributed to regulate plant growth and reproductive development [3, 12]. Additionally, these DMRs between QA and FA were distributed in different gene regions with great differences, which was also an important factor contributing to the regulation of related gene transcription (Fig. 3). Indeed, the DMRs happened on various genomic regions (i.e., promoter, gene-body, or intron) may lead to different gene transcript levels in multiple plants [3, 37], indicating that gene transcription is closely related to methylation. In our data, we found that methylation levels had a negative regulatory relationship with gene expression levels involved in the CG, CHG, and CHH contexts in the gene-body region (Fig. 6A), but there had a positive regulatory relationship when the methylation happened in promoter regions involved in CHH context (Fig. 6B), suggesting that DNA methylation may play a negative or positive roles in regulation of gene expression based on different gene regions, and the CG, CHG and CHH contexts. Additionally, similar results have been reported in other aspects of plants [3].

Whether plants flowering or not, ultimately depends on converge on the promoting expression of key flowering integration factors (i.e., FT, SOC1, SPLs, AP1 and LAFY ) involved in multiple flowering pathways under environmental and internal cues [24, 39, 41, 42]. For example, multi-layered regulation of SPL15 and cooperation with SOC1 Integrated gibberellin (GA) and age flowering pathways to regulate flowering at the SAM [43]. In our data, the key flowering integration factors (FT, SOC1, SPLs, AP1 and LEAFY) were significantly higher expression in axillary bud with profusely flowering ‘Qinguan’ than in no flowering ‘Nagafu No. 2’ associated with DNA methylation in different gene regions (Fig. 10). Additionally, studies showed that the expression levels of these key flowering genes were significantly affected by DNA methylation and epigenetic variations in different gene regions [3, 18, 23]. For example, targeted DNA methylation represses the Block C enhancers of FT gene (located 5-kilobases (kb) upstream of the transcriptional start site) contribute to regulating the floral transition in the long day photoperiod [23]. SIN3-LIKE genes mediate long-day induction of flowering through histone deacetylation to directly repress the expression of three FT repressors: AP2 family transcription factors [44]. Indeed, the upregulated of FT genes, but downregulated of AP2 genes in QA compared with FA (Fig. 12) were affected by their DNA methylation in different gene regions, contribute to ‘Qinguan’ axillary bud flowering (Fig. 1) through epigenetic modification. Chromatin-dependent repression of floral integrator genes (SOC1 and FT) by associated with plant specific PHD-containing proteins, played a key role in flowering [45]. Study showed that LFY and AP1 were conserved floral regulators in plants [39], as well as the floral transition was largely controlled by the meristem identity regulator LAFY gene that activating APETALA1 (AP1) gene [46]. In a word, flowering in axillary bud of ‘Qinguan’ was largely determined by the expression of the key flowering genes (i.e., FT, SOC1, LAFY and AP1), which were significantly regulated by transcription and epigenetic regulation in multiple pathways.

In our data, the expression of MIKC^c subgroups MADS-box family genes had significantly differences between in QA and FA (Fig. 8, 9), suggesting that they had a key role in flowering associated with epigenetic variations in different regions, similar results had been reported in other study [47]. For example, the different MADS-box family genes including floral regulators FLC and SOC1 that can directly regulate SVP gene expression play a key role in floral transition by antagonistic chromatin modifications [32]. The NF-Y interacts with CONSTANS directly regulated the SOC1 expression by modulating trimethylated H3K27 levels on promoter involved in epigenetic regulation [27]. Our data showed that multiple NF-YA/B/C family genes, as well as SOC1 were significantly higher levels in QA than in FA, accompanied by DNA methylation differences in different regions, contribute to relatively higher flowering rate in ‘Qinguan’ axillary bud (Fig. 10). However, some MADS-box genes (i.e., SVP, AGL15 and AGL18 as floral repressors) had an opposite result between QA and FA (Fig. 8, 10). Indeed, even for the MIKC^c-type family genes, there were subdivided into 14 monophyletic clades (i.e., AGL6, SEP, AP1, SOC1, AP3/PI, SVP, and AGL15) play positive or negative regulatory roles in floral induction [30]. The AGL15 and AGL18, along with SVP, were necessary to block initiation of floral induction [32]. AGL15 repressed flowering associated with recruitment of histone deacetylase complex components (SIN3/HDAC) [48]. SHORT VEGETATIVE PHASE (SVP) was a key flowering regulator that repressed the FT and SOC1 expression by binding to their regulatory sequences [49].

High germination rate and outgrowth was the guarantee of high flowering rate for multiple woody plants [3, 50]. Phytohormone including strigolactone (SL), auxin, cytokinin, abscisic acid (ABA), and gibberellin (GA) with the antagonistic or cooperative complex regulation network play key roles in outgrowth of axillary bud [51–55]. Our data showed that the cytokinin (CK) level and germination rate were significantly higher in axillary bud with profusely flowering ‘Qinguan’ than in no flowering ‘Nagafu No. 2’ (Fig. 11), but ABA and GA levels had an opposite result (Fig. 11), indicating that multiple hormones had negative or positive function on axillary bud germination and flowering. Studies showed that SLs repressed axillary bud outgrowth associated with SL signaling genes [53, 56], with a similar result that the significantly downregulated genes of two MAX2 and three D14 were in axillary bud with profusely flowering ‘Qinguan’ than in no flowering ‘Nagafu No. 2’, but no differences of genes belong to strigolactone biosynthesis pathway (Fig. 11). There was a contrasting behavior of auxin and cytokinin in control of axillary bud development and outgrowth in strawberry [55], and the model in Arabidopsis showed that high auxin content inhibited the activity of axillary bud, while cytokinin had an opposite effect [57]. Among them, auxin efflux PIN genes played a positively role in IAA-mediated axillary bud outgrowth by discharging transport of high auxin from the axillary bud meristem [55, 57]. Indeed, in axillary bud meristem, the eight PIN genes (i.e., two PIN1, one PIN3, and five PIN-like) were significantly higher expression in QA than in FA (Fig. 11) with relatively higher germination and flowering rates in ‘Qinguan’ axillary bud, indicating that whether IAA transported out or not from the axillary bud meristem was the key factor of axillary bud outgrowth. Researched that CsBRC1 inhibited axillary bud outgrowth by directly repressing CsPIN3 expression in cucumber [54]. Auxin transport by PIN1 gene, contributed to strigolactone-mediated shoot branching control process [58]. Additionally, BRC1 was upstream of ABA involved in regulation of bud axillary outgrowth by upregulating of ABA synthesizes key NCED3 genes expression, contributing to higher levels of ABA [51, 59], a similar result was showed in our data with upregulated of three NCED1, two NCED4, one ZEP, two ABA2 genes belong to ABA synthesizes pathway in FA with lower germination rate and no flowering than in QA, contributing to higher ABA level in FA (Fig. 11). Moreover, study showed that GA play key roles in regulation of axillary bud outgrowth involved the interaction between DELLA and BRC1 proteins [52, 60].

In a word, multi-layered regulation process, such as multiple key flowering genes with complex networks, MIKC^c-type MADS-box genes, as well as phytohormone with complex regulation network involved in axillary bud outgrowth by DNA methylation and epigenetic modification were contributed to regulate ‘Qinguan’ axillary buds flowering in apple (Fig. 12C).

During the research, a number of investigations have been conducted to confirm whether there are any epigenetic changes between axillary bud with profusely flowering ‘Qinguan’ (QA) and no flowering ‘Nagafu No. 2’ (FA) apple varieties, and whether the DNA methylation changes contribute to the expression patterns of genes involved in apple floral induction and development. We constructed whole-genome bisulfite sequencing (BS) libraries of QA and FA with diverse flowering capabilities to find whole-genome cytosine methylation profiles. Additionally, combining the RNA sequencing of QA and FA with diverse flowering capabilities, we identified the gene expression patterns and the correlation with the methylation states, which contributed to a better understanding of epigenetic regulation mechanisms of floral induction and formation in apple.

Plant materials and sample preparation

The 12-year-old ‘Qinguan’ and ‘Fuji’ apple trees were growing in the apple orchard (the experiment site of the Northwest A&F University) in Fufeng County, Baoji City, Shaanxi, China. The axillary bud with profusely flowering ‘Qinguan’ and no flowering ‘Nagafu No. 2’ were collected from 18 trees in July 2th, 2018, respectively (Fig. 1A, C), and then immediately frozen in liquid nitrogen. The axillary bud samples for both ‘Qinguan’ and ‘Fuji’ were future used for RNA sequencing and whole-genome BS libraries (to analyze DNA methylation) construction, as well as qRT-PCR to verify the genes expression with three biological replicates.

Morphological and phenotypic analysis of the ‘Qinguan’ and ‘Fuji’ axillary buds

We used a balance and a digital-display Vernier caliper to evaluate the size including length and width, and fresh weight of the axillary buds for both ‘Qinguan’ and ‘Fuji’ (Fig. 1C, D). Additionally, ~ 360 axillary buds from 12 annual apple shoots with similar length in QA and FA, respectively, were used for calculating flowering (April 10th, 2019) and germination rates (March 10th, 2019) with previously described method [3, 50], with twelve biological replicates (Fig. 1A, B). Longitudinal structural sections of ‘Qinguan’ and ‘Fuji’ axillary buds were observed with paraffin section as previously described [3].

Measurement of hormones contents

~ 0.5 g (fresh weight) from axillary buds of ‘Qinguan’ and ‘Fuji’ samples were used for hormones (i.e., IAA, CTK, ABA and GA) extracting with previously described method [61]. And then we used a high-performance liquid chromatography system (Waters 2498 UV-Visible detector, USA) for detecting and analysis the hormones levels with as previously described method [62], with three biological replicates.

Bisulfite sequencing library construction, sequencing, and data analysis

We used a modified CTAB DNA extraction method for BS libraries construction from the ‘Qinguan’ and ‘Fuji’ axillary buds [63]. Briefly, ~ 5.2 µg genomic DNA was fragmented to generate 200–300 bp fragments, and then underwent an end-repair and adenylation step for sequencing with the Illumina HiSeq 2500 platform (90 cycles) by Novogene (Beijing, China) [64, 65]. The clean reads by removing the adapters and low-quality reads from the raw data for both QA and FA were mapped to the Malus x domestica reference genome (https://www.rosaceae.org/species/malus/malus_x_domestica/genome_GDDH13_v1.1) [36].

The DNA methylation levels of the ‘Qinguan’ and ‘Fuji’ axillary buds were analyzed by a Python script (P < 0.0025) [3]. With the following formula: ML (C) = reads (mC)/reads (mC + umC). The calculated ML was then corrected as follows, according to a previous study: ML_corrected = (ML − r)/(1 − r), in which the ‘r’ is the bisulfite non-conversion rate for estimating the corrected ML [3].

Analysis of differentially methylated regions

The differentially methylated regions (DMRs) between the QA and FA were analyzed by the DSS software [66]. Briefly, the DMR distribution in different gene regions involved in the gene-body region [from the transcription start site (TSS) to the transcription termination site (TTS)] or the promoter region (2-kb upstream from the TSS) were analyzed with previously described method [67]. And we used Fisher’s exact test for calculating methylation levels between QA and FA with methylKit, with P-values. The identified DMRs based on the following four criteria: 1) at least three methylated cytosines in all the contexts; 2) more than 10× coverage; 3) bases (> 99.9th percentile coverage) were removed to minimize errors; and 4) q-values < 0.01.

RNA sequencing and identification of differentially expressed genes

Six ‘Qinguan’ and ‘Fuji’ axillary buds with three biological replicates were used for constructing RNA-seq libraries, and then sequenced with the Illumina HiSeq 2500 system by Novogene (Beijing, China) [68]. Through eliminating adapter sequences, mismatches, and low-quality reads, the clean reads were mapped to the Malus x domestica reference genome (https://www.rosaceae.org/species/malus/malus_x_domestica/genome_GDDH13_v1.1) [36] using the SOAPaligner/SOAP2 software (http://soap.genomi cs.org.cn/soapaligner.html) [69]. The gene expression levels were calculated according to the fragments per kilobase of transcript per million mapped reads (FPKM) values. The DEGs (Differentially expressed genes) with FPKM were identified according to two criteria: 1) log₂ fold-change > 0; 2) padj < 0.05 [70]. The Pearson correlation between the RNA-seq data and DNA methylation levels in each ‘Qinguan’ and ‘Fuji’ axillary buds were analyzed by previously described method [71]

GO and KEGG functional enrichment analyses of DEGs and DMRs

Significantly enriched gene ontology (GO) terms for both DEGs and DMR-associated genes between the ‘Qinguan’ and ‘Fuji’ axillary buds were identified (P < 0.05) (http://www.geneontology.org/) [72]. Additionally, the significantly enriched KEGG pathways involved in DEGs and DMRs between QA and FA were also identified (P < 0.05) (http://fttp.genome.jp/pub/kegg/pathway/) [73].

Identification of the related genes expression levels by qRT-PCR

We used a qRT-PCR assay to identify the relative gene expression levels including DEGs and DMR-associated genes between the QA and FA with three biological replicates as previously described method [2]. Additionally, the information of primers used in this study are listed in Table S6.

Statistical analysis

The DPS (v7.0) (Zhejiang University, Hangzhou, China) was used for calculating the statistical significance of the differences in the data with t-test, p < 0.05 between QA and FA. The different lowercase letters (a, b, or c) were used for indicating significant differences.

QA: Qinguan; FA: Fuji; DMRs: Different Dethylation Regions; DEGs: Different Expression genes; GA: Gibberellic Acid

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The datasets used and analysed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests.

Funding

This study was financially supported by the National Natural Science Foundation of China (31801813), the Science and Technology Innovative Engineering Project in Shaanxi province, China (2017ZDXM-NY-017), and the China Postdoctoral Science Foundation (2018M631207, 2017M623254). The funding agencies played no role in the design of the study, data collection, analysis, and interpretation or writing the manuscript.

Author contributions

The study was designed and conceived by LX. Experiments were conducted by LX, HZ, SQ, MH and NA. The data were analyzed by LX, AN, WL, WZ and HZ. The manuscript was written by LX. All the authors have read and approved the manuscript.

Acknowledgments

The authors thank the lab members for their assistance.

Xing L, Zhang D, Zhao C, et al. Shoot bending promotes flower bud formation by miRNA-mediated regulation in apple (Malus domestica Borkh.). Plant Biotechnology Journal. 14: 2016; 749–770.
Xing L, Zhang D, Qi S, et al. Transcription profiles reveal the regulatory mechanisms of spur bud changes and flower induction in response to shoot bending in apple (Malus domestica Borkh.). Plant Molecular Biology. 2019; 99: 45–66.
Xing L, Li Y, Qi S, et al. Comparative RNA-Sequencing and DNA Methylation Analyses of Apple (Malus domestica Borkh.) Buds with Diverse Flowering Capabilities Reveal Novel Insights into the Regulatory Mechanisms of Flower Bud Formation. Plant & Cell Physiology. 2019; 60: 1702–1721.
Chen X, Qi S, Zhang D, et al. Comparative RNA-sequencing-based transcriptome profiling of buds from profusely flowering “Qinguan” and weakly flowering “Nagafu no. 2” apple varieties reveals novel insights into the regulatory mechanisms underlying floral induction. BMC plant biology. 2018; 18: 370.
Haberman A, Ackerman M, Crane O, Kelner JJ, Costes E, Samach A. Different flowering response to various fruit loads in apple cultivars correlates with degree of transcript reaccumulation of a TFL1-encoding gene. Plant Journal. 2016; 87: 161–173.
Xing L, Zhang D, Song X, et al. Genome-Wide Sequence Variation Identification and Floral-Associated Trait Comparisons Based on the Re-sequencing of the ‘Nagafu No. 2’ and ‘Qinguan’ Varieties of Apple (Malus domestica Borkh.). Frontiers in Plant Science. 2016; 7.
Albani MC, Coupland G. Comparative analysis of flowering in annual and perennial plants In: Current topics in developmental biology. 2010; Elsevier, 323–348.
Xing L-B, Zhang D, Li Y-M, et al. Transcription Profiles Reveal Sugar and Hormone Signaling Pathways Mediating Flower Induction in Apple (Malus domestica Borkh.). Plant and Cell Physiology. 2015; 56: 2052–2068.
Kim J-Y, Oh JE, Noh Y-S, Noh B. Epigenetic control of juvenile-to-adult phase transition by the Arabidopsis SAGA-like complex. The Plant Journal. 2015; 83: 537–545.
Xu M, Hu T, Smith MR, Poethig RS. Epigenetic Regulation of Vegetative Phase Change in Arabidopsis. The Plant Cell. 2016; 28: 28–41.
Narsai R, Secco D, Schultz MD, Ecker JR, Lister R, Whelan J. Dynamic and rapid changes in the transcriptome and epigenome during germination and in developing rice (Oryza sativa) coleoptiles under anoxia and re-oxygenation. The Plant Journal. 2017; 89: 805–824.
Han Q, Bartels A, Cheng X, et al. Epigenetics Regulates Reproductive Development in Plants. Plants. 2019; 8: 564.
Seymour G, Poole M, Manning K, King GJ. Genetics and epigenetics of fruit development and ripening. Current Opinion in Plant Biology. 2008; 11: 58–63.
Liu R, How-Kit A, Stammitti L, et al. A DEMETER-like DNA demethylase governs tomato fruit ripening. Proceedings of the National Academy of Sciences of the United States of America. 112: 2015; 10804–10809.
Gallusci P, Hodgman C, Teyssier E, Seymour GB. DNA Methylation and Chromatin Regulation during Fleshy Fruit Development and Ripening. Frontiers in Plant Science. 2016; 7: 807.
Wang Z, Meng D, Wang A, et al. The Methylation of the PcMYB10 Promoter Is Associated with Green-Skinned Sport in Max Red Bartlett Pear1[C][W]. Plant Physiology. 2013; 162: 885–896.
Xu J, Xu H, Xu Q, Deng X. Characterization of DNA Methylation Variations During Fruit Development and Ripening of Sweet Orange. Plant Molecular Biology Reporter. 2015; 33: 1–11.
Agustí M, Mesejo C, Muñoz-Fambuena N, et al. Fruit-dependent epigenetic regulation of flowering in Citrus. New Phytologist. 2020; 225: 376–384.
Amasino RM, Michaels SD. The timing of flowering. Plant Physiology. 2010; 154: 516–520.
Andrés F, Coupland G. The genetic basis of flowering responses to seasonal cues. Nature Reviews Genetics. 2012; 13: 627.
Jung C. Flowering time regulation: Agrochemical control of flowering. Nature Plants. 2017; 3: 17045.
Dorca-Fornell C, Gregis V, Grandi V, Coupland G, Colombo L, Kater MM. The Arabidopsis SOC1-like genes AGL42, AGL71 and AGL72 promote flowering in the shoot apical and axillary meristems. The Plant Journal. 2011; 67: 1006–1017.
Johan, Zicola, Liangyu, et al. Targeted DNA methylation represses two enhancers of FLOWERING LOCUS T in Arabidopsis thaliana. Nature plants. 2019.
Qin Z, Wu J, Geng S, et al. Regulation of FT splicing by an endogenous cue in temperate grasses. Nature Communications. 2017; 8: 14320.
Jeong J-H, Song H-R, Ko J-H, et al. Repression of FLOWERING LOCUS T chromatin by functionally redundant histone H3 lysine 4 demethylases in Arabidopsis. PLoS one. 2009; 4: e8033.
He Y. Enabling photoperiodic control of flowering by timely chromatin silencing of the florigen gene. Nucleus. 2015; 6: 179–82.
Hou X, Zhou J, Liu C, Liu L, Shen L, Yu H. Nuclear factor Y-mediated H3K27me3 demethylation of the SOC1 locus orchestrates flowering responses of Arabidopsis. Nature Communications. 2014; 5: 4601.
Meijón M, Cañal MJ, Valledor L, Rodríguez R, Feito I. Epigenetic and physiological effects of gibberellin inhibitors and chemical pruners on the floral transition of azalea. Physiologia Plantarum. 2011; 141: 276–288.
Smaczniak C, Immink RGH, Angenent GC, Kaufmann K. Developmental and evolutionary diversity of plant MADS-domain factors: insights from recent studies. Development. 2012; 139: 3081–3098.
Chen F, Zhang X, Liu X, Zhang L. Evolutionary Analysis of MIKCc-Type MADS-Box Genes in Gymnosperms and Angiosperms. Frontiers in Plant Science. 2017; 8: 895.
Zhao T, Holmer R, Bruijn S de, Angenent GC, van den Burg HA, Schranz ME. Phylogenomic Synteny Network Analysis of MADS-Box Transcription Factor Genes Reveals Lineage-Specific Transpositions, Ancient Tandem Duplications, and Deep Positional Conservation. Plant Cell. 2017; tpc.00312.2017.
Fernandez DE, Wang C-T, Zheng Y, et al. The MADS-Domain Factors AGAMOUS-LIKE15 and AGAMOUS-LIKE18, along with SHORT VEGETATIVE PHASE and AGAMOUS-LIKE24, Are Necessary to Block Floral Gene Expression during the Vegetative Phase. Plant Physiology. 2014; 165: 1591–1603.
Sacharowski SP, Gratkowska DM, Sarnowska EA, et al. SWP73 Subunits of Arabidopsis SWI/SNF Chromatin Remodeling Complexes Play Distinct Roles in Leaf and Flower Development. Plant Cell. 2015; tpc.15.00233.
Pazhouhandeh M, Molinier J, Berr A, Genschik P. MSI4/FVE interacts with CUL4-DDB1 and a PRC2-like complex to control epigenetic regulation of flowering time in Arabidopsis. Proceedings of the National Academy of Sciences of the United States of America. 2011; 108: 3430–3435.
Li D, Zhang H, Mou M, et al. Arabidopsis Class II TCP Transcription Factors Integrate with the FT–FD Module to Control Flowering. Plant Physiology. 2019; 181: 97–111.
Daccord N, Celton J-M, Linsmith G, et al. High-quality de novo assembly of the apple genome and methylome dynamics of early fruit development. Nature Genetics. 2017; 49: 1099–1106.
Wang L, Kong D, Lv Q, et al. Tetrahydrofolate modulates Floral Transition through Epigenetic Silencing. Plant Physiology. 2017; 174: 1274–1284.
Huijser P, Schmid M. The control of developmental phase transitions in plants. Development. 2011; 138: 4117–4129.
Monniaux M, Mckim SM, Cartolano M, et al. Conservation vs divergence in LEAFY and APETALA1 functions between Arabidopsis thaliana and Cardamine hirsuta. New Phytologist. 2017; 216.
Xing M-Q, Zhang Y-J, Zhou S-R, et al. Global Analysis Reveals the Crucial Roles of DNA Methylation during Rice Seed Development. Plant Physiology. 2015; 168: 1417–1432.
Chahtane Hicham, Zhang Bo, Norberg Mikael. LEAFY activity is post-transcriptionally regulated by BLADE ON PETIOLE2 and CULLIN3 in Arabidopsis. New Phytologist. 2018.
Yao T, Park BS, Mao H-Z, et al. Regulation of flowering time by SPL10/MED25 module in Arabidopsis. New Phytologist. 2019; 224: 493–504.
Hyun Y, Richter R, Vincent C, Martinez-Gallegos R, Porri A, Coupland G. Multi-layered Regulation of SPL15 and Cooperation with SOC1 Integrate Endogenous Flowering Pathways at the Arabidopsis Shoot Meristem. Developmental Cell. 2016; S1534580716301964.
Fei Huang, Wenya Yuan, Shu Tian, Qijie Zheng, Yuehui He. SIN3 LIKE genes mediate longã࿽¾ay induction of flowering but inhibit the floral transition in short days through histone deacetylation in Arabidopsis. The Plant Journal. 2019.
Lopez-Gonzalez L, Mouriz A, Narro-Diego L, et al. Chromatin-Dependent Repression of the Arabidopsis Floral Integrator Genes Involves Plant Specific PHD-Containing Proteins. Plant Cell. 2014; 26: 3922–3938.
Pastore JJ, Limpuangthip A, Yamaguchi N, et al. LATE MERISTEM IDENTITY2 acts together with LEAFY to activate APETALA1. Development. 2011; 138: 3189–3198.
Rothkegel K, Sánchez E, Montes C, et al. DNA methylation and small interference RNAs participate in the regulation of MADS-box genes involved in dormancy in sweet cherry (Prunus avium L.). Tree Physiology. 2017; 37: 1739–1751.
Hill K, Huai Wang, Sharyn E. Perry. A transcriptional repression motif in the MADS factor AGL15 is involved in recruitment of histone deacetylase complex components. Plant Journal. 2008; 53: 172–185.
Shen L, Kang YGG, Liu L, Yu H. The J-Domain Protein J3 Mediates the Integration of Flowering Signals in Arabidopsis. The Plant Cell. 2011; 23: 499–514.
Yang W, Pallas B, Durand J-B, Martinez S, Han M, Costes E. The impact of long-term water stress on tree architecture and production is related to changes in transitions between vegetative and reproductive growth in the “Granny Smith” apple cultivar. Tree Physiology. 2016; 36: 1369–1381.
Gonzalez-Grandio E, Poza-Carrion C, Sorzano COS, Cubas P. BRANCHED1 Promotes Axillary Bud Dormancy in Response to Shade in Arabidopsis. Plant Cell. 2013; 25: 834–850.
Ni J, Gao C, Chen MS, Pan BZ, Xu ZF. Gibberellin Promotes Shoot Branching in the Perennial Woody Plant Jatropha curcas. Plant and Cell Physiology. 2015; 56.
Stes E, Depuydt S, De Keyser A, et al. Strigolactones as an auxiliary hormonal defence mechanism against leafy gall syndrome in Arabidopsis thaliana. Journal of Experimental Botany. 2015; 66: 5123–5134.
Shen J, Zhang Y, Ge D, et al. CsBRC1 inhibits axillary bud outgrowth by directly repressing the auxin efflux carrier CsPIN3 in cucumber. Proceedings of the National Academy of Sciences. 116: 2019; 17105–17114.
Qiu Y, Guan SC, Wen C, Li P, Gao Z, Chen X. Auxin and cytokinin coordinate the dormancy and outgrowth of axillary bud in strawberry runner. BMC Plant Biology. 2019; 19: 528.
Wang B, Smith SM, Li J. Genetic Regulation of Shoot Architecture. Annual Review of Plant Biology. 2018; 69: 437–468.
Waldie T, Leyser O. Cytokinin Targets Auxin Transport to Promote Shoot Branching. Plant Physiology. 2018; 177: 803–818.
van Rongen M, Bennett T, Ticchiarelli F, Leyser O. Connective auxin transport contributes to strigolactone-mediated shoot branching control independent of the transcription factor BRC1 (S Hake, Ed.). PLOS Genetics. 2019; 15: e1008023.
González-Grandío E, Pajoro A, Franco-Zorrilla JM, Tarancón C, Cubas P. Abscisic acid signaling is controlled by a BRANCHED1/HD-ZIP i cascade in Arabidopsis axillary buds. Proceedings of the National Academy of Sciences of the United States of America. 2016; 114: 201613199.
Davière J-M, Wild M, Regnault T, et al. Class I TCP-DELLA Interactions in Inflorescence Shoot Apex Determine Plant Height. Current Biology. 2014; 24: 1923–1928.
Dobrev PI, Vankova R. Quantification of abscisic acid, cytokinin, and auxin content in salt-stressed plant tissues In: Plant Salt Tolerance. Springer. 2012; 251–261.
Djilianov DL, Dobrev PI, Moyankova DP, et al. Dynamics of endogenous phytohormones during desiccation and recovery of the resurrection plant species Haberlea rhodopensis. Journal of Plant Growth Regulation. 2013; 32: 564–574.
Abdel-Latif A, Osman G. Comparison of three genomic DNA extraction methods to obtain high DNA quality from maize. Plant Methods. 2017; 13: 1.
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011; 27: 1571–1572.
Feng H, Conneely KN, Wu H. A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data. Nucleic Acids Research. 2014; 42: e69.
Park Y, Wu H. Differential methylation analysis for BS-seq data under general experimental design. Bioinformatics. 2016; 32: 1–10.
Wu H, Xu T, Feng H, et al. Detection of differentially methylated regions from whole-genome bisulfite sequencing data without replicates. Nucleic Acids Research. 2015; 43: e141.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-SEq. Nature Methods. 2008; 5: 621–628.
Li R, Yu C, Li Y, et al. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics (Oxford, England). 2009; 25: 1966–1967.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome biology. 11: 2010; R106.
Sánchez-Taltavull D, Ramachandran P, Lau N, Perkins TJ. Bayesian Correlation Analysis for Sequence Count Data. PLoS ONE. 2016; 11.
Du Z, Zhou X, Ling Y, Zhang Z, Su Z. agriGO: a GO analysis toolkit for the agricultural community. Nucleic acids research. 2010; 38: W64–W70.
Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Research. 2016; 44: D457–D462.

No competing interests reported.

Supplementarydata1.docx
Supplementary File 1 Table S1. Details regarding the bisulfite sequencing reads for the ‘Fuji’ and ‘Qinguan’ axillary buds. Table S2. Bisulfite sequencing data of mapped reads for the ‘Fuji’ and ‘Qinguan’ axillary buds. Table S3. The detail information of gene coverage for the ‘Fuji’ and ‘Qinguan’ axillary buds. Table S4. Details regarding the RNA-sequencing data for the ‘Fuji’ and ‘Qinguan’ axillary buds. Table S5. Details regarding the mapping of RNA-sequencing data for the ‘Fuji’ and ‘Qinguan’ axillary buds. Table S6. Primers used in this study. Figure S1. Detailed micrographs of longitudinal structural sections of ‘Fuji’ and ‘Qinguan’ axillary buds. Figure S2. Additional details regarding of epigenome identification of ‘Qinguan’ and ‘Fuji’ axillary buds. Figure S3. Bisulfite sequencing data for the mean depth and proportion of covered bases in each chromosome for the ‘Fuji’ and ‘Qinguan’ axillary buds. Figure S4. Methylation densities (A) levels (B) and in the CG, CHG, and CHH contexts in the ‘Fuji’ and ‘Qinguan’ axillary buds. Figure S5. DNA methylation patterns in various genomic regions in ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Figure S6. Additional details regarding of comparative analysis of DNA methylation levels in various genomic regions between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Circos plots comparing DNA methylation levels in all contexts between QA and FA. Figure S7. Go enrichment terms analysis of the DMR_genes in the CG, CHG, and CHH contexts between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Figure S8. Go enrichment terms analysis of the DMR_genes in the CG, CHG, and CHH contexts between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Figure S9. KEGG enrichment pathways of the DMR_genes in the CG, CHG, and CHH contexts between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Figure S10. KEGG enrichment pathways of the DMR_promoter_genes in the CG, CHG, and CHH contexts between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Figure S11. The gene expression distribution, as well as Pearson correlation among samples of ‘Fuji’ and ‘Qinguan’ axillary buds. Figure S12. Phylogenetic analysis of Apple and Arabidopsis thaliana MADS-box genes. Protein designations: Arabidopsis (At, blue circle) and Apple (Md, red square).
Supplementarydata2.xlsx
Supplementary File 2 Additional file 1. Detailed information regarding the DEGs associated with MADS-box family genes between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Additional file 2. Detailed information regarding the DEGs in multiple flowering pathways between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Additional file 3. Detailed information regarding the DEGs in autonomous, vernalization pathways, as well as different composite components between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud. Additional file 4. Detailed information regarding the DEGs in hormone pathways between ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud.
TOC.png

Download PDF

Version 1

posted

You are reading this latest preprint version

Epigenomic Regulatory Mechanism of Flowering in Apple Demonstrated in Two Varieties with Contrasted Flowering Behaviors

Status:

Version 1

Abstract

Figures

Background

Results

Phenotypes characteristics of ‘Fuji’ and ‘Qinguan’ axillary buds and their flowering ability

Genome-wide examination of cytosine methylation in ‘Qinguan’ and ‘Nagafu No. 2’ axillary bud

Distribution of DNA methylation in genic regions

Identification of differentially methylated regions (DMRs) between axillary bud with profusely flowering ‘Qinguan’ and no flowering ‘Nagafu No. 2’, and functional annotation analysis

Transcriptome profiling and functional analysis of differentially expressed genes between axillary bud with profusely flowering ‘Qinguan’ and no flowering ‘Nagafu No. 2’ apple

Association between DNA methylation and transcript abundance

Huge differences in transcription of MIKCc-Type MADS-box genes by changing DNA methylation, contributed to contrasted flowering behaviors

Differentially expression genes in multiple flowering pathways between axillary buds of ‘Qinguan’ and ‘Nagafu No. 2’ apple by changing DNA methylation, contributed to contrasted flowering behaviors

Significant differences in multiple hormone levels were due to DEGs and their DMRs in their synthesis pathways, leading to both contrasted axillary bud outgrowth and flowering behaviors.

Identification of related differentially expressed genes by a qRT-PCR assay

Discussion

Conclusion

Materials And Methods

Plant materials and sample preparation

Morphological and phenotypic analysis of the ‘Qinguan’ and ‘Fuji’ axillary buds

Measurement of hormones contents

Bisulfite sequencing library construction, sequencing, and data analysis

Analysis of differentially methylated regions

RNA sequencing and identification of differentially expressed genes

GO and KEGG functional enrichment analyses of DEGs and DMRs

Identification of the related genes expression levels by qRT-PCR

Statistical analysis

Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and materials

Competing interests

Funding

Author contributions

Acknowledgments

References

Additional Declarations

Supplementary Files

Status:

Version 1