Pathogen-induced alterations in methylome during Lr48-mediated APR in wheat (Triticum aestivum L.)

Differential DNA methylation due to APR gene Lr48 for leaf rust resistance in wheat cv. CSP44 was examined at pre-adult plant (P-AP) susceptible stage and adult plant (AP) resistant stage. As many as 52,872 differentially methylated regions (DMRs) carrying 897 differentially methylated genes (DMGs) and many intergenic regions were identied. DMGs included exons/introns, promoters and UTRs. Little less than half of DMGs, were also compared with transcriptome data based on RNA-sEq. Susceptibility (at P-AP) was found to be governed by activation (due to hypomethylation) of relatively more genes, whereas resistance (at AP) involved silencing of relatively large number of genes. Using GO terms, DMGs were found to belong to a variety of biological processes including transcription regulation, protein synthesis, signal transduction, defense, photosynthesis, lipid and carbohydrate metabolism. The results of the present study improved our understanding of the epigenetic control of leaf rust APR in wheat.


Introduction
Leaf rust caused by the fungal pathogen Puccinia triticina Erikss & Henn. is one of the most devastating diseases of wheat 1 . Incorporation of adult plant resistance (APR) is one of the major strategies in breeding for durable resistance against this fungal pathogen. The APR genes provide enhanced defense leading to reduced compatible host-pathogen interactions resulting in durable resistance. In the past, this kind of resistance has been achieved using a number of APR genes including Lr34 2 and Lr46 3 for leaf rust and Sr2 for stem rust 4 . Disease resistance including APR is known to be partly genetic and partly epigenetic in nature 5 .
APR genes belong to a separate class of resistance genes, which differ from race-speci c R genes, which generally encode proteins carrying NBS-LRR domain 6 . The APR genes cloned so far have been found to confer resistance due to proteins other than those carrying NBS-LRR. These other proteins include cytoplasmic kinases 7 , ABC transporters 2 and hexose transporters 8 . A number of expression studies for wheat-rust pathosystem in the past suggested that the host rapidly and precisely re-program a large number of downstream signalling and defense genes for providing partial durable resistance against a broad range of fungal pathogens , 9,10,11,12,13 . This reprogramming may also involve epigenetic modi cations as shown in our own studies involving epigenetic modi cations (methylation, chromatin modi cation and ncRNA) of some downstream defense genes involved in leaf rust resistance mediated by race-speci c gene Lr28 and APR gene Lr48 12,13,14,15,16 .
Among epigenetic modi cations, DNA methylation has received maximum attention, and has been shown to be involved in regulation of expression of genes involved in a variety of cellular processes, which control development and stress responses 18 . This DNA methylation can be examined either at individual gene level or at the whole genome level. Methylation patterns suggest conservation of methylation in exonic regions (gene body) of genes also in addition to the well-known conservation of transposon methylation 19 . In higher plants, 5-methyl-cytosines, which accounts for methylation of > 30% UTRs (Supplementary Table S1). The level of methylation was highest at the transcription start site (TSS) relative to other anking regions in all the four pairs of comparison (Fig. 2).
The frequencies of DMRs on 21 different chromosomes also varied but, the pattern in each of the four comparisons was similar in different chromosomes. It is interesting to notice that the frequencies of DMRs in each chromosome was high in two pairs (S0 vs R0 and R0 vs R96) relative to those in the two other pairs (S0 vs S96 and S96 vs R96). Chromosomes 3D and 5D had relatively higher frequencies of DMRs, and the chromosome 4D had relatively fewer DMRs. A large number of DMRs could not be assigned to any of the 21 chromosomes. These DMRs were described as unassigned (Un); the relative abundance of these unassigned DMRs among the four pairs was still the same as in individual chromosomes ( Fig. 3; Supplementary Table S10).
DMRs were further classi ed into unique for each pair of comparison and those, which occurred in more than one pairs of comparisons. The frequencies of unique DMRs (hypomethylated + hypermethylated) was highest in the pair S0 vs R0 (14,386), followed by R0 vs R96 (9,789), S96 vs R96 (2,433)  The frequencies of hypermethylated and hypomethylated DMRs were high in S0 vs R0 (12,234 and 12,090) and R0 vs R96 (9,514 and 11,123) relative to those in other two pairs. (Range: 1,068 to 2,216; see Fig 6).

