Transcriptome Analysis of Eucalyptus grandis Implicates Brassinosteroid Signaling in Defense Against Myrtle Rust (Austropuccinia psidii)

Eucalyptus grandis, in its native Australian range, varies in resistance to Austropuccinia psidii (syn. Puccinia psidii). The biotrophic rust fungus, A. psidii is the causal agent of myrtle rust and poses a serious threat to Australian biodiversity. The pathogen produces yellow pustules of urediniospores on young leaves and shoots, resulting in shoot tip dieback, stunted growth, and death. Dissecting the underlying mechanisms of resistance against this pathogen will contribute to improved breeding and control strategies to mitigate its devastating effects. The aim of this study was to determine the molecular dialogue between E. grandis and A. psidii, using an RNA-sequencing approach. Resistant and susceptible E. grandis seedlings grown from seed collected across its natural range were inoculated with the pandemic biotype of A. psidii. The leaf tissue was harvested at 12-h post inoculation (hpi), 1-day post inoculation (dpi), 2-dpi and 5-dpi and subjected to RNA-sequencing using Illumina 50 bp PE reads to a depth of 40 million reads per sample. Differential gene expression and gene ontology enrichment indicated that the resistant seedlings showed controlled, coordinated responses with a hypersensitive response, while the susceptible seedlings showed no systemic response against myrtle rust. Brassinosteroid signaling was apparent as an enriched term in the resistant interaction at 2-dpi, suggesting an important role of this phytohormone in defense against the pathogen. Brassinosteroid mediated signaling genes were also among the candidate genes within two major disease resistance loci (Puccinia psidii resistance), Ppr3 and Ppr5. While brassinosteroids have been tagged as positive regulators in other plant disease resistance interactions, this is the first report in the Eucalyptus – Austropuccinia psidii interaction. Furthermore, several putative resistance genes, underlying known resistance loci and implicated in the interaction have been identified and highlighted for future functional studies. This study provided further insights into the molecular interactions between E. grandis and A. psidii, contributing to our understanding of this pathosystem.


