Transcriptomic Analysis Reveal the Molecular Mechanisms of Wheat Higher-Temperature Seedling-Plant Resistance to Puccinia striiformis f. sp. tritici.

Stripe rust, caused by Puccinia striiformis f. sp. tritici (Pst), is a destructive disease of wheat worldwide. The disease is preferably controlled by growing resistant cultivars. Wheat cultivar Xiaoyan 6 (XY 6) has been resistant to stripe rust since its release. In the previous studies, XY 6 was found to have higher-temperature seedling-plant (HTSP) resistance. However, the molecular mechanisms of HTSP resistance were not clear. To identify differentially expressed genes (DEGs) involved in HTSP resistance, we sequenced 30 cDNA libraries constructed from XY 6 seedlings exposed to several temperature treatments. Compared to the constant normal (15°C) and higher (20°C) temperature treatments, 1395 DEGs were identified in seedlings exposed to 20°C for 24 h (to activate HTSP resistance) and then kept at 15°C. These DEGs were located on all 21 chromosomes, with 29.2% on A, 41.1% on B and 29.7% on D genomes, by mapping to the Chinese Spring wheat genome. The 1395 DEGs were enriched in ribosome, plant-pathogen interaction and glycerolipid metabolism pathways, and some of them were identified as hub proteins (phosphatase 2C10), resistance protein homologs, WRKY transcription factors and protein kinases. The majority of these genes were up-regulated in HTSP resistance. Based on the differential expression, we found that phosphatase 2C10 and LRR receptor-like serine/threonine protein kinases are particularly interesting as they may be important for HTSP resistance through interacting with different resistance proteins, leading to a hypersensitive response.

Introduction

Stripe (yellow) rust, caused by Puccinia striiformis f. sp. tritici (Pst), is a destructive disease affecting wheat production world wide (Wan, 2003; Chen, 2005). Breeding resistant cultivars is the best approach for controlling stripe rust (Zhang et al., 2001; Chen, 2007; Sui et al., 2009). Different types of stripe rust resistance have been identified and used for developing resistant cultivars (Chen, 2005, 2013). Based on specificity, resistance can be classified as race-specific and non-race-specific. Race-specific resistance is usually controlled by major genes and effective throughout plant development. However, new virulent Pst races can overcome race-specific resistance (Chen, 2005; Zheng et al., 2013). For example, the rapid development of Pst races that have overcome Yr2, Yr9, Yr17, and Yr27 resistance has led to destructive epidemics in many parts of the world (Wellings, 2011). In contrast, non-race-specific resistance is usually quantitative and often controlled by several genes (Coram et al., 2008a; Chen, 2013).

High temperature resistance to Pst is activated by changes in temperature and is believed to be non-race-specific. Use of such temperature induced resistance could thus be considered as a durable method for managing stripe rust (Shang, 1998; Ma and Shang, 2000; Chen, 2013; Zhou et al., 2014). Resistance in both seedling and adult plants can be induced by temperature changes. High-temperature adult plant (HTAP) resistance has been successfully used to develop durable resistant cultivars in the United States since the early 1960s (Chen and Line, 1995; Line, 2002; Chen, 2005, 2013). Cultivars with only HTAP resistance are susceptible in seedlings when temperatures are low (diurnal temperatures changing from 4 to 20°C), but gradually become more resistant when plants grow older and temperatures are higher (diurnal temperatures changing from 10 to 30°C; Chen, 2013). HTAP resistance usually becomes visible after the tillering stage and reaches to the highest level on the flag leaves (Qayoum and Line, 1985; Milus and Line, 1986; Chen, 2013). Numerous genes or quantitative trait loci conferring HTAP resistance have been identified and used to develop wheat cultivars with durable resistance. HTAP resistance is generally partial and can have a wide range of levels, depending on individual genes and the number of genes in a cultivar (Chen, 2013). Although HTAP resistance is influenced by temperature and growth stage, different HTAP resistance genes may have different sensitivities to temperature and/or plant growth stage. Similar to HTAP resistance, higher-temperature seedling-plant (HTSP) is also induced by higher-temperature. However, typical HTSP resistance is not affected much by plant growth stages. HTAP resistance is often reversible as plants become susceptible or less resistant when temperature changes from high to low (Qayoum and Line, 1985; Chen, 2013). In contrast, at least with the wheat cultivars studied, seedlings with HTSP resistance continue showing resistance after exposure to 18~21°C for only 24 h (Lu and Li, 1958; Lu, 1996; Ma and Shang, 2000; Hu X. P. et al., 2012; An et al., 2015).

Winter wheat cultivar Xiaoyan 6 (XY 6), developed from a cross between a wheat (Triticum aestivum) cultivar and Elytrigia elongatum (Li, 1986), has shown partial resistance to stripe rust, and the resistance has been characterized as HTSP resistance (Ma and Shang, 2000; Hu X. P. et al., 2012). An et al. (2015) found that the treatments of 18–24°C after inoculation of seedlings significantly reduced infection type and uredospore production compared to the seedlings grown at constant 16°C. At 8 days after inoculation when plants had the mosaic symptom without sporulation, the plants exposed to the optimal temperature of 20°C for 24 h showed incompatible reaction, in contrast to the compatible reaction on the plants without the higher temperature treatment. These results show that XY 6 has HTSP resistance to stripe rust and this type of resistance is induced by higher-temperatures.

