N Cui ,Xiofeng Chen ,Yn Shi ,Meirong Chi ,Jintun Hu ,Kunlong Li ,Zonghu Wng ,b,c,*,Hifeng Wng ,d,*

a State Key Laboratory of Ecological Pest Control for Fujian and Taiwan Crops,College of Plant Protection,Fujian Agriculture and Forestry University,Fuzhou 350002,Fujian,China

b Marine and Agricultural Biotechnology Laboratory,Institute of Oceanography,Minjiang University,Fuzhou 350108,Fujian,China

c Fujian Province Key Laboratory of Pathogenic Fungi and Mycotoxins,Fujian Agriculture and Forestry University,Fuzhou 350002,Fujian,China

d State Key Laboratory for Conservation and Utilization of Subtropical Agro-Bioresources,Guangxi Key Lab of Sugarcane Biology,College of Agriculture,Guangxi University,

Nanning 530004,Guangxi,China

ABSTRACT DNA methylation participates in regulating the expression of coding and non-coding regions in plants.To investigate the association between DNA methylation and pathogen infection,we used whole-genome bisulfite sequencing to survey temporal DNA methylation changes in rice after infection with the rice blast fungus Magnaporthe oryzae.In contrast to previous findings in Arabidopsis,global DNA methylation levels in rice increased slightly after rice blast infection.We identified over 38,000 differentially methylated regions (DMRs),and hypermethylated DMRs far outnumbered hypomethylated DMRs.Most DMRs were located in transposable element regions.Using transcriptome analysis,we identified 8830 differentially expressed genes (DEGs) after 1,3,and 5 days of infection.Over one-third of DEGs,most of which were CHH-type DMRs,were associated with DMRs.Functional analysis of the CHH DMR-DEGs indicated their involvement in many important biological processes,including cell communication and response to external stimulus.The transcription of many NBS-LRR family genes was affected by changes in DNA methylation,suggesting that DNA methylation plays essential roles in the response of rice to M.oryzae infection.More broadly,the DNA methylation analysis presented here sheds light on epigenomic involvement in plant defense against fungal pathogens.

Keywords:DNA methylation Epigenome Magnaporthe oryzae Rice Transcriptome

1.Introduction

DNA methylation plays roles in transcription regulation and the silencing of repetitive elements,such as transposable elements(TEs).In eukaryotes,many biological processes depend on regulation of DNA methylation levels.In contrast to animals,in which methylated cytosines are found predominantly in the CG context,methylated cytosines in plants are present in all three sequence contexts:CG,CHG(H being A,C or T),and CHH[1].Although most knowledge of plant DNA methylation comes from Arabidopsis,progress in elucidating the mechanism of DNA methylation has been made in rice (Oryza sativa),including characterization ofde novoand maintenance-methylation enzymes [2].As in Arabidopsis,the establishment of all three sequence contexts of DNA methylation in rice is also mediated by small RNAs via an RNA-directed DNA methylation(RdDM)pathway[3].In rice,OsDRM2,a homolog of Arabidopsis DRM2 (DOMAINS REARRANGED METHYLASE 2),is responsible forde novomethylation in CG and non-CG contexts in the rice genome.In Arabidopsis,MET1 (METHYLTRANSFERASE 1),is the major CG methylase responsible for CG methylation maintenance [4].In rice,twoMET1homologs (OsMET1aandOsMET1b)have been identified,andOsMET1bis the major methyltransferase gene,given that a loss-of-function mutation inOsmet1bled to a greater reduction in CG methylation than that inOsmet1acompared with the wild type[5].As in Arabidopsis,there are threeCMTgenes in the rice genome:OsCMT2,OsCMT3a,andOsCMT3b.OsCMT3ais the major gene and is involved in the maintenance of CHG sites.Compared with theOsCMT3amutation,theOscmt3bmutation does not lead to abnormal morphology and is expressed only in panicles [6].ForOsCMT2,although there are few studies of its function,based on its similarity with ArabidopsisCMT2it is speculated to play a critical role in CHH methylation maintenance.Compared with the Arabidopsis genome,the rice genome is larger and contains more repetitive elements,such as transposable elements(accounting for more than 40%of the genome).Transposons in rice are widely spread through the genome,with a majority associated with genes.Because transposons are the main target sites for DNA methylation,this difference may partly explain the stronger effects of rice DNA methyltransferases than those in Arabidopsis.

Rice is a staple food for half of the human population.The DNA methylation rate in rice is 24.3%of all covered cytosines,four times higher than that in Arabidopsis [7,8].Previous studies [9,10] indicated that rice could easily gain or lose DNA methylation during either tissue culture or regeneration.In rice endosperm,the maternal genome is locally hypomethylated,helping to establish specific gene expression patterns for seed development[11].Besides developmental process,DNA methylation also plays a role in rice response to abiotic stress.Environment-responsive genes such as drought-response genes are likely to be regulated by DNA methylation [12].Phosphate starvation causes gene transcription changes and DNA methylation variation,part of which can be transmitted through mitosis [13].Although there is much epigenetic evidence for the role of DNA methylation in rice [2],little research has dealt with rice under biotic stress.Only a recent study of RNA polymerase IV in rice suggested a possible role of DNA methylation in rice–virus interaction [14].In view of the agricultural impact of plant pathogens,a detailed understanding of DNA methylation mechanisms and changes in crop plants under pathogen attack is desirable.