INTRODUCTION
Austropuccinia psidii (G. Winter) Beenken (formerly Puccinia psidii, commonly myrtle rust) is a biotrophic fungal pathogen causing widespread devastation to the natural and commercial range of numerous plant species within the Myrtaceae family. Originating in South America, this pathogen jumped from native Myrtaceae to introduced Eucalyptus and has since spread across the world causing devastation to natural and introduced Myrtaceae species (Carnegie and Pegg, 2018). Currently, A. psidii has been discovered in South and Central America (Coutinho et al., 1998) and Florida, California and Hawaii in the United States, North America (Marlatt and Kimbrough, 1979;Rayachhetry et al., 1997;Uchida et al., 2006). The pathogen has also been reported in countries such as Japan (Kawanishi et al., 2009), Australia (Carnegie et al., 2010), China (Zhuang and Wei, 2011), South Africa (Roux et al., 2013), New Caledonia (Giblin, 2013), Indonesia (Mctaggart et al., 2016), Singapore, and New Zealand (Du Plessis et al., 2019). Myrtle rust has an unusually broad host range on various members of the Myrtaceae family, affecting approximately 480 species within 86 genera, producing yellow pustules of urediniospores on growing leaves and shoots (Coutinho et al., 1998;Soewarto et al., 2019). This results in dieback of shoot tips, stunted growth and in severe cases plant death (Carnegie et al., 2016;Pegg et al., 2017Pegg et al., , 2020Fensham et al., 2020). Eucalyptus spp. are amongst the most susceptible members of this family to myrtle rust, posing a risk to their natural and introduced ranges (Coutinho et al., 1998;Junghans et al., 2003a), with Eucalyptus grandis showing variable responses ranging from highly susceptible to resistant.
Plants have evolved complex and integrated defense mechanisms against pathogenic organisms, including preformed and induced defenses. Preformed defenses include factors such as cuticle waxes on the leaf surface, cell walls and secondary metabolites (Naidoo et al., 2014). Dos Santos et al. (2019) found the chemical composition of the cuticle waxes in E. grandis and E. phaeotricha contribute significantly to susceptibility against A. psidii, with various cuticular signals influencing the interaction between host and pathogen. Additionally, older leaves of E. grandis have been shown to comprise of greater amounts of waxes than younger leaves which may affect the ability of A. psidii urediniospores to adhere to or penetrate the leaf surface. This results in pre-penetration resistance in mature trees against this pathogen (Xavier et al., 2015). When a urediniospore successfully attaches to the surface of the leaf, the spore germinates and develops a germ tube. From the germ tube, an appressorium forms which gives rise to a penetration peg that infiltrates the cuticle (Coutinho et al., 1998). Following leaf infiltration, haustoria develop in the host cells (Coutinho et al., 1998;Glen et al., 2007). It is through these haustoria that the pathogen acquires nutrients from the host plant to ensure its survival (Hahn and Mendgen, 2001) and sustained release of proteins that suppress the host defense responses (Jones and Dangl, 2006). These proteins are known as effectors (or avirulence proteins) that contribute to the success of the pathogen (Jones and Dangl, 2006). Junghans et al. (2003a) identified a locus associated with A. psidii resistance in E. grandis, termed Puccinia psidii resistance 1 (Ppr1), as per the previous classification of A. psidii. Ppr1 was shown to contribute to a large portion of resistance against A. psidii. This locus was found to be associated with numerous nucleotide-binding site leucine-rich repeat (NBS-LRR) resistance genes (R-genes) (Thumma et al., 2013), which suggests resistance against A. psidii is controlled by a combination of these genes. The discovery of this locus allowed targeted breeding programs, selecting for rust resistance to combat this pathogen. Despite the significance of Ppr1 in rust resistance, clones carrying this locus were shown to be susceptible to certain biotypes of A. psidii, indicating the locus was failing at providing resistance against A. psidii (Graça et al., 2011). Four additional loci contributing to myrtle rust resistance were discovered in E. globulus (Butler et al., 2016), termed Ppr2, Ppr3, Ppr4, and Ppr5. These quantitative trait loci (QTLs) determine if the plant shows symptoms (Ppr2 and Ppr3) and if there is a hypersensitive response (HR, Ppr4 and Ppr5) (Butler et al., 2016). The identification of additional QTLs contributing to disease resistance highlights the complexity of rust resistance. Despite the significance of these resistance loci, few studies have investigated the expression patterns of genes underlying these loci and little is known about how these genes contribute to disease resistance.
A recent study aimed to investigate the transcriptome of Melaleuca quinquenervia upon infection by A. psidii, to identify differences in responses among individuals with varying responses to the pathogen (Hsieh et al., 2018). At 5-dpi, they found that susceptible plants had higher levels of fungal transcripts than resistant plants, suggesting A. psidii could grow more successfully in susceptible leaves. Additionally, the authors found numerous differentially expressed genes (DEGs) associated with defense. These include NBS-LRR genes, WRKY transcription factors and Pathogenesis-related (PR) genes, such as PR-1, PR-2, PR-3, PR-4, and PR-5, with expression varying among resistant individuals. The authors suggest a potential role of the phenylpropanoid pathway in resistance against A. psidii, however, expression of associated genes differed among resistant samples and further studies are required to confirm this involvement. Sekiya et al. (2021), studying the interaction between E. grandis and A. psidii, found that rust resistance largely depends on the control of phenylpropanoid pathway related enzymes, corroborating the involvement of this pathway in defense against A. psidii. Tobias et al. (2018) studied the transcriptional responses of Syzygium leuhmannii challenged by A. psidii. The authors found up-regulation of receptor-like kinase genes, NBS-LRR genes as well as genes encoding enzymes involved in secondary metabolism which may contribute to disease resistance. Furthermore, they identified two TIR-NBS-LRR-type genes differentially expressed at 2-dpi with E. grandis homologs predicted within the Ppr1 locus. Santos et al. (2020) investigated the transcriptomes of resistant and susceptible E. grandis upon infection with A. psidii after 24-hpi. The authors found constitutive overexpression of defense-related genes in resistant compared to susceptible genotypes. These responses included up-regulation of salicylic acid (SA) associated genes, signal transduction and protein kinase leucine-rich receptors (PK-LRR). These studies provide some of the first insights into the transcriptional responses of Myrtaceae species challenged by A. psidii and further studies of these interactions are required to gain a deeper understanding of the molecular mechanisms governing host resistance.
The present study aimed to investigate the molecular interaction of resistant and susceptible Eucalyptus grandis challenged by A. psidii over time. In susceptible E. grandis genotypes, mycelia and haustoria have reportedly been observed between 12-and 18-hpi and in resistant genotypes, these have been observed between 18-and 1-dpi. Tobias et al. (2018) found significant transcriptome changes between resistant and susceptible samples 2-dpi whereas Hsieh et al. (2018) reported the possibility of fungal parasitic feeding stage at 5-dpi. Therefore, we hypothesize that at these and earlier time points (12-hpi, 1-, 2-, and 5-dpi) investigation of the expression profiles of resistant and susceptible seedlings would reveal molecular responses governing resistance, including changes in phytohormones, reactive oxygen species (ROS) and resistance genes. Through this, we have identified putative pathways and genes that contribute to host defense in E. grandis that can be manipulated in breeding and control programs, to mitigate the effects of myrtle rust. This is the first study to look at the transcriptome of E. grandis upon infection by A. psidii over a time course and provides further insights into the molecular interactions.

Austropuccinia psidii Inoculation Trial
Eucalyptus grandis seedlings were grown from seed sourced from wild plants across its natural distribution in eastern Australia, ranging from Coffs Harbor in New South Wales to northern tropical regions of Queensland. Seed was germinated in glasshouse conditions, where temperatures ranged from 20 to 30 • C. The seedlings, all with at least four young leaves, were inoculated to initially determine the phenotypes. Inoculation were generally carried out as described by Pegg et al. (2014). Spores of A. psidii were however collected from Syzygium jambos growing as street trees in suburbs of Brisbane, Queensland rather than using spores bulked from single spore inoculations in the laboratory. Spores were shaken into paper bags and passed through a sieve in the laboratory to remove plant debris. The collected spores were then desiccated in silica gel for 48 h and stored at −80 • C. Urediniospores of the pandemic biotype were removed from −80 • C storage and left at room temperature before addition of sterile distilled water with 0.05% tween 20. Spore were counted using a hemocytometer, with the concentration adjusted to 1 × 10 5 spore mL −1 . Six-weekold E. grandis seedlings were spray inoculated with A. psidii urediniospores using a fine mist at a pressure of 2.9 kPa. Spores were applied to the upper and lower leaf surfaces, and the seedling trays were covered with a plastic sheet for 24 h in the dark at 18 • C. After 24 h, the plastic covers were removed, and the seedlings were assessed for resistance and susceptibility based on their responses at 2 weeks post inoculation. Leaves were rated on a scale from 1 to 5, where rating one was considered resistant (R-interaction) and rating five was considered susceptible (Sinteraction) based on the system used in Junghans et al. (2003b) (Table 1).
Seedlings rated R and S were selected and separated out from other seedlings, and diseased tissues were removed by pruning. The seedlings were left to reshoot and after 8 weeks inoculation trials were repeated to obtain samples for RNA-seq. Seedlings were either inoculated as described above (infected), or mock-inoculated with sterile distilled water with 0.05% tween 20 (control). The youngest two sets of leaves of the infected and control plants were sampled as these are the most susceptible to A. psidii infection. Samples were collected from the seedlings at four time intervals namely, 12-hpi, 1-, 2-, and 5-dpi with every time point sampled from the same individual. Three replicates per time point were collected, comprising 14 seedlings per replicate to account for the genetic diversity within the seedlings.

