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⁵ spore mL^–1. Six-week-old 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 (S-interaction) 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.

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.

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 two-fold change (> 1.00 or < –1.00 log₂ (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₂ (fold change) values.

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 down-regulated DEGs. Additionally, Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation enrichment was performed on DEGs to identify enriched KEGG terms using the same method.

To identify candidate genes underlying A. psidii resistance QTLs (Ppr1, Ppr2, Ppr3, Ppr4, and Ppr5), the physical positions of their flanking markers were obtained (Table 2; Junghans et al., 2003a; Mamani et al., 2010; Butler et al., 2016). Phytozome v12.1¹ was used to identify genes underlying these physical positions of each locus on its respective linkage group using the genome browser.

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 14-dpi. 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.

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 up-regulated 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, up- and down-regulated terms are represented. Terms enriched in the R-interaction involve up-regulation of genes associated 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 up-regulation of genes associated with JA and abscisic acid (ABA) pathways. Interestingly, the R-interaction showed up-regulation 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 1-dpi, 2-dpi, and 5-dpi. Interestingly, at 5-dpi, the S-interaction showed down-regulation of cutin, suberin and wax biosynthesis (Figure 4).

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), leucine-rich 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).

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). Defense-related 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₂ (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 5-dpi. 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).

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₂O₂ 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₂O₂, 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.

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 defense-related 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 defense-related 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 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 up-regulation 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).

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 receptor-like 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 receptor-like 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.

Plants recognize pathogen effectors through resistance proteins, most commonly the nucleotide-binding domains and leucine-rich 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 R-genes 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 R-genes 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.