Identification of Candidate Susceptibility Genes to Puccinia graminis f. sp. tritici in Wheat

Wheat stem rust disease caused by Puccinia graminis f. sp. tritici (Pgt) is a global threat to wheat production. Fast evolving populations of Pgt limit the efficacy of plant genetic resistance and constrain disease management strategies. Understanding molecular mechanisms that lead to rust infection and disease susceptibility could deliver novel strategies to deploy crop resistance through genetic loss of disease susceptibility. We used comparative transcriptome-based and orthology-guided approaches to characterize gene expression changes associated with Pgt infection in susceptible and resistant Triticum aestivum genotypes as well as the non-host Brachypodium distachyon. We targeted our analysis to genes with differential expression in T. aestivum and genes suppressed or not affected in B. distachyon and report several processes potentially linked to susceptibility to Pgt, such as cell death suppression and impairment of photosynthesis. We complemented our approach with a gene co-expression network analysis to identify wheat targets to deliver resistance to Pgt through removal or modification of putative susceptibility genes.


INTRODUCTION
Stem rust caused by Puccinia graminis f. sp. tritici (Pgt) is one of the most devastating foliar diseases of wheat (Triticum aestivum) and barley (Hordeum vulgare). The economic relevance of this pathogen to food security is demonstrated by the impact of historical and recent epidemics (Pretorius et al., 2000;Peterson, 2001;Olivera et al., 2015;Singh et al., 2015;Bhattacharya, 2017;Steffenson et al., 2017). Consistent with its biotrophic lifestyle, Pgt develops an intricate relationship with its host in order to acquire nutrients and survive. Early stages of infection involve the germination of urediniospores (asexual spores) and host penetration through the formation of appressoria over stomata (Staples and Macko, 1984). As the fungus reaches the mesophyll cavity of the plant, it develops infection hyphae which penetrate plant cell walls and differentiate into specialized feeding structures, known as haustoria. Haustorial development takes place during the first 24 h post-infection and is critical for colony establishment and sporulation that re-initiates the infection cycle (Harder and Chong, 1984). Similar to other plant pathogens, cereal rust infections involve the translocation of effectors to the plant cell as a mechanism to shut down basal defenses activated by PAMP triggered immunity (PTI) and manipulate host metabolism (Dodds and Rathjen, 2010;Couto and Zipfel, 2016). In rust fungi, the haustorium mediates the secretion of effectors, although the underlying molecular mechanism that facilitates this process is not known (Garnica et al., 2014;Petre et al., 2014). The plant targets of effectors and other plant genes that mediate compatibility and facilitate pathogen infection are often regarded as susceptibility (S) genes (Lapin and van den Ackerveken, 2013;van Schie and Takken, 2014;Lo Presti et al., 2015).
To avoid infection by adapted pathogens, plants employ effector-triggered immunity (ETI) which is mediated by the recognition of effectors by nucleotide-binding domain leucinerich repeat (NLR) receptors (Flor, 1971;Dodds and Rathjen, 2010;Garnica et al., 2014;Petre et al., 2014). These specific recognition events often induce localized cell death at infection sites (hypersensitive response, HR) which restrict pathogen growth. In wheat-rust interactions, ETI is manifested by the reduction or absence of fungal growth and sporulation (Periyannan et al., 2017). The use of NLR genes to provide crop protection was a critical component of the Green Revolution which diminished the impact of stem rust epidemics (Ellis et al., 2014). While this approach still contributes to the development of wheat cultivars with genetic resistance to stem rust, the durability of such resistant cultivars is hampered by the evolution of rust populations to avoid recognition by NLRs. Given the economic and environmental advantages of genetic disease control over chemical applications, the identification of alternative genetic sources of resistance are a priority for securing future wheat production. In this context, the discovery of S genes could have important translational applications for agriculture and potential durable disease control. Mutations in S genes, although often recessive, could shift a genotype to a non-suitable host due to alterations in initial recognition stages or loss of pathogen establishment requirements (van Schie and Takken, 2014;Lo Presti et al., 2015).
The genetic factors that contribute to wheat susceptibility to biotrophic pathogens such as rust fungi remain largely unknown. Numerous structural and physiological alterations have been observed in wheat-rust compatible interactions. At early infection stages, 4-6 days post-inoculation (dpi), the cytoplasm of infected mesophyll cells increases in volume and an extensive network of the endoplasmic reticulum is built near the haustorium (Bushnell, 1984). The nucleus of infected cells also increases in size and migrates toward the haustorium, and in some cases both structures appear in proximity. These observations suggest that plant cells undergo a massive transcriptional reprogramming to either accommodate rust colonization or initiate a cascade of plant defenses to prevent infection. In addition, many biotrophs are known to increase the ploidy of host cell nuclei near infection sites (Wildermuth, 2010). Advances in next generation sequencing and data mining bring new opportunities to deepen our understanding of plant-pathogen interactions and the relationship between plant metabolism and disease resistance or susceptibility. Several transcriptome profiling studies comparing compatible and incompatible wheat-rust interactions provide strong evidence for the complexity of these interactions (Bozkurt et al., 2010;Zhang et al., 2014;Chandra et al., 2016;Dobon et al., 2016;Yadav et al., 2016). Although S genes in rust pathosystems are largely unknown, several susceptibility factors to other plant pathogenic fungi have been identified in Arabidopsis thaliana, H. vulgare (barley), and solanaceous plants (van Schie and Takken, 2014;Zaidi et al., 2018).
To expand our knowledge of wheat-rust interactions and identify candidate S genes to direct future functional studies, we conducted a comparative RNA-seq analysis of the molecular responses to Pgt in compatible and incompatible interactions. We included a susceptible genotype (W2691) of Triticum aestivum (bread wheat) and the same genotype containing the resistance gene Sr9b, which confers race-specific responses to various Pgt isolates (McIntosh et al., 1995). We also included the related grass species Brachypodium distachyon, which is recognized as a nonhost to various cereal rust species (Kellogg, 2001;Ayliffe et al., 2013;Figueroa et al., 2013Figueroa et al., , 2015Bettgenhaeuser et al., 2014Bettgenhaeuser et al., , 2018Gilbert et al., 2018;Omidvar et al., 2018). As part of our analysis, we examined the expression profiles of T. aestivum and B. distachyon orthologs of several known S genes in Arabidopsis thaliana, H. vulgare, as well as other characterized S genes in T. aestivum, and identified groups of genes co-regulated with these S gene candidates. In conclusion, this study provides an overview of global expression changes associated with failure or progression of Pgt infection in T. aestivum and B. distachyon and insights into the molecular processes that define disease incompatibility.
Pgt Infection of T. aestivum and B. distachyon Genotypes Brachypodium distachyon seeds were placed in petri dishes with wet grade 413 filter paper (VWR International) at 4 • C for 5 days and germinated at room temperature for 3 days before sowing to synchronize growth with wheat plants which did not require stratification. Seeds of both wheat and B. distachyon were sown in Fafard R Germination Mix soil (Sun Gro Horticulture, Agawam, MA, United States). All plants were grown in growth chambers with a 18/6 h light/dark cycle at 21/18 • C light/dark and 50% relative humidity. Urediniospores of Pgt were activated by heatshock treatment at 45 • C for 15 min and suspended in Isopar M oil (ExxonMobil) at 10 mg/mL concentration. Inoculation treatments consisted of 50 µl of spore suspension per plant, whereas mock treatments consisted of 50 µl of oil per plant. Fungal and mock inoculations were conducted on 7-day old wheat plants (first-leaf stage) and 12-day old B. distachyon plants (three-leaf stage). After inoculations, plants were kept for 12 h in mist chambers with repeated misting for 2 min every 30 min and returned to growth chambers under the previously described conditions.

