Lijuan Wu1,2
Junqiang Wang
Xinkun Hu
94% of researchers rate our articles as excellent or good
Learn more about the work of our research integrity team to safeguard the quality of each article we publish.
Find out more
- 1Institute of Ecology, China West Normal University, Nanchong, Sichuan, China
- 2College of Agronomy, Sichuan Agricultural University, Chengdu, Sichuan, China
- 3College of Life Science, China West Normal University, Nanchong, Sichuan, China
Introduction: Fusarium head blight (FHB), caused primarily by Fusarium graminearum (Fg), poses a significant threat to wheat production. It is necessary to deeply understand the molecular mechanisms underlying FHB resistance in wheat breeding.
Methods: In this study, the transcriptomic responses of two Chinese wheat landraces—Wuyangmai (WY, resistant) and Chinese Spring (CS, susceptible)—to F. graminearum infection were examined using RNA sequencing (RNA-seq). Differential expression of mRNAs, long non-coding RNAs (lncRNAs), circular RNAs (circRNAs), and microRNAs (miRNAs) was analyzed at 3 and 5 days post-Fg inoculation (dpi).
Results: The results showed that WY exhibited a targeted miRNA response, primarily modulating defense-related pathways such as glutathione metabolism and phenylpropanoid biosynthesis, which are crucial for oxidative stress regulation and pathogen defense response. In contrast, CS displayed a broader transcriptional response, largely linked to general metabolic processes rather than immune activation. Notably, the up-regulation of genes involved in oxidative stress and immune defense in WY confirmed its enhanced resistance to FHB. The integrated analysis of miRNA-mRNA interactions highlighted miRNAs as central regulators of defense mechanisms in WY, particularly at later stages of infection. These miRNAs targeted genes involved in immune responses, while lncRNAs and circRNAs played a more limited role in the regulation of defense responses. The GO and KEGG pathway enrichment analyses further revealed that WY enriched for plant-pathogen interaction and secondary metabolite biosynthesis pathways, which are crucial for pathogen resistance. In contrast, CS prioritized metabolic homeostasis, suggesting a less effective defense strategy.
Discussion: Overall, this study underscores the critical role of miRNA-mediated regulation in FHB resistance in WY. These insights into miRNA-mediated regulatory mechanisms provide a molecular basis for breeding FHB-resistant wheat varieties and highlight miRNA-mRNA interactions as promising targets for enhancing disease resilience.
1 Introduction
As one of the vital food crops, wheat (Triticum aestivum L.) supplies approximately 20% of the caloric intake for the global population (Ma et al., 2020). Wheat Fusarium head blight (FHB), destructive fungal disease primarily caused by F. graminearum (Fg), leads to a significant yield losses and quality deteriorates during epidemic years (Zhu et al., 2019). Additionally, FHB contaminates grain with harmful mycotoxins, including nivalenol and deoxynivalenol (DON), which seriously endanger the health of humans and livestock (Chen et al., 2019; Hu et al., 2023). Recent climate changes and wheat farming practices have increased the frequency and severity of FHB outbreaks (Ma et al., 2020). Identification and utilization of resistant germplasms in breeding programs are one of the most sustainable and economical approach to manage FHB.
FHB resistance is a complex quantitative trait governed by multiple genes and influenced by genotype-environment interactions (Zhang et al., 2021). Approximately 500 quantitative trait loci (QTLs) related to FHB resistance were identified from wheat and its relatives, and distributed across all 21 wheat chromosomes (Buerstmayr et al., 2009, Buerstmayr et al., 2020). However, most of these QTLs have minor effects and limited breeding value in addition to seven widely recognized major resistance genes (Fhb1 to Fhb7). Recently, two additional major QTLs, Fhb8 and Fhb9, were identified and associated with specific resistance types (Wang et al., 2024a; Zhang et al., 2024).
Despite these advances, only Fhb1 and Fhb7 have been cloned using map-based techniques and have shown significant resistance. Fhb1 confers resistance through a histidine-rich calcium-binding protein (His or TaHRC), while the exact mechanism remains unclear (Su et al., 2019; Li et al., 2019a). Fhb7, derived from Thinopyrum species, encodes a glutathione-S-transferase (GST) that detoxifies trichothecene toxins (Wang et al., 2020a). However, wheat varieties incorporated with Fhb7 are not yet widely available, and those with Fhb1 displayed a moderate resistance, likely due to interactions with parent genotype (Li et al., 2019b). FHB resistance is generally inherited quantitatively and affected by multiple minor genes (Bai and Shaner, 2004). Other resistance genes, such as TaFROG and TaABCC3, have been identified and characterized for their roles in response to DON and F. graminearum (Perochon et al., 2015; Walter et al., 2015). To effectively combat FHB in wheat effectively, more effort is still needed to identify new resistant genes.
Transcriptome analysis is a powerful tool for studying FHB resistance mechanisms, offering insights into wheat-pathogen interactions. Previous studies utilized microarray techniques to analyze transcriptomic responses (Bernardo et al., 2007), while recent ones have applied RNA-seq to reveal new pathways and defense-related genes (Xiao et al., 2013; Brauer et al., 2019; Pan et al., 2018). Small RNAs, particularly miRNAs, have been found to regulate the response of wheat to F. graminearum, either by silencing fungal genes or modulating host defense pathways (MaChado et al., 2018; Jin et al., 2020; Biselli et al., 2018). However, little is known about the role of other non-coding RNAs, such as long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs), in response to FHB resistance.
Chinese wheat landraces are gene pools for FHB-resistant. Some of the well-known resistant resources including Wangshuibai, Sanyuehuang, Shuilizhan, and Wuyangmai (WY), originated in China (Wan et al., 1997). WY, in particular, is an FHB-resistant landrace from Yibin city in Sichuan, China. However, the loci against FHB in WY have not yet been genetically identified, and the mechanisms underlying FHB resistance remain unknown. In this study, RNA-seq analysis was performed on WY (an FHB-resistant genotype) and Chinese Spring (CS, an FHB-susceptible genotype) in response to F. graminearum at 3 and 5 dpi to identify lncRNAs, circRNAs, miRNAs, and mRNAs involved in host-pathogen interactions. Differentially expressed miRNAs and mRNAs were used as central elements to target other differentially expressed RNAs, including lncRNAs and host genes associated with circRNAs. Additionally, differential expression analysis, along with Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses, was conducted on significantly up-regulated genes in CS and WY at different time points. Finally, the regulatory relationships between differentially expressed miRNAs and mRNAs were characterized, leading to the construction of a miRNA-mRNA targeting network map. The results provide valuable information for gaining a deeper understanding the FHB resistance differences of the two Sichuan wheat landraces, WY and CS.
2 Materials and methods
2.1 Plant materials and cultivation conditions
The two Sichuan wheat landraces, WY and CS, which exhibit high resistance (HR) and high susceptibility (HS) to FHB, respectively (Wan et al., 1997), were used in this study. The two landraces were provided by the Triticeae research Institute of Sichuan Agricultural University of China, and showed similar plant heights, spikes, heading and flowering dates, minimizing the potential discrepancies due to variations in plant architecture and development periods. The experiment was conducted at the breeding field of China West Normal University (Nanchong, Sichuan Province, China; 30°48′N, 106°05′E) during the 2022-2023 growing season. A randomized block design was implemented, with each genotype sown in a plot containing ten rows. Each row was 1.5 meters long, with 20 seedlings spaced 0.3 meters apart. Field cultivation and fertilization practice were the same as local wheat cultivation standards, and no fungicide was applied.
2.2 Macroconidia preparation and F. graminearum inoculation
A highly virulent and 15-acetyldeoxynivalenol (15ADON) producing isolate of F. graminearum, F0609, supplied by the Jiangsu Academy of Agricultural Science (Nanjing, Jiangsu Province, China), was used to infect the wheat plants. To prepare the macroconidia, isolate F0609 was first cultured on potato dextrose agar (PDA) medium at room temperature for 10 days under fluorescent-UV lights. A small piece of mycelium was then transferred to a mung bean medium and shaken at 180 rpm at 28°C for 4 days. The mycelia were filtered with double-layered medical gauze, and the macroconidia concentration was measured by a blood cell counting plate under a microscope. The macroconidia suspension was then diluted to a concentration of 1 × 105 spores/mL.
At mid-anthesis, 10 µL of the F. graminearum macroconidial spore suspension (1 × 105 spores/mL) was point-inoculated into the two basal florets of fully developed spikelet between the lemma and palea using a micropipette (Hu et al., 2019). Each wheat genotype was inoculated with three biological replicates, each consisting of 18 spikes. The inoculated spikes were then covered with a transparent plastic bag to maintain humidity and sprayed with sterile water twice daily until harvest. Spikelets were harvested at 0 (sterile water mock inoculation control), 3, and 5 dpi, with three biological replicates for each time point. Those samples collected from CS and WY at 0, 3 and 5 dpi were designated as CS-0, CS-3, CS-5, WY-0, WY-3, and WY-5, respectively. A total of nine samples were collected for each genotype.
2.3 RNA extraction, library construction and sequencing
Total RNA of each sample was extracted with the TRIzol reagent (Invitrogen, CA, USA), and the residual DNA was removed using RNase-free DNase I (Invitrogen, CA, USA) following the manufacturer’s instructions. The integrity and quality of the RNA were preliminarily assessed with 1.5% formaldehyde denaturing agarose gel electrophoresis. The concentration and purity of RNA were measured by a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, MA, USA). The RNA integrity was further determined by using the RNA 6000 Nano Assay Kit on the Agilent Bioanalyzer 2100 System (Agilent, CA, USA).
For cDNA library construction, 1.5 μg of RNA per sample was used for rRNA removal with the Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, USA). Sequencing libraries were then prepared with the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, MA, USA) following the manufacturer’s guidelines, with index codes included for sample identification. For small RNA (sRNA) library construction, 2.5 μg of RNA per sample was used. Sequencing libraries were generated with the NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, MA, USA) based on the manufacturer’s instructions, and index codes were also added for sample identification (Zhang et al., 2020). The quality of all libraries was assessed using the Agilent Bioanalyzer 2100. Paired-end sequencing (PE150-bp) was performed on the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) by Biomarker Technologies (Qingdao, China).
2.4 Data processing
The raw data (raw reads) were filtered with in-house Perl scripts to remove low-quality reads and sequences containing adapters or poly-N reads. Clean data (clean reads) were collected, and analyzed for the Q20, Q30, GC content, and sequence duplication levels. For sRNA-seq data, sequences ranging from 15 to 35 nucleotides (nt) in length were retained after trimming. All downstream analyses were conducted with high-quality clean data. The clean reads were then mapped to the CS reference genome IWGSC_RefSeq_v2.1 (IWGSC; Alaux et al., 2023) using HISAT2 software version 2.2.1 (Kim et al., 2019). Only reads with a perfect match or a single mismatch were further analyzed and annotated based on the reference genome.
2.5 Identification of mRNA, lncRNA, circRNA, and miRNA
For mRNA and lncRNA identification, the transcriptome was assembled using StringTie (v2.2.0) (Kovaka et al., 2019) based on the reads mapped to the wheat CS reference genome (IWGSC_RefSeq_v2.1). The GffCompare program (v0.12.6) (Pertea and Pertea, 2020) was used to annotate the assembled transcripts. Unknown transcripts were screened as potential lncRNAs. Four computational tools—CPC, CNCI, Pfam, and CPAT—were combined to distinguish non-protein-coding RNA (ncRNA) candidates from the unknown transcripts. Putative protein-coding RNAs were filtered out using thresholds for minimum transcript length and exon number. LncRNA candidates were selected based on transcript lengths greater than 200 nucleotides (nt) and more than 2 exons, and further validated with the four computational tools mentioned above. Those tools effectively differentiate protein-coding genes from non-coding ones and classify various lncRNA types, such as lincRNAs, intronic lncRNAs, antisense, and sense lncRNAs (Zhang et al., 2020).
The circRNA identification tool “CIRI” (CircRNA Identifier) (Gao et al., 2015) was used to detect circRNAs from the transcriptome data. To collect sufficient information for circRNA identification and characterization, the SAM files were scanned twice with CIRI. The identified circRNAs were then output with annotation information. Target miRNAs of the circRNAs and the target genes of miRNAs were predicted using miRanda (Doron et al., 2008) and RNAhybrid for animals (Rehmsmeier et al., 2004), and TargetFinder for plants (Bo and Wang, 2005). During the predictions, FASTA sequences of these circRNAs and miRNAs were used as input files.
The clean reads were aligned against four bioinformatic databases, Silva, GtRNAdb, Rfam, and Repbase with Bowtie software (Langmead and Salzberg, 2012) to identify candidate miRNA. The alignments allowed for length variations at both the 5′ and 3′ ends, as well as one internal mismatch. Common RNA families, including rRNA, tRNA, snRNA, snoRNA, and other ncRNAs, along with repeats, were filtered out. The remaining reads were then subject to miRNA identification through comparison with the wheat reference genome. The known and novel miRNAs were detected with the miRBase v22 database (Kozomara et al., 2019) and miRDeep2 module (Friedländer et al., 2012), respectively.
2.6 Differential expression analysis
The expression levels of coding genes and lncRNAs in each sample were estimated with fragments per kilobase of transcript per million fragments mapped (FPKM) values (Trapnell et al., 2010). The expression of circRNAs was determined based on the number of junction reads identified by the CIRI tool. For miRNAs, expression levels of each sample were estimated through the following steps: (1) sRNAs were mapped back to the precursor sequence, and (2) the read count for each miRNA was obtained from the mapping results. The expression levels of miRNAs and circRNAs in each sample were calculated using transcripts per million (TPM) (Fahlgren et al., 2007). Differential expression analysis among three time points (0, 3, and 5 dpi) of two genotypes (WY and CS) was conducted with DESeq2 R package (v1.10.1) (Love et al., 2014). Differentially expressed genes (DEGs) were identified using an adjusted p-value (calculated using the Benjamini and Hochberg method) ≤ 0.01 to control the false discovery rate (FDR) and a |log2 (Fold change (FC))| > 1. Venny 2.1 (Oliveros, 2015) was then used for Venn diagram analysis, and SRplot web server (Tang et al., 2023; available at http://www.bioinformatics.com.cn/SRplot) was used for heatmap drawing.
2.7 Analysis of competing endogenous RNA
Candidate ceRNA relationship pairs were identified based on miRNA targeting relationships. Using the predicted miRNA-mRNA, lncRNA-miRNA, and circRNA-miRNA interaction pairs, groups of lncRNA-miRNA-mRNA or circRNA-miRNA-mRNA that shared the same miRNAs were collected. The screening criteria for ceRNAs were as follows: (1) The number of shared miRNAs among ceRNAs must be greater than or equal to 5; (2) The hypergeometric test p-value must be < 0.01, and the corrected FDR value must be < 0.01; (3) mRNAs and lncRNA/circRNA interaction pairs with a co-expression correlation (r > 0.9) were further selected for ceRNA network construction.
2.8 Gene functional annotation and GO and KEGG enrichment analysis
Gene function annotation was performed using multiple databases, including Nr (Non-redundant protein), Nt (Non-redundant nucleotide), Pfam (Protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and non-redundant protein sequence database), KEGG, and GO. GO enrichment analysis of DEGs was conducted with the GOseq R package, based on Wallenius’ non-central hypergeometric distribution (Young et al., 2010), which adjusts for gene length bias in DEGs. The enrichment analysis of DEGs in KEGG pathways was performed using KOBAS software (version 2.0) (Xie et al., 2011).
2.9 Real-time quantitative PCR validation
For RT-qPCR analysis, cDNA of all RNA samples (those used for Illumina RNA-sequencing) was synthesized with the PrimeScript RT reagent Kit (TaKaRa, Dalian, China) according to the manufacturer’s protocol. The reaction volume was 20 μL including random hexamers and 1 μg of each RNA sample. Two wheat genes—glyceraldehyde-3-phosphate dehydrogenase (GAPDH, TraesCS7A02G313100) and indole-3 acetaldehyde oxidase (IAAOx, TraesCS2A02G246300) (Hu et al., 2019)—were used as reference genes to normalize gene expression in the two wheat genotypes. Primers (Dataset S1) were designed using the PrimerQuest™ Tool from Integrated DNA Technologies (IDT) and synthesized by Sangon Biotech (Shanghai, China). For each reaction, 1 μL of 10×-diluted cDNA was added to a 10 μL reaction volume using the SYBR Green PCR kit (TaKaRa, Dalian, China). RT-qPCR was carried out in an ABI7500 Real-Time PCR System (Applied Biosystems, CA, USA). The cycling conditions were 5 minutes at 95°C, followed by 40 cycles of 30 seconds at 95°C, 30 seconds at the melting temperature (as indicated in Dataset S1), and 30 seconds at 72°C. The melting curve analysis was conducted from 55°C to 95°C, with readings taken every 1°C and held for 5 seconds. Four technical replicates were performed for each sample. The relative expression levels were calculated using the 2−ΔΔCt method (Livak and Schmittgen, 2001).
3 Results
3.1 Sequencing and identification of RNAs
After sequencing of the 18 samples, the raw data were filtered to obtain clean data, resulting in a total of 1,035,915,752 reads and 309.95 Gb of clean data. Each sample yielded approximately 15.95 Gb of clean data, with Q20 and Q30 base percentages exceeding 97.02% and 92.05%, respectively. The GC content ranged from 47.39% to 50.51% (Dataset S2). These metrics indicate that the quality of RNA sequencing in this study was high and suitable for further analysis. The clean reads were mapped to the reference wheat genome sequences (IWGSC_RefSeq_v2.1) of Chinese Spring (Alaux et al., 2023). After identification of RNAs, a total of 15,765 lncRNAs, 286,397 genes, 621 circRNAs, and 3,063 miRNAs were detected.
3.2 Differential expression analysis of mRNA, lncRNA, circRNA, and miRNA
To identify differentially expressed RNAs, pairwise comparisons were conducted between time points for the two wheat genotypes. In CS, the numbers of differentially expressed mRNAs, lncRNAs, circRNAs, and miRNAs for 0 vs. 3-dpi, 0 vs. 5-dpi, and 3 vs. 5-dpi, were 5,713, 174, 8, and 29; 20,596, 331, 7, and 172; and 16,028, 229, 10, and 69, respectively. In WY, the numbers of differentially expressed mRNAs, lncRNAs, circRNAs, and miRNAs for 0 vs. 3-dpi, 0 vs. 5-dpi, and 3 vs. 5-dpi were 2,196, 150, 5, and 30; 9,793, 288, 11, and 179; and 6,769, 188, 5, and 239, respectively (Figure 1). The expression profiles of those differentially expressed RNAs were visualized using circos maps, with the height of the profiles representing significance (-log10 (FDR)). The results showed that the F. graminearum invasion had a significant impact on the expression of mRNA and miRNA, and had moderate to low effects on the expression of lncRNAs and circRNAs, respectively (Figure 2) (Supplementary Figure 1).
Figure 1
Figure 1. Statistics of differentially expressed RNAs across various comparisons. The 0, 3, and 5-dpi represent samples collected at 0, 3, and 5 days post F. graminearum inoculation (dpi), respectively. CS and WY refer to Chinese Spring and Wuyangmai, respectively. LncRNA, long non-coding RNA; miRNA, microRNA; circRNA, circular RNA.
Figure 2
Figure 2. Expression profiles of differentially expressed RNAs at 5 dpi with F. graminearum. The outermost ring shows chromosome information, followed by mRNA (gene), lncRNA, circRNA, and miRNA. For each group of differentially expressed RNAs, red, blue, and the height indicate up-regulation, down-regulation, and significance [-log10 (FDR)], respectively.
Analysis of the differentially expressed RNAs showed that the numbers of differentially expressed genes in CS were 2.60, 2.10, and 2.37 times greater than those in WY for the comparison time points of 0 vs. 3-dpi, 0 vs. 5-dpi, and 3 vs. 5-dpi, respectively. The differentially expressed lncRNAs in CS were slightly higher than those in WY, with fold changes of 1.16, 1.15, and 1.22 for the above-mentioned time points. However, the numbers of differentially expressed circRNAs were very low in both CS and WY, ranging from 5 to 11. Interestingly, the number of differentially expressed miRNAs in WY was slightly higher than that in CS for the comparison time points of 0 vs. 3 dpi and 0 vs. 5 dpi, and 3.46 times higher for the comparison time points of 3 vs. 5 dpi (Table 1).
Table 1
Table 1. Comparative analysis of differentially expressed RNAs between Chinese Spring (CS) and Wuyangmai (WY).
3.3 Integrated analysis of differentially expressed RNA targeting relationships
Differentially expressed lncRNAs, circRNAs, miRNAs, and mRNAs were analyzed by taking each RNA type as the center and examining the targeted relationships involving other differentially expressed RNAs (including host genes corresponding to circRNAs). First, an integrated analysis was conducted using differentially expressed lncRNAs (DE_lncRNAs) as the center. In the comparisons of 0 vs. 3-dpi, 0 vs. 5-dpi, and 3 vs. 5-dpi of CS, the Venn diagram revealed that 0 (1), 11 (14), and 2 (4) lncRNAs, respectively, were located at the intersection of the three data sets: DE_lncRNA, DE_Cis.mRNA_Target lncRNA (or DE_Trans.mRNA_Target lncRNA), and DE_miRNA_Target lncRNA. Similarly, for the pairwise comparisons of the three time points in WY, the Venn diagram indicated that 0 (0), 1 (7), and 1 (4) lncRNAs were located at the intersection of the same three data sets, respectively (Supplementary Figures 2A, B). When circRNAs were taken as the center for the targeting relationship analysis, no circRNAs were found at the intersection of the three data sets: DE_circRNA, DE_Hostgene_circRNA, and DE_miRNA_Target circRNA, for any of the pairwise comparisons among the three time points in CS or WY (Supplementary Figure 3).
Subsequently, an integrated analysis of targeting relationships was conducted using differentially expressed miRNAs (DE_miRNA) as the center. In the pairwise comparisons of the three time points in CS, the Venn diagram indicated that 0 (3), 4 (9), and 1 (4) miRNA were present at the intersection of the three data sets: DE_miRNA, DE_circRNA_Target miRNA (or DE_lncRNA_Target miRNA), and DE_mRNA_Target miRNA, respectively. Similarly, for the pairwise comparisons of the three time points in WY, the Venn diagram showed that 0 (0), 1 (28), and 0 (28) miRNAs were located at the intersections of the same data sets, respectively (Figure 3). Additionally, the analysis revealed that the primary targeting relationships existed between the DE_miRNA and DE_mRNA_Target miRNA data sets. Specifically, 24 (18), 150 (182), and 61 (172) miRNAs were identified at their intersection in the pairwise comparisons of the three CS (or WY) time points, respectively (Figure 3).
Figure 3
Figure 3. Interaction of differentially expressed miRNAs with all miRNAs targeted by differentially expressed mRNAs, lncRNAs, and circRNAs. Here, DE_miRNA represents differentially expressed miRNAs; DE_mRNA_Target miRNA, DE_lncRNA_Target miRNA, and DE_circRNA_Target miRNA represent miRNAs targeted by differentially expressed genes, lncRNAs, and circRNAs, respectively.
Finally, differentially expressed mRNAs (DE_mRNA) were used as the center for analyzing targeting relationships. The Venn diagram showed that the primary targeting relationships occurred between the DE_mRNA data set and either DE_lncRNA_Target mRNA or DE_miRNA_Target mRNA. In the pairwise comparisons of the three CS time points, 39 (1,241), 283 (7,027), and 166 (5,608) mRNAs were detected at the intersection of DE_mRNA with DE_lncRNA_Target Cis.mRNA (or DE_lncRNA_Target Trans.mRNA), respectively. Similarly, for the pairwise comparisons of the three WY time points, 12 (167), 42 (1,551), and 36 (1,358) mRNAs were found at the corresponding intersections of the above two data sets (Supplementary Figures 4A, B). When considering the targeting relationships between DE_mRNA and DE_miRNA_Target mRNA, 72 (23), 2,152 (1,239), and 1,179 (783) mRNAs were identified at their intersections in the pairwise comparisons of the three CS (or WY) time points, respectively (Supplementary Figures 4A, B).
3.4 Identification of ceRNA
mRNAs, lncRNAs, and circRNAs can function as competing endogenous RNAs (ceRNAs) in a ceRNA network, where they regulate expression levels of each other by competitively binding to the same miRNA response elements (MREs) (Salmena et al., 2011). Based on the screening criteria for ceRNAs, co-expression analysis identified 98 circRNA-miRNA-mRNA relationship pairs and 49 lncRNA-miRNA-mRNA pairs (Dataset S3). However, no differentially expressed ceRNA relationship pairs were identified using a one-step nearest-neighbor network analysis, when each group of differentially expressed RNAs was extracted from the ceRNA relationship pairs.
3.5 Differential expression analysis for up-regulated genes
Differential gene expression analysis revealed that in CS, there were 3,989, 13,141, and 10,521 up-regulated DEGs, and 1,724, 7,455, and 5,507 down-regulated DEGs for 0 vs. 3-dpi, 0 vs. 5-dpi, and 3 vs. 5-dpi comparisons, respectively. The number of up-regulated DEGs was consistently much greater than that of down-regulated DEGs across all three comparisons. In WY, there were 1,060, 5,099, and 4,608 up-regulated DEGs, and 1,136, 4,694, and 2,161 down-regulated DEGs for the same comparisons. The number of up-regulated DEGs was higher than down-regulated DEGs only in the 3 vs. 5-dpi comparison, while for the 0 vs. 3-dpi and 0 vs. 5-dpi comparisons, the up-regulated DEGs were slightly fewer or greater than the down-regulated ones, respectively (Figure 4A). Additionally, the expression levels and significance of up-regulated DEGs in CS were consistently much higher than those in WY across all three comparisons (Figure 4B) (Supplementary Figures 5A, 6A).
Figure 4
Figure 4. Differential expression analysis for up-regulated genes. (A) Statistical analysis for all genes. (B) Volcano plot of DEGs obtained in the comparison of 3 vs. 5-dpi. (C) GO term and (D) KEGG pathway enrichment analysis for significantly up-regulated DEGs. (E) Gene expression heatmap of glutathione metabolism pathway.
GO enrichment analysis of up-regulated DEGs revealed distinct defense responses between the FHB-resistant wheat CS and the susceptible wheat WY. For the 0 vs. 3-dpi comparison, two of the top five GO terms differed between CS and WY. In CS, the top terms were “phenylalanine ammonia-lyase activity” and “heme binding”, while in WY, they were “hydrolase activity, hydrolyzing O-glycosyl compounds” and “polysaccharide binding” (Supplementary Figure 5B). In the 0 vs. 5-dpi comparison, CS showed enrichment for “glutathione transferase activity” and “iron ion binding”, while WY had “ATPase-coupled transmembrane transporter activity” and “ATPase activity” (Supplementary Figure 6B). In the 3 vs. 5-dpi comparison, three of the top five GO terms were different: CS showed “protein kinase activity”, “protein serine/threonine kinase activity”, and “phenylalanine ammonia-lyase activity”, while WY showed “glutathione transferase activity”, “heme binding”, and “UDP-glycosyltransferase activity” (Figure 4C).
KEGG pathway enrichment of up-regulated DEGs also highlighted distinct defense responses between CS and WY. In the 0 vs. 3-dpi comparison, only one of the top five KEGG pathways was different: “biosynthesis of amino acids” (ko01230) in CS, compared to “plant hormone signal transduction” (ko04075) in WY, although the enrichment of the top five pathways in WY was not significant. In the 0 vs. 5-dpi comparison, two of the top five pathways differed: “phenylpropanoid biosynthesis” (ko00940) and “carbon metabolism” (ko02010) in CS, and “glutathione metabolism” (ko00480) and “ABC transporters” (ko02010) in WY. For the 3 vs. 5-dpi comparison, four of the top five pathways were different between CS and WY: “plant-pathogen interaction” (ko04626), “biosynthesis of amino acids” (ko01230), “starch and sucrose metabolism” (ko00500), and “carbon metabolism” (ko01200) in CS, compared to “glutathione metabolism” (ko00480), “ABC transporters” (ko02010), “flavonoid biosynthesis” (ko00941), and “phenylalanine metabolism” (ko00360) in WY (Figure 4D) (Supplementary Figures 5C, 6C).
Additionally, a log2 (FC) > 2 threshold was applied to filter up-regulated genes enriched in the “glutathione metabolism” pathway in the 3 vs. 5-dpi comparison of WY, resulting in the identification of 93 DEGs. The expression heatmap of these up-regulated DEGs showed that 80 DEGs were shared between CS and WY. Among these, two DEGs (TraesCS2B03G1525600 and TraesCS6B03G0425000) in CS had expression levels more than two times higher than in WY. Furthermore, the heatmap also showed that 13 DEGs were differentially expressed only in WY (Figure 4E).
3.6 Integrated analysis of the targeting relationship between miRNA and mRNA
Using small RNA and transcriptome sequencing data, DE_miRNAs and DEGs were identified in the two sample groups. The regulatory relationships between miRNAs and mRNAs were then explored, focusing on the negative regulatory effects of miRNAs on mRNAs. First, DE_miRNAs were used as screening criteria to identify the mRNAs regulated by these miRNAs, followed by an analysis of the pairs of DE_miRNAs and mRNAs with negative regulatory relationships (Dataset S4). Similarly, DEGs were used as screening criteria to identify miRNAs that regulate these DEGs, and DEG-miRNA pairs with negative regulatory relationships were analyzed (Dataset S5). The targeting relationships between miRNAs and mRNAs were visualized using Cytoscape software (v3.10.2). The top 20 DE_miRNAs (or top 18 for 0 vs. 3-dpi comparison in WY) most relevant to target gene regulation were selected to construct the miRNA-mRNA targeting relationship network map (Dataset S6).
The results showed that only two miRNAs were shared between CS and WY in the 0 vs. 3-dpi comparison, whereas 16 and 15 miRNAs were shared in the 0 vs. 5-dpi and 3 vs. 5-dpi comparisons, respectively (Figure 5A) (Supplementary Figure 7A). In the 0 vs. 3-dpi comparison, ten DE_miRNAs in CS and eight in WY specifically targeted DEGs with negative regulatory relationships, with each DE_miRNA targeting one to four DEGs (Supplementary Figure 7B). In the 0 vs. 5-dpi and 3 vs. 5-dpi comparisons, all DE_miRNAs specifically expressed in CS and WY targeted DEGs with negative regulatory relationships, with each DE_miRNA targeting 14 to 30 DEGs. Notably, novel_miR_228, which was up-regulated in CS, and tae-miR1122a, which was down-regulated in WY, were present in both comparisons. Novel_miR_228 down-regulated 28 and 23 DEGs in CS, while tae-miR1122a up-regulated 27 and 23 DEGs in WY, respectively (Figure 5B) (Supplementary Figure 7B).
Figure 5
Figure 5. Targeting relationships between DE_miRNAs and DEGs. (A) Comparison of the top 20 DE_miRNAs most relevant to target gene regulation in CS and WY. (B) Targeting relationships between the DE_miRNAs specifically expressed in CS and WY and DEGs. (C) Expression heatmap of genes significantly down-regulated by novel_miR_228 and significantly up-regulated by tae-miR1122a with molecular function (MF) GO term. Heatmap represents log2 (FC) between F. graminearum infected and control.
GO analysis of these DEGs revealed that seven DEGs down-regulated by novel_miR_228 in CS had molecular function (MF) terms in the 0 vs. 5-dpi comparison, and nine DEGs in the 3 vs. 5-dpi comparison. Significant GO term enrichment analysis (using a q-value < 0.05 as the threshold) showed two DEGs, TraesCS1A03G0367900 (histone acetyltransferase activity; q = 1.52e-02) and TraesCS6B03G0592100 (double-stranded DNA binding; q = 1.67e-43), were significantly enriched in the 0 vs. 5-dpi comparison of CS, while six DEGs—TraesCS2A03G0587000 (protein kinase activity; q = 9.31e-06), TraesCS2D03G0569400 (protein kinase activity; q = 9.31e-06), TraesCS3D03G0550100 (hydrolase activity; q = 1.61e-02), TraesCS5A03G0080300 (microtubule binding; q = 1.12e-02), TraesCS7A03G0224800 (single-stranded DNA binding; q = 4.39e-04), and TraesCS7A03G0450300 (GTP binding; q = 4.47e-08)—were significantly enriched in the 3 vs. 5-dpi comparison of CS. The expression of these DEGs in CS was lower than in WY (Figure 5C) (Supplementary Figure 7C).
Likewise, GO analysis showed that 14 DEGs up-regulated by tae-miR1122a in WY had MF terms in both the 0 vs. 5-dpi and 3 vs. 5-dpi comparisons. Nine DEGs up-regulated by tae-miR1122a in WY were not differentially expressed in CS in the 0 vs. 5-dpi comparison, and three (TraesCS2B03G0537600, TraesCS2D03G0425900, and TraesCS6A03G0855700) were significantly enriched in the MF term for metal ion binding (q = 2.14e-03). In the 3 vs. 5-dpi comparison, nine DEGs up-regulated by tae-miR1122a in WY were not differentially expressed in CS, and four (TraesCS2B03G0537600, TraesCS2D03G0425900, TraesCS6A03G0855700, and TraesCS7D03G0982200) were significantly enriched in the MF term for metal ion binding (q = 2.94e-02). The expression of these DEGs in WY was higher than in CS (Figure 5C) (Supplementary Figure 7C).
3.7 KOG functional classification of DEGs regulated by miRNA
After analyzing the targeting relationships between miRNAs and mRNAs, functional classification analysis was performed on the DEGs identified. The KOG functional classification of DE_miRNA-regulated DEGs was conducted. Focusing on the top five functional classes, four classes—T (signal transduction mechanisms), R (general function prediction only), O (posttranslational modification, protein turnover, chaperones), and J (translation, ribosomal structure, and biogenesis)—were shared by CS and WY in the comparison of 0 vs. 3-dpi. The fifth class differed, with class S (function unknown) identified in CS and class A (RNA processing and modification) identified in WY. However, in the 0 vs. 5-dpi and 3 vs. 5-dpi comparisons, the top five functional classes for both CS and WY were identical: T, R, O, G (carbohydrate transport and metabolism), and K (transcription) (Supplementary Figure 8). Similarly, the KOG functional classification of miRNA-regulated DEGs was performed. The results showed that the top five functional classes—T, R, O, G, and Q (secondary metabolites biosynthesis, transport, and catabolism)—were consistent across all three comparisons for both CS and WY (Supplementary Figure 9).
3.8 GO analysis of DEGs regulated by miRNA
GO analysis of DE_miRNA-regulated DEGs revealed 22, 18, and 17 functional classifications in the biological process (BP), cellular component (CC), and molecular function (MF) categories, respectively, for both CS and WY. Among these, the percentage of genes in 16 BP, 15 CC, and 12 MF classifications exceeded 0.1%. The results showed that most DEGs in the BP category were associated with metabolic processes, cellular processes, and single-organism processes. In the CC category, the majority of DEGs were located in the membrane, cell, cell part, membrane part, and organelle. In the MF category, most DEGs were related to catalytic activity and binding functions (Supplementary Figure 10). Similar findings were obtained from the GO analysis of miRNA-regulated DEGs (Supplementary Figure 11).
GO enrichment analysis was conducted on the identified DEGs (Dataset S7, Dataset S8), with the top 20 terms for each comparison selected for display, using a q-value < 0.05 as the significance threshold. GO enrichment based on DE_miRNA-regulated DEGs revealed the most significant BP terms in CS were phenylalanyl-tRNA aminoacylation, RNA modification, and leucyl-tRNA aminoacylation for the 0 vs. 3-dpi, 0 vs. 5-dpi, and 3 vs. 5-dpi comparisons, respectively, whereas those in WY were phenylalanyl-tRNA aminoacylation, defense response, and mitochondrial RNA modification. The BP terms with the highest number of genes in CS included gene silencing by RNA, defense response, and protein folding, while those in WY were gene silencing by RNA, defense response, and nucleic acid phosphodiester bond hydrolysis. For the CC category, the most significant terms in CS were the THO complex (part of the transcription export complex), microtubule, and Ino80 complex, while those in WY were the ROC complex, Golgi membrane, and cell plate. The CC terms with the most genes in CS were cytoplasm (for the first two comparisons) and microtubule (for the third one), while in WY, the terms with the most genes were cytoplasm (for the first comparison) and Golgi membrane (for the subsequent two). In the MF category, the most significant terms in CS were phenylalanine-tRNA ligase activity, double-stranded DNA binding, and carbohydrate binding, while in WY they were calcium transmembrane transporter activity, phosphorylative mechanism, and double-stranded DNA binding. The terms with the most genes in CS were protein binding, catalytic activity, and protein kinase activity, whereas in WY, they were RNA binding, protein kinase activity, and RNA binding (Figure 6).
Figure 6
Figure 6. GO term enrichment analysis for DE_miRNA-regulated DEGs. GeneNum denotes the gene number.
GO enrichment analysis based on miRNA-regulated DEGs revealed that the most significant BP terms in CS were protein-chromophore linkage, trehalose biosynthetic process, and leucyl-tRNA aminoacylation in the 0 vs. 3-dpi, 0 vs. 5-dpi, and 3 vs. 5-dpi comparisons, respectively. In WY, the most significant BP terms were protein-chromophore linkage in both 0 vs. 3-dpi and 0 vs. 5-dpi, and trehalose biosynthetic process for 3 vs. 5-dpi. The BP terms with the most genes in CS included defense response, carbohydrate metabolic process, and recognition of pollen, while in WY, they were defense response, photosynthesis, and response to oxidative stress. In the CC category, significant terms were observed only in two comparisons of WY (0 vs. 3-dpi and 0 vs. 5-dpi), with the most significant CC term being photosystem I in both cases. The CC terms with the most genes were the chloroplast thylakoid membrane for both comparisons, though this term was not statistically significant in the second comparison. In the MF category, the most significant MF term in the comparison of 0 vs. 3-dpi in CS was calmodulin binding, with the term having the most genes being protein kinase activity. In the other two comparisons (0 vs. 5-dpi and 3 vs. 5-dpi), the most significant MF term and the term with the most genes in CS were both protein kinase activity. For WY, the most significant MF terms in the 0 vs. 3-dpi and 0 vs. 5-dpi comparisons were chlorophyll binding and protein serine/threonine kinase activity, respectively. In both cases, the MF term with the most genes was protein kinase activity. In the 3 vs. 5-dpi comparison in WY, both the most significant MF term and the MF term with the most genes were protein serine/threonine kinase activity (Supplementary Figure 12).
3.9 KEGG analysis of DEGs regulated by miRNA
KEGG analysis based on DE_miRNA-regulated DEGs indicated that most metabolic pathways were within the metabolism category, with the pathway containing the most DEGs being the plant-pathogen interaction pathway (Supplementary Figure 13). Similar results were obtained from the KEGG analysis based on miRNA-regulated DEGs (Supplementary Figure 14). KEGG enrichment analysis was conducted on these DEGs (Dataset S9, Dataset S10), and the top 20 pathways for each comparison were selected for display, using q-value < 0.05 as the significance threshold. KEGG enrichment of DE_miRNA-regulated DEGs showed that the most significantly enriched pathway with the highest number of DEGs was the plant-pathogen interaction pathway across all comparisons, except for 3 vs. 5-dpi in CS, where the starch and sucrose metabolism pathway was the most significantly enriched (Figure 7). Similarly, KEGG enrichment based on miRNA-regulated DEGs also identified the plant-pathogen interaction pathway as the most significantly enriched pathway (Supplementary Figure 15).
Figure 7
Figure 7. KEGG pathway enrichment analysis for DE_miRNA-regulated DEGs. GeneRatio indicates the gene ratio.
To compare KEGG pathway enrichment between CS and WY at different time points, the top five pathways from each comparison were selected for further analysis. First, the top five KEGG pathways enriched based on DE_miRNA-regulated DEGs were examined. Results showed that only two pathways were identical between CS and WY in the 0 vs. 3-dpi comparison: plant-pathogen interaction (ko04626) and RNA transport (ko03013). However, these pathways were not significant in WY. The three specific pathways in CS were ubiquitin-mediated proteolysis (ko04120), spliceosome (ko03040; not significant), and ABC transporters (ko02010). In contrast, the three specific pathways in WY were phenylpropanoid biosynthesis (ko00940), aminoacyl-tRNA biosynthesis (ko00970), and RNA polymerase (ko03020), and none of them were significant. Additionally, both CS and WY had the same top five pathways for the 0 vs. 5-dpi comparison. For the 3 vs. 5-dpi comparison, four of the top five pathways were shared by CS and WY. The only differing pathway in CS was other glycan degradation (ko00511), while in WY, it was plant-pathogen interaction (ko04626) (Table 2).
Table 2
Table 2. KEGG pathway enrichment analysis of DE_miRNA-regulated DEGs in CS and WY (top 5 pathways).
The KEGG pathways were enriched based on miRNA-regulated DEGs. In the 0 vs. 3-dpi comparison, four of the top five pathways were identical between CS and WY. However, the carbon metabolism pathway (ko02010) in CS and the starch and sucrose metabolism pathway (ko00500) in WY were not significant. The only differing pathway in CS was phenylpropanoid biosynthesis (ko00940), while in WY, it was photosynthesis (ko00195). For the other two comparisons (0 vs. 5-dpi, and 3 vs. 5-dpi), two of the top five pathways were different between CS and WY in each case, although the other three pathways—plant-pathogen interaction (ko04626), starch and sucrose metabolism (ko00500), and ABC transporters (ko02010)—were identical. In the 0 vs. 5-dpi comparison, the two different pathways in CS were flavonoid biosynthesis (ko00941) and glyoxylate and dicarboxylate metabolism (ko00630). In WY, the different pathways were carbon metabolism (ko01200) and photosynthesis (ko00195; not significant). In the 3 vs. 5-dpi comparison, the two different pathways in CS were flavonoid biosynthesis (ko00941) and cyanoamino acid metabolism (ko00460). In WY, the different pathways were phenylpropanoid biosynthesis (ko00940; not significant) and phenylalanine metabolism (ko00360) (Table 3).
Table 3
Table 3. KEGG pathway enrichment analysis of miRNA-regulated DEGs in CS and WY (top 5 pathways).
3.10 RT−qPCR validation of DEGs regulated by miRNA
Based on the integrated analysis of the targeting relationships between miRNA and mRNA, two DE_miRNAs—tae-miR1122a and novel_miR_228—were identified from the top 20 DE_miRNAs and were specifically expressed in WY and CS at 5-dpi. Six DEGs, three with a negative regulatory relationship with tae-miR1122a and the other three with novel_miR_228, were selected for RT-qPCR validation of their expression levels. The DEGs targeted by tae-miR1122a were TraesCS1B03G0240900, TraesCS2D03G0425900, and TraesCS3A03G0573700, associated with the KEGG pathway of glutathione metabolism (ko00480), phagosome pathway (ko04145), and autophagy-other (ko04136), respectively. The DEGs targeted by novel_miR_228 were TraesCS2A03G0587000, TraesCS2D03G0569400, and TraesCS4B03G1000000, and both of the first two were associated with the plant-pathogen interaction pathway (ko04626), and the third one had no KEGG pathway annotation. The relative expression data from RT-qPCR for the selected genes were compared with the RNA-seq analysis data, and the results revealed that the expression patterns were largely consistent between the two methods, confirming that the RNA-seq analysis in this study was reliable and suitable for further research (Figure 8).
Figure 8
Figure 8. Relative expression comparison for selected genes using RNA-seq and RT-qPCR data.
4 Discussion
4.1 Transcriptomic study for FHB resistance
Fusarium head blight (FHB) significantly threatens wheat production globally, leading to yield loss and grain contamination by mycotoxins. Breeding and utilizing FHB-resistant wheat cultivars provide effective and eco-friendly solutions to mitigate these effects. A comprehensive understanding of the molecular mechanisms underlying wheat-pathogen interactions has laid the genetic foundation for wheat breeding programs aimed at combating FHB. Recent transcriptomic studies advanced our knowledge of key molecular pathways associated with FHB resistance and susceptibility (Kazan and Gardiner, 2018), and will aid in identifying key genes and mechanisms involved in the plant’s defense response. For example, previous studies have shown that hormone biosynthesis and signal transduction pathways actively respond to FHB infection, and contribute to plant defense reaction (Biselli et al., 2018). Specifically, salicylic acid (SA) and jasmonic acid (JA) positively affect FHB resistance, while auxin and abscisic acid (ABA) are linked to susceptibility of FHB. Ethylene has dual roles in the interaction with F. graminearum (Wang et al., 2018). Additionally, resistant genotypes exhibit early and intense expression of defense-related genes, including those involved in redox homeostasis and secondary metabolite biosynthesis (Xiao et al., 2013). Differentially expressed miRNAs and lncRNAs play roles in regulating gene expression related to biotic and abiotic stress responses, respectively (Biselli et al., 2018; Soresi et al., 2023). Here, a comprehensive transcriptomic analysis of two Sichuan wheat landraces with contrasting FHB resistance (WY was resistant, while CS was susceptible) was conducted based on the prior research (Wan et al., 1997). Analysis of the expression and regulatory networks of lncRNAs, circRNAs, mRNAs, and miRNAs provides a detailed view of the regulatory differences between resistant and susceptible genotypes.
4.2 Differential expression of lncRNAs, circRNAs, and mRNAs
Our findings reveal significant differences in RNA expression profiles between WY and CS at different stages post-inoculation (Figures 1, 2). The differential expression of lncRNAs is particularly noteworthy. Previous studies have indicated that lncRNAs modulate the expression of defense-related genes, such as those involved in the JA pathway, which is crucial for plant disease resistance (Chen et al., 2023). Specific lncRNAs, like ALEX1 in rice, enhance resistance to bacterial blight by activating JA signaling and increasing JA content (Yu et al., 2019). Furthermore, lncRNAs can function as ceRNAs, decoying miRNAs to regulate the expression of target genes involved in immune responses. For instance, lncRNAs in tomato modulate MYB transcription factors by decoying miR159, enhancing resistance to Phytophthora infestans (Cui et al., 2020). LncRNAs have also been shown to regulate pathogen resistance and often peak early during pathogen invasion (Duan et al., 2020). Our results confirm this, with WY and CS displaying substantial lncRNA activity at 3 and 5-dpi, suggesting that early lncRNA activity is crucial for an effective defense against F. graminearum. Moreover, the greater number of differentially expressed lncRNAs observed in the field conditions compared to greenhouse settings suggests that environmental factors significantly influence lncRNA-mediated responses (Soresi et al., 2023).
In addition to lncRNAs, which are known for their regulatory roles, the circRNAs have emerged as another crucial class of regulatory RNAs. The role of circRNAs in responding to FHB is still not well-defined. Our data showed limited differential expression of circRNAs in both genotypes (Figures 1, 2), aligning with findings that indicate a reduction in circRNA expression following F. graminearum infection, possibly due to their involvement in fine-tuning gene expression during initial infection stages (Yin et al., 2022). Together, these results suggest that lncRNAs actively mediate gene regulation during early infection, while circRNAs may play a more subdued role, potentially involved in more sustained immune responses.
The differentially expressed mRNAs was more pronounced in CS than in WY (Table 1), consistent with the previous studies that susceptible genotypes often exhibit a broader transcriptional response. This reflects a less-targeted response, activating a wide array of genes rather than focusing on specific defense-related pathways (Erayman et al., 2015). In contrast, resistant genotypes like WY likely mount a more efficient and focused response, resulting in fewer DEGs but potentially more effective pathogen resistance mechanisms (Walker et al., 2024).
4.3 Integrated analysis of targeting relationships among RNAs
To further understand the regulatory roles of these RNAs, we analyzed their interactions and targeting relationships. Focusing on miRNAs as central regulators, we observed primary targeting relationships with mRNAs (Figure 3), particularly at 3 and 5 dpi. This miRNA-mRNA interaction network suggests that miRNAs play a vital role in silencing pathogen-related genes, curbing the spread of infection and enhancing plant resistance (Jiao and Peng, 2018; MaChado et al., 2018; Fu et al., 2023). In contrast, fewer interactions were observed when lncRNAs or circRNAs were analyzed as the central regulators (Supplementary Figures 2, Supplementary Figures 3). This lack of significant ceRNA relationships at later infection stages indicates that while lncRNAs and circRNAs may be important during initial pathogen recognition, their influence diminishes as the infection progresses. This temporal division of labor highlights miRNAs as the primary regulatory players in sustained immune responses, while lncRNAs and circRNAs contribute more to the early stages of infection.
4.4 Differential expression of up-regulated genes and their role in wheat FHB resistance
The differential expression analysis showed that CS had more up-regulated DEGs than WY across all comparisons (Figures 4A, B). In CS, the number of up-regulated genes exceeded the down-regulated ones at all time points, indicating a stronger overall transcriptional response to infection. However, this response was likely non-specific and may contribute to the higher susceptibility of CS to FHB. In contrast, WY displayed more balanced gene expression, with up-regulated DEGs outnumbering down-regulated ones only at 3 vs. 5 dpi. This suggests that while both genotypes mount a defense response to F. graminearum, a broader transcriptional activation of CS may lead to an inflammatory response that fails to efficiently control pathogen spread (Pan et al., 2018; Erayman et al., 2015).
The GO term enrichment analysis revealed differences in the molecular functions (MF) of the up-regulated DEGs between CS and WY (Figure 4C). At 0 vs. 3 dpi, CS showed enrichment for phenylalanine ammonia-lyase activity and heme binding, both of which are important for phenylpropanoid biosynthesis and managing oxidative stress. These pathways help protect the plant from damage and activate secondary metabolites involved in defense against fungal pathogens (Ramaroson et al., 2022). In contrast, WY was enriched for hydrolase activity and polysaccharide binding, suggesting a more localized, mechanical defense, likely related to altering the cell wall to prevent fungal invasion. At 0 vs. 5 dpi, CS enriched for glutathione transferase activity and iron ion binding, both involved in detoxifying reactive oxygen species (ROS) and maintaining redox balance. These processes are crucial during pathogen infection (Van Breusegem et al., 2008). WY, on the other hand, showed enrichment for ATPase-coupled transmembrane transporter activity and ATPase activity, which are involved in ion transport and energy metabolism, possibly indicating a different or less efficient stress response. At 3 vs. 5 dpi, CS showed enrichment for protein kinase activity and protein serine/threonine kinase activity, both crucial for signalling immune responses (Pan et al., 2018; Jiang et al., 2022). The presence of phenylalanine ammonia-lyase activity in CS also suggests ongoing activation of pathways important for producing lignin and antimicrobial compounds. WY, by contrast, was enriched for glutathione transferase activity, heme binding, and UDP-glycosyltransferase activity, which are involved in detoxification, ROS management, and modifying secondary metabolites (Wang et al., 2020a; Gharabli et al., 2023; Walker et al., 2024). These differences highlight the more specialized defense mechanisms in WY, which may explain its greater resistance to FHB.
The KEGG pathway enrichment analysis highlighted further differences in the defense mechanisms of CS and WY (Figure 4D). At 0 vs. 3 dpi, CS showed enrichment in the biosynthesis of amino acids pathway, which is important for producing amino acids needed for protein synthesis and secondary metabolite formation. In contrast, WY enriched for plant hormone signal transduction, suggesting that hormonal regulation (such as JA and SA) plays a key role in its early defense response (Pieterse et al., 2012; Pan et al., 2018). This supports the idea that the resistance of WY to FHB may involve a more coordinated immune response, with specific hormones activated by fungal infection. At 0 vs. 5 dpi, CS enriched for pathways related to phenylpropanoid biosynthesis and carbon metabolism, both crucial for stress responses, secondary metabolite production, and energy balance during infection. In contrast, WY enriched for glutathione metabolism and ABC transporters, which help detoxify oxidative stress and remove harmful substances (Wang et al., 2020b; Gullner et al., 2018). This suggests that the resistance of WY may depend more on detoxification and transporting stress-related compounds. At 3 vs. 5 dpi, the differences were even more pronounced. CS enriched for plant-pathogen interaction and amino acid biosynthesis, which are important for immune recognition and defense protein synthesis. WY, however, enriched for glutathione metabolism and ABC transporters, helping maintain cellular integrity during stress (Wang et al., 2020b; Gullner et al., 2018). Additionally, WY showed enrichment in flavonoid biosynthesis and phenylalanine metabolism, which suggests a shift towards producing secondary metabolites as a defense strategy to limit pathogen growth.
A detailed analysis of the up-regulated genes in the glutathione metabolism pathway during the 3 vs. 5 dpi comparison in WY identified 93 DEGs, of which 80 were shared between CS and WY (Figure 4E). These genes are involved in detoxifying ROS and maintaining cellular redox balance, which are essential for plant defense against oxidative stress caused by fungal infection. Interestingly, 13 DEGs were only differentially expressed in WY, suggesting a unique resistance mechanism in this genotype. The heatmap analysis showed that two genes, TraesCS2B03G1525600 and TraesCS6B03G0425000, had significantly higher expression levels in CS compared to WY, which may indicate a larger-scale response in CS, but this doesn’t necessarily mean a more effective defense. The differences in gene expression highlight the distinct defense strategies between the two genotypes, with WY showing a more localized, focused response to oxidative stress.
4.5 Integrated analysis of miRNA-mRNA targeting relationships
The results show a clear difference in miRNA profiles between the two wheat genotypes, especially at different time points after infection. At 0 vs. 3 dpi, only two miRNAs were common between CS and WY, indicating that each genotype activates different miRNA pathways early in the infection. However, by 0 vs. 5 dpi and 3 vs. 5 dpi, 16 and 15 miRNAs were shared between the genotypes, suggesting a more similar response as the infection progresses. Notably, novel_miR_228 and tae-miR1122a were key miRNAs with different roles in the two genotypes. Novel_miR_228, which was up-regulated in CS, down-regulated 28 and 23 genes related to defense in CS at 0 vs. 5 dpi and 3 vs. 5 dpi, respectively (Figure 5B) (Supplementary Figure 7B). This suggests that it may weaken the defense response of host, making CS more susceptible to FHB. On the other hand, tae-miR1122a, which was down-regulated in WY, up-regulated 27 and 23 defense-related genes in WY at the same time points, suggesting that it helps enhance resistance in WY. These findings support the idea that miRNAs play a role in regulating the defense mechanisms of host in response to fungal stress, contributing to FHB resistance (MaChado et al., 2018; Jin et al., 2020).
The GO analysis of the DEGs regulated by the identified miRNAs revealed important molecular functions linked to FHB resistance. In CS, novel_miR_228 down-regulated genes related to key functions like histone acetyltransferase activity and DNA binding, which are essential for transcription and DNA repair. This down-regulation could weaken the ability of host to manage gene expression and repair DNA, making it more vulnerable to fungal damage. Genes like TraesCS1A03G0367900 and TraesCS6B03G0592100 showed strong associations with these functions, highlighting how F. graminearum might exploit these weaknesses to infect the host. In contrast, in WY, tae-miR1122a up-regulated genes involved in metal ion binding, such as TraesCS2B03G0537600, TraesCS2D03G0425900, and TraesCS6A03G0855700 (Figure 5C) (Supplementary Figure 7C). These genes help maintain cellular stability under stress, suggesting that they play a role in the defense of plants (Wang et al., 2020b). The up-regulation of these genes in WY may support antioxidant and detoxification pathways, helping the plant cope with the oxidative damage caused by fungal infection (Gullner et al., 2018). The significant enrichment of these genes in WY points to their potential role in strengthening the ability of host to resist fungal stress by binding toxic ions released by the pathogen, thus helping maintain cellular integrity and defense against FHB.
The comparison between WY and CS shows a clear difference in how they respond to F. graminearum infection. In WY, the expression of defense-related genes is stronger, with down-regulation of miRNAs like tae-miR1122a and up-regulation of their target genes involved in metal ion binding and stress responses. This suggests that the resistance of WY to FHB may be due to better regulation of metal ion balance, oxidative stress, and immune signaling, all of which are essential for protecting cells during infection. In CS, the up-regulation of novel_miR_228 and the down-regulation of key defense genes indicate a weaker response to the fungal infection. This may explain why CS, a more susceptible genotype, struggles to mount a strong defense against F. graminearum, making it more prone to disease.
4.6 GO term and KEGG pathway enrichment of DEGs regulated by miRNA
GO and KEGG pathway analyses further illustrate the functional divergence between WY and CS in response to FHB. In GO classification, the DEGs of WY were enriched in defense-related categories, including gene silencing by RNA and defense response (Figure 6). This enrichment suggests that WY is genetically predisposed to actively respond to pathogen invasion, prioritizing immune functions over general metabolic processes. In the cellular component (CC) category, the DEGs of WY were primarily localized to membrane and cytoplasmic components, critical sites for pathogen recognition and signaling. In contrast, CS displayed a more generalized distribution across organelles, indicative of a less-targeted defense approach.
KEGG pathway enrichment revealed that immune pathways were central to the response of WY, particularly the plant-pathogen interaction pathway (Figure 7). This pathway was consistently enriched across all time points, reflecting a robust immune response in WY, channeling resources toward defense mechanisms upon pathogen detection (Zaynab et al., 2018). Conversely, the response of CS emphasized starch and sucrose metabolism pathways (Table 2), suggesting a baseline metabolic stability rather than a targeted immune response. This finding is consistent with previous studies that susceptible genotypes often allocate resources to general metabolism at the expense of pathogen-specific defenses (Brauer et al., 2019).
Additionally, pathways like phenylpropanoid biosynthesis and phenylalanine metabolism were enriched in WY but not in CS (Table 3). These pathways are known to produce secondary metabolites that strengthen structural defenses and stress responses, underscoring the preparedness of WY for pathogen attack (Zaynab et al., 2018). The divergence in pathway enrichment between WY and CS highlights how WY allocates resources more strategically, channeling efforts toward immune functions while CS appears to prioritize metabolic maintenance. This strategic resource allocation likely contributes to enhancing the resistance of WY to FHB.
4.7 Role of miRNAs in regulating defense mechanisms
This study highlighted the critical role of miRNAs in response to FHB resistance. miRNAs are known to target specific mRNAs, modulating gene expression in response to stressors, including pathogens. For example, miRNAs are integral to the defense mechanisms of plants against pathogens. They regulate the expression of genes involved in phytohormone signaling, ROS production, and nucleotide-binding site-leucine-rich repeat (NBS-LRR) gene expression, which are critical for mounting an effective defense response (Yang et al., 2021). The interaction between miRNAs and lncRNAs complicated this regulatory network, further fine-tuning the immune response of plants (Song et al., 2021). We found that WY exhibited a notably stronger miRNA response than CS, particularly at later infection stages (Table 1). miRNAs in WY were enriched in defense-related pathways such as glutathione metabolism and plant-pathogen interactions, crucial for oxidative stress regulation—an essential factor in plant defense that limits pathogen spread (Dorion et al., 2021). The enrichment of miRNAs in plant-pathogen interaction pathways supports their role in enhancing immune responses, allowing WY to recognize and respond more effectively to F. graminearum invasion. In contrast, miRNA pathways in CS were predominantly linked to general metabolic processes rather than defense-specific responses, suggesting that the miRNAs of CS may not be primed for targeted pathogen defense. This aligns with previous studies that susceptible genotypes often prioritize maintaining metabolic homeostasis over rapid immune activation, leaving them more vulnerable to pathogens (Brauer et al., 2019). These differences indicate that miRNAs in resistant genotypes like WY are more adept at regulating specific pathways linked to pathogen defense, facilitating a faster and more efficient immune response.
In summary, this study provides valuable insight into the molecular mechanisms underlying FHB resistance in two Sichuan wheat landraces WY (HR) and CS (HS). Through a comprehensive transcriptomic analysis, we identified distinct differences in how these two genotypes respond to F. graminearum. The resistant genotype, WY, showed a targeted miRNA response, particularly at later infection stages, which effectively regulated pathways important for defense, such as glutathione metabolism and phenylpropanoid biosynthesis. In contrast, the susceptible genotype, CS, exhibited a broader transcriptional response focused more on general metabolism than specific pathogen defense. Our findings emphasize the critical role of miRNA-mRNA interactions in enhancing the resistance of WY. Although lncRNAs and circRNAs were also involved in resistance, their contributions were limited in comparison with miRNAs. This research not only emphasized the importance of miRNA-mediated defense mechanisms but also suggested potential molecular targets to improve FHB resistance in breeding programs. Future studies should focus on validating these pathways and incorporating deep learning method and genomic selection (GS) (Chen et al., 2024; Wang et al., 2024b; Ma et al., 2024) techniques to develop more resilient wheat varieties, ensuring sustainable wheat production amid ongoing challenges posed by FHB disease.
Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: NCBI, PRJNA1183537.
Author contributions
LW: Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. JW: Investigation, Writing – review & editing. SS: Investigation, Writing – review & editing. ZY: Writing – review & editing. XH: Conceptualization, Resources, Writing – review & editing, Formal analysis, Funding acquisition, Methodology, Project administration, Supervision.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. Funding was provided by “the Key Research and Development Program of Science and Technology Department of Sichuan Province, grant number 2021YFS0342” and “Fundamental Research Funds of China West Normal University, grant number 18Q053 and 20A018”.
Acknowledgments
We would like to thank Dr. Youping Xu and Ms. Tingting Wang for preparing the inoculum, and thank Professor Zehong Yan for providing wheat materials and revising the article.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2025.1537605/full#supplementary-material
References
Alaux, M., Dyer, S., Sen, T. Z. (2023). “Wheat data integration and FAIRification: IWGSC, grainGenes, ensembl and other data repositories,” in The wheat genome. Eds. Appels, R., Eversole, K., Feuillet, C., Gallagher, D. (Compendium of Plant Genomes. Springer, Cham). doi: 10.1007/978-3-031-38294-9_2
Bai, G., Shaner, G. (2004). Management and resistance in wheat and barley to fusarium head blight. Annu. Rev. Phytopathol. 42, 135–161. doi: 10.1146/annurev.phyto.42.040803.140340
Bernardo, A., Bai, G., Guo, P., Xiao, K., Guenzi, A. C., Ayoubi, P. (2007). Fusarium graminearum-induced changes in gene expression between Fusarium head blight-resistant and susceptible wheat cultivars. Funct. Integr. Genomics 7, 69–77. doi: 10.1007/s10142-006-0028-1
Biselli, C., Bagnaresi, P., Faccioli, P., Hu, X., Balcerzak, M., Mattera, M. G., et al. (2018). Comparative transcriptome profiles of near-isogenic hexaploid wheat lines differing for effective alleles at the 2DL FHB resistance QTL. Front. Plant Sci. 9. doi: 10.3389/fpls.2018.00037
Bo, X., Wang, S. (2005). TargetFinder: a software for antisense oligonucleotide target site selection based on MAST and secondary structures of target mRNA. Bioinformatics 21, 1401–1402. doi: 10.1093/bioinformatics/bti211
Brauer, E. K., Rocheleau, H., Balcerzak, M., Pan, Y. L., Fauteux, F., Liu, Z. Y., et al. (2019). Transcriptional and hormonal profiling of Fusarium graminearum-infected wheat reveals an association between auxin and susceptibility. Physiol. Mol. Plant Pathol. 107, 33–39. doi: 10.1016/j.pmpp.2019.04.006
Buerstmayr, H., Ban, T., Anderson, J. A. (2009). QTL mapping and marker assisted selection for Fusarium head blight resistance in wheat: a review. Plant Breed. 128, 1–26. doi: 10.1111/j.1439-0523.2008.01550.x
Buerstmayr, M., Steiner, B., Buerstmayr, H. (2020). Breeding for Fusarium head blight resistance in wheat-progress and challenges. Plant Breed. 139, 429–454. doi: 10.1111/pbr.12797
Chen, D., Zhang, Z., Chen, Y., Li, B., Chen, T., Tian, S. (2023). Transcriptional landscape of pathogen-responsive lncRNAs in tomato unveils the role of hydrolase encoding genes in response to Botrytis cinerea invasion. Plant Cell Environ. 47, 651–663. doi: 10.1111/pce.14757
Chen, Y., Kistler, H. C., Ma, Z. (2019). Fusarium graminearum trichothecene mycotoxins: Biosynthesis, regulation, and management. Annu. Rev. Phytopathol. 57, 15–39. doi: 10.1146/annurev-phyto-082718-100318
Chen, Y., Wang, W., Yang, Z., Peng, H., Ni, Z., Sun, Q., et al. (2024). Innovative computational tools provide new insights into the polyploid wheat genome. Abiotech 5, 52–70. doi: 10.1007/s42994-023-00131-7
Cui, J., Jiang, N., Hou, X., Wu, S., Zhang, Q., Meng, J., et al. (2020). Genome-wide identification of lncRNAs and analysis of ceRNA networks during tomato resistance to Phytophthora infestans. Phytopathology 110, 456–464. doi: 10.1094/PHYTO-04-19-0137-R
Dorion, S., Ouellet, J. C., Rivoal, J. (2021). Glutathione metabolism in plants under stress: beyond reactive oxygen species detoxification. Metabolites 11, 641. doi: 10.3390/metabo11090641
Doron, B., Manda, W., Aaron, G., Marks, D. S., Chris, S. (2008). The microRNA.org resource: targets and expression. Nucleic Acids Res. 36, 149–153. doi: 10.1093/nar/gkm995
Duan, X., Song, X., Wang, J., Zhou, M. (2020). Genome-wide identification and characterization of Fusarium graminearum-responsive lncRNAs in Triticum aestivum. . Genes 11, 1135. doi: 10.3390/genes11101135
Erayman, M., Turktas, M., Akdogan, G., Gurkok, T., Inal, B., Ishakoglu, E., et al. (2015). Transcriptome analysis of wheat inoculated with Fusarium graminearum. Front. Plant Sci. 6. doi: 10.3389/fpls.2015.00867
Fahlgren, N., Howell, M. D., Kasschau, K. D., Chapman, E. J., Sullivan, C. M., Cumbie, J. S., et al. (2007). High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PloS One 2, e219. doi: 10.1371/journal.pone.0000219
Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688
Fu, K., Wu, Q., Jiang, N., Hu, S., Ye, H., Hu, Y., et al. (2023). Identification and expressional analysis of siRNAs responsive to Fusarium graminearum infection in wheat. Int. J. Mol. Sci. 24, 16005. doi: 10.3390/ijms242116005
Gao, Y., Wang, J., Zhao, F. (2015). CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 16, 4. doi: 10.1186/s13059-014-0571-3
Gharabli, H., Della Gala, V., Welner, D. H. (2023). The function of UDP-glycosyltransferases in plants and their possible use in crop protection. Biotechnol. Adv. 67, 108182. doi: 10.1016/j.bioteChadv.2023.108182
Gullner, G., Komives, T., Király, L., Schröder, P. (2018). Glutathione S-transferase enzymes in plant-pathogen interactions. Front. Plant Sci. 9. doi: 10.3389/fpls.2018.01836
Hu, W. J., Fu, L. P., Gao, D. R., Li, D. S., Liao, S., Lu, C. B. (2023). Marker-assisted selection to pyramid Fusarium head blight resistance loci Fhb1 and Fhb2 in the high-quality soft wheat cultivar Yangmai 15. J. Integr. Agr. 22, 360–370. doi: 10.1016/j.jia.2022.08.057
Hu, X., Rocheleau, H., McCartney, C., Biselli, C., Bagnaresi, P., Balcerzak, M., et al. (2019). Identification and mapping of expressed genes associated with the 2DL QTL for fusarium head blight resistance in the wheat line Wuhan 1. BMC Genet. 20, 47. doi: 10.1186/s12863-019-0748-6
Jiang, L., Zhang, S., Su, J., Peck, S. C., Luo, L. (2022). Protein kinase signaling pathways in plant-Colletotrichum interaction. Front. Plant Sci. 12. doi: 10.3389/fpls.2021.829645
Jiao, J., Peng, D. (2018). Wheat microRNA1023 suppresses invasion of Fusarium graminearum via targeting and silencing FGSG_03101. J. Plant Interact. 13, 514–521. doi: 10.1080/17429145.2018.1528512
Jin, X., Jia, L., Wang, Y., Li, B., Sun, D., Chen, X. (2020). Identification of Fusarium graminearum-responsive miRNAs and their targets in wheat by sRNA sequencing and degradome analysis. . Funct. Integr. Genomics 20, 51–61. doi: 10.1007/s10142-019-00699-8
Kazan, K., Gardiner, D. M. (2018). Transcriptomics of cereal-Fusarium graminearum interactions: what we have learned so far. Mol. Plant Pathol. 19, 764–778. doi: 10.1111/mpp.12561
Kim, D., Paggi, J. M., Park, C., Bennett, C., Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907–915. doi: 10.1038/s41587-019-0201-4
Kovaka, S., Zimin, A. V., Pertea, G. M., Razaghi, R., Salzberg, S. L., Pertea, M. (2019). Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol. 20, 278. doi: 10.1186/s13059-019-1910-1
Kozomara, A., Birgaoanu, M., Griffiths-Jones, S. (2019). miRBase: from microRNA sequences to function. Nucleic Acids Res. 47, D155–D162. doi: 10.1093/nar/gky1141
Langmead, B., Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods. 9, 357–359. doi: 10.1038/nmeth.1923
Li, G., Zhou, J., Jia, H., Gao, Z., Fan, M., Luo, Y., et al. (2019a). Mutation of a histidine-rich calcium-binding-protein gene in wheat confers resistance to Fusarium head blight. Nat. Genet. 51, 1106–1112. doi: 10.1038/s41588-019-0426-7
Li, T., Zhang, H., Huang, Y., Su, Z., Deng, Y., Liu, H., et al. (2019b). Effects of the Fhb1 gene on Fusarium head blight resistance and agronomic traits of winter wheat. Crop J. 7, 799–808. doi: 10.1016/j.cj.2019.03.005
Livak, K. J., Schmittgen, T. D. (2001). Analysis of relative gene expression data using Real-Time quantitative PCR and the 2–ΔΔCt method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi: 10.1186/s13059-014-0550-8
Ma, X., Wang, H., Wu, S., Han, B., Cui, D., Liu, J., et al. (2024). DeepCCR: large-scale genomics-based deep learning method for improving rice breeding. Plant Biotechnol. J. 22 (10), 2691–2693. doi: 10.1111/pbi.14384
Ma, Z., Xie, Q., Li, G., Jia, H., Zhou, J., Kong, Z., et al. (2020). Germplasms, genetics and genomics for better control of disastrous wheat Fusarium head blight. Theor. Appl. Genet. 133, 1541–1568. doi: 10.1007/s00122-019-03525-8
MaChado, A. K., Brown, N. A., Urban, M., Kanyuka, K., Hammond-Kosack, K. E. (2018). RNAi as an emerging approach to control Fusarium head blight disease and mycotoxin contamination in cereals. Pest Manage. Sci. 74, 790–799. doi: 10.1002/ps.4748
Oliveros, J. C. (2015). “Venny,” in An interactive tool for comparing lists with Venn’s diagrams. Available at: https://bioinfogp.cnb.csic.es/tools/venny/index.html (Accessed May 25, 2024).
Pan, Y., Liu, Z., Rocheleau, H., Fauteux, F., Wang, Y., McCartney, C., et al. (2018). Transcriptome dynamics associated with resistance and susceptibility against fusarium head blight in four wheat genotypes. BMC Genomics 19, 642. doi: 10.1186/s12864-018-5012-3
Perochon, A., Jianguang, J., Kahla, A., Arunachalam, C., Scofield, S. R., Bowden, S., et al. (2015). TaFROG encodes a pooideae orphan protein that interacts with SnRK1 and enhances resistance to the mycotoxigenic fungus Fusarium graminearum. Plant Physiol. 169, 2895–2906. doi: 10.1104/pp.15.01056
Pertea, G., Pertea, M. (2020). GFF utilities: gffRead and gffcompare. F1000Research 9, 304. doi: 10.12688/f1000research.23297.1
Pieterse, C. M., van der Does, D., Zamioudis, C., Leon-Reyes, A., Van Wees, S. C. (2012). Hormonal modulation of plant immunity. Annu. Rev. Cell Dev. Bi. 28, 489–521. doi: 10.1146/annurev-cellbio-092910-154055
Ramaroson, M. L., Koutouan, C., Helesbeux, J. J., Le Clerc, V., Hamama, L., Geoffriau, E., et al. (2022). Role of phenylpropanoids and flavonoids in plant resistance to pests and diseases. Molecules 27, 8371. doi: 10.3390/molecules27238371
Rehmsmeier, M., Steffen, P., Hochsmann, M., Giegerich, R. (2004). Fast and effective prediction of microRNA/target duplexes. RNA 10, 1507–1517. doi: 10.1261/rna.5248604
Salmena, L., Poliseno, L., Tay, Y., Kats, L., Pandolfi, P. P. (2011). A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell 146, 353–358. doi: 10.1016/j.cell.2011.07.014
Song, L., Fang, Y., Chen, L., Wang, J., Chen, X. (2021). Role of non-coding RNAs in plant immunity. Plant Commun. 2, 100180. doi: 10.1016/j.xplc.2021.100180
Soresi, D., Díaz, M., Bagnaresi, P., Gallo, C., Carrera, A. (2023). Gene and lncRNA expression patterns associated with Qfhs.ndsu-3AS Fusarium head blight resistance QTL in durum wheat. Can. J. Plant Pathol. 45, 41–56. doi: 10.1080/07060661.2022.2091663
Su, Z., Bernardo, A., Tian, B., Chen, H., Wang, S., Ma, H., et al. (2019). A deletion mutation in TaHRC confers Fhb1 resistance to Fusarium head blight in wheat. Nat. Genet. 51, 1099–1105. doi: 10.1038/s41588-019-0425-8
Tang, D., Chen, M., Huang, X., Zhang, G., Zeng, L., Zhang, G., et al. (2023). SRplot: A free online platform for data visualization and graphing. PloS One 18, e0294236. doi: 10.1371/journal.pone.0294236
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., Van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
Van Breusegem, F., Bailey-Serres, J., Mittler, R. (2008). Unraveling the tapestry of networks involving reactive oxygen species in plants. Plant Physiol. 147, 978–984. doi: 10.1104/pp.108.122325
Walker, P. L., Belmonte, M. F., McCallum, B. D., McCartney, C. A., Randhawa, H. S., Henriquez, M. A. (2024). Dual RNA-sequencing of Fusarium head blight resistance in winter wheat. Front. Plant Sci. 14. doi: 10.3389/fpls.2023.1299461
Walter, S., Kahla, A., Arunachalam, C., Perochon, A., Khan, M. R., Scofield, S. R., et al. (2015). A wheat ABC transporter contributes to both grain formation and mycotoxin tolerance. J. Exp. Bot. 66, 2583–2593. doi: 10.1093/jxb/erv048
Wan, Y. F., Yen, C., Yang, J. L. (1997). Sources of resistance to head scab in Triticum. Euphytica 94, 31–36. doi: 10.1023/A:1002982005541
Wang, R., Huang, J., Liang, A., Wang, Y., Mur, L. A. J., Wang, M., et al. (2020b). Zinc and copper enhance cucumber tolerance to fusaric acid by mediating its distribution and toxicity and modifying the antioxidant system. Int. J. Mol. Sci. 21, 3370. doi: 10.3390/ijms21093370
Wang, X., Li, G., Jia, H., Cheng, R., Zhong, J., Shi, J., et al. (2024). Breeding evaluation and precise mapping of Fhb8 for Fusarium head blight resistance in wheat (Triticum aestivum). Plant Breed. 143, 26–33. doi: 10.1111/pbr.13113
Wang, H., Sun, S., Ge, W., Zhao, L., Hou, B., Wang, K., et al. (2020a). Horizontal gene transfer of Fhb7 from fungus underlies Fusarium head blight resistance in wheat. Science 368, 822–823. doi: 10.1126/science.aba5435
Wang, L., Wang, L., Li, Q., Liu, Z., Surendra, A., Pan, Y., et al. (2018). Integrated transcriptome and hormone profiling highlight the role of multiple phytohormone pathways in wheat resistance against fusarium head blight. PloS One 13, e0207036. doi: 10.1371/journal.pone.0207036
Wang, H., Yan, S., Wang, W., Cheng, Y., Hong, J., He, Q., et al. (2024). Cropformer: An interpretabl deep learning framework for crop genome prediction. Plant Commun. doi: 10.1016/j.xplc.2024.101223
Xiao, J., Jin, X. H., Jia, X. P., Wang, H. Y., Cao, A. Z., Zhao, W. P., et al. (2013). Transcriptome-based discovery of pathways and genes related to resistance against Fusarium head blight in wheat landrace Wangshuibai. BMC Genom. 14, 197. doi: 10.1186/1471-2164-14-197
Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., et al. (2011). KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322. doi: 10.1093/nar/gkr483
Yang, X., Zhang, L., Yang, Y., Schmid, M., Wang, Y. (2021). miRNA mediated regulation and interaction between plants and pathogens. Int. J. Mol. Sci. 22, 2913. doi: 10.3390/ijms22062913
Yin, J., Han, X., Zhu, Y., Fang, Z., Gao, D., Ma, D. (2022). Transcriptome profiles of circular RNAs in common wheat during Fusarium head blight disease. . Data 7, 121. doi: 10.3390/data7090121
Young, M. D., Wakefield, M. J., Smyth, G. K., Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11, R14. doi: 10.1186/gb-2010-11-2-r14
Yu, Y., Zhou, Y., Feng, Y., He, H., Lian, J., Yang, Y., et al. (2019). Transcriptional landscape of pathogen-responsive lncRNAs in rice unveils the role of ALEX1 in jasmonate pathway and disease resistance. Plant Biotechnol. J. 18, 679–690. doi: 10.1111/pbi.13234
Zaynab, M., Fatima, M., Abbas, S., Sharif, Y., Umair, M., Zafar, M. H., et al. (2018). Role of secondary metabolites in plant defense against pathogens. Microb. Pathogenesis 124, 198–202. doi: 10.1016/j.micpath.2018.08.034
Zhang, S., Shen, S., Yang, Z., Kong, X., Liu, F., Zhen, Z. (2020). Coding and Non-coding RNAs: Molecular basis of forest-insect outbreaks. Front. Cell Dev. Biol. 8. doi: 10.3389/fcell.2020.00369
Zhang, Y. D., Yang, Z. B., Ma, H. C., Huang, L. Y., Ding, F., Du, Y. Y., et al. (2021). Pyramiding of Fusarium head blight resistance quantitative trait loci, Fhb1, Fhb4, and Fhb5, in modern Chinese wheat cultivars. Front. Plant Sci. 12. doi: 10.3389/fpls.2021.694023
Zhang, F., Zhang, H., Liu, J., Ren, X., Ding, Y., Sun, F., et al. (2024). Fhb9, a major QTL for Fusarium head blight resistance improvement in wheat. J. Integr. Agr. doi: 10.1016/j.jia.2024.03.045
Keywords: Sichuan wheat landraces, Fusarium head blight, lncRNAs, circRNAs, miRNAs, RNA sequencing
Citation: Wu L, Wang J, Shen S, Yang Z and Hu X (2025) Transcriptomic analysis of two Chinese wheat landraces with contrasting Fusarium head blight resistance reveals miRNA-mediated defense mechanisms. Front. Plant Sci. 16:1537605. doi: 10.3389/fpls.2025.1537605
Received: 01 December 2024; Accepted: 13 February 2025;
Published: 28 February 2025.
Edited by:
Hao Zhou, Sichuan Agricultural University, ChinaReviewed by:
Shenghui Zhou, Chinese Academy of Agricultural Sciences (CAAS), ChinaYongming Chen, Peking University, China
Copyright © 2025 Wu, Wang, Shen, Yang and Hu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xinkun Hu, aHV4aW5rdW4xMjNAMTYzLmNvbQ==