RNA Extraction and Sequencing
Leaf samples collected from infected and control groups were subsequently frozen in liquid nitrogen. Total RNA was extracted as described by Zeng and Yang (2002) and treated with Qiagen RNase-free DNase 1 enzyme (Qiagen Inc.). Following successful RNA extraction, the samples were purified using Qiagen RNeasy Mini Kit as per the manufacturer's instructions. To determine the concentration and quality of the extracted RNA, the samples were analyzed using Bio-Rad Experion analyzer. Purified RNA from three biological replicates were submitted to the Beijing Genomics Institute (BGI) for mRNA-sequencing using 50 bp paired-end Illumina HiSeq 2,500 with an insert size of 300 bp and a sequencing depth of 40 million reads per sample.

Mapping and Analysis of Reads Against
Eucalyptus grandis V2.0 Genome Assembly RNA-seq data were analyzed using the Galaxy workspace (Goecks et al., 2010) and FASTQC v0.11.3 was used to assess read quality. Reads from the R-and S-interactions at each time point (both infected and control samples) were mapped to the E. grandis v2.0 reference genome (Myburg et al., 2014) using Bowtie 2 (Langmead and Salzberg, 2012), Tophat2 v2.0.14 (Trapnell et al., 2012) and featureCounts v1.5.0-p3 (Liao et al., 2013). Pathogen transcripts were removed for the purpose of this study and will be investigated in future studies using the A. psidii reference genome (Tobias et al., 2021). Data was normalized and variance stabilized transformation was applied in DESeq2 (Love et al., 2014) and transcripts per kilobase million (TPM) in R v3.6.0 (R Core Team, 2018). DEGs were analyzed using DESeq2 with filtering based on a twofold change (> 1.00 or < -1.00 log 2 (fold change)) and a false discovery rate of P < 0.001. We compared changes in gene expression between resistant, control relative to resistant, infected (R-interaction) and susceptible, control relative to susceptible, infected (S-interaction) to determine the effect pathogen infection would have on each interaction. Additionally, we compared susceptible, control relative to resistant, control (S-C vs. R-C) to determine the background gene expression. Heatmaps were generated using the R package ComplexHeatmap v2.0.0 (Gu et al., 2016) using log 2 (fold change) values.

Enrichment Analysis
Identified DEGs were analyzed for Gene Ontology (GO) enrichment over-representation in the category biological processes (BP), to identify biological pathways potentially involved in the interaction. Enrichment was calculated using Fisher's exact test and p-values were adjusted using the Benjamini and Hochberg (BH) method with a false discovery rate (FDR) of 0.05, with all calculations performed in R v3.6.0. Enrichment was determined separately for all time points as well as up-and downregulated DEGs. Additionally, Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation enrichment was performed on DEGs to identify enriched KEGG terms using the same method.

Eucalyptus grandis Disease Progression in Response to Austropuccinia psidii
Disease progression was monitored routinely for the presence of symptoms. Disease assessments were performed on the primary and secondary foliage (youngest two expanding pairs of leaves) to avoid the influence leaf age may have on disease severity levels.
Based on symptom development, seedlings were assessed 14dpi. Different responses were observed in E. grandis seedlings following inoculations as shown in Figure 1. Symptoms were similar to those observed in previous inoculations performed by Junghans et al. (2003b). This included light chlorotic flecking on seedlings rated one. Visible HR and necrotic lesion formation on seedlings rated two. Progression of necrosis and formation of uredinia pustules were observed in ratings three and four. Severe pustules were observed in seedlings rated five, with disease commonly occurring on the petiole and stem. For the purpose of this study, plants rated one were considered resistant (R) and those rated five on the severity scale considered susceptible (S), with samples harvested from these ratings for further differential gene expression analysis.

