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 cm³) 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.

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.

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).

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.

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_i 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; α, β_I, and β_T are the batch effect, inoculation effect and temperature effect, respectively; IT is the interaction term for inoculation and temperature; β_IT is the IT interaction effect; and ε_i is the random error following a normal distribution with mean of 0 and variance of σ². 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 log₂-fold change > 1 or < −1 and a logCPM (log₂ counts per million) > −2 were used as additional criteria to select DEGs (McCarthy et al., 2012; Liu Y. et al., 2014).

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).

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).

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).

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).

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.

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.

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).

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).

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).

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.

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).

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).

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.

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 (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.

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.

The CDPK gene was highly up-regulated in the induction process of HTSP. CDPK is an important protein associated with Ca²⁺ 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 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.

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 Ca²⁺ 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.

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 Ca²⁺ 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.

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.