HTAP resistance has been found to involve different mechanisms. The HTAP resistance controlled by Yr36 can be observed at the seedling stage under high temperatures, but the highest level of resistance is expressed at the adult-plant stage at high temperatures (Uauy et al., 2005; Chen, 2013). Originally from Triticum dicoccoides, Yr36 encodes a predicted kinase and a steroidogenic acute regulatory protein-related lipid transfer (START) domain (Fu et al., 2009). Both the kinase and START domains are necessary for the resistance function. Temperature and Pst inoculation consistently up-regulates expression of the resistance alleles, but down-regulates the susceptible alleles. The START domain has the ability to bind lipids from stripe rust fungus at high temperature and change its conformation, which may cause the kinase domain to initiate a signaling cascade leading to programmed cell death. Yr18 (also known as Lr34) is considered as a HTAP resistance gene as the level of resistance is increased by high temperatures (Chen, 2013). This gene encodes a putative ATP-binding cassette (ABC) transporter (Krattinger et al., 2009a,b). This drug resistance gene product contains two cytosolic nucleotide binding domains (NBD) and two hydrophobic transmembrane domains. The pleotropic resistance gene in wheat confers non-race-specific resistance to stripe rust, leaf rust, stem rust and powdery mildew, and also confers resistance to other diseases when transferred into barley, corn, rice and other plant species (Krattinger et al., 2011, 2013, 2016; Risk et al., 2013). Yr46 is also a pleiotropic gene, providing adult-plant resistance to stripe rust, leaf rust (Lr67), stem rust (Sr55) and powdery mildew (Pm46) (Herrera-Foessel et al., 2011; Chhetri et al., 2016). This gene encodes a hexose transporter that differs from the susceptible form of the same protein by just two conserved amino acids (Moore et al., 2015). The susceptible allele functions as a higher affinity glucose transporter, while the resistant allele has a dominant-negative effect through heterodimerization with functional transporters to reduce glucose uptake. These cloned adult-plant or HTAP resistance genes do not have LRR domains and do not show race specificity. However, nine NBS-LRR genes are involved in Yr39-controlled HTAP resistance based on transcriptomics analyses (Lin and Chen, 2007; Coram et al., 2008a). In addition to HTAP resistance, transcriptomics analyses have been used to study mechanisms of race-specific all-stage resistance and numerous genes with diverse functions have been found to be involved in this type of resistance (Coram et al., 2008b, 2010; Chen et al., 2013; Zhang et al., 2014; Hao et al., 2016). The effects of temperature on plant defense have been studied for other diseases (Wang et al., 2017). In contrast, there were no reports on molecular mechanisms on HTSP resistance to stripe rust, and the identity of genes and biochemical pathways involved in HTSP resistance were unknown before the present study.

The main objective of this study was to identify co-regulated genes that show significant changes in expression patterns related to HTSP resistance. We confirmed that exposure of XY 6 seedlings to 20°C for 24 h was sufficient for activating the resistance to Pst, and used this exposure regime to study gene expression during the activation of HTSP resistance in comparison with the inoculated seedlings grown at constant temperatures of 15 and 20°C. Through the comparison, we identified a large number of differentially expressed genes (DEGs) induced by the higher-temperature treatment. These DEGs allowed us to infer the mechanisms underlying HTSP resistance.

Materials and methods

Plant materials, growth conditions, and temperature treatments

Chinese yellow rust race 32 (CYR32) was used to inoculate wheat cv. XY 6 (susceptible but possessing HTSP) and Mingxian 169 (MX 169, susceptible without HTSP). Seeds (10–15) were sown in plastic pots (10 × 10 × 10 cm3) at a seed-to-seed distance of ca. 1.5 cm. Urediniospores of CYR32 were added to sterile water at a ratio of ~1:6–9 (v/v) and stirred with a vaccination needle; urediniospores were floating at the top of water surface and this spore suspension was then used to inoculate seedlings. For each cultivar, a total of 120 pots (90 inoculated and 30 not inoculated) were used for each of three biological replicate experiments over time. At the one-leaf stage (~10 days after sowing), seedlings were inoculated with the urediniospore suspension using a paint brush and then kept in a growth chamber (Percival E-30B, Perry, IA, USA) in dark at 10 ± 1°C and 80% relative humidity for 24 h, as described previously (Wang et al., 2017). Thereafter, the inoculated seedlings were divided into three groups for exposure to different temperature regimes. The first group was maintained at 15 ± 1°C [the normal temperature treatment (N)]; the second group was for the normal-higher-normal temperature treatment (NHN) – seedlings were kept at 15 ± 1°C from two to eight days after inoculation (dpi) and then transferred at 20 ± 1°C for 24 h, and finally moved back to 15 ± 1°C. This regime was shown previously to activate HTSP in XY 6 (An et al., 2015). The third group was maintained at 20 ± 1°C (H) 24 h after Pst inoculation. Wheat plants inoculated with sterile water were used as controls (Figure S1). The first sampling time (0 h) corresponded to the beginning of NHN treatment (at 8 dpi, i.e., 192 h after Pst inoculation). Leaf tissue from seedlings under three temperature treatments was sampled at 0, 12, 24, 48, and 120 h after temperature treatment was imposed (i.e., 192~312 h after Pst inoculation). For each cultivar, leaf tissue samples were collected at each time point for each biological replicate under each treatment.

Histopathological analysis

Pst inoculated leaves (XY 6 and MX 169) sampled at 12, 24, 48, and 120 h were assessed for the number of necrotic cells per infection site, the length of uredinium and the number of uredinia per leaf under a microscope. The number of necrotic cells per infection site and the length of uredinium were measured based on a published method (Wang et al., 2017). To measure the number of uredinia per leaf, 10 leaves were randomly selected for each of the three biological replicates. Microscopic observations were performed using an Olympus BX-51 microscope (Olympus Corporation, Tokyo, Japan) or an Olympus SZ-PT anatomical lens (Olympus Corporation, Tokyo, Japan); the data were measured using Cell Sens Entry software (v.1.6). The length and number of uredinia per leaf as well as the number of necrotic cells were analyzed using analysis of variance (ANOVA) and multiple comparison tests, which were performed using a generalized linear model with a Poisson distribution and the Tukey test by the glm and glht functions of R software (v.3.2.3), respectively.

Establishment and sequencing quality evaluation of cDNA libraries

The leaves of XY 6 sampled at 0 and 24 h (8 and 9 dpi, respectively) were used for RNA-Seq (Figure S1). Samples at 0 h from N and NHN treatments were the same, as the samples had not yet been subjected to higher temperature. For each of the three repetitions, one sample was selected and RNA was extracted under each combination of treatment and time point, a total of 30 samples from XY 6 were collected to extract RNA. RNA extracting was conducted using the PureLink® Plant RNA Reagent (Invitrogen, Carlsbad, CA, USA) and then treated with DNase I (Thermo Fisher, Waltham, MA, USA) at a concentration of 1 U/mg. The quality and concentration of extracted RNA were checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany). Thirty paired-end (PE) cDNA libraries were sequenced separately and PE reads were generated separately with Q30 as a base phred quality score threshold (Each library > 4 Gb and single plexing). Sequencing was performed on each library from each sample to generate 100 bp PE reads for transcriptome sequencing on an Illumina High-Seq 2000 platform. Library construction was accomplished using commercial products (Illumina); sequencing was done by Macrogen (Seoul, South Korea).

Optimization and evaluation of CDMC assembly strategy