Analysis of Fungal Colonization and Growth
At 2, 4, and 6 dpi T. aestivum and B. distachyon leaves were sampled and cut into 1 cm sections before staining with Wheat Germ Agglutinin Alex Fluor R 488 conjugate (WGA-FITC; ThermoFisher Scientific) following previously described procedures (Omidvar et al., 2018). Time points to represent stages of Pgt infection were selected based on previous characterization (Figueroa et al., 2013(Figueroa et al., , 2015. To determine the level of fungal colonization, the percentage of urediniospores that germinated (GS), formed an appressorium (AP), established a colony (C), and differentiated a sporulating colony (SC) were visualized using a fluorescence microscope (Leica model DMLB; 450-490 nM excitation). The progression of fungal growth was recorded for 100 infection sites for each of the three biological replicates. Genomic DNA was extracted from T. aestivum (three infected primary leaves) and B. distachyon (three infected secondary leaves) using the DNeasy Plant Mini Kit (Qiagen) and were standardized to a 10 ng/µl concentration. The ITS regions were amplified by qPCR using ITS-specific primers provided by the Femto TM Fungal DNA Quantification Kit (Zymo Research) to quantify the relative abundance of fungal DNA following the manufacturer's recommendations for the three biological replicates. The GAPDH housekeeping gene from each species was used as an internal control to normalize fungal DNA quantities (Omidvar et al., 2018).

RNA Isolation, Purification, and Sequencing
Infected and mock treated primary leaves from W2691 and W2691+Sr9b and secondary leaves from Bd21-3 were collected at 2, 4, and 6 dpi. For each of the three biological replicates, three infected leaves were pooled for RNA extraction using the RNeasy Plant Mini Kit (Qiagen). Subsequently, stranded-RNA libraries were constructed, and 125 bp paired-end reads were sequenced on an Illumina HiSeq TM 2500 instrument at the University of Minnesota Genomics Center. On average, more than 10 million reads were generated per time point in each of the previously listed plant-rust interactions (Supplementary Table 1).

Expression Profiling and Identification of Differentially Expressed Genes
Reads were mapped to T. aestivum and B. distachyon gene features using htseq v.0.11.0 to obtain count values (Anders et al., 2015). Normalized read counts and differential expression (DE) analysis were performed with DESeq2 v1.28.1 (Love et al., 2014). Genes with a | log2 fold change| ≥1.5 and a p-value < 0.05 were identified as differentially expressed genes (DEGs).

Gene Ontology Analysis
Gene ontology (GO) terms were obtained from GOMAP track data for T. aestivum (Alaux et al., 2018) and previously published data for B. distachyon (Brachypodium distachyon Bd21-3 v1.2 DOE-JGI, http://phytozome.jgi.doe.gov/) annotation files. GO terms in wheat and B distachyon were mapped to the GOslim plant subset using OWLTools with the command owltools -map2slim 1 . GO enrichment analysis for DEGs was performed using the topGO R package using the "weight01" algorithm and fisher test statistic (Alexa and Rahnenfuhrer, 2020). Enriched terms were considered significant with a Fisher test p-value < 0.01 (Supplementary Table 2). Enrichment analyses using the GOslim subset were performed on all differentially expressed wheat and B. distachyon genes, as well as on genes within the S-gene orthologs clusters. Enrichment analysis with the full GO set was only performed on the differentially expressed T. aestivum and B. distachyon genes using the same methods described above.

Orthology Analysis
Protein sequences from S genes of interest (Supplementary Table 3) as cited in original publications as reviewed by van Schie and Takken (2014) were cross-checked using gene name and synonym information and the Basic Local Alignment Search Tool (BLAST) functions in the TAIR gene search database 2 , EnsemblPlants 3 , UniPro 4 , and the IPK blast server 5 . OrthoFinder version 2.4.0 (Emms and Kelly, 2019) was used to identify orthologs between A. thaliana 6 , H. vulgare 7 , T. aestivum (annotation version 1.1 https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_ Annotations/v1.1/iwgsc_refseqv1.1_genes_2017July06.zip), and B. distachyon (annotation version 1.2 https://phytozomenext.jgi.doe.gov/info/BdistachyonBd21_3_v1_2) proteins. For genes in A. thaliana, H. vulgare, and T. aestivum with multiple isoforms, perl scripts for each species were used to retain only the longest representative transcript for use in the orthology analysis 8 . The longest transcript file for B. distachyon (BdistachyonBd21_3_537_v1.2.protein_primaryTranscriptOnly. fa) was obtained from Phytozome. The default settings of OrthoFinder were used, and orthologs of the four species were obtained in a single run. The URGI BLAST tool 9 was used to identify candidates for missing subgenome representatives.

Protein Sequence Phylogenetic Analysis
Using the longest protein sequence from known S genes in A. thaliana (Lamesch et al., 2011) and H. vulgare (Howe et al., 2019), as well as the longest protein sequences from the orthologous candidate S genes in T. aestivum and B. distachyon (Supplementary Table 3), phylogenetic trees were constructed to examine the relationship of ortholog families using the web-based tool NGPhylogeny (Lemoine et al., 2019). Default parameters for the FastME one-click workflow were used for MAFFT alignment, BMGE curation, and FASTME tree inference 10 . A R script using the packages ggplot2, ggtree, and ape was used to generate visualizations of the generated phylogenetic trees (Wickham, 2016;Yu et al., 2017;Paradis and Schliep, 2018).

Gene Co-expression Network Analysis
Individual gene co-expression networks (GCNs) were constructed and analyzed for T. aestivum W2691, W2691+Sr9b, and B. distachyon Bd21-3 genotypes using the python package Camoco (Schaefer et al., 2018). To build each network, all three independent RNA-seq replicates from all three time points (2, 4, and 6 dpi) of infected and mock-inoculated treatments were used. HTSeq read counts were converted to FPKM values for Camoco compatibility, and subjected to inverse hyperbolic sine transformation normalized against median FPKMs across all samples. Genes with coefficient of variation <0.1 across all samples or without a single sample having an expression above 0.5 FPKM were removed from analysis. Additionally, genes with a FPKM value > 0.001 across 60% of samples were included in network analyses. Pearson correlation metrics between all gene pairs were calculated and subjected to Fisher transformation to generate Z-scores with a cutoff of Z > = 3 to allow comparisons between networks (Huttenhower et al., 2006). Finally, correlation metrics were used to build weighted GCNs. Clusters containing susceptibility gene orthologs were visualized using ggplot2 (Wickham, 2016), ggnetwork (Briatte, 2020), sna (Butts, 2019), and network (Butts, 2015) R packages.

Data Availability
Sequence data was deposited in NCBI under BioProject PRJNA483957 (Supplementary Table 1). Unless specified otherwise, Supplementary Tables, scripts and files for analysis and visualizations are available at https://github.com/henni164/ stem_rust_susceptibility. Due to large file sizes Supplementary Tables 8, 9 can only be found the github page mentioned above.

T. aestivum and B. distachyon Differ in Susceptibility to Pgt
We compared the infection and colonization of Pgt in two T. aestivum isogenic lines that were susceptible (W2691) or resistant (W2691+Sr9b) to Pgt as well as in the non-host B. distachyon Bd21-3. The Bd21-3 line was selected to ease future adoption of reverse genetics approaches as an extensive collection of T-DNA insertional mutants is available in the same background (Bragg et al., 2012). Symptom development upon infection was consistent with previous observations reporting susceptibility of W2691 and W2691+Sr9b mediated resistance (intermediate) to race (pathotype) SCCL (Figures 1A,B; Zambino et al., 2000). Susceptibility was manifested by formation of large sporulating pustules in W2691, while small pustules surrounded by a chlorotic halo were characteristic of Sr9b mediated-resistance at 6 days post-inoculation (dpi). Susceptibility differences between W2691 and W2691+Sr9b were evident at 6 dpi as formation of fungal colonies was present in both genotypes, but colony sizes were larger in W2691 than W2691+Sr9b (Figure 1). B. distachyon supports the formation of colonies that are smaller than those in the resistant T. aestivum line W2691+Sr9b with no visible macroscopic symptoms observed at 6 dpi ( Figure 1C). To monitor the progression of fungal growth and colonization, we quantified the percentage of GS, and interaction sites displaying the formation of AP, colony formation (C), and colony sporulation (SC) at 2, 4, and 6 dpi using microscopy ( Figure 1D). The germination frequency (∼95%) was similar between all three genotypes tested (ANOVA test, p > 0.05). The percentage of interaction sites showing appressorium formation (AP) was higher in wheat than in B. distachyon at 4 dpi (ANOVA test, p ≤ 0.035). The genotype W2691 displayed the highest percentage of interaction sites showing colony formation at 4 and 6 dpi (ANOVA test, p ≤ 0.002), and sporulation at 6 dpi (ANOVA test, p ≤ 0.0015). In contrast, a smaller number of rust colonies formed in B. distachyon, and these colonies did not show signs of sporulation. To estimate rust colonization levels on T. aestivum and B. distachyon, we quantified the abundance of fungal DNA in infected leaves at 2, 4, and 6 dpi ( Figure 1E). Fungal DNA abundance (rust colonization) among all genotypes was not significantly different at 2 dpi (ANOVA test, p > 0.05); however, there was a trend at 4 and 6 dpi for higher rust colonization in W2691 than in W2691+Sr9b and B. distachyon (ANOVA test, p > 0.05).

Putative Biological Processes Associated With in planta Responses to Pgt
The transcriptome profiles of T. aestivum (W2691 and W2691+Sr9b) and B. distachyon (Bd21-3) in response to Pgt infection at 2, 4, and 6 dpi were examined using RNA-seq expression profiling (Supplementary Table 1). DE analysis was used to compare responses to rust infection relative to the baseline mock treatments. Overall, the number of DEGs increased in W2691, W2691+Sr9b, and Bd21-3 over the course of infection (Table 1). Between 11 and 12.9% of T. aestivum genes were differentially expressed at 6 dpi, whereas in Bd21-3 only 6.2% were differentially expressed. We conducted a GOslim enrichment analysis on up-and down-regulated DEGs for each interaction at the infection time points (Figure 2). At 2 dpi, W2691 and W2691+Sr9b had only a few GOslim terms enriched in either up-or down-regulated DEGs. At 4 dpi, greater similarities between the T. aestivum genotypes emerged with very similar enrichment patterns in GOslim terms. The similarity of GOslim term enrichment continued at 6 dpi, with W2691 and W2691+Sr9b having nearly identical enrichment patterns. W2691+Sr9b had one additional term enriched in both up-regulated (cytoplasm, GO:0005737) and down-regulated (chromatin binding, GO:0003682) genes. Compared to the two T. aestivum genotypes, Bd21-3 had fewer terms enriched across all three timepoints and only a few terms were in common with W2691 and W2691+Sr9b, i.e., extracellular region (GO:0005576), DNA-binding transcription factor activity (GO:0003700). Bd21-3 had several unique terms in both upand down-regulated categories, among them mitochondrion (GO:0005739), transporter activity (GO:0005215), catalytic activity (GO:0003824), and DNA binding (GO:0003677) were up-regulated, while intracellular (GO:0005622), DNAbinding transcription factor activity (GO:0003700), catalytic activity (GO:0003824), and DNA binding (GO:0003677) were down-regulated. The full GO set also demonstrated clear differences between the T. aestivum genotypes and Bd21-3. Photosynthesis-related terms such as chloroplast photosystem I and II (GO:0030093 and GO:0030095), photosystem II antenna complex (GO:0009783), and PSII associated light-harvesting complex II (GO:0009517) were overrepresented at 4 and 6 dpi in W2691 and W2691+Sr9b, but not in Bd21-3 (Supplementary Table 2). In addition, Bd21-3 only had enrichment in 11 terms across the cellular component (CC), biological process (BP), and molecular function (MF) categories compared to the terms enriched 741 across the three categories in W2691 and W2691+Sr9b (Supplementary Table 2). Overall, this analysis highlights how the molecular and genetic responses of Bd21-3  Chen et al., 2007Chen et al., , 2010Low et al., 2020), and this knowledge has allowed us to further understand molecular plantmicrobe interactions. With an interest in identifying potential S genes in T. aestivum as well as creating resources to enable future studies, we designed an experimental workflow based on the identification of known S gene orthologs, gene expression comparisons and co-expression network analysis (Figure 3). A curated set of previously characterized or postulated S genes as summarized by van Schie and Takken (2014) Table 4). These genes from T. aestivum  and B. distachyon were selected as S gene orthologs. A total of 29,767 genes (orthogroup OG0029421 to OG0059187) were assigned groups with only one member (singleton orthogroups; Supplementary Table 5).
The gene expression patterns of S gene orthologs in T. aestivum and B. distachyon were used to identify which orthologs may act as susceptibility factors (Figure 3 and Supplementary Table 6). The selection criterion was applied to include DEGs that showed a progressive increase in log2 fold change (mock vs infected, | log2 fold change| ≥1.5 and a p-value < 0.05) in W2691 or in both W2691 and W2691+Sr9b, but the corresponding orthologs in B. distachyon and/or W2691+Sr9b showed a decrease or no change, as observed in various systems Pessina et al., 2014). The assumption is that S genes will be up-regulated during infection when the pathogen reaches the sporulation stage (e.g., in a susceptible or intermediate resistant host represented by W2691 and W2691+Sr9b, respectively) but with a low or no regulatory change in a non-host (Bd21-3). Expression data for all genes can be found in Supplementary Table 6 in association with orthogroup number. Most genes in the 70 orthogroups did not demonstrate major changes in expression over the course of the experiment (Supplementary Figure 1), including the orthogroup OG0001703, which contains the Mlo (Mildew locus O) alleles and orthologous sequences. Eight orthogroups that demonstrated these expression patterns were chosen for further analysis; these included ortholog genes for AGD2 (aberrant growth and death 2), BI-1 (BAX inhibitor-1), DMR6 (downy mildew resistance 6), DND1 (defense, no death), FAH1 (fatty acid hydroxylase 1), IBR3 (IBA response 3), VAD1 (vascular associated death 1), and WRKY25 (WRKY DNA binding protein 25; Figure 4, Table 2 and Supplementary Table 7). Among the eight susceptibility orthogroups, T. aestivum orthologs of BI-1, DMR6, and WRKY25 showed the greatest increase in fold change (Supplementary Table 7) in either W2691 or W2691+Sr9b, particularly at 6 dpi (Figure 4). The gene ortholog of DND1 displayed a higher fold change in W2691 than in W2691+Sr9b.
The phylogenetic relationships of the orthogroups to known S genes were confirmed using NGphylogeny (Supplementary  Figure 2). A phylogenetic tree for DND1 was not generated since the orthogroup (OG0018857) only contains three genes (TraesCS5D01G404600, BdiBd21-3.1G0110600, and AT5G15410). Complete sets of T. aestivum homeologs from the three subgenomes were found in four out of the eight examined orthogroups. There were only two of three expected T. aestivum homeologs in the DMR6 orthogroup, with TraesCS4B02G346900 and TraesCS4D02G341800 representing the B and D subgenomes, respectively. A tblastn of these sequences to chromosome 4A revealed TraesCS4A02G319100, a partial match of 30-31% identity (1e-42 to 1e-44). This gene has low expression and is found in orthogroup OG0006808, which contains two other T. aestivum genes, one B. distachyon gene, two A. thaliana genes, and one H. vulgare gene (Supplementary Tables 4, 6). Despite the low sequence similarity, TraesCS4A02G319100 and TraesCS4B02G346900 are at more similar positions (4A:608043459 and 4B:640532917, respectively) to each other than to TraesCS4D02G341800 (4D:498572979). A tblastn to the entire genome revealed 20 other matches for the two DMR6 orthologs with 31-74% identity. Thus, it does not seem that the T. aestivum genome reference (Chinese Spring) contains a homeolog of DMR6 in the A genome. To further examine if a DMR6 ortholog is present in the A subgenome, the B and D subgenome DMR6 homeologss were BLASTed to the wheat pangenome CDS sequences 11 . Similar results were obtained; there were 30 total hits ranging from 36.765 to 39.564% identity (1.07e-61 to 1.99e-69) on chromosome 4A of the 10 genomes. Based on the low sequence identity, it is likely that there is not an ortholog of DMR6 on chromosome 4A. However, hits with high identity (97.619-98.214, e = 0) were found on chromosome 5A in all 10 genomes.
Gene ontology enrichment tests using GOslim annotations were conducted on the clusters to investigate functional processes. Across all eight S gene orthogroups, at least one gene from each is in a cluster with GO enrichment in at least one genotype (Figure 5). DMR6, FAH1, and WRKY25 are the only candidates to have enrichment in all three genotypes, AGD2 and DND1 only has enrichment in W2691, and BI-1, IBR3, and VAD1 have enrichment in both W2691 and W2691+Sr9b. Terms commonly enriched in the T. aestivum genotypes include the Golgi apparatus (GO:0005794), endosome (GO:0005768), endoplasmic reticulum (GO:0005783), protein binding (GO:0005515), transporter activity (GO:0005215), vacuole (GO:0005773), and peroxisome (GO:0005777; Figure 5). Only one GO term, catalytic activity (GO:0003824) is unique to Bd21-3, with other terms like DNA-binding transcription factor activity (GO:0003700) being enriched in the Bd21-3 and T. aestivum genotypes.

DISCUSSION
Susceptibility (S) genes are an essential component of compatible plant pathogen interactions (Engelhardt et al., 2018). The opportunity to genetically manipulate such genes to engineer disease resistance in important crops such as T. aestivum has captured significant scientific interest in recent years. However, our understanding of the genetic basis of disease susceptibility in cereals is limited to a few examples (van Schie and Takken, 2014;Engelhardt et al., 2018). Thus, important questions regarding the biological functions of these genes and their activation remain to be answered. As a first step to uncover putative stem rust S genes, we conducted a comparative RNA-seq experiment coupled with gene co-expression network analysis to determine transcriptional responses in T. aestivum genotypes and B. distachyon Bd21-3. We compared a compatible interaction (W2691) with an incompatible interaction controlled by the racespecific resistance gene Sr9b in the same genetic background (W2691+Sr9b). Sr9b restricts pathogen growth; however, it also allows the development of small sporulating colonies of a Pgt isolate which belongs to the race SCCL (Zambino et al., 2000). A more stringent incompatibility scenario is given by Bd21-3 genotype of B. distachyon, which allows restricted colony formation of Pgt without sporulation. These observations were consistent with previous descriptions of B. distachyon as a nonhost to rust pathogens (Figueroa et al., 2013, 2015, Omidvar et al., 2018. Thus, a strength of this study is the survey of molecular responses associated with increasing levels of susceptibility. Consistent with findings from other transcriptomic studies of wheat-rust interactions (Manickavelu et al., 2010;Zhang et al., 2014;Chandra et al., 2016;Dobon et al., 2016;Yadav et al., 2016), major transcriptional changes were detected in response to infection in both T. aestivum and B. distachyon, which reflect the complexity of these plant-microbe interactions. A significantly FIGURE 6 | Network diagrams for clusters containing orthologs of (A) DND1, (B) VAD1, and (C) DMR6 with corresponding plots showing log2 fold change of all nodes across 2, 4, and 6 dpi. Only connections with Z > = 3 are shown. Red lines, points, and counts represent T. aestivum and B. distachyon orthologs of S genes. Cluster identifiers (IDs) and gene names presented, left to right: DND1: 4 (TraesCS5D02G404600), 122 (TraesCS5D02G404600), and 652 (BdiBd21-3.1G0110600); VAD1: 0 (TraesCS2D02G236800), 0 (TraesCS2D02G236800), and 3085 (BdiBd21-3.1G0357000); and DMR6: 0 (TraesCS4B02G346900), 0 (TraesCS4D02G341800), and 4 (BdiBd21-3.1G1026800).
higher number of up-or down-regulated genes were found in T. aestivum than B. distachyon. The greater fungal colonization of T. aestivum as indicated by in planta fungal growth assays of Pgt is likely a result of the pathogen's failure to effectively manipulate the metabolism of B. distachyon. GOslim term analyses indicated an enrichment for Golgi apparatus, peroxisome, vacuole, and cell wall related functions in up-regulated genes in T. aestivum. These results are not surprising as a large proportion of immune receptors and plant defense signaling components play a role in plant-microbe interactions (Dodds and Rathjen, 2010;Couto and Zipfel, 2016). The plant Golgi apparatus and peroxisomes have been reported as targets of effectors from various pathogenic filamentous fungi (Robin et al., 2018). The enrichment of these GO terms in up-regulated genes in T. aestivum suggests that these CCs may be direct or indirect targets for effectors derived from Pgt. Analyses with the full GO term set revealed many enriched terms among downregulated genes related to photosynthesis in W2691 and W2691+Sr9b; a decrease in chlorophyll and photosynthetic activity has been previously reported in wheat infected with Pgt (Berghaus and Reisener, 1984;Moerschbacher et al., 1994).
Several S genes to diverse pathogens have been identified or postulated in various plant species (van Schie and Takken, 2014;Engelhardt et al., 2018). While this area of research for cereal rust pathogens is in its infancy, positive results from other pathosystems make a strong case to consider the modification of S genes as an approach to deliver durable and broadspectrum disease resistance. So far, only a few host-delivered avirulence effectors, AvrSr50 (Chen et al., 2017), AvrSr35 (Salcedo et al., 2017), and AvrSr27 (Upadhyaya et al., in press) from any cereal rust fungi have been isolated. These were identified in Pgt and how these effectors disrupt defense responses in compatible interactions remains unknown. Future research seeking to identify which plant proteins these effectors target will help elucidating S genes or processes required for stem rust susceptibility.
Here, expression patterns of gene orthologs in T. aestivum and B. distachyon corresponding to previously characterized S genes in H. vulgare and A. thaliana were examined to develop a framework to study S genes in wheat. A key focus of this study was to develop a workflow to extract orthologs with high expression in stem rust susceptible T. aestivum, but low expression in either T. aestivum with intermediate resistance, or B. distachyon. We note that differential gene expression between T. aestivum genotypes (W2691 and W2691+Sr9b) can also provide an opportunity to discovery S genes since rust infection in both genotypes differs. To link these candidate S genes with the biological pathways in T. aestivum and B. distachyon, we constructed GCNs, which can be explored to determine the role of components of these pathways and the complex interplay toward regulation of susceptibility in Pgt-T. aestivum interactions.
The biological functions of S genes in compatible-plant microbe interactions are diverse, as these genes play roles in a wide array of events that are critical for pathogen accommodation and survival (Engelhardt et al., 2018). Some of these susceptibility genes can act as negative regulators of immune responses, such as PTI, cell death, and phytohormone-related defense.
Our study determined that T. aestivum orthologs of the BAX inhibitor-1 (BI-1) gene in H. vulgare are candidate S genes, as these were up-regulated in W2691 (6 dpi) and W2691+Sr9b (4-6 dpi) whereas their expression in B. distachyon was not affected. BI-1 is an endoplasmic reticulum membrane-localized cell death suppressor in A. thaliana, and its wheat ortholog TaBI-1 (accession GR305011) is proposed to contribute to susceptibility in T. aestivum to the biotrophic pathogen Puccinia striiformis f. sp. tritici (Wang et al., 2012). Interestingly, the highest upregulation of the BI-1 was detected in the W2691+Sr9b genotype where it is necessary to regulate a HR upon Pgt recognition. Given this result it should be examined if BI-1 may be a conserved plant S factor to wheat rust fungi. Various orthologs of FAH1, which encodes a ferulate 5-hydroxylase in A. thaliana, were up-regulated in the T. aestivum genotypes upon Pgt infection (Mitchell and Martin, 1997). According to studies in A. thaliana FAH1 plays a role in BI-1-mediated cell death suppression through interaction with cytochrome b 5 and biosynthesis of very-long-chain fatty acids (Nagano et al., 2012). Additional findings further suggest that Pgt can also interfere with cell death signaling by altering VAD1 expression. The VAD1 gene encodes a putative membrane-associated protein with lipid binding properties and it is proposed to act as negative regulator of cell death (Lorrain et al., 2004;Khafif et al., 2017). Transcriptional activation of VAD1 has been shown to occur in advanced stages in plant pathogen interactions (Bouchez et al., 2007). We detected an upregulation of VAD1 orthologs in T. aestivum at 6 dpi, which is considered a late infection stage in the establishment of rust colonies.
Salicylic acid (SA) is a key phytohormone required to orchestrate responses to many pathogens (Ding and Ding, 2020). Similar to VAD1 whose function as a S factor is SA-dependent, we also uncovered other up-regulated genes that may also participate in defense suppression. The orthologs of the DMR6 are highly up-regulated in T. aestivum at 4 and 6 dpi in both compatible and incompatible interactions. As characterized in A. thaliana, DMR6 encodes a putative 2OG−Fe(II) oxygenase that is defense−associated and required for susceptibility to downy mildew through regulation of the SA pathway (van Damme et al., 2008;Zhang et al., 2017). The role of DMR6 in disease susceptibility holds significant promise to control diverse pathogens. For instance, mutations in DMR6 confer resistance to hemibiotrophic pathogens Pseudomonas syringae and Phytophthora capsici (Zeilmaker et al., 2015) and silencing of DMR6 in potato increases resistance to the potato blight causal agent, P. infestans (Sun et al., 2016). It has also been shown that the H. vulgare ortholog genetically complements DMR6 knock-out A. thaliana lines and restores susceptibility to Fusarium graminearum (Low et al., 2020). Gene orthologs of DND1 were also identified as up-regulated in both T. aestivum genotypes. The gene DND1 encodes a cyclic nucleotide-gated ion channel and its activity is also related to SA regulation . Mutations in A. thaliana DND1 display enhanced resistance to viruses, bacteria and fungal pathogens Jurkowski et al., 2004;Genger et al., 2008;Sun et al., 2017). We also noted that several wheat orthologs of the A. thaliana gene IBR3 also increased in expression as Pgt infection advanced. The role of IBR3 in susceptibility to P. syringae in A. thaliana has been confirmed by mutations and overexpression approaches (Huang et al., 2013). Consistent with our results, IBR3 is up-regulated in A. thaliana upon infection by P. syringae.
Plant transcriptional reprogramming triggered by pathogen perception is often mediated by WRKY transcription factors through activation of the MAP kinase pathways (Eulgem and Somssich, 2007;Rushton et al., 2010). Here, we detected an upregulation of the expression of WRKY25 orthologs that was most prominent at 6 dpi in the W2691+Sr9b genotype. The Arabidopsis gene AtWRKY25 is induced in response to the bacterial pathogen Pseudomonas syringae and the SA-dependent activity of AtWRKY25 is also linked to defense suppression (Zheng et al., 2007). According to results from this study, the contribution of orthologs of AGD2 to stem rust susceptibility in T. aestivum should also be examined. AGD2 encodes an aminotransferase and participates in lysine biosynthesis at the chloroplast (Song et al., 2004). Given that several oomycete and fungal effectors target the chloroplast (Kretschmer et al., 2020), effector research in cereal rust pathogens will be crucial to determine if these pathogens also target this organelle.
A classic example of S genes in barley is given by the Mlo gene (Jørgensen, 1992) in which a recessive mutation results in broad spectrum resistance to Blumeria graminis f. sp. hordei, the causal agent of powdery mildew. The Mlo gene family is highly conserved across monocot and dicot plants and gene editing of Mlo homeologs in wheat confers resistance to powdery mildew (Acevedo-Garcia et al., 2017). Interestingly, Mlo genes in T. aestivum have not been reported to provide protection against cereal rust diseases. Consistent with this, this study did not detect a significant change in the expression of Mlo alleles in T. aestivum genotypes (W2691 and W2691+Sr9b) over the course of the experiment.
One caveat of this study is that some S genes in T. aestivum for Pgt may not be found in model species like A. thaliana or detected using other pathogens. However, this is a first step to identify candidates to guide functional studies. While in this study we focused on orthologous S genes, the GCNs presented here are excellent resources to identify additional candidate S factors. It is possible that some of the genes included in clusters of these networks are part of the regulatory process that control expression of S genes or are part of essential pathways although their function may not be characterized yet in other systems. Future functional studies are required to validate the function of these genes in T. aestivum as S factors for rust infection and determine if these can be exploited for agricultural practice. A key aspect for the success of these novel approaches is the absence of plant developmental defects resulting from mutations of S genes. In some cases, the loss-of-function of negative regulators leads to constitutive activation of plant defense responses that manifest as poor growth or lesion-mimic phenotypes among other pleiotropic effects (Büschges et al., 1997;Ge et al., 2016). VIGS-mediated transient gene silencing (Lee et al., 2015), RNAi-mediated silencing Waterhouse and Helliwell, 2003;Sun et al., 2016), TILLING populations include some of the approaches to explore the potential use of these S gene candidates. Gene editing technologies through Zinc Finger nucleases, TALENs, CRISPR/Cas9 systems also offer options to generate transgene free plants (Urnov et al., 2010;Gaj et al., 2013;Wang et al., 2014;Jia et al., 2017;Nekrasov et al., 2017;Kim et al., 2018;Luo et al., 2019). In conclusion, as the demands for multi-pathogen durable disease resistance rise, our ability to target S genes may serve as a sound approach to harness genetic diversity and maximize the resources to meet critical these grand challenges.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, Bioproject PRJNA483957.

AUTHOR CONTRIBUTIONS
MF, CH, CM, and SK conceived and designed the study. VO, MM, FL, and MF conducted the experiments. EH, EG, J-MM, RD, CH, SG, JV, and MM contributed to data analysis. EH, BS, SK, CH, and MF interpreted results. EH, VO, CH, and MF wrote the manuscript. All authors contributed to manuscript editing, revisions and approved the submitted version.