Genome-wide analysis of day/night DNA methylation differences in Populus nigra

DNA methylation is an important mechanism of epigenetic modification. Methylation changes during stress responses and developmental processes have been well studied; however, their role in plant adaptation to the day/night cycle is poorly understood. In this study, we detected global methylation patterns in leaves of the black poplar Populus nigra ‘N46’ at 8:00 and 24:00 by methylated DNA immunoprecipitation sequencing (MeDIP-seq). We found 10,027 and 10,242 genes to be methylated in the 8:00 and 24:00 samples, respectively. The methylated genes appeared to be involved in multiple biological processes, molecular functions, and cellular components, suggesting important roles for DNA methylation in poplar cells. Comparing the 8:00 and 24:00 samples, only 440 differentially methylated regions (DMRs) overlapped with genic regions, including 193 hyper- and 247 hypo-methylated DMRs, and may influence the expression of 137 downstream genes. Most hyper-methylated genes were associated with transferase activity, kinase activity, and phosphotransferase activity, whereas most hypo-methylated genes were associated with protein binding, ATP binding, and adenyl ribonucleotide binding, suggesting that different biological processes were activated during the day and night. Our results indicated that methylated genes were prevalent in the poplar genome, but that only a few of these participated in diurnal gene expression regulation.


Introduction
DNA methylation is a mechanism of epigenetic modification in eukaryotic organisms. In plants, DNA methylation generally occurs at cytosine bases in three different sequence contexts: CG, CHG and CHH contexts (H = A, T, or C) [1,2]. The mechanisms responsible for the establishment and maintenance of cytosine methylation have been best studied in Arabidopsis [3], where de novo methylation is catalyzed by Domains-Rearranged Methyltransferases2 (DRM2), and the maintenance of DNA methylation can be classified into two types: CG and CHG methylation are catalyzed by DNA Methyltransferase1 (MET1) and Chromomethylase3 (CMT3), respectively; CHH methylation is methylated by DRM2 [2,3] In both animals and plants, DNA methylation status might have effects on gene expression, splicing, and polyadenylation [4][5][6], and can influence DNA repair, recombination and meiotic crossover in euchromatic regions [7][8][9]. In plants, studies have indicated that DNA methylation plays an important role in parental imprinting, floral symmetry, flowering time, pigmentation, fruit ripening, sex determination, and stomatal development [10][11][12][13][14][15][16][17]. Therefore, variation in DNA-methylation may be functional and result in phenotypic variability in plants.
In recent years, increasing evidence has shown that DNA methylation is sensitive to extrinsic signals. For example, in the dandelion (Taraxacum officinale), alternative phenotypes may be induced by environmental stress due to DNA methylation changes [18]. Therefore, external environmental factors might influence gene expression through DNA methylation. Light can alter the promoter DNA methylation level in the suprachiasmatic nucleus (SCN) of the hypothalamus and reprogram circadian behavior in mice [19]. In plants, cold stress induced hypomethylation in Antirrhinum and growth temperature has been shown to positively correlate with the DNA methylation degree of Tam3 transposon [20]. In Arabidopsis elevated temperature have been shown to influence the expression of gene At3g50770 by DNA methylation changes in promoter region [21]. Additionally, different methylation levels were observed between plants grown in light and dark [22]. In the natural environment, plants are growing under constantly changing environmental conditions (light intensity, temperature and humidity etc.) over the course of a day, daily DNA methylation change are therefore expected in plant genomes.
Poplars (Populus L.) are model species in tree biology, and typically produce genetically identical clones quickly. Therefore, poplars are ideal plant materials for the study of DNA methylation changes under different environmental conditions. In a previous study, we didn't detect consistent DNA methylation variations over the course of a day in mature leaves of the P. nigra clone 'N46' by methylation-sensitive amplification polymorphism method (MSAP) [23]. The MSAP can only detect the methylation changes of the CCGG sites, but not the CHG and CHH sites which are prevalent in plant genome. In this study, using the same clone, we detected changes in DNA methylation of the whole genome between day and night using the methylated-DNA immune precipitation sequencing (MeDIP-Seq) method. We found that methylated genes were prevalent in the poplar genome, but that only a few of these participated in diurnal gene expression regulation. We will discuss the prevalence of DNA methylation in the poplar genome and its participation in diurnal regulation.

