Next Article in Journal
Process-Based Modeling of Phenology and Radial Growth in Pinus tabuliformis in Response to Climate Factors over a Cold and Semi-Arid Region
Next Article in Special Issue
Dissection of Common Rust Resistance in Tropical Maize Multiparent Population through GWAS and Linkage Studies
Previous Article in Journal
Long-Lasting Silver Nanoparticles Synthesized with Tagetes erecta and Their Antibacterial Activity against Erwinia amylovora, a Serious Rosaceous Pathogen
Previous Article in Special Issue
Genome-Wide Association Analysis Uncovers Genes Associated with Resistance to Head Smut Pathotype 5 in Senegalese Sorghum Accessions
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Plants
Volume 13
Issue 7
10.3390/plants13070982
plants-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Transcriptome-Wide N6-Methyladenosine (m6A) Methylation Analyses in a Compatible Wheat–Puccinia striiformis f. sp. tritici Interaction

by
Elif Naz Cerav
†,
Nan Wu
† and
Mahinur S. Akkaya
*
School of Bioengineering, Dalian University of Technology, No. 2 Linggong Road, Dalian 116024, China
*
Author to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Plants 2024, 13(7), 982; https://doi.org/10.3390/plants13070982
Submission received: 25 February 2024 / Revised: 21 March 2024 / Accepted: 27 March 2024 / Published: 29 March 2024
(This article belongs to the Special Issue Molecular Biology and Genomics of Plant-Pathogen Interactions)
Browse Figures
Versions Notes

Abstract

:
N6-methyladenosine (m6A) is a prevalent internal modification in eukaryotic mRNA, tRNA, miRNA, and long non-coding RNA. It is also known for its role in plant responses to biotic and abiotic stresses. However, a comprehensive m6A transcriptome-wide map for Puccinia striiformis f. sp. tritici (Pst) infections in wheat (Triticum aestivum) is currently unavailable. Our study is the first to profile m6A modifications in wheat infected with a virulent Pst race. Analysis of RNA-seq and MeRIP-seq data revealed that the majority of differentially expressed genes are up-regulated and hyper-methylated. Some of these genes are enriched in the plant–pathogen interaction pathway. Notably, genes related to photosynthesis showed significant down-regulation and hypo-methylation, suggesting a potential mechanism facilitating successful Pst invasion by impairing photosynthetic function. The crucial genes, epitomizing the core molecular constituents that fortify plants against pathogenic assaults, were detected with varying expression and methylation levels, together with a newly identified methylation motif. Additionally, m6A regulator genes were also influenced by m6A modification, and their expression patterns varied at different time points of post-inoculation, with lower expression at early stages of infection. This study provides insights into the role of m6A modification regulation in wheat’s response to Pst infection, establishing a foundation for understanding the potential function of m6A RNA methylation in plant resistance or susceptibility to pathogens.

1. Introduction