Rice blast,caused by the fungal pathogenM.oryzae,is a major disease that can reduce annual rice harvests by 10%–30% [15,16].Many studies of host–pathogen interactions have investigated the molecular mechanisms of eitherM.oryzaevirulence or rice resistance [17].Most current work concentrates on specific virulence genes or resistance gene families,with few studies focusing on rice epigenomic changes during pathogen infection.

The object of this study was to measure dynamic DNA methylation changes in rice(O.sativassp.indica)after infection withM.oryzaeby using whole-genome bisulfite sequencing,and surveyed temporal gene expression across different infection stages.This comparative approach revealed local specific DNA methylation changes associated with rewiring of the transcriptional networks in response to pathogeninfection in rice.This work demonstrates a critical role of DNA methylation in pathogen infection response of rice.

2.Materials and methods

2.1.Plant and fungal strains

TheM.oryzaestrain Guy11 was used for infection assays.It was cultured on starch–yeast medium (SYM) containing 0.2% (w/v)yeast extract,1.0% (w/v) starch,0.3% sucrose and 1.5% agar.For conidiation,the strain was cultured on rice bran medium containing 2.0% rice bran and 1.5% agar (pH 6.0).Conidial suspensions(1 × 105conidia/mL in 0.02% Tween 20 solution) of Guy11 were used for inoculation assays on 3-week-old rice seedlings of the susceptibleO.sativassp.indicacultivar CO39 as described previously[18].Rice leaves inoculated withM.oryzae0,1,3,and 5 days were collected as four samples and each sample had two biological replicates.

2.2.Whole genome bisulfite sequencing and analysis

Genomic DNA was extracted from all samples before bisulfite treatment.After library construction,DNA was sequenced using the Illumina HiSeq2500 platform.Low-quality reads were removed from the raw data to leave clean reads.TheO.sativassp.indicaASM465v1 assembly was downloaded from Ensembl Plants(http://plants.ensembl.org/Oryza_indica/Info/Index?db=core)as the reference genome.Clean reads were mapped to the reference genome using BSMAP v2.90 [19] with a mismatch rate of up to 4%,and only uniquely mapped reads were kept for subsequent analysis.Cs covered by over four reads were used to compute individual cytosine methylation ratios according to the following formula:Cs/(Cs+Ts).To estimate the repeatability,the genome was divided into 2-kb windows,the methylation level of each window was determined,and the Pearson correlation coefficient between the two replicates were calculated.The average methylation levels were calculated as weighted methylation levels in terms of CG,CHG,and CHH[20].Circos-0.69-6(http://circos.ca/software/download/) was used for drawing circos plot.For DMR detection,the R package ‘‘methylKit” (http://bioconductor.org/packages/methylKit/)was employed,and the genome sequence was divided into 100-bp bins.The absolute methylation level differences between test and control were required to be at least 0.4,0.2,and 0.1 for CG,CHG,and CHH,respectively.To identify DMRs that overlapped genes or TEs,bedtools (https://github.com/arq5x/bedtools2) was used with an overlap region of >90% as a cutoff value.

2.3.qPCR validation

To study the expression patterns of methylation-related genes in rice during infection byM.oryzae,rice leaf samples were collected at 0,1,3,and 5 days post-inoculation (dpi).Total RNA was extracted using the Eastep Super Total RNA Extraction Kit (Promega,Madison,[WI] U.S.A.),and first-strand cDNA was synthesized using the PrimeScript RT Reagent Kit with gDNA Eraser(TaKaRa Bio Inc.,Otsu,Shiga,Japan) according to the manufacturer’s instructions.Quantitative real-time RT-PCR (qRT-PCR) and SYBR Premix Ex Taq(TaKaRa Bio Inc.,Otsu,Shiga,Japan)were used to quantify the relative transcript levels,and the actin gene was used as internal control for normalization.All reactions were performed for three independent biological experiments with at least three replicates,and the relative transcript level was calculated using the 2-ΔΔCt method [21].The gene-specific primers used for qRT-PCR are listed in Table S1.

2.4.TE annotation and analysis

TEs were re-annotated with RepeatMasker v4.0.5 (http://www.repeatmasker.org)on the reference genome with the repeat library from RepBase (v20170127).For miniature inverted-repeat transposable element (MITE) annotation,968 TE seed sequences with target-site duplications(TSDs)inO.sativassp.japonicacv.Nipponbare were downloaded from the P-MITE database (http://pmite.hzau.edu.cn/).These seed sequences were queried against the reference genome using RepeatMasker.Only repeats <15% divergent from the seed sequences were marked as MITEs.These MITEs were excluded from other TE superfamilies using bedtools with >80%overlap.For each TE superfamily,the expected value was calculated as previously [22] described.

2.5.RNA sequencing and analysis

Total RNA was extracted as described above.First-strand cDNAs were synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (New England BioLabs,Ipswich,MA,USA)and sequenced on the Illumina HiSeq2500 platform.After removal of low-quality reads,the clean reads were mapped to the reference genome using Hisat2 v2.1.0 (https://daehwankimlab.github.io/hisat2/).Only uniquely mapped paired-end reads were kept for read counting with HTSeq-count v0.7.2(https://htseq.readthedocs.io/en/release_0.9.1/index.html).Reads with alignment quality lower than 5 were removed.For expression values and the identification of differentially expressed genes (DEGs),the R package‘‘NOISeq”v2.24.0 (http://www.bioconductor.org/packages/devel/bioc/html/NOISeq.html) was used with default parameters.Trimmed Mean of M (TMM) values were calculated to represent gene expression levels.To estimate the repeatability of the transcriptome,TMM values were log2transformed and a Pearson correlation coefficient were calculated based on log2TMM values across two biological replicates.DEG identification was based on the following cutoff values:fold change >2 and probability >0.95.The probability of 0.95 in NOISeq software is equivalent to an FDR of 0.05.R package ComplexHeatmap (http://www.bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html) was employed for heatmaps and K-means clustering.The line plots were drawn with dplyr (https://www.rdocumentation.org/packages/dplyr/versions/0.7.8).Gene Ontology (GO) annotation and enrichment analysis of the DEGs was performed using AgriGO v2(http://systemsbiology.cau.edu.cn/agriGOv2/index.php),and aPvalue <0.05 was used to identify enriched GO terms.

3.Results

3.1.Global DNA methylation patterns in rice in response to pathogen infection

To investigate whether global DNA methylation is altered in rice in response to pathogen infection,we performed whole genome bisulfite sequencing(termed BS-seq or Methyl-seq)of rice infected withM.oryzaefor 1,3,and 5 days and uninfected rice as a control.DNA were extracted from leaves of each stage sample and sequenced with two biological replicates per sample to assess repeatability.Each library was independently mapped to the cultivar 93–11 reference genome with BSMAP [19].More than 80% of the clean reads were uniquely mapped with a depth of more than 10× and a conversion rate of at least 99% (Table S2).The repeatability of replicates in different stages was high,with Pearson correlation coefficient >0.96 between replicates (Fig.S1).As shown in Fig.1A,the blast fungus grew asymptomatically at dpi1 and dpi3 and showed overt lesions at dpi5.In contrast to the changes observed in Arabidopsis in response to pathogen attack,global DNA methylation levels were similar between untreated and infected rice,suggesting species-specific mechanisms in response to pathogen,bacterial,or oomycete infection.Only a slight increase in DNA methylation levels was observed in any of the three sequence contexts in infected rice (Fig.1B;Table S2).

In agreement with previous studies,DNA methylation was abundant in transposable element regions,which commonly occur at the pericentromeric regions of rice chromosomes (Fig.1C)[23,24].The TE information used in this study was re-annotated based on theO.sativassp.indicaASM465v1 assembly.A total of 447,331 TEs were annotated,more than the 400,176 TEs predicted by BGI ten years ago [25].Most TEs (49.2%) were located in intergenic regions,followed by 2 kb up-or downstream of the gene body region (41.2%) and the gene body region (9.5%) (Fig.S2).For all stages of pathogen infection,similar patterns and profiles of DNA methylation in all sequence contexts were observed across each chromosome.

3.2.Transcription levels of methylation-associated genes in response to pathogen infection

Although global DNA methylation patterns and levels were similar across the different stages after infection,there was a slight increase in DNA methylation in the three sequence contexts after infection.These changes prompted us to investigate whether genes involved in DNA methylation or demethylation showed altered expression levels in response to pathogen infection.We used RTqPCR to quantify the transcript levels of DNA methyltransferase genes includingMET1,CMT2,CMT3,DDM1,andDRM2,and demethylase genesROS1,DML2,andDML3.The transcript levels of most of the analyzed DNA methyltransferase genes showed less-than-twofold changes after infection compared to the uninfected control (Fig.1D).However,the DNA demethylase genesDML3aandROS1bshowed changes of greater than twofold in transcript abundance in response to infection,indicating that their expression was differentially regulated in response to the pathogen.Considering the significantly increased expression levels of these demethylation factors together with the largely similar DNA methylation levels across the different stages of infection,we hypothesized that both DNA methylation and demethylation pathway genes participate in the response of rice plants to pathogen infection.Moreover,in contrast to the DNA methylation changes induced by pathogens in Arabidopsis,there may be local rather than global DNA methylation changes in this response.

3.3.DNA methylation patterns in genic and intergenic regions after infection

DNA methylation is involved in silencing transposons and maintaining genome stability,gene methylation is always associated with gene expression regulation [26],and environmental stimuli such as abiotic and biotic stresses have been linked to transposon reactivation [27].To further investigate the methylation patterns of different genomic structures,we measured the DNA methylation profiles of both genic and TE regions.Similar to the patterns in untreated rice,CG-context methylation occurred exclusively in gene body regions,a pattern known to be conserved in many other plant species (Fig.2A).We observed striking differences of all three contexts of DNA methylation in flanking regions of genes,and these regions are important as binding sites for transcriptional regulatory elements (Fig.2A–C).Genic regions showed slightly increased methylation of all sequence contexts after pathogen infection,findings consistent with those observed for the global pattern of DNA methylation.

For TE regions,the main targets of DNA methylation,DNA methylation levels in all three sequence contexts were much higher than those in genic regions (Fig.2D–F).The discernible peaks in CHH methylation near both the start and end sites of TE body regions(Fig.2F)were similar to the pattern of CHH methylation enrichment described in the maize genome [28].Fig.2D,E shows that CG and CHG methylation were comparable in TE regions across different stages of infection.However,significant differences are seen in inset panels of Fig.2D,E after rescaling.Moreover,there were striking differences in CHH-context methylation across the stages of infection.The increased CHH methylation in infected rice suggests that the RdDM pathway is involved in pathogen response.

To further quantify the differences of DNA methylation in gene body and TE regions at the different infection stages,we examined the DNA methylation level of each gene and TE region.As shown in Fig.2G–I,there were significant increases in DNA methylation levels of gene body regions in all three sequence contexts after infection.For TE regions,the methylated levels of all methylation contexts were also significantly higher in infection stages than in the uninfected stage except for CG methylation at dpi3.We found that both global average DNA methylation analysis (using a metaplot)and collective gene/TE methylation examination(using a boxplot) were necessary to distinguish DNA methylation differences across different stages of infection.Collectively,our results showed that increased DNA methylation occurred in both genic and TE regions after infection,especially DNA methylation in flanking regions of genes and CHH methylation of TEs.Our subsequent differential methylation region (DMR) analysis supported this finding:most DMRs were enriched in TE regions and flanking regions of genes,and hypermethylated DMRs far outnumbered hypomethylated DMRs.

Fig.1.The rice leaf phenotype and DNA methylation changes after rice blast inoculation.(A)Rice leaf phenotypes after inoculation with M.oryzae.0,Control leaf inoculated with inoculation buffer;1,3,and 5 represent respectively one,three,and five days after M.oryzae inoculation.(B)Genome-wide average CG,CHG,and CHH DNA methylation levels after M.oryzae inoculation.(C)Circos plot of the O.sativa ssp.indica 93–11 chromosomes.Track order from outside to inside:TE density,gene density,CG methylation level for dpi0 to dpi5,CHG methylation level for dpi0 to dpi5,and CHH methylation level for dpi0 to dpi5.The heatmaps display average cytosine methylation in 100-kb windows.Green color represents low gene/TE density and low methylation levels.Red color represents high gene or TE density and high methylation levels.The scale indicates chromosome length in million bp (Mb) units.(D) qRT-PCR transcript fold changes of selected DNA methylation-related genes in infected rice compared to the untreated control.

3.4.Identification of differentially methylated regions (DMRs) in the rice genome after infection

Considering the subtle changes in the average DNA methylation level throughout the genome among the different samples,we focused instead on DMRs in the untreated control (dpi0) and the other three samples (dpi1,dpi3,and dpi5).Using the R package‘‘methylKit”,we obtained over 38,000 DMRs (Table S3,see details in Section 2).Most of the DMRs were of the CHH type,and the CG type was the least common.Among the different infection periods,the number of DMRs was greatest for dpi5,followed by dpi1 and dpi3.These findings suggest strong local DNA methylation changes on the first day after infection,followed by a reduction and then an increase to an even higher level on the fifth day.The faster reaction of DNA methylation in contrast to transcription may result in the higher number of DMRs in dpi1 than in dpi3.The subsequent gene expression changes in dpi3 may have contributed to the DMR maximum in dpi5.Although the number of CG DMRs was relatively stable across the different stages,the number of CHG and CHH DMRs fluctuated sharply with the infection process (Fig.3A).Except for dpi3/dpi0 (dpi3 compared to dpi0)CG and CHG DMRs,hypermethylated DMRs far outnumbered hypomethylated DMRs at all stages.This finding is in clear contrast with the hypomethylation induced by pathogen infection in Arabidopsis or phosphate starvation stress in rice,suggesting that the pathway induced by rice blast fungus is distinct from the biological pathway induced by phosphate starvation.

Fig.2.DNA methylation of gene and transposable element regions in rice after infection.(A–C)DNA methylation patterns of genes.(D–F)DNA methylation patterns of TEs.In the figure panels,‘‘-2 kb”denotes the 2 kb region upstream of the transcription start site(TSS)or TE body;‘‘2kb”denotes the 2 kb downstream of the transcription end site(TES)or TE body.The upstream region,gene/TE body and downstream region were divided into 20 proportional bins to draw the metaplot.(G–I) Box plots of CG,CHG,and CHH methylation of gene and TE regions.Asterisks indicate significant differences between dpi1/3/5 and dpi0 based on the Mann–Whitney test (**, P <0.01).

To determine whether DMRs occurred in specific genomic locations after pathogen infection,we examined the overlap between different genomic structures and DMRs.For CG,CHG,and CHH DMRs,the percentage of DMRs in gene body regions was very low.Most CG and CHH DMRs were located in the 2-kb upstream and downstream regions of genes,and most CHG DMRs occurred in intergenic regions.The genomic distribution of DMRs differed slightly among infection stages.For instance,the number of CG and CHG DMRs within the 2-kb upstream and downstream regions of genes decreased from dpi1 to dpi3 and then increased from dpi3 to dpi5.CHH DMRs showed an entirely different pattern than that observed for CHG DMRs (Fig.3B).Considering the importance of TEs in plant physiology and pathology,we also analyzed the overlap of DMRs and TEs using bedtools[29].Although only 38%of CG DMRs were located in TEs,over 60%of CHG DMRs and 80%of CHH DMRs overlapped with TEs.The The majority of CHH DMRassociated TEs were DNA-type (85.5%),followed by LTR (5.65%),SINE(2.51%),RC/Helitron(1.63%)and other TEs(Fig.3C).Hypermethylated CHH DMRs were significantly enriched in the MITE,SINE/tRNA and DNA/PIF-Harbinger families,whereas hypomethylated CHH DMRs were found in the MITE,DNA/TcMar-Stowaway and DNA/PIF-Harbinger families (Fig.3D).CHG DMR-associated TEs showed a different pattern:the most abundant family was LTR(49.21%),followed by DNA (38.39%),LINE (2.76%),RC/Helitron(1.58%)and other TEs(Fig.S3A).Hypermethylated CHG DMRs were markedly enriched in the LTR/Gypsy,Copia and rRNA families,whereas hypomethylated CHG DMRs were enriched in the LTR/Gypsy and DNA/CMC-EnSpm families (Fig.S3B).

3.5.Global transcriptional changes in response to pathogen infection

To investigate the global transcriptomic response to pathogen infection,we performed mRNA-seq simultaneously with BS-seq using the same plant samples.To evaluate the reproducibility of the RNA-seq data,two biological replicates were compared for each sample,with Pearson correlation coefficients between replicates ranging from 0.94 to 0.97 (Fig.S4).More than 80% of clean reads were uniquely mapped to the reference genome (Table S4).A total of 8830 DEGs were identified in at least one of the infection stages compared to the untreated control,and K-means analysis clustered these DEGs into four groups (G1,G2,G3,and G4)(Fig.4A,B).Most of the 8830 DEGs were differentially expressed at dpi5 (7036/8830) compared to dpi0,followed by dpi3(5808/8830) and then dpi1 (100/8830) (Table S5).The numbers of DEGs in the dpi3 and dpi5 datasets were about 58 and 70 times greater,respectively,than that in dpi1 (Fig.S5A),suggesting greater genome-wide transcriptional changes at the later infection stages.Across the different stages,the number of up-regulated DEGs was similar to the number of down-regulated DEGs.Most of the DEGs(4568)were identified in both the dpi3 and dpi5 datasets,with the expression change in the same direction(Fig.S5B,C).This finding indicates similar transcriptional regulation between these two stages.Comparing these findings to those of previous studies [30–32],we found many more new DEGs in our study and established a relationship between DEGs and DMRs (Fig.S6).

Fig.3.Genomic distribution of differentially methylated regions.(A) Hyper-and hypomethylated differentially methylated region (DMR) numbers for CG,CHG,and CHH after rice blast infection.(B)DMR genomic distribution.Gene body regions comprise the CDS and introns.Up-and downstream regions correspond respectively to the 2 kb upstream of the transcription start site(TSS)and the 2 kb downstream of the transcription end site(TES).The rest of the genome makes up the intergenic regions.(C)DMR distribution in TE and non-TE regions.The pie chart on the right is the distribution of TE-region CHH DMRs among the indicated TE classes.(D)Detailed values for hyper-and hypomethylated CHH DMRs for different TE classes.Significant differences between observed and expected numbers of DMRs are based on Fisher’s exact test (**,false discovery rate <0.01).

Fig.4.Hierarchical clustering and GO enrichment analysis of DEGs.(A)Heatmap of DEG expression levels across the different infection stages.The TMM value of each DEG was normalized by z-score,with the color gradient corresponding to the expression level.Black and white color represent respectively higher and lower expression,with shades of red between them.(B)Line plots of the G1,G2,G3,and G4 DEG groups across the four time points.The z-scores calculated in(A)were used to draw the lines.Each black line represents one DEG,and the red line represents the mean value of the group.N indicates the number of DEGs in each group.(C)Heatmap of enriched GO terms for G1,G2,G3,and G4 DEGs.Yellow,orange/red,and black indicate respectively no significant enrichment,significant enrichment,and extremely significant enrichment.

Numerous pathogen response genes were induced afterM.oryzaeinfection,including chitin binding protein AT2G43570(BGIOSGA01480).Genes belonging to several transcription factor(TF)families also showed significantly altered levels after pathogen infection (Table S6).Considering the well-established roles of TFs in regulating gene expression in many biological processes,we further analyzed the TF DEGs identified from the RNA-seq data.These 423 TF genes were divided into two groups by hierarchical clustering (Fig.S7A),and the expression patterns of these genes suggested that transcriptional regulation was dramatically altered during the pathogen infection process.The group I TF genes were highly expressed at dpi0 and dpi1 compared to the levels at dpi3 and dpi5,indicating these TF genes were active at the early stage of infection but showed decreased expression in the following days.Group II TF genes were induced at the late stage of pathogen infection (dpi3 and dpi5),suggesting that these TF genes are not immediately sensitive to pathogen infection and may be involved in a later response process.In some cases,TF genes belonging to the same family showed distinct expression patterns across infection stages.For example,some members of the NAC family,which is associated with biotic and abiotic stress responses,showed expression patterns similar to those of group I genes,while others showed expression patterns similar to those of group II genes(Fig.S7B).The same phenomenon was also observed for ERF and WRKY family genes (Fig.S7C,D).Overall,the varied expression patterns of these TF genes in response to pathogen infection suggested that plants use many different regulators to achieve a comprehensive response to pathogens.

We also performed GO analysis to identify enriched terms among the genes in the different DEG groups(Table S7).In all four groups (G1,G2,G3,and G4),genes related to cell communication,response to external stimulus and anatomical structure development were enriched,suggesting that genes with related functions show entirely different expression patterns during infection(Fig.4C).The analysis also suggested that genes associated with response to biotic stimulus,cell death,signal transduction,regulation of cellular processes,and stress response were downregulated after infection (G1 and G2),indicating that these pathways were suppressed byM.oryzae.In contrast,in the upregulated groups (G3 and G4),there was an enrichment of genes associated with translation,gene expression,cellular component organization,macromolecule modification and many biosynthetic and metabolic processes.These diverse functions hint at the dramatic changes occurring in plants during the late stages of infection.Genes involved in the developmental processes of reproduction,photosynthesis,generation of precursor metabolites and energy,embryo development,and lipid metabolic processes were enriched in the G3 and G4 groups.These functions are suggestive of accelerated energy compensation and reproduction processes at the late infection stages.The terms nitrogen compound metabolic process and response to stress were enriched among the G2,G3,and G4 genes.These included numerous NBS-LRR (R)genes,which are known to contribute to stress resistance in plants.For example,the expression levels ofBGIOSGA034857(disease resistance geneRGA2),BGIOSGA027329(disease resistance geneRPM1) andBGIOSGA033823(disease resistance geneRGA5)increased more than 20-fold after infection compared to untreated plants (Table S5).These DEGs may thus be involved in transcription regulation and stress response during pathogen infection.

3.6.Transcriptional changes associated with differential methylation

Pathogen infection induced global transcriptional changes in a temporal manner.To investigate the relationship between differential methylation and transcription in response to pathogen infection,we assessed the correspondence between DEGs and DMRs.We identified DMR-associated genes by assigning DMRs to nearby genes (flanking 2 kb region or gene body region) (Table S8).The numbers of DEGs associated with DMRs were as follows:2703 DEGs associated with 3428 CHH DMRs,163 DEGs associated with 170 CHG DMRs,and 139 DEGs associated with 140 CG DMRs(Fig.5A).Approximately 75%(3381)of the DMRs among the three infection stages were hypermethylated,whereas only 25% (1176)of the DMRs were hypomethylated.Most of the CHH DMRs associated with DEGs were located in TE regions,followed by the 2-kb upstream region,the 2-kb downstream region,introns,and the CDS of the DEGs.MITEs represented more than half of the total TEs containing CHH DMRs.Compared to the distribution of hypermethylated CHH DMRs,there was a slightly increased distribution of hypomethylated CHH DMRs in the upstream 2 kb region,the downstream 2 kb region,and introns(Fig.5B).We also performed GO enrichment analysis of these DMR-associated DEGs (DMRDEGs) (Table S9).Many enriched GO terms were found for CHH DMR-DEGs in dpi3/dpi0 and dpi5/dpi0,including cell communication,response to external and abiotic stimulus,anatomical structure development,and cellular component organization,suggesting that these processes are modulated by DNA methylation at late stages of infection (Fig.5C).In dpi3,genes involved in many synthetic and metabolic processes,translation,signal transduction,and photosynthesis processes were associated with CHH DMR variation.In dpi5,the analysis suggested that genes associated with response to stress and biotic stimulus and multiorganism processes are also under regulation by CHH methylation.

Given the substantial overlap of CHH DMRs with TEs(Figs.3C,D and 5B),we wondered whether the CHH DMRs broadly upregulated or down-regulated TE expression levels.By analyzing the RNA-seq data with TE annotation predicted by RepeatMasker,we found 537 and 522 DMR-related differentially expressed TEs(DETEs) in dpi3/dpi0 and dpi5/dpi0 (Table S10).Only 11 DETEs were associated with DMR-DEGs,indicating that the overwhelming majority of CHH DMR-DEGs (2692/2703) were unrelated to DETEs.However,the method we used to detect DETEs provided 2977 candidate TEs with proper expression levels,while other TEs (44,354/447,331,99.3%) were lowly expressed and can’t meet the screening criteria.Thus,whether pathogen infection influences the DNA methylation of TEs to activate them awaits further study.One representative gene isBGIOSGA016941,which encodes a heat shock factor (HSF) family member that specifically binds heat shock promoter elements (HSE).Both the expression level of the HSF protein gene and the upstream CHH methylation level increased along with the infection stages (Fig.5D).The upstream CHH DMR was located in a DNA/TcMar-Stowaway TE,but the increased methylation did not activate the expression of this TE(Table S10).One possibility is that the hypermethylated DNA/Tar-Stowaway TE effectively stabilized the upstream sequence,allowing transcription of the HSF protein gene to proceed more smoothly.Interestingly,the HSF family DMR-DEGs identified in the analysis were all up-regulated whether they were hypermethylated or hypomethylated in an upstream,downstream or intron region.In other words,the activation of HSFs not only was caused by changes of DNA methylation,but also involved other factors.Alternatively,the activation of HSF was independent of DNA methylation changes.Along withBGIOSGA016941,133 DMR-DEGs were annotated as TF genes belonging to 35 families (Table S8),including WRKY,MYB,ERF,and other TF families closely associated with stress resistance.Most MYB DMR-DEGs were downregulated,suggesting that plant resistance was influenced by DNA methylation that suppressed the expression of these MYB genes.

To identify the possible role of NBS-LRR genes during the infection process,we retrieved 494 NBS-LRR gene sequences ofO.sativassp.japonicacv.Nipponbare from Uniprot and identified 405O.sativassp.indicaputative NBS-LRR homologous genes using a similarity search (Table S11).Among those genes,83 were differentially expressed after pathogen infection,71% (59) of which were up-regulated.Of the 83 NBS-LRR DEGs,20 were associated with DMRs.As shown in Fig.S8,three representative NBS-LRR genes(BGIOSGA014719,BGIOSGA019764,andBGIOSGA030493) were upregulated in dpi5 in contrast to dpi0.We found at least one hypermethylated CHH-DMR located in TE at up-stream 2 kb or intron regions of them.Referring to the study of Deng and colleagues[33],we also found thePigm-Rhomologous gene inO.sativassp.indica(BGIOSGA022714) as DEG in our study.The expression level ofPigm-Rwas up-regulated in dpi3 and dpi5 samples,a finding consistent with that ofPigm-Rinduced resistance in rice [33].Many blast resistance genes have been identified in previous studies,and we retrieved the sequences of 28 major rice blastresistance genes from the China Rice Data Center (http://www.ricedata.cn/gene/gene_pi.htm,2012-6-20) to query them against theO.sativassp.indicagene set using BLAST.Seven of these genes were differentially expressed in dpi3 or dpi5,and three were DMRDEGs (Table S12).These findings indicate that pathogen infection induced many NBS-LRR genes and a handful of known rice blastresistance genes and suggest that DNA methylation participates in the expression regulation of these genes.

4.Discussion

Fig.5.Characterization of DMR-associated DEGs.(A)The number of differentially methylated region(DMR)-associated DEGs(DMR-DEGs).The bars at left and right indicate the number of DMRs and DEGs,respectively.(B) The distribution of hyper-and hypomethylated CHH DMRs in different regions of DMR-DEGs.MITEs and other TEs are excluded from upstream 2 kb,downstream 2 kb,CDS and intron regions of DMR-DEGs.(C)Enriched GO terms for CHH DMR-DEGs in dpi3 and dpi5.Yellow,orange/red and black indicate respectively no significant enrichment,significant enrichment,and extremely significant enrichment.(D) Transcription and methylation levels of HSF transcription factor gene BGIOSGA016941 and the upstream region on chromosome 4 (27,419,000–27,421,000 bp).The red box indicates DMRs.

DNA methylation is known to play crucial roles in many important biological processes,but few studies have focused on its role in plants under pathogen infection,especially in crop species and non-model plants.In rice infected withM.oryzae,many stressresponsive genes and disease-resistance genes were modulated by DNA methylation during pathogen infection.The DMRassociated DEGs included many TF genes,such as those belonging to the WRKY,MYB,NAC,ERF and bZIP families,which are associated with abiotic and biotic stress responses.In a previous study,the activation of PAMP-triggered immunity (PTI) was found to induce mitogen-activated protein (MAP) kinase signaling,transcriptional recoding mediated by WRKY TFs and the production of reactive oxygen species (ROS) [34].These defense reactions effectively prevent pathogen infection.Antagonizing members of the plant-specific WRKY TF family mediate transcriptional reprogramming,which is influenced by PTI.WRKY TF domains bind to the promoter W-box of defense genes to regulate their transcription.As a result,WRKY TFs act as regulators that help provide resistance to biotic stimuli [35].Our analysis identified three DMR-DEGs annotated as WRKY TF genes 1,47,and 53.Transgenic wild tomato (Solanum arcanumPeralta) overexpressingWRKY1showed higher levels of early blight resistance,whereas RNAisilenced lines showed decreased resistance toAlternaria solani[36].Thus,the observed up-regulation ofWRKY1may have helped rice respond toM.oryzae.Together with WRKY46 and WRKY70,WRKY53 promotes resistance toPseudomonas syringaein Arabidopsis,possibly enhancing the expression of salicylic acid (SA)-dependent genes to suppress jasmonic acid (MeJA)-induced expression of PDF1.2 [37].Based on this activity,downregulation ofWRKY53,as observed in the present study,may decrease plant resistance to bacterial pathogens but increase resistance to fungal pathogens.In addition to the WRKY family,other DMR-DEGs encoding different MYB,NAC,and ERF TFs also showed expression changes in different directions (i.e.,up-or down-regulated).All these diverse changes indicate that the rice blast response is a complex process involving epigenetic modifications and consequent diversified expression of various TF genes.

Rice blast infection induced slight but statistically significant hypermethylation in gene regions and the upstream and downstream regions of TEs.Hypermethylated DMRs far outnumbered hypomethylated DMRs,with varying degrees of hypermethylation both throughout the genome and at specific sites.While the observed hypermethylation could be a response of rice plants to rice blast infection,it may also increase the susceptibility of the plants by suppressing the expression of major rice blastresistance genes.One example is the rice blast-resistance factor RGA4 (BGIOSGA034263).Working with RGA5,RGA4 activates plant effector-triggered immunity (ETI) by directly interacting with theM.oryzaeeffectors AVR-Pia and AVR1-CO39 [38].In our study,down-regulation of anRGA4gene (BGIOSGA034263) was detected along with CHH hypermethylation in its downstream region.Likewise,a gene encoding Pi-d2 (BGIOSGA021200),annotated as a receptor-like serine/threonine protein kinase in the Uniprot database,was also down-regulated with CHH hypermethylation in its downstream region.Mutation of Pi-d2 is known [39] to increase rice susceptibility toM.oryzae.Similarly,hypermethylation has been associated with reduced resistance to pathogens in Arabidopsis mutants[40,41].Our genome-wide analysis provided evidence of hypermethylation,in some cases affecting stress-response genes,in rice infected with rice blast.

Our findings also suggest a link between hypermethylation induced by rice blast infection and TEs.The vast majority of CHG and CHH DMRs were located in TEs,and many TEs were associated with DEGs.It has been shown [42] that TE-associated genes have lower transcriptional activities than non TE-related genes.DNA methylation levels were also lower in transcribed TE-related genes[42].In rice,DNA methylation may influence genome stability and regeneration via TEs [43,44].Thus,hypermethylated TEs may induce rice susceptibility to rice blast fungus by affecting the expression of defense genes.As most DMRs were located in MITEs and the majority of CHH DMR-DEGs were also associated with MITEs,this TE class may provide clues to the role of DNA methylation in the susceptibility or defense of rice plants against rice blast fungus.Previous studies in rice have established that MITEs play roles in species diversity and gene expression.For example,genes near MITEs show lower expression levels than genes far from MITEs [45].In rice,a group of MITEs,mPing,regulates cold-and salt stress responsiveness by affecting nearby genes,possibly creating new stress-inducible regulatory networks [46].MITEs are also often found around R genes,suggesting their role in R gene evolution [47–49].Based on these reported functions,rice blast infection may lead to altered R gene expression by affecting the methylation of nearby MITEs to overcome the defense system of rice.As described above for the rice blast resistance factor RGA4,the decreased expression of oneRGA4gene may have been due to the CHH hypermethylation in a downstream MITE.Another major rice blast resistance genePitwas upregulated significantly,possibly as a result of the hypermethylation in an upstream MITE.One possibility is that the site of a hypermethylated MITE also influences the expression of a nearby gene.

During rice blast infection,epigenetic activation preceded expression variation of genes involved in the defense responses.Considering CHH hypermethylation of both genes and TEs,the activation of DNA methylation-related enzymes and the number of DMRs were both relatively high at dpi1,in contrast to the few DEGs at this stage.At dpi3,the number of DEGs surged,and the number of DMRs was slightly reduced.Thereafter,the number of DMRs and DEGs increased at a slower rate.These patterns may be indicative of the infection process,with DNA methylation playing a key part from beginning to end.We hypothesize that DNA methylation is used by rice blast to reduce plant resistance at the very early stage,with stable and maintained DNA hypermethylation at the later stages facilitating rice blast infection.However,how the rice blast pathogen could induce DNA hypermethylation in rice awaits future studies.In view of the large roles of other epigenetic modifications in plant response to biotic stress [50],investigating the roles of multiple epigenetic modifications,such as siRNA,miRNA,and DNA methylation,in rice response to blast fungus infection will provide epigenetic clues to breeding highly resistant rice cultivars.

CRediT authorship contribution statement

Haifeng Wang and Zonghua Wang:conceived and designed this research.Na Cui and Yan Shi:conducted methylation and transcriptome analysis.Xiaofeng Chen:performed the experiments and collected tissues.Meirong Chi,Jiantuan Hu,and Kunlong Lai:provided comments and discussion.Na Cui and Haifeng Wang:wrote the manuscript.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the Natural Science Foundation of Fujian Province(2018J06006),National Key Research and Development Program of China (2016YFD0300700) and National Natural Science Foundation of China (31770156),State Key Laboratory of Ecological Pest Control for Fujian and Taiwan Crops(SKL2018006),the Program for New Century Excellent Talents of Fujian Province University,and the Pre-eminent Youth Fund and Distinguished Young Scholars of Fujian Province.This study was supported by the Supercomputing Center at the College of Plant Protection of Fujian Agriculture and Forestry University.

Data availability

The data reported in the present paper were deposited into the Short Read Archive (SRA) database of NCBI (accession number SRP230674 and SRP230704).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2020.10.002.