Dynamic changes in methylome and transcriptome patterns in response to methyltransferase inhibitor 5-azacytidine treatment in citrus.

DNA methylation is known to play an important role in various developmental processes in plants. However, there is a general lack of understanding about the possible functions of DNA methylation in fruit trees. Using callus as a model, methylome, transcriptome and metabolite changes were assessed after treatment with the DNA methyltransferase inhibitor 5-azacytidine (5azaC). Genome-wide methylome analysis revealed the demethylation of a diverse of genes, including many genes encoding transcription factors (TFs), genes involved in biological processes, and the up-regulation of a wide range of transposable elements (TEs). Combined with the RNA-seq data, we observed no obvious genome-wide correlation between the changes in methylation status and expression levels. Furthermore, 5azaC treatment induced carotenoid degradation along with strong activation of carotenoid cleavage dioxygenases 1 (CpCCD1). Functional complementation analysis in bacterial system showed that CpCCD1 exhibited strong catalytic activities toward zeaxanthin, β-carotene and lycopene. In summary, 5azaC treatments induced carotenoid degradation by CpCCD1 activation and led to a genome-wide demethylation effect.

1. Introduction

In higher plants, DNA can be methylated at cytosine bases in all sequence contexts (CG, CHG and GHH, where H is A, C or T). Cytosine methylation is regulated by a diversity of known factors, including DNA methyltransferases (MET1, CMTs and DRMs), demethylases (DME, ROS1, DML2 and DML3) and the chromatin-remodelling factor DDM1. Proper regulation of cytosine methylation is essential for silencing transposable elements (TEs), moderating appropriate gene expression and maintaining proper chromatin structure. Much of our understanding of the role of DNA methylation in plant processes derives from experiments utilizing 5-azacytidine (5azaC), a chemical analogue of cytidine that inhibits DNA methyltransferases and reduces genomic methylation levels by direct incorporation into DNA and substitutions for cytosine. For example, 5azaC treatment accelerated the flowering process in the long-day plant Silene armeria and the short-day plant Perilla frutescens under non-inductive photoperiodic conditions. Treatment of tomato fruit could accelerate the ripening process by enhancing the transcriptional abundances of the hallmark ripening genes pg2a (encoding the pectinase polygalacturonase) and psy1 (phytoene synthase). Further methylome analysis of the tomato fruit ripening process revealed that DNA methylation triggered the fruit ripening process and regulated the quality of fruit formation, such as the cell wall, colour and aroma. Moreover, the involvements of DNA methylation regulating secondary metabolism has also been described in the plant kingdom. For example, in apple and pear, DNA methylation acts as a regulator at the MYB promoter region to determinate peel colour via different levels of anthocyanin accumulation., Similarly, different DNA methylation status effects on the chalcone synthase gene determine the anthocyanin pigmentation in the floral tissues of two Oncidium orchid cultivars. In Vitis amurensis, 5azaC treatment of the rolB-transgenic callus, which accumulates higher resveratrol contents, further increased the resveratrol production, indicating that DNA methylation could also influence the secondary metabolism.

Carotenoids are a diverse group of natural products that are widely distributed throughout plants, algae, fungi and bacteria. Carotenoids not only comprise colourful pigments in flowers and fruits but they also play important roles in photosynthesis, antenna assembly and photoprotection. The precursor of carotenoid biosynthesis is geranylgeranyl diphosphate (GGPP), which is synthesized via the 2-C-methyl-d-erythritol (MEP) pathway., Two molecules of GGPPs are synthesized by phytoene synthase (PSY) to generate phytoene. This procedure is widely accepted as a rate-limiting step in the carotenoid pathway. Four additional enzymes act to convert phytoene to all trans-lycopene: phytoene desaturase (PDS), ζ-carotene isomerase (Z-ISO), ζ-carotene desaturase (ZDS), and carotenoid isomerase (CRTISO)., The synergistic actions of lycopene β-cyclase (LCYB) and lycopene ε-cyclase (LCYE) then determine the lycopene flux that produces lutein, β-carotene, and the xanthophyll cycle carotenoids. Carotenoids can be enzymatically cleaved to produce apocarotenoids, which are commonly known as abscisic acid (ABA), strigolactones and volatile compounds . Two types of carotenoid cleavage enzymes have been identified in plants: 9-cis carotenoid cleavage dioxygenases (NCEDs), which catalyze the first step in ABA biosynthesis, and carotenoid cleavage dioxygenases (CCDs), which cleave double bonds at different positions in various carotenoid substrates, giving rise to volatile and non-volatile apocarotenoids. CCD1 cloned from maize, rice, Rosa damascene and other plants can cleave a variety of carotenoid substrates in vitro, indicating that the cleavage positions of CCD1 are not regioselective. In citrus, CitCCD1 has been reported to cleave at the 9–10 and 9′–10′ positions of β-cryptoxanthin, zeaxanthin and all-trans-violaxanthin and at the 9′–10′ position of 9-cis-violaxanthin., However, another member of the CCD family in citrus, CitCCD4, has shown cleaving specificity for the 7–8 or 7′–8′ positions of β-cryptoxanthin and zeaxanthin, leading to the formation of β-citraurin.

Citrus fruits accumulate an abundance of carotenoids and present various colours with different carotenoid accumulation patterns. To understand the basis of carotenoid biosynthesis, we expressed the phytoene synthase gene CrtB from Erwinia herbicola in citrus callus in vitro, resulting in an accumulation of high levels of carotenoids., To explore the possible involvement of DNA methylation in carotenoid metabolism, we treated the CrtB-transgenic callus (named M33) with incremental concentrations of the DNA methyltransferase inhibitor 5azaC. The 5azaC treatments triggered carotenoid degradation and thus decreased carotenoid contents. Furthermore, analysis of the methylome and transcriptome patterns was performed to identify the global changes in DNA methylation associated with 5azaC treatment. The results revealed that a diverse set of genes was demethylated, including genes encoding transcription factors (TFs), and 5azaC treatment led to the up-regulation of a wide range of TEs. These data provide insight into whole-genome changes in biological processes in response to 5azaC (demethylation) treatment and improve our understanding of the regulatory mechanism underlying DNA methylation.

2. Materials and methods

2.1. Plant materials

The wild-type callus was propagated by abortive ovule embryogenic calli derived from Pink Marsh grapefruit (Citrus paradisi Macf.) and, designated as WT. The CrtB (CrtB, a bacterial phytoene synthase gene, GenBank No. M90698) transgenic callus line M33 used in this study was established by transformation of the suspended grapefruit callus (WT) with Agrobacterium tumefaciens strain EHA105 containing the PMV vector with the chimeric gene (tp-rbcS-CrtB) as described previously. The previous study showed that the total carotenoid content was markedly increased in CrtB-transgenic callus M33 compared with the non-coloured WT callus. The specimens of WT and M33 were preserved and supplied by the Key Laboratory of Horticultural Plant Biology of Ministry of Education at Huazhong Agricultural University (HZAU), Wuhan, China.

2.2. 5azaC treatments

The 5azaC reagent (Sigma-Aldrich) was dissolved in sterile double-distilled water at 40 mM, and the solution was filtered with a sterile 0.22-μm filters. After filtration, the sterile 5azaC solution was added to liquid MT to the designated concentrations. The 3 g of input grapefruit callus was added to liquid MT medium with the designated 5azaC concentrations of 0 μM (control), 10, 50, 100, 150, 200 and 400 μM with shaking (130 rpm) for 15 days. Then, the 400 μM-treated callus was cultured on normal solid MT medium for one month and transferred to liquid MT medium with shaking (130 rpm) for 15 days to maintain similar culture conditions. After treatment, the callus was harvested as the TR line.

2.3. Methylated DNA immunoprecipitation sequencing (MeDIP-seq) analysis

The genomic DNAs of CK, 400 and TR callus were used for global methylation analysis. Each sample consisted of two biological replicates for independent library construction and sequencing. The genomic DNA was sheared with Bioruptor UCD 200 (Diagenode) to 100–500 bp. The library constructions of shared DNA were performed with the DNA-seq Library Preparation Kit-Illumina Compatible (Gnomegen K02422). The Methylated-DNA IP Kit (ZYMO D5101) was applied to the methylated DNA immunoprecipitation. The library products were sequencing using an Illumina HiSeqTM 2500. After filtration, clean reads were mapped to the reference genome using BWA (Burrows-Wheeler Aligner) software. Peak-calling was performed with MACS 1.4.0 software. The differentially methylated peaks were analyzed with the edgeR package. The significant differently methylated peaks were located near different elements according to the reference genome annotation.

2.4. RNA-seq and global detection of differentially expressed genes (DEGs)

The RNAs extracted from M33 control callus (CK), 400 μM 5azaC treated callus (400) and treatment recovery callus (TR) were applied to DEG analysis. Each sample had three biological replicates for independent library construction and sequencing. The library products were sequencing by Illumina HiSeqTM 2000. After the filtration, clean reads were mapped to reference genome using SOAPaligner/SOAP2. The gene expression level was calculated by using RPKM method (Reads Per kb per Million reads). We used ‘FDR (False Discovery Rate) ≤ 0.001 and the absolute value of log2Ratio ≥ 1’ as the threshold to judge the significance of gene expression difference. After getting Gene Ontology (GO) annotation for DEGs, WEGO software was applied to analyze the GO functional classification for DEGs. Furthermore, pathway enrichment analysis (KEGG) was applied to identify significantly enriched metabolic pathways or signal transduction pathways in DEGs comparing with the whole genome background.

2.5. Carotenoid profile analysis

About 0.3 g of lyophilized callus sample was applied to extract the total carotenoids and each sample extracted three independent replications. The absorbance of total carotenoids content was analyzed by a Microplate spectrophotometer (TECAN, infinite M200). The carotenoid extraction and analysis using reversed-phase high-performance liquid chromatography (RP-HPLC) was applied as previously described., The analysis of carotenoid profile was performed on a Waters liquid chromatography system (600E solvent pump system, 2996 photodiode array detection, 717 plus autosampler). A C30 carotenoid column (150 × 4.6 mm; YMC, Japan) was used to elute the carotenoids. Each sample was performed with three independent replications and SPSS software was applied to statistical analysis (P <0.05 (*), P <0.01 (**)).

2.6. Quantitative analysis of gene expression

Total RNA was extracted from the callus according to the method described previously. First strand cDNA was synthesized from 2.0 μg of total RNA using the ReverAid first strand cDNA synthesis KIT (Fermentas). The real-time PCR primers used in this study were as listed in the previous study or designed by Primer Express software (Applied Biosystems, Foster City, CA, USA) based on the sequenced genome data. The primer sequences were listed in Supplementary Table S5. The primers were diluted in Power SYBR® Green PCR Master Mix (Applied Biosystems) and the amplification mixture volume was 10 μL per reaction. Reaction conditions were an initial incubation for 2 min at 50 °C and 95 °C for 1min, and then followed by 40 cycles of 95 °C/15 s and 60 °C/1 min. Reactions were run on a 7900HT Fast Real-Time PCR System with 384-Well Block Module (Applied Biosystems). The β-actin gene was used as an endogenous control and comparative Ct method (2−ΔΔCt) was adopted to calculate the expression data according to the previous study. SPSS software was applied to the statistical analysis (P <0.01 (**)).