Ethics statement
The plants used in this study were propagated under permission from the state forestry administration of China. The location is not privately-owned.

Plant materials and growth conditions
We used the P. nigra clone 'N46' in this experiment. The production of 'N46'remets and growth conditions were as previously described [23]. After four months of culture in the greenhouse, four homogeneous plants selected and the 4 th -6 th leaves from the top were collected from each plant at 8:00 and 24:00, yielding two plants for each time point. Based on the photoperiod at the sample collection day (5:19~19:20), samples collected at 8:00 and 24:00 were used to present the day and night sample respectively. All leaves were immediately frozen in liquid nitrogen and stored at -80˚C.

Genomic DNA extraction and MeDIP-Seq
Genomic DNA was isolated from mature leaves of P. nigra clone 'N46' using the DNeasy Plant Mini Kit (Qiagen), and purified according to the manufacturer's instructions.
The MeDIP-seq DNA libraries were constructed according to the instruction manual for the MagMeDIP kit (Diagenode). Sonication was used to fragment genomic DNA to the size of 100-600 bp. A 3' single-base dA overhang on the end of each DNA fragment was ligated and then Illumina sequencing adapters were added to the ends by the Paired-End DNA Sample Prep Kit (Illumina). Denatured double-stranded DNA fragments were immunoprecipitated with 5-methylcytosine antibody (Diagenode). The quality of immunoprecipitated fragments were validated by quantitative real-time polymerase chain reaction (qPCR). DNA fragments of 200-300 bp were gel-excised and purified using the Qiagen gel extraction kit. The Quant-iT dsDNA High Sensitivity Assay Kit (Invitrogen) was used to quantify the extracted fragments on an Agilent 2100 Analyzer (Agilent Technologies, Germany). Then, DNA libraries were paired-end sequenced using the Illumina HiSeq 2000 platform (Illumina). The raw image files produced were analyzed by Illumina Real-Time Analysis (RTA) and base calling.
The adapters and low-quality reads were filtered out by Fastx (version: 0.0.13). The clean reads obtained were then aligned to the Populus trichocarpa reference genome (JGI 2.0.15). Genome mapping was performed using the Bowtie software (version: 0.12.8) in local alignment mode with default parameters to generate the BAM files.

Identification of differentially methylated regions (DMRs)
Methylation-enriched regions (or peaks) were detected from the BAM files using ChIP-seq (MACS, version 1.4), and peaks with a p-value less than 1 × 10 −5 were defined as passed peaks. The passed peaks were retained for further analyses.
DMRs were identified using the method developed by Wang et al. [24]. If the peak detected in one sample overlapped with peaks in another sample, the genomic region covering these was selected as a candidate DMR region. The log 2 ratio of 24:00 read numbers vs. those for 8:00 was calculated for each candidate DMR region, normalized, test p-value and using the DEGseq R package [25]. When p < 0.001 for the difference between the read numbers of each methylated region for two samples, we determined that the result was significant. The normalized log 2 value and p-value were combined to determine whether a region was differentially methylated. DMRs with a greater than twofold difference in read numbers were classified as hypo-or hyper-methylated regions.