Differentially methylated genes (DMGs)
DMRs were used for identi cation of associated genes, leading to identi cation of 897 differentially methylated genes (DMGs). The frequencies of DMRs as well as DMGs in S0 vs R0 and R0 vs R96 was fairly high relative to those in the other two pairs of treatments (Fig. 5). Also, the relative frequencies of hyper-and hypomethylated DMGs in these two pairs were just opposite to each other, which is an important observation (Fig. 5). The pattern of DMGs in R0 vs R96 (442 hypermethylated and 154 hypomethylated) was just the opposite of that in S0 vs R0 (206 hypermethylated and 491 hypomethylated). However, among the remaining two treatment pairs, the difference in the number of hyper-and hypomethylated DMGs largely corresponded that in the DMRs. For instance, in S0 vs S96, 29 hypermethylated and 20 hypomethylated DMGs were identi ed. Similarly, in S96 vs R96, 54 hypermethylated and 32 hypermethylated DMGs were identi ed. (Fig. 5; Supplementary Tables S11-S18).

Transposable elements (TEs) in DMRs
Within DMRs, transposable elements (TEs) were identi ed, which included DNA transposons (DNA/CMC-EnSpm) representing the most widely occurring TE family, followed by LTR retrotransposons (gypsy and Copia); majority of the other classes were largely absent, barring some exceptions ( Table 3). The number of TEs in DMRs for S0 vs R0 and R0 vs R96 was relatively very high. The frequencies of TEs in different regions of the genome largely corresponded the frequencies of DMRs, such that a majority of TEs were present in intergenic regions, although a signi cant number of TEs was also present in the protein coding genes (Supplementary Tables S2-S5).

Functional annotation of DMGs
DMRs were used for identi cation of associated genes, leading to identi cation of 897 protein coding differentially methylated genes (DMGs). These genes were assigned to 12 different classes of functions ( Table 4) that were assigned on the basis of speci c domains in their encoded proteins (Supplementary Tables S11-S18). Further functional annotation was done using GO terms, which were available for only 651 of 897 DMGs, and were classi ed into three well-known broad classes, namely cellular components, molecular functions and biological processes. In each class, DMGs were also further assigned to different functions (Fig. 6).
The relative proportions of hypermethylated and hypomethylated genes differed not only among three broad classes, but also among different categories of genes in each of the three classes (Fig. 6a, 6b). For instance, while hypermethylated DMGs in large number belonged to incompatible interaction at APR involving the three treatment pairs, namely S0 vs R0, S96 vs R96 and R0 vs R96 (Fig 6a), the hypomethylated DMGs largely belonged to both the compatible (S0 vs S96) and incompatible interactions (R0 vs R96) and to some extent also to S96 vs R96. Very few hypomethylated DMGs were available in the pair S0 vs R0 (Fig. 6b).
In the class cellular components, generally high proportion of genes were hypermethylated in S96 vs R96 except the genes involving membrane and membrane part which were relatively more abundant in S0 vs R0. In contrast, relatively high proportion of hypomethylated genes in the class cellular components belonged to S0 vs S96 and R0 vs R96. In the class molecular functions, the proportion of hypermethylated genes was generally low except that for structural molecular activity, which was relatively high in case of in S96 vs R96. The proportion of hypomethylated genes belonging to membrane part was high in S0 vs S96, while that of genes belonging to catalytic activity and binding was high in R0 vs R96. For biological processes, the proportions of hypermethylated and hypomethylated DMGs were generally low except DMGs which were hypermethylated and hypomethylated in almost equal proportions for metabolic and cellular processes, although the identity of genes belonging to these two categories involved in hyper-and hypomethylation may differ.