2.7. Cloning of CpCCD1 and expression of CpCCD1 in carotenoid accumulating Escherichiacoli

The Gateway Cloning system was applied to the vector construction for the CpCCD1 bacterial expression analysis. The whole CDS sequence of citrus CCD1 (genome ID: Cs7g01710) was amplified by the primer pairs with the attB adapter (Forward: 5′-AAAAAGCAGGCT TC ATGGTGGAGAAAGAAAAGC-3′; Reverse: AGAAAGCTGGGT G TCACAATTTTGCTTGCTCT). The amplified fragment was used as the template to perform the second round PCR reaction with the attB-adapter primers (attB1-adapter: GGGGACA AGT TTGTACAAAAAAGCAGGCT; attB2-adapter: GGGGAC CACTTTGTACAAGAA AGCTGGGT). Then, the PCR product was cloned to the entry vector pDONR221 (Kanamycin) by performing the BP reaction and sequenced. After confirming the CpCCD1 sequence, the CpCCD1 was cloned to the final expression vector pDEST-17 (Ampicillin) by LR reaction. E. coli transformants harbouring plasmids pACCRT-EIB, pACCAR16△crtX and pACCAR25△crtX which accumulated lycopene, β-carotene and zeaxanthin, respectively, with CpCCD1-pDEST-17 were transformed to E. coli (BL21). The culture conditions and carotenoid analysis were according to the previous study.

2.8. Bisulfite sequencing analysis of CpCCD1 gene

