Integration of multi-omics data accelerates molecular analysis of common wheat traits

Plant materials and growth conditions

We selected a common wheat, Yunong 268, released in 2023, as a middle-strong gluten wheat cultivar owing to its excellent comprehensive agronomic traits and great resistance to disease and stress. Yunong 268 was grown in pot culture at 25/15 °C day/night before the heading stage and at 32/26 °C day/night until the filling stage with a 16/8 h light/dark photoperiod in a greenhouse. The positions of plants in the growth chamber were randomly shuffled during plant growth to minimize the microenvironmental effect. The plants with normal and consistent growth were sampled for each experiment. We collected 20 different developmental stages across wheat vegetative and reproductive phases, including roots and seedlings at the seedling stage (SS_RT, SS_SG), leaves, roots, and the basal part of stems at the jointing stage (JS_LF, JS_RT, JS_BPST), leaves, roots, basal part of stems, and spikes at the booting stage (BS_LF, BS_RT, BS_BPST, BS_SP), flag leaves, roots, basal part of stems, stalks, sheaths, and spikes at the heading stage (HS_LF, HS _RT, HS _BPST, HS_ SL, HS_ SH, HS _SP), and roots, basal part of stems, and seeds at 14 days post anthesis (DPA), 28 DPA, and 35 DPA at the filling stage (FS _RT, FS _BPST, FS_14d_seed, FS_28d_seed, FS_35d _seed). After the samples were collected, immediately flash frozen in liquid nitrogen, and stored at −80 °C for protein and RNA extraction.

Fusarium species isolation and inoculation on wheat

The wheat cultivar Yunong 268 was planted in a greenhouse and inoculated with a prevalent Chinese F. pseudograminearum (F. pg) isolate, namely FP-Ta5. FP-Ta5 is a highly aggressive strain isolated from infected wheat crowns in northwest of the Yellow and Huai wheat regions of China. The strain was routinely cultured on potato dextrose agar (PDA) at 25°C. When seedlings grew to 3 cm long, 8 ~ 10 seedlings per cultivar that grew uniformly were inoculated, then 0.4 g colonized millet grains were scattered to each plastic container⁴⁸. All seedlings were grown in a greenhouse under 16/8 h day/night conditions at 25/20°C day/night with 60–80% relative humidity. After F. pg inoculation for 0, 24, 48, 72, and 96 hours at the one-leaf stage, basal stems with three biology replicates were collected, and frozen in liquid nitrogen, and stored at −80 °C for further use.

Protein extraction of 20 samples from wheat developmental stages

To extract tissue proteins, lysis buffer (10 mM dithiothreitol, 1% TritonX-100, 50 μM PR-619, 50 mM NAM, 3 μM TSA, 2 mM EDTA, and 1% Protease Inhibitor Cocktail) was used⁵⁷. The proteins were re-dissolved in 8 M urea, and a BCA Protein Assay Kit (Catalog No. P0012, Beyotime, Beijing, China) was used to determine the protein concentration according to the manufacturer’s instructions.

Trypsin digestion

For initial digestion, trypsin was added at a ratio of 1:50 trypsin to protein mass overnight. Next, the sample was reduced with 5 mM dithiothreitol for 60 min at 37 °C and alkylated with 11 mM iodoacetamide for 45 min at room temperature in the dark. Finally, Strata X SPE column was used to desalt the peptides.

High-performance liquid chromatography (HPLC) fractionation

The tryptic peptides were fractionated using high-pH reverse-phase HPLC with Agilent 300 Extend C18 column (5 μm particles, 4.6 mm ID, 250 mm length). The process involved separating peptides into 60 fractions using a gradient of 2% to 60% acetonitrile in 10 mM ammonium bicarbonate (pH 10) over 60 min. The peptides were then combined into 24, 10, and four fractions for proteome, phosphoproteome, and acetylproteome, respectively, and dried by vacuum centrifuging.

Phospho-peptide enrichment

Peptide mixtures were subjected to incubation with a suspension of immobilized metal affinity chromatography (IMAC) microsphere, accompanied by vibration in the loading buffer (50% acetonitrile/0.5% acetic acid). In order to eliminate non-specifically adsorbed peptides, a sequential washing procedure was employed suing 50% acetonitrile/0.5% acetic acid and 30% acetonitrile/0.1% trifluoroacetic acid solutions on the IMAC microspheres. After adding an elution buffer containing 10% NH₄OH and shaking the mixture, the phospho-peptides supernatant was collected and lyophilized for LC-MS/MS analysis.

Acetyl-lysine enrichment

The tryptic peptides were dissolved in NETN buffer (100 mM NaCl, 1 mM EDTA, 50 mM Tris-HCl, 0.5% NP-40, pH 8.0) and incubated with anti-acetyl-lysine antibody agarose conjugated antibody beads (PTM Biolabs, Lot number PTM-104) at 4 °C overnight with gentle shaking to enrich acetylation-modified peptides. The bound peptides were eluted from the beads using 0.1% trifluoroacetic acid, and then the beads were washed four times in NETN buffer and twice in deionized water. The eluted fractions were combined and vacuum-dried. The resulting peptides were cleaned with C18 ZipTips (Millipore) as per the manufacturer’s instructions.

Tandem mass spectroscopy (LC-MS/MS) analysis

The tryptic peptides were dissolved in solvent A (0.1% formic acid, 2% acetonitrile in water) and then directly loaded onto a reversed-phase analytical column (25 cm length, 100 μm i.d.). For the proteome, peptides were separated using a nanoElute UHPLC system (Bruker Daltonics) with a 6% to 24% gradient of solvent B (0.1% formic acid in acetonitrile) over 42 min, 24% to 32% over 12 min, climbing to 80% over 3 min, and held at 80% for the last 3 min, with a constant flow rate of 450 nL/min the whole time. For the phosphoproteome, peptides were separated with a 2% to 22% gradient of solvent B over 44 min, 22% to 35% over 10 min, climbing to 90% over 3 min, and held at 90% for the last 3 min, with a constant flow rate of 450 nL/min on a nanoElute UHPLC system (Bruker Daltonics). For the acetylproteome, peptides were separated with a 6% to 24% gradient of solvent B over 40 min, 24% to 32% over 12 min, climbing to 80% over 4 min, and held at 80% for the last 4 min, all at a constant flow rate of 450 nL/min on a nanoElute UHPLC system (Bruker Daltonics).

The peptides were subjected to a capillary source followed by timsTOF Pro (Bruker Daltonics) mass spectrometry (4D Mass Spectrometer). For the proteome, the electrospray voltage applied was 1.70 kV. Precursors and fragments were analyzed with a TOF detector, with an MS/MS scan range of 100–1700 m/z. The timsTOF Pro was operated in parallel accumulation serial fragmentation (PASEF) mode. Precursors with charge states 0 to 5 were selected for fragmentation, with 10 PASEF-MS/MS scans acquired per cycle. The dynamic exclusion was set to 30 s. For the phosphoproteome, the electrospray voltage applied was 1.75 kV, with the dynamic exclusion set to 24 s and the other parameters in line with the proteome. For the acetylproteome, the dynamic exclusion was set to 24 s, with the other parameters in line with the proteome.

Database search

The resultant data from MS/MS were processed using the MaxQuant search engine (v.1.6.15.0). Tandem mass spectra were searched against the latest version (IWGSC RefSeq v2.1)¹⁹ of the International Wheat Genome Sequencing Consortium (IWGSC) concatenated with reverse decoy database. Trypsin/P was specified as the cleavage enzyme, allowing up to two missing cleavages. The mass tolerance for precursor ions was set as 20 ppm for both the initial and subsequent searches, with the mass tolerance for fragment ions set as 0.02 Da. Peptide length was restricted to no less than 7, and no more than 5 peptide modifications were allowed. A fixed modification for the proteome included carbamidomethyl on cysteine (Cys), whereas variable modifications included acetylation on the protein N-terminal and oxidation on methionine (Met). For the phosphoproteome, carbamidomethyl on Cys was specified as a fixed modification, with acetylation on the protein N-terminal, oxidation on Met, and phosphorylation on serine (Ser), threonine (Thr), and tyrosine (Tyr) specified as variable modifications. For the acetylproteome, a fixed modification included a carbamidomethyl on Cys and three variable modifications included acetylation on the protein N-terminal, oxidation on Met, and acetylation on lysine (Lys). The false discovery rate (FDR) was adjusted to <1%. When calculating total protein, if unique peptides = 0, the total number of peptides identified should be at least 2, and a protein should be determined by two peptides to ensure the reliability of the identification; If unique peptides ≥ 1, we did not screen total peptides.

Fusarium crown rot responsive multi-omics

To find Fusarium crown rot (FCR) responsive proteins (proteins with ≥ 2 unique peptides), the wheat differential multi-omics were conducted with three biological replicates to examine the differentially expressed proteins (DEPs), phosphorylation proteins (DEPPs), and acetylation proteins (DEAPs) after F. pg inoculation for 0, 24, 48, 72, and 96 hours using the above method. Fold changes were determined using the abundance ratio of these proteins identified at 24, 48, 72, and 96 hours relative to those at 0 hour. Moreover, significant differences were determined between means of control (0 hour) and each of four inoculation times (24, 48, 72, 96 hours), respectively, using a two-sample t-test. DEPs, DEPPs, and DEAPs were defined based on thresholds of >1.5- or <0.67-fold change ratios and P < 0.05 at 24, 48, 72, and 96 hours compared to those at 0 hour.

RNA sequencing

Trizol® Reagent Plant RNA Reagent (Invitrogen) was used to extract total RNA from 20 tissue and FCR assay samples with three biological replicates. The sequencing depth is 10 G of a single sample. After the DNase I (TaKara) treatment, RNA-seq libraries were constructed using a TruSeq^TM RNA sample preparation Kit (Illumina, USA), with 150 bp paired-end sequences generated using an Illumina DNBSEQ-T7 sequencer. SeqPrep and Sickle with default parameters were used to trim and quality control the raw paired-end reads. The clean reads were aligned to the reference genome (IWGSC RefSeq v2.1)¹⁹ with orientation mode using HISAT2 software⁵⁸. StringTie was used to assemble the mapped reads of each sample in a reference-based approach⁵⁹. According to the transcripts per million reads (TPM) method, the expression level of each transcript was calculated by RSEM⁶⁰ to quantify gene abundance. DESeq2⁶¹ provided statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg’s approach to control the FDR. Genes with |log₂FC | >1 and FDR < 0.05 were considered differentially expressed genes (DEGs).

Transcript abundance distributions

For all 132,570 transcripts, a kernel density estimation of log₂ (TPM) distribution was calculated. Similarly, other kernel density estimates were calculated for the subsets of transcripts with 30,772 proteins, 18,912 phosphoproteins, and 11,751 acetylproteins. The distributions were plotted using R.

Tissue-wise mRNA-to-protein correlations

Pearson’s correlation coefficient (r) was calculated for the gene expression of all pairwise tissue combinations to assess tissue similarity. Correlations were displayed as separate heat maps at the protein level and transcript level. For the shared data, including 27,149 shared_transcript, 4002 shared_protein, and 2751 shared_intersection data, the Pearson’s correlation values between protein and transcript were calculated using genes from the 2,751 shared_intersection dataset.

Gene expression analysis

Gene expression analysis was conducted at transcript and protein levels¹⁰. We used proteins (intensity-based absolute quantification, iBAQ) and transcripts (TPM) from highest to lowest for the cumulative abundance calculation of eight representative tissues. The six most abundant transcripts or proteins are listed. The contribution of each gene representing the most abundant protein in at least one tissue was determined by taking the associated protein abundance value (iBAQ) and dividing it by the total number of protein abundance values (∑iBAQ) in that tissue.

Bioinformatics

We constructed the circular proteome map (CPM) using Circos⁶² normalizing the organ LFQ values (Lt) for each protein, as follows:

where \({gen}{e}_{i,j}\) represents the organ LFQ value of the protein i in organ j or time j, and n represents the number of organs or times.

Shannon entropy was used to assign organ specificity to the proteome, phosphoproteome, and acetylproteome^63,64. We calculated a relative expression value for each protein, as follows:

where m represents the number of proteins.

For each protein, the specificity score was calculated based on its relative expression, as follows:

where i represents the protein i.

Each protein has a specificity score, with the density curve made to determine the specificity threshold and screen out the specific protein. Phosphoproteome and acetylproteome applied a five-tenth percentile cutoff. All statistical analyses were performed using R. Gene ontology (GO) enrichment was carried out, with only GO terms with P < 0.05 (p.adjust/p-value < 0.05) considered significant. We analyzed the model of amino acid sequences in specific positions of phospho-12-mers or acetyl-20-mers in all protein sequences using the Motif-X software⁶⁵. The wheat proteome was used as the background database, with the parameters set as occurrences = 20, Bonferroni corrected P = 0.000001, and other parameters set to default values. Subcellular localization was predicted using WoLFPSORT (http://www.genscript.com/wolf-psort.html). Protein pathways were annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database⁶⁶; P < 0.05 (p.adjust/FDR < 0.05) was the threshold for pathway significance. The protein domain functional description was annotated using InterPro⁶⁷ (http://www.ebi.ac.uk/interpro/). Protein-protein interactions of FCR responsive proteins were determined from the protein (STRING) database v10.5, with the network visualized in Cytoscape⁶⁸.

Prediction of kinase-substrate regulations

The prediction of kinase-substrate regulations using proteome data was conducted by Group-based Prediction System (GPS) 5.0⁶⁹ (http://gps.biocuckoo.cn), which indemnified the primary specificity based on the theory of short linear motifs around phosphorylation sites (p-sites). The corresponding kinase proteins in the kinase family were acquired through a comparison with the kinase sequence in the IEKPD2.0 database. A “medium” threshold was chosen in GPS 5.0. One kinase can regulate multiple substrates, and a phosphorylation site can be regulated by more than one kinase. Therefore, a regulatory network diagram of protein kinase-substrate in the same tissue was constructed.

Protein/mRNA relation

We determined the relationship between protein and mRNA according to the previously reported method¹⁰. In brief, for more than eight tissues (n = 18,194), we calculated the r between median protein (median iBAQ of 20 tissues) and median transcript (median TPM of 20 tissues abundance. Protein-to-mRNA ratio (PTR) values are calculated for each tissue individually by constructing the ratio between the corresponding transcript abundance (TPM) and protein abundance (iBAQ) for each gene. The Tissue PTR values are then utilized to determine the median and median absolute deviation (MAD) of PTR across all 20 developmental stages. Further, quantiles (Q1–Q5) were created using the PTR MAD values for the median PTR distribution. Both phosphoproteome and acetylproteome were developed in a similar way.

Co-expression analysis

Co-expression networks were constructed using the weighted gene co-expression network analysis (WGCNA) R package³³ as described in the prior study¹¹. We constructed a set of four networks to capture all available information and use all available genes detected in at least three tissues. These networks comprised 58,263 mRNA, 8,231 proteins, 15,930 phospho-sites, and 11,033 acetyl-sites. Due to excessive transcriptome data flux, one parameter used for transcript (minModuleSize = 150) differed from other omics (minModuleSize = 30), with all other parameters used for each network identical. To get rid of differences in sequencing depth between studies, the count expression level of each gene was normalized with the variance stabilizing transformation from DESeq2⁶¹. The soft power threshold was calculated to exceed a scale-free topology fit index of 0.9 for each network separately. The soft power threshold SEDs were: transcript = 16, proteome =18, phosphoproteome = 24, and acetylproteome = 1. The topographical overlap matrices (TOM) were calculated block-wise. The module function using TOM Type = “signed”, with the minimum module size set to 30. Co-expression modules were constructed through hierarchical clustering of the TOM distance (1-TOM) using the hclust function with method = “average”. A parameter merge Cut Height = 0.15 was used to merge similar modules. Using these four omics data sets, we identified 172 modules (including 39,558 hub genes) by WGCNA (Supplementary Data 7–9).

Identifying highly connected hub genes

The signedKME function of the WGCNA R package was employed to compute hub genes for each module. Genes that exhibited a stronger absolute correlation with the eigengene (absolute correlation coefficient > 0.8) were classified as hub genes.

Gene regulatory network

For the gene regulatory networks (GRNs) inference, we generated three individual networks from transcriptome, proteome and phospho-proteome in our atlas, respectively. Transcription factors (TFs) were predicted using PlantTFDB (Plant Transcription Factor Database)³⁰. A set of 24,955 transcripts, 24,442 proteins, and 7187 phospho-proteins (Supplementary Data 4) were used as potential targets regulated by the predicted TFs. These transcripts, proteins, and phospho-proteins were observed in at least 2 omics (TPM > 0, iBAQ > 0, intensity > 0) in at least six samples. TFs, which were composed of 2414 transcripts, 2729 proteins, and 1049 phospho-proteins detected in at least 2 omics in at least six samples, were selected potential regulators (Supplementary Data 4). Then, three full GRNs were constructed by GRNboost2⁷⁰ (version 0.1.6) with default parameters. For each inferred GRN, we retained only edges with an importance score ≥1. Based on previous study¹¹, regulatory edges presented in more than two GRNs were incorporated to the total GRN, where its importance score was determined by the average. GO functional enrichment (p. value < 0.05) was performed for target genes of some TFs by using an R Bioconductor package cluster profiler (version 4.4.4). To investigate the tissue-specific regulatory mechanism, we employed a sample-specific network (SSN) method⁷¹ to construct tissue-specific GRN based on transcriptome data (TFs) and the total GRN. The regulatory edges with P value ≤ 0.01 were regarded as tissue-specific. These huge data of GRNs were deposited in our WheatPro (https://www.csuligroup.com/WheatPro/#/).

Quantitative contribution of post-translational modifications and transcript

The data subset with genes detected in at least eight tissues on a protein level (n = 2,2941) was used for all further calculations involving the quantitative contribution of post-translational modifications (PTM) and transcripts (PTC). For the phosphoproteome and acetylproteome, Pearson’s correlation coefficients (PCC) for phosphoprotein and protein, transcript and protein, and acetyl-protein and protein were calculated; proteins with at least one kind of modification were retained (Fig. 3c). To standardize the relative ratio between PCCs, we normalized the PCC for each protein with absolute values, defining PTC as follows:

For proteins modified by phosphorylation and acetylation simultaneously, the type with the higher PTC was selected as the dominant modification type. Three linearization categories were defined based on PTC values: balanced (BL), with PTC from 40–60%, transcript dominant (TD), with PTC from 0–40%, and PTM dominant (PD), with PTC from 60–100%. GO enrichment was undertaken, with only GO terms with P < 0.05 considered significant.

Principal component analysis and correlation coefficient analysis

We used categorized protein abundance values to visualize sub-genome and tissue expression-relatedness. These expression values (iBAQ) were averaged across tissues and log₂-transformed. The principal component analysis employed singular value decomposition via the prcomp function in R. Pearson’s correlation coefficients were determined with the Euclidean distance and complete linkage methods used to implement hierarchical clustering.

Relative expression levels of sub-genome homologs across triads

The analysis focused exclusively on the gene triads with a 1: 1: 1 correspondence across the three homologous sub-genomes, including 22,097 syntenic triads (66,291 genes) (http://wheat.cau.edu.cn/TGT/download.html). For relative expression levels of sub-genome homologs across triad, the transcripts/proteins with gene ID.1 (for instance, TraesCS1A02G002100.1) were selected. Protein abundance was normalized (riBAQ) by the ratios of the corresponding protein abundance value (iBAQ) in all protein abundance values (∑iBAQ). A protein triad was defined when the sum of the A, B, and D sub-genome homologs was > 0.01 riBAQ. Using this criterion, we defined 10,695 proteins (3565 triads) that were considered expressed (Supplementary Data 16). To standardize the relative expression of each homolog across the triad, we normalized the absolute riBAQ for each protein within the triad as follows:

where A, B, and D represent the gene corresponding to the A, B, and D homologs in the triad. The normalized expression was calculated for each tissue, and the average across all tissues. In addition, we defined 56,757 genes (18,919 triads) considered expressed at the mRNA level (Supplementary Data 16)³⁵.

For the phosphoproteome and acetylproteome, modification abundance was normalized (rintensity) by the ratios of their corresponding modification abundance value (intensity) in all modification abundance values (∑intensity), with 5358 phosphoproteins (1786 triads) and 1275 acetylproteins (425 triads) considered expressed. Protein abundance was normalized (rintensity) by the ratios of their corresponding modification abundance value (intensity) in all modification abundance values.

Definition of homolog expression bias categories

The ideal normalized expression bias for the seven categories was defined in Supplementary Data 16. The homolog expression bias category for each triad by selecting the shortest distance³⁵ for each of the intermediate tissue as well as for the average across all expressed tissues (combined analysis).

RNA isolation and quantitative real-time PCR (qRT-PCR)

After isolating total RNA with TRIZOL isolation reagent (Yeasen), we utilized the one-step ReScript II RT All-in-One Mix (with dsDNase) (Nobelab) to conduct reverse transcription procedures in accordance with the manufacturer’s instructions. The endogenous control was the β-actin gene (GenBank No. AB181991) (Supplementary Data 29).

Subcellular localization in wheat protoplasts

pJIT163-Ubi: delta-1-pyrroline-5-carboxylate synthase (TaP5CS1, TraesCS1B03G0812300)-hGFP, pJIT163-Ubi: histone deacetylase 9 (TaHDA9, TraesCS2B02G309700.1) -hGFP (Supplementary Data 29) and pJIT163-Ubi:hGFP were transfected into wheat mesophyll protoplast cells using a PEG-mediated method⁷². After incubation for 20 h in the dark at 22 °C, confocal microscopy (Carl Zeiss, LSM710) was used to examine the protoplasts.

Yeast two-hybrid (Y2H) assays

The CDS of the TaP5CS1 gene was cloned into the bait vector pGBKT7. A cDNA library was constructed from seedlings of wheat variety Aikang 58⁷³ to screen interaction proteins of TaP5CS1 using a yeast two-hybrid (Y2H) assay according to the manufacturer’s protocol (Clontech, USA). Yeast colonies were amplified by PCR and sequenced by T7-specific primers (Supplementary Data 29). After sequencing, the client genes were identified. The CDS of a histone deacetylase (TaHDA9) was amplified and fused to the prey vector pGADT7 to confirm protein interaction with the bait in the Y2H screening (Supplementary Data 29). The bait and prey vectors were co-transformed into yeast strain AH109 and screened on SD/Leu-Trp and SD/Trp-Leu-His-Ade + X-a-Gal medium.

Luciferase complementation imaging assay (LCI)

The CDSs of TaHDA9 and TaP5CS1 were subcloned into pCAMBIA1300-cLUC and pCAMBIA1300-nLUC to form TaHDA9-cLUC and TaP5CS1-nLUC, respectively. The recombinant vectors were transformed into Agrobacterium strain GV3101. Different vector combinations were co-transfected into the lower epidermis of N. benthamiana leaves. TaP5CS1-nLUC and cLUC, nLUC and TaHDA9-cLUC, cLUC and nLUC were used as negative controls. Luminescence images were captured using a plant-living imaging system (Berthold, Night Shade LB 985).

Prokaryotic expression of TaP5CS1

The TaP5CS1 was cloned into the prokaryotic expression vector pGEX-6P-1 using BamH1 and Not1 as restriction enzyme cutting sites (Supplementary Data 29). E. coli BL21 (DE3) was used to express recombinant TaP5CS1. TaP5CS1 protein was purified using a GST-tag Protein Purification Kit (Beyotime, China, Catalog No. P2262). The collected TaP5CS1 recombinant protein was analyzed by immunoblotting using anti-GST (PTMbio, PTM 5046) and pan anti-acetyl-lysine (anti-Kac antibody, Cell Signaling Technology, Catalog No. 9681) antibodies.

Transient expression assays

35S: TaP5CS1-GFP-Flag and 35S: TaHDA9-GFP-His or 35S: GFP were injected into four-week-old tobacco (N. benthamiana) leaf cells. Then, the total proteins extracted in lysis buffer [20 mM Tris-HCl (pH 8.0), 2 mM DTT, 1% Triton X-100, 800 μM PMSF, 250 mM sucrose] containing protease inhibitor cocktail. After washing⁷³, we performed immunoprecipitation using anti-Flag antibody (Abcam, M20008) and Protein A/G PLUS-Agarose (Santa Cruz, Lot number sc-2003) to purify TaP5CS1-GFP-Flag proteins from the above extracts. Eluted proteins were detected by immunoblotting with anti-Flag (Abcam, M20008) and anti-Kac (Cell Signaling Technology, Catalog No. 9681) antibodies using the ECL SuperSignal system.

Protein extracts of tobacco leaves with transient expression of 35S: TaP5CS1-GFP-Flag and 35S: TaHDA9-GFP-His or 35S: GFP were analyzed by immunoblotting with anti-Flag antibody (Abcam, M20008) to examine whether TaHDA9 could regulate the expression of TaP5CS1 protein. Anti-actin antibody (Abbkine, Lot number A01050) was used as a loading control.

Co-immunoprecipitation (Co-IP) assays

We transfected 35S:TaP5CS1-GFP-Flag and 35S: TaHDA9-GFP-His or 35S: GFP into tobacco leaf cells to extract the total proteins. TaP5CS1-GFP-Flag proteins from the above extracts were IP with anti-Flag antibody (Abcam, M20008) and Protein A/G PLUS-Agarose (Santa Cruz, Lot number sc-2003). The eluted proteins were detected by immunoblotting with anti-Flag antibody (Abcam, M20008) and anti-His antibody (Abmart, M20001), respectively.

Pull-down protein interaction assays

TaHDA9-His and TaP5CS1-GSTor GST protein were expressed in E. coli. Purification of TaP5CS1-GST or GST proteins was accomplished with a GST-tag Protein Purification Kit (Beyotime, China, Catalog No. P2262) and beyoGold^TM GST-tag Purification Resin in buffer (50 mM Tris, 150 mM NaCl, 10 mM GSH, pH 8.0). The pulled-down proteins were mixed with the SDS sample buffer and detected by immunoblotting with anti-GST (PTMbio, PTM5046) and anti-His (Abmart, M20001) antibodies.

In vitro lysine deacetylase enzyme activity assay

To determine whether TaHDA9 can deacetylate TaP5CS1 in vitro, E. coli was transformed with TaP5CS1-GST and TaHDA9-His or a His vector. The TaP5CS1-GST recombinant protein was isolated from E. coli cultures that were grown both with and without the deacetylation inhibitors 5 nM trichostatin A (TSA, MedChemExpress, Catalog No. HY-15144) and 5 mM β-nicotinamide (NAM, MedChemExpress, Catalog No. HY-B0150). Anti-GST (PTMbio, PTM5046) and pan anti-Kac (Cell Signaling Technology, Catalog No. 9681) antibodies were used for immunoblotting detection of the isolated proteins.

Site-directed mutagenesis

The 35S: TaP5CS1-GFP-Flag plasmid was subjected to site-directed mutagenesis of K221, K231, and K634 (K221R, 231R, K634R, AAA-AGA, AAG-AGG, AAA-AGA) using the Fast Site-Directed Mutagenesis Kit (Tiangen, China, Catalog No. KM101). PCR was conducted using site-specific primers (Supplementary Data 29). Mutations were verified by DNA sequencing. Four-week-old tobacco (N. benthamiana) leaf cells were used to express recombinant wild-type and three mutants of TaP5CS1, respectively. Proteins were purified by IP and analyzed by immunoblotting with anti-Flag (Abcam, M20008) and anti-pan-Kac (Cell Signaling Technology, Catalog No. 9681) antibodies. In addition, TaP5CS1 and TaP5CS1^K634R were transiently expressed in tobacco, respectively, and analyzed by immunoblots with an anti-Flag antibody (Abcam, M20008) to examine the expression of wild-type and three mutants of TaP5CS1. An anti-actin antibody (Abbkine, Lot number A01050) was used as the loading control.

Barley strip mosaic virus (BSMV) induced gene silencing

Virus-induced gene silencing (VIGS) experiments were conducted on the middle-sensitivity to BSMV wheat cultivar Yunong 268. The primers for peroxidase (TaPOD, TraesCS2B03G1538300.1), sucrose synthase (TaSUS, TraesCS7D03G0082900.2), and TaP5CS1 were generated from 312, 259, and 232 base pair fragments, respectively, using Primer 3.0 (Supplementary Data 29). The α, β, and γ RNAs of the BSMV genome were synthesized from linearized plasmids, using RiboMAX^TM Large Scale RNA Production System-T7 (Promega, P1300) and Ribo m7G Cap Analog (Promega, P1712). The positive control used plants inoculated with BSMV: PDS (phytoene desaturase). Transcripts of each vector (α, β, γ, or recombinant γ-gene) were mixed in a 1:1:1 ratio for inoculating the wheat cultivar. Ten days after the wheat plants were infected with F. pg, the mixture RNA virus was added to the FES buffer for inoculation of the second leaves of wheat seedlings, and after darkness for 24 h, all plants were grown at 23 °C in 60–80% relative humidity. After two weeks, the FCR resistance of wheat plants was investigated, with the third and fourth leaf tissues were collected to determine the efficiency of silencing of TaPOD, TaSUS, and TaP5CS1.

Ethyl methanesulfonate mutants of TaHDA9

Seeds of the tetraploid wheat Kronos were mutagenized by Jorge Dubcovsky’s team using the chemical mutagen ethyl methanesulfonate (EMS)⁷⁴. The mutant lines P5CS1-1B (K3196: C/T, mutation effect = stop gained), HDA9-2A (K4505: C/T, mutation effect = splice_acceptor_variant), HDA9-2A (K4236: C/T, mutation effect = stop gained), HDA9-2B (Kronos733:A497G, mutation effect = splice_acceptor_variant), and HDA9-2B (Kronos2217, C/T, mutation effect = stop gained) were sequenced using specific primers (Supplementary Data 29). The positive mutant individual plants selected for backcrossing with the wild-type to obtain BC₂ mutants for functional verification.

Overexpressed transgenic plants of TaP5CS1 and TaHDA9

The CDSs of TaP5CS1 and TaHDA9 were cloned into the LGY-OE3 vector with the Ubi promoter to produce TaP5CS1-OE and TaHDA9-OE plants, respectively. The vectors containing targeted genes were transformed by Agrobacterium-mediated infection into immature embryos of Fielder (FCR-susceptible) to obtain TaP5CS1-OE and TaHDA9-OE lines. T₀ transgenic plants detected as positive by PCR using specific primers (Supplementary Data 29) were further self-pollinated. The T₂ non-segregating plants of three lines with high expression levels were selected for further analysis.

Disease index (DI) investigation

F. pg was routinely cultured on potato dextrose agar (PDA) at 25°C. Millet grains were used as the pathogen medium to inoculate the wheat seedlings. All seedlings were grown in a greenhouse under 16 h/8 h day/night conditions at 25/20°C day/night temperatures with 60-80% relative humidity. The severity of FCR was evaluated at the fourth week after being inoculated with F. pg and the resistance of the plants was estimated with the DI, which was calculated on a scale from 0 to 9 for every plant. The detailed method is based on our previous study⁴⁸. Briefly, a total of 12 infection stems were counted for each experiment. Five independent biological repeats were performed. The field test of the FCR resistance was using inoculated millet grains as a pathogen medium concurrent with the planting material, and the severity of the disease was investigated at wheat grain filling stage, which was calculated by the DI with a scale from 0 to 4⁷⁵.

Microscopic analysis of pathogen infection

To observe F. pg development on wheat plants, we sampled the stems of the overexpressed transgenic wheat and wild type after F. pg inoculation for four weeks, then stained them with 2.5% glutaraldehyde. The F. pg fungals were observed with a scanning electron microscope (TEM, Hitachi, HT7700). A total of 12 infection stems were counted from for each experiment. Three independent biological repeats were performed.

We measured the amount of F. pg biomass on wheat plants by taking samples from the stems of overexpressed transgenic wheat and wild type wheat that had been infected with F. pg for four weeks. We then extracted the total DNA. The semi-quantitative PCR was conducted using tubulin gene from wheat (GenBank No. MG852130.1) and F. pg (GenBank No. CP102997.1), respectively. A total of 12 infection stems were counted for each experiment. Three independent biological repeats were performed.

To assess the cellular response of the overexpressed transgenic wheat and wild type to F. pg infection, cell necrotic death at the infection site was stained by trypan blue³⁹. The inoculated stems were detached from the plants and stained in a 0.5% solution of trypan blue in lactophenol for 2 days at room temperature and mounted in a solution of 2.5 mg/mL chloral hydrate. The stained segments were observed and photographed by a stereomicroscope (Olympus, Tokyo, Japan). A total of 12 infection stems were counted for each experiment. Three independent biological repeats were performed.

Proline content measurement

Proline content was quantified in the overexpressed transgenic wheat and wild type using the acid ninhydrin assay. Briefly, proline was extracted by boiling 0.5 g of plant material in 2 mL of distilled water before adding 500 μL of 0.2 mM sodium citrate (pH 4.6) and 2 mL of 1% ninhydrin to 0.5 μL of plant extract. The mixture was boiled for 1 h before adding 2 mL toluene for extraction and centrifuging. A standard curve was prepared by measuring the absorbance of the standard proline (Himedia, India). Proline content was measured at 520 nm using a spectrophotometer (Shimadzu, Japan). All experiments were conducted in triplicate.

Assay of proline inhibition against F. pg growth

To evaluate the effect of exogenous proline on F. pg growth, 0, 1, 3, 5, 7, 9, and 10 mM proline were separately added to potato dextrose agar (PDA) medium. Since valine has a similar structure to proline, it was selected as a control. Subsequently, F. pg mycelia were inoculated into each pore of PDA and incubated at 25°C for 6 days. From the first to sixth day post-inoculation (DPI), we recorded the diameter of the mycelial colony and hyphal morphology. The mycelium growth inhibition rate was determined by measuring the colony diameter⁷⁶:

A total of five treated PDA were counted for each experiment. Three independent biological repeats were performed.

Effect of exogenous proline on wheat Fusarium crown rot resistance

After wheat germination, 10, 20, and 30 mM proline were separately added to wheat plot culture and an equal amount water was used as a control. After F. pg inoculation for four weeks, disease phenotype of wheat seedlings was evaluated, with DI calculated as above. Twelve infection stems were counted for each experiment. Five independent biological repeats were performed.

Quantification and statistical analysis

Two-tailed Student’s t-test (*P < 0.05, **P < 0.01) was used for statistical analysis of experimental and control groups. As for multiple samples, one-way analysis of variance (ANOVA) with Tukey’ multiple comparisons test was conducted.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.