A comparison of DMGs with DEGs
For a comparison of DMGs with DEGs, only 340 of the 897 DMGs were available in RNA-seq data.
As expected, for majority of DMGs, methylation had a negative relationship with expression level. In three treatment pairs, namely S0 vs S96, R0 vs R96 and S96 vs R96, hypermethylation was associated with reduced gene expression, whereas in the remaining treatment pair (S0 vs R0), reverse was true (hypomethylation associated with enhanced expression of genes). These genes showing expected relationship between methylation and gene expression, mainly encoded the following seven types of proteins: (i) ribosomal proteins, (ii) photosystem genes, (iii) cytochrome genes, (iv) ATP synthase, (v) Clp protease, (vi) NADH oxidoreductase and (vii) ribosomal proteins. There were also DMGs, which had unexpected positive relationship with gene expression level (Table 5;Supplementary Tables S19-22).

Quantitative real time PCR
In order to validate the relationship between methylation of genes and their expression, a set of 15 DMGs was selected for qRT-PCR involving three of the four treatment pairs (except S0 vs R0). Out of the 15 genes for which qRT-PCR analysis was performed, RNA-seq data was also available for 10 genes. Eleven among the above 15 genes used for qRT-PCR analysis followed the expected inverse pattern of methylation and gene expression (Fig. 7). All these eleven genes were methylated in the unassigned regions. However, six out of these eleven genes were such which were also methylated in the promoter regions besides unassigned regions whereas one gene was methylated in the UTR regions.
Brie y, the results of the eleven methylated genes showing inverse relationship with gene expression using qRT-PCR are described as follows: (i) R0 vs R96: ve hypermethylated genes encoding the following proteins had lower expression in R96 relative to R0; photosystem I, ATP dependent Clp protease, Acylcoenzyme-A oxidase and SAP domain, double stranded DNA binding protein; one hypomethylated gene encoding an uncharacterised protein showing higher abundance in R96. (ii) S0 vs S96: three hypermethylated genes encoding proteins for pentatricopeptide domain, ribosomal protein and ATP synthase had low expression in S96 relative to S0; and (iii) S96 vs R96: One hypermethylated gene encoding S4 RNA binding DNA protein and one hypomethylated gene encoding NADH quinone oxidoreductase showing reduced and higher expression in R96. (iv) The remaining four genes (out of 15) which showed inverse relationship of DNA methylation with gene expression data but not with qRT-PCR encoded the following proteins: NADH dehydrogenase, uncharacterised protein, ATP dependent Clp protease and cytochrome C-asm proteins.

Discussion
Disease resistance in plants including leaf rust resistance in wheat can be either all stage resistance (ASR) or adult plant resistance (APR). ASR provides race-speci c resistance during seedling stage to the adult stage, whereas APR provides resistance only at the adult stage. In wheat, APR genes for leaf rust have been classi ed by some in two groups, one group including 7 genes (Lr12, Lr13, Lr22a/b, Lr35, Lr37, Lr48 and Lr49) that are race-speci c and hypersensitive and the other group having 8 genes (Lr34, Lr46, Lr67, Lr68, Lr74, Lr75, Lr77 and Lr78) that are race-nonspeci c and non-hypersensitive , 41,42 . Other workers recognized only the 8 non-race-speci c non-hypersensitive genes as APR genes 43,44 .
The race non-speci c APR genes are also described as 'slow rusting' genes or 'partial resistance' genes, which typically reduce the disease development in adult plants. It has also been shown that while ASR is qualitative in nature following Flor's gene-for-gene relationship, APR is quantitative in nature, where an individual APR gene or QTL perhaps functions along with a number of minor genes. As many as ~80 QTLs were known for APR for leaf rust as in the year 2014 42 ; 170 more such genes must have been reported since then (manuscript under preparation). After recognizing the pathogen, like ASR, APR also initiates a defense response involving signal transduction and spatial/temporal differential expression of a large number of defense genes 43 . This aspect has not been studied in detail, although limited information is available 3,8,9,10,13,14 .
The gene Lr48 involved in the present study is a race-speci c hypersensitive APR gene 33 . In our earlier study, a candidate gene for Lr48 was also inferred, which was shown to encode a receptor like wall associated kinase (WAK) and was shown to be an orthologue of rice gene OsWAK8 13 . In order to understand the Lr48-mediated molecular events involved in downstream signal transduction and defense responses, transcriptome studies have also been conducted 8,9,13 . We believe that the expression of each APR gene as well as downstream signal transduction and defense responses might be partly regulated through epigenetic modi cations 44,45,46 . In particular, the role of DNA methylation in expression of APR genes for bacterial blight in rice was earlier examined, where a correlation was reported between DNA methylation (examined using MSAP) and repression of gene activity examined using Northern hybridization 44 .
The present study was focused on an analysis of the effect of the presence of the pathogen on genomewide DNA methylation in Lr48-mediated APR for leaf rust in wheat genotype CSP44. This allowed us to identify changes in differential methylome in different genic and intergenic regions. We assume that majority of the alteration in DNA methylation occurs of downstream genes, some of which encode proteins for LRR, photosystem I and II, Clp protease, peroxidases, glutaredoxin, kinases, transferases, etc.
While analysing the data, we noticed that the methylation level (hypo-as well hypermethylation) during normal growth and development (S0 vs R0) is rather high (Fig. 5, Table 2), which means that a fairly large number of genes undergo hypomethylation (Total 491, unique 130) and another set of relatively smaller number of genes undergo hypermethylation (total 206 DMGs, unique 99) during transition from P-AP to AP stage in the absence of pathogen. This methylation level is affected rather drastically due to the presence of the pathogen. A large number of differentially methylated regions (DMRs) were also available in each of the remaining three treatment pairs (S0 vs S96, R0 vs R96, and S96 vs R96), which suggest that the number of hyper-and hypomethylated DMRs is drastically reduced in the presence of the pathogen during susceptible P-AP stage but not during resistant AP stage ( Table 2; Fig.5). Similar results were also observed earlier in case of rice under drought stress, where the proportion of DMRs induced due to drought stress were remarkably low in sensitive line relative to a tolerant line 28 . However, the situation of DMGs differs from the situation of DMRs. It is apparent that very few genes are differentially methylated during P-AP stage (S0 vs S96) relative to AP stage (R0 vs R96) (Fig. 5). At the resistant AP stage (R0 vs R96), the methylation status of DMGs is just the opposite of the condition observed during normal growth and development (S0 vs R0). At the resistant AP stage, the frequencies of hypomethylated DMGs is greatly reduced while that of hypermethylated DMGs is greatly increased. This appears to be a signi cant observation suggesting that fewer DMGs are activated relative to a large number of DMGs, which are silenced. Some of the important genes which are activated due to hypomethylation include Fbox, wall associated kinase (WAK), glutaredoxin, thioredoxin, armadillo type fold, NB-ARC, LRR, etc. Similarly, the genes which are silenced due to hypermethylation encode proteins for photosystem I and II, ATP synthase, NADH-quinone oxidoreductase, cytochrome C oxidase, etc.
The genes that are activated due to hypomethylation may be perhaps involved in positive regulation of The genes that are silenced or downregulated due to hypermethylation may be perhaps involved in negative regulation of Lr48 mediated APR so that silencing of these genes will facilitate APR. There are very few examples of genes that are involved in negative regulation of disease resistance. Therefore, this seems to be a new information, which may need validation in future. For instance, the genes involved in photosynthesis like cyt b6/f of photosystem II are earlier reported to induce resistance through negative regulation. These genes exhibited reduced expression in response to P. syringeae infection in soybean leading to resistance involving hypersensitive reaction. This resistance was induced due to damage of PSII system, which leads to oxidative burst leading to expression of other defense related genes 54 . In fact, the repression or downregulation of genes involved in photosynthesis during biotic stress like leaf rust in the present study was earlier shown to be an adaptive response to biotic stresses 55 . Our results are also largely in agreement with the results of an earlier study in Arabidopsis, where DMGs associated with a large number of DMRs were shown to express in response to the bacterial pathogen, Pseudomonas syringae pv. tomato DC3000 (Pst) 18 .
In the present study, 229 unique DMGs (130 hypo-and 99 hyper-) were observed in the treatment pair S0 vs R0 indicating the role of these genes mostly in transition from P-AP to AP stage during normal growth and development (Fig. 4). The methylation status of majority of these genes changes during susceptible P-AP stage (S0 vs S96); these genes encode proteins for peptidases, transferases, zinc nger proteins, glutaredoxins and thioredoxins, SWEET sugar transporters, expansins, heat shock proteins, FAD binding domains, etc. The role of some of these genes encoding proteins like expansins, FAD binding domains, zinc nger proteins in growth and development has been widely recognized 56,57,58 .
There are only a few genes (20 hypo-and 29 hypermethylated) which are still differentially methylated during pathogen invasion. These genes include those encoding proteins ferredoxins, aldehyde oxidase, ribosomal proteins, methionine-S-methyl transferase, protein kinases, LRR, etc. The role of these genes during susceptibility needs further examination.
The above results and interpretation receive support from GO analysis. For instance, during susceptibility (S0 vs S96), hypo-and hypermethylated DMGs were found to be different involving completely different functions. These included multi-organ process, reproduction, signalling and biological regulation for hypermethylated DMGs and cellular and metabolic processes for hypomethylated DMGs. On the other hand, the hypermethylated genes during the resistant AP stage were found to be largely involved in structural molecule activity, metabolic and cellular processes, whereas the hypomethylated DMGs during resistance were involved in reproduction, signalling, response to stimulus, etc.
RNA-seq data was also available for approximately 38% (340/897 genes) of all the DMGs identi ed in the present study 13 . For 312/340 DMGs (Supplementary Tables S19-22), the expression pattern matched with the level of methylation (assuming that hypermethylation was associated with repression of gene expression and hypomethylation was associated with enhanced gene expression). The results of methylation induced changes in some representative genes were also validated through the results of qRT-PCR analysis, where the expression of 11 out of the 15 genes examined showed inverse relationship with DNA methylation. In an earlier study in rice also, under drought stress, level of methylation in more than half of the DMGs showed negative relationship with gene expression 28 . Further, in the present study, 76 of the above 312 DMGs, exclusive methylation of the promoter region suppressed gene expression. It is known that DNA methylation in promoter directly represses transcription in one or more of the following ways: (i) inhibition of the binding of transcription activators, (ii) promotion of the binding of transcription repressors, (iii) indirect repression of transcription by promoting repressive histone modi cations such as H3K9me2, and (iv) inhibiting permissive histone modi cations such as histone acetylation 59,60 .
Another interesting observation of the present study is an abundance of DMRs in intergenic regions, relative to DMRs in genic regions, which included exons/introns, promoters and other regulatory regions (Table 2). This also received support from two recent studies where >70% of the DMRs were found in the intergenic regions. These included one of our own study involving MeDIP analysis for Lr28-mediated ASR 16 and a study involving WGBS during lead (Pb) stress in maize 31 . The intergenic regions are known to carry many enhancers and other regulatory regions like non coding RNA genes and TEs. Therefore, it is possible that intergenic DNA methylation involves methylation of enhancers and other regulatory sequences encoding non-coding RNAs for trans-regulation as observed in case of genes controlling precursor B-cell acute lymphoblastic leukaemia in mammals 61 . Similar results were available in our earlier study, where we examined the role of histone modi cations in leaf rust resistance 17 . It is also known that the effect of methylation in gene body and the associated 3' UTR and 5' UTR regions differs from the effect of methylation in the promoter sequences. DNA methylation in the gene body region is known to promote transcription involving alternate splicing 62 , while that in the promoter regions inhibits transcription 48,63 . The results of the present study are in agreement with this general observation.
During the present study, TEs were also examined in the DMRs, although speci c methylation of TEs could not be examined. The number of TEs (located in the gene body, promoter and intergenic regions) was higher in the DMRs at AP stage (R0 vs R96) than the P-AP stage (S0 vs S96) ( Table 3; Supplementary  Tables S2-S5). This suggested increasing role of possible methylation of TEs in silencing of genes with passage of time from P-AP susceptible stage to AP resistant stage. These results are in agreement with the results of an earlier study 64 . Among different types of TEs, the LTR retrotransposons like GYPSY, COPIA and Em-spm were abundant among DMRs relative to other types of transposons, particularly in S0 vs R0 and R0 vs R96. Occurrence of TEs and their activity has also been shown to be related with methylation patterns in genomes of Arabidopsis and cassava 65,66 ; GYPSY and COPIA were more heavily methylated than other TEs in cassava 67 .
The results of DNA methylation at the P-AP and AP stages involving APR gene Lr48 in wheat genotype CSP44 were also compared with methylation status in susceptible and resistant NILs for ASR gene Lr28 in wheat cultivar HD2329 16,17 . The results revealed a similar pattern of methylation change, i.e. hypermethylation mediated gene silencing in resistant reaction indicating some common epigenetic mechanism involving methylation during ASR and APR. Further observations also revealed two common genes encoding NADPH oxidoreductase and cyt C assembly which were hypomethylated in the present study and also showed differential binding sites for H3K27me3 (a repression mark) in an earlier ChIP-seq analysis 17 . Similarly, while comparing the differentially methylated regions obtained during ASR due to Lr28 in our earlier study and APR due to Lr48 in the present study, 17 hypermethylated genes and 12 hypomethylated genes were found to be common. These genes encoded proteins for photosystem I, cyt C, cyt b559, different ribosomal proteins, etc. suggesting the role of regulated expression of these genes in disease resistance.

Summary And Conclusions
In the present study, the role of DNA methylation in APR was examined in hexaploid wheat genotype CSP44, which carries an APR gene Lr48. The status of methylation after inoculation was examined at two P-AP susceptible stages (S0 and S96) and two AP resistance stages (R0 and R96). Genome-wide methylation was examined using MeDIP approach, which is cost-effective. Methylation data was used for identi cation of differentially methylated regions (DMRs) of the genome, which were then utilized for identi cation of differentially methylated genes (DMGs), belonging to different classes determined on the basis of gene ontology (GO). The data on DMRs and DMGs was also compared with data on differentially expressed genes (DEGs) identi ed through RNA-Seq mediated transcriptome analysis. The results suggested that susceptibility involves activation (due to hypomethylation) and repression (due to hypermethylation) of relatively fewer but different sets of genes, while resistance involves hypermethylation-mediated silencing of many genes and hypomethylation-mediated activation of fewer genes. In particular, during susceptibility, some genes (e.g. genes encoding clp protease, acyl co-A oxidase, photosystem I and II, etc.) became activated following hypomethylation; there were other genes (glyceraldehyde 3 phosphate, ATPase, serine threonine like kinase, etc.), which were silenced following hypermethylation. Similarly, during resistance, some genes (wall associated kinase, leucine rich repeats, peroxidases, cytP450, etc.) which are known to be involved in resistance are activated whereas some others (NADH ubquinone oxidoreductase, photosystem II, pyridoxal phosphate, etc.) were silenced. The methylation results were compared with results of expression studies involving transcriptomics using RNA-Seq data, which largely con rmed that high expression of genes involves hypomethylation and vice versa. The results of the study were validated using qRT-PCR, although there were examples, where qRT-PCR results deviated from expectation. However, DNA methylation is not the only epigenetic phenomenon and there are other genetic and epigenetic mechanisms for gene regulation. We believe that the present study will improve our understanding of the genetics and epigenetics of APR.