Many biological processes are controlled by RNA molecules, which play a vital role in transferring genetic information. To perform specific molecular functions, RNA transcripts can be chemically modified. Based on the MODOMICS database (https://iimcb.genesilico.pl/modomics/ (accessed on 1 October 2023)), there have been more than 180 RNA chemical modification discoveries to date, mainly in transfer RNAs (tRNAs) and ribosomal RNAs (rRNAs), where they play important roles in RNA function [1]. The most common post-transcriptional modifications to RNA are the following: N6-methyladenosine (m6A), 5-methylcytosine (m5C), N1-methyladenosine (m1A), pseudouridine (Ψ), inosine (I), and N6,2′-O-dimethyladenosine (m6Am) [2].
In eukaryotic mRNA, m6A is the most abundant methylation type which is also found in bacteria and RNA viruses [3]. The m6A RNA modification is a preserved regulatory system across highly diverse organisms; the system encompasses the methyltransferase complexes known as “writers” that add methyl groups to mRNA, demethylases termed “erasers” that remove these modifications for dynamic regulation, and m6A-binding proteins or “readers” which recognize and interact with methylated RNA sites, collectively playing a crucial role in mRNA function by impacting processes such as stability, splicing, localization, and translation, thereby illustrating the evolutionary significance of m6A modifications [4,5]. There are several m6A methyltransferases, called m6A writers, including mRNA adenosine methylases (MTA and MTB), VIRILIZER (VIR), FK506-binding protein 12 (FKBP12)-interacting protein 37 (FIB37), and the E3 ubiquitin ligase HAKAI [6]. In addition, m6A demethylases or erasers, such as the alkylated DNA repair proteins AlkB homologs (e.g., AlkBH4b, AlkBH8, AlkBH9B, AlkBH10B, AlkBH11B) [7,8] and fat mass and obesity-associated protein (FTO) [9], can dynamically modulate m6A modification. Moreover, there is a series of m6A-binding proteins for recognizing m6A-containing mRNAs named m6A readers like the EVOLUTIONARILY CONSERVED C-TERMINAL REGION (ECT) family (e.g., ECT2, ECT3, ECT4). It has been demonstrated that ECT2 is involved in mRNA stability [10]. YT521-B homology (YTH) family proteins, whose YTH domain can recognize and bind m6A-containing RNA, are the main readers of m6A [11]. The m6A modification very much plays a pivotal role in regulating plant transcription and gene expression, growth, and development [12].
The most canonical consensus motif of the m6A modification site is RRACH (R = G or A; H = U, A or C) [13], which accounts for more than 90% of the m6A sites in wheat mRNAs [14]. Numerous studies in plants also identified this dominant consensus motif [15,16,17]. Alternative motifs of GGAU and URUAY (Y = C or U) were also found to be enriched in plant species by epitranscriptome analysis [18,19]. The distribution of m6A modification within the mRNA is not uniform and is mainly enriched close to the stop codon and 3’ untranslated regions (3’ UTRs) [20]. Despite preliminary research on m6A methylation in the plant response to abiotic stress [21,22,23], biotic stress has received far less attention. To highlight a few studies, in Arabidopsis thaliana, AtALKBH9B demethylase was reported to erase m6A from RNAs of alfalfa mosaic virus (AMV), and its interaction with the CP (coat protein) of AMV resulted in successful viral infection. The suppression of AtALKBH9B led to the hyper-methylation of AMV RNAs and attenuated AMV infectivity, whereas the infection of cucumber mosaic virus (CMV) was unaffected, which might correspond to the absence of interaction between AtALKBH9B and the CMV CP [24]. The mRNAs of certain genes of pear (Pyrus bretschneideri) involved in defense mechanisms were m6A-modified to enhance the abundance of defense-related transcripts against the inoculation of fire blight pathogen Erwinia amylovora [25].
The plant defense regulation mechanism against pathogen assaults has been studied over the years. Many genes were identified based on the changes in expression levels, and the biological and biochemical roles were inferred. The involvement of NBS-LRR genes in plant–pathogen interactions, MAPK signaling, and the biosynthesis of secondary metabolites have been highlighted, and more detailed information is accumulating [26]. Similarly, the pathogenesis-related protein 1 (PR1), along with calcium-dependent protein kinases and calmodulin genes, play crucial roles in plant–pathogen interactions, in signal transduction. The biosynthesis pathways of phenylpropanoids and flavonoids have also been demonstrated [27]. In a study, Wang et al., 2016, emphasized the significance of mitogen-activated protein (MAP) kinase genes in the biosynthesis of phenylpropanoids, flavonoids, stilbenoids, and defense-related pathways [28]. 3-ketoacyl-CoA synthase genes are vital in transcription regulation, signal transduction, plant–pathogen interactions, and MAPK signaling pathways [29]. Additionally, WRKY transcription factor genes, as reported by Chen et al., 2014, are key players in carbohydrate metabolism, fatty acid metabolism, oxidative phosphorylation, and defense and signal transduction pathways [30]. It has also been reported that respiratory burst oxidase genes are involved in the regulation of transcription, transport of organic acids and metal cations, carbohydrate metabolism, and photosynthesis-antenna proteins [31]. Heat shock protein genes were shown to be critical in transcription regulation, metabolism, and signal transduction in stress and defense-related pathways [32]. In interactions between plants and pathogens, the significance of cysteine protease genes in the biosynthesis of phenylpropanoids and defense mechanisms was explored [33] through the extensive identification and characterization of the differential expression during biotic and abiotic stresses [34].
One of the most harmful diseases to affect wheat (Triticum aestivum) worldwide is known as stripe rust or yellow rust disease and is brought on by Puccinia striiformis f. sp. tritici (Pst). Pst leads to a substantial decrease in grain quality and in yield [35]. For long-term field resistance, one aspect is to thoroughly understand the molecular-level interaction mechanisms between wheat and Pst, elucidating how wheat responds to stripe rust infection at the molecular level. This study, from the perspective of m6A RNA methylation, combined RNA-seq (RNA sequencing) and MeRIP-seq (methylated RNA immunoprecipitation sequencing) analyses to understand the changes in gene expression levels and m6A modifications in wheat–Pst compatible interaction. Expression level changes in m6A regulators of wheat at various time points were also investigated, exploring how m6A-mediated post-transcriptional regulation affects the Pst-infected wheat, emphasizing a further complexity to the interactions in the pathosystem. We identified the genes going through m6A methylation transcriptional regulations in many of the above pathways, introducing an additional level of understanding to the changes in the expression levels.

2. Results

2.1. Transcriptome Profile of Wheat under Compatible Pst Interactions

To understand how the wheat transcriptome responds to attack by Puccinia striiformis f. sp. tritici (Pst) and how it is modified in such an interaction, this study used the Chinese yellow rust race, CYR32 (a virulent Pst race on many wheat cultivars in China [36]), to inoculate Avocet/Yr7, creating a compatible interaction model. The total RNA extracted from 7 dpi inoculated and CK samples were utilized for next-generation sequencing (NGS) as input samples, followed by an in-depth analysis of the sequencing data (Figure 1A). The phenotype of inoculated wheat leaves compared to CK showed slight chlorosis with the infection sites exhibiting necrotic, chlorotic flecks and blotches and, only in a few, sporadic sporulation (Figure 1B,C). The successful infection was verified with the expression, and the accumulation of the protein synthesis elongation factor 1 gene, PsEF1, belonging to Pst was verified with consistency across two biological replicates (Figure 1D and Figure S1A).
RNA-seq data from two biological replicates each from CK and infected samples, after filtering for adapters (Table S1), quality control (Table S2), and removal of rRNA (Table S3), were aligned to the wheat reference genome (IWGSC CS RefSeq v2.1, https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_018294505.1 (accessed on 10 September 2023)) (Table S4), producing approximately 37–61 million reads for each of the four samples. Principal component analysis (PCA) of the four samples, considering the distribution of sample points along PC1 and PC2, showed that I-1_In and I-2_In were closely related, as were CK-1_In and CK-2_In (Figure 2A). Pearson correlation analysis of the four samples, shown in a heatmap (Figure 2B), revealed that both the correlation coefficient between I-1_In and I-2_In and between CK-1_In and CK-2_In exceeded 0.9. Integrating the PsEF1 expression data from the two biological replicates (Figure 1D), the sequencing quality control, and genome alignment results (Tables S1–S4), along with the PCA and Pearson correlation analysis (Figure 2A,B), we demonstrated the reproducibility of the data, which are suitable for further analysis.

2.2. Most Differentially Expressed Genes (DEGs) in Wheat Were Induced by Pst Infection, while a Significant Decrease Was Observed in Genes Related to Photosynthesis

Expression levels were first presented as a raw read count, which indicates the number of reads corresponding to a transcript. However, the raw read count is influenced by sequencing depth and gene length, making it unsuitable for comparing differential gene expression across samples. To ensure accuracy in subsequent analyses, we normalized the sequencing depth, gene, or transcript length to obtain the TPM (transcript per kilobase per million mapped reads) for each gene in four samples before further analysis (Table S5). Without considering the FDR (false discovery rate), comparing TPM values of infected and CK samples revealed a fold change (FC), with log2(FC) > 0 indicating 69,083 genes with up-regulation in the infected sample and log2(FC) < 0 indicating 28,424 genes with down-regulation. After further filtering for genes with |log2(FC)| > 1 and FDR < 0.05, the differentially expressed genes (Table S6) were identified, 4100 up-regulated and 1814 down-regulated genes (Figure 2C), showing that most differentially expressed genes in wheat were induced by the pathogen.
Hierarchical clustering was performed on the expression patterns of 5914 differentially expressed genes (Table S7), categorizing them into 20 clusters, and the clustering results are displayed in the heatmap (Figure 2D). These genes, with similar expression patterns, may have common functions or be involved in the same metabolic pathways and signaling pathways. Among them, cluster No. 20, which is enriched with the highest number of differentially expressed genes, involves 1010 genes that are in the pathogen-infected samples. Within these 1010 genes, 13 genes are annotated in the KEGG pathway as a “plant-pathogen interaction”, including Calmodulin (TraesCS1B03G1255400, TraesCS5B03G0150400, TraesCS7B03G1109900, TraesCS7D03G1164400), NBS-LRR disease resistance protein (TraesCS1A03G0059400, TraesCS7D03G0011300), pathogenesis-related protein 1 (TraesCS5B03G1087600, TraesCS7D03G0450000), 3-ketoacyl-CoA synthase (TraesCS4A03G0012800), calcium-dependent protein kinase (TraesCS4A03G0728000), cyclic nucleotide-gated channel (TraesCS5D03G0900900), protein kinase family protein (TraesCS6B03G0271200), and chaperone protein htpG family protein (TraesCS7A03G1288100), indicating genes related to calcium ion channel regulation in plant immunity, showing a degree of a coordinated response to combat Pst invasion.
The significantly expressed genes with up- and down-regulations in comparison to CK are presented in the volcano plot with FDR values. Many of the genes appear within the fold change values of 1 < |log2(FC)| < 5, having relatively low FDR values, with highly credible differential changes. The three most significantly up-regulated genes are TraesCS4B03G0119300, TraesCS5B03G1087600, and TraesCS5D03G0710500LC; the most significantly down-regulated genes are TraesCS1A03G0499300, TraesCS5D03G0016200LC, and TraesCS6D03G0396800 (Figure 2E). Among the up-regulated genes, TraesCS5B03G1087600, a pathogenesis-related protein 1 (PR1), is considered an important defense protein. The production and accumulation of plant PR proteins are key responses to various biotic and abiotic stresses [37]. The two most significantly down-regulated genes are related to photosynthesis. In summary, wheat RNA expression was mostly up-regulated by Pst.
Differential genes were mapped to terms in the Gene Ontology (GO) database (http://www.geneontology.org/ (accessed on 11 October 2022)) as a molecular function, cellular component, and biological process. This mapping allowed for the determination of the number of differentially expressed genes per term, yielding lists of differential genes associated with specific GO functions. A hypergeometric test was then applied to identify GO terms significantly enriched among the differentially expressed genes in comparison to the background, resulting in 3204 terms enriched in biological process, 1611 in molecular function, and 540 in cellular component (Tables S8–S10). The top twenty enriched GO terms are presented in a GO enrichment circle plot (Figure 3A), with fifteen items enriched in biological process and five in molecular function. Each GO term predominantly featured more up-regulated than down-regulated genes, consistent with the overall changes observed in differentially expressed genes. Specifically, the biological process GO terms GO:0016054 (organic acid catabolic process), GO:0046395 (carboxylic acid catabolic process), and GO:0044282 (small molecule catabolic process) and the molecular function GO terms GO:0016853 (isomerase activity) and GO:0016829 (lyase activity) showed a lower q-value and a higher rich factor, contributing to high significance. Due to fewer entries in the cellular component, this GO term was not noticeable in the enrichment circle plot. Therefore, the top 20 enriched items under cellular component were further highlighted (Figure 3B), showing clear enrichment in terms such as “plastid”, “chloroplast”, and “thylakoid”, indicating a relation to photosynthesis. To compare the ratio of up-regulated to down-regulated genes within the main GO terms of biological process, molecular function, and cellular component (Table S11), a GO enrichment bar chart was created (Figure 3C). Integrating RNA-Seq data analysis, 24 annotations for 35 genes were identified in the GO term related to the “immune system process”, elucidating the interplay of genes involved in plant immune/defense mechanisms against pathogens (Table 1. In most GO terms, the number of up-regulated differential genes exceeded that of down-regulated ones, except for a few terms where down-regulated genes were more numerous: GO:0032991 (macromolecular complex), GO:0044422 (organelle part), GO:0099512 (supramolecular fiber), GO:0009295 (nucleoid), GO:0031012 (extracellular matrix), GO:0044420 (extracellular matrix component), GO:0005198 (structural molecule activity), GO:0045735 (nutrient reservoir activity), and GO:0045182 (translation regulator activity).
Pathway analysis based on KEGG (Kyoto Encyclopedia of Genes and Genomes), a major public database for pathways, aids in further understanding the biological functions of genes. Pathway significance enrichment analysis uses KEGG pathways as units and applies a hypergeometric test to identify pathways significantly enriched among the differential genes in comparison to the entire background. This analysis can determine the primary biochemical metabolic and signal transduction pathways in which the genes are involved (Table S12), enriching a total of 121 KEGG A class pathways, with 97 annotated as metabolism. The list of the “plant-pathogen interaction pathway” through RNA-Seq data and KEGG revealed 18 annotations for 80 genes, highlighting the dynamic battle where increased pathogen virulence exerts selective pressure on plants to enhance or alter their immune/defense mechanisms (Table 2). The top 20 enriched KEGG pathways are presented in a KEGG enrichment circle plot (Figure 3D), and those with the smallest q-values are depicted (Figure 3E). Notably, the pathway KO00196 (photosynthesis-antenna proteins) had the highest rich factor, and all the enriched genes were significantly down-regulated. In summary, following infection by Pst, wheat’s transcriptome was predominantly up-regulated, with fewer pathways of down-regulated genes. In particular, genes related to photosynthesis were notably and comprehensively down-regulated.

2.3. The Landscape of m6A Methylome in Pst-Infected Wheat Leaves

Methylated RNA immunoprecipitation (MeRIP) is based on the principle of recognizing and binding m6A-modified RNA with m6A-specific antibodies, enriching methylated fragments through RNA immunoprecipitation, and combining this with high-throughput sequencing to study regions of RNA undergoing m6A methylation modification across the entire transcriptome (Figure 1A). Partial RNA from four samples for RNA-seq, CK-1, CK-2, I-1, and I-2, underwent enrichment and elution with m6A-specific antibodies for MeRIP-seq analysis. After filtering low-quality data and adapters (Table S13), quality control (Table S14), and removal of rRNA (Table S15), the data were aligned to the wheat reference genome (Table S16), resulting in approximately 41–54 million reads for each of the four samples. Peak analysis identified regions and sites of RNA with m6A modification across the whole genome. Using the R package exomePeak2 [38], peak scanning (peak calling) was conducted on a group basis across the genome with a default threshold of p-value < 0.05. This was followed by the analysis of the genomic location of peaks, peak region sequence information, and annotation of genes (Tables S17 and S18), leading to the screening and annotation of peak-associated genes (Tables S19 and S20).
Principal component analysis (PCA) was performed using uniquely mapped reads, specifically by dividing the genome into 10 kb windows, counting the number of reads in each window, and reducing the read counts across all windows to a few principal components. PCA of the four MeRIP-seq samples, considering the distribution of sample points along PC1 and PC2, showed that I-1_In and I-2_In were closely related, as were CK-1_In and CK-2_In (Figure 4A). A heatmap was generated with the pairwise sample correlation using uniquely mapped reads by dividing the genome into 10 kb windows, counting the number of reads in each window, and then calculating the Pearson correlation coefficient for read counts across all windows. The Pearson correlation coefficients for both I-1_In vs. I-2_In and CK-1_In vs. CK-2_In were greater than 0.9. A higher correlation coefficient indicated a higher similarity in RNA methylation patterns between samples (Figure 4B). Integrating the PsEF1 expression data from the two biological replicates (Figure 1D), the sequencing quality control, genome alignment, and peak calling results (Tables S13–S20), along with the PCA and Pearson correlation analysis (Figure 4A,B), demonstrated reproducibility among biological replicates, high sequencing quality, and alignment accuracy, making the MeRIP-seq sequencing data suitable for further analysis.
In CK and infected samples, the number of peaks was 15,249 and 17,864, respectively (Table S21); the infected ones not only showed a higher number of peaks in comparison to CK but also a greater chromosomal peak coverage (Figure S2). Further analysis of the distribution of peaks in the 5′ UTR (untranslated region), start codon, CDS (coding sequence), stop codon, and 3′ UTR (Figure 4C and Figure S3) showed that the number of peaks in infected samples was higher across different gene functional elements in comparison to CK, with the largest increases observed at the 5′ UTR and start codon, which were 40.34% and 30.37%, respectively. The smallest increase was in the 3′ UTR at 11.78%, but the changes in the 3′ UTR peaks were concentrated in the middle position of the 3′ UTR. Meanwhile, the peak distribution in this study shows a similar m6A modification distribution to other plant species, regardless of whether they are monocots or dicots, whether in fruits or other tissues or under some stress conditions [39].
The m6A modification is the most prevalent internal modification in eukaryotic mRNA, with m6A methylation being highly conserved, typically embedded in conserved consensus sequences like 5’-DRACH-3’ and 5’-RRACH-3’ (D = G/A/U, R = G/A, H = A/U/C). The frequency distribution of specific motifs (RRACH, DRACH) within all peaks across all samples was analyzed using HOMER to construct average base frequency matrices. In CK sample two, in the infected sample, three significant motifs were identified, two of which were very similar to the two motifs found in the control sample. A particular motif emerged as a unique methylation site motif, GTACAA, which might specifically be determined by Pst (Figure 4D). We found that 6174 genes with this motif were methylated in our experiments. These genes associated with disease resistance include RPP13-like, RGA2-like EDR2L, putative RGA1, and RPM1-like. Similarly, genes related to the immune response, including those involved in the immune system and effector processes, such as WD repeat-containing protein 7, RRP44A, MACPF domain-containing protein, and numerous others like Histone acetyltransferase and IAN9 (immune-associated nucleotide-binding protein 9), are essential for the robust functioning of the immune system. These genes were identified in a sample infected with Pst having the GTACAA motif, which is a novel methylation site linked to this pathogen. The discovery of this unique methylation motif indicates a specific epitranscriptional response to Pst infection, suggesting a sophisticated mechanism of plant–pathogen interaction and defense regulation (Figure 4D, Table S22).

2.4. The Majority of Wheat Genes Associated with Differentially Methylated Peaks (DMPs) Showed Hyper-Methylation in Response to Pst

RNA-seq data from “CK” and “Infected” samples (Table S5) and peak data (Tables S19 and S20) were merged, resulting in 22,629 merged peaks and their peak-related gene expression counts and annotations across the four sequencing samples (Table S23). TPM and p-values were further calculated. The filtering criteria of |log2(FC)| > 1 and p-value < 0.05 identified differentially expressed genes associated with differentially methylated peaks (DMP-related genes) (Table S24), totaling 1804 DMP-related genes. Of these, 1167 were up-regulated DMP-related genes (indicating hyper-methylation) and 518 were down-regulated DMP-related genes (indicating hypo-methylation) (Figure 4E), suggesting that most DMP-related genes exhibit a trend towards hyper-methylation upon Pst infection.
Volcano plot analysis was conducted for DMP-related genes comparing infected samples to CK samples (Figure 4F). The volcano plot displays the situation of differentially expressed DMP-related genes, with genes closer to edges showing greater differences. DMP-related genes around |log2(FC)| = 10 had relatively lower FDR values, indicating more credible and significant changes. The three most significantly hyper-regulated genes were TraesCS3B03G1139000, TraesCS2A03G0079700, and TraesCS7D03G0606900; the most significantly hypo-regulated genes were TraesCS1A03G0058100, TraesCS7D03G0910600, and TraesCS1A03G0605600 (Figure 4F). Among the significantly hypo-regulated genes, TraesCS1A03G0058100 is annotated as a putative NBS-LRR disease resistance protein. To our knowledge, this may be the first report of an NBS-LRR gene showing hypo-methylation after Pst–wheat compatible interaction, suggesting that the m6A mechanism could directly regulate key genes, such as resistance genes, in plant immunity.
DMP-related genes were enriched in the GO categories of biological process, cellular component, and molecular function, with respective counts of 2838, 571, and 1274 (Tables S25–S27). The top 20 enriched GO terms, all from the biological process category (Table S25), were presented in a GO enrichment circle plot (Figure 5A), with a relatively high rich factor for GO:0006767 (water-soluble vitamin metabolic process), GO:1901607 (alpha-amino acid biosynthetic process), and GO:0008652 (cellular amino acid biosynthetic process) and a relatively low q-value for GO:1901564 (organonitrogen compound metabolic process), GO:0044711 (single-organism biosynthetic process), and GO:0044283 (small molecule biosynthetic process). The top 20 entries under the cellular component category (Table S26) showed enrichment as in the RNA-seq cellular component results, with DMP-related genes enriched in GO cellular component annotations related to photosynthesis, such as “chloroplast”, “plastid”, “thylakoid”, and “stroma” (Figure 5B). Further analysis of the quantity and annotations of DMP-related genes in the main GO terms of biological process, cellular component, and molecular function (Table S28) led to the creation of a GO enrichment bar chart (Figure 5C).
Significant enrichment in KEGG pathways identifies the primary biochemical metabolic and signal transduction pathways involved by DMP-related genes (Table S29), enriching a total of 115 KEGG A class pathways, with 84 annotated as metabolism. The top 20 enriched KEGG pathways are presented in the KEGG enrichment circle plot (Figure 5D); those with the smallest q-values are depicted in a bubble chart (Figure 5E). In the KEGG enrichment circle plot, sixteen entries were enriched in metabolism and four in genetic information processing. The bubble chart shows that aside from the first 14 pathways with q-values less than 0.05 indicating significance, the q-values of other enriched pathways were greater than 0.05, suggesting lower significance. The pathways with higher enrichment factors and smaller q-values include KO00860 (porphyrin metabolism), KO00920 (sulfur metabolism), KO00030 (pentose phosphate pathway), KO00410 (beta-alanine metabolism), KO00280 (valine, leucine, and isoleucine degradation), and KO00300 (lysine biosynthesis). This suggests that in wheat, upon compatible interaction with stripe rust, these pathways of secondary metabolism were regulated by the m6A mechanism, thereby affecting wheat’s growth and development.

2.5. Association Analysis of DEGs and DMP-Related Genes in Pst-Infected Leaves in Comparison to CK

To further clarify the relationship between gene expression levels and m6A modifications, an analysis was conducted based on the number of genes (Table S5) and the number of peaks (Table S23) detected in both CK and infected samples from combined RNA-seq and MeRIP-seq data, leading to Figure 6A. Comparing the number of m6A modification peaks between CK and infected samples, there were 11,685 common peaks, 3565 unique peaks in CK samples, and 6179 unique peaks in infected samples, showing an approximate 73% increase in peak types, indicating that m6A modifications increase due to infection (Figure 6A, left panel). Further analysis of the number of genes associated with these peaks revealed 11,727 common genes, 2928 unique genes in CK samples, and 5313 unique genes in infected samples (Figure 6A, right panel). Typically, one or more peaks map to a gene, so the number of genes associated with peaks is less than the number of peaks. In this study, the total number of peaks in CK and infected samples was 21,428, with 19,968 genes associated with peaks, aligning with the expected mapping ratio between peaks and genes (Figure 6A). However, the number of common genes is slightly higher than the number of common peaks, possibly because most peaks are not regulated upon Pst infection associated with genes having only one m6A modification peak, leading to a similar count of common genes and peaks. Also, a gene may have multiple transcripts, and if a peak maps to an overlapping part of a transcript, it could be associated with multiple transcripts, thus increasing the count of common genes over common peaks. Interestingly, the number of genes associated with peaks in both CK and infected samples is significantly lower than the number of peaks, suggesting that genes associated with unique peaks are more likely to be regulated by more than one m6A modification. Genes modified by multiple peaks are more likely to be involved in response to the Pst pathogen, and their numbers are substantial, indicating that the transcriptome is extensively regulated by the m6A mechanism by the pathogen.
To further correlate m6A modification with gene expression, a co-differential gene analysis was conducted. This involved combining differentially expressed genes from RNA-seq (Table S6) with differential m6A peaks from MeRIP-seq (Table S24) together (Table S30). By using a filter criterion that required |log2(FC)| > 1 in both RNA-seq and MeRIP-seq, an FDR < 0.05 in RNA-seq data, and a p-value < 0.05 in MeRIP-seq data, a total of 1242 significant co-differential genes were filtered (Table S31). A quadrant volcano plot allows for the presentation of the expression level and the methylation level changes together (Figure 6B). The largest group consisted of 642 genes with significant hyper-methylation and up-regulated gene expressions. Following this, 556 genes showed significant hypo-methylation with down-regulated expression. Additionally, 34 genes had significant hypo-methylation with up-regulated expression, and 10 genes had significant hyper-methylation with down-regulated expression.
To focus deeper on the annotation information and function of markedly co-differentiated genes, the filtering criteria were further tightened by adjusting |log2(FC)| > 1.5 in both RNA-seq and MeRIP-seq data, FDR < 0.05 in RNA-seq data, and a p-value < 0.05 in MeRIP-seq data. Under these new settings, 346 genes were identified in the hyper-up category, 212 in hypo-down, 10 in hypo-up, and 2 in hyper-down. For genes in the hyper-up and hypo-down categories, the ones with the lowest FDR and p-value were further selected. A total of 48 genes are presented in Table 3, along with gene IDs, annotations, expression patterns, peak distributions, peak start positions, IP_log2(FC) values, and Input_log2(FC) values. The expression levels of these 48 genes and the |log2(FC)| from MeRIP-seq are displayed as histograms, and |log2(TPM + 1)| from RNA-seq are shown as heatmaps in Figure 6C. Unlike the genes in the hyper-up category, the overall expression level of significantly up-regulated genes in the hypo-up category was lower across the transcriptome. In the hypo-down category, except for TraesCS1D03G0954900 and TraesCS2D03G0286600LC, the other 12 significant co-differential genes were annotated as chlorophyll a/b binding proteins. Overall, this study found that by stripe rust invasion, most genes in wheat responded in a hyper-up manner, consistent with the observation that most genes in RNA-seq were up-regulated and most genes in MeRIP-seq showed hyper-methylation. A notable discovery was that photosynthesis-related genes, possibly due to the secretion of effectors or other pathogenic factors from Pst, might be involved in the regulation of the wheat m6A mechanism, leading to hypo-methylation modifications, decreased expression of photosynthesis-related genes, and impaired photosynthesis function, thus affecting basic biomolecular metabolism and secondary metabolite production, further reducing the wheat resistance to stripe rust and resulting in successful invasion.

2.6. Divergent Expression Dynamics of m6A Regulators in Response to Pst Invasion

Since previous research rarely focused on the expression changes in m6A regulators in the wheat–Pst pathosystem, we further analyzed the expression patterns of m6A regulators in RNA-seq and MeRIP-seq data. According to a review by Yue et al., (2019) [40], which analyzed and annotated 90 wheat m6A regulator genes, including 20 m6A writers, 29 m6A erasers, and 41 m6A readers (Table S32), 28 of these m6A regulator genes were detectable in our RNA-seq and MeRIP-seq data, consisting of 12 m6A writers, 5 m6A erasers, and 11 m6A readers (Figure 7A). We further combined the average of TPM + 1 from two biological replicates of CK and infected samples in RNA-seq and MeRIP-seq and created a heatmap after converting to log2(TPM + 1) (Figure 7A). The heatmap shows that these 28 m6A regulator genes exhibit a variety of different expression patterns, whether they are m6A writers, m6A erasers, or m6A readers, or in terms of their expression patterns in RNA-seq and MeRIP-seq, indicating that the m6A regulatory mechanism in wheat–Pst interactions is complex and diverse. However, it is still possible to observe that some m6A regulator genes have relatively similar expression patterns in both RNA-seq and MeRIP-seq data, as well as in terms of expression levels, such as TaHAKAI1-A and TaHAKAI1-B among m6A writers, TaALKBH3B and TaALKBH4B among m6A erasers, and TaECT21 and TaECT31 among m6A readers. To further understand the temporal changes of m6A regulators after virulent Pst infection, we selected four writers, two erasers, and three readers for time-dependent expression analysis (Figure 1A). To test the changes in expression levels of m6A regulators after Pst infection, wheat materials without infection were collected as CK (also noted as 0 dpi) and at 0.25, 1, 4, 7, and 10 dpi after infection, with three biological replicates each. sqRT-PCR (semi-quantitative reverse transcription and polymerase chain reaction) products of the PsEF1 gene in each sample showed that the infection was successful, with the intensity of the bands increasing over time and consistent across the three biological replicates (Figure 7B). To determine the relative expression levels of the PsEF1 gene, the intensity of the bands in Figure 7B was analyzed using ImageJ (Figure 7C and Figure S1B). In the early stages of infection (0.25, 1 dpi), the expression levels of PsEF1 were lower, indicating no significant proliferation of Pst, but to some extent, transcriptomic activity was already present due to the formation of haustoria. As the infection progressed, the amplification of PsEF1 significantly increased, reflecting the clear increase in the biomass of stripe rust with the duration of infection, indicating the successful colonization in the host by overcoming its immune system.
Figure 7D–L show the relative expression levels of nine m6A regulators over time as the infection progresses using the qRT-PCR (Real-Time Quantitative Reverse Transcription PCR) assay, revealing that different m6A regulators exhibit various dynamic response patterns during the wheat response to Pst infection. For the two eraser genes, TaALKBH11B and TaALKBH4B (Figure 7H,I), it was observed that TaALKBH11B showed a clear response and induced expression in the early and initial stages of infection (0.25, 1, 4 dpi), while TaALKBH4B was clearly responsive and induced expression during the colonization phase (7, 10 dpi) of Pst inoculation. The three m6A reader genes detected all showed lower relative expression levels at 1 dpi (Figure 7J–L), which may be due to the lower response of m6A writer genes at 1 dpi (Figure 7D,F,G), indicating that the expression of m6A readers is subject to feedback regulation by the function of m6A writers. Unlike other writer genes, TaFIP37-2D continuously increased its induced expression throughout the infection stages (Figure 7E).

3. Discussion

This study combines RNA-seq and MeRIP-seq to analyze genes with significant differences in mRNA and m6A modifications in wheat after compatible Pst interaction transcriptome-wide. We also examined the temporal changes in wheat m6A regulators after Pst infection through qRT-PCR analysis (Figure 1A). From the perspective of m6A regulation, this study elucidates how wheat responds to Pst infection and how stripe rust interferes with wheat growth, metabolism, and immune pathways to successfully colonize. This research provides an important theoretical basis for understanding the disruption of m6A regulatory mechanisms and the genes regulated by the m6A mechanism during the infection process of biotrophic fungi in plants.
Regarding the omics studies on the interaction between Pst and wheat, with the development of sequencing technologies, improved genome assembly and annotation for wheat [41] and Pst (https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=27350 (accessed on 15 February 2024)) have led to the accumulation of data. Our research previously reviewed some of the transcriptome, microarray, and proteome omics analyses of wheat and Pst, analyzing potential Pst effector candidates [42,43]. Recently, the study of the pan-genome of Pst has laid the foundation for the population genetics and comparative genomics research of the Pst population [44]. Currently, dual RNA sequencing seems more suitable for samples that cannot be separated from the host wheat after Pst infection. This sequencing technology uses unique molecular identifier (UMI) tags for sequencing, which involves several more rounds of amplification of the reverse-transcribed templates than other sequencing methods. This approach optimizes the issue where plant RNA is abundant and pathogen RNA is relatively less abundant and harder to detect in plant samples after pathogen infection [45]. It is known that RNA has over 100 modifications, with m6A methylation being the most common modification in eukaryotic mRNA. Studies have found that m6A modification is involved in plant responses to biotic or abiotic stress, but there have been no reports on wheat response to Pst infection in this aspect. Therefore, this study combines RNA-seq and MeRIP-seq data to analyze the RNA response and m6A modification in wheat after virulent infection with stripe rust.
Research on photosynthesis and its related genes in the interaction between stripe rust and wheat has been reported from various perspectives. Chen et al. (2015) [46] found that after infection by stripe rust, resistant wheat had higher levels of D1 protein and light-harvesting complex II (LHCII) accumulation compared to susceptible wheat. The light-harvesting antenna protein CP29 in both resistant and susceptible wheat was phosphorylated under stripe rust infection, with more significant phosphorylation in resistant wheat. The thylakoid membranes in susceptible wheat suffered more extensive damage [46]. In this study, we also found a significant decrease in the expression of many light-harvesting antenna protein genes. Wang et al. (2019) discovered that the peripheral protein of photosystem II, PsbO, could be regulated by phosphorylation by a protein kinase WKS1 encoded by YR36, inhibiting ascorbate peroxidase, leading to the accumulation of substrate hydrogen peroxide, initially inhibiting the growth of the pathogen. Subsequently, photosystem II, without the protection of PsbO, acts as a source of hydrogen peroxide production in the broad-spectrum disease resistance process, inhibiting the reproduction of stripe rust and demonstrating plant resistance [47]. During the normal photosynthesis process, zeaxanthin epoxidase 1 (ZEP1) can transfer singlet oxygen to zeaxanthin, producing violaxanthin, avoiding the accumulation of reactive oxygen species and protecting the photosynthetic system. Chang et al. (2023) found that ZEP1 mutation could inhibit the growth of stripe rust, promote the accumulation of reactive oxygen species at the invasion site, and significantly enhance wheat resistance [48]. Xu et al. (2019) discovered that the effector protein Pst_12806 could be transported to chloroplasts, where it interacts with the Rieske structure of the plant cytochrome b6f protein subunit TaISP (iron–sulfur protein), weakening the electron transfer capability of the TaISP Rieske domain during photosynthesis, thereby disrupting plant photosynthesis, inhibiting chloroplast-mediated reactive oxygen species, and further promoting pathogen colonization in the plant [49]. In this study, we found three genes of Rieske-type iron–sulfur protein enriched in the level 2 GO term annotated as an immune system process, and their m6A modifications showed significant changes after stripe rust infection. Andac et al. (2020) also identified a stripe rust effector, PstCTE1, which, despite lacking a clear conserved transport signal peptide as predicted by TargetP and ChloroP, was found to localize in chloroplasts regardless of whether the N-terminus was tagged [50]. The discovery of Pst effectors Pst_12806 and PstCTE1 suggests that disrupting host wheat photosynthetic machinery proteins is one of the strategies for Pst to infect the host. The above research on photosynthesis-related mechanisms in the interaction between stripe rust and wheat mostly focused on post-translational modifications. In contrast, this study provides new insights into the level of post-transcriptional regulation, specifically m6A modification of RNA. We found that photosynthesis-related genes with m6A modifications mostly showed hypo-methylation after stripe rust infection (Table S24), particularly among the most significantly hypo-methylated genes, most of which are related to photosynthesis (Figure 5B, Table 1 and Table S28). Thus, on one hand, the expression of photosynthesis-related genes may decrease, leading to a reduction in related proteins, and on the other hand, the regulation of m6A modifications may also interfere with the translation level of related proteins, ultimately leading to reduced photosynthesis, decreased substance accumulation, and increased plant susceptibility. Further research on photosynthesis and plant immunity is needed at the molecular level. Photosynthesis is a complex and precise regulatory mechanism, and genetic modification of a few photosynthesis-related genes might enhance disease resistance. However, the balance of the photosynthesis mechanism itself is easily disturbed and broken. Since photosynthesis is closely related to yield and specific agricultural practices do not wish to sacrifice yield for increased disease resistance, seeking a balance between the two is necessary.

4. Materials and Methods

4.1. Plant Materials and Pst Inoculation

The wheat cultivar MingXian169 (MX169) was cultivated in a growth chamber with a light intensity of 200 μmol·m−2·s, temperature maintained at 18 °C, and humidity at 80%. The light–dark cycle was set at 16 h of light and 8 h of darkness. When MX169 reached the two-leaf stage with fully expanded leaves, it was used for inoculation with Pst spores. Spores of Chinese yellow rust race CYR32 stored at −80 °C were used, and they were incubated in a 42 °C water bath for 10 min. Wheat leaves were sprayed with 2 mg/mL of Novec 7100, followed by the inoculation of spores. Infected MX169 were kept in darkness at 12 °C with 100% humidity for 1 d, then transferred to a growth chamber under normal conditions. After 2 weeks, the leaves showed abundant fresh spores, which were used to inoculate wheat cultivar Avocet/Yr7 (Avocet/Yr7 the stripe rust Avocet differential lines were originally generated by C.R. Wellings of Plant Breeding Institute, Narrabri, Australia). Inoculation conditions were the same as those for MX169. Infected Avocet/Yr7 were collected at 0.25, 1, 4, 7, and 10 dpi. Avocet/Yr7 only sprayed with 2 mg/mL of Novec 7100 solution without inoculation were used as control (CK) samples. Three biological replicates were collected for both CK and infected samples with different time points and then stored at −80 °C for subsequent RNA extraction.

4.2. RNA Extraction and Library Construction for RNA-Seq and MeRIP-Seq

Total RNA was extracted from samples stored at −80 °C using the TRIzol reagent (Life Technologies, Carlsbad, CA, USA, cat15596026) with the user guide of the manufacturer involving the homogenization of leaf tissue, phase separation with chloroform, RNA precipitation, and purification. Extracted RNA quality was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), ensuring suitability for library preparation and sequencing. Then, mRNA enriched by Oligo(dT) magnetic beads was fragmented into short fragments using a fragmentation buffer. Fragmented RNA was divided into two parts, one of which was used as the input (no immunoprecipitation experiment was performed). The other RNA enriched with an m6A-specific antibody was used as the IP. The input RNA and IP RNA were reversely transcribed into cDNA with random primers using NEBNext® UltraTM RNA Library Prep Kit for Illumina (New England Biolabs, MA, USA catE7530). Next, the cDNA fragments were end-repaired, “A” base-tailed, and ligated to Illumina sequencing universal adapters. Finally, the qualified cDNA library with PCR amplification was sequenced using Illumina Novaseq6000 by GENE DENOVO Biotechnology Co., Ltd. (Guangzhou, China).

4.3. Data Analysis of Sequencing

In preprocessing RNA-seq and MeRIP-seq raw data, FASTQ v0.18.0 [51] was employed as a trimming adapter, filtering reads exceeding 10% N content, removing poly A sequences, and discarding low-quality reads. High-quality clean reads were aligned to the species ribosome database using Bowtie2 [52], with ribosomal RNA reads subsequently removed, and HISAT2 [53] was used for genome alignment. Differential gene expression was identified using criteria of FDR < 0.05 and |log2(FC)| > 1, classifying genes with significant changes in expression. PCA was performed with the R package models (http://www.r-project.org/ (accessed on 8 October 2022)). Differential RNA methylation was analyzed with DiffBind v1.20.0 [54], identifying peaks as significantly differential with p-value < 0.05 and |log2(FC)| > 1. StringTie v1.3.1 was utilized to assemble RNA-seq data from each sample against a known genome [55,56]. RSEM calculated TPM (transcript per kilobase per million mapped reads) values to quantify RNA transcript expression levels, enabling accurate across-sample comparisons [57]. GO [58] (http://www.geneontology.org/ (accessed on 11 October 2022)) and KEGG [59] enrichment was analyzed to clarify the function of differentially expressed genes. Using MACS2 [60] (version 2.1.2), read-enriched regions in MeRIP-seq data were identified through a dynamic Poisson distribution, defining peaks with a q-value < 0.05 and retaining those with over 50% overlap in at least two of the biological duplicates for further analysis. By analyzing the genomic location and gene annotation of peaks, related genes were confirmed, and the distribution of peaks across functional regions like 5′ UTR, CDS, and 3′ UTR was assessed [61]. Motif analysis was performed by HOMER [62], and significant motif identification in peak-associated transcripts was carried out by the MEME Suite (http://meme-suite.org/meme/ (accessed on 6 October 2022)) and DREME (http://meme-suite.org/tools/dreme (accessed on 6 October 2022)). The clustering heatmap (http://cloud.oebiotech.com/task/detail/heatmap/?id=31 (accessed on 11 February 2024)) and quadrant image (www.omicstudio.cn/tool/31 (accessed on 21 January 2024)) were formed using online software.

4.4. sqRT-PCR and qRT-PCR

After adjusting the concentration of RNA extracted from CK and 0.25, 1, 4, 7, and 10 dpi infected samples were ensured to be the same as 350 ng/μL. The genomic DNA was removed and synthesis of the first cDNA strand was performed according to the specifications of SweScript RT II First Strand cDNA Synthesis Kit (Servicebio, Wuhan China catG3333). A 10 μL genome removal reaction mixture was prepared with 8 μL RNA adjusted to 350 ng/μL, 1 μL 10×gDNA remover buffer, and 1 μL gDNA remover, lightly mixed and briefly centrifuged, then incubated at 37 °C for 2 min and placed on ice for later use. Then, 20 μL reverse transcription reaction mixture was prepared, including 10 μL genome removal mixture, 3 μL nuclease-free water, 4 μL 5 × reaction buffer, 1 μL 100 μmol/L oligo (dT)18 primer, 1 μL 100 μmol/L random hexamer primer, and 1 μL SweScript RT II enzyme mix, lightly mixed and briefly centrifuged, then incubated at 25 °C for 5 min, 55 °C for 15 min, and 85 °C for 5 s. The synthesized cDNA was stored at −20 °C.
A 20 μL sqRT-PCR mixture was prepared, which included 1 μL cDNA of each sample as a template, 10 μL 2 × Fast sTaq PCR master mix (Servicebio, Wuhan China catG3304), 0.8 μL 10 μM forward primer, 0.8 μL 10 μM reverse primer, and 7.4 μL nuclease-free water. The primers for gene PsEF1 (Puccinia striiformis elongation factor 1, accession number: KNE93481) are shown in Table S33. After being gently mixed and briefly centrifuged, the sqRT-PCR mixture was denatured at 95 °C for 2 min, then 30 cycles were amplified according to the setting of “denaturation at 95 °C for 10 s, annealing at 57 °C for 30 s, extension at 72 °C for 10 s” and finally extended at 72 °C for 10 min. Then, 5 μL of PCR product was loaded on 2.5% agarose gel containing 4S Green Plus Nucleic Acid Stain, and electrophoresis was performed using the setting “100 V, 45 min”. Grey-scale analysis of electrophoresis images was performed by ImageJ to indicate the relative expression level of PsEF1 in each sample.
A 15 μL qRT-PCR mixture was prepared, which contained 1 μL cDNA of each sample as a template, 7.5 μL universal blue SYBR green qPCR master mix (Servicebio, Wuhan China G3326), 0.3 μL of 10 μM forward primer, 0.3 μL 10 μM reverse primer, and 5.9 μL nuclease-free water. The primers of reference gene TaEF1α (Triticum aestivum elongation factor 1 alpha, accession number: Q03033) and m6A regulator genes are shown in Table S33. After the qRT-PCR mixture was gently mixed and briefly centrifuged, it was placed in the CFX Connect fluorescent quantitative PCR detection system. The detection mode was adjusted to SYBR, and the denaturation was set at 95 °C for 30 s, and then 45 cycles were amplified according to the program of “denaturation at 95 °C for 15 s, annealing and extension at 60 °C for 30 s”. Finally, the melting curve was made at intervals of 0.5 °C for 5 s between 65 °C and 95 °C to determine the specificity of PCR amplification. Then, the amplified gene name, sample name, and operation repetition for each reaction well were set in CFX Maestro v2.2, and TaEF1α was set as the reference gene. The relative expression levels of each tested gene were calculated by the 2−ΔΔCt method [63]. Each test was repeated five times.

5. Conclusions

In the differential gene analysis of RNA-seq, we found that among the significantly up-regulated 4100 genes (Figure 2C), 1010 wheat genes had similar expression patterns after stripe rust infection (Figure 2D, Table S7), including 18 genes annotated in the KEGG pathway as “Plant-pathogen interaction”, indicating that disease-resistance-related genes synergistically respond to resist the stripe rust assault. Among these eighteen genes, two are NBS-LRR disease resistance proteins (TraesCS1A03G0059400, TraesCS7D03G0011300), with TraesCS1A03G0059400 annotated in UniProt as disease resistance protein RPM1 and TraesCS7D03G0011300 as Rx N-terminal domain-containing protein. Two of them are pathogenesis-related protein 1 (TraesCS5B03G1087600, TraesCS7D03G0450000), and also, a significantly up-regulated PR1 gene was found among the overall up-regulated genes (Figure 2E), indicating that the TaPR1 gene is indeed an important gene in response to stripe rust, worthy of further research [37]. Six genes are annotated as “Calmodulin”, “Calcium-dependent protein kinase”, and “Cyclic nucleotide-gated channel”, indicating that calcium regulation plays a significant role in plant immune responses to affect the invasion of pathogens. Additionally, among the genes with the most significant decrease in m6A modification, we found that the most significantly down-regulated TraesCS1A03G0058100 is annotated as an NBS-LRR disease-resistance protein. NBS-LRR genes play a crucial role in plant resistance against pathogens. This might be the first report of NBS-LRR genes showing hypo-methylation after infection in a Pst–wheat compatible interaction, suggesting that the m6A mechanism can directly regulate key genes in the plant immune pathway and affect the plant response to pathogens.
We observed how host wheat actively responds to stripe rust infection through the perspective of up-regulated genes suppressing the host wheat immune system for the successful invasion and colonization of Pst. We also turned our attention to the analysis of down-regulated genes. Among these genes, the two most significantly down-regulated are related to photosynthesis: TraesCS1A03G0499300 (protein PAM68, chloroplast) and TraesCS5D03G0016200LC (photosystem I assembly protein Ycf4). When we further analyzed the overall significantly up-regulated and down-regulated genes in the enrichment of GO terms and KEGG pathways (Figure 3), we found that in the cellular component category of GO terms, there were strikingly multiple photosynthesis-related GO terms (Figure 3B, Table S10). Similarly, in the KEGG pathway enrichment, KO00196 (photosynthesis-antenna proteins) had the highest rich factor and the lowest q-value, with all 80 enriched genes being down-regulated (Figure 3D,E). This indicates that after compatible Pst infection, the photosynthesis-related genes in host wheat were significantly disrupted, and their regulatory levels were disturbed. The reduction in photosynthesis affects the accumulation of substances within the plant and the regulatory level of secondary biomass metabolism, as we also found that the regulatory levels of pathway genes in basic metabolism and secondary metabolism were reduced (for example, KO00710: carbon fixation in photosynthetic organisms; KO00630: glyoxylate and dicarboxylate metabolism) (Figure 3D,E). This may lead to a reduction in wheat resistance during the early stage of stripe rust infection. In the later stages of infection, the reduction in photosynthesis and substance accumulation forces the Pst to absorb nutrients from the host nearby leaf tissue, further expanding the invasion and colonization of Pst.
Additionally, this study also focused on the expression changes in m6A regulators in wheat after stripe rust infection. Among the 90 wheat m6A regulator genes summarized and annotated by Yue et al. (2019) [40], 28 m6A regulator genes were detected in this study (Figure 7A), indicating that these m6A regulator genes are closely associated with the response to stripe rust infection and are involved in m6A regulation. At the same time, the m6A modification levels of some m6A regulator genes themselves also showed significant changes after infection, such as hyper-methylation in TaECT21 and TaECT31 and hypo-methylation in TaHAKAI1-A and TaHAKAI1-B. This suggests that the regulation of m6A modification by m6A regulators is also under the control of the m6A mechanism, implying that the m6A mechanism is a precisely regulated biological process. The regulatory genes show time-point-dependent expressions; different genes have divergent dynamic response patterns, suggesting that the m6A methylation mechanism dynamically responds to the infection process. It is noticed that at 1 dpi, the expression levels of all tested m6A readers were lower. The observation is in accordance with the lower expression levels of m6A writers (Figure 7D–L), which is supported by previous results [40], indicating that most m6A writer genes also had lower expression levels in the early stages of stripe rust infection. Using the rust expression browser (http://www.rust-expression.com (accessed on 24 February 2024)) [64], we examined the expression levels of eight searchable m6A writer genes (TaVIR-D, TaVIR-A, TaMTA-D, TaMTB-D, TaMTA-B, TaMTB-B, TaFIP37-2D, TaFIP37-2A) in the wheat variety “Vuka” after infection with Pst isolates 87/66 at 0, 1, 2, 3, 5, 7, 9, and 11 dpi [65] (Figure S4) and found that most m6A writer genes also had lower expression levels in the early stages of stripe rust infection (2, 3 dpi).
Overall, we have produced an immense amount of data determining the genes undergoing m6A methylation for the first time in the wheat–Pst pathosystem. There is no doubt that these data will facilitate further in-depth analyses in system biology.

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/plants13070982/s1, Figure S1. Uncropped agarose gel (2.5%) electrophoresis pattern of semi-quantitative reverse transcription PCR (sqRT-PCR) product for PsEF1; Figure S2. Distribution and coverage of m6A peaks from CK sample (A) and infected sample (B) on wheat chromosome; Figure S3. Number of m6A peaks distributed within 5′UTR (untranslated region) (A), start codon (B), CDS (coding sequence) (C), stop codon (D), and 3′UTR (E) of genes in CK sample and infected sample; Figure S4. The heatmap visualization of the expression levels of m6A writer genes in the wheat variety “Vuka” when infected with Pst isolate 87/66, detected at 0, 1, 2, 3, 5, 7, 9, and 11 dpi, as logged in the rust expression browser. Table S1. Data filtering statistics of RNA-seq; Table S2. Base information statistics of RNA-seq; Table S3. Ribosomal RNA mapping statistics of RNA-seq; Table S4. Reference genome mapping statistics of RNA-seq; Table S5. All gene annotations of RNA-seq; Table S6. Annotation of all genes with significant differential expression in RNA-seq; Table S7. Heatmap clustering of differentially expressed genes (DEGs); Table S8. Gene Ontology (GO) terms enrichment in biological processes of differentially expressed genes (DEGs); Table S9. Gene Ontology (GO) terms enrichment in molecular function of differentially expressed genes (DEGs); Table S10. Gene Ontology (GO) term enrichment in cellular component of differentially expressed genes (DEGs); Table S11. Level 2 Gene Ontology (GO) terms of differentially expressed genes (DEGs); Table S12. KEGG enrichment analysis of differentially expressed genes (DEGs); Table S13. Data filtering statistics of MeRIP-seq; Table S14. Base information statistics of MeRIP-seq; Table S15. Ribosomal RNA mapping and filtering statistics of MeRIP-seq; Table S16. Reference genome mapping statistics of MeRIP-seq; Table S17. CK sample peak identification of MeRIP-seq; Table S18. Infected sample peak identification of MeRIP-seq; Table S19. CK sample peak-related gene annotation of MeRIP-seq; Table S20. Infected sample peak-related gene annotation of MeRIP-seq; Table S21. Statistics of peak number in each group; Table S22. Pst-induced methylation sites and associated genes in infected sample motif; Table S23. All peak annotations of MeRIP-seq after merging peaks; Table S24. Annotation of differentially m6A-modified peaks; Table S25. Gene Ontology (GO) terms enrichment in biological processes of differentially m6A-modified peaks (DMPs); Table S26. Gene Ontology (GO) terms enrichment in cellular component of differentially m6A-modified peaks (DMPs); Table S27. Gene Ontology (GO) terms enrichment in molecular function of differentially m6A-modified peaks (DMPs); Table S28. Level 2 Gene Ontology (GO) terms of differentially m6A-modified peaks (DMPs); Table S29. KEGG enrichment analysis of differentially m6A-modified peaks (DMPs); Table S30. Joint statistics on fold changes of m6A-modified peak-related genes and transcripts in infected leaves in comparison to CK; Table S31. Significant fold changes of both m6A-modified peak-related genes and transcripts in infected leaves in comparison to CK; Table S32. Accession numbers, amino acid sequences, and name codes of m6A writers, erasers, and readers of Triticum aestivum; Table S33. The primers used in this work.

Author Contributions

Conceptualization, M.S.A.; Methodology, E.N.C.; Investigation, E.N.C.; Writing—Original Draft, E.N.C. and N.W.; Writing—Review and Editing, M.S.A., N.W. and E.N.C.; Supervision, M.S.A.; Project Administration, M.S.A. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by a grant from Dalian University of Technology (DUT18RC(3)050).

Data Availability Statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: National Center for Biotechnology Information (NCBI) BioProject, PRJNA1072674.

Acknowledgments

We thank Zhensheng Kang from Northwest A&F University for providing Pst race CYR32 and wheat cultivar MingXian169. We thank Taiguo Liu from the Chinese Academy of Agricultural Sciences for providing wheat Avocet/Yr7 used in this work.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Boccaletto, P.; Stefaniak, F.; Ray, A.; Cappannini, A.; Mukherjee, S.; Purta, E.; Kurkowska, M.; Shirvanizadeh, N.; Destefanis, E.; Groza, P.; et al. MODOMICS: A database of RNA modification pathways. 2021 update. Nucleic Acids Res. 2021, 50, D231–D235. [Google Scholar] [CrossRef] [PubMed]
  2. Zhang, Y.; Lu, L.; Li, X. Detection technologies for RNA modifications. Exp. Mol. Med. 2022, 54, 1601–1616. [Google Scholar] [CrossRef] [PubMed]
  3. Bi, Z.; Liu, Y.; Zhao, Y.; Yao, Y.; Wu, R.; Liu, Q.; Wang, Y.; Wang, X. A dynamic reversible RNA N6-methyladenosine modification: Current status and perspectives. J. Cell. Physiol. 2019, 234, 7948–7956. [Google Scholar] [CrossRef] [PubMed]
  4. Lin, H.; Shi, T.; Zhang, Y.; He, C.; Zhang, Q.; Mo, Z.; Pan, W.; Nie, X. Genome-wide identification, expression and evolution analysis of m6A writers, readers and erasers in Aegilops_tauschii. Plants 2023, 12, 2747. [Google Scholar] [CrossRef] [PubMed]
  5. Shen, L.; Ma, J.; Li, P.; Wu, Y.; Yu, H. Recent advances in the plant epitranscriptome. Genome Biol. 2023, 24, 43. [Google Scholar] [CrossRef] [PubMed]
  6. Shinde, H.; Dudhate, A.; Kadam, U.S.; Hong, J.C. RNA methylation in plants: An overview. Front. Plant Sci. 2023, 14, 1132959. [Google Scholar] [CrossRef] [PubMed]
  7. Rajecka, V.; Skalicky, T.; Vanacova, S. The role of RNA adenosine demethylases in the control of gene expression. Biochim. Biophys. Acta Gene Regul. Mech. 2019, 1862, 343–355. [Google Scholar] [CrossRef] [PubMed]
  8. Shen, L.; Liang, Z.; Wong, C.E.; Yu, H. Messenger RNA modifications in plants. Trends Plant Sci. 2019, 24, 328–341. [Google Scholar] [CrossRef] [PubMed]
  9. Yu, Q.; Liu, S.; Yu, L.; Xiao, Y.; Zhang, S.; Wang, X.; Xu, Y.; Yu, H.; Li, Y.; Yang, J.; et al. RNA demethylation increases the yield and biomass of rice and potato plants in field trials. Nat. Biotechnol. 2021, 39, 1581–1588. [Google Scholar] [CrossRef]
  10. Bai, B. An update on principles of m6A targeting. Trends Plant Sci. 2022, 27, 224–226. [Google Scholar] [CrossRef]
  11. Ma, W.; Cui, S.; Lu, Z.; Yan, X.; Cai, L.; Lu, Y.; Cai, K.; Zhou, H.; Ma, R.; Zhou, S.; et al. YTH domain proteins play an essential role in rice growth and stress response. Plants 2022, 11, 2206. [Google Scholar] [CrossRef]
  12. Tang, J.; Chen, S.; Jia, G. Detection, regulation, and functions of RNA N6-methyladenosine modification in plants. Plant Commun. 2023, 4, 100546. [Google Scholar] [CrossRef]
  13. Wang, K.; Peng, J.; Yi, C. The m6A consensus motif provides a paradigm of epitranscriptomic studies. Biochemistry 2021, 60, 3410–3412. [Google Scholar] [CrossRef] [PubMed]
  14. Huang, T.; He, W.J.; Li, C.; Zhang, J.B.; Liao, Y.C.; Song, B.; Yang, P. Transcriptome-wide analyses of RNA m6A methylation in hexaploid wheat reveal its roles in mRNA translation regulation. Front. Plant Sci. 2022, 13, 917335. [Google Scholar] [CrossRef] [PubMed]
  15. Arribas-Hernández, L.; Brodersen, P. Occurrence and functions of m6A and other covalent modifications in plant mRNA. Plant Physiol. 2019, 182, 79–96. [Google Scholar] [CrossRef] [PubMed]
  16. Zhao, J.; Zhang, C.; Li, S.; Yuan, M.; Mu, W.; Yang, J.; Ma, Y.; Guan, C.; Ma, C. Changes in m6A RNA methylation are associated with male sterility in wolfberry. BMC Plant Biol. 2023, 23, 456. [Google Scholar] [CrossRef] [PubMed]
  17. Chuong, N.N.; Doan, P.P.T.; Wang, L.; Kim, J.H.; Kim, J. Current insights into m6A RNA methylation and its emerging role in plant circadian clock. Plants 2023, 12, 624. [Google Scholar] [CrossRef] [PubMed]
  18. Anderson, S.J.; Kramer, M.C.; Gosai, S.J.; Yu, X.; Vandivier, L.E.; Nelson, A.D.L.; Anderson, Z.D.; Beilstein, M.A.; Fray, R.G.; Lyons, E.; et al. N6-methyladenosine inhibits local ribonucleolytic cleavage to stabilize mRNAs in Arabidopsis. Cell Rep. 2018, 25, 1146–1157.e3. [Google Scholar] [CrossRef] [PubMed]
  19. Miao, Z.; Zhang, T.; Qi, Y.; Song, J.; Han, Z.; Ma, C. Evolution of the RNA N6-methyladenosine methylome mediated by genomic duplication. Plant Physiol. 2020, 182, 345–360. [Google Scholar] [CrossRef]
  20. Meyer, K.D.; Saletore, Y.; Zumbo, P.; Elemento, O.; Mason, C.E.; Jaffrey, S.R. Comprehensive analysis of mRNA methylation reveals enrichment in 3′ UTRs and near stop codons. Cell 2012, 149, 1635–1646. [Google Scholar] [CrossRef]
  21. Hu, J.; Manduzio, S.; Kang, H. Epitranscriptomic RNA methylation in plant development and abiotic stress responses. Front. Plant Sci. 2019, 10, 500. [Google Scholar] [CrossRef] [PubMed]
  22. Zhou, L.; Gao, G.; Tang, R.; Wang, W.; Wang, Y.; Tian, S.; Qin, G. m6A-mediated regulation of crop development and stress responses. Plant Biotechnol. J. 2022, 20, 1447–1455. [Google Scholar] [CrossRef] [PubMed]
  23. Rudy, E.; Grabsztunowicz, M.; Arasimowicz-Jelonek, M.; Tanwar, U.K.; Maciorowska, J.; Sobieszczuk-Nowicka, E. N6-methyladenosine (m6A) RNA modification as a metabolic switch between plant cell survival and death in leaf senescence. Front. Plant Sci. 2023, 13, 1064131. [Google Scholar] [CrossRef] [PubMed]
  24. Martínez-Pérez, M.; Aparicio, F.; López-Gresa, M.P.; Bellés, J.M.; Sánchez-Navarro, J.A.; Pallás, V. Arabidopsis m6A demethylase activity modulates viral infection of a plant virus and the m6A abundance in its genomic RNAs. Proc. Natl. Acad. Sci. USA 2017, 114, 10755–10760. [Google Scholar] [CrossRef] [PubMed]
  25. Han, C.; Zhang, F.; Qiao, X.; Zhao, Y.; Qiao, Q.; Huang, X.; Zhang, S. Multi-Omics analysis reveals the dynamic changes of RNA N6-Methyladenosine in pear (Pyrus bretschneideri) defense responses to Erwinia amylovora pathogen infection. Front. Microbiol. 2022, 12, 803512. [Google Scholar] [CrossRef]
  26. Rathod, V.; Hamid, R.; Tomar, R.S.; Patel, R.; Padhiyar, S.; Kheni, J.; Thirumalaisamy, P.P.; Munshi, N.S. Comparative RNA-Seq profiling of a resistant and susceptible peanut (Arachis hypogaea) genotypes in response to leaf rust infection caused by Puccinia arachidis. 3 Biotech. 2020, 10, 284. [Google Scholar] [CrossRef] [PubMed]
  27. Yang, Y.; Chen, T.; Dai, X.; Yang, D.; Wu, Y.; Chen, H.; Zheng, Y.; Zhi, Q.; Wan, X.; Tan, X. Comparative transcriptome analysis revealed molecular mechanisms of peanut leaves responding to Ralstonia solanacearum and its type III secretion system mutant. Front. Microbiol. 2022, 13, 998817. [Google Scholar] [CrossRef] [PubMed]
  28. Wang, H.; Lei, Y.; Wan, L.; Yan, L.; Lv, J.; Dai, X.; Ren, X.; Guo, W.; Jiang, H.; Liao, B. Comparative transcript profiling of resistant and susceptible peanut post-harvest seeds in response to aflatoxin production by Aspergillus flavus. BMC Plant Biol. 2016, 16, 54. [Google Scholar] [CrossRef]
  29. Jiang, C.; Zhang, H.; Ren, J.; Dong, J.; Zhao, X.; Wang, X.; Wang, J.; Zhong, C.; Zhao, S.; Liu, X.; et al. Comparative Transcriptome-based mining and expression profiling of transcription factors related to cold tolerance in peanut. Int. J. Mol. Sci. 2020, 21, 1921. [Google Scholar] [CrossRef]
  30. Chen, Y.; Ren, X.; Zhou, X.; Huang, L.; Yan, L.; Lei, Y.; Liao, B.; Huang, J.; Huang, S.; Wei, W.; et al. Dynamics in the resistant and susceptible peanut (Arachis hypogaea L.) root transcriptome on infection with the Ralstonia solanacearum. BMC Genomics. 2014, 15, 1078. [Google Scholar] [CrossRef]
  31. Xiao, D.; Li, X.; Zhou, Y.Y.; Wei, L.; Keovongkod, C.; He, H.; Zhan, J.; Wang, A.Q.; He, L.F. Transcriptome analysis reveals significant difference in gene expression and pathways between two peanut cultivars under Al stress. Gene 2021, 781, 145535. [Google Scholar] [CrossRef] [PubMed]
  32. Tang, Y.; Wang, C.; Yang, G.; Feng, T.; Gao, H.; Wang, X.; Chi, X.; Xu, Y.; Wu, Q.; Chen, D. Identification of chilling-responsive transcripts in peanut (Arachis hypogaea L.). Electron. J. Biotechnol. 2011, 14, 5. [Google Scholar]
  33. Kumar, K.R.R.; Kirti, P.B. Differential gene expression in Arachis diogoi upon interaction with peanut late leaf spot pathogen, Phaeoisariopsis personata and characterization of a pathogen induced cyclophilin. Plant Mol. Biol. 2011, 75, 497–513. [Google Scholar] [CrossRef]
  34. Kumar, D.; Kirti, P.B. The genus Arachis: An excellent resource for studies on differential gene expression for stress tolerance. Front. Plant Sci. 2023, 14, 1275854. [Google Scholar] [CrossRef] [PubMed]
  35. He, C.; Zhang, Y.; Zhou, W.; Guo, Q.; Bai, B.; Shen, S.; Huang, G. Study on stripe rust (Puccinia striiformis) effect on grain filling and seed morphology building of special winter wheat germplasm Huixianhong. PLoS ONE 2019, 14, e0215066. [Google Scholar] [CrossRef]
  36. Zheng, W.; Huang, L.; Huang, J.; Wang, X.; Chen, X.; Zhao, J.; Guo, J.; Zhuang, H.; Qiu, C.; Liu, J.; et al. High genome heterozygosity and endemic genetic recombination in the wheat stripe rust fungus. Nat. Commun. 2013, 4, 2673. [Google Scholar] [CrossRef] [PubMed]
  37. Liu, R.; Lu, J.; Xing, J.; Xue, L.; Wu, Y.; Zhang, L. Characterization and functional analyses of wheat TaPR1 genes in response to stripe rust fungal infection. Sci. Rep. 2023, 13, 3362. [Google Scholar] [CrossRef]
  38. Meng, J.; Lu, Z.; Liu, H.; Zhang, L.; Zhang, S.; Chen, Y.; Rao, M.K.; Huang, Y. A protocol for RNA methylation differential analysis with MeRIP-Seq data and exomePeak R/Bioconductor package. Methods 2014, 69, 274–281. [Google Scholar] [CrossRef] [PubMed]
  39. Prall, W.; Ganguly, D.R.; Gregory, B.D. The covalent nucleotide modifications within plant mRNAs: What we know, how we find them, and what should be done in the future. Plant Cell 2023, 35, 1801–1816. [Google Scholar] [CrossRef] [PubMed]
  40. Yue, H.; Nie, X.; Yan, Z.; Weining, S. N6-methyladenosine regulatory machinery in plants: Composition, function and evolution. Plant Biotechnol. J. 2019, 17, 1194–1208. [Google Scholar] [CrossRef]
  41. Appels, R.; Eversole, K.; Feuillet, C.; Gallagher, D. The Wheat Genome, 1st ed.; Springer: Cham, Switzerland, 2023. [Google Scholar]
  42. Ozketen, A.C.; Andac-Ozketen, A.; Dagvadorj, B.; Demiralay, B.; Akkaya, M.S. In-depth secretome analysis of Puccinia striiformis f. sp. tritici in infected wheat uncovers effector functions. Biosci. Rep. 2020, 40, BSR20201188. [Google Scholar] [CrossRef]
  43. Wu, N.; Ozketen, A.C.; Cheng, Y.; Jiang, W.; Zhou, X.; Zhao, X.; Guan, Y.; Xiang, Z.; Akkaya, M.S. Puccinia striiformis f. sp. tritici effectors in wheat immune responses. Front. Plant Sci. 2022, 13, 1012216. [Google Scholar] [CrossRef]
  44. Wang, J.; Peng, Y.; Xu, Y.; Li, Z.; Zhan, G.; Kang, Z.; Zhao, J. Pan-genome analysis reveals a highly plastic genome and extensive secreted protein polymorphism in Puccinia striiformis f. sp. tritici. J. Genet. Genom. 2023. [Google Scholar] [CrossRef]
  45. Qu, Z.; Zhang, H.; Wang, Q.; Zhao, H.; Liu, X.; Fu, Y.; Lin, Y.; Xie, J.; Cheng, J.; Li, B.; et al. Exploring the symbiotic mechanism of a virus-mediated endophytic fungus in its host by dual unique molecular identifier-RNA sequencing. mSystems 2021, 6, e0081421. [Google Scholar] [CrossRef] [PubMed]
  46. Chen, Y.E.; Cui, J.M.; Su, Y.Q.; Yuan, S.; Yuan, M.; Zhang, H.Y. Influence of stripe rust infection on the photosynthetic characteristics and antioxidant system of susceptible and resistant wheat cultivars at the adult plant stage. Front. Plant Sci. 2015, 6, 779. [Google Scholar] [CrossRef]
  47. Wang, S.; Li, Q.P.; Wang, J.; Yan, Y.; Zhang, G.L.; Yan, Y.; Zhang, H.; Wu, J.; Chen, F.; Wang, X.; et al. YR36/WKS1-mediated phosphorylation of PsbO, an extrinsic member of photosystem II, inhibits photosynthesis and confers stripe rust resistance in wheat. Mol. Plant 2019, 12, 1639–1650. [Google Scholar] [CrossRef]
  48. Chang, C.Y.; Yang, S.X.; Zhang, M.Q.; Guo, Y.T.; Li, X.M.; Yan, Y.; Ding, C.H.; Niu, K.X.; Wang, M.L.; Li, Q.Q.; et al. Suppression of ZEAXANTHIN EPOXIDASE 1 restricts stripe rust growth in wheat. Plant Commun. 2023, 4, 100608. [Google Scholar] [CrossRef] [PubMed]
  49. Xu, Q.; Tang, C.; Wang, X.; Sun, S.; Zhao, J.; Kang, Z.; Wang, X. An effector protein of the wheat stripe rust fungus targets chloroplasts and suppresses chloroplast function. Nat. Commun. 2019, 10, 5571. [Google Scholar] [CrossRef] [PubMed]
  50. Andac, A.; Ozketen, A.C.; Dagvadorj, B.; Akkaya, M.S. An effector of Puccinia striiformis f. sp. tritici targets chloroplasts with a novel and robust targeting signal. Eur. J. Plant Pathol. 2020, 157, 751–765. [Google Scholar] [CrossRef]
  51. Chen, S.; Zhou, Y.; Chen, Y.; Gu, J. fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 2018, 34, i884–i890. [Google Scholar] [CrossRef]
  52. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef]
  53. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef]
  54. Stark, R.; Brown, G. DiffBind: Differential Binding Analysis of ChIP-Seq Peak Data; Bioconductor: Boston, MA, USA, 2011; 100p. [Google Scholar]
  55. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 2015, 33, 290–295. [Google Scholar] [CrossRef] [PubMed]
  56. Pertea, M.; Kim, D.; Pertea, G.M.; Leek, J.T.; Salzberg, S.L. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 2016, 11, 1650–1667. [Google Scholar] [CrossRef] [PubMed]
  57. Li, B.; Dewey, C.N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, 323. [Google Scholar] [CrossRef] [PubMed]
  58. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T.; et al. Gene ontology: Tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 2000, 25, 25–29. [Google Scholar] [CrossRef] [PubMed]
  59. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28, 27–30. [Google Scholar] [CrossRef] [PubMed]
  60. Zhang, Y.; Liu, T.; Meyer, C.A.; Eeckhoute, J.; Johnson, D.S.; Bernstein, B.E.; Nusbaum, C.; Myers, R.M.; Brown, M.; Li, W.; et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008, 9, R137. [Google Scholar] [CrossRef] [PubMed]
  61. Duan, H.C.; Wei, L.H.; Zhang, C.; Wang, Y.; Chen, L.; Lu, Z.; Chen, P.R.; He, C.; Jia, G. ALKBH10B is an RNA N6-methyladenosine demethylase affecting Arabidopsis floral transition. Plant Cell 2017, 29, 2995–3011. [Google Scholar] [CrossRef] [PubMed]
  62. Heinz, S.; Benner, C.; Spann, N.; Bertolino, E.; Lin, Y.C.; Laslo, P.; Cheng, J.X.; Murre, C.; Singh, H.; Glass, C.K. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 2010, 38, 576–589. [Google Scholar] [CrossRef]
  63. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCt method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
  64. Adams, T.M.; Olsson, T.S.G.; Ramírez-González, R.H.; Bryant, R.; Bryson, R.; Campos, P.E.; Fenwick, P.; Feuerhelm, D.; Hayes, C.; Henriksson, T.; et al. Rust expression browser: An open source database for simultaneous analysis of host and pathogen gene expression profiles with expVIP. BMC Genom. 2021, 22, 166. [Google Scholar] [CrossRef] [PubMed]
  65. Dobon, A.; Bunting, D.C.; Cabrera-Quio, L.E.; Uauy, C.; Saunders, D.G. The host-pathogen interaction between wheat and yellow rust induces temporally coordinated waves of gene expression. BMC Genom. 2016, 17, 380. [Google Scholar] [CrossRef] [PubMed]
Plants 13 00982 g001
Figure 1. An overview of the experimental setup in wheat after Pst infection, integrating MeRIP-seq, RNA-seq, and qRT-PCR with phenotype analysis. (A) Schematic of the experimental workflow for MeRIP-seq, RNA-seq, and qRT-PCR in CK and Pst-infected wheat leaves. This diagram illustrates the step-by-step procedure adopted for analyzing both mock inoculation wheat leaves (CK) and Pst-infected wheat leaves (infected) at 7 dpi using NGS technology. The workflow highlights the incorporation of the MeRIP technique to selectively isolate m6A-modified RNA fractions for subsequent sequencing. (B) The phenotype of CK and infected wheat seedlings with Pst infection at 7 dpi. The CK wheat leaves showed no visible symptoms, but the infected wheat leaves showed chlorotic patches and blotches. (C) The magnification of CK and infected wheat leaves in the dotted box in (B). (D) sqRT-PCR assay of PsEF1 in CK and infected wheat leaves at 7 dpi.
Figure 1. An overview of the experimental setup in wheat after Pst infection, integrating MeRIP-seq, RNA-seq, and qRT-PCR with phenotype analysis. (A) Schematic of the experimental workflow for MeRIP-seq, RNA-seq, and qRT-PCR in CK and Pst-infected wheat leaves. This diagram illustrates the step-by-step procedure adopted for analyzing both mock inoculation wheat leaves (CK) and Pst-infected wheat leaves (infected) at 7 dpi using NGS technology. The workflow highlights the incorporation of the MeRIP technique to selectively isolate m6A-modified RNA fractions for subsequent sequencing. (B) The phenotype of CK and infected wheat seedlings with Pst infection at 7 dpi. The CK wheat leaves showed no visible symptoms, but the infected wheat leaves showed chlorotic patches and blotches. (C) The magnification of CK and infected wheat leaves in the dotted box in (B). (D) sqRT-PCR assay of PsEF1 in CK and infected wheat leaves at 7 dpi.
Plants 13 00982 g001
Plants 13 00982 g002
Figure 2. Comprehensive analysis of RNA-seq data from CK and infected wheat samples at 7 dpi. (A) Principal component analysis (PCA) of CK groups (green dots) and infected groups (yellow triangles) with two biological replicates each in RNA-seq. (B) Heatmap of Pearson correlation coefficients across four samples. The gradient of colors, transitioning from white to dark blue, represents the correlation strength between each pair of samples, with darker shades indicating stronger correlations. (C) Numbers of differentially expressed genes (DEGs) in wheat after Pst infection with |Log2 (FC)| > 1 and FDR (false discovery rate) < 0.05. (D) Heatmap of DEGs with 20 clusters. Rows, individual mRNA transcripts; columns, individual CK and infected samples. Red and blue represent up-regulation and down-regulation of mRNA levels in CK and infected samples, respectively. (E) Volcano plot showing up-regulated genes (red) and down-regulated genes (blue) in wheat with Pst infection. The gene ID and corresponding gene descriptions for the three most significantly up-regulated genes as well as the three most significantly down-regulated genes are provided.
Figure 2. Comprehensive analysis of RNA-seq data from CK and infected wheat samples at 7 dpi. (A) Principal component analysis (PCA) of CK groups (green dots) and infected groups (yellow triangles) with two biological replicates each in RNA-seq. (B) Heatmap of Pearson correlation coefficients across four samples. The gradient of colors, transitioning from white to dark blue, represents the correlation strength between each pair of samples, with darker shades indicating stronger correlations. (C) Numbers of differentially expressed genes (DEGs) in wheat after Pst infection with |Log2 (FC)| > 1 and FDR (false discovery rate) < 0.05. (D) Heatmap of DEGs with 20 clusters. Rows, individual mRNA transcripts; columns, individual CK and infected samples. Red and blue represent up-regulation and down-regulation of mRNA levels in CK and infected samples, respectively. (E) Volcano plot showing up-regulated genes (red) and down-regulated genes (blue) in wheat with Pst infection. The gene ID and corresponding gene descriptions for the three most significantly up-regulated genes as well as the three most significantly down-regulated genes are provided.
Plants 13 00982 g002
Plants 13 00982 g003
Figure 3. GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of differentially expressed genes (DEGs) identified based on the RNA-seq data from CK and Pst-infected wheat leaves. (A) GO enrichment circle diagram of differentially expressed genes (DEGs). From outer circle to inner circle: Circle 1: The top 20 enriched GO terms of DEGs, and the coordinate scale outside the circle is the number of DEGs. Yellow and blue represent the GO terms of biological process and molecular function, respectively. Circle 2: The background of DEGs enriched in each GO term. The greater the number of DEGs, the longer the bar, and the smaller the q-value, the redder the color. Circle 3: The bar of the proportion of up-regulated (dark purple) and down-regulated (light purple) DEGs. The specific values are shown below. Circle 4: The rich factor value of each GO term (the number of DEGs divided by the total number of genes in the GO term); each grid line of the background grid represents 0.1. (B) Bubble chart illustrating the top 20 enriched GO terms in the cellular component category of DEGs based on q-value, with GO term ID and annotation on the y-axis and rich factor on the x-axis. (C) The histogram of level 2 GO terms enrichment classification of DEGs in infected samples in comparison to CK samples. The histogram depicts the number of up-regulated (in red) and down-regulated (in green) genes in each GO term. (D) KEGG pathway enrichment circle diagram of DEGs. From outer circle to inner circle: Circle 1: The top 20 enriched A class KEGG pathways of DEGs, and the coordinate scale outside the circle is the number of DEGs. Yellow represents the KEGG pathway of metabolism. Circle 2: The background of DEGs enriched in each pathway. The greater the number of DEGs, the longer the bar, and the smaller the q-value, the redder the color. Circle 3: The bar of the proportion of up-regulated (dark purple) and down-regulated (light purple) DEGs enriched in each pathway. The specific values are shown below. Circle 4: The rich factor value of each pathway (the number of DEGs divided by the total number of genes in the pathway); each grid line of the background grid represents 0.1. (E) Bubble chart illustrating the top 20 enriched KEGG pathways of DEGs based on q-value, with KEGG term annotation on the y-axis and rich factor on the x-axis.
Figure 3. GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of differentially expressed genes (DEGs) identified based on the RNA-seq data from CK and Pst-infected wheat leaves. (A) GO enrichment circle diagram of differentially expressed genes (DEGs). From outer circle to inner circle: Circle 1: The top 20 enriched GO terms of DEGs, and the coordinate scale outside the circle is the number of DEGs. Yellow and blue represent the GO terms of biological process and molecular function, respectively. Circle 2: The background of DEGs enriched in each GO term. The greater the number of DEGs, the longer the bar, and the smaller the q-value, the redder the color. Circle 3: The bar of the proportion of up-regulated (dark purple) and down-regulated (light purple) DEGs. The specific values are shown below. Circle 4: The rich factor value of each GO term (the number of DEGs divided by the total number of genes in the GO term); each grid line of the background grid represents 0.1. (B) Bubble chart illustrating the top 20 enriched GO terms in the cellular component category of DEGs based on q-value, with GO term ID and annotation on the y-axis and rich factor on the x-axis. (C) The histogram of level 2 GO terms enrichment classification of DEGs in infected samples in comparison to CK samples. The histogram depicts the number of up-regulated (in red) and down-regulated (in green) genes in each GO term. (D) KEGG pathway enrichment circle diagram of DEGs. From outer circle to inner circle: Circle 1: The top 20 enriched A class KEGG pathways of DEGs, and the coordinate scale outside the circle is the number of DEGs. Yellow represents the KEGG pathway of metabolism. Circle 2: The background of DEGs enriched in each pathway. The greater the number of DEGs, the longer the bar, and the smaller the q-value, the redder the color. Circle 3: The bar of the proportion of up-regulated (dark purple) and down-regulated (light purple) DEGs enriched in each pathway. The specific values are shown below. Circle 4: The rich factor value of each pathway (the number of DEGs divided by the total number of genes in the pathway); each grid line of the background grid represents 0.1. (E) Bubble chart illustrating the top 20 enriched KEGG pathways of DEGs based on q-value, with KEGG term annotation on the y-axis and rich factor on the x-axis.
Plants 13 00982 g003
Plants 13 00982 g004
Figure 4. Comprehensive analysis of MeRIP-seq data from CK and infected wheat samples at 7 dpi. (A) Principal component analysis (PCA) of CK groups (green dots) and infected groups (yellow triangles) with two biological replicates each in MeRIP-seq. (B) Heatmap of Pearson correlation coefficients across four samples in MeRIP-seq. The color from blue to red represents the correlation from weak to strong. (C) Distribution and density of m6A peaks of CK sample and infected sample within 5′UTR (untranslated region), CDS (coding sequence), and 3′UTR. (D) Motif analysis of enriched RRACH and DRACH (R is A/G, H is A/U/C and D is A/G/U) conserved motifs for m6A peaks in wheat from CK and Pst-infected samples. (E) Numbers of differentially m6A-modified peaks (DMPs) in wheat after Pst infection with |log2(FC)| > 1 and p-value < 0.05. (F) Volcano plot showing up-regulated (red) and down-regulated (blue) m6A peak-related genes in wheat after Pst infection. The gene ID and corresponding gene descriptions for the three most significantly up-regulated peak-related genes as well as the three most significantly down-regulated peak-related genes are provided.
Figure 4. Comprehensive analysis of MeRIP-seq data from CK and infected wheat samples at 7 dpi. (A) Principal component analysis (PCA) of CK groups (green dots) and infected groups (yellow triangles) with two biological replicates each in MeRIP-seq. (B) Heatmap of Pearson correlation coefficients across four samples in MeRIP-seq. The color from blue to red represents the correlation from weak to strong. (C) Distribution and density of m6A peaks of CK sample and infected sample within 5′UTR (untranslated region), CDS (coding sequence), and 3′UTR. (D) Motif analysis of enriched RRACH and DRACH (R is A/G, H is A/U/C and D is A/G/U) conserved motifs for m6A peaks in wheat from CK and Pst-infected samples. (E) Numbers of differentially m6A-modified peaks (DMPs) in wheat after Pst infection with |log2(FC)| > 1 and p-value < 0.05. (F) Volcano plot showing up-regulated (red) and down-regulated (blue) m6A peak-related genes in wheat after Pst infection. The gene ID and corresponding gene descriptions for the three most significantly up-regulated peak-related genes as well as the three most significantly down-regulated peak-related genes are provided.
Plants 13 00982 g004
Plants 13 00982 g005
Figure 5. GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of differentially m6A-modified peak (DMP)-related genes identified based on the MeRIP-seq data from Pst-infected wheat leaves in comparison to CK samples. (A) GO enrichment circle diagram of differentially m6A-modified peak (DMP)-related genes. From outer circle to inner circle: Circle 1: The top 20 enriched GO terms of DMP-related genes, and the coordinate scale outside the circle is the number of DMP-related genes. Yellow represents the GO term of biological process. Circle 2: The background of DMP-related genes enriched in each GO term. The greater the number of DMP-related genes, the longer the bar, and the smaller the q-value, the redder the color. Circle 3: The number of DMP-related genes enriched in each GO term. Circle 4: The rich factor value of each GO term (the number of DMP-related genes divided by the total number of genes in the GO term); each grid line of the background grid represents 0.1. (B) Bubble chart illustrating the top 20 enriched GO terms in the cellular component category of DMP-related genes based on q-value, with GO term ID and annotation on the y-axis and rich factor on the x-axis. (C) The histogram of level 2 GO terms enrichment classification of DMP-related genes in infected samples in comparison to CK samples. (D) KEGG pathway enrichment circle diagram of DMP-related genes. From outer circle to inner circle: Circle 1: The top 20 enriched A class KEGG pathways of DMP-related genes, and the coordinate scale outside the circle is the number of DMP-related genes. Yellow and pink represent the KEGG pathways of metabolism and genetic information processing, respectively. Circle 2: The background of DMP-related genes enriched in each pathway. The greater the number of DMP-related genes, the longer the bar, and the smaller the q-value, the redder the color. Circle 3: The number of DMP-related genes enriched in each pathway. Circle 4: The rich factor value of each pathway (the number of DMP-related genes divided by the total number of genes in the pathway); each grid line of the background grid represents 0.1. (E) Bubble chart illustrating the top 20 enriched KEGG pathways of DMP-related genes based on q-value, with KEGG term annotation on the y-axis and rich factor on the x-axis.
Figure 5. GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of differentially m6A-modified peak (DMP)-related genes identified based on the MeRIP-seq data from Pst-infected wheat leaves in comparison to CK samples. (A) GO enrichment circle diagram of differentially m6A-modified peak (DMP)-related genes. From outer circle to inner circle: Circle 1: The top 20 enriched GO terms of DMP-related genes, and the coordinate scale outside the circle is the number of DMP-related genes. Yellow represents the GO term of biological process. Circle 2: The background of DMP-related genes enriched in each GO term. The greater the number of DMP-related genes, the longer the bar, and the smaller the q-value, the redder the color. Circle 3: The number of DMP-related genes enriched in each GO term. Circle 4: The rich factor value of each GO term (the number of DMP-related genes divided by the total number of genes in the GO term); each grid line of the background grid represents 0.1. (B) Bubble chart illustrating the top 20 enriched GO terms in the cellular component category of DMP-related genes based on q-value, with GO term ID and annotation on the y-axis and rich factor on the x-axis. (C) The histogram of level 2 GO terms enrichment classification of DMP-related genes in infected samples in comparison to CK samples. (D) KEGG pathway enrichment circle diagram of DMP-related genes. From outer circle to inner circle: Circle 1: The top 20 enriched A class KEGG pathways of DMP-related genes, and the coordinate scale outside the circle is the number of DMP-related genes. Yellow and pink represent the KEGG pathways of metabolism and genetic information processing, respectively. Circle 2: The background of DMP-related genes enriched in each pathway. The greater the number of DMP-related genes, the longer the bar, and the smaller the q-value, the redder the color. Circle 3: The number of DMP-related genes enriched in each pathway. Circle 4: The rich factor value of each pathway (the number of DMP-related genes divided by the total number of genes in the pathway); each grid line of the background grid represents 0.1. (E) Bubble chart illustrating the top 20 enriched KEGG pathways of DMP-related genes based on q-value, with KEGG term annotation on the y-axis and rich factor on the x-axis.
Plants 13 00982 g005
Plants 13 00982 g006
Figure 6. Conjoint analysis of MeRIP-seq and RNA-seq data. (A) Dual Venn diagrams of m6A-modified peaks and their corresponding genes in CK and infected samples. Left: numbers of unique peaks in CK and infected samples, along with common peaks. Right: numbers of these peaks represent genes of the two groups. (B) A four-quadrant diagram illustrating significant changes in DEGs and DMP-related genes, comparing infected samples to CK. |log2(FC)| > 1 in both MeRIP-seq and RNA-seq; FC is fold change. (C) Bar chart and heat map of the hyper-up, hyper-down, hypo-up, and hypo-down genes are shown in Table 3.
Figure 6. Conjoint analysis of MeRIP-seq and RNA-seq data. (A) Dual Venn diagrams of m6A-modified peaks and their corresponding genes in CK and infected samples. Left: numbers of unique peaks in CK and infected samples, along with common peaks. Right: numbers of these peaks represent genes of the two groups. (B) A four-quadrant diagram illustrating significant changes in DEGs and DMP-related genes, comparing infected samples to CK. |log2(FC)| > 1 in both MeRIP-seq and RNA-seq; FC is fold change. (C) Bar chart and heat map of the hyper-up, hyper-down, hypo-up, and hypo-down genes are shown in Table 3.
Plants 13 00982 g006
Plants 13 00982 g007
Figure 7. Expression patterns of m6A writers, erasers, and readers in wheat after Pst inoculation. (A) Heatmap illustrating the expression profiles of m6A writers, erasers, and readers from MeRIP-seq and RNA-seq data in CK and infected samples. The color scale within the dendrogram represents the log2(TPM + 1) of each m6A regulator gene. Gene names in bold were further analyzed in (D–L). (B) Agarose gel (2.5%) electrophoresis displaying the sqRT-PCR products of PsEF1 in Pst-infected wheat leaves collected from 0 (CK), 0.25, 1, 4, 7, 10 dpi with three biological replicates. (C) Quantitative analysis of band intensity for sqRT-PCR products of PsEF1, as depicted in (B). Intensity values of all samples are normalized against the intensity value of the 0.25 dpi sample. (D–L) Relative expression levels of m6A regulators in Pst-infected wheat determined by RT-PCR. Each panel corresponds to a different m6A regulator gene: (D) TaVIR-D, (E) TaFIP37-2D, (F) TaVIR-A, (G) TaHAKAI1-A, (H) TaALKBH11B, (I) TaALKBH4B, (J) TaECT25, (K) TaECT31, and (L) TaECT21. Expression levels are normalized to those detected in the CK (0 dpi). The small case letters on the bars of the figures refer to the relative expression levels.
Figure 7. Expression patterns of m6A writers, erasers, and readers in wheat after Pst inoculation. (A) Heatmap illustrating the expression profiles of m6A writers, erasers, and readers from MeRIP-seq and RNA-seq data in CK and infected samples. The color scale within the dendrogram represents the log2(TPM + 1) of each m6A regulator gene. Gene names in bold were further analyzed in (D–L). (B) Agarose gel (2.5%) electrophoresis displaying the sqRT-PCR products of PsEF1 in Pst-infected wheat leaves collected from 0 (CK), 0.25, 1, 4, 7, 10 dpi with three biological replicates. (C) Quantitative analysis of band intensity for sqRT-PCR products of PsEF1, as depicted in (B). Intensity values of all samples are normalized against the intensity value of the 0.25 dpi sample. (D–L) Relative expression levels of m6A regulators in Pst-infected wheat determined by RT-PCR. Each panel corresponds to a different m6A regulator gene: (D) TaVIR-D, (E) TaFIP37-2D, (F) TaVIR-A, (G) TaHAKAI1-A, (H) TaALKBH11B, (I) TaALKBH4B, (J) TaECT25, (K) TaECT31, and (L) TaECT21. Expression levels are normalized to those detected in the CK (0 dpi). The small case letters on the bars of the figures refer to the relative expression levels.
Plants 13 00982 g007
Table 1. List of 24 annotations for 35 genes in the GO term related to “immune system process” using RNA-Seq data.
Table 1. List of 24 annotations for 35 genes in the GO term related to “immune system process” using RNA-Seq data.
Annotations/Description/Gene SymbolGene IDRegulation
UDP-glucose: glycoprotein glucosyltransferase (UGGT)TraesCS1A03G0409400
TraesCS1D03G0377700		Up
Lectin receptor kinase (LECRK91)TraesCS1B03G1249300
TraesCS1D03G1014400		Up
Receptor-like protein kinase, putative, expressed (CRK6)TraesCS2B03G0583500Up
GRAM domain family protein (VAD1)TraesCS2B03G0620400Up
Cyclic nucleotide-gated channel (CNGC4)TraesCS3A03G0774300Up
Calreticulin (CRT3)TraesCS3A03G0936700Up
Pathogenesis-related protein PR-4 (PR4A, PR4B)TraesCS3A03G1217600
TraesCS3B03G1460900
TraesCS3D03G1160600		Up
Glutamate decarboxylase (SPL)TraesCS3B03G0229400
TraesCS3D03G0168800		Up
MACPF domain protein (CAD1)TraesCS3D03G0628000Up
Pheophorbide a oxygenase, chloroplastic (PAO)TraesCS4A03G1017200
TraesCS4B03G0809700		Up
Receptor protein kinase, putative (RLK7)TraesCS5A03G0031900Up
Cysteine protease (CCP2)TraesCS5B03G0553900Up
Aminotransferase-like protein (OAT)TraesCS5D03G0861400Up
Pleiotropic drug resistance ABC transporter (ABCG48)TraesCS5D03G0992800Up
E3 ubiquitin-protein ligase (ATL6)TraesCS6A03G0290200Up
Alcohol dehydrogenase, putative (CAD2)TraesCS6A03G0429300
TraesCS6B03G0520600
TraesCS6D03G0363400		Up
Allene oxide cyclase (AOC)TraesCS6B03G1030800
TraesCS6D03G0731500		Up
BZIP transcription factor protein (BZIP50)TraesCS7A03G0964400Up
MACPF domain-containing protein (MACPF)TraesCS7B03G0282300Up
WRKY transcription factor (WRKY28)TraesCS7D03G1174900Up
GTPase Der (FZL)TraesCS1A03G0630100Down
Apolipoprotein D (CHL)TraesCS2D03G1058400Down
NAD(P)H-quinone oxidoreductase subunit N (ndhN)TraesCS3A03G0952700
TraesCS3B03G1089400
TraesCS3D03G0886900		Down
Telomere-associated protein (RIF1)TraesCS7A03G0315600Down
Table 2. List of 18 annotations for 80 genes in the KEGG “plant-pathogen interaction pathway” using RNA-Seq data.
Table 2. List of 18 annotations for 80 genes in the KEGG “plant-pathogen interaction pathway” using RNA-Seq data.
Annotations/Description/Gene SymbolGene IDRegulation
Glycerol kinase (GLPK)TraesCS2B03G1367900, TraesCS2D03G1154600Up
Cysteine proteinase (CCP1)TraesCS2A03G0676600, TraesCS2D03G0630100, TraesCS5A03G0791200, TraesCS5B03G0823900, TraesCS5D03G0747200Up
Calmodulin, putative (CML27, CML37, CML4)TraesCS1B03G1255000, TraesCS1B03G1255400, TraesCS5A03G0918000Up
Heat shock protein 90 (HSP81-1)TraesCS7A03G0555500, TraesCS7D03G0537900Up
Cyclic nucleotide-gated channel (CNGC2)TraesCS5D03G0900900Up
Chaperone protein htpG family protein (HSP90-5, HSP90)TraesCS5B03G0646100, TraesCS5D03G0595900, TraesCS7A03G1288100, TraesCS7B03G1204400, TraesCS7D03G1221800Up
Calcium-dependent protein kinase (CPK15, CPK9, CPK10, CPK8, CPK4)TraesCS1D03G0995500, TraesCS4A03G0728000, TraesCS5A03G1009700, TraesCS5B03G1055900, TraesCS5D03G0955400, TraesCS5B03G1159600, TraesCS5D03G1051900, TraesCS6A03G0203200, TraesCS6B03G0271200Up
WRKY transcription factor (WRKY24)TraesCS1A03G0172500, TraesCS1A03G0743700, TraesCS1B03G0235500, TraesCS1B03G0848500, TraesCS1D03G0166900, TraesCS1D03G0711800Up
Kinase family protein (PTI13)TraesCS2D03G0794500Up
Respiratory burst oxidase (RBOHA, RBOHC, RBOHE)TraesCS3A03G0705100, TraesCS4B03G0848700, TraesCS5A03G1171300, TraesCS5D03G0694900, TraesCS5B03G0560500, TraesCS5D03G0519000, TraesCS6D03G0366700Up
Calcium-binding protein (CML19, CML27, CML36, CML5, CML21, CML8, CML25)TraesCS3B03G1378700, TraesCS4B03G0485600, TraesCS5A03G0382400LC, TraesCS5B03G0150400, TraesCS5B03G0370400, TraesCS5D03G0789100, TraesCS7B03G1109900, TraesCS7B03G1175200, TraesCS7D03G1164400Up
Pathogenesis-related protein 1-1 (PRMS)TraesCS5A03G0484700, TraesCS5A03G1037300, TraesCS5A03G1037400, TraesCS5B03G0483900, TraesCS5B03G1087600, TraesCS5D03G0980700, TraesCS5D03G0980800, TraesCS7D03G0362800, TraesCS7D03G0450000Up
RPM1-interacting protein 4 (RIN4)TraesCS7A03G0593800Up
NBS-LRR disease resistance protein, putative, expressed (RPM1, RPP13L4, PIK-2, PIK6-NP)TraesCS1A03G0059400, TraesCS1D03G0032300, TraesCS2D03G0030700, TraesCS2D03G1070500, TraesCS3D03G0966800, TraesCSU03G0241100, TraesCS1D03G0036500, TraesCS7B03G0033000, TraesCS7D03G0011300Up
Disease resistance protein (NBS-LRR class) (RPS2)TraesCS2A03G0930000, TraesCS2B03G1021500, TraesCS2D03G0868500Up
3-ketoacyl-CoA synthase (FDH, KCS20, CUT1, KCS11)TraesCS4A03G0012800, TraesCS4A03G0046100, TraesCS4B03G0728500, TraesCSU03G0063000LC, TraesCS7D03G0827800Up
Mitogen-activated protein kinase (MPK5)TraesCS4A03G0218300Up
Elongation factor Tu (TUFA)TraesCS6B03G0694200, TraesCS6D03G0480200Down
Table 3. List of 48 genes that exhibit a significant differentiation in both m6A-modified peaks and transcript abundance in Pst-infected wheat leaves in comparison to CK.
Table 3. List of 48 genes that exhibit a significant differentiation in both m6A-modified peaks and transcript abundance in Pst-infected wheat leaves in comparison to CK.
Gene ID, Strand aDescription bPattern cPeak
Annotation d		Peak
Start e		Peak
End f		IP_
L2FC g		In_
L2FC h
TraesCS6B03G0561000, +Cytochrome P450 family proteinHyper-up5′, start, CDS2840106392840108145.011.64
TraesCS1D03G0829400, −Chaperone DnaJHyper-upstop, 3′4407010514407011014.644.02
TraesCS2B03G0935100, −Malate synthaseHyper-upstart, CDS5261030005261031254.484.57
TraesCS5B03G0999600, +ChitinaseHyper-upstart, CDS5836007145836007644.503.12
TraesCS7A03G1159700, +Ferredoxin--NADP reductase 2Hyper-upstart, CDS6765264766765265764.232.61
TraesCS3A03G0941700, −Tropinone reductase-like proteinHyper-up5′, start, CDS6477147366477148614.873.22
TraesCS4D03G0669200, −DUF1997 family proteinHyper-upstart, CDS4536709304536709804.483.28
TraesCS5A03G0935200, −Cysteine protease, putativeHyper-up5′, start, CDS5885195705885197336.035.66
TraesCS3B03G0133200, −D-Ala-D/L-Ala epimeraseHyper-upstart, CDS39875470398755455.023.99
TraesCS5B03G0752100, −Serine aminopeptidase S33Hyper-up5′, start, CDS4798158614798160114.644.80
TraesCS1D03G0951700, +GDP-mannose-3’,5’-epimeraseHyper-up5′4745017544745018294.193.12
TraesCS2D03G0725900, +Kinesin-like proteinHyper-upstop, 3′4078281854078283736.573.53
TraesCS3D03G0726600, +3-hydroxyacyl-CoA dehydrogenaseHyper-upstop, 3′4349588444349589194.511.83
TraesCS4A03G0035700, −plant/proteinHyper-upstart, CDS13580848135809234.771.69
TraesCS5D03G0513000, +Jasmonate zim-domain proteinHyper-upstart, CDS3315983153315984154.732.48
TraesCS2A03G0387200, +vacuolar sorting-associated proteinHyper-upCDS1532360441532361194.142.20
TraesCS3B03G0090100, +Nuclease S1Hyper-up5′, start, CDS24977420249776204.132.69
TraesCS1B03G0217600, +AAA+ ATPaseHyper-upstart, CDS75657827756579774.381.63
TraesCS2B03G1234100, −Amine oxidase family proteinHyper-up5′6955093496955095244.681.63
TraesCS2D03G0901800, −SPX domain-containing proteinHyper-upstop, 3′5099637905099639154.272.57
TraesCS2D03G0880500, −Aldehyde dehydrogenaseHyper-upstop, 3′4939751954939753204.192.57
TraesCS1A03G0652200, +Pyruvate phosphate dikinaseHyper-upstop, 3′4465854674465855674.434.06
TraesCS7D03G0606900, +Peptidylprolyl isomeraseHyper-downstart, CDS2541929412541929913.17−1.75
TraesCS3B03G1139000, +Alpha/beta hydrolaseHyper-downCDS7206052937206054232.39−1.57
TraesCS3B03G1515300LC, +P-loop NTPaseHypo-upCDS, stop847160574847160674−1.531.80
TraesCS7D03G0910600, +ABC transporterHypo-upstart, CDS502724870502724945−3.192.70
TraesCS3A03G0325200, +GDSL esterase/lipaseHypo-upCDS, stop131460529131460623−3.151.80
TraesCS3A03G0354000, +Cytochrome P450, putativeHypo-upCDS155007091155007221−1.501.98
TraesCS1D03G0777100, +Cysteine protease family proteinHypo-up5′, start421851652421851802−2.473.29
TraesCS1A03G0652400, −Pyruvate phosphate dikinaseHypo-upstart, CDS447084657447084809−1.591.71
TraesCS6B03G0574100, −CASP-like proteinHypo-up5′, start, CDS296054675296054825−2.691.97
TraesCS5A03G0517600, −Homeobox protein, putativeHypo-up5′404317960404318085−2.372.24
TraesCS1B03G0741000, +Pyruvate phosphate dikinaseHypo-up5′, start471557711471557786−1.852.48
TraesCS1B03G0574200, −Glutathione S-transferaseHypo-up5′, start, CDS355708499355708549−2.561.52
TraesCS1D03G0954900, +Chlorophyll a/b binding proteinHypo-down5′, start, CDS475689826475689876−2.17−3.10
TraesCS2D03G0959100, −Cold-responsive proteinWCOR15Hypo-down5′, start, CDS540101549540101599−3.47−2.90
TraesCS5D03G1027900, −Chlorophyll a/b binding proteinHypo-down5′, start, CDS511061211511061286−1.67−3.89
TraesCS1D03G0729600, −Chlorophyll a/b binding proteinHypo-down5′, start, CDS405463983405464083−2.37−3.47
TraesCS1B03G1053300, +Chlorophyll a/b binding proteinHypo-downstart, CDS629357113629357163−3.85−3.69
TraesCS6A03G0217100, +Chlorophyll a/b binding proteinHypo-downstart, CDS6426736264267412−2.64−3.37
TraesCS5B03G1135900, −Chlorophyll a/b binding proteinHypo-down5′, start, CDS638614728638614778−2.47−3.18
TraesCS1B03G1053200, +Chlorophyll a/b binding proteinHypo-downstart, CDS629262631629262681−3.97−3.86
TraesCS5D03G1027800, −Chlorophyll a/b binding proteinHypo-down5′, start, CDS511058147511058297−3.14−3.25
TraesCS1A03G0988100, +Chlorophyll a/b binding proteinHypo-down5′, start, CDS569868615569868715−4.36−3.34
TraesCS1D03G0873500, −Chlorophyll a/b binding proteinHypo-downCDS454772861454772911−4.66−3.19
TraesCS1A03G0761900, −Chlorophyll a/b binding proteinHypo-down5′, start, CDS500100318500100462−5.01−3.65
TraesCS1B03G0871200, −Chlorophyll a/b binding proteinHypo-down5′, start549194399549194449−3.22−3.48
TraesCS2D03G0286600LC, −RING/U-box superfamily proteinHypo-downstart, CDS8235692082356989−2.24−3.05
a The numbers and letters in bold in wheat gene names represent the chromosome where the gene is located. Wheat genes come from the sense strand (+) or antisense strand (−) of DNA; b gene descriptions come from WheatOmics (http://wheatomics.sdau.edu.cn/ (accessed on 13 February 2024)) or UniProtKB (https://www.uniprot.org/uniprotkb (accessed on 13 February 2024)); c The pattern indicates whether the m6A-modified peak of the differential gene is hyper-methylated or hypo-methylated and whether the transcript is up-regulated or down-regulated. The background colors are pink for hyper-up, purple for hyper-down, yellow for hypo-up, and green for hypo-down. d Distribution of m6A-modified peaks on different gene functional elements. “5′ “is 5′ untranslated region (UTR), “start” is start codon, “CDS” is coding sequence, “stop” is stop codon, and “3′ “is 3′ UTR. e The starting position of the m6A-modified peak in the chromosome; f the termination position of the m6A-modified peak in the chromosome. g In MeRIP-seq (IP) analysis, the fold change (FC) of m6A-modified peaks in infected samples relative to the CK is calculated and then subjected to log2 transformation (L2); h in RNA-seq (Input, In) analysis, the fold change (FC) of transcripts in infected samples relative to the CK is calculated and then subjected to log2 transformation (L2). The colors distinguish the genes showing differential patterns of regulation: pink for hyper-up, purple for hyper-down, yellow for hypo-up, and green for hypo-down, the same color scheme was used in Figure 6C.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Cerav, E.N.; Wu, N.; Akkaya, M.S. Transcriptome-Wide N6-Methyladenosine (m6A) Methylation Analyses in a Compatible Wheat–Puccinia striiformis f. sp. tritici Interaction. Plants 2024, 13, 982. https://doi.org/10.3390/plants13070982

AMA Style

Cerav EN, Wu N, Akkaya MS. Transcriptome-Wide N6-Methyladenosine (m6A) Methylation Analyses in a Compatible Wheat–Puccinia striiformis f. sp. tritici Interaction. Plants. 2024; 13(7):982. https://doi.org/10.3390/plants13070982

Chicago/Turabian Style

Cerav, Elif Naz, Nan Wu, and Mahinur S. Akkaya. 2024. "Transcriptome-Wide N6-Methyladenosine (m6A) Methylation Analyses in a Compatible Wheat–Puccinia striiformis f. sp. tritici Interaction" Plants 13, no. 7: 982. https://doi.org/10.3390/plants13070982

APA Style

Cerav, E. N., Wu, N., & Akkaya, M. S. (2024). Transcriptome-Wide N6-Methyladenosine (m6A) Methylation Analyses in a Compatible Wheat–Puccinia striiformis f. sp. tritici Interaction. Plants, 13(7), 982. https://doi.org/10.3390/plants13070982

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop