Preliminary Study of Differential Expression
We conducted a preliminary RNA-seq study by inoculating reference Salix genotypes with M. americana to determine the optimum time post-inoculation to observe differential expression. We inoculated S. purpurea hosts ‘Fish Creek’ and 94006 with uredospores of M. americana isolate R15-033-03 and then extracted RNA at 0, 18, 42, 66, 90, and 114 HPI from inoculated leaves and un-inoculated control leaves. These two host genotypes were selected because 94006 is one of the grandparents and ‘Fish Creek’ is a parent of the F2 mapping population. A direct contrast between the inoculated and control treatment for each genotype-by-time was performed to generate a total number of differentially expressed genes (DEGs) up-regulated and down-regulated for each host genotype (Additional Fig. 1).
The total number of DEGs (p ≤ 0.05) for ‘Fish Creek’ were 0 (0 HPI), 0 (24 HPI), 5,589 (48 HPI), 562 (72 HPI), 1637 (96 HPI), and 3061 (120 HPI) (Additional Fig. 1), whereas DEGs for parent 94006 was 0, 0, 3796, 948, 597, and 1,293 for each ascending time-point, respectively (Additional Fig. 1). Neither parent displayed signs nor symptoms of infection during the experiment, however, signs of rust were visibly detected at 210 HPI. While uredospore sporulation appeared greater on ‘Fish Creek’ by 258 HPI, both genotypes were susceptible to the pathogen. The greatest number of DEGs was observed in both genotypes around 48 HPI (Additional Fig. 1). Thus, time-points 42 and 66 HPI were selected for the full experiment to capture the maximum host and pathogen response after inoculation.
Greenhouse inoculation of selected resistant and susceptible F2 genotypes
Based on field ratings of rust severity conducted in 2015 and 2017 in a replicated trial of an F2 S. purpurea QTL mapping population (25), 28 resistant and 28 susceptible genotypes were selected for controlled inoculation and 3' RNA-sEq. At 42 and 66 HPI, leaf discs were collected from six leaves of the 56 F2 genotypes, the two parents, and two grandparents of the population, as well as from uninoculated control plants. This experiment was conducted twice in separate greenhouses. Leaf rust severity was assessed in the inoculated treatment at 9 DPI as total percent leaf area coverage of uredospore pustules. The greenhouse ratings were moderately correlated with the 2015 and 2017 field ratings, with Pearson’s correlation values of 0.48 (p-value = 9.4 × 10− 5) and 0.53 (p-value = 1.6 × 10− 5), respectively. The susceptible genotypes had a significantly greater mean rust severity (44.8% - CV: 17%) than the resistant genotypes (28.1% - CV: 54%) based on a t-test (CI = 95%) despite considerably more variability among the resistant genotypes (Fig. 1).
Differential Expression Analysis of S. purpurea Transcripts
Two separate contrasts in DESeq2 were used to identify differentially expressed genes in this study. In the direct contrast between inoculated susceptible and resistant groups, there were 19 and 105 differentially expressed genes at time points 42 HPI and 66 HPI, respectively (Fig. 2A). Of the 19 DEGs at 42 HPI, six were up-regulated in the resistant genotypes, including a polyubiquitin protein (UBQ10), a plasma membrane intrinsic protein (PIP2;8), a phosphoglycerate kinase 1 (PGK1), a chaperone DnaJ-domain superfamily protein, and two genes of unknown function (DUF). The remaining 13 differentially expressed genes at 42 HPI were up-regulated in the susceptible genotypes and included several genes associated with the flavanone synthesis pathway. The 105 DEGs at 66 HPI consisted of 35 genes up-regulated in the resistant group, while the remaining 70 were up-regulated in the susceptible group. Genes up-regulated at 66 HPI in the resistant group include several involved in defense response such as: wall-associated kinase 2 (WAK2), WRKY DNA-binding protein 51, CAP superfamily protein, cytochrome P450, and chitinase A, but as a group, were not significantly enriched for any GO terms. Gene enrichment of the up-regulated susceptible genes were response to heat, stress, and reactive oxygen species (Additional Table 1).
The contrast of inoculated treatments versus uninoculated controls highlighted the response to infection. By performing separate paired analyses for both the resistant and susceptible groups, then intersecting DEGs, variable responses to inoculation were identified at each time-point. We classified DEGs as susceptible-specific, resistant-specific, and not-type-specific (common response between the resistant and susceptible groups). At both time points, the largest group of DEGs was the not-type-specific, positive log2-fold change (LFC) group, with 990 and 1862 genes at 42 HPI and 66 HPI respectively (Fig. 2B). All groups of DEGs that were up-regulated after inoculation were enriched for defense response at 42 HPI. However, only the resistant-specific and not-type-specific groups retained enrichment of upregulated defense response genes at 66 HPI. At 66 HPI, the susceptible-specific group lacked genes associated with defense response, but instead displayed upregulation of heat response genes (Additional Table 2). The resistant-specific and the susceptible-specific groups that were down-regulated at 42 HPI were both enriched for chloroplast components, with the susceptible-specific category also enriched for down-regulated ‘response to heat’ genes. There was no significant GO term enrichment at 66 HPI for genes down-regulated in the susceptible-specific category, while both the down-regulated resistant-specific and not-type-specific categories were enriched for genes associated with photosynthesis (Additional Table 2).
Network Analysis of S. purpurea Transcripts
A comparison between transcriptome-wide expression in the inoculated resistant and susceptible groups was performed in WGCNA, which defined co-expression modules based on correlated gene expression. Each module was randomly assigned a color name by the R package and is only relevant in distinguishing modules within networks, not in making comparisons between them. In this study, modules are referred to either as ‘R-module’ or ‘S-module’ to distinguish between those associated with resistant (R-) or susceptible (S-) plant networks. After removal of outlier samples and genes with low counts, the resistant network retained 75 samples and 16,410 genes, while the susceptible network retained 73 samples and 16,427 genes.
Of the 16,410 genes expressed in the resistant network, 10,176 genes were assigned to 14 modules, while the other 6,234 genes were assigned to the ‘grey’ module (unassigned genes). Modules sizes ranged from 33 to 5,085 genes, of which nine modules were correlated with time-point (Additional Fig. 2). The largest module ‘R - turquoise’ (n = 5,085) was positively correlated with time point (r = 0.92) and was the only module enriched for defense-related GO terms in the resistant network (Additional Table 3). The ‘R-blue’ module (n = 1,853) was inversely correlated with time point (r = − 0.89) and enriched for photosynthesis-related GO terms. A total of 10,977 genes in the susceptible network were placed into 15 modules, with the remaining 5,450 placed within the ‘grey’ module. Co-expression modules ranged in size from 25 to 4,661 genes, of which 12 were correlated with time point (Additional Fig. 2).
A hypergeometric test (p ≤ 0.05) facilitated a direct comparison between the resistant and susceptible networks to identify significant representation of the susceptible network modules within the ‘R-turquoise’ and ‘R-blue’ resistant modules. The ‘R-turquoise’ and ‘R-blue’ modules shared significant portions of four and six modules, respectively (Fig. 3A). Two modules correlated with time point in the susceptible network with significant ‘R-turquoise’ module representation were ‘S-turquoise’ (n = 4,661, r = 0.88) and ‘S-salmon’ (n = 89, r = 0.51), and were the only susceptible modules enriched for defense-related GO terms (Additional Table 4). Concomitantly, among the six susceptible modules represented within the ‘R-blue’ module and correlated with time point, only the ‘S-brown’ (n = 1,258, r = − 0.83) and ‘S-red’ (n = 264, r = − 0.57) modules were enriched for photosynthetic genes.
To gain insight in the role of hub genes in module composition, hub gene analysis was performed on the ‘R-turquoise’ and ‘R-blue’ modules, in addition to the ‘S-turquoise’, ‘S-salmon’, ‘S-brown’, and ‘S-red’ modules from the susceptible network (34, 35). Significant differences in mean expression of each module's hub genes and genes commonly co-expressed across networks were determined using Fisher’s least significant difference (p < 0.05). The ‘R-turquoise’ and ‘S-turquoise’ modules had 3,572 genes in common, yet at 42 HPI and 66 HPI the mean expression of these genes was greater among resistant genotypes (Fig. 3B). This trend persisted at 42 and 66 HPI among their respective hub genes, whose expression exceeded that of the shared genes. There were only 57 genes shared between the ‘R-turquoise’ and ‘S-salmon’ modules and were not differentially expressed throughout the experiment. However, the expression of ‘S-salmon’ hub genes did not significantly increase until 66 HPI, while ‘R-turquoise’ hub gene expression increased over time.
The ‘R-blue’ module from the resistant network was enriched for photosynthesis-related GO terms and shared commonly co-expressed genes with the ‘S-brown’ and ‘S-red’ modules from the susceptible network, similarly, enriched for photosynthesis (Fig. 3C). There were 773 shared genes in ‘R-blue’ and ‘S-brown’ modules with similar patterns of decreased expression over time, however, the mean expression of corresponding ‘R-blue’ hub genes was lower at each time point. The genes commonly co-expressed in ‘R-blue’ and ‘S-red’ only accounted for 188 genes that gradually decreased expression through time. Their hub genes, however, show that while the ‘R-blue’ genes decreased after 0 HPI and were beginning to level off by 42 HPI, the ‘S-red’ genes held similar expression throughout.
eQTL Analysis of S. purpurea Transcripts
Mapping of eQTL was performed using 22,068 SNPs and 16,270 genes to interrogate eQTL associated with the response to inoculation, removing those that were detected either at T0 or within the control treatment at the same time point. A total of 38,480 cis and 9,460 trans eQTL were identified at 42 HPI, 45,148 cis and 10,638 trans eQTL at 66 HPI, and 13,860 cis and 1,839 trans eQTL at both time points (Fig. 4A). Any SNP with more than 14 eQTL, the 95% confidence threshold identified through permutation, was identified as an eQTL hotspot. A hotspot was considered to be a locus influencing the regulation of multiple genes related to allelic phase. Simple correlation analysis (p < 0.05) condensed the significant eSNPs into eight eQTL hotspots at 42 HPI and six at 66 HPI (Fig. 4B). Hotspot sizes ranged from 14 to 55 eQTL associations and only three hotspots were enriched for any GO terms (Additional Table 5). The chr03 hotspot at 42 HPI (C3) was enriched for cell communication and signaling while the chr06 hotspot at 42 HPI (C6A) was enriched for chloroplast components. The only hotspot at 66 HPI showing GO enrichment was located on chr16 for photosynthesis and chloroplast components.
Candidate Genes for S. purpurea Resistance to M. americana
Candidate genes which potentially determine a compatible interaction (successful infection) between S. purpurea and M. americana were identified using the intersection of network analysis, differential expression, and eQTL mapping. Candidate genes were defined as the hub genes of modules found to be enriched for plant defense-related terms and differentially expressed either between resistant and susceptible genotypes or between the inoculated and control treatments. While associations with an eQTL hotspot for response to inoculation were not required for identification as candidate genes, it does aid in prioritization for further research. We identified candidate genes associated with the defense response enriched ‘R-turquoise’ module at 42 HPI (n = 31) and 66 HPI (n = 69) (Additional Table 6), of which 18 and 20 genes were correlated with leaf rust severity, respectively. Hub genes from the ‘R-blue’ module were associated with a reduction in photosynthesis through GO enrichment analysis. From these hub genes only 3 (42 HPI) and 21 (66 HPI) met our criteria for candidate gene selection, with all three genes at 42 HPI and one gene at 66 HPI having a significant correlation with leaf rust severity (Additional Table 6).
Differential Expression Analysis of M. americana Transcripts
Total raw reads of the inoculated treatments for each of the 60 willow genotypes (two replicates) were aligned to the M. americana reference genome R15-033-03 v1.0 (Crowell et al. 2021 - submitted). A direct contrast between genotypes previously identified as resistant and susceptible was performed at each time point (42 HPI and 66 HPI). A total of 22 M. americana genes were differentially expressed (FDR = 0.1) between the resistant and susceptible willow genotypes at 42 HPI, yet none at 66 HPI (Fig. 5). The majority of differentially expressed genes were up-regulated in the resistant group (20 genes) as compared to the susceptible (2 genes) (Additional Table 7). A BLAST search of these 22 DEGs was queried against the NCBI nt database (36). One transcript sequence (CDS_5062) was homologous to a known effector ubiquitin carboxyl extension protein in the plant parasitic nematode Globodera rostochiensis (37).
The in-silico proteome of the M. americana reference genome was analyzed using SignalPv5.0 using default settings to generate an in silico secretome, which resulted in 1,779 predicted secreted proteins (Additional Table 8) and analyzed for effector prediction using EffectorPv2.0 using the default settings (Additional Table 9). These proteins were then cross-referenced to the list of differentially expressed fungal transcripts between resistant and susceptible host groups. One (CDS_12834) of the 22 transcripts differentially expressed between the resistant and susceptible hosts was identified as a potential effector.