RNA-Seq data trimming and adapter clipping were performed using Trimmomatic (v.0.33; Bolger et al., 2014). Reads were assembled by de novo, combined with reference-based mapping using CD-HIT-EST (CDMC) to remove redundancy (Li et al., 2013). The de novo process was conducted by reconstructing a transcript library with reads using Trinity (v.2.0.6; Haas et al., 2013). The mapping process was carried out by mapping the reads to the Chinese Spring wheat genome (T. aestivum ftp://ftp.ensemblgenomes.org/pub/release-8/plants/fasta/triticum_aestivum/dna/) using Tophat (v.2.1.0). Then, a transcript library was reconstructed using Cufflinks (v.2.2.1) and Coffmerge (v.2.2.1; Trapnell et al., 2012). The CDMC process was conducted by combining the above two libraries constructed from de novo and mapping strategies using CD-HIT-EST (v.4.6.4) with a similarity threshold of 95% identity (Fu et al., 2012) to generate a new transcript library.

Analysis of differentially expressed genes

The normalization factors were calculated using the trimmed mean of M-values (TMM) method of RSEM (v.1.1.17; Li and Dewey, 2011). The TMM-FPKM value was the expression level of each transcript expressed as the fragments per transcript kilobase per million fragments mapped value by the TMM normalization of the RNA-Seq data. In this study, we focused on DEGs during the induction process of HTSP resistance to Pst. Thus, significant DEGs over time were not analyzed. To eliminate the time effect, the 0 h samples were used as background. TMM-FPKM values for 24-h samples were divided by the corresponding values for 0 h samples, to avoid division by zero, 1 was added to each value (Oono et al., 2013). To validate whether data met the requirement for a normal distribution, the Kolmogorov-Smirnov test was used (SPSS software, v.20; Table S1). HTSP resistance to Pst was induced by two factors (higher temperature and inoculation with Pst); thus the interaction between these two factors should be considered. Therefore, the DEGs were evaluated using the following model:

where y is the expression value for gene i; μ is the population mean expression; B, I and T are the indicator variables that describe the batch, inoculation and temperature treatment, respectively; α, β, and β are the batch effect, inoculation effect and temperature effect, respectively; IT is the interaction term for inoculation and temperature; β is the IT interaction effect; and εi is the random error following a normal distribution with mean of 0 and variance of σ2. The model terms were tested by ANOVA, and P-values for all genes were adjusted using a false discovery rate (FDR) of α < 0.05 (Benjamini and Hochberg, 1995). In addition, a log2-fold change > 1 or < −1 and a logCPM (log2 counts per million) > −2 were used as additional criteria to select DEGs (McCarthy et al., 2012; Liu Y. et al., 2014).

Identification of chromosomes of DEGs

Identified DEGs were mapped to the hexaploid wheat genome of cv. Chinese Spring (CS) (ftp://ftp.ensemblgenomes.org/pub/release-8/plants/fasta/triticum_aestivum/dna/). These mapped DEGs were categorized into two groups: (1) DEGs with chromosome location information from the reference wheat genome of CS and (2) DEGs without location information mapped from the de novo assembly were used as probes for use with BLASTN (with an E-value < 1E−50) against the predicted mRNA database of the CS genome to search location information on chromosomes (Laudencia-Chingcuanco et al., 2006). A circles-plot was made using R software (v.3.2.3).

Protein-protein interaction (PPI) network analysis

Proteins encoded by identified DEGs were used to analyse protein-protein interactions (PPIs) based on the STRING database of the model plant Arabidopsis thaliana (with an E-value < 1E−10) (http://string-db.org/). PPIs with combined confidence scores greater than 0.7 were selected (Franceschini et al., 2013; Liu et al., 2016) and visualized with CYTOSCAPE (v.2.8, http://cytoscape.org/; Shannon et al., 2003).

Functional annotation and enrichment analysis

For functional annotation, transcripts were subjected to BLASTX (v.2.2.28, E-value < 1E−5) analysis against several protein databases, including the Nr, Swiss-Prot and KEGG databases. The transcripts were named according to the annotation in the Nr database. BLAST2GO (Conesa et al., 2005) was then used to obtain GO annotations. To investigate the metabolic pathways of those annotated transcripts, the transcripts were aligned to the KEGG database. PPI network, GO terms and KEGG pathways with FDR(BH adjustment) corrected P-values smaller than 0.05 were considered statistically significant. R genes in DEGs were predicted (E-value < 1E−5) based on the information from the Plant Resistance Gene Database (Sanseverino et al., 2013). Transcription factors (TFs) in DEGs were predicted (E-value < 1E−5) according to the Plant Transcription Factor Database (Jin et al., 2014).

Quantitative reverse-transcription-PCR analysis

Twelve transcripts were randomly selected for qRT-PCR analysis (Table S2). UltraSYBR Mixture (Kangwei, Beijing, China) and iQ™ 5 (Bio-Rad, Hercules, CA, USA) were used for qRT-PCR analysis of all reactions according to the manufacturer's instructions. Data were collected from three replicate experiments—the samples used for qRT-PCR were the same as those used for RNA-Seq, each consisting of at least three technical repeats. Negative controls that lacked templates were also included. The amplification efficiency of primers was determined (Figure S2) using LinReg PCR (Ramakers et al., 2003). The wheat ATP-dependent 26S proteasome regulatory subunit (26S) and cell division control (CDC) genes were chosen as internal reference genes for each qRT-PCR assay (Paolacci et al., 2009; Scholtz and Visser, 2013). The relative expression of selected transcripts was calculated using the 2−ΔΔCT method (Livak and Schmittgen, 2001).

RNA-seq data submission

The raw data used in the present study for transcriptome assembly and gene expression analysis have been submitted to the NCBI Sequence Read Archive (SRA) database under accession numbers from SRR5580869 to SRR5580898.

Results

Histopathological observation of Pst infections

Among the normal temperature (N), normal-higher-normal temperature (NHN), and higher temperature (H) treatments, the necrosis on leaves at 16 dpi was observed only in cv. XY6 from the NHN treatment. There was no change in the infection type in cv. MX 169 in all three treatments (Figure 1A). After exposure to 20°C for 24 h, necrosis of host cells (NC) was observed around secondary hyphae (SH) at the infection site in the NHN treatment of XY 6 (Figure 1B). XY 6 had more (P < 0.05) necrotic cells per infection site, shorter (P < 0.05) uredinial length, and fewer (P < 0.05) uredinia per leaf than MX 169 in the NHN treatment at all-time points (Figures 2A–C,E). In addition, the number of necrotic cells per infection site and number of uredinia per leaf in XY 6 were greater (P < 0.05) in the NHN 24 h treatment than in the NHN 12 h (Figures 2D,F). These results confirmed that HTSP resistance was activated and that the hypersensitive response (HR) of XY 6 to Pst was induced by exposure to 20°C for 24 h.

Figure 1

Histopathology observation during Puccinia striiformis f. sp. tritici infection in XY 6 (susceptible but possessing HTSP) and MX 169 (susceptible without HTSP) under different temperature treatments. (A) The infection types on XY 6 and MX 169 under different temperature treatments were observed at 8, 12 and 16 days post-inoculation (dpi). (B) Leaves of XY 6 infected by Puccinia striiformis f. sp. tritici at normal temperature 0, 24 h (I-N-0 and I-N-24), normal-higher-normal temperature 24 h (I-NHN-24), are examined under an epifluorescence microscope. I-N-0: The substomatal vesicle (SV), infection hypha (IH), secondary hyphae (SH), and haustorial mother cells (HMC) formed at the infection site at 0 h under the N treatment. I-N-24: SH formed at the infection site, which extends rapidly and formed larger colonies, and then further produce a large number of uredinia (U) at 24 h under the N treatment. I-NHN-24: SH formed at the infection site at 24 h and further induces the necrosis of host cell (NC) under the NHN treatment. 0~24 h means the hours post-temperature treatment which represents 192~216 h after inoculation. Bars, 20~100 μm. The samples of I-N-0 and I-NHN-0 are same before the temperature treatment. The results of H treatment have not been shown, because the results are similar to the N treatment.

Figure 2

The influence of Puccinia striiformis f. sp. tritici development and the wheat response among different temperature treatments (black, red and blue represent normal temperature (N), normal-higher-normal temperature (NHN) and higher temperature (H) treatment, respectively) during the initial symptom expression stage of stripe rust development in XY 6 (Filled histograms) and MX 169 (Open histograms). 12~120 h means the hours post-temperature treatment which represents 204~312 h after inoculation. (A) The number of necrosis of host cell, (B) Linear length of uredinia and (C). The number of uredinia per leaf. The lowercase letters means significant difference between different temperatures treatment in the different cultivars at the same time. The (D–F) are necrosis number of host cell, linear length of uredinia and the number of uredinia per leaf at different time points under NHN treatment in XY 6, respectively. For each treatment, the error bars are the mean ± SE of three replicates.

Identification of DEGs and their chromosome locations

A total of 65.95 million 101-bp paired-end clean reads were obtained, 91.23% of which had quality scores at or greater than the Q30 level (Table S3). The CDMC assembly strategy yielded 445226 transcripts with an N50 length of 1,849 bp (Table S4). Approximately 81.1 and 2.8% of the transcripts were longer than 400 bp and 4,000 bp, respectively (Figure S3A). About 37.0 and 3.4% of the transcripts with a predicted open reading frame (ORF) were longer than 300 and 1,000 bp, respectively (Figure S3B).

To identify genes with differential expression during the process of HTSP resistance induction, a linear model (Equation 1) was used to analyse the RNA-Seq data assembled using the CDMC strategy. Significant differences between NHN and N treatments, and between NHN and H treatments were found for 3596 (1788 up-regulated, 1808 down-regulated) and 5379 (2278 up-regulated, 3101 down-regulated), respectively (Table S5). Functions for ca. 20% of the DEGs were unknown based on the Nr database. Among the 8975 DEGs, 1395 were identified in both NHN-N and NHN-H comparisons. Thus, there were 7580 unique DEGs, of which 2201 and 3984 expressed in the N and H treatments, respectively (Figure 3A).

Figure 3

The distribution features of differentially expressed genes (DEGs) under different temperature treatment models. (A) The distribution features of DEGs (I*T) under normal (N), normal-higher-normal (NHN) and higher (H) temperature treatments. (B) Circles-plot of DEGs location in Chinese spring wheat genome. The circles from outside going in represent chromosomes (different colors represent different chromosomes) and each black short line at outside circle represents one DEGs gene. The first outside imaginary line circle is the cluster of false discover rate (FDR) values of DEGs for NHN vs. H (I*T) treatment, the second is the cluster of FDR values of DEGs for NHN vs. N (I*T) treatment. The FDR values of DEGs in inner circles are higher than those in the outer circles. The black dotted line represents threshold (FDR = 0.01).

Chromosome location information was obtained for 3731 of the 7580 unique DEGs. They were located on B chromosomes (41.1%), D chromosomes (29.7%) and A chromosomes (27.2%). Among the B chromosomes, 12.3, 8.7, 8.0, 4.4%, 4.3, 3.0, and 1.9% of the DEGs were on 3B, 5B, 2B, 4B, 1B, 7B, and 6B, respectively (Data S1; Figure 3B).

Verification of RNA-seq analysis by qRT-PCR

To verify the gene expression profiles from the RNA-Seq analysis, 12 transcripts were randomly selected for qRT-PCR analysis. The I-N-0 sample was used as a control when calculating relative expression levels. RNA-Seq (TMM-FPKM) and qRT-PCR results are shown in Figure 4. The relative expression levels of the transcripts from qRT-PCR were nearly identical to those from the RNA-Seq data (Figure 4A), the correlation coefficient was 0.80 (P < 0.0001; Figure 4B).

Figure 4

Verification of RNA-Seq analysis by qRT-PCR. (A) Relative expression levels of 12 randomly selected transcripts are verified through qRT-PCR. The gray histograms represent the relative gene expression levels analyzed using qRT-PCR. The black histograms represent TMM-FPKM values which are the relative expression levels of RNA-Seq data. The error bars are the mean ± SE of three replications. (B) Comparison between the log2 of expression ratios of DGEs obtained from RNA-Seq and qRT-PCR.

Functional annotation and enrichment analysis of DEGs

For the NHN treatment, 1395 DEGs were enriched in KEGG pathways. Notably, DEGs in ribosome metabolism (KO 03010), plant-pathogen interaction (KO 04626) and glycerolipid metabolism (KO 00561) were significantly enriched (Q-value < 0.05; Table 1). Ten DEGs were enriched in the plant-pathogen interaction pathway under the NHN treatment. When the NHN treatment was compared with the N treatment, calcium-dependent protein kinase 1 (CDPK, TCONS_00113300), heat shock cognate protein 80 (Hsp80, TR82962|c0_g1_i1), calcium-binding protein CML31 (CaMCML31, TCONS_00079635) and disease resistance protein RPS2 (RPS2, TR216110|c0_g1_i2) were up-regulated by 11.6, 2.5, 9.3, and 4.0-fold, respectively. When the NHN treatment was compared with the H treatment, CDPK, Hsp80, CML31, and RPS2 were up-regulated by 2.4, 3.7, 5.8, and 1.7-fold, respectively. Those 1395 DEGs under the NHN treatment were enriched based on the GO database (P < 0.05), and included 9, 25, and 20 terms for cellular component, molecular function and biological process categories, respectively (Table S6). GO enrichment results showed that HTSP-related DEGs were mainly involved in membrane proteins, ribonucleoside binding proteins, protein kinase activity, serine family amino acid metabolic processes, phosphotransferase activity, oxylipin metabolism and cell surface receptor signaling pathway processes.

Table 1

Significant KEGG pathways under normal-higher-normal (HNH) treatment for 24 h.

Type	Pathway	DEGs with pathway annotation	All genes with pathway annotation	P-value	Q-value	Pathway ID
NHN (I*T)	Ribosome	39 (17.89%)	2309 (4.38%)	7.50E-66	7.28E-64	ko03010
	Plant-pathogen interaction	10 (4.59%)	619 (1.17%)	0.000391	0.0190	ko04626
	Glycerolipid metabolism	7 (3.21%)	392 (0.74%)	0.001317	0.0329	ko00561

The values of correct-p (Q-value) < 0.05 are considered.

Putative R genes and TFs involved in HTSP resistance to Pst

Sixty-four putative R genes (paralogs and spliceforms) mainly belonging to the RLP (23, eLRR-TM-S/TPK domain), NL (18, NBS-LRR domain), and CNL (9, NB-ARC domain) classes were identified (Table 2). The relative expression levels of 58 putative R genes were higher in the NHN than in the other treatments (Figure 5). Products of nine, six, three, and two putative R gene were homologous to the RPM1, RGA3, RPP13, and RPS2 proteins of A. thaliana, respectively. Twelve putative R gene products were homologous to the LRR receptor-like serine/threonine protein kinase (Ser/Thr PK) of Aegilops tauschii, which belong to RLP. Approximately 21% of those up-regulated putative R genes were located on B chromosomes (Table 2).

Table 2

Summary of predicted HTSP response associated R genes homologous under normal-higher-normal (HNH) treatment for 24 h.

Transcript ID	Type	E-value	Fold change (NHN vs. N)	Fold change (NHN vs. H)	Nr functional annotation	Chromosome location
TCONS_00131852	CNL	0	2.818552	2.703682	Disease resistance protein RPM1 [Aegilops tauschii]	5DL
TCONS_00124527	CNL	3E-120	3.800193	4.137525	Disease resistance protein RPP13 [Triticum urartu]
TCONS_00200075	CNL	0	7.181612	6.031909	Disease resistance protein RGA3 [Brachypodium distachyon]
TCONS_00200311	CNL	0	5.330094	3.185791	Disease resistance protein RGA3 [Brachypodium distachyon]
TR151221_c4_g1_i1	CNL	0	5.715421	5.109299	Disease resistance protein RGA3 [Brachypodium distachyon]
TCONS_00272084	CNL	0	1.029699	0.62448	Disease resistance protein RPM1 [Aegilops tauschii]	7DS
TCONS_00169146	CNL	4E-95	−1.76117	7.233904	Disease resistance protein RGA3 [Aegilops tauschii]
TR216110_c0_g1_i2	CNL	0	4.02049	1.720958	Disease resistance protein RPS2 [Triticum urartu]
TR187245_c0_g1_i1	CNL	2E-38	6.77577	6.946225	Disease resistance protein RPM1 [Triticum urartu]
TR141582_c0_g1_i3	Mlo-like	5E-25	5.27466	4.725616	MLO protein homolog 1	4DL
TR145148_c0_g2_i1	Mlo-like	1E-166	−4.37887	−5.69556	Predicted protein [Hordeum vulgare subsp. vulgare]	2DS
TR101843_c0_g1_i2	Mlo-like	4E-81	5.305871	4.665338	MLO protein-1-like protein [Aegilops tauschii]	4DL
TR202262_c1_g3_i3	NL	0	1.866966	0.266813	Disease resistance protein RPP13 [Aegilops tauschii]	6BL
TCONS_00023192	NL	8E-99	7.252431	5.384763	Disease resistance protein RGA3 [Brachypodium distachyon]
TR200236_c1_g1_i1	NL	0	4.637216	3.008804	Disease resistance protein RGA3 [Brachypodium distachyon]
TR211806_c2_g1_i2	NL	0	−5.03394	−5.63732	Disease resistance protein RGA2 [Triticum urartu]	2AL
TR215734_c0_g2_i1	NL	9E-73	3.415211	2.013215	Predicted protein [Hordeum vulgare subsp. vulgare]	4AS
TCONS_00047832	NL	0	3.178319	3.785673	Disease resistance RPP13-like protein 4 [Aegilops tauschii]
TCONS_00049108	NL	0	6.898312	1.918458	Disease resistance protein RPS2 [Aegilops tauschii]	2DL
TR181960_c2_g1_i7	NL	0	4.130234	1.477969	Resistance protein [Triticum aestivum]	2BL
TCONS_00154991	NL	4E-146	3.296805	3.403045	Disease resistance protein RPM1 [Triticum urartu]	7AS
TCONS_00079413	NL	0	2.55712	1.504546	Disease resistance protein At4g27190-like [Brachypodium distachyon]
TR211589_c1_g1_i5	NL	1E-73	−3.12754	12.24328	Disease resistance protein RPM1 [Aegilops tauschii]	2DS
TR211589_c1_g1_i6	NL	1E-73	−1.36936	9.992023	Disease resistance protein RPM1 [Aegilops tauschii]	2DS
TCONS_00161767	NL	0	4.441952	2.579976	Disease resistance protein RPM1 [Triticum urartu]	7BS
TR230052_c0_g1_i2	NL	0	4.911197	2.643646	Disease resistance protein RPM1 [Triticum urartu]	7BS
TCONS_00003885	NL	0	2.55409	3.007257	Disease resistance protein RPM1 [Triticum urartu]	1AS
TR192331_c2_g4_i3	NL	0	1.397204	1.296923	NBS-LRR disease resistance protein homolog [Hordeum vulgare]	7AS
TCONS_00038730	NL	0	5.263287	4.754892	NBS-LRR disease resistance protein homolog [Hordeum vulgare]	2BS
TCONS_00116383	NL	1E-35	6.941874	4.123755	Predicted protein [Hordeum vulgare subsp. vulgare]	5BL
TCONS_00209942	RLK-GNK2	0	3.438375	4.586322	Cysteine-rich receptor-like protein kinase 10 [Aegilops tauschii]	3AS
TR172572_c3_g1_i2	RLK-GNK2	9E-173	4.135715	7.999334	Cysteine-rich receptor-like protein kinase 10 [Aegilops tauschii]	3AS
TCONS_00188007	RLK-GNK2	0	7.227389	3.927101	Cysteine-rich receptor-like protein kinase 10 [Aegilops tauschii]	1DL
TR151596_c1_g1_i2	RLK-GNK2	2E-143	7.46698	4.99814	Cysteine-rich receptor-like protein kinase 10 [Aegilops tauschii]	1DL
TR151596_c1_g1_i3	RLK-GNK2	0	7.686276	4.382317	Predicted protein [Hordeum vulgare subsp. vulgare]	1DL
TR238545_c0_g1_i2	RLK-GNK2	0	6.931735	3.217689	Cysteine-rich receptor-like protein kinase 25 [Aegilops tauschii]	1DL
TR207091_c2_g1_i1	RLK-GNK2	2E-153	1.643478	5.98336	Predicted protein [Hordeum vulgare subsp. vulgare]
TR238545_c0_g1_i1	RLK-GNK2	0	5.823815	3.253725	Cysteine-rich receptor-like protein kinase 25 [Aegilops tauschii]	1DL
TCONS_00128885	RLP	6E-70	6.362442	4.4303	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]
TCONS_00012063	RLP	0	4.023625	7.780404	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	1BL
TCONS_00022549	RLP	0	3.748471	5.147999	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	1DL
TCONS_00185542	RLP	0	2.675208	6.200027	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]
TCONS_00160106	RLP	0	6.053778	8.606126	Probable LRR receptor-like serine/threonine-protein kinase [Brachypodium distachyon]	7AL
TCONS_00048175	RLP	7E-162	3.243099	4.127629	Unnamed protein product [Triticum aestivum]
TCONS_00199321	RLP	0	4.359545	5.702405	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	2BS
TCONS_00032886	RLP	0	3.010608	3.337234	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	2BS
TR189161_c0_g1_i8	RLP	0	1.329894	2.968987	LRR receptor-like serine/threonine-protein kinase EFR [Aegilops tauschii]	5BL
TCONS_00023422	RLP	0	2.083794	3.047996	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	2AS
TCONS_00018367	RLP	3E-32	2.83137	3.420183	F-box/WD-40 repeat-containing protein [Aegilops tauschii]
TCONS_00153886	RLP	0	2.703369	3.923887	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	6DL
TR237425_c1_g4_i5	RLP	0	−3.33041	5.691968	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	2BL
TCONS_00084905	RLP	0	3.093912	8.190439	Predicted protein [Hordeum vulgare subsp. vulgare]	3DL
TCONS_00197961	RLP	6.00E-165	6.101454	4.94249	Probable inactive receptor kinase At1g27190 [Brachypodium distachyon]	2BL
TR77612_c0_g1_i1	RLP	2.00E-52	4.406603	4.290932	Probable inactive receptor kinase At1g27190 [Brachypodium distachyon]	2AL
TR82468_c0_g1_i1	RLP	7.00E-124	2.351209	3.951302	LRR receptor-like serine/threonine-protein kinase At3g47570 [Oryza brachyantha]	4BL
TCONS_00147950	RLP	0	3.837951	3.667534	Tyrosine-sulfated glycopeptide receptor 1 [Aegilops tauschii]	6DS
TCONS_00031899	RLP	2E-73	7.067972	6.183877	Sulfotransferase 17 [Aegilop stauschii]
TCONS_00161220	RLP	6E-88	3.897112	6.287682	Sulfotransferase 17 [Aegilops tauschii]
TR118780_c0_g1_i3	RLP	6E-38	3.485217	3.081509	Somatic embryogenesis receptor kinase 3, partial [Commiphora wightii]
TCONS_00006438	RLP	0	4.191838	2.755574	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	1AL
TR163704_c0_g1_i5	Other	0	3	2.312676	G-type lectin S-receptor-like serine/threonine-protein kinase SD2-5 [Triticum urartu]
TCONS_00019317	RLP	0	8.087051	3.460008	LRR receptor-like serine/threonine-protein kinase [Aegilops tauschii]	1DL
TCONS_00055825	Other	0	3.562162	2.970888	Zeamatin [Aegilops tauschii]	3AS
TR171566_c1_g1_i1	N	4E-111	1.951928	5.056486	ABC transporter, ATP-binding protein [Galdieria sulphuraria]

Mlo-like, mlo-like resistant proteins; CNL, contains a central nucleotide-binding (NB) subdomain as part of a larger entity called the NB-ARC domain; RLP, receptor like proteins, consist of a LRR-like repeat, a transmembrane region of ~25 AA, and a short cytoplasmic region, with no kinase domain; NL, contains NBS at N-terminal and LRR at C-terminal, and lacks the CC domain; RLK, class with additional domain GNK2; N, contains NBS domain only; Other, has resistance function but does not fit the known classes. The E-values are calculated by plant resistance gene database. Chromosome location information is from support data .

Figure 5

Hierarchical clustering of 64 putative R proteins encoded by differentially expressed genes (DEGs) from the Table 2. The signal ratios are shown in a red-green color scale, where red represents up-regulation and green represents down-regulation. Each column represents the mean expression value (log2TMM-FPKM, 24 h sample are divided by values of 0 h samples) of the RNA-Seq data obtained from three biological replicates. Each row represents a DEG.

There were 227 DEGs that were putatively identified as TFs belonging to different families (WRKY, NAC, and MYB, etc.): 177 were up-regulated and 50 were down-regulated in the NHN treatment when compared with the N treatment (Figure 6A). Among these TFs, WRKY was the family with the largest number of genes, but further domain alignment analysis indicated that only 17 (including alternative splicing) of these genes contained a complete WRKY domain and the zinc finger motif type, including WRKY41, WRKY70, WRKY55, WRKY53, WRKY51, WRKY50, WRKY48, WRKY46, WRKY45 and WRKY15; most of these WRKY genes were up-regulated under the NHN treatment compared with those under the other treatments. Among the WRKY genes, the relative expression level of WRKY 41 was the highest; it had five homologs (paralogs or spliceforms) under the NHN treatment (Figure 6B).

Figure 6

Identification of differentially expressed genes (DEGs) by plant transcription factor database. (A) Identification of the distribution feature of differentially expressed transcription factors under the NHN treatment. Red bars represent up-regulated genes and green bars represent down-regulated genes. (B) Heatmap of the confirmed WRKY transcription factors (containing complete WRKY domain) responding to HTSP from the NHN treatment.

Protein interaction network in the HTSP response to Pst

A total of 135 proteins, mostly consisting of lignin and fatty acid synthesis related proteins, ribosome proteins, protein kinases, heat shock proteins, WRKY TFs, and R proteins homologous, were involved in the main network based on the STRING database of A. thaliana (Data S2; Figure 7). In particular, compared to other proteins, phosphatase 2C10 (PP2C10, TCONS_00197067) was the hub protein with the highest degree of 52 (Data S2); this protein was predicted to interact with 35 protein kinases. In addition, most of R proteins interacted with Hsp80. For example, Hsp80 interacted with Ser/Thr PK and RPS2. Additionally, lignin and fatty acid biosynthesis related proteins formed a separate interaction network. For instance, 4-coumarate ligase 2 (4CL2, TR211614|c0_g1_i10), elicitor-activated gene 3-1 (ELI3-1, TR172862|c0_g1_i2), and elicitor-activated gene 3-2 (ELI3-2, TR171566|c1_g1_i1) were associated with lignin biosynthesis. When the NHN treatment was compared with both the N and H treatments, the majority of the ribonucleoproteins (RPs) encoded by DEGs interacted with each other and all these DEGs were down-regulated.

Figure 7

Protein interaction network in higher-temperature seedling-plant (HTSP) responding to Puccinia striiformis f. sp. tritici under the normal-higher-normal (NHN) temperature treatment. Different colors represent fold changes of differentially expressed genes (DEGs) (NHN vs. N). Each node is a DEG. Each size of the node represents a false discover rate (FDR) value and the smaller node was more significant than a bigger one. The FDR values ranged from 1.39E-59 to 0.049924. Protein interaction network was constructed by CYTOSCAPE software.

Discussion

Identification of differentially expressed genes under the NHN treatment involved in HTSP resistance to Pst

Studies on the inheritance and ultrastructural analysis of the HTSP resistance in XY 6 have been reported previously (Wang and Shang, 2003; Yao et al., 2006). These studies showed that constant relatively high temperature cannot induce resistance to Pst in XY 6 (Shang and Wang, 1997) and that only temperature changes can activate the resistance. Similarly, HTAP resistance to Pst becomes effective when wheat plants were grown under a night/day cycle of 10~12°C/25~30°C after inoculation (Chen, 2005, 2013; Coram et al., 2008a; Bryant et al., 2014). In fact, wheat crops grown in the fields are almost all under fluctuating temperature conditions every day and changing temperatures throughout the growth season, and therefore, screening germplasm for stripe rust resistance should be carried out at diurnally changed temperatures (Chen, 2013). Wheat plants at 18°C facilitated resistance to the Triticum mosaic virus (TriMV), which was controlled by temperature dependent Wsm1 and Wsm2 genes (Tatineni et al., 2016). The present study confirmed that the HTSP resistance in XY 6 was activated by exposure to 20°C for 24 h. Also, the expression level trends of 1395 DEGs in both the H and N treatments were similar, consistent with the histological data on Pst development. These results showed that sudden changes in temperature play an active role in the defense response to Pst in XY 6. Coram et al. (2008a) identified 99 transcripts involved in the Yr39-mediated HTAP resistance to stripe rust, including R protein homologs, pathogenesis-related (PR) proteins, protein kinases and phenylpropanoid biosynthesis. In contrast, we identified a total of 1395 DEGs in the HTSP resistance induced by the NHN treatment, including genes coding phosphatase 2C10, protein kinases, R protein homologs, TFs and RPs as well function unknown proteins were specific resistance. The high number of DEGs was achieved by taking the advantage of RNA-Seq over the previously used microarray technique (Coram et al., 2008a).

Phosphatase 2C proteins may play a positive role in HTSP resistance to Pst

Phosphatase 2C (PP2C) has been reported to be involved in the regulation of plant development and the adaptation to environmental stresses (Schweighofer et al., 2004; Bhatnagar et al., 2017). Recently, there has been an increasing focus on the role of PP2C in plant stress signaling: cold, drought, high salt, etc. (Hu et al., 2010; Liu et al., 2012; Arshad and Mattsson, 2014), indicating that PP2C could receive signals rapidly under abiotic stresses. However, reports about the function of PP2C in disease resistance are limited. In an alien substitution line of Triticum aestivum-Elytrigia elongatum, a phosphateserine aminotransferase was shown to be involved in defense to powdery mildew (He and Wang, 2005). Protein kinases such as Ser/Thr PK can alter the functions of proteins by phosphorylating the OH group of serine or threonine residues, and protein phosphorylation plays an important role in disease resistance (Cao et al., 2011). In the present study, PP2C10 in XY 6 was up-regulated and predicted to directly interact with 35 protein kinases during the induction process of HTSP resistance based on the analysis of the database of Arabidopsis protein interactions. This suggests that PP2C10 could be a hub protein and may play a pivotal role, such as signal switch in HTSP resistance against Pst by regulating the activity of protein kinases.

Chromosomal locations and predicted functions of R genes

Wang and Chen (2017) summarized a total of 451 genes and QTL with chromosomal locations for resistance to stripe rust in wheat identified through molecular mapping, of which 49% are on B chromosomes while only 31% on A chromosomes and 20% on D chromosomes, indicating that the B genome is more involved in stripe rust resistance than either A or D genomes. In the present study, we found 41% of the 3731 DEGs with mapped chromosomal locations were on B chromosomes, more than either A chromosomes (27.2%) or D chromosomes (29.7%). Furthermore, 21% of the 64 R genes were on B chromosomes. Our results also suggest that the B genome of wheat contains more genes for defense to Pst.

There have been many studies on temperature sensitive R genes against stripe rust, especially non-race-specific HTAP resistance, such as Yr36 (Fu et al., 2009; Bryant et al., 2014); Yr52 (Ren et al., 2012); Yr59 (Zhou et al., 2014), Yr62 (Lu et al., 2014), LrZH22 (Wang et al., 2016), and Yr79 (Feng et al., 2018). Previously cloned non-race-specific resistance genes to stripe rust and/or leaf rust, such as Yr18/Lr34 (Krattinger et al., 2009a), Yr36 (Fu et al., 2009), and Yr46/Lr67 (Moore et al., 2015), do not contain NBS-LRR domains (Chen, 2013), in contrast to many race-specific resistance to stripe rust and leaf rust, such as Lr10 (Feuillet et al., 2003), Lr21 (Huang et al., 2003); Lr1 (Cloutier et al., 2007); Yr10 (Liu W. et al., 2014). Base on the association of the types of genes to the types of resistance, temperature-sensitive and non-NBS-LRR genes have been connected to the non-race specificity and therefore durability of stripe rust resistance (Chen, 2013; Chen et al., 2013; Wang and Chen, 2017). Although we did not study genes in XY 6 for controlling the HTSP resistance to stripe rust as in those studies mentioned above, we identified 23 putative R genes encoding an eLRR-TM-S/TPK domain and most of them were up-regulated during the HTSP induction process. The finding of the numerous R genes homologs involved in HTSP resistance is similar to nine R genes contributing to the Yr39-controlled HTAP resistance (Coram et al., 2008a). The LRR domain has been implicated in PPIs (Martin et al., 2003). PPIs related to disease resistance can be complex. For example, the RPM1 protein of A. thaliana contains NBS and LRR domains, and RPM1-interacting protein 4 (RIN4), a membrane receptor protein, recognizes the avrRpt2 effector (Axtell and Staskawicz, 2003). The cleavage of RIN4 results in RPS2-mediated elicitor-triggered immunity, which can be induced by relatively high temperature (Wang et al., 2009; Zhu et al., 2010). It will be interesting to further study the putative R genes identified in the present study to determine how they interact to each other and to other genes contributing to the HTSP resistance.

Hsp80, which is involved in a variety of regulatory and defense responses, has a molecular chaperone function in which the protein interacts with different domains of R proteins (Takahashi et al., 2003; Maimbo et al., 2007; Wang et al., 2011). The present study predicts that Hsp80 interacts with different domains of R proteins based on the analysis of the database of Arabidopsis protein interactions, such as the NBS-LRR (RPM1, RPS2) and eLRR-TM-S/TPK (Serine threonine-protein kinase) domains. Moreover, the Pto and PBS1 proteins are members of the R protein family containing a S/TPK domain and have been shown to require the NBS-LRR-domain R proteins Prf and RPS5 respectively for their resistance functions (Brueggeman et al., 2008). Based on the up-regulation of Hsp80 and putative R genes during the HTSP induction, we hypothesize that the eLRR-TM-S/TPK domains of the putative R proteins could interact with the NBS-LRR-domain proteins via the Hsp80 protein, leading to the HTSP resistance in XY 6 against Pst.

CDPK and ribosomal proteins are associated with important Ca2+ signaling components involved in HTSP

The CDPK gene was highly up-regulated in the induction process of HTSP. CDPK is an important protein associated with Ca2+ signaling components in immune and stress signaling networks (Bolton, 2009). Previous reports have shown that calcium-dependent CDPK4 and CDPK5 regulate ROS production by phosphorylating NADPH oxidase in potato (Kobayashi et al., 2007). However, ROS is important not only for signaling mechanisms for defense (Eckardt, 2017) but also for regulating programmed cell death through the establishment of the HR (Tamas et al., 2010). Therefore, the CDPK protein identified in the present study may function in ROS accumulation and cell death under HTSP resistance against Pst. In addition, a ribosome translocon complex mediates calcium leakage from endoplasmic reticulum stores, which regulate many physiological processes, including apoptosis (Coppenolle et al., 2004; Garcia et al., 2017). Hence, we suggest that during the initial stage of higher temperature treatment, the down-regulation of ribosomal genes in XY 6 may be part of an emergency reaction against temperature stress. This process could change the composition of ribosomes (Wang et al., 2013), involved in the regulation of calcium leakage, and increase resistance to Pst in XY 6.

WRKY TFs are positively involved in the HTSP resistance of XY 6 to Pst

WRKY TFs belong to a large gene family, are regulated by MAPKs and mediate plant defense responses (Chen et al., 2012). A large number of studies have shown that most of WRKY TFs are involved in the salicylic acid (SA) signaling pathway in defense responses (Dong et al., 2003; Xing et al., 2008; Hu Y. et al., 2012; Shimono et al., 2012; Wang et al., 2017). A putative WRKY5 gene was found to be up-regulated in the HTAP resistance controlled by Yr39 (Coram et al., 2008a). In the present study, 17 WRKY genes, including WRKY41 and WRKY70, were identified to be up-regulated during the induction process of HTSP resistance. Overexpression of WRKY41 and WRKY70 leads to the constitutive expression of PR5 and PR1 genes in A. thaliana, which increased resistance to P. syringae pv. syringae and Erysiphe cichoracearum, respectively (Li et al., 2006; Higashi et al., 2008). Silencing TaWRKY70 leads to enhanced susceptibility to Pst when subjected to higher temperature during the initial Pst incubation stage (Wang et al., 2017). Moreover, WRKY proteins interact with not only PR proteins and receptor-like kinases but also other members of the WRKY family in disease resistance (Yu et al., 2001; Peng et al., 2008; Hu Y. et al., 2012). Based on the criterion that PPI confidence scores were greater than 0.7, the interactions between WRKY15 (FDR < 0.02) and WRKY70 (FDR < 0.00002), WRKY15 (FDR < 0.02) and WRKY53 (FDR < 0.03) were identified. These results indicate that WRKY TFs may work together with other members of WRKY proteins and play an important role in the crosstalk of SA signaling pathways in HTSP resistance to Pst.

DEGs are involved in signaling pathways during the HTSP response of XY 6 to Pst

A complex network of signaling pathways induced by phytohormones such as SA regulates local and systemic resistance to invasive pathogens (Panstruga et al., 2009). We identified several major signal transduction pathways that are likely involved in the HTSP response of XY 6 to Pst (Figure 8). First, HTSP resistance to Pst is induced by temperature changes and Pst inoculation, which may induce changes in the phosphorylation status of Ser/Thr PK via PP2C10 or membrane associated proteins such as RIN4. These changes could then lead to R genes directly or indirectly recognizing Avr elicitors of Pst by interacting with Ser/Thr PK or RIN4, activating the HR of XY 6 with the help of Hsp proteins. Second, R gene-mediated resistance pathways may activate SA signaling pathways, which regulate downstream MAPK proteins and transcription factors (WRKY); this process results in a defense response in XY 6. Third, phosphorylation of membrane-associated proteins may activate H+-ATPases and act as an important primary sensor of Ca2+ leakage, leading to the activation of CDPK or RPs. This process would subsequently regulate ROS production by phosphorylating NADPH oxidase, leading to the HR in XY 6.

Figure 8

A summary of the molecular pathways and cellular processes in higher-temperature seedling-plant (HTSP) resistance to Puccinia striiformis f. sp. tritici in XY 6. The red dotted lines indicate the proposed pathways and the green lines denote the protein-protein interactions based on the KEGG and PPIs database of Arabidopsis thaliana. A question mark indicates the interaction of proteins that needs functional verification.

In summary, HTSP resistance to Pst in XY 6 is induced by the exposure to 20°C for 24 h during the early Pst incubation stage and may be controlled by several major genes in conjunction with other minor-effect genes. Functions of many identified DEGs remain unknown although some candidate genes were identified based on the current NCBI database. For example, PP2C10 and LRR receptor-like serine/threonine protein kinases may play important roles in the processes of SA and Ca2+ signal transduction in HTSP resistance to Pst. Identified DEGs were located on A chromosomes (29.2%), B chromosomes (41.1%), and D chromosomes (29.7%); most of the defense related DEGs were located on the B chromosome group. Together, these results constitute a strong base for future research on HTSP resistance to Pst in wheat.

Author contributions

XH, XX, and JY planned and designed the research. FT, JW, and ZG performed RNA-Seq experiments. FT and JH analyzed the data. FT, XH, XC, and XX wrote the manuscript.

Conflict of interest statement

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.