Differential Gene Expression in Resistant and Susceptible Eucalyptus grandis
To determine the effects pathogen infection would have on infected seedlings, we compared infected and control seedlings for both the R-and S-interactions across the time series. There were 4,681 significantly DEGs over the time course in the R-interaction, while there were 2,800 DEGs in the S-interaction over the time course (data not shown). During early infection at 12-hpi, there were 50 down-and 96 up-regulated genes in the R-interaction, while the S-interaction had 39 down-and 309 upregulated genes (Figure 2A). The number of up-regulated genes began increasing at 1-dpi in both the interactions. Up-regulation of genes peaked in the R-interaction at 2-dpi with 2,221 genes and 5-dpi in the S-interaction, with 2,146 genes (Figure 2A). Basal gene expression in resistant and susceptible plants might contribute to disease phenotypes when a pathogen is introduced, as previously observed in the A. psidii-E. grandis pathosystem (Santos et al., 2020). To investigate this, we compared susceptible, control relative to resistant, control (S-C vs. R-C) ( Figure 2B).
There was a total of 209 DEGs across the time course (data not shown). These results suggest that differences in basal expression exist.
To identify pathways putatively involved in resistance and susceptibility against A. psidii, GO enrichment was performed. No over-represented GO terms were identified in S-C vs. R-C, which is expected due to the low numbers of DEGs. Figure 3 shows a heatmap of over-represented defense-related terms in both the R-and S-interaction over the time series. Defense-related terms were separated into terms associated with phytohormone signaling and responses, oxidative burst, secondary metabolism and defense responses. Additionally, upand down-regulated terms are represented. Terms enriched in the R-interaction involve up-regulation of genes associated FIGURE 1 | Disease progression of Eucayptus grandis in response to Austropuccinia psidii. The top leaves show the adaxial while the bottom shows the abaxial view of the leaves. Rating one was considered resistant while ratings two to five were considered susceptible to A. psidii. Rating one showed slight chlorotic flecking, progressing in rating two and three where the flecking became necrotic and uredinia pustules began to form. Uredinia pustules became more apparent in rating four and severe pustules began forming by rating five. Rating one and rating five were selected for further analysis.
with brassinosteroid (BR) signaling, at 2-dpi, genes associated with a range of salicylic acid (SA) associated pathways as well as jasmonic acid (JA) and ethylene (ET) associated pathway. Enriched terms in the S-interaction involve upregulation of genes associated with JA and abscisic acid (ABA) pathways. Interestingly, the R-interaction showed upregulation of secondary metabolic pathways including terpenoid transport, while these responses with absent/down-regulated in the S-interaction. Moreover, defense responses such as callose deposition were up-regulated in the R-interaction, while these responses were limited in the S-interaction.
The DEGs identified in the R-and S-interaction were analyzed with KEGG to further elucidate the putative role these genes play in defense pathways. KEGG was used to determine pathways involved in defense against this pathogen. KEGG enrichment revealed up-regulation of "photosynthesis-antenna proteins" at 12-hpi, up-regulation of "terpenoid backbone biosynthesis" at 1-dpi and up-regulation of "amino sugar and nucleotide sugar metabolism" at 2-dpi unique to the R-interaction. Shared responses included plant-pathogen interaction at 1dpi, 2-dpi, and 5-dpi. Interestingly, at 5-dpi, the S-interaction showed down-regulation of cutin, suberin and wax biosynthesis (Figure 4).

Austropuccinia psidii Infection
In the R-interaction, BR signaling was observed at 2-dpi, with these responses absent in the S-interaction (Figure 3). Further investigations into the genes contributing to BR signaling revealed DEGs in both interactions, with the R-interaction expressing more of these genes at greater intensities. Signaling genes included BRI1-associated receptor kinase (BAK1), leucinerich repeat protein kinase family proteins and BZR1 family proteins. These genes showed different expression patterns across the R-and S-interaction over time, with the greatest expression observed in the R-interaction at 2-dpi. Similar responses were observed at 5-dpi in the S-interaction, although these terms were not enriched in the GO analysis ( Figure 5).

Expression of Candidate Genes Underlying Austropuccinia psidii Resistance Loci
To identify candidate genes underlying known resistance loci, the physical positions of the flanking markers were identified ( Table 2) and the genes between the markers were selected for further investigation (Supplementary Figure 1). Table 3 shows the total number of genes identified within each locus, with the number of DEGs underlying these loci. Additionally, the number of R-genes and differentially expressed (DE) R-genes underlying these loci are shown (Table 3 and Supplementary  Table 1). To determine the putative biological pathways that underlie the A. psidii resistance QTLs, the candidate genes identified were investigated for GO biological processes functional characterization (Supplementary Figure 1). Defenserelated enriched terms included "response to salicylic acid stimulus, " "defense response, " and "salicylic acid metabolic process." The involvement of these Ppr loci in defense-related processes highlighted throughout literature, the results of the GO BP enrichment analysis as well as the identification of numerous R-genes prompted further investigations into the expression profiles of the candidate genes underlying these QTLs.
To determine the roles the A. psidii resistance loci play in defense, DEGs underlying the QTLs were identified in the R-, S-interactions (Figure 6) and S-C vs. R-C (Figures 7A-I). Across all comparisons, there was significant differential expression of R-genes, namely NBS-LRR genes ( Table 3 and Supplementary  Table 2). Figure 6 is a heatmap representing log 2 (fold change) values of the DE R-genes underlying these QTLs in the R-and S-interactions. Differential expression of these R-genes is greater in the R-interaction, with up-regulation of numerous genes earlier in the interaction between 1 and 2-dpi when compared to the S-interaction. The S-interaction shows similar differential expression patterns later in the time series, between 2 and 5dpi. The DE R-genes underlying Ppr1, Ppr2, Ppr3, and Ppr5 had low TPM values, despite showing significant differential expression. Interestingly, Ppr4, a locus determining if the plant exhibits HR (Butler et al., 2016), showed increased expression of these R-genes, with expression of these genes greater than those observed at the other disease resistance QTLs. This corroborates the findings in Figure 3, where HR-associated responses were significantly enriched in both the R-and S-interactions.
To determine the differences in constitutive expression of the R-genes underlying the disease resistance loci, we analyzed S-C relative to R-C and found nine DE R-genes underlying these loci in the control samples (Table 3 and Supplementary Table 2). Figures 7A-I show the expression (TPM) of these genes over the time series. Eucgr.F01132 ( Figure 7G) and Eucgr.F01139 (Figure 7H), underlying Ppr4, are highly expressed, this is similar to the findings in the R-and S-interactions, in which R-genes were highly expressed. Eucgr.I00210 (Figure 7I), underlying Ppr5, showed greater constitutive expression in R-C at 12-hpi, 1-and 2-dpi, suggesting a role in constitutive resistance against A. psidii. Interestingly, Eucgr.C01921 (Figure 7B) underlying Ppr1 showed significantly greater DE in R-C at 2-dpi than in S-C. Gene Eucgr.C02650 (Figure 7F) showed greater expression in S-C at 12-hpi, with the expression of this gene increasing R-C at 5-dpi. The timing and intensity of basal expression of these genes may contribute to the phenotypes observed when resistant and susceptible E. grandis is challenged by A. psidii. Further studies into the functional role of these genes are required to determine if they contribute to rust resistance and how these genes differ at the molecular level between resistant and susceptible seedlings.
Genes involved in BR mediated signaling was significantly over-represented in the R-interaction at 2-dpi (Figures 3, 5).
To determine if the Ppr loci harbor genes associated with BR signaling, genes associated with this pathway identified through GO (GO:0009742) were compared to genes underlying the resistance loci. No BR mediated signaling genes were identified within Ppr1, Ppr2, and Ppr4. Two and four BR mediated signaling genes were identified within Ppr3 and Ppr5, respectively. These genes include a signaling kinase (BSK1), leucine-rich receptors (BRI1), protein kinases and a BR receptor (Supplementary Table 3).

DISCUSSION
The aim of this study was to compare the transcriptional responses of resistant and susceptible E. grandis with A. psidii over time. Through this, we identified pathways that may be involved in the E. grandis defense response against myrtle rust, such as the hypersensitive response, oxidative burst, and a range of phytohormones, with brassinosteroids signaling identified in resistant seedlings. Additionally, the timing of these responses varied between R-and S-seedlings, with delays and absences of responses in susceptible hosts. This suggests the importance in a coordinated, rapid response in resistance against A. psidii. Analysis of loci associated with resistance revealed these resistance loci harbor a plethora of R-genes, and we have identified candidate R-genes that are DE in the interactions, showing potential involvement in the disease outcome. Overall, we suggest a combination of these responses contribute to the resistance observed in E. grandis.

Selection of Resistant and Susceptible Seedlings
Despite the absence of physical symptoms of HR in the resistant seedlings, HR was significantly enriched in the resistant RNA-seq data, suggesting that this response is occurring at a molecular level. The plant material sourced for this study were seedlings rather than clones. This introduces genetic diversity within the each rating class, affecting the results obtained in downstream data analysis (Tobias et al., 2018). It is possible that the molecular responses observed are representative of the genetic diversity within the samples, with the HR varying among seedlings within the resistant rating class. This accounts for the variation seen in seedlings rated one, where few showed no visible infection and most showed light chlorotic flecking ( Table 1). A study conducted by Wang et al. (2018) showed that HR is involved in resistance against anthracnose (Colletotrichum) in the tea plant (Camellia sinensis). Through transcriptional analysis and histochemistry, the authors found the significant involvement of the HR, where cell walls were reinforced in resistant cultivars as a result of H 2 O 2 accumulation. Despite the involvement of HR in resistant cultivars, resistant leaves did not present with physical symptoms of disease or show signs of HR. Despite this, future studies will need to investigate the levels of HR associated compounds, such as H 2 O 2 , in resistant E. grandis to determine if this mechanism is associated with plant resistance at a molecular level, despite the lack of physical evidence of the HR.
Seedlings rated one were considered resistant while seedlings rated five were considered susceptible. This was characterized by either showing few symptoms of infection or large pustules on leaves, shoots and petioles ( Table 1). Seedlings rated two can also be considered resistant, as the presence of HR is indicative of host resistance against infection. However, the presence of necrosis suggests cell death as the result of disease rather than resistance, as described by Morel and Dangl (1997). Additionally, we observed in some instances that seedlings rated two developed pustules, despite the presence of an HR. For this reason, despite an HR observed in seedlings rated two, we only selected seedlings rated one and five for RNA-seq analysis to investigate the transcriptional responses of the two extremes of host resistance and susceptibility. This is a criteria also selected by Hsieh et al. (2018) when investigating the responses of M. quinquenervia upon A. psidii infection.

Improper Timing and Coordination May Lead to Host Susceptibility
Gene ontology biological processes functional characterization of DEGs in both the R-and S-interactions over time revealed interesting patterns of responses against A. psidii. In both interactions, there was significant up-regulation of defenserelated terms over time, with few or no down-regulation observed. These responses appeared to begin earlier in the S-interaction, with significant over-representation of defenserelated terms observed at 12-hpi. These included responses associated with SA biosynthesis and signaling, hydrogen peroxide metabolism and MAPK cascades (Figure 3). A possible explanation for the earlier increased expression of genes may be due to a lack of effective preformed defenses. A study by Hsieh et al. (2018) suggests that resistance in some individuals of M. quinquenervia against myrtle rust may be influenced by preformed defenses such as physical barriers and chemical constituents. The R-interaction may have effective preformed defenses that prevent appressorium formation at 12-hpi, while these barriers may be absent more susceptible hosts. Additionally, Xavier et al. (2001) reported that primary hyphae were observed at 12-hpi in susceptible E. grandis hosts, while only observed in the resistant hosts at 18-hpi. This is further supported by a study conducted by Dos Santos et al. (2019) in which susceptibility in E. grandis and E. phaeotricha was greatly influenced by the chemical composition of the cuticle waxes. Susceptible hosts lack effective preformed foliar compounds to prevent fungal penetration and as a result, respond to the presence of the pathogen earlier than in resistant hosts.  Despite this observation, the R-interaction shows a rapid increase in responses over 1-and 2-dpi, while similar responses are muted or absent in the S-interaction. Our study suggests that the S-interaction can induce a response against A. psidii, as evident by the changes in the transcriptome upon infection. However, it cannot implement adequate and sustained mechanisms to combat disease shown by the absence of critical defense responses by 2-dpi. A recent study showed that susceptible genotypes of S. luehmannii had no coordinated systemic response to A. psidii, with low host transcript numbers detected (Tobias et al., 2018). This is further supported by Santos et al. (2020), in which susceptible E. grandis genotypes were able to activate defense responses after 1-dpi with A. psidii, but the responses did not result in successful defense due to a lack of effective signal transduction. This suggests that susceptible E. grandis hosts and similarly of other Myrtaceae species, lack coordination in regulation of defense-related genes, which may contribute to disease susceptibility. A previous study on the interaction between A. psidii and M. quinquenervia showed susceptible FIGURE 7 | Differentially expressed (DE) R-genes (A-I) underlying the Austropuccinia psidii resistance loci in the susceptible, control compared to resistant, control seedlings over the time series, representing TPM values of DE R-genes in all biological replicates in the controls of the R-and S-interactions i.e., constitutive expression of the R-genes are being compared, where black represents resistant and light gray represents susceptible seedlings and asterisks represent significant DE. plants had various up-regulated defense-related responses at 5-dpi (Hsieh et al., 2018). This is similar to the results obtained in this study, suggesting that Myrtaceae species are able to deploy defenses against A. psidii, but delays in these responses may contribute to susceptibility. Over-compensation later in the interaction at 5-dpi may be too late to prevent spread of disease.
Interesting KEGG over-represented terms included "photosynthesis-antenna proteins" up-regulated in the R-interaction at 12-hpi, with this term absent in the S-interaction over the time series. Photosynthetic processes are known to play roles in plant defense against pathogenic invasion. Santos et al. (2020) found the constitutive involvement of this pathway in the defense response in resistant E. grandis genotypes. Furthermore, "amino sugar and nucleotide sugar metabolism" was unique to the R-interaction at 2-dpi. Sugar signaling has been implicated in disease resistance, with high levels of sugars in plants contributing to defense responses. Increased plant sugars contribute to resistance by inducing expression of PR genes, interacting with the plant phytohormone biosynthesis and signaling pathways (Morkunas and Ratajczak, 2014) and increasing production of ROS (Morkunas et al., 2008). Exogenous applications of sucrose in Vitis vinifera induces the production of secondary metabolites such as phenylpropanoid, which may contribute to disease resistance modulated by sugar signaling (Ferri et al., 2011). The upregulation of phenylpropanoid biosynthesis identified in this study (Figure 4), as well as identification of this pathway in other studies investigating interactions with myrtle rust, suggests the importance of phenylpropanoids in resistance, with sugar metabolism putatively contributing to the induction of this pathway. Moreover, sugar transport, mediated by the sugar transport protein TaHTP, was found to contribute to resistance against the biotrophic rust pathogen Puccinia triticina in wheat (Savadi et al., 2018).

