Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Unravelling the epigenomic interactions between parental inbreds resulting in an altered hybrid methylome in pigeonpea

Unravelling the epigenomic interactions between parental inbreds resulting in an altered hybrid... DNA methylation is an important heritable landmark conferring epigenetic changes in hybrids and has fascinated biologists and plant-breeders over the years. Although epigenetic changes have been documented in rice and maize hybrids, such investigations have not been reported in pigeonpea. Here, we report genome-wide methylation profiles of pigeonpea sterile and fertile inbred lines and their fertile F1 hybrid at single base resolution. We found that pigeonpea genome is relatively enriched in CG methylation. Identification of differentially methylated re- gions (DMRs) in the sterile and fertile parent revealed remarkable differences between their methylation patterns. Investigation of methylation status of parental DMRs in hybrid revealed non-additive methylation patterns resulting from trans-chromosomal methylation and trans- chromosomal demethylation events. Furthermore, we discovered several DMRs negatively associated with gene expression in the hybrid and fertile parent. Interestingly, many of those DMRs belonged to transposable elements and genes encoding pentatricopeptide repeats asso- ciated proteins, which may mediate a role in modulating the genes impacting pollen fertility. Overall, our findings provide an understanding of two parental epigenomes interacting to give rise to an altered methylome in pigeonpea hybrids, from genome-wide point of view. Key words: DNA methylation, TcM, TcDM, pigeonpea, hybrid 1. Introduction C-5 position of the cytosine ring by DNA methyltransferases and is Epigenetic regulation has remained of fundamental interest in under- transmitted across generations. DNA methylation acts to regulate standing gene regulation in eukaryotes. Epigenetic regulation in- transposons and repetitive genome sequences and often results in the 4,5 volves differential methylation at cytosine residues and post- suppression of gene activity. A major difference between DNA translational modifications of histones. One of the most comprehen- methylation in animals and plants is that in plant genome, DNA sively investigated mechanism of chromatin modification in plants, methylation occurs at CG, CHH and CHG sequence contexts (where DNA methylation involves covalent transfer of a methyl group to H could be C, T or A). It has been shown in Arabidopsis that V C The Author(s) 2018. Published by Oxford University Press on behalf of Kazusa DNA Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] 361 Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 362 The epigenomic interactions between parental inbreds different enzymes regulate methylation at different sequence con- variation in pigeonpea hybrids. This particular property has become texts. Methyltransferase 1 (MET1) causes and maintains methylation popular among breeders to produce varieties with improved fertility 6–8 at CG sites and chromomethylase 3 (CMT3) produces and main- and been utilized through cytoplasmic male sterility (CMs) system. 6,9,10 tains methylation at CHG sites and this maintenance is further CMS is a condition in which the ability of a plant to produce fertile required to preserve methylation patterns over generations. On the pollen is severely impaired. In this condition, rearrangements in mito- other hand, domains rearranged methyltransferase 2 (DRM2) and chondrial genome result in chimeric mitochondrial open reading 24 nt siRNAs, which guide methyltransfrase to specific DNA seg- frames (ORFs), generating non-functionality in anthers. However, the 6,10,11 ments, regulate de novo methylation CHH sites, which is main- expression of such chimeric mitochondrial ORFs can be suppressed by tained by CMT2. fertility restorer gene in F1 hybrid, thus restoring fertility. This kind of De novo methylation of all three different sites involves 24 nt hybrid breeding technology involves three parents— a male sterile par- smRNAs (small RNAs), DRM2 and the Dnmt3 homologue involved ent, its maintainer and a fertility restorer line (commonly called—A, B in RNA-directed DNA methylation (RdDM). In RdDM, RNA is and R lines, respectively; reviewed in detail in the literature). Few transcribed by the RNA polymerase POLIV, which is then made CMS associated mitochondrial orfs have been identified in double stranded by RNA-dependent RNA polymerases (RDRs), pigeonpea. cleaved to 24 nt by Dicer-like enzymes and loaded onto argonaute In this study, we have investigated whether there are any epige- proteins. This effector complex then interacts with other proteins netic changes in F1 progeny of cross between 11 A (male sterile) and and Polymerase V transcripts to target homologous regions of the ge- 303 R (fertile) inbred lines and whether those changes could poten- nome for methylation by the recruitment of DRM2. The methylation tially mediate an impact on the gene expression patterns including state of CG is maintained during cell replication by MET1, which those involved in fertility restoration in the hybrid. We conducted recognizes hemi methylated CG sites and methylate the newly whole genome bisulfite sequencing of 11 A, 303 R parental genotypes formed DNA strand. CMT3 maintains CHG methylation, while and their hybrid to generate whole genome cytosine methylation pro- CHH methylation is recognized by a histone methyltransferase files and found that overall, hybrid methylome was remarkably al- kryptonite (KYP), which can then methylate the H3K9 residues. tered compared with the parental methylome due to the events of These two enzymes work in a feedback loop enabling the mainte- transchrosomal methylation and demethylation in hybrid. We next nance of both epigenetic marks, with the absence of either CMT3 or carried out transcriptome profiling of the same plants and compared KYP leading to a reduction in CHG methylation. DNA methylation methylation status of the differentially expressed genes between pa- plays an important role during different aspects of plant development rental inbreds and the F1 hybrids. We found differences in methyla- and growth. Of these heterosis or hybrid vigour is of vital impor- tion patterns in the three genotypes, suggesting the involvement of tance in agriculture. Heterosis is a phenomenon where F1 progeny of DNA methylation in regulation of gene expression and possibly a genetic cross exhibits superior characteristics comparedwith its in- influencing sterility and fertility. bred parental lines and has been documented nearly 300 years ago. In recent years, several studies have focused on the role of epigenetic 2. Materials and methods regulation in mediating hybrid vigour in model plant 1,14–17 Arabidopsis. Previous studies in Arabidopsis have suggested 2.1. Plant material and growth conditions that interaction between similar genomes and different epigenomes Pigeonpea lines 11 A and 303 R were crossed by hand pollination of the two different parental ecotypes could result in new transcrip- using 303 R as pollen donor to generate hybrid. Seeds were collected tional profiles in hybrids that might contribute to hybrid vigour in and F1 seed resulting from the cross along with seeds from parental 18,19 the progeny. Also, it has been proposed that epigenetic regula- lines were sown again together for further experimental work. Seeds tion via changes in methylation patterns might play a role in mediat- were directly sown in the soil in green house under natural condi- ing differences in gene expression levels that may confer increased tions. All the plants were grown at the same time under exactly simi- biomass in heterotic hybrid offspring of Arabidopsis ecotypes C24 lar growth conditions to reduce the impact of any environmental and Landsberg erecta. It has also been documented that both DNA variations. Fertility status of all three genotypes was verified by pol- methylation patterns and sRNA (small RNA) levels are altered in in- len analysis. Pollens were fixed on a glass slide using 2% acetocar- 15,20–22 23 24 traspecific hybrids of Arabidopsis, rice and maize. mine. Pollens stained with acetocarmine were considered to be Recent advances in DNA methylation profiling techniques have en- fertile. After pollen verification, young leaves of same age from the abled in understanding the correlation between methylation patterns three lines were harvested, frozen in liquid nitrogen and used to ex- and hybrid vigour in several crop species as well such as rice and tract DNA and RNA for further work. 25,26 Brassica but knowledge about interaction of methylation patterns and gene expression in legume plants such as pigeonpea has still re- 2.2. Whole genome bisulfite sequencing library mained elusive. Pigeonpea (Cajanus cajan (L) Millsp) is a important construction and data analysis grain crop with high protein (20–22%) grain legume, mainly culti- High molecular weight genomic DNA was isolated from the leaves vated in the tropical and subtropical parts of the world. Pigeonpea using cetyl trimethylammonium bromide method. DNA was frag- contains 11 pairs of chromosomes (2n¼ 22) and has a genome size of 858 Mbp. Globally pigeonpea is cultivated on 5.4 Mha, with an an- mented into 100–300 bp using Covaris followed by end repair by nual production of 4.4 MT and a yield of 8299 hg/ha. Pigeonpea is adding ‘A’ and True seq methylated adapters were ligated. The li- grown in about 50 countries with India occupying 66% of its area gated DNA fragments were subjected to bisulfite treatment using EZ (FAO, 2016, http://www.fao.org/faostat/en/#data/QC). India is also DNA Methylation-Gold Kit, according to manufacturer’s instruc- the largest producer and consumer of pigeonpea with an annual pro- tions. The BS libraries were amplified, quality tested using Bio- duction of 2.51 MT. analyzer (Agilent 2100) and subjected to Ilumina Hi-seq2000 for 90 Pigeonpea is characterized by significant genetic variability. Genetic cycles in paired-end mode. The adaptors and poor quality reads were studies have also revealed the occurrence of non-additive genetic removed from the raw data using TrimGalore (https://github.com/ Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 363 FelixKrueger/TrimGalore). High quality reads were mapped to the Table S2) using Bowtie alignment program, then relative abundance 27,32 33 39 reference genome using BSMAP at default parameters. of transcripts was calculated using RSEM. Transcripts with zero Methylation information was extracted using Python script. Reads Fragments Per Kilobase of transcript per Million mapped reads were mapped on chloroplast genome to calculate methylation effi- (FPKM) value were excluded from analysis. Raw reads count data ciency and error rate. A bisulfite conversion rate of >99% was ob- were imported into R environment and data table containing fragment tained for all three samples. In-order to determine true methylated counts for each transcript were combined and matrix was formed to cytosines, P-values were calculated using binomial distribution identify differentially expressed gene using edgeR. P¼ binomial (m, x, error rate), where m¼ number of methylated reads, n¼ number of non-methylated reads and x¼ m þ n (coverage 2.6. GO enrichment analysis on a single C). Cytosines with P-value< 0.0025 were considered to To identify and categorize cluster of differentially expressed genes be truly methylated. Fraction and mean methylation level was calcu- (DEGs) and genes enriched for methylation in all the three categories lated as described previously. (11 A vs. 303 R, 11 A vs. hybrid and 303 R vs. hybrid), GO enrich- ment analysis was performed using python based program 2.3. Identification of differentially methylated regions GOATOOLS (https://github.com/tanghaibao/goatools) and Wego. To identify differentially methylated regions (DMRs) at each se- Genes with Bonferroni-corrected P-values< 0.05 using GOATOOLS quence context (CG, CHG and CHH), genomes were tiled using a were considered to be enriched. sliding window approach with 100 bp window sliding at 100 bp in- terval for CG and CHG and 150 bp for CHH. Methylation informa- 2.7. Identification of transposable elements tion for each window was calculated using Methyl kit. A denovo repeat library was generated using repeatmodeller (http:// P-value for each window was calculated using Fisher’s exact test. www.repeatmasker.org/RepeatModeler/; version 1.0.11).This cus- To determine q-values, P-values for each window were corrected us- tomized library was used to run RepeatMasker (http://www.repeat ing sliding linear model . Following criteria were used for identifica- masker.org/; version- 4.0.7) to identify and annotate transposable el- tion of DMRs. ements (TEs) in pigeonpea at default parameters. TE shorter than 100 bp were excluded from the analysis and the distance between 1. Windows contained at-least three methylated cytosines for CG TE, less than 50 bp were conconated (Supplementary Table S3). and CHG; and six for CHH DMRs. 2. Bases having more than 10X coverage were considered. 3. Bases having more than 99.9th percentile coverage were discarded 2.8. Allele specific methylation analysis using single to minimize errors arising due to PCR bias. nucleotide polymorphism 4. Windows with methylation difference greater than 20% for CG, Single nucleotide polymorphisms (SNPs) from parental lines and hy- CHG and 25% for CHH were selected. 42 brid were identified using BIS–SNP, a variant caller from BS reads. 5. Fisher’s exact test corrected q-values < 0.01. For this, the mapped methylated reads were used as input files along with pigeonpea reference genome. SNPs detection was performed 2.4. Validation of methylation levels using bisulfite using the following parameters: ‘-T BisulfiteGenotyper -stand_call_- sequencing PCR conf 20 -stand_emit_conf 0 -mmq 30 -mbq 17 -minConv 20’. Spurious SNP with quality score less than 20, more than 120x read Two different regions were chosen one for validating methylation coverage, strand bias more than—0.02, quality score by depth less levels, one in the intergenic region and the other one in the exonic re- than 1.0, mapping quality zero reads, fraction more than 0.1 and 2 gion. For bisulfite PCR, genomic DNA was extracted from 11 A SNPs within the 20 bp bin were further filtered out using VCF post and F1 hybrid and subjected to bisulfite conversion using EZ TM process command. DNA Methylation-Gold Kit. Primers were designed using METHPRIME (Supplementary Table S1). Bisulfite converted DNA was amplified and the fragments were cloned into pGEMT vector. 3. Results Ten colonies were picked randomly for sequencing. The results were analyzed using QUMA. 3.1. Genome wide methylation analysis in pigeonpea To understand the possible roles of epigenetic regulation and the dy- 2.5. RNA sequencing library construction and analysis namics of DNA methylation during hybrid formation in pigeonpea, we crossed 303 R (Fertile, pollen donor) and 11 A(sterile) lines. Since Total RNA was extracted from the same plants used for BS Seq. Equal their F1 progeny restored fertility (Fig. 1A), we hypothesised that it amount of total RNA from three biological replicates were pooled. About 1 mg of the total pooled RNA was used to purify poly-A con- could possibly partially attributed to the new epigenome resulting taining mRNA molecules. The purified mRNA was fragmented into from the interaction between the two parental epigenomes. Thus, in- smaller pieces and reverse transcribed into cDNA. The cDNA frag- order to gain deep insights into the epigenome of the parental lines ments were subjected to end repair by adding single ‘A’ base. True seq and their F1 hybrid, we deciphered their respective whole genome adaptors were ligated and cDNA libraries were amplified for 15 cy- methylation pattern at single base resolution. For this high molecular weight genomic DNA was isolated from the leaves of 10 days after cles. Quality of the library was assessed using Bioanalyzer and the vali- flowering (DAF) old plants from three different genotypes, namely dated library was sequenced by Illumina Hi-Seq 1000 according to 11 A, 303 R and F1 hybrid progeny of 11 A and 303 R. The DNA manufacturers’ instructions. The quality of the raw reads was deter- samples were subjected to bisulfite treatment to generate genomic li- mined using FastQc (https://www.illumina.com/products/by-type/infor matics-products/basespace-sequence-hub/apps/fastqc.html). Adaptors braries for high throughput bisulfite sequencing. The reads obtained and poor quality reads were cleaned using Trim Galore. High quality were processed for removal of adapter and low quality sequences. paired end reads were aligned to reference genome (Supplementary About 80–100 million cleaned reads (in each of the three genotypes) Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 364 The epigenomic interactions between parental inbreds Figure 1. DNA methylation profiles of pigeonpea inbred lines (11 A, 303 R) and their hybrid. (A) Acetocarmine stained pollens of 11A (sterile parent), 303R (fertile parent) and fertile hybrid. Note the density of pollen is less in 11A compared to 303R and hybrid. Also note viable pollens in 303R and hybrid.(B) Genome-wide relative proportions of CG (blue), CHG (red) and CHH (green) methyl cytosine in 11A, 303R and F1 hybrid. (C and D) Fractional methylation levels in the three lines of pigeonpea in CG, CHG and CHH contexts (B) and comparisons with published reports in Soybean, Rice, Populus, Arabidopsis and Brassica (D). (E) A cir- cos plot representing chromosome wise distribution of mCG, mCHG, mCHH, gene, transposons density and FPKM (outer to inner) in F1 hybrid. Yellow to red color relates to lower to higher methylation density in a 100 Kb window (Colored figure is available in the online version of the article). Scale bars 50mm (A). 27,32 were mapped on to the reference genome. We found that more cytosine methylation at single base resolution. Further, to calculate than 80% reads got mapped on to the reference genome. Thus even- the frequency of cytosine conversion, reads were mapped to un- tually, 46 408 326, 54 030 096 and 43 763 469 reads were obtained methylated chloroplast genome. We obtained a conversion rate for 11 A line, 303 R line and their F1 hybrid progeny, respectively, of more than 99% in all the three genotypes, indicating high effi- that mapped uniquely to the reference genome and providing a cov- ciency of cytosine conversion in our experiments (Supplementary erage more than 20X. These were subjected to further analysis of Table S4). Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 365 We then deciphered the fractional proportion of DNA methyla- 3.2. Patterns of DNA methylation in genic and tion and found that CG methylation was highest in 11 A (41.7%) transposable elements while CHH methylation was highest in F1 hybrid (28.82%). To precisely understand the methylation patterns of the protein cod- Methylation in CHG context were similar in all the three genotypes ing genes, genes with inserted TE elements were excluded from (34%) (Fig. 1B). Highest mean methylation was observed in CG analysis. We calculated the methylation levels in the genic region and context (74%-79%), closely followed by CHG (56%–63%) while 2 kb upstream and downstream regions of the gene bodies in the two mean methylation in CHH context (10–15%) was lowest in both the parental lines and their hybrid and found significant methylation parents and the hybrid (Fig. 1C). Comparative analysis of fractional differences between the three lines (Wilcoxon rank sum test, P-val- methylation levels in pigeonpea, soybean, brassica, oryza and ue< 2.2e-16, n ¼ 26,774). These differences reflect differential meth- populus revealed that pigeonpea genome was hypermethylated com- ylation patterns in the genic regions of parental lines and their pared with the other crops. Interestingly, CHH methylation levels in hybrid. Furthermore, upon a detailed investigation of methylation pigeonpea genome were significantly higher than those observed in pattern, we observed the occurrence of a single CG methylation peak 25,43 soybean, brassica, oryza and populus (Fig. 1D). It has been doc- within the gene body and two troughs near transcriptional start and umented that the DNA repeats in the genomes are highly methylated termination sites in pigeonpea, invariably to earlier reports in other 9,44,45 43,48,49 48 in symmetric contexts. Pigeonpea geneome is highly repetitive plants. However, contrary to a recent report in rice, we 46 47 (47%) compared with soyabean and Arabidopsis. Thus, such found over all low CG methylation levels within the coding regions high levels of DNA methylation in pigeonpea may arise from the occ- compared with their flanking regions. A similar trend was observed curance highly repetive genome. for CHG methylation, however, the average CHG methylation both Visualization of chromosome wide distribution of CG, CHG and within the coding regions and their flanking regions were compara- CHH methylation revealed higher CG and CHG methylation levels tively lower than CG methylation. CHH methylation however compared with CHH while relatively even distribution pattern of showed a different trend to both CG and CHG methylation in terms CHH methylation (Fig. 1E, Supplementary Fig. S1A and B). High of a plateau instead of a single peak within the gene body. Overall CG and CHG methylation were correlated with transposon rich re- CHH methylation both within the genic regions and their flanking gions across the genome (Fig. 1E, Supplementary Fig. S1A and B). regions were significantly attenuated than CG and CHG Figure 2. Distribution of methyl cytosines in genes and transposable elements in pigeonpea. (A) Line plots showing methylation level in 2Kb flanking (upstream and downstream) and within the gene body in 11A, 303R and F1 hybrid in CG, CHG and CHH contexts. (B–D) Line plots showing methylation level in 2Kb flank- ing (upstream and downstream) and within the gene body of transposons- long terminal repeats (B); LINE (C) and DNA Elements (D) in 11A, 303R and F1 hybrid in CG, CHG and CHH context. Note higher mCs throughout the transposons compared to genes (Colored figure is available in the online version of the article). Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 366 The epigenomic interactions between parental inbreds methylation. However, the hybrid exhibited mid parental range for different criteria were applied for CHH DMR identification. To as- CG and CHG but a higher CHH methylation compared with paren- sign a DMR to CHH category, each window had at-least eight cyto- tal lines, suggesting altered CHH methylation patterns in the genic sines with minimum of 10X coverage and a minimum of 25% regions in pigeonpea hybrids (Fig. 2A). difference in CHH methylation. Eventually, we identified much more Eukaryote genomes are abundant in TEs, also called mobile ele- CG DMRs (11,300) than CHG (6988) and CHH (8061) DMRs be- ments. TEs are structurally divergent, epigenetically controlled and are tween 303 R and 11 A parental lines (Supplementary Table S5). We mostly involved in genome evolution and gene activity regulation. also found that the differential methylation between the two parents Thus, we sought to gain a comprehensive insight of the influence of was more enriched in upstream flanking sequences of the genes than DNA methylation on the transposable elements within the pigeonpea the down stream regions (Fig. 3A). We next assessed relative methyl- genome. Since no complete dataset on pigeonpea TEs is available, we ation percentages of the obtained DMRs between the 11 A and used repeat masker to identify and classify them. TEs shorter than 303 R parents (11 A vs. 303 R or Sterile vs. fertile) and found enrich- 100 bp were excluded from analysis. In total, we found 105 733 TEs, ment for hypermethylated CHH DMRs and hypomethylated CG which we classified into retro transposons and DNA elements. We DMRs (Fig. 3B). This suggests that in 11 A line, genes involved in next calculated their CG, CHG and CHH methylation levels in all the regulation of CG methylation are less active than the RdDM path- three genotypes and found that irrespective of their categories, methyl- way genes and CMT2 that are involved in the regulation CHH meth- ation was highly enriched in all the contexts in TEs compared with the ylation (Supplementary Fig. S2). gene body. However, interestingly, we found relatively lower CG and Interaction of two dissimilar parental epigenomes can result in an CHG methylation levels in the TEs flanking regions than those of the altered methylome in their hybrid progeny. Since fertility is restored genic regions. By comparing the methylation levels in long terminal in the F1 hybrid and earlier studies in Arabidopsis showed higher 15,20 repeats within the three genotypes, we found relatively higher methyl- methylation levels in hybrid compared with parental lines, we ation levels in F1 hybrid in all the three contexts, compared with the therefore investigated whether there is differential methylation be- parental lines except for CG methylation in flanking regions where tween individual parental lines and the fertile hybrid (Supplementary 11 A exhibited higher methylation levels (Fig. 2B). While, in the case Table S5). Consistent with previous studies, we found a relative of LINE elements, 11 A displayed hyper methylation compared with increase in the percentage of methylated cytosines in CG, CHG and the F1 hybrid and 303 R in CG and CHG contexts (Fig. 2C). DNA CHH DMRs in the hybrid compared with parental lines (Fig. 3B) in- elements were relatively hypermethylated than the class 1 transposons. dicative of occurrence of non-additive methylation in pigeonpea hy- F1 hybrid displayed uniformly higher methylation irrespective of the brids. We randomly chose hyper and hypo methylated DMRs mC contexts (Fig. 2D). Overall we observed comparatively lower between sterile parent and F1 hybrid in intergenic region and coding methylation levels in 303 R, compared with 11 A and F1 hybrid regions. Analysis of their methylation levels using bisulfite sequenc- (Fig. 2B–D), suggesting differences in the regulatory mechanisms of ing PCR revealed relatively higher methylation levels in intergenic re- TE silencing via methylation in pigeonpea. gion associated DMR than in the coding region associated DMR (Supplementary Fig. S3) consistent with our in-silico data. In-order to assess the functional categories of the genes associated 3.3. Fertile and sterile pigeonpea inbred lines display with the DMRs between parents and hybrids and to investigate locus specific differences in cytosine methylation whether DMRs were associated with genes involved in DNA methyla- Earlier work in crop plants such as rice and maize have documented tion, epigenetic processes or fertility regulation, gene ontology enrich- differences in methylation levels between a male sterile plant and a ment analysis was performed. We found genes ontology terms for fertile plant and revealed hyper-methylation in sterile lines compared regulation of flowering time, DNA methylation, regulation of gene si- 50,51 with fertile lines. This suggests differential cytosine methylation lencing by miRNA, developmental processes, Histone H4 K4 methyla- could possibly influence fertility and sterility in plants by regulating tion, DNA methyltransferase, transcription factors, proton gene expression. Therefore to investigate whether pigeonpea 11 A transporting ATP synthase activity, regulation of post transcriptional and 303 R inbreds exhibit differences in their methylation levels, we gene silencing were enriched (Fig. 3C). Interestingly, we also observed compared the methylome of the two lines. Several studies have docu- a strong enrichment of DMRs in gene is associated with aldehyde de- mented that viewing differential methylation in windows of multiple hydrogenase activity. A previous study have reported that RF2 cytosines seems to reflect more functional consequences rather than encoded aldehyde dehydrogenase activity is required to restore male 52,53 viewing at single sites. Also it is known that CG, CHG and fertility in maize and also plays a role in anther development. Our CHH methylation levels are maintained through different regulatory results thus suggest a possible correlation between differential methyl- mechanisms. Thus, to identify DMRs between 11 A and 303 R lines, ation and fertility or sterility in 11 A, 303 R pigeonpea inbred and we decided to investigate and categorize differential methylation for their hybrid. Also, genes involved in cellular processes, biosynthetic all the three contexts, CG, CHG and CHH, separately. We averaged processes and ATP binding were significantly enriched (Bonferroni the DNA methylation for all the three contexts independently across corrected P-value< 0.05; Fig. 3C, Supplementary data S1). non-overlapping 100 bp regions (windows), as done previously. Overall, our results indicate that differential methylation could be In-order to categorize a DMR window as CG, CGH or CHH, differ- involved in mediating and regulating several cellular and develop- ent conditions were applied. To assign a DMR to CG or CHG cate- mental processes in pigeonpea. gory, windows having at-least three cytosines with a minimum of 10X coverage were used. A minimum difference of 20% in CG and 3.4. Trans-chromosomal methylation events alter in CHG methylation each was chosen. Previously, it has been re- pigeonpea hybrid methylome ported that CHH methylation patterns appeared to be more variable than CG and CHG. Since CHH methylation sites are much abun- Our results so far show altered methylome in pigeonpea hybrid com- dant in genomes and still, the CHH methylation levels were observed pared with parents. The changes in hybrid methylome occur at to be much lower than the CG and CHG methylation levels (Fig. 1), DMRs between the parents and are mediated by either increase or Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 367 Figure 3. Comparative analyses of DMRs between the inbred lines and their hybrid. (A) Pie charts showing the distribution of DMRs in the different genomic re- gions in 11A vs. 303R, 11A vs. hybrid and 303R vs. hybrid. (B) A bar plot showing number of hypo and hypermethylated CG, CHG and CHH DMRs in 11A vs. 303R, 11A vs. hybrid and 303R vs. hybrid. (C) A heat map showing enrichment of DMR associated genes in 11A vs. 303R, 11A vs. hybrid and 303R vs. hybrid. Color code represents low to high P-values (0.05 to 8.09e-10) (Colored figure is available in the online version of the article). decrease in DNA methylation termed as non-additive methyla- of 303 R fertile parent and hybrid were much higher than the 11 A 15,20 tion. This increase or decrease occurs via altering methylation sterile parent. In contrast, CHH methylation levels of 303 R fertile parent and hybrid were lower than that of the 11 A sterile parent status of one parental allele to resemble that of the other parental al- lele, events known as trans-chromosomal methylation (TCM) and (Fig. 4A). Additionally, to decipher the number of direct TcM and trans-chromosomal demethylation (TCdM). We investigated non- TcDM non additive mC clusters, we used an alternative approach of SNPs as described previously. For this, we tracked 75,340 mC clus- additive methylation in F1 hybrid to understand its methylome dy- ters of parental origin and followed the methylation pattern of pa- namics. To understand this mechanism, DMRs between the two pa- rental lines were identified and mapped onto the F1 hybrid genome. rental alleles with differences in methylation level (3,222 out of 75,340) in hybrid using SNPs. We observed 27% mC clusters were Average methylation of the mapped loci was calculated in F1 hybrid. stably inherited in hybrids indicative of additive methylation. Among Mean parental value (MPV) of the DMRs between parents were the clusters showing non additive methylation, we observed, TCM obtained and compared with observed methylation in hybrid. Loci events in 18.5% cases where methylation of epi-alleles from low having methylation difference of atleast 20% were considered as parent increased in the hybrid, whereas methylation of high parent non-additive methylation sites. Subsequently, we found non-additive epi-alleles remained unchanged in the hybrid, resulting in an over all methylation in 28% DMRs suggesting methylation interaction at increase in methylation in the hybrid. While in 24.1% clusters, these regions in F1 hybrid. We further analyzed the nature of this TCdM was observed where methylation of the epi-allele from low non-additive methylation and observed patterns of Trans- parent remained unchanged but methylation of the epi-allele from chromosomal methylation at 3,911 DMRs and Trans-chromosomal high parent decreased in hybrid, resulting in an overall decrease in demethylation in 5,377 DMRs (Fig. 4A and B). By investigating the methylation in hybrid (Fig. 4C). genomic features of the DMRs showing non-additive methylation, we found a higher proportion resided in TE (25%) as compared with the genic regions (15%). These results are in accord with previous 3.5. Transcriptome profiling in pigeonpea findings revealing the roles of TE in genome stability and their regu- In order to assess the transcript abundance and differential gene ex- lation by methylation. Furthermore we found that in all the three pression between the three genotypes, their whole transcriptome pro- contexts, the methylome of F1 hybrid resembled more of the 303 R filing was performed by RNA sequencing. We first compared the fertile parent than the 11 A sterile. CG and CHG methylation levels transcription between different genotypes (11 A vs. 303 R, 11 A vs. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 368 The epigenomic interactions between parental inbreds Figure 4. (A) Heat maps showing the methylation status of DMRs of the parents and methylation pattern at those loci in the F1 hybrid. Note in all the mCs con- text, methylation patterns of F1 hybrid are similar to the fertile parent (303 R) and different to the sterile parent. (B) Bar plot showing non-additive methylation in F1 hybrid. Note increased average methylation in hybrid than the MPV in the upper graph indicating the occurrence of TCM. Similarly, note a decreased aver- age methylation in hybrid than the MPV in the lower upper graph indicating TCdM. (C) Bar plots showing changes in parental allelic methylations in the hybrid. Majority of the non-additive methylation occurred through TCdM (Blue). Red bar represent TcM. Methylation levels of most parental alleles remained unchanged in hybrid (gray) (Colored figure is available in the online version of the article). HP: high parent; LP: low parent. F1 hybrid and 303 R vs. F1 hybrid) to compute genes that were dif- that gene ontology terms included reproductive processes, develop- ferentially expressed between these lines. The transcripts showing a mental processes, protein–DNA complexes and nucleic acid binding minimum of 2-fold change and P-value<0.05 were considered dif- proteins that were significantly enriched (Fisher’s exact test P-val- ferentially expressed. We found 1,621 (460 up, 1,161 down), 1,790 ue<0.05) in all the three comparative scenarios 11 A vs. 303 R, 11 A (1,052 up, 738 down) and 1,380 (963 up, 417 down) genes that vs. hybrid and 303 R vs. hybrid. showed differential expression in 11 A vs. 303 R, 11 A vs. hybrid and 303 R vs. hybrid, respectively. We found that in all the three cases of 3.6. Correlation between DMRs and gene expression in comparison, the differentially expressed genes were well represented pigeonpea inbred and hybrid line in gene categories like pentatricopeptide repeats (PPRs)- associated proteins involved in fertility restoration, methylases, demethylases We next asked if the differential methylation in the three different ge- notypes could potentially influence the gene expression in these and histone modifying genes (that play a role in epigenetic regula- tion) and flowering-associated genes. We performed gene ontology plants. For this, first DEGs with DMRs localized in 2 kb flanking (up- stream and downstream) and gene body regions were identified as enrichment analysis to identify functional categories of significantly enriched genes (Supplementary Figs S4–S6). Interestingly, we found DMR-associated DEGs. Then the correlation between whole genome Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 369 differential methylation and gene expression was examined investigated the correlation between methylation status and expres- (Fig. 5A–C). We next determined the genes that were enriched for sion in the DEGs associated with DMRs and found that 127 (16 DMRs between the three genotypes and found that DMRs were local- hypomethylated and up-regulated; 111 hyper-methylated and down- ized in 3,019, 1,996 and 2,233 genes in 11 A vs. hybrid, 11 A vs. regulated), 75 (29 hypo-methylated and up-regulated; 46 hyperme- 303 R and 303 R vs. hybrid, respectively (Fig. 5A–C). Among these, thylated and down-regulated) and 229 (197 hypo-methylated and up- 290, 586 and 578 genes were found to be both differentially ex- regulated; 22 hyper-methylated and down-regulated) genes showed a pressed and associated with DMRs in 11 A vs. F, 11 A vs. 303 R and negative correlation between their expression and methylation in 11 303 R vs. F, respectively (Fig. 5A–C). Furthermore, we also A vs. R, 11 A vs. Hybrid and 303 R vs. Hybrid, respectively (Fig. 5D Figure 5. Characterization of DMRs associated DEGs in pigeonpea. (A–C) Venn diagrams showing unique and common DEGs (differentially expressed genes) to the genes associated with DMRs in 11A vs. hybrid (A), 11A vs. 303R (B) and 303R vs. hybrid (C). (D) Bar plots showing number of DEGs that were differentially methylated in 11A vs. 303R, 11A vs. hybrid and 303R vs. hybrid. (E–G) Heat maps displaying comparative methylation levels of selective genes that showed neg- ative correlation between their expression and methylation in 11A vs. 303R (E), 11A vs. hybrid (F) and 303R vs. hybrid (G). Color code indicates methylation level (0–1) (Colored figure is available in the online version of the article). Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 370 The epigenomic interactions between parental inbreds and Supplementary data S2). In order to gain insights into the biologi- and downstream regulatory sequences of the genes, we also sought cal significance of the observed negative association between DEGs to examine methylation levels in the gene flanking sequences. Similar and DMRs, we investigated the functional categories of those genes. methylation levels were observed in the sequences upstream and Interestingly we found many DMRs showing negative association downstream of genes with high, moderate and low expression. with gene expression in fertile hybrid and 303 R lines, compared with However, contrary to this, genes with no expression exhibited high 11 A (sterile) line belonged to transposable elements and genes that methylation in their flanking regions suggesting that a high methyla- coded for PPR-associated proteins (Fig. 5E–G). In addition, we also tion in upstream and downstream regulatory sequences might be found genes coding for histone deacytylase and methyltransferase critical for suppressing transcript levels (Fig. 6A and Supplementary (Fig. 5E–G). This suggests that the transcription of those genes might Fig. S7). We also found that the presence of TE elements within the possibly be influenced by methylation. genes or in the vicinity resulted in a comparatively reduced expres- sion relative to the genes without TE elements within or in the neigh- boring regions (Fig. 6B). 3.7. Association between DNA methylation and transcript abundance Previous studies have reported a negative association between gene 4. Discussion expression and CHG and CHH methylation in plants like cotton. Rapid advancements in genome sequencing techniques have enabled Therefore to investigate the occurrence of any relation between genomic comparisons to understand genetic diversities between differ- changes in methylation levels and gene expression in pigeonpea, we ent plant species or different ecotypes of the same species. However, classified genes into four quartiles based on their FPKM values— understanding of epigenetic variation among different ecotypes has be- undetectable expression (FPKM 0-1), low expression (FPKM 1-10), come equally important to decipher its role in generating phenotypic moderate expression (FPKM 10–100) and high expression variations. A fundamental question of great importance is that how (FPKM> 100) according to the criteria described in earlier stud- 58,59 ies and correlated with changes in their methylation levels. the parental genome combination is regulated in hybrids to give rise to Consistent with previous studies, we found CHG and CHH methyla- significantly different hybrid genome which results in altered genetic traits. It is widely speculated that many agriculturally important phe- tion were negatively correlated with gene expression levels. nomena such as hybrid vigour and fertility restoration may arise due However, additionally, we also found negative association between CG methylation and gene expression (Fig. 6A and Supplementary to epigenetic changes between inbred parents and hybrid lines. Thus Fig. S7). Since transcription is importantly controlled by upstream the field of epigenetic research has expanded its horizon by adopting Figure 6. Correlation between gene expression and DMRs in pigeonpea. (A) Lineplots showing correlation between CG, CHG and CHH DMRs and the gene expres- sion which was divided into four quartiles—no expression (purple); low expression (green); moderate expression (red) and high expression (cyan) in F1 hybrid. (Also see Figure S7 for correlation in 11A and 303R parental lines). Note a negative correlation between gene expression and methylation in all the three contexts. (B) Box plots showing expression levels of genes inserted with and without TE within their gene bodies and 2kb upstream and downstream flanking regions. Gene þTE indicates insertion of transposon. Gene-TE indicates no insertion of transposon (Colored figure is available in the online version of the article). Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 371 many crop species as model systems. However, there are still agricul- Transposable elements can influence neighboring gene expression, turally important crops such as pigeonpea, in which existence of such contribute to inheritance of epi-alleles and thus represent the most mechanisms have remained unexplored. In this study, we investigated striking and common example of epigenetic inheritance in plants. and compared epigenetic changes in three different pigeonpea geno- Consistently, our data also revealed low expression of genes with in- types. We have generated whole genome methylation profiles of sterile serted TEs suggesting similar roles of TEs in pigeonpea. (11 A), fertile parents (303 R) and their fertile F1 hybrid, which were We also found a strong enrichment of DMRs in a wide spectrum unavailable till date. Our analysis of relative fractions of methyl- of genes that are associated with chromatin modification and have cytosines in each methylation context in the three pigeonpea lines re- been reported to regulate flowering, heterosis and reproduction in vealed highest fraction of mCs in CG context, followed by CHG and rice. We also found genes involved in DNA methylation that are CHH. Similar pattern of distribution of methylation fractions has known to play a role in transposon repression and genes involved 60 61 been observed previously in rice seedlings and Soybean leaves. We in histone methylation that help regulate flowering. also found that average methylation levels of the fertile hybrid and the DNA methylation is known to influence gene expression. fertile parental line (303 R) were lower than the sterile parent (11 A) However, as mentioned above, it has been proposed that DNA meth- suggesting that comparative low DNA methylation of genes may play ylation and H3K9 methylation in the gene body regions do not affect a role in restoring the fertility in the F1 hybrid. gene abundance. Despite increased CHG and H3K9 methylation in Mean methylation has been reported to be highest in CG contexts loss of function mutants of a H3K9 demethylase gene increased in 60 15 in other plants such as Rice and Arabidopsis, and followed bonsai methylation 1 (ibm1), the expression of their target genes re- closely by CHG and CHH. However, our results show that in mained unaffected. Also, absence of CG gene body methylation in pigeonpea, mean methylation levels in CHH context are much lower the loss of function MET1 mutants did not alter gene expression. than CG and CHG methylation levels. This suggests a possibility of Surprisingly, we found a negative correlation between transcript overall low expression of genes involved in regulation of CHH meth- abundance and methylation levels in all the three contexts in the ylations such as those involved in RdDM pathway that target homol- gene body suggesting a more prominent role of DNA methylation in ogous regions of the genome for methylation in pigeonpea regulating gene expression in pigeonpea. Previous studies in compared with Rice and Arabidopsis. Arabidopsis have reported that methylation in the promoter regions 71,72 Previously, it has been reported that transcription of genes is af- resulted in lower gene expression. Furthermore, another study in fected by methylation in the promoter region but not in the internal soybean revealed negative correlation between gene expression and regions. However, transposon silencing is regulated by methylation DNA methylation in their promoters. We also found that high both in the promoter region as well as internal regions. Our data methylation in the 5’ upstream regions of the genes were negatively demonstrate that fraction of mCs was higher in the transcribed re- correlated with gene expression levels, consistent with the idea that gions of transposons than the transcribed regions of the genes. As DNA methylation in the promoter regions affects gene abundance. shown in Figure 2, the fraction of methyl cytosines in upstream se- Our results in pigeonpea together with previous results in soybean quences of the transcribed regions of transposon was almost two open new avenues for understanding the relationship between gene times higher than fraction of mCs in upstream sequences of genes. expression and gene body DNA methylation and whether they play Additionally, gene bodies of TEs exhibited higher methylation than a direct role in restoring fertility in pigeonpea hybrids requires fur- their upstream and downstream regions. Contrary to this, genes had ther investigation. lower methylation in the transcribed regions than their flanking re- The raw reads of bisulfite and transcriptomic study were submit- 1,64 gions. As suggested previously for other crops, this could poten- ted to NCBI Sequence Read Archive under bioproject accession tially be a consequence of reduced RdDM pathway activity within PRJNA435649 and SRP133358 study. the gene body of the transcribed genes compared with TE. Furthermore, we also observed that the fraction of mCs in upstream Acknowledgements regions of transposons in fertile 303 R parent was slightly less than sterile 11 A parental line and the fertile hybrid. These observations We thank Council of Scientific and Industrial Research, Government of India suggest that one possible epigenetic mechanism of regulating trans- for the award of a senior research fellowship, ICAR-IARI, ICAR-NRCPB and poson inactivation in pigeonpea could be via methylation and this NRCPB in house project, New Delhi for supporting A.J. and providing other mechanism could be comparatively weaker in the fertile 303 R line funds and facilities. We also thank Dr N. Bhatia, The University of Sydney for than in the sterile parent and fertile hybrid. Over all, our results com- helpful feedback and ideas on the manuscript. ply with previously documented findings that CHH methylation lev- els are usually lower than CG and CHG, both on a genomic scale 65,66 Authors’ contributions and at specific loci. Further, our findings also reveal the existence of non-additive A.J and K.G conceived and designed the study. A.J. performed all the experi- methylation in F1 hybrid due TCM and TCdM events. TcM medi- mental work and data analysis. H.K. helped in data analysis and interpreting the results. A.J., K.G. wrote the manuscript. A.R. and N.K.S. helped editing ated increase in DNA methylation were more frequently observed in the manuscript. CG context and TcDM mediated decrease in DNA methylation were more frequently observed in CHH context, consistent with previous findings in Arabidopsis. Furthermore, we also dissected the mC Conflict of interest changes of differentially methylated parental epi-alleles in the hybrid, None declared. using SNPs for tracking the parental epi-alleles in the hybrid. Our analysis revealed that majority of the mC clusters having SNPs showed non- additive methylation in the hybrid when parents had Supplementary data different levels of methylation (HP or LP), consistent with our gen- eral mC cluster analysis. Supplementary data are available at DNARES online. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 372 The epigenomic interactions between parental inbreds 26. Chodavarapu, R. K., Feng, S., Ding, B., et al. 2012, Transcriptome and References methylome interactions in rice hybrids, Proc. Natl. Acad. Sci. U.S.A, 109, 1. Henderson, I. R. and Jacobsen, S. E. 2007, Epigenetic inheritance in 12040–5. plants, Nature, 447, 418–24. 27. Singh, N. K., Gupta, D. K., Jayaswal, P. K., et al. 2012, The first draft of the 2. Niederhuth, C. E. and Schmitz, R. J. 2014, Covering your bases: inheri- pigeonpea genome sequence, J. Plant Biochem. Biotechnol., 21, 98–112. tance of DNA methylation in plant genomes, Mol. Plant., 7, 472–80. 0 28. Saxena, K. B. and Singh, L. 1990, Low natural outcrossing in cleistoga- 3. Robertson, K. D. 2005, DNA methylation and human disease, Nat. Rev. 0 mous pigeonpea mutant, Int. Pigeonpea Newslett., 11, 9–10. ISSN Genet., 6, 597–610. 0255-786X. 4. Suzuki, M. M. and Bird, A. 2008, DNA methylation landscapes: provoca- 29. Chen, L. and Liu, Y. G. 2014, Male sterility and fertility restoration in tive insights from epigenomics, Nat Rev Genet., 9, 465–76. crops, Annu. Rev. Plant Biol., 65, 579–606. 5. Bird, A. 2007, Perceptions of epigenetics, Nature, 447, 396–8. 30. Tuteja, R., Saxena, R. K., Davila, J., et al. 2013, Cytoplasmic male 6. Law, J. A. and Jacobsen, S. E. 2010, Establishing, maintaining and modi- sterility-associated chimeric open reading frames identified by mitochon- fying DNA methylation patterns in plants and animals, Nat. Rev. Genet., drial genome sequencing of four Cajanus genotypes, DNA Res., 20, 11, 204–20. 485–95. 7. Vongs, A., Kakutani, T., Martienssen, R. A. and Richards, E. J. 1993, 31. Meshram, M. P. and Patil, A. N. 2018, Genetics of fertility restoration in Arabidopsis thaliana DNA methylation mutants, Science, 260, 1926–8., A2 cytoplasm based hybrids of pigeonpea. Int. J. Cur. Microbiol. Appl. 8. Kankel, M. W., Ramsey, D. E., Stokes, T. L., et al. 2003, Arabidopsis Sci., ISSN: 2319-7706 565-571. MET1 cytosine methyltransferase mutants, Genetics, 163, 1109–22. 32. Varshney, R. K., Chen, W., Li, Y., et al. 2011, Draft genome sequence of 9. Zemach, A., McDaniel, I. E., Silva, P. and Zilberman, D. 2010, pigeonpea (Cajanus cajan), an orphan legume crop of resource-poor farm- Genome-wide evolutionary analysis of eukaryotic DNA methylation, ers, Nat. Biotechnol., 30, 83–9. Science, 328, 916–9. 33. Xi, Y. and Li, W. 2009, BSMAP: whole genome bisulfite sequence 10. Cao, X. and Jacobsen, S. E. 2002, Role of the arabidopsis DRM methyl- MAPping program, BMC Bioinformatics, 10, 232. transferases in de novo DNA methylation and gene silencing, Curr. Biol., 0 0 34. Schultz, M. D., Schmitz, R. J. and Ecker, J. R. 2012, Leveling the playing 12, 1138–44. field for analyses of single-base resolution DNA methylomes, Trends 11. Mosher, R. A. and Melnyk, C. W. 2010, siRNAs and DNA methylation: Genet., 28, 583–5. seedy epigenetics, Trends Plant Sci., 15, 204–10. 35. Akalin, A., Kormaksson, M., Li, S., et al. 2012, methylKit: a comprehen- 12. Stroud, H., Do, T., Du, J., et al. 2014, Non-CG methylation patterns sive R package for the analysis of genome-wide DNA methylation profiles, shape the epigenetic landscape in Arabidopsis, Nat. Struct. Mol. Biol., 21, Genome Bio.l, 13,R87. 64–72. 36. Li, L. C. and Dahiya, R. 2002, MethPrimer: designing primers for methyl- 13. Reed, H. S. 1942, A short history of the plant sciences. https://www. ation PCRs, Bioinformatics, 18, 1427–31. scribd.com/document/190746815/Reed-Howard-S-1942-a-Short-History- 37. Kumaki, Y., Oda, M. and Okano, M. 2008, QUMA: quantification tool of-the-Plant-Sciences. for methylation analysis, Nucleic Acids Res., 36, W170–5. 14. Ng, S. Y., Lin, L., Soh, B. S. and Stanton, L. W. 2013, Long noncoding 38. Langmead, B. and Salzberg, S. L. 2012, Fast gapped-read alignment with RNAs in development and disease of the central nervous system, Trends Bowtie 2, Nat. Methods , 9, 357–9. Genet., 29, 461–8. 39. Li, B. and Dewey, C. N. 2011, RSEM: accurate transcript quantification 15. Greaves, I. K., Groszmann, M., Ying, H., Taylor, J. M., Peacock, W. J. from RNA-Seq data with or without a reference genome, BMC and Dennis, E. S. 2012, Trans chromosomal methylation in Arabidopsis Bioinformatics, 12, 323. hybrids, Proc. Natl. Acad. Sci. US. A, 109, 3570–5. 40. Robinson, M. D., McCarthy, D. J. and Smyth, G. K. 2010, edgeR: a bio- 16. Greaves, I. K., Groszmann, M., Wang, A., Peacock, W. J. and Dennis, E. conductor package for differential expression analysis of digital gene ex- S. 2014, Inheritance of trans chromosomal methylation patterns from pression data, Bioinformatics, 26, 139–40. Arabidopsis F1 hybrids, Proc. Natl. Acad. Sci. USA, 111, 2017–22. 41. Ye, J., Fang, L., Zheng, H., et al. 2006, WEGO: a web tool for plotting 17. Greaves, I. K., Gonzalez-Bayon, R., Wang, L., et al. 2015, Epigenetic GO annotations, Nucleic Acids Res , 34, W293–7. changes in hybrids, Plant Physiol., 168, 1197–205. 42. Liu, Y., Siegmund, K. D., Laird, P. W. and Berman, B. P. 2012, Bis-SNP: 18. Vaughn, M. W., Tanurdzic, M., Lippman, Z., et al. 2007, Epigenetic nat- combined DNA methylation and SNP calling for Bisulfite-seq data, ural variation in Arabidopsis thaliana, PLoS Biol., 5, e174. Genome Biol., 13,R61. 19. Schmitz, R. J., Schultz, M. D., Urich, M. A., et al. 2013, Patterns of popu- 43. Feng, S., Cokus, S. J., Zhang, X., et al. 2010, Conservation and diver- lation epigenomic diversity, Nature, 495, 193–8. gence of methylation patterning in plants and animals, Proc. Natl. Acad. 20. Shen, H., He, H., Li, J., et al. 2012, Genome-wide analysis of DNA meth- Sci. U.S.A., 107, 8689–94. ylation and gene expression changes in two Arabidopsis ecotypes and 44. Lisch, D. 2009, Epigenetic regulation of transposable elements in plants, their reciprocal hybrids, Plant Cell, 24, 875–92., Annu. Rev. Plant Biol., 60, 43–66. 21. Li, Y., Varala, K., Moose, S. P. and Hudson, M. E. 2012, The inheritance 45. Cokus, S. J., Feng, S., Zhang, X., et al. 2008, Shotgun bisulphite sequenc- pattern of 24 nt siRNA clusters in arabidopsis hybrids is influenced by ing of the Arabidopsis genome reveals DNA methylation patterning, proximity to transposable elements, PLoS One, 7, e47043. Nature, 452, 215. 22. Groszmann, M., Greaves, I. K., Albertyn, Z. I., Scofield, G. N., Peacock, 46. Qi, X., Li, M.-W., Xie, M., et al. 2014, Identification of a novel salt toler- W. J. and Dennis, E. S. 2011, Changes in 24-nt siRNA levels in ance gene in wild soybean by whole-genome sequencing, Nat. Commun., Arabidopsis hybrids suggest an epigenetic contribution to hybrid vigor, 5, 4340. Proc. Natl. Acad. Sci. U.S.A, 108, 2617–22. 47. The Arabidopsis Genome, I. 2000, Analysis of the genome sequence of the 23. He, G., Zhu, X., Elling, A. A., et al. 2010, Global epigenetic and tran- flowering plant Arabidopsis thaliana, Nature, 408, 796. scriptional trends among two rice subspecies and their reciprocal hybrids, 48. Zhang, J., Liu, Y., Xia, E. H., Yao, Q. Y., Liu, X. D. and Gao, L. Z. Plant Cell, 22, 17–33. 2015, Autotetraploid rice methylome analysis reveals methylation varia- 24. Barber, W. T., Zhang, W., Win, H., et al. 2012, Repeat associated small tion of transposable elements and their effects on gene expression, Proc. RNAs vary among parents and following hybridization in maize, Proc. Natl. Acad. Sci. U.S.A, 109, 10444–9. Natl. Acad. Sci. U.S.A, 112, E7022–9. 25. Parkin, I. A., Koh, C., Tang, H., et al. 2014, Transcriptome and methyl- 49. Li, X., Zhu, J., Hu, F., et al. 2012, Single-base resolution maps of culti- ome profiling reveals relics of genome dominance in the mesopolyploid vated and wild rice methylomes and regulatory roles of DNA methylation Brassica oleracea, Genome Biol., 15, R77. in plant gene expression, BMC Genomics, 13, 300. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 373 50. Lu, Y., Liu, Y., Wang, J., Cao, M. and Rong, T. 2010, Variation and pat- 62. Matzke, M., Kanno, T., Daxinger, L., Huettel, B. and Matzke, A. J. 2009, terns of DNA methylation in maize C-type CMS lines and their main- RNA-mediated chromatin-based silencing in plants, Curr. Opin. Cell Biol., 21, 367–76. tainers, J. Plant Biochem. Biotechnol., 19, 43–50. 63. Inagaki, S. and Kakutani, T. 2012, What triggers differential DNA meth- 51. Chen, X., Hu, J., Zhang, H. and Ding, Y. 2014, DNA methylation ylation of genes and TEs: contribution of body methylation? Cold Spring changes in photoperiod-thermo-sensitive male sterile rice PA64S under Harb. Symp. Quant. Biol., 77, 155–60. two different conditions, Gene, 537, 143–8. 64. Zakrzewski, F., Schmidt, M., Van Lijsebettens, M. and Schmidt, T. 2017, 52. Schmitz, R. J., Schultz, M. D., Lewsey, M. G., et al. 2011, DNA methylation of retrotransposons, DNA transposons and genes in Transgenerational epigenetic instability is a source of novel methylation sugar beet (Beta vulgaris L.), Plant J., 90, 1156–75. variants, Science, 334, 369–73. 65. West, P. T., Li, Q., Ji, L., et al. 2014, Genomic distribution of H3K9me2 53. Becker, C., Hagmann, J., Muller, J., et al. 2011, Spontaneous epigenetic and DNA methylation in a maize genome, PLoS One, 9, e105267. variation in the Arabidopsis thaliana methylome, Nature, 480, 245–9. 66. Gent, J. I., Ellis, N. A., Guo, L., et al. 2013, CHH islands: de novo DNA 54. Eichten, S. R., Briskine, R., Song, J., et al. 2013, Epigenetic and genetic in- methylation in near-gene chromatin regulation in maize, Genome Res., fluences on DNA methylation variation in maize populations, Plant Cell, 23, 628–37. 25, 2783–97. 67. Martienssen, R., Barkan, A., Taylor, W. C. and Freeling, M. 1990, 55. Liu, F., Cui, X., Horner, H. T., Weiner, H. and Schnable, P. S. 2001, Somatically heritable switches in the DNA modification of Mu transpos- Mitochondrial aldehyde dehydrogenase activity is required for male fertil- able elements monitored with a suppressible mutant in maize, Genes ity in maize, Plant Cell, 13, 1063–78. Dev., 4, 331–43. 56. Huang, W., Yu, C., Hu, J., et al. 2015, Pentatricopeptide-repeat family 68. Shi, J., Dong, A. and Shen, W. H. 2014, Epigenetic regulation of rice protein RF6 functions with hexokinase 6 to rescue rice cytoplasmic male flowering and reproduction, Front. Plant Sci., 5, 803. sterility, Proc. Natl. Acad. Sci. U.S.A, 112, 14984–9. 69. Higo, H., Tahir, M., Takashima, K., et al. 2012, DDM1 (decrease in 57. Song, Q., Guan, X. and Chen, Z. J. 2015, Dynamic roles for small RNAs DNA methylation) genes in rice (Oryza sativa), Mol. Genet. Genomics, and DNA methylation during ovule and fiber development in allotetra- 287, 785–92. ploid cotton, PLoS Genet., 11, e1005724. 70. Choi, S. C., Lee, S., Kim, S. R., et al. 2014, Trithorax group protein 58. Xu, W., Yang, T., Dong, X., Li, D. Z. and Liu, A. 2016, Genomic DNA Oryza sativa Trithorax1 controls flowering time in rice via interaction methylation analyses reveal the distinct profiles in castor bean seeds with with early heading date3, Plant Physiol., 164, 1326–37., persistent endosperms, Plant Physiol., 171, 1242–58. 71. Zhang, X., Yazaki, J., Sundaresan, A., et al. 2006, Genome-wide 59. Lu, X., Wang, W., Ren, W., et al. 2015, Genome-wide epigenetic high-resolution mapping and functional analysis of DNA methylation in regulation of gene transcription in maize seeds, PLoS One, 10, arabidopsis, Cell, 126, 1189–201., e0139582. 72. Zilberman, D., Gehring, M., Tran, R. K., Ballinger, T. and Henikoff, S. 60. Garg, R., Narayana Chevala, V., Shankar, R. and Jain, M. 2015, 2007, Genome-wide analysis of Arabidopsis thaliana DNA methylation Divergent DNA methylation patterns associated with gene expression in uncovers an interdependence between methylation and transcription, Nat rice cultivars with contrasting drought and salinity stress response, Sci. Genet., 39, 61–9. Rep., 5, 14922. 73. Schmitz, R. J., He, Y., Valdes-Lopez, O., et al. 2013, Epigenome-wide in- 61. Song, Q. X., Lu, X., Li, Q. T., et al. 2013, Genome-wide analysis of DNA heritance of cytosine methylation variants in a recombinant inbred popu- methylation in soybean, Mol. Plant., 6, 1961–74. lation, Genome Res., 23, 1663–74. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png DNA Research Oxford University Press

Unravelling the epigenomic interactions between parental inbreds resulting in an altered hybrid methylome in pigeonpea

DNA Research , Volume 25 (4) – Aug 1, 2018
Download PDF
13 pages

Loading next page...
 
/lp/ou_press/unravelling-the-epigenomic-interactions-between-parental-inbreds-Q6wRQBQxtg

References (86)

Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of Kazusa DNA Research Institute.
ISSN
1340-2838
eISSN
1756-1663
DOI
10.1093/dnares/dsy008
pmid
29566130
Publisher site
See Article on Publisher Site

Abstract

DNA methylation is an important heritable landmark conferring epigenetic changes in hybrids and has fascinated biologists and plant-breeders over the years. Although epigenetic changes have been documented in rice and maize hybrids, such investigations have not been reported in pigeonpea. Here, we report genome-wide methylation profiles of pigeonpea sterile and fertile inbred lines and their fertile F1 hybrid at single base resolution. We found that pigeonpea genome is relatively enriched in CG methylation. Identification of differentially methylated re- gions (DMRs) in the sterile and fertile parent revealed remarkable differences between their methylation patterns. Investigation of methylation status of parental DMRs in hybrid revealed non-additive methylation patterns resulting from trans-chromosomal methylation and trans- chromosomal demethylation events. Furthermore, we discovered several DMRs negatively associated with gene expression in the hybrid and fertile parent. Interestingly, many of those DMRs belonged to transposable elements and genes encoding pentatricopeptide repeats asso- ciated proteins, which may mediate a role in modulating the genes impacting pollen fertility. Overall, our findings provide an understanding of two parental epigenomes interacting to give rise to an altered methylome in pigeonpea hybrids, from genome-wide point of view. Key words: DNA methylation, TcM, TcDM, pigeonpea, hybrid 1. Introduction C-5 position of the cytosine ring by DNA methyltransferases and is Epigenetic regulation has remained of fundamental interest in under- transmitted across generations. DNA methylation acts to regulate standing gene regulation in eukaryotes. Epigenetic regulation in- transposons and repetitive genome sequences and often results in the 4,5 volves differential methylation at cytosine residues and post- suppression of gene activity. A major difference between DNA translational modifications of histones. One of the most comprehen- methylation in animals and plants is that in plant genome, DNA sively investigated mechanism of chromatin modification in plants, methylation occurs at CG, CHH and CHG sequence contexts (where DNA methylation involves covalent transfer of a methyl group to H could be C, T or A). It has been shown in Arabidopsis that V C The Author(s) 2018. Published by Oxford University Press on behalf of Kazusa DNA Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] 361 Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 362 The epigenomic interactions between parental inbreds different enzymes regulate methylation at different sequence con- variation in pigeonpea hybrids. This particular property has become texts. Methyltransferase 1 (MET1) causes and maintains methylation popular among breeders to produce varieties with improved fertility 6–8 at CG sites and chromomethylase 3 (CMT3) produces and main- and been utilized through cytoplasmic male sterility (CMs) system. 6,9,10 tains methylation at CHG sites and this maintenance is further CMS is a condition in which the ability of a plant to produce fertile required to preserve methylation patterns over generations. On the pollen is severely impaired. In this condition, rearrangements in mito- other hand, domains rearranged methyltransferase 2 (DRM2) and chondrial genome result in chimeric mitochondrial open reading 24 nt siRNAs, which guide methyltransfrase to specific DNA seg- frames (ORFs), generating non-functionality in anthers. However, the 6,10,11 ments, regulate de novo methylation CHH sites, which is main- expression of such chimeric mitochondrial ORFs can be suppressed by tained by CMT2. fertility restorer gene in F1 hybrid, thus restoring fertility. This kind of De novo methylation of all three different sites involves 24 nt hybrid breeding technology involves three parents— a male sterile par- smRNAs (small RNAs), DRM2 and the Dnmt3 homologue involved ent, its maintainer and a fertility restorer line (commonly called—A, B in RNA-directed DNA methylation (RdDM). In RdDM, RNA is and R lines, respectively; reviewed in detail in the literature). Few transcribed by the RNA polymerase POLIV, which is then made CMS associated mitochondrial orfs have been identified in double stranded by RNA-dependent RNA polymerases (RDRs), pigeonpea. cleaved to 24 nt by Dicer-like enzymes and loaded onto argonaute In this study, we have investigated whether there are any epige- proteins. This effector complex then interacts with other proteins netic changes in F1 progeny of cross between 11 A (male sterile) and and Polymerase V transcripts to target homologous regions of the ge- 303 R (fertile) inbred lines and whether those changes could poten- nome for methylation by the recruitment of DRM2. The methylation tially mediate an impact on the gene expression patterns including state of CG is maintained during cell replication by MET1, which those involved in fertility restoration in the hybrid. We conducted recognizes hemi methylated CG sites and methylate the newly whole genome bisulfite sequencing of 11 A, 303 R parental genotypes formed DNA strand. CMT3 maintains CHG methylation, while and their hybrid to generate whole genome cytosine methylation pro- CHH methylation is recognized by a histone methyltransferase files and found that overall, hybrid methylome was remarkably al- kryptonite (KYP), which can then methylate the H3K9 residues. tered compared with the parental methylome due to the events of These two enzymes work in a feedback loop enabling the mainte- transchrosomal methylation and demethylation in hybrid. We next nance of both epigenetic marks, with the absence of either CMT3 or carried out transcriptome profiling of the same plants and compared KYP leading to a reduction in CHG methylation. DNA methylation methylation status of the differentially expressed genes between pa- plays an important role during different aspects of plant development rental inbreds and the F1 hybrids. We found differences in methyla- and growth. Of these heterosis or hybrid vigour is of vital impor- tion patterns in the three genotypes, suggesting the involvement of tance in agriculture. Heterosis is a phenomenon where F1 progeny of DNA methylation in regulation of gene expression and possibly a genetic cross exhibits superior characteristics comparedwith its in- influencing sterility and fertility. bred parental lines and has been documented nearly 300 years ago. In recent years, several studies have focused on the role of epigenetic 2. Materials and methods regulation in mediating hybrid vigour in model plant 1,14–17 Arabidopsis. Previous studies in Arabidopsis have suggested 2.1. Plant material and growth conditions that interaction between similar genomes and different epigenomes Pigeonpea lines 11 A and 303 R were crossed by hand pollination of the two different parental ecotypes could result in new transcrip- using 303 R as pollen donor to generate hybrid. Seeds were collected tional profiles in hybrids that might contribute to hybrid vigour in and F1 seed resulting from the cross along with seeds from parental 18,19 the progeny. Also, it has been proposed that epigenetic regula- lines were sown again together for further experimental work. Seeds tion via changes in methylation patterns might play a role in mediat- were directly sown in the soil in green house under natural condi- ing differences in gene expression levels that may confer increased tions. All the plants were grown at the same time under exactly simi- biomass in heterotic hybrid offspring of Arabidopsis ecotypes C24 lar growth conditions to reduce the impact of any environmental and Landsberg erecta. It has also been documented that both DNA variations. Fertility status of all three genotypes was verified by pol- methylation patterns and sRNA (small RNA) levels are altered in in- len analysis. Pollens were fixed on a glass slide using 2% acetocar- 15,20–22 23 24 traspecific hybrids of Arabidopsis, rice and maize. mine. Pollens stained with acetocarmine were considered to be Recent advances in DNA methylation profiling techniques have en- fertile. After pollen verification, young leaves of same age from the abled in understanding the correlation between methylation patterns three lines were harvested, frozen in liquid nitrogen and used to ex- and hybrid vigour in several crop species as well such as rice and tract DNA and RNA for further work. 25,26 Brassica but knowledge about interaction of methylation patterns and gene expression in legume plants such as pigeonpea has still re- 2.2. Whole genome bisulfite sequencing library mained elusive. Pigeonpea (Cajanus cajan (L) Millsp) is a important construction and data analysis grain crop with high protein (20–22%) grain legume, mainly culti- High molecular weight genomic DNA was isolated from the leaves vated in the tropical and subtropical parts of the world. Pigeonpea using cetyl trimethylammonium bromide method. DNA was frag- contains 11 pairs of chromosomes (2n¼ 22) and has a genome size of 858 Mbp. Globally pigeonpea is cultivated on 5.4 Mha, with an an- mented into 100–300 bp using Covaris followed by end repair by nual production of 4.4 MT and a yield of 8299 hg/ha. Pigeonpea is adding ‘A’ and True seq methylated adapters were ligated. The li- grown in about 50 countries with India occupying 66% of its area gated DNA fragments were subjected to bisulfite treatment using EZ (FAO, 2016, http://www.fao.org/faostat/en/#data/QC). India is also DNA Methylation-Gold Kit, according to manufacturer’s instruc- the largest producer and consumer of pigeonpea with an annual pro- tions. The BS libraries were amplified, quality tested using Bio- duction of 2.51 MT. analyzer (Agilent 2100) and subjected to Ilumina Hi-seq2000 for 90 Pigeonpea is characterized by significant genetic variability. Genetic cycles in paired-end mode. The adaptors and poor quality reads were studies have also revealed the occurrence of non-additive genetic removed from the raw data using TrimGalore (https://github.com/ Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 363 FelixKrueger/TrimGalore). High quality reads were mapped to the Table S2) using Bowtie alignment program, then relative abundance 27,32 33 39 reference genome using BSMAP at default parameters. of transcripts was calculated using RSEM. Transcripts with zero Methylation information was extracted using Python script. Reads Fragments Per Kilobase of transcript per Million mapped reads were mapped on chloroplast genome to calculate methylation effi- (FPKM) value were excluded from analysis. Raw reads count data ciency and error rate. A bisulfite conversion rate of >99% was ob- were imported into R environment and data table containing fragment tained for all three samples. In-order to determine true methylated counts for each transcript were combined and matrix was formed to cytosines, P-values were calculated using binomial distribution identify differentially expressed gene using edgeR. P¼ binomial (m, x, error rate), where m¼ number of methylated reads, n¼ number of non-methylated reads and x¼ m þ n (coverage 2.6. GO enrichment analysis on a single C). Cytosines with P-value< 0.0025 were considered to To identify and categorize cluster of differentially expressed genes be truly methylated. Fraction and mean methylation level was calcu- (DEGs) and genes enriched for methylation in all the three categories lated as described previously. (11 A vs. 303 R, 11 A vs. hybrid and 303 R vs. hybrid), GO enrich- ment analysis was performed using python based program 2.3. Identification of differentially methylated regions GOATOOLS (https://github.com/tanghaibao/goatools) and Wego. To identify differentially methylated regions (DMRs) at each se- Genes with Bonferroni-corrected P-values< 0.05 using GOATOOLS quence context (CG, CHG and CHH), genomes were tiled using a were considered to be enriched. sliding window approach with 100 bp window sliding at 100 bp in- terval for CG and CHG and 150 bp for CHH. Methylation informa- 2.7. Identification of transposable elements tion for each window was calculated using Methyl kit. A denovo repeat library was generated using repeatmodeller (http:// P-value for each window was calculated using Fisher’s exact test. www.repeatmasker.org/RepeatModeler/; version 1.0.11).This cus- To determine q-values, P-values for each window were corrected us- tomized library was used to run RepeatMasker (http://www.repeat ing sliding linear model . Following criteria were used for identifica- masker.org/; version- 4.0.7) to identify and annotate transposable el- tion of DMRs. ements (TEs) in pigeonpea at default parameters. TE shorter than 100 bp were excluded from the analysis and the distance between 1. Windows contained at-least three methylated cytosines for CG TE, less than 50 bp were conconated (Supplementary Table S3). and CHG; and six for CHH DMRs. 2. Bases having more than 10X coverage were considered. 3. Bases having more than 99.9th percentile coverage were discarded 2.8. Allele specific methylation analysis using single to minimize errors arising due to PCR bias. nucleotide polymorphism 4. Windows with methylation difference greater than 20% for CG, Single nucleotide polymorphisms (SNPs) from parental lines and hy- CHG and 25% for CHH were selected. 42 brid were identified using BIS–SNP, a variant caller from BS reads. 5. Fisher’s exact test corrected q-values < 0.01. For this, the mapped methylated reads were used as input files along with pigeonpea reference genome. SNPs detection was performed 2.4. Validation of methylation levels using bisulfite using the following parameters: ‘-T BisulfiteGenotyper -stand_call_- sequencing PCR conf 20 -stand_emit_conf 0 -mmq 30 -mbq 17 -minConv 20’. Spurious SNP with quality score less than 20, more than 120x read Two different regions were chosen one for validating methylation coverage, strand bias more than—0.02, quality score by depth less levels, one in the intergenic region and the other one in the exonic re- than 1.0, mapping quality zero reads, fraction more than 0.1 and 2 gion. For bisulfite PCR, genomic DNA was extracted from 11 A SNPs within the 20 bp bin were further filtered out using VCF post and F1 hybrid and subjected to bisulfite conversion using EZ TM process command. DNA Methylation-Gold Kit. Primers were designed using METHPRIME (Supplementary Table S1). Bisulfite converted DNA was amplified and the fragments were cloned into pGEMT vector. 3. Results Ten colonies were picked randomly for sequencing. The results were analyzed using QUMA. 3.1. Genome wide methylation analysis in pigeonpea To understand the possible roles of epigenetic regulation and the dy- 2.5. RNA sequencing library construction and analysis namics of DNA methylation during hybrid formation in pigeonpea, we crossed 303 R (Fertile, pollen donor) and 11 A(sterile) lines. Since Total RNA was extracted from the same plants used for BS Seq. Equal their F1 progeny restored fertility (Fig. 1A), we hypothesised that it amount of total RNA from three biological replicates were pooled. About 1 mg of the total pooled RNA was used to purify poly-A con- could possibly partially attributed to the new epigenome resulting taining mRNA molecules. The purified mRNA was fragmented into from the interaction between the two parental epigenomes. Thus, in- smaller pieces and reverse transcribed into cDNA. The cDNA frag- order to gain deep insights into the epigenome of the parental lines ments were subjected to end repair by adding single ‘A’ base. True seq and their F1 hybrid, we deciphered their respective whole genome adaptors were ligated and cDNA libraries were amplified for 15 cy- methylation pattern at single base resolution. For this high molecular weight genomic DNA was isolated from the leaves of 10 days after cles. Quality of the library was assessed using Bioanalyzer and the vali- flowering (DAF) old plants from three different genotypes, namely dated library was sequenced by Illumina Hi-Seq 1000 according to 11 A, 303 R and F1 hybrid progeny of 11 A and 303 R. The DNA manufacturers’ instructions. The quality of the raw reads was deter- samples were subjected to bisulfite treatment to generate genomic li- mined using FastQc (https://www.illumina.com/products/by-type/infor matics-products/basespace-sequence-hub/apps/fastqc.html). Adaptors braries for high throughput bisulfite sequencing. The reads obtained and poor quality reads were cleaned using Trim Galore. High quality were processed for removal of adapter and low quality sequences. paired end reads were aligned to reference genome (Supplementary About 80–100 million cleaned reads (in each of the three genotypes) Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 364 The epigenomic interactions between parental inbreds Figure 1. DNA methylation profiles of pigeonpea inbred lines (11 A, 303 R) and their hybrid. (A) Acetocarmine stained pollens of 11A (sterile parent), 303R (fertile parent) and fertile hybrid. Note the density of pollen is less in 11A compared to 303R and hybrid. Also note viable pollens in 303R and hybrid.(B) Genome-wide relative proportions of CG (blue), CHG (red) and CHH (green) methyl cytosine in 11A, 303R and F1 hybrid. (C and D) Fractional methylation levels in the three lines of pigeonpea in CG, CHG and CHH contexts (B) and comparisons with published reports in Soybean, Rice, Populus, Arabidopsis and Brassica (D). (E) A cir- cos plot representing chromosome wise distribution of mCG, mCHG, mCHH, gene, transposons density and FPKM (outer to inner) in F1 hybrid. Yellow to red color relates to lower to higher methylation density in a 100 Kb window (Colored figure is available in the online version of the article). Scale bars 50mm (A). 27,32 were mapped on to the reference genome. We found that more cytosine methylation at single base resolution. Further, to calculate than 80% reads got mapped on to the reference genome. Thus even- the frequency of cytosine conversion, reads were mapped to un- tually, 46 408 326, 54 030 096 and 43 763 469 reads were obtained methylated chloroplast genome. We obtained a conversion rate for 11 A line, 303 R line and their F1 hybrid progeny, respectively, of more than 99% in all the three genotypes, indicating high effi- that mapped uniquely to the reference genome and providing a cov- ciency of cytosine conversion in our experiments (Supplementary erage more than 20X. These were subjected to further analysis of Table S4). Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 365 We then deciphered the fractional proportion of DNA methyla- 3.2. Patterns of DNA methylation in genic and tion and found that CG methylation was highest in 11 A (41.7%) transposable elements while CHH methylation was highest in F1 hybrid (28.82%). To precisely understand the methylation patterns of the protein cod- Methylation in CHG context were similar in all the three genotypes ing genes, genes with inserted TE elements were excluded from (34%) (Fig. 1B). Highest mean methylation was observed in CG analysis. We calculated the methylation levels in the genic region and context (74%-79%), closely followed by CHG (56%–63%) while 2 kb upstream and downstream regions of the gene bodies in the two mean methylation in CHH context (10–15%) was lowest in both the parental lines and their hybrid and found significant methylation parents and the hybrid (Fig. 1C). Comparative analysis of fractional differences between the three lines (Wilcoxon rank sum test, P-val- methylation levels in pigeonpea, soybean, brassica, oryza and ue< 2.2e-16, n ¼ 26,774). These differences reflect differential meth- populus revealed that pigeonpea genome was hypermethylated com- ylation patterns in the genic regions of parental lines and their pared with the other crops. Interestingly, CHH methylation levels in hybrid. Furthermore, upon a detailed investigation of methylation pigeonpea genome were significantly higher than those observed in pattern, we observed the occurrence of a single CG methylation peak 25,43 soybean, brassica, oryza and populus (Fig. 1D). It has been doc- within the gene body and two troughs near transcriptional start and umented that the DNA repeats in the genomes are highly methylated termination sites in pigeonpea, invariably to earlier reports in other 9,44,45 43,48,49 48 in symmetric contexts. Pigeonpea geneome is highly repetitive plants. However, contrary to a recent report in rice, we 46 47 (47%) compared with soyabean and Arabidopsis. Thus, such found over all low CG methylation levels within the coding regions high levels of DNA methylation in pigeonpea may arise from the occ- compared with their flanking regions. A similar trend was observed curance highly repetive genome. for CHG methylation, however, the average CHG methylation both Visualization of chromosome wide distribution of CG, CHG and within the coding regions and their flanking regions were compara- CHH methylation revealed higher CG and CHG methylation levels tively lower than CG methylation. CHH methylation however compared with CHH while relatively even distribution pattern of showed a different trend to both CG and CHG methylation in terms CHH methylation (Fig. 1E, Supplementary Fig. S1A and B). High of a plateau instead of a single peak within the gene body. Overall CG and CHG methylation were correlated with transposon rich re- CHH methylation both within the genic regions and their flanking gions across the genome (Fig. 1E, Supplementary Fig. S1A and B). regions were significantly attenuated than CG and CHG Figure 2. Distribution of methyl cytosines in genes and transposable elements in pigeonpea. (A) Line plots showing methylation level in 2Kb flanking (upstream and downstream) and within the gene body in 11A, 303R and F1 hybrid in CG, CHG and CHH contexts. (B–D) Line plots showing methylation level in 2Kb flank- ing (upstream and downstream) and within the gene body of transposons- long terminal repeats (B); LINE (C) and DNA Elements (D) in 11A, 303R and F1 hybrid in CG, CHG and CHH context. Note higher mCs throughout the transposons compared to genes (Colored figure is available in the online version of the article). Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 366 The epigenomic interactions between parental inbreds methylation. However, the hybrid exhibited mid parental range for different criteria were applied for CHH DMR identification. To as- CG and CHG but a higher CHH methylation compared with paren- sign a DMR to CHH category, each window had at-least eight cyto- tal lines, suggesting altered CHH methylation patterns in the genic sines with minimum of 10X coverage and a minimum of 25% regions in pigeonpea hybrids (Fig. 2A). difference in CHH methylation. Eventually, we identified much more Eukaryote genomes are abundant in TEs, also called mobile ele- CG DMRs (11,300) than CHG (6988) and CHH (8061) DMRs be- ments. TEs are structurally divergent, epigenetically controlled and are tween 303 R and 11 A parental lines (Supplementary Table S5). We mostly involved in genome evolution and gene activity regulation. also found that the differential methylation between the two parents Thus, we sought to gain a comprehensive insight of the influence of was more enriched in upstream flanking sequences of the genes than DNA methylation on the transposable elements within the pigeonpea the down stream regions (Fig. 3A). We next assessed relative methyl- genome. Since no complete dataset on pigeonpea TEs is available, we ation percentages of the obtained DMRs between the 11 A and used repeat masker to identify and classify them. TEs shorter than 303 R parents (11 A vs. 303 R or Sterile vs. fertile) and found enrich- 100 bp were excluded from analysis. In total, we found 105 733 TEs, ment for hypermethylated CHH DMRs and hypomethylated CG which we classified into retro transposons and DNA elements. We DMRs (Fig. 3B). This suggests that in 11 A line, genes involved in next calculated their CG, CHG and CHH methylation levels in all the regulation of CG methylation are less active than the RdDM path- three genotypes and found that irrespective of their categories, methyl- way genes and CMT2 that are involved in the regulation CHH meth- ation was highly enriched in all the contexts in TEs compared with the ylation (Supplementary Fig. S2). gene body. However, interestingly, we found relatively lower CG and Interaction of two dissimilar parental epigenomes can result in an CHG methylation levels in the TEs flanking regions than those of the altered methylome in their hybrid progeny. Since fertility is restored genic regions. By comparing the methylation levels in long terminal in the F1 hybrid and earlier studies in Arabidopsis showed higher 15,20 repeats within the three genotypes, we found relatively higher methyl- methylation levels in hybrid compared with parental lines, we ation levels in F1 hybrid in all the three contexts, compared with the therefore investigated whether there is differential methylation be- parental lines except for CG methylation in flanking regions where tween individual parental lines and the fertile hybrid (Supplementary 11 A exhibited higher methylation levels (Fig. 2B). While, in the case Table S5). Consistent with previous studies, we found a relative of LINE elements, 11 A displayed hyper methylation compared with increase in the percentage of methylated cytosines in CG, CHG and the F1 hybrid and 303 R in CG and CHG contexts (Fig. 2C). DNA CHH DMRs in the hybrid compared with parental lines (Fig. 3B) in- elements were relatively hypermethylated than the class 1 transposons. dicative of occurrence of non-additive methylation in pigeonpea hy- F1 hybrid displayed uniformly higher methylation irrespective of the brids. We randomly chose hyper and hypo methylated DMRs mC contexts (Fig. 2D). Overall we observed comparatively lower between sterile parent and F1 hybrid in intergenic region and coding methylation levels in 303 R, compared with 11 A and F1 hybrid regions. Analysis of their methylation levels using bisulfite sequenc- (Fig. 2B–D), suggesting differences in the regulatory mechanisms of ing PCR revealed relatively higher methylation levels in intergenic re- TE silencing via methylation in pigeonpea. gion associated DMR than in the coding region associated DMR (Supplementary Fig. S3) consistent with our in-silico data. In-order to assess the functional categories of the genes associated 3.3. Fertile and sterile pigeonpea inbred lines display with the DMRs between parents and hybrids and to investigate locus specific differences in cytosine methylation whether DMRs were associated with genes involved in DNA methyla- Earlier work in crop plants such as rice and maize have documented tion, epigenetic processes or fertility regulation, gene ontology enrich- differences in methylation levels between a male sterile plant and a ment analysis was performed. We found genes ontology terms for fertile plant and revealed hyper-methylation in sterile lines compared regulation of flowering time, DNA methylation, regulation of gene si- 50,51 with fertile lines. This suggests differential cytosine methylation lencing by miRNA, developmental processes, Histone H4 K4 methyla- could possibly influence fertility and sterility in plants by regulating tion, DNA methyltransferase, transcription factors, proton gene expression. Therefore to investigate whether pigeonpea 11 A transporting ATP synthase activity, regulation of post transcriptional and 303 R inbreds exhibit differences in their methylation levels, we gene silencing were enriched (Fig. 3C). Interestingly, we also observed compared the methylome of the two lines. Several studies have docu- a strong enrichment of DMRs in gene is associated with aldehyde de- mented that viewing differential methylation in windows of multiple hydrogenase activity. A previous study have reported that RF2 cytosines seems to reflect more functional consequences rather than encoded aldehyde dehydrogenase activity is required to restore male 52,53 viewing at single sites. Also it is known that CG, CHG and fertility in maize and also plays a role in anther development. Our CHH methylation levels are maintained through different regulatory results thus suggest a possible correlation between differential methyl- mechanisms. Thus, to identify DMRs between 11 A and 303 R lines, ation and fertility or sterility in 11 A, 303 R pigeonpea inbred and we decided to investigate and categorize differential methylation for their hybrid. Also, genes involved in cellular processes, biosynthetic all the three contexts, CG, CHG and CHH, separately. We averaged processes and ATP binding were significantly enriched (Bonferroni the DNA methylation for all the three contexts independently across corrected P-value< 0.05; Fig. 3C, Supplementary data S1). non-overlapping 100 bp regions (windows), as done previously. Overall, our results indicate that differential methylation could be In-order to categorize a DMR window as CG, CGH or CHH, differ- involved in mediating and regulating several cellular and develop- ent conditions were applied. To assign a DMR to CG or CHG cate- mental processes in pigeonpea. gory, windows having at-least three cytosines with a minimum of 10X coverage were used. A minimum difference of 20% in CG and 3.4. Trans-chromosomal methylation events alter in CHG methylation each was chosen. Previously, it has been re- pigeonpea hybrid methylome ported that CHH methylation patterns appeared to be more variable than CG and CHG. Since CHH methylation sites are much abun- Our results so far show altered methylome in pigeonpea hybrid com- dant in genomes and still, the CHH methylation levels were observed pared with parents. The changes in hybrid methylome occur at to be much lower than the CG and CHG methylation levels (Fig. 1), DMRs between the parents and are mediated by either increase or Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 367 Figure 3. Comparative analyses of DMRs between the inbred lines and their hybrid. (A) Pie charts showing the distribution of DMRs in the different genomic re- gions in 11A vs. 303R, 11A vs. hybrid and 303R vs. hybrid. (B) A bar plot showing number of hypo and hypermethylated CG, CHG and CHH DMRs in 11A vs. 303R, 11A vs. hybrid and 303R vs. hybrid. (C) A heat map showing enrichment of DMR associated genes in 11A vs. 303R, 11A vs. hybrid and 303R vs. hybrid. Color code represents low to high P-values (0.05 to 8.09e-10) (Colored figure is available in the online version of the article). decrease in DNA methylation termed as non-additive methyla- of 303 R fertile parent and hybrid were much higher than the 11 A 15,20 tion. This increase or decrease occurs via altering methylation sterile parent. In contrast, CHH methylation levels of 303 R fertile parent and hybrid were lower than that of the 11 A sterile parent status of one parental allele to resemble that of the other parental al- lele, events known as trans-chromosomal methylation (TCM) and (Fig. 4A). Additionally, to decipher the number of direct TcM and trans-chromosomal demethylation (TCdM). We investigated non- TcDM non additive mC clusters, we used an alternative approach of SNPs as described previously. For this, we tracked 75,340 mC clus- additive methylation in F1 hybrid to understand its methylome dy- ters of parental origin and followed the methylation pattern of pa- namics. To understand this mechanism, DMRs between the two pa- rental lines were identified and mapped onto the F1 hybrid genome. rental alleles with differences in methylation level (3,222 out of 75,340) in hybrid using SNPs. We observed 27% mC clusters were Average methylation of the mapped loci was calculated in F1 hybrid. stably inherited in hybrids indicative of additive methylation. Among Mean parental value (MPV) of the DMRs between parents were the clusters showing non additive methylation, we observed, TCM obtained and compared with observed methylation in hybrid. Loci events in 18.5% cases where methylation of epi-alleles from low having methylation difference of atleast 20% were considered as parent increased in the hybrid, whereas methylation of high parent non-additive methylation sites. Subsequently, we found non-additive epi-alleles remained unchanged in the hybrid, resulting in an over all methylation in 28% DMRs suggesting methylation interaction at increase in methylation in the hybrid. While in 24.1% clusters, these regions in F1 hybrid. We further analyzed the nature of this TCdM was observed where methylation of the epi-allele from low non-additive methylation and observed patterns of Trans- parent remained unchanged but methylation of the epi-allele from chromosomal methylation at 3,911 DMRs and Trans-chromosomal high parent decreased in hybrid, resulting in an overall decrease in demethylation in 5,377 DMRs (Fig. 4A and B). By investigating the methylation in hybrid (Fig. 4C). genomic features of the DMRs showing non-additive methylation, we found a higher proportion resided in TE (25%) as compared with the genic regions (15%). These results are in accord with previous 3.5. Transcriptome profiling in pigeonpea findings revealing the roles of TE in genome stability and their regu- In order to assess the transcript abundance and differential gene ex- lation by methylation. Furthermore we found that in all the three pression between the three genotypes, their whole transcriptome pro- contexts, the methylome of F1 hybrid resembled more of the 303 R filing was performed by RNA sequencing. We first compared the fertile parent than the 11 A sterile. CG and CHG methylation levels transcription between different genotypes (11 A vs. 303 R, 11 A vs. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 368 The epigenomic interactions between parental inbreds Figure 4. (A) Heat maps showing the methylation status of DMRs of the parents and methylation pattern at those loci in the F1 hybrid. Note in all the mCs con- text, methylation patterns of F1 hybrid are similar to the fertile parent (303 R) and different to the sterile parent. (B) Bar plot showing non-additive methylation in F1 hybrid. Note increased average methylation in hybrid than the MPV in the upper graph indicating the occurrence of TCM. Similarly, note a decreased aver- age methylation in hybrid than the MPV in the lower upper graph indicating TCdM. (C) Bar plots showing changes in parental allelic methylations in the hybrid. Majority of the non-additive methylation occurred through TCdM (Blue). Red bar represent TcM. Methylation levels of most parental alleles remained unchanged in hybrid (gray) (Colored figure is available in the online version of the article). HP: high parent; LP: low parent. F1 hybrid and 303 R vs. F1 hybrid) to compute genes that were dif- that gene ontology terms included reproductive processes, develop- ferentially expressed between these lines. The transcripts showing a mental processes, protein–DNA complexes and nucleic acid binding minimum of 2-fold change and P-value<0.05 were considered dif- proteins that were significantly enriched (Fisher’s exact test P-val- ferentially expressed. We found 1,621 (460 up, 1,161 down), 1,790 ue<0.05) in all the three comparative scenarios 11 A vs. 303 R, 11 A (1,052 up, 738 down) and 1,380 (963 up, 417 down) genes that vs. hybrid and 303 R vs. hybrid. showed differential expression in 11 A vs. 303 R, 11 A vs. hybrid and 303 R vs. hybrid, respectively. We found that in all the three cases of 3.6. Correlation between DMRs and gene expression in comparison, the differentially expressed genes were well represented pigeonpea inbred and hybrid line in gene categories like pentatricopeptide repeats (PPRs)- associated proteins involved in fertility restoration, methylases, demethylases We next asked if the differential methylation in the three different ge- notypes could potentially influence the gene expression in these and histone modifying genes (that play a role in epigenetic regula- tion) and flowering-associated genes. We performed gene ontology plants. For this, first DEGs with DMRs localized in 2 kb flanking (up- stream and downstream) and gene body regions were identified as enrichment analysis to identify functional categories of significantly enriched genes (Supplementary Figs S4–S6). Interestingly, we found DMR-associated DEGs. Then the correlation between whole genome Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 369 differential methylation and gene expression was examined investigated the correlation between methylation status and expres- (Fig. 5A–C). We next determined the genes that were enriched for sion in the DEGs associated with DMRs and found that 127 (16 DMRs between the three genotypes and found that DMRs were local- hypomethylated and up-regulated; 111 hyper-methylated and down- ized in 3,019, 1,996 and 2,233 genes in 11 A vs. hybrid, 11 A vs. regulated), 75 (29 hypo-methylated and up-regulated; 46 hyperme- 303 R and 303 R vs. hybrid, respectively (Fig. 5A–C). Among these, thylated and down-regulated) and 229 (197 hypo-methylated and up- 290, 586 and 578 genes were found to be both differentially ex- regulated; 22 hyper-methylated and down-regulated) genes showed a pressed and associated with DMRs in 11 A vs. F, 11 A vs. 303 R and negative correlation between their expression and methylation in 11 303 R vs. F, respectively (Fig. 5A–C). Furthermore, we also A vs. R, 11 A vs. Hybrid and 303 R vs. Hybrid, respectively (Fig. 5D Figure 5. Characterization of DMRs associated DEGs in pigeonpea. (A–C) Venn diagrams showing unique and common DEGs (differentially expressed genes) to the genes associated with DMRs in 11A vs. hybrid (A), 11A vs. 303R (B) and 303R vs. hybrid (C). (D) Bar plots showing number of DEGs that were differentially methylated in 11A vs. 303R, 11A vs. hybrid and 303R vs. hybrid. (E–G) Heat maps displaying comparative methylation levels of selective genes that showed neg- ative correlation between their expression and methylation in 11A vs. 303R (E), 11A vs. hybrid (F) and 303R vs. hybrid (G). Color code indicates methylation level (0–1) (Colored figure is available in the online version of the article). Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 370 The epigenomic interactions between parental inbreds and Supplementary data S2). In order to gain insights into the biologi- and downstream regulatory sequences of the genes, we also sought cal significance of the observed negative association between DEGs to examine methylation levels in the gene flanking sequences. Similar and DMRs, we investigated the functional categories of those genes. methylation levels were observed in the sequences upstream and Interestingly we found many DMRs showing negative association downstream of genes with high, moderate and low expression. with gene expression in fertile hybrid and 303 R lines, compared with However, contrary to this, genes with no expression exhibited high 11 A (sterile) line belonged to transposable elements and genes that methylation in their flanking regions suggesting that a high methyla- coded for PPR-associated proteins (Fig. 5E–G). In addition, we also tion in upstream and downstream regulatory sequences might be found genes coding for histone deacytylase and methyltransferase critical for suppressing transcript levels (Fig. 6A and Supplementary (Fig. 5E–G). This suggests that the transcription of those genes might Fig. S7). We also found that the presence of TE elements within the possibly be influenced by methylation. genes or in the vicinity resulted in a comparatively reduced expres- sion relative to the genes without TE elements within or in the neigh- boring regions (Fig. 6B). 3.7. Association between DNA methylation and transcript abundance Previous studies have reported a negative association between gene 4. Discussion expression and CHG and CHH methylation in plants like cotton. Rapid advancements in genome sequencing techniques have enabled Therefore to investigate the occurrence of any relation between genomic comparisons to understand genetic diversities between differ- changes in methylation levels and gene expression in pigeonpea, we ent plant species or different ecotypes of the same species. However, classified genes into four quartiles based on their FPKM values— understanding of epigenetic variation among different ecotypes has be- undetectable expression (FPKM 0-1), low expression (FPKM 1-10), come equally important to decipher its role in generating phenotypic moderate expression (FPKM 10–100) and high expression variations. A fundamental question of great importance is that how (FPKM> 100) according to the criteria described in earlier stud- 58,59 ies and correlated with changes in their methylation levels. the parental genome combination is regulated in hybrids to give rise to Consistent with previous studies, we found CHG and CHH methyla- significantly different hybrid genome which results in altered genetic traits. It is widely speculated that many agriculturally important phe- tion were negatively correlated with gene expression levels. nomena such as hybrid vigour and fertility restoration may arise due However, additionally, we also found negative association between CG methylation and gene expression (Fig. 6A and Supplementary to epigenetic changes between inbred parents and hybrid lines. Thus Fig. S7). Since transcription is importantly controlled by upstream the field of epigenetic research has expanded its horizon by adopting Figure 6. Correlation between gene expression and DMRs in pigeonpea. (A) Lineplots showing correlation between CG, CHG and CHH DMRs and the gene expres- sion which was divided into four quartiles—no expression (purple); low expression (green); moderate expression (red) and high expression (cyan) in F1 hybrid. (Also see Figure S7 for correlation in 11A and 303R parental lines). Note a negative correlation between gene expression and methylation in all the three contexts. (B) Box plots showing expression levels of genes inserted with and without TE within their gene bodies and 2kb upstream and downstream flanking regions. Gene þTE indicates insertion of transposon. Gene-TE indicates no insertion of transposon (Colored figure is available in the online version of the article). Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 371 many crop species as model systems. However, there are still agricul- Transposable elements can influence neighboring gene expression, turally important crops such as pigeonpea, in which existence of such contribute to inheritance of epi-alleles and thus represent the most mechanisms have remained unexplored. In this study, we investigated striking and common example of epigenetic inheritance in plants. and compared epigenetic changes in three different pigeonpea geno- Consistently, our data also revealed low expression of genes with in- types. We have generated whole genome methylation profiles of sterile serted TEs suggesting similar roles of TEs in pigeonpea. (11 A), fertile parents (303 R) and their fertile F1 hybrid, which were We also found a strong enrichment of DMRs in a wide spectrum unavailable till date. Our analysis of relative fractions of methyl- of genes that are associated with chromatin modification and have cytosines in each methylation context in the three pigeonpea lines re- been reported to regulate flowering, heterosis and reproduction in vealed highest fraction of mCs in CG context, followed by CHG and rice. We also found genes involved in DNA methylation that are CHH. Similar pattern of distribution of methylation fractions has known to play a role in transposon repression and genes involved 60 61 been observed previously in rice seedlings and Soybean leaves. We in histone methylation that help regulate flowering. also found that average methylation levels of the fertile hybrid and the DNA methylation is known to influence gene expression. fertile parental line (303 R) were lower than the sterile parent (11 A) However, as mentioned above, it has been proposed that DNA meth- suggesting that comparative low DNA methylation of genes may play ylation and H3K9 methylation in the gene body regions do not affect a role in restoring the fertility in the F1 hybrid. gene abundance. Despite increased CHG and H3K9 methylation in Mean methylation has been reported to be highest in CG contexts loss of function mutants of a H3K9 demethylase gene increased in 60 15 in other plants such as Rice and Arabidopsis, and followed bonsai methylation 1 (ibm1), the expression of their target genes re- closely by CHG and CHH. However, our results show that in mained unaffected. Also, absence of CG gene body methylation in pigeonpea, mean methylation levels in CHH context are much lower the loss of function MET1 mutants did not alter gene expression. than CG and CHG methylation levels. This suggests a possibility of Surprisingly, we found a negative correlation between transcript overall low expression of genes involved in regulation of CHH meth- abundance and methylation levels in all the three contexts in the ylations such as those involved in RdDM pathway that target homol- gene body suggesting a more prominent role of DNA methylation in ogous regions of the genome for methylation in pigeonpea regulating gene expression in pigeonpea. Previous studies in compared with Rice and Arabidopsis. Arabidopsis have reported that methylation in the promoter regions 71,72 Previously, it has been reported that transcription of genes is af- resulted in lower gene expression. Furthermore, another study in fected by methylation in the promoter region but not in the internal soybean revealed negative correlation between gene expression and regions. However, transposon silencing is regulated by methylation DNA methylation in their promoters. We also found that high both in the promoter region as well as internal regions. Our data methylation in the 5’ upstream regions of the genes were negatively demonstrate that fraction of mCs was higher in the transcribed re- correlated with gene expression levels, consistent with the idea that gions of transposons than the transcribed regions of the genes. As DNA methylation in the promoter regions affects gene abundance. shown in Figure 2, the fraction of methyl cytosines in upstream se- Our results in pigeonpea together with previous results in soybean quences of the transcribed regions of transposon was almost two open new avenues for understanding the relationship between gene times higher than fraction of mCs in upstream sequences of genes. expression and gene body DNA methylation and whether they play Additionally, gene bodies of TEs exhibited higher methylation than a direct role in restoring fertility in pigeonpea hybrids requires fur- their upstream and downstream regions. Contrary to this, genes had ther investigation. lower methylation in the transcribed regions than their flanking re- The raw reads of bisulfite and transcriptomic study were submit- 1,64 gions. As suggested previously for other crops, this could poten- ted to NCBI Sequence Read Archive under bioproject accession tially be a consequence of reduced RdDM pathway activity within PRJNA435649 and SRP133358 study. the gene body of the transcribed genes compared with TE. Furthermore, we also observed that the fraction of mCs in upstream Acknowledgements regions of transposons in fertile 303 R parent was slightly less than sterile 11 A parental line and the fertile hybrid. These observations We thank Council of Scientific and Industrial Research, Government of India suggest that one possible epigenetic mechanism of regulating trans- for the award of a senior research fellowship, ICAR-IARI, ICAR-NRCPB and poson inactivation in pigeonpea could be via methylation and this NRCPB in house project, New Delhi for supporting A.J. and providing other mechanism could be comparatively weaker in the fertile 303 R line funds and facilities. We also thank Dr N. Bhatia, The University of Sydney for than in the sterile parent and fertile hybrid. Over all, our results com- helpful feedback and ideas on the manuscript. ply with previously documented findings that CHH methylation lev- els are usually lower than CG and CHG, both on a genomic scale 65,66 Authors’ contributions and at specific loci. Further, our findings also reveal the existence of non-additive A.J and K.G conceived and designed the study. A.J. performed all the experi- methylation in F1 hybrid due TCM and TCdM events. TcM medi- mental work and data analysis. H.K. helped in data analysis and interpreting the results. A.J., K.G. wrote the manuscript. A.R. and N.K.S. helped editing ated increase in DNA methylation were more frequently observed in the manuscript. CG context and TcDM mediated decrease in DNA methylation were more frequently observed in CHH context, consistent with previous findings in Arabidopsis. Furthermore, we also dissected the mC Conflict of interest changes of differentially methylated parental epi-alleles in the hybrid, None declared. using SNPs for tracking the parental epi-alleles in the hybrid. Our analysis revealed that majority of the mC clusters having SNPs showed non- additive methylation in the hybrid when parents had Supplementary data different levels of methylation (HP or LP), consistent with our gen- eral mC cluster analysis. Supplementary data are available at DNARES online. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 372 The epigenomic interactions between parental inbreds 26. Chodavarapu, R. K., Feng, S., Ding, B., et al. 2012, Transcriptome and References methylome interactions in rice hybrids, Proc. Natl. Acad. Sci. U.S.A, 109, 1. Henderson, I. R. and Jacobsen, S. E. 2007, Epigenetic inheritance in 12040–5. plants, Nature, 447, 418–24. 27. Singh, N. K., Gupta, D. K., Jayaswal, P. K., et al. 2012, The first draft of the 2. Niederhuth, C. E. and Schmitz, R. J. 2014, Covering your bases: inheri- pigeonpea genome sequence, J. Plant Biochem. Biotechnol., 21, 98–112. tance of DNA methylation in plant genomes, Mol. Plant., 7, 472–80. 0 28. Saxena, K. B. and Singh, L. 1990, Low natural outcrossing in cleistoga- 3. Robertson, K. D. 2005, DNA methylation and human disease, Nat. Rev. 0 mous pigeonpea mutant, Int. Pigeonpea Newslett., 11, 9–10. ISSN Genet., 6, 597–610. 0255-786X. 4. Suzuki, M. M. and Bird, A. 2008, DNA methylation landscapes: provoca- 29. Chen, L. and Liu, Y. G. 2014, Male sterility and fertility restoration in tive insights from epigenomics, Nat Rev Genet., 9, 465–76. crops, Annu. Rev. Plant Biol., 65, 579–606. 5. Bird, A. 2007, Perceptions of epigenetics, Nature, 447, 396–8. 30. Tuteja, R., Saxena, R. K., Davila, J., et al. 2013, Cytoplasmic male 6. Law, J. A. and Jacobsen, S. E. 2010, Establishing, maintaining and modi- sterility-associated chimeric open reading frames identified by mitochon- fying DNA methylation patterns in plants and animals, Nat. Rev. Genet., drial genome sequencing of four Cajanus genotypes, DNA Res., 20, 11, 204–20. 485–95. 7. Vongs, A., Kakutani, T., Martienssen, R. A. and Richards, E. J. 1993, 31. Meshram, M. P. and Patil, A. N. 2018, Genetics of fertility restoration in Arabidopsis thaliana DNA methylation mutants, Science, 260, 1926–8., A2 cytoplasm based hybrids of pigeonpea. Int. J. Cur. Microbiol. Appl. 8. Kankel, M. W., Ramsey, D. E., Stokes, T. L., et al. 2003, Arabidopsis Sci., ISSN: 2319-7706 565-571. MET1 cytosine methyltransferase mutants, Genetics, 163, 1109–22. 32. Varshney, R. K., Chen, W., Li, Y., et al. 2011, Draft genome sequence of 9. Zemach, A., McDaniel, I. E., Silva, P. and Zilberman, D. 2010, pigeonpea (Cajanus cajan), an orphan legume crop of resource-poor farm- Genome-wide evolutionary analysis of eukaryotic DNA methylation, ers, Nat. Biotechnol., 30, 83–9. Science, 328, 916–9. 33. Xi, Y. and Li, W. 2009, BSMAP: whole genome bisulfite sequence 10. Cao, X. and Jacobsen, S. E. 2002, Role of the arabidopsis DRM methyl- MAPping program, BMC Bioinformatics, 10, 232. transferases in de novo DNA methylation and gene silencing, Curr. Biol., 0 0 34. Schultz, M. D., Schmitz, R. J. and Ecker, J. R. 2012, Leveling the playing 12, 1138–44. field for analyses of single-base resolution DNA methylomes, Trends 11. Mosher, R. A. and Melnyk, C. W. 2010, siRNAs and DNA methylation: Genet., 28, 583–5. seedy epigenetics, Trends Plant Sci., 15, 204–10. 35. Akalin, A., Kormaksson, M., Li, S., et al. 2012, methylKit: a comprehen- 12. Stroud, H., Do, T., Du, J., et al. 2014, Non-CG methylation patterns sive R package for the analysis of genome-wide DNA methylation profiles, shape the epigenetic landscape in Arabidopsis, Nat. Struct. Mol. Biol., 21, Genome Bio.l, 13,R87. 64–72. 36. Li, L. C. and Dahiya, R. 2002, MethPrimer: designing primers for methyl- 13. Reed, H. S. 1942, A short history of the plant sciences. https://www. ation PCRs, Bioinformatics, 18, 1427–31. scribd.com/document/190746815/Reed-Howard-S-1942-a-Short-History- 37. Kumaki, Y., Oda, M. and Okano, M. 2008, QUMA: quantification tool of-the-Plant-Sciences. for methylation analysis, Nucleic Acids Res., 36, W170–5. 14. Ng, S. Y., Lin, L., Soh, B. S. and Stanton, L. W. 2013, Long noncoding 38. Langmead, B. and Salzberg, S. L. 2012, Fast gapped-read alignment with RNAs in development and disease of the central nervous system, Trends Bowtie 2, Nat. Methods , 9, 357–9. Genet., 29, 461–8. 39. Li, B. and Dewey, C. N. 2011, RSEM: accurate transcript quantification 15. Greaves, I. K., Groszmann, M., Ying, H., Taylor, J. M., Peacock, W. J. from RNA-Seq data with or without a reference genome, BMC and Dennis, E. S. 2012, Trans chromosomal methylation in Arabidopsis Bioinformatics, 12, 323. hybrids, Proc. Natl. Acad. Sci. US. A, 109, 3570–5. 40. Robinson, M. D., McCarthy, D. J. and Smyth, G. K. 2010, edgeR: a bio- 16. Greaves, I. K., Groszmann, M., Wang, A., Peacock, W. J. and Dennis, E. conductor package for differential expression analysis of digital gene ex- S. 2014, Inheritance of trans chromosomal methylation patterns from pression data, Bioinformatics, 26, 139–40. Arabidopsis F1 hybrids, Proc. Natl. Acad. Sci. USA, 111, 2017–22. 41. Ye, J., Fang, L., Zheng, H., et al. 2006, WEGO: a web tool for plotting 17. Greaves, I. K., Gonzalez-Bayon, R., Wang, L., et al. 2015, Epigenetic GO annotations, Nucleic Acids Res , 34, W293–7. changes in hybrids, Plant Physiol., 168, 1197–205. 42. Liu, Y., Siegmund, K. D., Laird, P. W. and Berman, B. P. 2012, Bis-SNP: 18. Vaughn, M. W., Tanurdzic, M., Lippman, Z., et al. 2007, Epigenetic nat- combined DNA methylation and SNP calling for Bisulfite-seq data, ural variation in Arabidopsis thaliana, PLoS Biol., 5, e174. Genome Biol., 13,R61. 19. Schmitz, R. J., Schultz, M. D., Urich, M. A., et al. 2013, Patterns of popu- 43. Feng, S., Cokus, S. J., Zhang, X., et al. 2010, Conservation and diver- lation epigenomic diversity, Nature, 495, 193–8. gence of methylation patterning in plants and animals, Proc. Natl. Acad. 20. Shen, H., He, H., Li, J., et al. 2012, Genome-wide analysis of DNA meth- Sci. U.S.A., 107, 8689–94. ylation and gene expression changes in two Arabidopsis ecotypes and 44. Lisch, D. 2009, Epigenetic regulation of transposable elements in plants, their reciprocal hybrids, Plant Cell, 24, 875–92., Annu. Rev. Plant Biol., 60, 43–66. 21. Li, Y., Varala, K., Moose, S. P. and Hudson, M. E. 2012, The inheritance 45. Cokus, S. J., Feng, S., Zhang, X., et al. 2008, Shotgun bisulphite sequenc- pattern of 24 nt siRNA clusters in arabidopsis hybrids is influenced by ing of the Arabidopsis genome reveals DNA methylation patterning, proximity to transposable elements, PLoS One, 7, e47043. Nature, 452, 215. 22. Groszmann, M., Greaves, I. K., Albertyn, Z. I., Scofield, G. N., Peacock, 46. Qi, X., Li, M.-W., Xie, M., et al. 2014, Identification of a novel salt toler- W. J. and Dennis, E. S. 2011, Changes in 24-nt siRNA levels in ance gene in wild soybean by whole-genome sequencing, Nat. Commun., Arabidopsis hybrids suggest an epigenetic contribution to hybrid vigor, 5, 4340. Proc. Natl. Acad. Sci. U.S.A, 108, 2617–22. 47. The Arabidopsis Genome, I. 2000, Analysis of the genome sequence of the 23. He, G., Zhu, X., Elling, A. A., et al. 2010, Global epigenetic and tran- flowering plant Arabidopsis thaliana, Nature, 408, 796. scriptional trends among two rice subspecies and their reciprocal hybrids, 48. Zhang, J., Liu, Y., Xia, E. H., Yao, Q. Y., Liu, X. D. and Gao, L. Z. Plant Cell, 22, 17–33. 2015, Autotetraploid rice methylome analysis reveals methylation varia- 24. Barber, W. T., Zhang, W., Win, H., et al. 2012, Repeat associated small tion of transposable elements and their effects on gene expression, Proc. RNAs vary among parents and following hybridization in maize, Proc. Natl. Acad. Sci. U.S.A, 109, 10444–9. Natl. Acad. Sci. U.S.A, 112, E7022–9. 25. Parkin, I. A., Koh, C., Tang, H., et al. 2014, Transcriptome and methyl- 49. Li, X., Zhu, J., Hu, F., et al. 2012, Single-base resolution maps of culti- ome profiling reveals relics of genome dominance in the mesopolyploid vated and wild rice methylomes and regulatory roles of DNA methylation Brassica oleracea, Genome Biol., 15, R77. in plant gene expression, BMC Genomics, 13, 300. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018 A. Junaid et al. 373 50. Lu, Y., Liu, Y., Wang, J., Cao, M. and Rong, T. 2010, Variation and pat- 62. Matzke, M., Kanno, T., Daxinger, L., Huettel, B. and Matzke, A. J. 2009, terns of DNA methylation in maize C-type CMS lines and their main- RNA-mediated chromatin-based silencing in plants, Curr. Opin. Cell Biol., 21, 367–76. tainers, J. Plant Biochem. Biotechnol., 19, 43–50. 63. Inagaki, S. and Kakutani, T. 2012, What triggers differential DNA meth- 51. Chen, X., Hu, J., Zhang, H. and Ding, Y. 2014, DNA methylation ylation of genes and TEs: contribution of body methylation? Cold Spring changes in photoperiod-thermo-sensitive male sterile rice PA64S under Harb. Symp. Quant. Biol., 77, 155–60. two different conditions, Gene, 537, 143–8. 64. Zakrzewski, F., Schmidt, M., Van Lijsebettens, M. and Schmidt, T. 2017, 52. Schmitz, R. J., Schultz, M. D., Lewsey, M. G., et al. 2011, DNA methylation of retrotransposons, DNA transposons and genes in Transgenerational epigenetic instability is a source of novel methylation sugar beet (Beta vulgaris L.), Plant J., 90, 1156–75. variants, Science, 334, 369–73. 65. West, P. T., Li, Q., Ji, L., et al. 2014, Genomic distribution of H3K9me2 53. Becker, C., Hagmann, J., Muller, J., et al. 2011, Spontaneous epigenetic and DNA methylation in a maize genome, PLoS One, 9, e105267. variation in the Arabidopsis thaliana methylome, Nature, 480, 245–9. 66. Gent, J. I., Ellis, N. A., Guo, L., et al. 2013, CHH islands: de novo DNA 54. Eichten, S. R., Briskine, R., Song, J., et al. 2013, Epigenetic and genetic in- methylation in near-gene chromatin regulation in maize, Genome Res., fluences on DNA methylation variation in maize populations, Plant Cell, 23, 628–37. 25, 2783–97. 67. Martienssen, R., Barkan, A., Taylor, W. C. and Freeling, M. 1990, 55. Liu, F., Cui, X., Horner, H. T., Weiner, H. and Schnable, P. S. 2001, Somatically heritable switches in the DNA modification of Mu transpos- Mitochondrial aldehyde dehydrogenase activity is required for male fertil- able elements monitored with a suppressible mutant in maize, Genes ity in maize, Plant Cell, 13, 1063–78. Dev., 4, 331–43. 56. Huang, W., Yu, C., Hu, J., et al. 2015, Pentatricopeptide-repeat family 68. Shi, J., Dong, A. and Shen, W. H. 2014, Epigenetic regulation of rice protein RF6 functions with hexokinase 6 to rescue rice cytoplasmic male flowering and reproduction, Front. Plant Sci., 5, 803. sterility, Proc. Natl. Acad. Sci. U.S.A, 112, 14984–9. 69. Higo, H., Tahir, M., Takashima, K., et al. 2012, DDM1 (decrease in 57. Song, Q., Guan, X. and Chen, Z. J. 2015, Dynamic roles for small RNAs DNA methylation) genes in rice (Oryza sativa), Mol. Genet. Genomics, and DNA methylation during ovule and fiber development in allotetra- 287, 785–92. ploid cotton, PLoS Genet., 11, e1005724. 70. Choi, S. C., Lee, S., Kim, S. R., et al. 2014, Trithorax group protein 58. Xu, W., Yang, T., Dong, X., Li, D. Z. and Liu, A. 2016, Genomic DNA Oryza sativa Trithorax1 controls flowering time in rice via interaction methylation analyses reveal the distinct profiles in castor bean seeds with with early heading date3, Plant Physiol., 164, 1326–37., persistent endosperms, Plant Physiol., 171, 1242–58. 71. Zhang, X., Yazaki, J., Sundaresan, A., et al. 2006, Genome-wide 59. Lu, X., Wang, W., Ren, W., et al. 2015, Genome-wide epigenetic high-resolution mapping and functional analysis of DNA methylation in regulation of gene transcription in maize seeds, PLoS One, 10, arabidopsis, Cell, 126, 1189–201., e0139582. 72. Zilberman, D., Gehring, M., Tran, R. K., Ballinger, T. and Henikoff, S. 60. Garg, R., Narayana Chevala, V., Shankar, R. and Jain, M. 2015, 2007, Genome-wide analysis of Arabidopsis thaliana DNA methylation Divergent DNA methylation patterns associated with gene expression in uncovers an interdependence between methylation and transcription, Nat rice cultivars with contrasting drought and salinity stress response, Sci. Genet., 39, 61–9. Rep., 5, 14922. 73. Schmitz, R. J., He, Y., Valdes-Lopez, O., et al. 2013, Epigenome-wide in- 61. Song, Q. X., Lu, X., Li, Q. T., et al. 2013, Genome-wide analysis of DNA heritance of cytosine methylation variants in a recombinant inbred popu- methylation in soybean, Mol. Plant., 6, 1961–74. lation, Genome Res., 23, 1663–74. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/361/4939498 by Ed 'DeepDyve' Gillespie user on 22 August 2018

Journal

DNA Research – Oxford University Press

Published: Aug 1, 2018

There are no references for this article.