Bisulfite sequencing of differentially methylated genes
To validate the MeDIP-Seq results, DMRs overlapping in five poplar genes were confirmed using bisulfite PCR (BS-PCR) and Sanger sequencing. Bisulfite conversions were completed using the DNA Bisulfite Conversion Kit (Tiangen). Bisulfite-modified PCR primers were designed using the online program BiSearch (http://bisearch.enzim.hu/?m=search). Information on the amplified methylated regions and primers are listed in Table 1. PCR was performed using GoTaq Hot Start Polymerase (Promega) in 20-μl reaction volumes containing 100 ng template DNA and 10 ng of primers. The PCR amplification conditions were: 7 min at 95˚C, followed by 40 cycles of 40 s at 95˚C, 40 s at 62˚C and 1 min at 72˚C, and further elongation at 72˚C for 5 min using an ABI Veriti 96-Well Thermal Cycler (Applied Biosystems). The PCR products were gel separated, purified with an AxyPrep DNA gel extraction kit (Axygen). Purified amplicons were cloned into the pGEM-T Easy vector (Promega) and sequenced. The Sanger sequencing data were analyzed using MethTools 2.0 (http://194.167.139.26/methtools/MethTools2_submit.html) and the methylation rate was calculated by dividing the number of non-converted (methylated) cytosines by the total number of cytosines in the amplified fragment.

RNA isolation and qPCR analyses
Total RNA was extracted from mature leaves using RNeasy Plant Mini Kit (Qiagen) for gene expression validation experiments. After the treatment of RNase-Free DNase (Qiagen), each RNA sample was quantified using NanoDrop ND-1000. The cDNA was synthesized from 2 ug total RNA using moloney murine leukemia virus reverse transcriptase with an RNase inhibitor (Promega), at an annealing temperature of 37˚C for 1 hour, then maintained at 65˚C for 10 min to terminate the reaction.
We performed qPCR analysis in the LightCycler 480 System platform (Roche, Switzerland). Each 20 μl PCR system contained 1 μl of first-strand cDNA, 200 nM of primers and 1× SYBR PCR mixture (TaKaRa, Japan). The amplification conditions were: 10 s at 95˚C, followed by 45 cycles of 10 s at 95˚C, 10 s at 60˚C and 10 s at 72˚C. We performed three replicates for each sample. Relative quantification values were calculated using the 2 -⊿⊿CT method [26]. Amplification lengths and primers are listed in Table 2.

Day/night methylation profiles in P. nigra leaves
The MeDIP-seq data were generated from samples collected at 8:00 and 24:00. Following the removal of low-quality data, we obtained 37,039,173 reads (average length 100 bp) from the
We detected a total of 74,565 methylation peaks from the 8:00 samples and 79,586 methylation peaks from the 24:00 samples (Table 3). More DNA methylation peaks indicate that more loci in the genome methylated [24]. These peaks were distributed across all 19 chromosomes, and the number of peaks in each chromosome was correlated with its length, indicating the prevalence of methylation in the poplar genome. Comparing the number of peaks in the samples for each time point, we found a slight increase in methylation level in the 24:00 samples.
Genomic regions that overlapped with DNA methylation peaks were considered methylated [24]. Differential distribution of the identified peaks among genomic components was observed, with most methylation peaks located in the intergenic regions (approximately 81.5%). Within the genic regions, the methylation peaks were most often located in the promoters (47.06%) (Fig 1). The methylation peak distribution patterns of P. nigra in genic regions were similar to those of Arabidopsis thaliana [30] and P. trichocarpa [27]. Based on Transposable Elements Database of P.trichocarpa [31], distribution of peaks in TEs regions was analyzed. The methylation peaks were seen in all the TE types, and the LTR_Gypsy, LTR_Copia and Hilitron tended to have much more methylated peaks than other TEs, suggesting that the activity control of these types of transposons may be more important in poplars (Fig 2).
Based on poplar genome database transcript annotation and the UniProtKB/TrEMBL database, the methylation peaks in the genic regions were function-annotated using the BLASTX program. The encoding protein of the methylated genes was further compared with proteins in the Kyoto Encyclopedia of Genes and Genomes (KEGG). In total, 10,027 and 10,242 methylated genes were annotated in the 8:00 and 24:00 samples, respectively. Gene ontology (GO) analysis of these methylated genes was performed by agriGO to detect the biological nature of the DNA methylation in the poplar genome. The methylated genes were annotated with up to 2,351and 2,443 GO terms by agriGO (http://bioinfo.cau.edu.cn/agriGO/), with 30.66% biological process terms, 21.02% cellular component terms, and 48.32% molecular function terms. The most enriched GO terms in the 8:00 and 24:00 samples were similar, such as protein phosphorylation and metabolic processes in the biological process GO domain; membrane, plasma membrane, chloroplast, and nucleus in the cellular component GO domain; and protein binding, ATP binding, nucleic acid binding, zinc ion binding, nucleotide binding, protein serine/ threonine kinase activity, protein kinase activity, and catalytic activity in the molecular function GO domain. GO term enrichment analysis showed that significantly enriched methylated genes (corrected p value < 0.05, FDR < 0.05) appeared to be involved in death, cell death, programmed cell death, apoptosis, and protein binding (Fig 3). KEGG analysis showed that these methylated genes were involved in 128-129 pathways, including metabolic pathways, biosynthesis of secondary metabolites, plant hormone signal transduction, penylpropanoid biosynthesis, plant-pathogen interactions, carbon metabolism, and fatty acid metabolism. These results indicate that DNA methylation was prevalent in the poplar genome and may influence the expression of thousands of genes involved in a wide range of biological processes. Methylation changed between day and night in P. nigra In the present study, a total of 2,120 DMRs were identified between the 24:00 and 8:00 samples. Among the identified 2,120 DMRs, 952 and 1,168 regions were classified as hyper-methylated and hypo-methylated, respectively. To better describe the corresponding biological function of these DMRs, we identified DMRs that overlapped with gene regions; 193 hyper-methylated DMRs and 247 hypo-methylated DMRs were obtained. The distribution of DMRs in genic regions was similar to that of the peaks (Fig 4). For TEs regions, about 90.27% DMRs located in the LTR_Gypsy, LTR_Copia and Hilitron transposons and no DMRs were seen in high levels of methylation were seen in LTR_ERVL, TIR_MuLE and TIR_TcMar (Fig 5).
In total, 137 genes overlapped with DMRs. The DMR-overlapped genes were annotated and subjected to GO analysis to detect the biological nature of the DNA methylation changes between day and night in the poplar genome. The DMR-overlapped genes were annotated to 56 GO terms by agriGO, with 4 in the cellular component, 26 in molecular function, and 26 in biological processes. No significantly enriched GO terms were found. The most enriched GO term in the cell component was the membrane; in molecular function the most enriched terms were transferase activity, kinase activity, phosphotransferase activity, ATP binding, protein binding, and adenyl ribonucleotide binding; in biological processes, they were apoptosis, cell death, protein amino acid phosphorylation, and post-translational protein modification (Fig 6).  Additionally, we found that the most enriched GO terms differed between hypo-methylated genes and hyper-methylated genes. Among the hyper-methylated genes, the most enriched GO terms were transferase activity, kinase activity, and phosphotransferase activity; among the hypo-methylated genes, they were protein binding, ATP binding, and adenyl ribonucleotide binding. These results suggest that different biological processes are activated during the day and night. The KEGG pathway analysis indicated that DMR-overlapped genes were mapped to carbohydrate metabolism processes, nucleotide metabolism processes, genetic information processing, environmental adaptation, and signal transduction pathways. Thus, although only a small amount of genes changed DNA methylation status between day and night, they were important in the diurnal response of poplar to the daily environment.

Verification of certain DMRs and their effects on gene expression
The accuracy of MeDIP-seq was validated by bisulfite sequencing PCR (BSP) using the DNA samples for MeDIP-seq. Five genes that overlapped with DMRs in promoter regions were selected, including Potri.012G138800 (chalcone synthase-like protein), Potri.007G044000 (Glutamate-gated kainate-type ion channel receptor subunit GluR5), Potri.002G180100 (metal tolerance protein, PTRMTP2), Potri.001G189900 (histidine phosphotransfer protein) and Potri.003G099100 (Nbs-lrr resistance protein) ( Table 1). The bisulfite sequencing data were consistent with the MeDIP-seq results (Table 4). For example, the MeDIP-seq data revealed that the methylation level in the 367-bp promoter region of Potri.012G138800 was hyper-methylated when 24:00 and 8:00 samples were compared. The bisulfite sequencing result showed that this region in the 24:00 samples contained 50 methylated cytocines ( m C), whereas the number of m C in 8:00 samples was 21, indicating that the 367-bp promoter region of Potri.012G138800 was hyper-methylated in the 24:00 samples, confirming the MeDIP-seq results.
The expression levels of the five genes with DMRs in the promoter region were detected using qPCR. The expression levels of two hyper-methylated genes (CHSL2 and the histidine phosphotransfer protein-coding gene) were significantly down-regulated in the 24:00 samples compared to the 8:00 samples. Conversely, the other three hypo-methylated genes exhibited the opposite expression patterns (Fig 7).

DNA methylation was prevalent in the poplar genome and only a small portion of genes may be involved in diurnal regulation
The characterization of genome-wide methylation patterns in plant systems, particularly DNA methylation profile changes in response to various environmental conditions, has been a hot research area in recent years [32]. Many studies have focused on the altered DNA methylation induced by dramatic environmental changes, such as abiotic stress and chemical treatment [33]. However, there is limited information regarding methylation changes in response to minor environmental changes, such as altered light intensity and temperature. The present study is the first to systematically compare genome-wide methylation profiles in P. nigra mature leaves between day and night. Considering the average day length during the experimental period, we chose 8:00 and 24:00 as the two sample collection time points, to represent DNA methylation during the day and night of one natural day. Using the MeDIP-seq method, we obtained DNA methylation profiles for both time points, and numerous genes were found  to be methylated in different genic regions, especially in the promoter regions, suggesting the potential role of DNA methylation in downstream gene regulation. These genes were shown to be involved in multiple biological processes, molecular functions, and cellular components, indicating the importance of the roles of DNA methylation in poplar cells. Additionally, 2,120 DMRs were identified by comparing the 24:00 and 8:00 samples, with only 440 DMRs located in the genic region that can influence the expression of overlapped genes. Therefore, although several genes were found to be methylated in the P. nigra genome, only a tiny fraction of these changed methylation status during the day/night cycle. In hybrid aspen (P. tremula × P. tremuloides), approximately 1,695 were identified as diurnally regulated genes [34]. In the current study, only 137 DMR-overlapped genes were identified by comparing the 24:00 and 8:00 samples, indicating that only a small proportion of diurnally expressed genes might be regulated by DNA methylation in poplars.
Previous studies indicated that plants have evolved a range of gene regulatory mechanisms to adapt to various environmental changes during the day, such as changes in temperature and light intensity [35]. Due to their functional role in transcriptional regulation, together with our RT-qPCR results, the methylation changes in the small amount of genes detected in this study did appear to cause time-specific expression of these genes between day and night, but the exact functional role of DMRs in the day night cycle till needs further study. Therefore, although DNA methylation is not a main regulation pathway in plant diurnal rhythm, it might be a minor compensation to other regulation pathways in plant diurnal rhythm.

The relationship between the methylation level and expression level was not linear
Previous studies have shown the relationship between methylation levels and corresponding gene expression in plants, with promoter methylation generally negatively correlated with gene expression, whereas gene body methylation positively correlated with gene expression [36,37]. In our study, BSP results of five genes with promoter regions overlapping with DMRs confirmed that the methylation levels in promoter regions repressed the expression of downstream genes. However, the effects of the DNA methylation levels of promoter regions on the expression of downstream genes were different. For example, the methylation level of promoter regions of the three genes (Nbs-lrr resistance protein coding gene, Potri.003G099100; metal tolerance protein coding gene, Potri.002G180100; and glutamate-gated kainate-type ion channel receptor subunit GluR5 coding gene, Potri.007G044000) decreased in 24:00 samples compared to those in the 8:00 samples. However, there was a significant difference in the expression levels of the three genes, with two genes (Potri.002G180100, Potri.007G044000) increasing by approximately 2.5-3.0 times, whereas one gene increased by only a small amount (Fig 5). Therefore, the DMRs in promoter regions have diverse effects on gene expression, and there is no direct linear relationship between the methylation level and expression level. This result may also explain those of several previous studies that showed only a small proportion of the DMR-overlapped genes exhibiting significantly different expression levels [38][39][40]. In most cases, the threshold to determine differentially expressed genes (DEGs) between samples was equal to or greater than a 1.5-fold change, and the expression differences in a large proportion of DMR-overlapped genes could not meet this criterion. We suggest that in such cases, the threshold for DEGs should be set lower because of the complicated relationship between DNA methylation and gene expression.