Eucalyptus grandis Phytohormone Responses to Austropuccinia psidii
Phytohormones contribute significantly to plant defense responses, with crosstalk existing between them to facilitate defense against pathogenic invasion. In this study, genes associated with phytohormones were extensively enriched in both resistant and susceptible E. grandis seedlings upon infection with A. psidii (Figure 3). These include genes relating to SA, JA, ET, ABA, and auxin. Our results suggests putative involvement of SA biosynthesis and signaling, with downstream responses such as systemic acquired resistance identified (Figure 3) in both the resistant and susceptible hosts. Santos et al. (2020) found the up-regulation of SA signaling and responsive genes when comparing mock-inoculated myrtle rust resistant and susceptible E. grandis genotypes. These findings suggest the importance of SA in E. grandis defense against A. psidii, both at a constitutive and induced level. Interestingly, GO analysis revealed the involvement of BR at 2-dpi in the R-interaction, a pathway absent in susceptible hosts.
BRs are polyhydroxylated steroidal phytohormones found in plants that have versatile roles in many plant physiological processes. These processes include those such as cell growth and development, seed germination and response to abiotic stress (Saini et al., 2015;Yu et al., 2018;Kim and Russinova, 2020). In recent years, studies have emerged implicating BRs in resistance to various pathogens, with complex molecular interactions driving host responses against pathogenic invasion (Choudhary et al., 2012). BRs have been shown to confer resistance against biotrophic bacteria (Huang et al., 2014) and exogenous application confers systemic resistance to biotrophic fungi such as Oidium sp. and Magnaporthe grisea (Nakashita et al., 2003).
Many studies on the roles BR play in modulating immune responses have focused on the leucine-rich repeat receptorlike kinase (LRR-RLK) and BRI1 (BR-insensitive 1)-associated kinase (BAK1). BAK1 is a co-receptor of BRI1 and plays a role in BR signal transduction as well as interacting with pattern recognition receptors (PRR) to induce a PTI response against invading pathogens (Li et al., 2002;Vert, 2008;De Bruyne et al., 2014). Brassinolide (BL), an active form of BR, binds to BRI1 leading to autophosphorylation and the activation of BAK1 (De Bruyne et al., 2014). This ultimately leads to downstream transcriptional changes associated with BR-signaling. BAK1 has been implicated in disease resistance against fungal pathogens including Sclerotinia sclerotiorum and Verticillium dahliae (Fradin et al., 2009;Chaparro-Garcia et al., 2011) and the identification of nine DE BAK1 genes in the present study suggests a putative role of BR and associated receptors in the defense response against A. psidii. The higher number of DE BAK1 genes in response to A. psidii in the R-interaction may contribute to disease resistance. Similarly, Santos et al. (2020) found constitutive over expression of BAK1 when comparing mock-inoculated A. psidii resistant and susceptible E. grandis at 1-dpi, corroborating the findings of the present study. Additionally, in a transcriptome study between banana and Fusarium oxysporum f. sp. cubense, up-regulation of BAK1 was identified in relatively resistant banana plants, potentially implicating this gene in the innate immune response in banana (Li et al., 2013).
In addition to BAK1, the somatic embryogenesis receptorlike protein kinase-1 (SERK1) are implicated in disease resistance in tomato against the foliar fungal pathogen, Cladosporium fulvum (Fradin et al., 2011). The SERK protein family, which includes SERK3/BAK1 has been extensively shown to contribute to disease resistance. In this study, five SERK1 genes were DE in E. grandis upon challenge by A. psidii (Figure 4). The expression of these genes is similar to that of BAK1. Further analyses of the genes involved in BR signaling revealed DE of gene BZR1 encoding brassinazole-resistant 1. This is a transcription factor responsible for repressing expression of defense-related genes (Robert-Seilaniantz et al., 2011). Further studies into the role this gene is playing during the interaction of E. grandis are required. However, up-regulation of this gene was low or absent over the time course within both interactions. This suggests defense-related responses involved in PTI may not be repressed.
Despite our results suggesting a role of BR signaling in defense against A. psidii, studies on the effector AVR2 have shown that increased BR signaling in transgenic potato plants leads to enhanced susceptibility to Phytophthora infestans. This effector manipulates the crosstalk that exists between plant growth and immunity to facilitate enhanced pathogen colonization (Turnbull et al., 2017). This highlights the complexities of BR signaling and the roles this phytohormone pathway plays in plant defense. Different pathogens employ different virulence mechanisms to elicit disease in host plants, and this may explain the different outcomes observed in BR signaling between P. infestans (oomycete) and A. psidii (fungus). Moreover, the lifestyle of the pathogen contributes to the outcome of BR in plant defense. Despite AVR2 being up-regulated during the biotrophic phase in P. infestans, this oomycete is a hemi-biotroph while A. psidii is a biotroph and studies have shown that BR increases plant defense in response to biotrophs while decreasing plant defense against hemibiotrophs and necroptrophs (De Vleesschauwer et al., 2012;Turnbull et al., 2017;Yu et al., 2018), this contributes further to the complexities of this phytohormone in plant responses. For this reason, more studies are required to understand how BRs contribute to the responses of E. grandis to A. psidii.
The identification of BR mediated signaling genes underlying A. psidii resistance loci may implicate this pathway in the resistance conferred by these loci. Ppr3 harbors two genes that are associated with BR signaling and BR receptors, while Ppr5 harbors four genes that code leucine-rich receptor proteins and protein kinases. Eucgr.G01832 (BRI1-interacting signaling kinase, BSK1) underlies Ppr3, while Eucgr.I00302 and Eucgr.I00303 (BRI1) underlies Ppr5 (Supplementary Table 3). BRs associates with BRI1, causing phosphorylation and subsequent dissociation of BRI1 from its associated inhibitors. This cascade event leads to the activation of BAK1 and the phosphorylation of BSK1 (De Bruyne et al., 2014;Yu et al., 2018). This causes activation of downstream BR signaling responses and changes in gene expression. Despite these loci harboring BR signaling genes, these were not DE in this study. However, extensive polymorphisms may exist within these genes that affect the expression as well as the final protein product, which may contribute to BR signaling and resistance against A. psidii. Further studies are required to identify single nucleotide polymorphisms (SNPs) that may exist between resistant and susceptible seedlings to determine how these genes may contribute to disease resistance.