Plant material and inoculation
The genotype CSP44, which is a selection from Australian cv. Condor (WW80/2 * WW15 selection) and carried a recessive APR gene (Lr48) was used in the present study. For inoculation, a single spore culture of Puccinia triticina pathotype 77-5 (121R63-1), which is the most predominant and devastating pathotype in all parts of the Indian subcontinent, was obtained from Regional Station, ICAR-Indian

Institute of Wheat & Barley Research, Flowerdale, Shimla.
Seeds of CSP44 were sown in pots containing sterilized medium consisting of agropeats, vermiculite and sand in the ratio of 2:1:1. The plants were raised in a growth chamber under controlled conditions of 16h light (450μmol m -2 s -1 light intensity) with 25 o C and 8h dark with 18 o C at the National Phytotron Facility, ICAR-IARI, New Delhi. Inoculation of leaves was carried out following the procedure described elsewhere 34 . Brie y, the plants were inoculated by spraying urediniospores suspended in water forti ed with Tween-20 (0.75mL mL -1 ) at an average concentration of 20 urediniospores/microscopic eld (10X x 10X). Four sets of plants were used for collecting the material at four different stages 36 (Fig. 1). Two sets were inoculated; one at S0 (P-AP; GS39) and the other at R0 (AP; GS49). Material was collected as follows: (i) material for S0 was collected from one of the two sets which were not inoculated; (ii) material for R0 was collected from the other set which was also not inoculated; (iii) material for S96 was collected from the set where inoculation was done at S0; (iv) material for R96 was collected from the fourth set in which inoculation was done at R0. The infection type was recorded during 12-14 days after inoculation using the 1-4 rust infection scale of Stackman et al. (1962) 37 .