The genomic DNA was extracted as Cheng et al. and treated with the ZYMO EZ DNA Methylation-Gold Kit (D5005). The converted gDNA was amplified with specific primers of three selected regions (M1, M2 and M3) which located in promoter and first exon of the CpCCD1. The locations of these three primer pairs were as follows: M1 (from −1689 to −1286) Forward: TTGAAGTTGGGGGAAATGTTAAAGG, Reverse: CAAAATTTARTTACACTAACCTCT; M2 (from −636 to −140) Forward: TAGATTTTGTGAGAATTATTTGAGG, Reverse: TCCTTTACRACCCACARCAAAA; M3: (+105 to +589) Forward: TGTGGAGAAAYTGATAGTGAAGTGG, Reverse: CCACAATCTCTRTACATAACTCTCT. The positive clone sequences combined with the untreated sequences were analyzed by CyMATE tool (http://www.cymate.org/).

2.9. Measurement of ABA content

About 0.2 mg callus of each sample was applied to hormone extractions. Each sample was preformed with four replications with independent extraction. 10 ng of [2H6]ABA (2-cis, 4-trans-abscisic acid-[2H6], olchemim, Czech Republic) was used as internal standards for ABA detection. The hormone extractions were separated by HPLC (high-performance liquid chromatography; Agilent 1100, Agilent Technologies, Palo Alto, CA, USA) and measured with HPLC-electrospray ionization-mass spectrometry (ESI-MS) (API3000 mass spectrometer, Applied Biosystems, Foster City, CA, USA). The MS/MS conditions were set according to Pan et al.

3. Results

3.1. Citrus methylome patterns and their variations in response to 5azaC treatment

The responses of wild-type (WT) and CrtB-transgenic (M33 line) calli to various concentrations of 5azaC (control-CK, 0 μM; 10, 50, 100, 150, 200 and 400 μM) and the treatment recovery (TR) line were investigated (Fig. 1). To evaluate global DNA methylation distributions and variations in response to 5azaC treatment, CK, 400 and TR callus lines were analyzed using methylated DNA immunoprecipitation sequencing (MeDIP-seq) (red box in Fig. 1). For each sample, two biological replicates were applied for library construction and sequencing. Methylation status of samples including CK, 400 and TR, gene and TE density were displayed in Fig. 2A. The results showed that TEs presented a hypermethylation status and were mainly concentrated in pericentromeric heterochromatin regions. In contrast, methylation in gene-rich regions were characterized by hypomethylation and mainly concentrated in euchromatic regions. Moreover, the distributions of DNA methylation level in genes and TEs were shown in Fig. 2B. We found that gene body (middle region) presented a low methylation level and the up-stream and down-stream regions of gene body showed a relatively higher methylation level. In contrast, TEs showed a higher methylation level in TE body and lower methylation levels in up-stream and down-stream regions (Fig. 2B). Concerning the methylation changes in response to 5azaC, the methylation level of 400 callus was higher than it in CK and TR lines in gene body region, but lower in up-stream and down-stream regions, corresponding to the results that the demethylation effects caused by 5azaC were mainly occurred in intergenic areas (Fig. 3A). In TE and its flanking regions, the methylation level of 400 callus was consistently lower than CK and TR lines (Fig. 2B).

Figure 1

Seven 5-azactytine (5azaC) concentration gradients were assessed: control (CK, 0 μM), 10, 50, 100, 150, 200 and 400 μM on citrus callus. The 400 μM 5azaC-treated callus (400) was cultured on normal MT medium for one month as the TR line.

Figure 2

The citrus methylome. (A) Circos plots of nine citrus chromosomes in CK, 400 and TR samples. Track order (from inner to outer rim): gene density (a); TE density (b); methylation fold enrichment of CK (c), 400 (d) and TR (e); chromosome ideograms of citrus genome (f). The chromosome name and scale are indicated on the outer rim. (B) Distributions of DNA methylation status in genes (left) and TEs (right). Flanking regions are the same length as the gene or TE body (middle region). TSS, transcription start site; TTS, transcription termination site.

Figure 3

Global methylome analysis by MeDIP-seq. (A) Numbers of differentially methylated elements (intergenic, upstream2k, CDS, intron, 3′UTR, 5′UTR, downstream2k). Differentially methylated regions were first divided into intergenic and genic region. The latter was further divided into different gene elements, including upstream2k, CDS, intron, 3′UTR, 5′UTR and downstream2k. (B) Numbers of differentially methylated genes with promoter and/or gene body. (C) Venn diagrams of differentially methylated genes among CK vs. 400, 400 vs. TR and CK vs. TR. (D) Distributions and functional categories of DMGs analyzed by MapMan.

The distribution of regions showing changes in methylation in response to 5azaC is shown in Fig. 3A. The largest number of DNA methylation variations altered by 5azaC treatment occurred in intergenic areas. To investigate whether TE genes were activated by 5azaC, 3 TE genes with high methylation enrichment and 12 additional TE genes that were randomly selected from the genome were analyzed for changes in expression. The results showed that most TE genes were up-regulated following 5azaC treatment and re-silenced in the TR lines (Supplementary Fig. S1). As expected, most of the gene elements were hypomethylated (down) in response to 5azaC treatment (CK vs. 400) and hypermethylated (up) in the TR callus (400 vs. TR). Moreover, Fig. 3B shows that variations in the methylation of genes occurred most frequently in the promoter region. In total, 3082 (682 up; 2400 down), 3516 (2752 up; 764 down) and 2086 (1171 up; 915 down) differentially methylated genes (DMGs) (including the gene body and promoter regions) were identified among CK vs. 400, 400 vs. TR and CK vs. TR, respectively (Fig. 3C). Among CK vs. 400, 1821 (76%) of the down-regulated DMGs affected by 5azaC treatment were up-regulated in the TR callus, indicating that these genes were potentially demethylated by 5azaC treatment and recovered in the TR callus.

The distributions and functional categories of all DMGs were visualized using MapMan software (Fig. 3D). DMGs affected by 5azaC treatment were mainly categorized as proteins (Proteolysis), RNA (especially TFs), signalling, stress, secondary metabolism and hormone metabolism. Most of the DMGs involved in these pathways were demethylated by 5azaC treatment, while the majority of DMGs among the 400 vs. TR exhibited hypermethylation in these pathways (Supplementary Figure S2A). In the secondary metabolism category in particular, the DMGs were mainly involved in the pathways of ‘Phenlypropanoids’, ‘Lignin and lignans’ and ‘Flavonoids’ and included several carotenoid genes such as PSY (phytoene synthase) (Supplementary Fig. S2B). Since the methylome data showed that the methylation status of TFs was disrupted by 5azaC treatment, DMGs annotated as TFs (e-value <1 e−10, identity >70%) were further identified and are detailed in Supplementary Table S1. One hundred and thirty-five TFs (18 up, 117 down; belonging to 39 families), 167 TFs (146 up, 21 down; 43 families) and 95 TFs (60 up, 35 down; 29 families) were identified among the DMGs from CK vs. 400, 400 vs. TR and CK vs. TR, respectively. Mainly, TF families that were disturbed by 5azaC treatment are shown in Supplementary Fig. S2C. Particularly, 11 bHLHs (14 bHLHs in total), 8 NACs (9), 7 ERFs (9), 8 WRKYs (8) and 8 MYBs (8) were hypomethylated in response to 5azaC treatment, while 18 bHLHs (21), 12 MYBs (13), 11 NACs (11), 7 C3Hs (9), 6 C2H2s (9), 5 ERFs (7) and 6 WRKYs (6) were hypermethylated in the TR callus.

DNA methylation is involved in various biological processes by regulating gene expressions. To explore the relationship between DNA methylation and gene expression on a genome-wide scale, genes were first classified into non-expressed (FPKM value < 0.1) and expressed genes (FPKM value > 0.1). Expressed genes were further rank-ordered based on degree of gene expression and divided into quintiles. The first quintile is the 20% of genes with the lowest expression level and the fifth group is the highest. Then, the methylation levels were calculated for each quintiles. As shown in Supplementary Fig. S3, the first quintile shows the highest methylation level in control callus (CK), which means that genes with low expression levels present hyper-methylation status in citrus genome. In 400 μM 5azaC-treated callus, methylation pattern shows the opposite way compared with CK, which the first quintile presents the lowest methylation level and the second quintile is the second lowest. These findings also suggested that the 5azaC treatment presented the biggest demethylation effect on the genes with low expression level.

3.2. Transcriptome changes in response to 5azaC treatments

The gene expression of the citrus DNA methyltransferase genes CsMET1, CsCMT1 and CsDRM1; the DNA demethylase gene CsDME1; and the chromatin-remodeling gene CsDDM1, was investigated in response to 5azaC treatment. As shown in Supplementary Fig. S4, CsMET1, CsCMT1 and CsDDM1 were significantly up-regulated following 5azaC treatment independently of the WT and CrtB-transgenic calli. Additionally, the demethylase gene CsDME1 was slightly induced in both WT and M33 calli. These findings suggested that the methylome pattern in the callus was dynamically altered by 5azaC treatment.

To investigate global gene expression in response to 5azaC treatment, we utilized RNA-seq technology to explore the expression patterns of CK, 400 mM treatment and TR. The results showed that the global gene expression patterns were altered by 5azaC treatment and recovered in the TR callus (Fig. 4A). A summary of sequencing samples and results is shown in Table 1. A total of 3407 differentially expressed genes (DEGs) were identified between CK and 400, with 2074 up-regulated and 1333 down-regulated genes. A total of 792 genes were significantly up-regulated and 1592 genes were down-regulated in 400 vs. TR (Fig. 4A). Indeed, a large numbers of activated genes (1397) in response to 5azaC treatment were also identified in the TR down-regulated gene group (Fig. 4B and Supplementary Table S2). To validate the RNA-seq data, 15 randomly selected DEGs, including 8 up-regulated and 7 down-regulated genes, were analyzed by qRT-PCR. As shown in Supplementary Fig. S5, the fold change values obtained by qRT-PCR were strongly correlated with the RNA-seq data (R2=0.79), indicating that the RNA-seq data were reliable to reflect the expression patterns. These data suggest that global expression patterns are strongly altered by 5azaC treatment.

Figure 4

Global DEG analysis by RNA-seq. (A) Numbers of DEGs among CK vs. 400, 400 vs. TR and CK vs. TR. (B) Venn diagrams of DEGs among CK vs. 400, 400 vs. TR and CK vs. TR; (C) heat map of DEGs among CK, 400 and TR samples.

Table 1

Summary of RNA-seq reads mapped to the reference genome

Sample	Raw reads	Clean reads	Total mapped (%)	Unique mapped (%)	Multiple mapped (%)	Perfect match (%)	< =2 bp Mismatch (%)
CK1	12497778	12380722	9,935,538 (80.25)	9,191,654 (74.24)	743,884 (6.01)	6,918,384 (55.88)	3,017,154 (24.37)
CK2	11880144	11772060	9,416,689 (79.99)	8,710,220 (73.99)	706,469 (6.00)	6,392,327 (54.30)	3,024,362 (25.69)
CK3	12498558	12385033	9,894,235 (79.89)	9,148,753 (73.87)	745,482 (6.02)	6,662,223 (53.79)	3,232,012 (26.10)
400-1	11900480	11800960	9,538,966 (80.83)	8,803,426 (74.60)	735,540 (6.23)	6,639,095 (56.26)	2,899,871 (24.57)
400-2	11940577	11841416	9,558,922 (80.72)	8,818,951 (74.48)	739,971 (6.25)	6,577,332 (55.55)	2,981,590 (25.18)
400-3	11952924	11842571	9,543,781 (80.59)	8,809,653 (74.39)	734,128 (6.20)	6,556,128 (55.36)	2,987,653 (25.23)
TR1	12069960	11944069	9,583,414 (80.24)	8,897,339 (74.49)	686,075 (5.74)	6,939,245 (58.10)	2,644,169 (22.14)
TR2	12249152	12161431	9,791,014 (80.51)	9,085,493 (74.71)	705,521 (5.80)	6,733,171 (55.36)	3,057,843 (25.14)
TR3	12703092	12619477	10,069,849 (79.80)	9,279,193 (73.53)	790,656 (6.27)	6,773,587 (53.68)	3,296,262 (26.12)

To understand the functions of DEGs, they were categorized into gene ontology (GO) categories. Those GOs showing a significant difference from the expected distribution were cellular process, metabolic process, response to stimulus of biological process, cell, cell part and organelle of cellular component and binding, and binding and catalytic activity of molecular function (Fig. 5A). Furthermore, assignments of DEGs in pathways were visualized using MapMan tools. Fig. 5B shows that many DEGs were involved in ‘Proteolysis’, ‘Signalling’, ‘Cell wall’ and ‘Secondary metabolites’. In particular, an abundance of genes involved in secondary metabolism were enriched in the ‘Phenlypropanoids’ and ‘Lignin and lignans’ categories. In the ‘Proteolysis’ category, the F-box and RING genes involved in proteasome displayed altered expression levels in response to 5azaC treatment, implying that 5azaC treatment modified the Proteolysis process and induced the ubiquitin-dependent protein degradation (Supplementary Fig. S6). The functions of genes involved in hormone signalling were also disrupted by 5azaC treatment (Fig. 5). In detail, most genes involved in auxin, brassinosteroid, ethylene and Jasmonic acid (JA) signalling were up-regulated, while ABA signalling-related genes were down-regulated, consistent with the observed decrease in ABA content (Fig. 8).

Figure 5

(A) GO analysis of DEGs. (B) The MapMan bin of DEGs involved in the stress response.

Figure 8

Transcriptional analysis of genes involved in isoprenoid and carotenoid metabolism and the changes in ABA content in response to 5azaC. The transcript levels of the control (CK) were used as the calibrator for expression analysis. ** indicates a significant differences at P <0.01.

Interestingly, many genes assigned to the ‘Peroxidases’ and ‘Redox state’ categories were up-regulated in response to 5azaC treatment, suggesting that this treatment induced the callus peroxidation process. We further verified the changes in expression of all differentially expressed peroxidase genes and peroxidase (POD) activity in response to 5azaC. As shown in Supplementary Fig. S7, the expression levels of 17 peroxidase genes were strongly induced in 5azaC-treated calli and restored in the TR lines. Accordingly, POD activity was enhanced by 5azaC and recovered in the TR callus. These results suggested that 5azaC induced callus stress reactions, which involved strong activation of POD genes and enzyme activities to reduce peroxide levels. Furthermore, DEGs were subjected to KEGG pathway analysis (http://www.genome.jp/kegg/) to understand the biological functions of DEGs. Many DEGs were enriched in ‘metabolites pathways’, particularly in phenylalanine and phenylpropanoid metabolism, biosynthesis of secondary metabolites, amino acid metabolism, plant hormone signal transduction, sugar metabolism and others (Supplementary Table S3). Since genes involved in sugar metabolism were disrupted by 5azaC, GC-MS analysis of callus extraction was further performed to analyze the basic primary metabolites in callus (Supplementary Fig. S8). The results showed that the contents of most types of sugars were increased, while various types of acids were decreased in response to 5azaC treatment.

As shown in Fig. 5, TFs were also activated by 5azaC (demethylation) treatment. To analyze the changes in expression of TFs under 5azaC treatment, DEGs were identified as TFs using BLASTP software (e-value <1 e−10, identity >70%) based on planttfdb database (http://planttfdb.cbi.pku.edu.cn/index.php). As a result, 231 TFs (belonging to 43 families), 153 TFs (37 families) and 57 TFs (23 families) were identified among the 3407, 2384 and 623 DEGs between CK and 400, 400 and TR and CK and TR, respectively (Supplementary Table S4). The transcriptional changes in the top 16 TF families are shown in Fig. 6A. In detail, 21 (81% of total 26 MYBs) MYBs, 24 WRKYs (100%; 24), 9 NACs (43%; 21), 11 ERFs (55%; 20) and 9 bHLHs (82%; 11) were activated by 5azaC treatment (CK vs. 400), while 17 MYBs (89%; 19), 18 WRKYs (100%; 18), 8 ERFs (67%; 12), 5 NACs (45%; 11) and 6 bHLHs (75%; 8) were down-regulated in 400 vs. TR (Supplementary Table S4). The heat map also shows that many TFs were activated in response to 5azaC treatment and restored in the TR callus (Fig. 6B). Furthermore, qRT-PCR was performed to verify the expression levels of selected TFs belonging to the MYB, WRKY, NAC and ERF families (Fig. 6C). Interestingly, all the WRKY members were activated by 5azaC treatment and re-silenced in the TR lines, indicating that the WRKY family played a role in the response to 5azaC treatment.

Figure 6

(A) Numbers of DEGs annotated to different TF families. (B) The heatmap of differentially expressed TFs in CK, 400 and TR. (C) qRT-PCR analysis of selected MYB, WRKY, NAC and ERF genes. ** indicates a significant difference at P <0.01.

To investigate the possible association of methylation changes and gene expression alterations, we analyzed the intersection of DMGs and DEGs. As shown in Supplementary Fig. S9, for CK vs. 400, the DMGs and DEGs shared 343 common genes. Among these genes, 187 of the hypo-methylated DMGs were up-regulated and 17 of hyper-methylated DMGs were down-regulated, counting that 59.5% (204 genes) of the total common genes methylation changes showed a negative correlation to their expression variations (Supplementary Fig. S9). Similarly, 286 genes were both identified in DMGs and DEGs among 400 vs. TR. Among these genes, 189 genes (66.1%), including 12 hypo-methylated/up-regulated genes and 177 hyper-methylated/down-regulated genes, were presented a negative correlation between methylation changes and expression variations (Supplementary Fig. S9). These findings revealed that gene expression alterations were not always correlated to the corresponding DNA methylation changes.

3.3. Enhanced CpCCD1 expression contributes to carotenoid degradation in response to 5azaC

The carotenoids were degraded upon 5azaC-treatment in the CrtB-transgenic callus line. The 400 μM 5azaC-treated callus (400) was cultured on normal MT medium for one month as a treatment recovery line (TR) (Fig. 1). As shown in Fig. 1, the carotenoid pigmentation of the M33 callus gradually decreased with increasing 5azaC concentrations, and was re-established after 1 month of growth on normal MT medium in the TR lines. Indeed, the total carotenoid contents decreased with increasing 5azaC concentrations and were recovered in the TR lines (Fig. 7A). Furthermore, HPLC analysis demonstrated that the M33 callus mainly accumulated β-carotene, lutein, and α-carotene, with these three carotenoids showing a 64%, 88% and 79% reduction in 400 μM 5azaC-treated callus, respectively (Fig. 7B). Accordingly, the ABA content also decreased in the 5azaC-treated callus (Fig. 8), indicating that carotenoid degradation resulted in lack of ABA biosynthesis. In WT callus, the total carotenoid contents were significantly reduced in WT callus under the 100, 200 and 400 μM 5azaC-treatments (Fig. 7A). However, the ABA contents increased with an increase in 5azaC concentrations (Supplementary Fig. S10), suggesting that the observed increase in ABA content was not caused by carotenoid metabolism.

Figure 7

The carotenoid profiles of citrus callus in response to 5azaC treatments. (A) The total carotenoid contents of CtrB-transgenic callus (M33) and wild-type callus (WT) treated with 5azaC. (B) HPLC analysis of carotenoid profiles in M33 callus affected by 5azaC. * indicates a significant difference at P <0.05, ** indicates a significant difference at P <0.01.

To assess the possibility that the reduced carotenoid concentration observed in the 5azaC-treated CrtB-transgenic callus resulted from enhanced carotenoid metabolism, we analyzed the expression of genes involved in isoprenoid and carotenoid metabolism (Fig. 8). First, we tested the expression level of exogenous gene CrtB, and the results showed that CrtB was up-regulated in response to 5azaC treatment. Additionally, 5azaC treatment resulted in a significant increase in expression of the endogenous gene DXR (1-deoxy--xylulose 5-phosphate reductoisomerase) in the isoprenoid pathway. For carotenoid pathway genes, the crucial rate-controlling gene, PSY, presented an increasing expression level with the increase in 5azaC concentrations (Fig. 8), which correlated with the demethylation of PSY in response to 5azaC (Supplementary Fig. S2B). Additionally, PDS and LCYE (lycopene ε-cyclase) were induced by 5azaC treatments. The most prominent change in expression related to 5azaC was the ∼8–12-fold increase in the expression of CpCCD1 (Fig. 8). Our phylogenetic analyses indicated that the CCD family in citrus comprised eight members. However, the expression levels of these genes were not induced by 5azaC treatment (Supplementary Fig. S11). Given the important role of CCD1 in carotenoid metabolism and the prominent effect of 5azaC on CpCCD1 expression, we focused on this gene.

We isolated a 1.6-Kb CDS sequence of CpCCD1 (GenBank KT183033) from citrus callus. This sequence would encode a 547-amino-acid protein with an estimated molecular mass of 61.8 kDa. Phylogenetic analysis of CCD proteins showed that CpCCD1 clustered with other typical CCD1 proteins and with the highest similarities (87%, 77% and 80% identity) to grape (VvCCD1), maize (ZmCCD1) and Arabidopsis (AtCCD1) CCD1 proteins, respectively (Supplementary Fig. S12A and B). To confirm the carotenoid degradation ability of CpCCD1, we expressed the full-length CDS of CpCCD1 in E. coli strains harbouring plasmids pACCAR25△crtX, pACCAR16△crtX and pACCRT-EIB, which were engineered to accumulate zeaxanthin, β-carotene and lycopene. All three carotenoid-accumulating E. coli strains showed a sharp reduction in carotenoid pigmentation when CpCCD1 was expressed (Fig. 9). HPLC analysis also indicated that CpCCD1 had a strong ability to degrade zeaxanthin, β-carotene and lycopene (Fig. 9). Taken together, these results suggest that enhanced CpCCD1 expression is at least a strong contributor to the effect of 5azaC on CrtB-transgenic callus.

Figure 9

Functional analysis of CpCCD1 protein. Escherichia coli strains engineered to accumulate zeaxanthin, β-carotene and lycopene were cotransformed with the CpCCD1 expression vector. Further carotenoid analysis of Escherichia coli strains was performed using HPLC.

To investigate whether CpCCD1 induction was due to the 5azaC-dependent change in methylation of the CpCCD1 gene, we performed bisulfite-sequencing analysis targeting three representative regions (M1, M2 and M3) located in CpCCD1. M1 and M2 were located in promoter regions and M3 was in the first exon. As shown in Fig. 10, the M1 region was hypomethylated (∼10%) in CrtB-transgenic callus. Oddly, the percentage of methylation level in 5azaC-treated callus increased to ∼23–45%, indicating that M1 region was hypermethylated in response to 5azaC (Fig. 10). This finding suggested that although the 5azaC treatment induced extensive demethylation processes in treated callus, methylation statuses were also increased at some loci. Further methylome analysis confirmed this conclusion (Fig. 3). However, the M2 and M3 regions located near the start site were barely methylated (Supplementary Fig. S13), possibly to assure the initiation of transcription. Clearly, the M1, M2 and M3 regions of CpCCD1 were not demethylated by the 5azaC effect. This finding suggests the occurrence of either mechanisms other than DNA methylation of CpCCD1 induce CpCCD1 expression, or that regulatory changes in methylation occurred outside of the tested region.

Figure 10

Methylation analysis of CpCCD1 using bisulfite sequencing. M1, M2 and M3 located in the promoter and first exon of CpCCD1 were selected.

4. Discussion

The ongoing development of high-throughput sequencing technologies provides the possibility of exploring the DNA methylation patterns on a genome-wide scale. MeDIP-seq is a powerful technology and is widely applied in plant DNA methylome analyses. To our knowledge, this study presents the first global DNA methylation analysis in citrus species based on high-throughput sequencing. The 5azaC, which is known to block DNA methylation, has been widely used in DNA methylation studies. For example, 5azaC-treated Arabidopsis seedlings display an altered gene expression profile and up-regulated expression levels of many TEs. Moreover, in tomato, 5azaC treatment on the tomato fruit activated PSY expression and triggered the tomato fruit ripening process.Vitis amurensis cultured cells treated with 5azaC showed activation of the resveratrol biosynthesis gene STS (stilbene synthase) and increased resveratrol biosynthesis. Here, we used MeDIP-seq to characterize genome-wide DNA methylation patterns in citrus and combined the methyltransferase inhibitor 5azaC to explore the regulatory roles of DNA methylation in the biological processes. The results revealed an overall demethylation effect across the whole genome, which was accompanied by TE activation (Supplementary Fig. S2). These data will provide a basis for understanding the regulatory roles of DNA methylation in various biological processes in citrus species.

DNA methylation has generally been considered to be a ‘silencing’ epigenetic mark for gene expression., At present, the improved high-throughput sequencing technology provides new insight into the functions of DNA methylation, and an increasing number of genome-scale methylome analyses have revealed that the relationship between DNA methylation and transcription is more nuanced than initially realized. For example, a genome-wide analysis of the interplay among DNA methylation, histone methylation and gene expression in rice indicated that ‘promoter region methylation may not affect transcription’, in contrast to previous studies. However, a more recent single-base resolution analysis of the rice methylome revealed that the promoter methylation negatively correlated to gene transcripts, but the negative correlation was only significant for these heavily-methylated genes. In contrast, gene body methylation was usually positively correlated with gene expression. Consistent with these findings, genome-scale analysis of DNA methylation variation associated with gene expression in Arabidopsis also verified the positive correlation between DNA methylation and gene expression and presented a more explicit point that DNA methylation only marginally contributed to gene transcripts. Correspondingly, our results also presented a complicated relationship between DNA methylation and gene expression and no obvious correlations between methylation and gene transcripts at the genome-wide level (Supplementary Fig. S3). Moreover, our results also showed that abundant gene expression alterations were not negatively related to the changes in DNA methylation (Supplementary Fig. S9). Similarly, methylome analysis of tomato fruit under chilling stress indicated that only 37–67% of different GO term genes presented a negative correlation between methylation level and gene expression. Genome-wide DNA methylation profiling of Arabidopsis plants exposed to a pathogen also demonstrated that abundant DEGs were not directly targeted by DNA methylation, but DNA methylation was at least partially responsible for the transcriptional control of these genes. Correspondingly, the exploration of possible causal relationships among genetic variation, methylation variation and expression variation in Arabidopsis also showed that DNA methylation could directly and indirectly regulate gene expression. Taken together, the poor correlation between methylation variations and expression alterations in this study were mainly caused by two reasons: (1) the regions in which methylation changes did not affect gene expression, like some gene body regions; (2) gene expression changes that were not directly controlled by DNA methylation.

Interestingly, carotenoid degradation was observed in 5azaC-treated callus. Further analysis revealed that the up-regulation of CpCCD1 contributed to carotenoid degradation in response to 5azaC treatment. We analyzed the methylation levels of the CpCCD1 promoter and exon regions by bisulfite sequencing and tried to link the phenotype to changes in methylation. However, the results indicated that the activation of CpCCD1 was not directly regulated by methylation changes (Fig. 10). As mentioned above, DNA methylation could indirectly regulate expression alterations. For example, chilling stress on tomato fruit increased the methylation levels of RIN (RIPENING INHIBITOR) TF and decreased the levels of RIN transcripts, resulting in the down-regulation of RIN binding target genes. This finding indicated that DNA methylation could affect transcriptional networks by regulating up-stream TFs. Similarly, Paper bagging altered DNA methylation and histone modifications on MdMYB1 and activated the MdMYB1 gene and its downstream genes, leading to differential anthocyanin pigmentation. The differential methylation of the MYB promoter also altered the MYB transcripts, resulting in differential anthocyanin accumulation in apple and pear.6, Excluding TFs, the indirect regulation of DNA methylation could be mediated by multiple mechanisms, like hormone pathways. Methylome analysis of Arabidopsis exposed to a pathogen revealed that DNA methylation activated transcripts of genes involved in the SA pathway and affected the expression levels of abundant SA downstream genes. Compared with the regulation of structural genes, it is understandable that plants prefer a more highly efficient and energy-saving way of controlling TFs or hormone pathways to further regulate biological processes. Thus, the activation of CpCCD1 was not directly controlled by methylation but was rather the consequence of methylation-dependent alterations in transcription networks. Concerning the mechanism of transcriptional regulation in CCDs, a recent study revealed that the TF SPL15 strongly and specifically interacted with the CCD gene promoter to regulate carotenoid accumulation in Arabidopsis. Moreover, transcripts of CCD1 were activated by hormonal treatments including ABA and strigolactone, and abiotic stresses including salt, drought and cold in Brassica rapa and Brassica oleracea. In this study, the RNA-seq data show that 5azaC treatment as a stimulant induced callus stress reactions (Fig. 5) and activated peroxidase gene expression and enzyme activities (Supplementary Fig. S7). These findings implied that the stress reaction caused by 5azaC could also contribute to the up-regulation of CpCCD1. Taken together, although our results indicated that DNA methylation was not directly responsible for CpCCD1 activation, it might be controlled by up-stream TFs or interaction factors (IFs), or hormonal pathways, among others, to affect CpCCD1 expression.

In total, a model was established to explain the effect of 5azaC on carotenoid degradation (Fig. 11). Transcriptional activation of CpCCD1 was potentially regulated by upstream TFs or IFs (interaction factors) that were induced by 5azaC treatment. However, 5azaC as a stimulator triggered the stress response in callus, and stress enhanced the transcriptional abundance of CpCCD1. The up-regulated CpCCD1 reacted with carotenoids and induced carotenoid degradation. The decrease in upstream carotenoid also blocked ABA biosynthesis. Moreover, the carotenoid degradation could also be associated with 5azaC stress. Correspondingly, 5azaC treatment stimulated callus peroxidation processes, and POD genes and enzyme activities were strongly activated to reduce the levels of peroxide. Furthermore, 5azaC induced the transcriptional profiles of ‘Phenlypropanoids’, ‘lignin and lignans’ and ‘Flavonoids’ pathways in response to 5azaC stress (Fig. 11).

Figure 11

A model was established to explain the effect of 5azaC on carotenoid degradation. The 5azaC treatment potentially activated upstream TFs or IFs (interaction factors) and then induced CpCCD1 expression levels. Additionally, 5azaC as a stimulator triggered the stress response in callus, and the stress enhanced the transcriptional abundance of CpCCD1. The up-regulated CpCCD1 induced carotenoid degradation. The decrease in carotenoids also blocked ABA biosynthesis. Moreover, carotenoid degradation could also be associated with stress induced by 5azaC treatment.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.