Constitutive and Induced Expression of Candidate Resistant Genes
Plants recognize pathogen effectors through resistance proteins, most commonly the nucleotide-binding domains and leucinerich repeats (NLRs) (Jones and Dangl, 2006). Effectors are small secreted proteins (SSPs) that the pathogen deploys to facilitate the manipulation of host responses, allowing the pathogen to acquire host nutrients and reduce plant defense (Fudal et al., 2018). The recognition of these pathogenic effector proteins leads to effector-triggered immunity (ETI), resulting in oxidative burst of ROS as well as hypersensitive responses (HR) (Jones and Dangl, 2006). These responses lead to extensive modifications in defense-related gene expression and prevention of pathogen spread from the site of infection (Dodds and Rathjen, 2010).
In Eucalyptus, numerous R-gene loci conferring resistance to A. psidii have been discovered (Junghans et al., 2003a;Butler et al., 2016). Ppr1 confers resistance to myrtle rust in E. grandis while Ppr2-5 has been shown to confer resistance to myrtle rust in E. globulus, a close relative of E. grandis (Butler et al., 2016). Little is known about the changes in expression over time of genes associated with these loci and how they contribute to disease resistance. To further investigate the genes underlying these loci, we identified various DEGs in the interaction between A. psidii and E. grandis, many of which were R-genes (Figures 6, 7A-I and Supplementary Table 2). In general, there were more significantly DE Rgenes in the R-interaction underlying these loci than in the S-interaction. Additionally, the R-interaction showed significant DE predominately between 1-and 2-dpi, while significant expression was only observed in the S-interaction during later stage infection. A recent study by Santos et al. (2020) identified two significantly DE R-genes underlying Ppr1 when investigating mock-inoculated resistant and susceptible E. grandis as well as susceptible infected and control E. grandis, namely Eucgr.C02650 and Eucgr.C02749. Eucgr.C02650 was down-regulated in both, while Eucgr.C02749 was up-regulated in both. In our study, we found significant DE of Eucgr.C02650 when investigating resistant, infected relative to control as well as investigations between resistant and susceptible controls.
Ppr2 and Ppr3 determine if the plant shows disease symptoms (Butler et al., 2016). In the R-interaction, these loci showed significant up-regulation at 1-and 2-dpi. In contrast, while the expression patterns are similar, the S-interaction showed fewer significant DE of these R-genes. This may contribute to the disease symptoms observed in susceptible seedlings. Ppr4 and Ppr5 determine if the plant exhibits a hypersensitive response (Butler et al., 2016). Both hosts showed down-regulation of these genes in Ppr4, suggesting the pathogen is potentially manipulating the HR that are controlled by these R-genes (Figure 6). Despite this, there was significant expression of some of these R-genes in both interactions at Ppr4 showing high TPM values. Similarly, at Ppr5, the R-interaction shows more DE than in the S-interaction, with significant up-regulation in the R-interaction at 2-dpi. This suggests an involvement of these Rgenes in the HR observed in the resistant hosts. This observation correlates with those made in the GO responses in which there is a delay in HR in the S-interaction.
It is important for plant defense responses to be suppressed in the absence of plant pathogens, as constitutive expression of defense-related compounds may be detrimental to plant growth and development (Marone et al., 2013). However, constitutive expression of defense-related genes, including R-genes, has been shown to contribute to disease resistance against pathogens. For example, in Arabidopsis ADR1 codes a specific CC-NB-LRR resistance gene. Mutations in this gene led to the constitutive expression of this R-gene, conferring broad-spectrum disease resistance against biotrophic pathogens and drought tolerance (Grant et al., 2003;Chini et al., 2004). Furthermore, constitutive expression of a TIR-NB-ARC-LRR disease resistance protein conferred resistance to phytopathogenic fungi and bacteria (Wen et al., 2017). Our results suggest that the constitutive expression of R-genes underlying resistance QTLs may putatively contribute to disease resistance. These results corroborate those from a recent study investigating the transcriptome of E. grandis, which found constitutive expression of genes encoding R proteins, such as Eucgr.C02650 underlying Ppr1, that was found to be significantly constitutively DE in the present study. The authors suggest this constitutive expression of defense-related genes potentially contributes to disease resistance observed against A. psidii (Santos et al., 2020).
The identification of DE R-genes when investigating the E. grandis transcriptomes, both at the induced and constitutive levels highlights candidate genes putatively involved in the defense responses and further investigations into the molecular differences that exist within these genes between genotypes is required to elucidate their roles.

CONCLUSION
This study investigated the molecular interactions of E. grandis in response to A. psidii to investigate mechanisms underlying host resistance. We identified pathways and genes potentially involved in the interaction. This includes brassinosteroid signaling, a pathway not previously reported in the A. psidii-host interactions and various R-genes that have different expression patterns, both at the constitutive and induced levels. Moreover, our findings suggest that resistant and susceptible hosts share similar responses, but timing and intensity determines disease severity. Due to the bioinformatic nature of RNA-seq analysis, further functional studies are required to characterize the observed responses to validate the results in the present study. This includes studies that investigate the brassinosteroid levels between resistant and susceptible seedlings to confirm the involvement of this pathway in host resistance. Moreover, investigations into the identified resistance genes are required to elucidate the roles they play during the interaction. Future studies should also aim to include investigations of the pathogen responses to identify mechanisms A. psidii is employing to manipulate susceptible hosts, as it will provide insight into novel targets of pathogen control, facilitated by the use of the A. psidii reference genome (Tobias et al., 2021). In conclusion, the results from this study will contribute to improved selection and breeding of E. grandis by manipulating pathways and genes tagged as important contributors to resistance.

DATA AVAILABILITY STATEMENT
The data have been uploaded to the NCBI under bioproject PRJNA763498.

AUTHOR CONTRIBUTIONS
LSS and SN conceived the study. SN obtained funding for the study. LSS and GSP performed the experiments and collected sample material. SS and CNO analyzed the data. SS wrote the first draft of the manuscript with input from SN. SS, CNO, LSS, GSP, and SN commented and revised the manuscript. All authors approved the submitted version.

FUNDING
This work was supported by the South African National Research Foundation (NRF) Y-rated Grant to SN (UID105767) and post doctoral support to LSS. The authors acknowledge funding from the Technology Innovation Agency of South Africa through the Forest Molecular Genetics Cluster Program to SN, CNO, and SS. Opinions and conclusions arrived at by the author(s) are not necessarily attributed to the NRF and other funding agencies.