DNA isolation and library preparation for MeDIP-seq
Leaf samples for DNA isolation were collected at P-AP and AP stages, each at 0 h (mock-inoculated samples) and 96 h after inoculation (hai), thus constituting four treatments (with three biological replicate per treatment). These four treatments were named as follows: (i) susceptible, mock (S0), (ii) susceptible, 96 hours after inoculation (hai) at P-AP stage (S96), (iii) resistant, mock (R0) and (iv) resistant, 96 hours after inoculation at AP stage (R96). Genomic DNA (50-100 ng) from all the four treatments was used for immunoprecipitation according to manufacturer's protocol (Zymo Research; EZ DNA Methylation-Gold kit). In order to reduce the sequencing cost, DNA samples of three biological replicates for each of the four treatments were pooled in equal concentrations. These pooled DNA samples were used for preparation of MeDIP-seq libraries using the protocol of Li et al. (2010) 38 with a few modi cations. Brie y, following steps were involved: (i) Genomic DNA (1 μg) from each of the four treatments was sheared independently and end-repaired using End Repair Mix with 3' to 5' exonuclease activity to obtain sheared fragments with blunt ends. (ii) A single 'A' nucleotide was added to the blunt ends of DNA fragments, to prevent the ligation of 3′-ends of the blunt fragments to one another during the adapter ligation reaction. Similarly, a single 'T' nucleotide was added to the 3′-end of the adapter to generate a complementary overhang for ligating the adapter to the fragment. This strategy ensured a low rate of chimera (concatenated template) formation. (iii) The multiple indexing adapters were then ligated to the ends of the DNA fragments to prepare them for hybridization on to the ow cell. (iv) In the next step, immunoprecipitation of the DNA was carried out using 5mC speci c monoclonal antibody and MeDIP buffer (Methylamp Methylated DNA capture kit, Epigentek, USA). (v) PCR was performed with primers having the adaptor barcodes/indices. The number of PCR cycles was kept to the minimum to retain library representation. (vi) Qubit facility was used for quality check of the library; this was done by loading 1 μl of sample on a 2100 Bioanalyzer (Agilent Technologies) using a DNA-speci c chip (e.g., Agilent DNA 1000 or DNA HS chip). The MeDIP libraries thus obtained were outsourced to SciGenome Pvt. Ltd. India for sequencing on Illumina HiSeq2000 platform to obtain 2x100bp paired end reads. The data (in the form of reads) was obtained in the form of FASTA les. The MeDIP-seq raw reads of P-AP susceptible and AP resistant stages were deposited in Sequence Read Archive (SRA), NCBI and can be accessed through Bioproject PRJNA693113 under the accession numbers SAMN17377447, SAMN17377526, SAMN17377527, and SAMN17377578.

MeDIP-seq analysis
The raw sequence reads, procured as above, were processed to obtain clean reads, which were assessed for their quality using FastQC tool using the following parameters: sequence quality Phred score (>30); read length distribution, 100 bp and GC distribution of reads (>40-60%). The clean reads were then aligned to the reference genome of wheat (Triticum_aestivum IWGSC_ensembl_release 48) using bowtie program (version 2.0) with default parameters (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml).
The uniquely mapped reads (.bam les) were analysed using Bioconductor package involving quality control of quantitative data derived from sequences. In order to normalize CpG density bias due to CpG coupling, a control sample was used as a reference 24 . Saturation analysis of the mapped reads was performed to determine su cient coverage pro le of the reference genome. The uniquely aligned reads were used for estimating methylation coverage using MEDIPS R package. The fractions of methylations (CpG contexts) were extracted using MethylDackel software (https://github.com/dpryan79/MethylDackel) with fraction option value set between 0 and 1. The extracted data was transformed to a bigWig le using bedGraph to BigWig script. ComputeMatrix option was used to calculate scores at reference point TSS using bigwig score le. The matrix thus generated was employed to draw TSS plot that evaluates methylation fraction across all transcription start sites.
The TSS plots were also used to compute the methylation coverage for a region covering 2000 bp downstream and upstream of the TSS.

Identi cation of differentially methylated regions (DMRs)
The differentially methylated regions (DMRs) in the four pairs of treatments (S0 vs S96, S0 vs R0, R0 vs R96 and S96 vs R96) were identi ed using Bioconductor edgeR Package at a genome wide sliding window of 100bp resolution, using a p-value ≤ 0.05 and log2FC ≥ 2 or ≤ −2. This analysis was based on CpG count that was treated to represent methylation events. In all these four comparisons, the second treatment was compared to the rst. The DMRs were assigned to different genomic regions, which included genes (exons, introns), intergenic regions, regulatory regions which included regions up to 2kb upstream (5' UTRs and promoters) and 2kb downstream (including 3' UTRs) and also regions, which could not be assigned to any region and were therefore described as unassigned.
Functional annotation, Gene ontology (GO) analysis and identi cation of transposable elements.
The differentially methylated genes (DMGs) were annotated using InterPro database and categorized into different classes. GO analysis was carried out using Blast2GO software 39 and differentially methylated genes (DMGs) were assigned to three well-known regions including cellular localizations, molecular functions and biological processes. WEGO 2.0 version in native text format was used for plotting histograms 40 . In order to identify transposable elements (TEs), including both LTR (long terminal repeats) and non-LTRs, the nucleotide sequences of all the DMRs were used as input in the online tool RepeatMasker (http://www.repeatmasker.org).

Comparison of DMGs with DEGs (based on RNA-seq data)
The presence of methylation in the differentially expressed genes (DEGs), was examined using RNA-seq data from our earlier study 13 . Fold changes in methylation of the genic regions (intron, exon, promoter and TTS) of DMGs were compared with fold changes in expression of the corresponding DEGs. The fold change in methylation levels of genes showing differential methylation in more than one out of the above four genic regions were merged; in each case, average fold change in methylation was compared with corresponding fold change in expression that was available from RNA-Seq data.
For expression analysis using qRT-PCR, 15 DMGs were selected; for 10 of these 15 genes, RNA seq data was also available. Following criteria were used for selection of these 15 genes: (i) occurrence of domains involved in biotic/abiotic stress tolerance in the proteins encoded by these genes; (ii) ≥ 2 or ≤ −2-fold difference in methylation level; (iii) hypo/hyper-methylated genes from at least three out of the four treatment pairs (10 genes for R0 vs R96, 2 genes for S96 vs R96 and 3 genes for S0 vs S96). For S0 vs R0, no genes could be studied as the three primers that were synthesised for this treatment pair either failed to show expected amplicon or gave ambiguous results.
Leaf samples for qRT-PCR were collected at S0, S96, R0 and R96 stages, as detailed above. The collected leaf samples were frozen in liquid nitrogen and stored at −80 °C. For each treatment in each experiment, three replications were used. Total RNA was isolated from frozen leaf tissue using the TRIzol Reagent (Ambion) as per the manufacturer's speci cations and treated with RNase-free DNase I (Invitrogen) for 15 min to eliminate any residual genomic DNA. First-strand cDNA was synthesized from DNase I-treated total RNA using Verso cDNA Synthesis Kit (Thermo Scienti c Inc., USA) according to the manufacturer's instructions.
The primers for qRT-PCR analysis were designed using Primer Quest program from IDT Technologies, which are listed in Supplementary Table S1. The qRT-PCR reactions were carried out in three technical replicates per biological replicate according to the following conditions: (i) one cycle for initial denaturation at 95°C for 30 sec, (ii) 40 cycles, each with 95°C for 5 sec, 60°C for 34 sec; the gene expression was calculated using 2 -ΔΔCT method. Wheat 18S RNA gene was used as an endogenous control. For each primer pair, a negative control (water) was also included. The signi cance of difference between the expressions of genes in the two treatments was tested using t-test.  *The regions where the methylation in genes could not be assigned to any of the ve regions (exon, intron, 3' UTR, 5' UTR and promoter) are designated as unassigned.