Introduction

Sugarcane is currently cultivated in 27 million hectares worldwide, 52% of which are located in Americas, mostly in South America (11 million ha). Brazil is the leading producer, responsible for 758.55 million tonnes (FAOSTAT, 2017). The term sugarcane includes at least six species of the genus Saccharum (Daniels and Roach 1987; Kim et al. 2014). The ancestor species are Saccharum robustum and S. spontaneum, which present variable ploidy levels with chromosome number ranges of 2n = 60–170 and 2n = 36–128, respectively (Irvine 1999; Zhang et al. 2012). The octoploid (2n = 80) S. officinarum is an ancient cultivated sugarcane with high sugar content that formed the group of noble cultivars  (Grivet et al. 2004; Irvine 1999; Paterson et al. 2013). Sugarcane plants currently used for industrial purposes are modern cultivars derived from interspecific crosses, mainly between S. officinarum and S. spontaneum. As a consequence of successive backcrosses, the majority of chromosomes in homology groups are from S. officinarum, with a lower frequency of S. spontaneum chromosomes, along with some recombinant chromosomes and frequent aneuploidy (Grivet and Arruda 2002; Piperidis et al. 2010).

Sugarcane is important for sugar and ethanol production and, moreover, new perspectives are being emerged because of its capacity for energy production and due to new uses for the industrial residues (Carpio and Simone de Souza 2017; Manochio et al. 2017; Moraes et al. 2015; Rodrigues Reis and Hu 2017). Despite its economic importance, sugarcane production is affected by insect and nematode attacks, weed competition and diseases, leading to losses in productivity (Cheavegatti-Gianotto et al. 2011). Among the diseases currently affecting sugarcane crops, three important biotrophic fungi that stand out are: (i) Puccinia melanocephala, responsible for brown rust; (ii) Puccinia kuehnii, causal agent of orange rust; and (iii) Sporisorium scitamineum, which causes smut. Both rusts affect leaf tissue, while smut affects meristem development. Particularly, more attention was given to orange rust following an outbreak in Australia roughly eighteen years ago (Magarey et al. 2001). Occurrences were then reported in America in the last 10 years (Comstock et al. 2008; Funes et al. 2016), but noticeably with smaller diversity of isolates than that observed in the Eastern hemisphere, suggesting a spread of particular strains from Asia and Australia (Glynn et al. 2010). The impacts of P. kuehnii on sugarcane yield have been associated with variation in photosynthesis, because there are considerable changes in leaf physiological measures in susceptible cultivars where symptoms are evident (Zhao et al. 2011; Chapola et al. 2016).

Biotrophic pathogens require living host tissues to capture nutrients for their development. The establishment of biotrophic fungi consists in spore adhesion, germination on plant surface and formation of appressoria, which are the penetration structures. Later, the spread of fungal cells occurs by different modes of colonization. Obligate biotrophs, like the species of the genus Puccinia, initially grow hyphae inside the host tissues, moving to intracellular colonization by haustoria, side branches from these hyphae (O’Connell and Panstruga 2006; Spanu and Panstruga 2017). Haustoria act as an interface to acquire nutrients from the host and to deliver effectors (O’Connell and Panstruga 2006; Voegele and Mendgen 2003). Effectors are essential molecules secreted by pathogens to control the plant immune system. They may be recognized by plant resistance proteins, which is generally followed by an oxidative burst and promotion of the hypersensitive response. The relevance of the latter is that it induces cell death to restrict pathogen development, ensuring the non-capture of nutrients and water. Resistance gene products also activate hormone-dependent signaling, by salicylic (SA) or jasmonic acid (JA), which regulates the expression of other defense response genes (Glazebrook 2005; Spanu and Panstruga 2017).

Plant infection is a dynamic process and, at the RNA level, the expressed genes of host and pathogen reflect the so-called interaction transcriptome (Birch and Kamoun 2000). Through this concept, it is possible to capture the profiles of transcripts with relevant roles, such as the members of the signaling pathways which prevent disease establishment. Techniques that measure changes in RNA levels such as RNA Sequencing (RNA-Seq) are used to detect differentially expressed genes. This methodology is based on the development of Next-Generation Sequencing (NGS) technologies, enabling the study of variations in RNA composition between samples from distinct experimental conditions. RNA-Seq was used in sugarcane to provide molecular information about transcripts from a range of superior genotypes (Cardoso-Silva et al. 2014), understanding changes in lignin content (Vicentini et al. 2015) and the development of organs (Mattiello et al. 2015). With regard to pathogenic interactions, RNA-Seq has the sensitivity to capture both pathogen and host transcripts, not only providing catalogs of transcripts, but also changes in gene expression levels due to the infection process (Westermann et al. 2012). In plants, transcriptome profiling studies are covering pathosystems in annual (Cho et al. 2015; Dobon et al. 2016; Loarce et al. 2016; Tremblay et al. 2013; Xia et al. 2017) and perennial species (Gervasi et al. 2018; Kan et al. 2017; Páez et al. 2015; Teixeira et al. 2014; Zhu et al. 2017). Recently, studies concerning the interaction of sugarcane with the fungus S. scitamineum (McNeil et al. 2018; Que et al. 2014; Schaker et al. 2016), with the bacteria Acidovorax avenae subsp. avenae  (Santa Brigida et al. 2016) and with the bacteria Leifsonia xyli subsp. xyli were published  (Cia et al. 2018).

Knowledge about the molecular mechanics of infection in sugarcane is important and there are efforts in identifying molecular markers for disease resistance  (Glynn et al. 2013; dos Santos et al. 2017; Yang et al. 2018). Efforts were taken for developing molecular markers associated with resistance against both sugarcane rusts. For brown rust, markers associated with the Bru1 QTL have been used in sugarcane research  (Glynn et al. 2013; Parco et al. 2017; Yang et al. 2017), mainly for breeding purposes  (Glynn et al. 2013), as well as for investigations about its occurrence in genotype panels  (Parco et al. 2017). With regard to orange rust, cultivars have shown distinct responses when infected with P. kuehnii strains  (Moreira et al. 2018). Primarily, studies in orange rust focused on molecular detection of P. kuehnii using primers that amplify ITS regions  (Glynn et al. 2010), the characterization of pathogen developmental structures during infection  (Gomez 2013) or the assessment of infection degree in sugarcane cultivars  (Chapola et al. 2016; Zhao et al. 2011). Yang et al. (2018) detected three QTL associated with resistance against orange rust. They also identified putative resistance genes based on sugarcane transcripts with orthologs in the sorghum genome. For doing so, they searched for homology with leucine-rich-repeat (LRR) and nucleotide binding (NB) sites using five combined databases. The authors identified a marker associated with orange rust resistance as a possible candidate to be used for marker assisted selection. In addition, non-inoculated susceptible and resistant genotypes were grouped based on metabolite fingerprinting  (Leme et al. 2014). However, there are no previous studies gauging gene expression changes as a result of rust infections in sugarcane. To our knowledge, there is no information about high-throughput molecular analyses during orange rust establishment in sugarcane. Our work thus aims to study a time series of transcriptional changes in susceptible sugarcane infected by P. kuehnii.

Materials and methods

Plant material

Plants of the sugarcane cultivar SP89-1115, susceptible to orange rust  (Barbasso et al. 2010; Chapola et al. 2016; Moreira et al. 2018), were grown in 500 mL plastic cups containing sterilized substrate in a greenhouse during 40 days and fertilized with an ammonium sulfate solution (30 g/L) ten days before inoculation. One day before inoculation, spores of P. kuehnii were obtained using a vacuum pump (Prismatec 101) from the abaxial side of infected leaves from Centro de Tecnologia Canavieira (Piracicaba, São Paulo, Brazil) that contained open pustules. A solution was prepared by homogenizing the spores for 30 min in distilled water, with a final concentration of \(10^3\) spores/mL. Approximately 4 mL of this solution was sprayed over the plants, preferably reaching the abaxial side of leaves, and they were kept at 22 °C for 24 h in a dew chamber. After this period, they were moved to a growth chamber at 25 °C with a photoperiod of 12 h light/12 h dark. A more detailed description of these experimental procedures is available at Gomez (2013). For molecular analyses, five biological replicates of +1 leaves were taken at 12 h after inoculation (hai), 24 hai, 48 hai, 5 days after inoculation (dai) and 12 dai. Biological replicates of uninoculated control leaves were also collected, with an experimental design similar to other works (Dobon et al. 2016; Li et al. 2015; Lysøe et al. 2011; Yadav et al. 2016; Zhang et al. 2014).

RNA extraction and sequencing

Total RNA from three biological replicates of each time point was extracted with the RNeasy Plant Mini Kit (Qiagen) in a final volume of 60 \(\mu\)L. We assessed the quality and quantity of the RNA via electrophoresis on 2% agarose gel, ND-1000 spectrophotometer (NanoDrop) and 2100 Bioanalyzer (Agilent Technologies). Next we used aliquots to prepare mRNA-Seq indexed-libraries with the TruSeq Stranded kit (Illumina). A total of 18 libraries were sequenced in three lanes of an Illumina HiSeq 2500 to obtain 125 bp paired-end reads. Some libraries were sequenced in multiple lanes, such that we could assess the occurrence of batch effects. The sequencing data has been deposited at DDBJ/EMBL/GenBank under the BioProject ID PRJEB31605.

Data processing, de novo assembly and functional annotation

Following a strategy similar to Hoang and colleagues  (Hoang et al. 2017), adapter sequences and low quality bases were removed with Trimmomatic v.0.38 (Bolger et al. 2014), using windows with a minimum average Phred quality score of 20. The first 12 bases of each read were removed and only reads with at least 75 bases were kept. We used sortMeRNA v.2.1 (Kopylova et al. 2012) to filter contaminant rRNA, by aligning the reads to prokaryotic and eukaryotic rRNA databases, then selecting the remaining reads. The processed reads of the 18 libraries were used for transcriptome assembly (Supplementary Table 1).

Transcriptome de novo assembly was performed using Trinity v.2.5.1  (Grabherr et al. 2011) using the default k-mer size of 25. We set the maximum number of reads to combine into a single path (max_reads_per_graph) to 3,000,000. Minimum percent identity (min_per_id_same_path) and maximum differences between two paths (max_diffs_same_path) to merge into one were set to 90 and 10, respectively. Finally, minimum contig length (min_contig_length) used for the assembly was set to 300. We established these values based on our previous analyses to minimize transcript fragmentation. We also evaluated the integrity of assembled transcripts by comparing with Sorghum bicolor coding sequences, the distribution of contig lengths and N50 statistics. The final transcriptome comprised the longest isoform of each transcript.

To remove potential transcripts from P. kuehnii, we compared the assembled transcripts with Viridiplantae and Fungi protein sequences in the Swiss-Prot database using BLASTX  (Altschul et al. 1990). Transcripts matching Fungi sequences but not plants were discarded. Next, the remaining transcriptome was compared with the non-redundant (nr) protein database using BLASTX with an e-value cutoff of \(10^{-6}\), considering up to 20 homologous hits. The annotation was complemented with domain and protein family assignments using InterProScan and all InterPro Consortium databases. We also used Blast2GO v.5.0.8 (Conesa et al. 2005) to map and assign Gene Ontology (GO) terms to transcripts. An additional functional annotation was performed with Trinotate v.3.1.1 (Bryant et al. 2017), comparing the assembled sequences with UniProt proteins and Pfam protein domains  (Finn et al. 2016).

Differential expression and enrichment analysis

We aligned the high quality reads of each sample to the de novo assembled transcriptome using HISAT2 v.2.1.0  (Kim et al. 2015), with the option no-spliced-alignment to disable spliced alignments. We quantified gene expression levels by assigning reads to transcripts and then counting fragments with featureCounts v.1.6.0 (Liao et al. 2014), independently for each sample. Next, differential expression analysis was conducted with the edgeR package v.3.20.9 (Robinson et al. 2010). We began by selecting genes with count-per-million (CPM) greater than or equal to one in at least three samples, to filter out the lowly expressed genes. Using the full factorial model with block and time effects, we observed no significant effect of the sequencing runs. Therefore, we fitted a model considering only the time effect for all subsequent analyses. Quasi-likelihood F tests were employed to detect differentially expressed genes (DEGs)  (Lun et al. 2016). We conducted five different tests, each between a pair of adjacent time points. We also compared control plants with samples from 24 hai, checking whether the results were similar to those at 12 hai, to verify the possible influence of circadian rhythm. To minimize false positives, we used the Benjamini–Hochberg method at a specified False Discovery Rate (FDR) of 5% for each time point comparison.

We utilized the GOseq package v.1.30.0  (Young et al. 2010) to identify enriched GO terms, using genes kept after filtering for low expression as the background set. We carried out functional enrichment analyses separately for each test performed in the previous step. For all adjacent comparisons the criterion to consider a term as enriched was an FDR-adjusted p-value less than 0.05. This method enabled us to find GO terms enriched either as overrepresented or as underrepresented, i.e., terms for which the number of DEGs was greater or less than that expected by chance alone.

Graphical representation of temporal expression of annotated DEGs in groups of processes or specific pathways

Our functional annotation included information about the Enzyme Code (EC) for some transcripts. For each pairwise comparison, we built a file containing all up and downregulated transcripts that had an associated EC number. Next, we used KEGG Mapper (https://www.genome.jp/kegg/mapper.html) to search for global metabolic maps containing DEGs.

To visualize processes and pathways as informative diagrams we used MapMan v.3.6.0 (Thimm et al. 2004). First, we obtained the appropriate BIN ontologies from Mercator v.3.6 (Lohse et al. 2014), which are based on information contained in specific plant annotation databases. For doing so, we used only the genes that were differentially expressed in at least one pairwise comparison. Only DEGs with an absolute log2 Fold Change (LFC) greater than two were considered to facilitate the visualization of functional BIN members, using the expression levels of these genes in the control samples as a baseline. Finally, we selected for visualizing the diagram of biotic stress, which contains processes likely involved with the response to infection by P. kuehnii.

We browsed the functional annotation description fields to look for transcripts with potential relevance for plant signaling pathways or direct defense responses.

Expression analysis with qRT-PCR

We used qRT-PCR to validate the results obtained with RNA-Seq. The primers were designed with Primer3Plus  (Untergasser et al. 2012) and Beacon Designer™Free Edition (http://www.premierbiosoft.com), with the following parameters: (i) fragment size between 90 and 200 bp; (ii) primer size ranging from 18 to 23 bp; (iii) melting temperature between 50 and 60 \(^\circ\)C; and (iv) GC content from 40 to 60%. To verify that these primers did not spuriously amplify other transcripts, we mapped them to the transcriptome and checked that the alignments were unique.

RNA samples were treated with DNAse I (Sigma–Aldrich) and GoScript™Reverse Transcription System (Promega) was used for reverse transcription according to manufacturer instructions. qRT-PCR reactions were performed in a 7500 Fast Real-Time PCR system (Applied Biosystems) with a GoTaq® qPCR master mix kit (Promega). We used two technical replicates for each biological replicate. Reactions consisted in 6.25 \(\mu\)L of SuperMix amplification buffer, 50 nM internal standard (ROX), 0.3 \(\mu\)L of each primer at a 10 \(\mu\)M, 2 \(\mu\)L of cDNA and nuclease-free water to a final reaction volume of 12.5 \(\mu\)L. Amplification consisted in 95 \(^\circ\)C for five minutes, 40 cycles of 95 \(^\circ\)C for 10  s and 58 \(^\circ\)C for 30 s. Product specificity was evaluated by dissociation curves. As internal controls we used a known reference gene, glyceraldehyde-3-phosphate dehydrogenase (GAPDH)  (Iskandar et al. 2004), and a Proteasome subunit beta type that was discovered among the non-DEGs. Reaction efficiency and Ct values were determined with LinReg  (Ramakers et al. 2003). The relative changes in expression level were calculated with REST  (Pfaffl et al. 2002) using the same comparisons used for RNA-Seq differential expression analyses.

Results and discussion

Rapid changes in gene expression in the early hours following inoculation

We sequenced libraries from 18 leaf samples, totaling 1.17 billion reads. After quality control and trimming, 845 million high quality reads passed the filtering steps (Supplementary Table 1). The de novo assembly with filtered reads produced 451,462 transcripts with an N50 of 1201 bp, and the final assembled transcriptome consisted of the 259,804 longest isoforms, with an N50 of 672 bp. On average, 67.76% of our reads were mapped to this reference. Although we expected only a minor proportion of transcript fragments from Puccinia in comparison with sugarcane, it is possible that their frequency increased with time as in other pathosystems  (Dobon et al. 2016). Because our interest was in the expression of sugarcane genes, we removed potential P. kuehnii transcripts. We applied homology search (BLASTX) against Swiss-Prot proteins of Viridiplantae and Fungi, identifying that 92.33% of the transcripts were probable plant transcripts. They presented an N50 of 682 bp and ranged in length from 301 bp to 14.5 kbp.

For differential expression analysis, to remove lowly expressed genes, we only used those with a CPM greater than or equal to one in at least three samples. The analysis resulted in a total of 51,439 expressed genes. Next, to evaluate the presence of DEGs over time, we tested contrasts between adjacent time points. Those with the largest number of DEGs were: (i) 12 hai compared with control: 24,359 DEGs; (ii) 48 hai against 24 hai: 18,227; and (iii) 12 dai compared with 5 dai: 13,660 (Supplementary Fig. 1). These particular time points are of biological relevance, because spore germination occurs at 12 hai, and microscopic observations confirmed fungal penetration at 24 hai  (Gomez 2013). The time frame in the first 48 hours after inoculation includes the events of spore germination, establishment of the infection and hyphae development, which were crucial to understand the initial plant defense responses. The later time points examined (days after inoculation) contributed to the understanding of further metabolic responses related to the fungal colonization to allow fungal sporogenesis. Indeed, the last time point sampled (12 dai) was marked by the confirmation of sporulation through scanning electron microscopy  (Gomez 2013).

Our first strategy to compare DEGs in adjacent time points was to use the Blast2GO enrichment annotation tool to identify relevant molecular events associated with different biological moments of disease progression. We thus found over and underrepresented GO terms in enrichment tests for each comparison (Supplementary File 1). Additional enrichment tests performed with the Trinotate annotation yielded similar results (Supplementary File 2), such that we will discuss them jointly, pointing occasional differences when necessary.

Modulation of pathways related to biotic stress responses promoted by the Puccinia kuehnii infection

P. kuehnii infection caused significant changes in the expression of sugarcane genes involved in the known molecular processes in response to biotic stresses. However, our data suggest that the immune system was activated lately, two days after the inoculation (Supplementary File 3-A). Supposedly, P. kuehnii was not immediately perceived by this particular susceptible sugarcane genotype as a threat.

Sets of repressed DEGs were found in almost all of these processes at 12 hai (Fig. 1). At this time point transcripts associated with basal defenses, such as those encoding \(\beta\)-1,3-glucanases, were repressed. \(\beta\)-1,3-glucanases in combination with other antifungal proteins are involved in degradation of fungal cell walls (Supplementary File 3-H and Fig. 1). As the infection progressed more transcripts of this enzyme were detected above the basal level. The same profile was observed for transcripts involved in host cell wall reinforcements. Several genes encoding enzymes of the the phenylpropanoid biosynthesis pathway leading to lignin formation were identified as downregulated until 24 hai. Only later they were detected as upregulated (Supplementary File 3 H-K). For the GO term response to stress (GO:0006950) the transcripts with higher absolute LFC at 12 hai were downregulated (Supplementary File 3-C). The innate immune system was positively regulated two days after the inoculation (Supplementary File 3-A). GO term defense response to bacterium (GO:0042742) was among the overrepresented GO terms at 48 hai (Supplementary File 1) and response to stress was overrepresented at 5 dai. The latter included transcripts encoding proteins involved in the biosynthesis of thiamine upregulated at 12 dai (Supplementary File 3-C). Thiamine is a known activator of plant disease resistance by inducing transient expression of pathogenesis-related (PR) genes  (Ahn et al. 2005).

Assessment of plant signaling indicated that genes related to hormones commonly involved with resistance against biotrophs, abscisic acid (ABA) and SA, were repressed at 12 hai (Fig. 1). The biosynthesis of JA and ethylene were also modulated, however, some genes were up and others downregulated. Waves of expression, marked by downregulation followed by upregulation (or vice-versa) of these genes commonly associated to biotic stress response were evident during disease progression. Many coding transcripts relevant to the redox state (Figs. 1 and  2d) and most of the photosynthesis/chloroplast-related (Fig. 2c) showed changes in opposite directions in the samples collected at 12 and 48 hai. Sets of genes encoding ROS scavenging enzymes and others in response to oxidative stress were mostly variable in expression, with subgroups of similar expression (Supplementary File 3 K-N). It has been demonstrated that smut-susceptible sugarcane plants have a delay in initiating the oxidative burst after fungal infection when compared to resistant plants  (Peters et al. 2017).

Fig. 1
figure 1

Expression profiles of genes mapped to processes of the MapMan biotic stress pathway. For each functional class (indicated in bold in each ellipsis), results for 12, 24 and 48 hai are presented by adjacent heatmaps separated by dashed lines. Each square in each of the heatmaps represents a different gene mapped to the same process. The color scale indicates the \(\hbox {log}_2\) Fold Change of each gene in comparison to the baseline expression level in the control samples. Blue, white and red represent the downregulated, non-differentially expressed and upregulated genes, respectively. This figure was adapted from the MapMan representation of the biotic stress (Thimm et al. 2004). Genes within the grey area are those with experimental indication of involvement in biotic stress, while those outside the grey area are putatively involved in biotic stress pathway

Full size image

Possible late pathogen recognition by sugarcane cultivar SP89-1115

The molecular basis of plant defense mechanisms against pathogens is often explained by a two-branched immune system. The first one is based on transmembrane receptors that recognize pathogen-associated molecular patterns (PAMPs), called PAMP-triggered immunity (PTI), and the second is based on the recognition of pathogen effectors, also called effector-triggered immunity (ETI)  (Fesel and Zuccaro 2016; Jones and Dangl 2006; Hammond-Kosack and Jones 1997). Pathogen perception by the host involves the expression of resistance genes which encode proteins that typically have an NB domain and an LRR motif  (McDowell and Woffenden 2003; Jones and Dangl 2006), followed by the induction of the hypersensitive response, production of reactive oxygen species, cell wall reinforcement and hormone signaling  (Jones and Dangl 2006; Hammond-Kosack and Jones 1996, 1997).

An early fungal perception by the sugarcane plants was likely hampered. None of the GO terms usually associated with hypersensitive response (GO:0009626), defense response (GO:0006952) or response to stress (GO:0006950) were enriched until 24 hai (Supplementary File 1). In fact, enrichment analysis with the Trinotate annotation showed that plant-type hypersensitive response underrepresentation started at 12 hai and persisted until 48 hai (Supplementary File 2), and defense response was significantly underrepresented with DEGs at 48 hai (Supplementary File 2). It is possible that sugarcane cultivar SP89-1115 either was late to recognize the molecular patterns of P. kuehnii to activate defense responses, or fungal effectors may have suppressed this form of defense. Later as disease progressed, defense response to bacterium (GO:0042742) transcripts were overrepresented at 48 hai (Supplementary File 1). Although transcripts related to the immune response regulation (GO:0050776) were upregulated at 48 h, the term was not significantly enriched (Supplementary File 2 and Supplementary File 3-A). Curiously, comparing 12 dai with 5 dai, the positive regulation of defense response was underrepresented, with no DEG annotated with this GO term (Supplementary File 2).

Manual inspection of the 36 DEGs annotated with the term defense response to bacterium showed that in fact they represent general defense mechanisms that may act against fungi as well. They were related to protein-like kinases, heat-shock, ubiquitin and others (Supplementary File 3-B). Expression profiles of these genes did not show any clear pattern, except for a restoration to basal level at 48 hai. A small set of upregulated genes encoded two peptidyl-prolyl isomerase activity (PPIase) enzymes and a peroxiredoxin at this time point. Moreover, 62.82% of the transcripts annotated with the PPIase term (GO:0003755) were differentially expressed, resulting in a significant overrepresentation (Supplementary File 1). Examining the profile of these DEGs, it is clear that they were stimulated at 48 hai (Fig. 2b). PPIases were initially associated with protein folding, but some groups are immunophilins which are involved in diverse processes  (Galat 2003; Vasudevan et al. 2015). One cyclophilin, a protein with enzymatic PPIase activity, was accumulated in potato from 14 to 72 h after a combination of wounding and Fusarium solani inoculation. In that study, the expression of this gene was influenced by the presence of ABA and JA and by abiotic stresses, in a plant tissue dependent manner  (Godoy et al. 2000), which agrees with the wide range of functions of the PPIases  (Vasudevan et al. 2015). Peroxiredoxin, in turn, is a thiol peroxidase that controls redox homeostasis, photosynthesis and the transmission of redox signal  (Liebthal et al. 2018).

Fig. 2
figure 2

Time series profiles of differentially expressed genes (DEGs) in different Gene Ontology categories. The y-axis represents expression levels in normalized counts-per-million, and each solid black line corresponds to one DEG. The x-axis breaks indicate each sampled time point—T0: Control; T12: 12 hai; T24: 24 hai; T48: 48 hai; T5D: 5 dai; and T12D: 12 dai. Gene ontology categories are: (a) plant-type hypersensitive response (GO:0009626), (b) peptidyl-prolyl cis-trans isomerase activity (GO:0003755) (c) photosynthesis (GO:0015979), (d) cellular response to oxidative stress (GO:0034599)

Full size image

The heatmap of the defense response to bacterium transcripts encoding a Barwin homolog revealed upregulation from 12 to 48 hai. The Barwin domain in plant pathogenesis-related (PR) proteins is known to be involved in a response to stimulus against phytopathogenic fungi  (Li et al. 2010; Marra et al. 2006; Svensson et al. 1992). Another transcript, the 70 kDa heat shock protein (Hsp), showed significant downregulation at all time points. Members of this family are stimulated when plants are submitted to mechanical, heat or biotic stress  (Piterková et al. 2013). The chaperones Hsp70 and Hsp90 trigger the expression of PR genes and defense-related transcription factors  (Kanzaki et al. 2003).

Regarding response to stress, most differentially expressed calmodulin-binding encoding transcripts (CaM-binding) were repressed at 12 hai, with an increasing expression over time (Supplementary File 3-C). Calmodulins act in transduction of \(\hbox {Ca}^+\) signaling, playing a role in basal defense  (Reddy et al. 2003; Takabatake et al. 2007). Loss of a CaM-binding domain in an MLO (Mildew resistance locus) protein increased susceptibility of barley to powdery mildew fungus  (Kim et al. 2002). An upregulated CaM-binding gene was a member of the enriched immune system biological process in the transgenic papaya resistant to ringspot virus  (Fang et al. 2016).

Five transcripts annotated with terms thiamine and thiazole were repressed in most of the time points compared (GO:0006950 in Supplementary File 3-C). However, samples collected at 24 and 48 hai showed slightly higher expression, followed by an intense repression 5 dai, which resulted in the significant enrichment of thiazole biosynthetic terms (Supplementary File 1). Thiamine has an important role in systemic acquired resistance, which is shown to reduce disease severity over time  (Ahn et al. 2005; Kamarudin et al. 2017). In fact, both overexpression and silencing of a rice resistance gene induced a gene involved in thiamine synthesis and increased the content of this compound (Wang et al. 2016). This work revealed that resistance genes can also balance the pathways that drive thiamine biosynthesis in chloroplasts. Our results suggest that the late expression of thiamine biosynthesis genes can indicate a failure of timely resistance signaling.

We also investigated the DEGs that composed the hypersensitive and defense response terms (GO:0009626 and GO:0006952), where transcripts homologous to resistance genes were present. They were annotated as RGA (Resistance Genes Analogs), RPP (resistance to Hyaloperonospora parasitica) and RPM1 (Resistance to Pseudomonas syringae pv. maculicola 1). The proteins encoded by these genes have domains in common with those found in the NB-LRR protein family  (McDowell and Woffenden 2003). Activation of resistance genes such as NB-LRR can trigger immunity by increasing PR proteins  (Wang et al. 2016). Of the 614 transcripts annotated with defense response, 104 were differentially expressed at 12 hai. At later time points the number of DEGs for this term did not exceed 70 (Supplementary File 3-D). Additionally, half of the DEGs involved in plant-type hypersensitive response were upregulated at 12 hai, although with a low LFC (Fig. 2f and Supplementary File 3-F). Only a small portion of transcripts from these groups showed expression changes with high LFC, and the presence of many non-DEGs explained the significant underrepresentation. Given these results, we suggest that defense systems may not have been adequate to initiate a coordinated plant response against P. kuehnii.

WRKY transcription factors, cell wall reinforcement and stress-regulated transcripts altered in the first two days following inoculation

Several transcripts encoding transcription factors, physical barrier precursors or proteins that operate directly on pathogen structure can be involved in avoiding pathogen establishment. We investigated a group of transcription factors belonging to the WRKY family (Supplementary File 3-I). They are key regulators in many processes in plants and are also involved in responses to abiotic and biotic stresses  (Rushton et al. 2010). They are regulators of the defense response via the mitogen-activated protein kinase (MAPK) cascade, SA and JA pathways  (Eulgem and Somssich 2007; Pandey and Somssich 2009; Rushton et al. 2010). Genes encoding the transcription factors WRKY-75, WRKY-90, WRKY-48 and WRKY-33 were among the most upregulated in resistant plants of a papaya pathosystem  (Fang et al. 2016). In our study, the contact with the fungus repressed one of the genes encoding a WRKY-33 at 12 hai and the other two showed increased expression over time. We observed a different behavior for the genes encoding WRKY-24, WRKY-45, WRKY-51 and an uncharacterized WRKY transcription factor. The upregulation occurred since the first contact and expression levels were increased as the disease progressed (24 and 48 hai) (Supplementary File 3-I). An allele of WRKY-45 was reported as a negative regulator of bacterial infection in rice, while the other allele conferred resistance  (Tao et al. 2009). Interestingly, both alleles positively regulated the resistance against the fungus Magnaporthe grisea  (Tao et al. 2009).

Disease resistance can also be initially triggered by genes that act over pathogen structure after PAMP recognition, such as \(\beta\)-1,3-glucanase and chitinase  (Balasubramanian et al. 2012; Fesel and Zuccaro 2016). We found two groups of \(\beta\)-1,3-glucanases with different expression patterns, particularly during the first 48 hai. The first group contained downregulated genes, with some members reaching a minimum LFC of − 2.0 at 12 hai and 24 hai. The second group showed LFC values up to 5.0 at 24 hai, indicating a strong stimulus to expression (Supplementary File 3-H). These enzymes are hydrolases and are classified as PR because they act by directly degrading \(\beta\)-1,3-glucans, the most abundant polysaccharides present in fungal cell walls, and that serve as PAMPs  (Balasubramanian et al. 2012; Fesel and Zuccaro 2016). A \(\beta\)-1,3-glucanase 5 was highly upregulated at 3 hai in a peach cultivar resistant to Xanthomonas arboricola pv. pruni  (Gervasi et al. 2018). During smut infection in sugarcane, RNA-Seq analysis revealed an induced \(\beta\)-1,3-glucanase gene at 48 hai (McNeil et al. 2018). Additionally, their qRT-PCR analyses revealed \(\beta\)-1,3-glucanase repression in a susceptible sugarcane cultivar, while a cultivar with external resistance exhibited upregulation.

A widely studied fungal PAMP is chitin, which is a structural polymer that composes part of the cell wall and is degraded by another kind of hydrolase, the chitinase  (Fesel and Zuccaro 2016). Supplementary File 3-J shows that eight chitinase-annotated transcripts presented a pattern similar to that of the second group of \(\beta\)-1,3-glucanases, with expression levels always above the basal level. Only three genes showed expression levels below that observed in the control samples. Others presented varying expression profiles across time points. Chitinase genes were reported as differentially expressed even in sugarcane cultivars susceptible to smut (Que et al. 2014; Su et al. 2015). A similar mixed expression pattern was also described in these cultivars, while in the case of resistance they were upregulated in a tissue-dependent manner at 24 or 48 hai  (Su et al. 2015). Other studies reported chitinases as primary temporal responses in apple and peach plants infected with Alternaria alternata and Xanthomonas arboricola pv. pruni, respectively (Socquet-Juglard et al. 2013; Zhu et al. 2017), differentially expressed in coconut infected with phytoplasma  (Nejat et al. 2015) and in the cacao-Moniliophthora perniciosa pathosystem  (Teixeira et al. 2014). Genes involved in the synthesis of this enzyme were first induced 1 dai and most prominently 11 dai in wheat susceptible to yellow rust  (Dobon et al. 2016).

Among the physical barriers against infection, lignin deposition in the cell wall is considered a localized response that can be hampered in susceptible plants as a result of pathogen modulation  (Nicholson and Hammerschmidt 1992). Transcripts annotated with the lignin biosynthetic process term (GO:0009809) showed a substantial decrease in expression in the first time point after inoculation, with some reestablishment at 48 hai (Fig. 3a). Lignin deposition in the epidermis is an important event to prevent sugarcane infection by S. scitamineum at 24 hai  (Marques et al. 2018). We observed that almost all genes of cinnamyl-alcohol dehydrogenase (CAD), a precursor of lignin biosynthesis, returned to their basal expression level after marked downregulation at 12 hai. In contrast, it is noticeable that three genes were highly upregulated from 12 hai to 12 dai (Fig. 3b). To achieve an overview of the phenylpropanoid biosynthesis, we selected DEGs annotated as enzymes of this pathway to visualize their expression profiles. Among them there were a CAD (EC:1.1.1.195) and a caffeic acid 3-O-methyltransferase (EC:2.1.1.68) upregulated in the transition from 12 to 24 hai (Fig. 3d). Also, Puccinia induced all phenylalanine/tyrosine ammonia-lyase (PTAL, EC:4.3.1.25), four 4-coumarate-CoA ligase (EC:6.2.1.12) and all cinnamoyl-CoA reductase (EC:1.2.1.44) 48 hai. Curiously, upregulation of the bifunctional PTAL instead of the phenylalanine ammonia-lyase (PAL, EC:4.3.1.24) and tyrosine accumulation after whip emission in the smutted cane is a hallmark of whip emission, probably associated with tissue lignification  (Schaker et al. 2017). Various studies describe the enrichment of differentially expressed genes involved in the pathway of phenylpropanoids with diverse outcomes  (McNeil et al. 2018; Santa Brigida et al. 2016; Xia et al. 2017; Zhu et al. 2017). For instance, poplar infected by a necrotroph shows a temporary activation of 4-coumarate-CoA ligase and cinnamoyl-CoA reductase at 36 hai  (Zhang et al. 2018). Other authors observed that upregulation of these genes could be associated with resistance  (Fang et al. 2016; Que et al. 2014). Otherwise, the compatible interaction of an apple cultivar and A. alternata resulted in phenylpropanoid biosynthesis activation at 12 hai  (Zhu et al. 2017). Our results indicated that the last component of the pathway, peroxidase (EC:1.11.1.7), was represented by groups of genes with variable expression, as well as some with a more constant expression pattern along the investigated time interval (Fig. 3d and Supplementary File 3-K). This pattern of expression can be partly due to the many functions attributed to the enzyme  (Passardi et al. 2005).

The phenylpropanoid pathway is commonly associated with plant defense, not only by providing molecules for lignification but also antimicrobial compounds such as coumarins and flavonoids  (Naoumkina et al. 2010). One well-studied enzyme that is at the onset of the biosynthesis of these components is PAL. High expression of PAL-encoding genes is a common attribute of resistant cultivars  (Loarce et al. 2016; Gervasi et al. 2018), including those of sugarcane  (McNeil et al. 2018; Que et al. 2014; Santa Brigida et al. 2016). In our results, we identified that P. kuehnii repressed initially (12 hai) PAL DEGs in sugarcane; three of them were induced in the following hours, returning to the basal level at 5 dai (Fig. 3c, d). We suggest the involvement of both PAL and PTAL in infected sugarcane plants, which although activated, was not timely enough to control the Puccinia colonization.

Fig. 3
figure 3

Genes in the phenylpropanoid biosynthesis pathway. (a) Expression profile of genes annotated with term lignin biosynthetic process (GO:0009809); (b) Heatmap of genes annotated as cinnamyl-alcohol dehydrogenase (CAD); (c) Heatmap of genes annotated as phenylalanine ammonia-lyase (PAL); (d) Simplified representation of the KEGG phenylpropanoid biosynthesis pathway with differentially expressed Enzyme Codes (ECs). Only ECs containing differentially expressed genes are shown in the pathway. The expression profiles of the genes annotated with these ECs are shown in the heatmaps on the sides, with colors matching the annotation in the pathway. Heatmap cells in (b), (c) and (d) present the \(\hbox {log}_2\) Fold Change for each time point—indicated in the columns—compared to control samples. In (d), the ECs correspond to: EC:4.3.1.24— phenylalanine ammonia-lyase; EC:4.3.1.25— phenylalanine/tyrosine ammonia-lyase; EC:6.2.1.12— 4-coumarate-CoA ligase; EC:1.2.1.44—cinnamoyl-CoA reductase; EC:1.1.1.195—cinnamyl-alcohol dehydrogenase; EC:3.2.1.21—beta-glucosidase; EC:1.11.1.7— peroxidase; EC:2.1.1.68—caffeic acid 3-O-methyltransferase; EC:2.3.1.133—shikimate O-hydroxycinnamoyltransferase; EC:1.14.14.91— trans-cinnamate 4-monooxygenase; EC:2.3.1.99—quinate O-hydroxycinnamoyltransferase

Full size image

Examination of other enriched terms revealed additional information relevant to stress responses. We identified at 48 hai terms related to transposable elements (TE) and the regulation of immune response as highly induced when compared to any other time point. The transposition related term with the lowest p-value was GO:0032199 (reverse transcription involved in RNA-mediated transposition), detected in homologs of LINE-1 retrotransposons, endonucleases, exonucleases and reverse transcriptases. Trinotate also detected TE as relevant at this same time after inoculation (48 hai) revealing terms such as RNA-mediated (GO:0032197) and viral processes (GO:0044826 and GO:0075713) (Supplementary File 2). TEs compose approximately 50% of the sugarcane genome (de Setta et al. 2014), are tightly regulated (De Araujo et al. 2005) and were previously shown as differentially expressed in sugarcane interacting with symbionts, pathogenic virus and other fungi (Belarmino et al. 2013; Medeiros et al. 2014; Thokoane and Rutherford 2001). Plants commonly activated TE gene expression in stress conditions, including when infected with fungus (Melayah et al. 2001). In rice, for instance, more than 10% of the transposable elements showed changes in expression in seedlings infected by the stripe virus (Cho et al. 2015). The infection by and development of P. kuehnii in sugarcane cultivar SP89-1115 did not seem to lead to a coordinated expression of these genes until 24 hai, which was remarkably changed 48 hai (Supplementary File 3-G).

Expression waves of genes involved in photosynthesis and oxidative processes

Puccinia infection also modulated sugarcane genes associated with photosynthesis and metabolism of reactive oxygen species (ROS) (Supplementary File 1). Changes in the primary plant metabolism of plants, including those related to carbohydrate partition, are a consequence of biotic stresses and necessary to meet the nutritional requirements of the pathogen  (Berger et al. 2007). All the genes associated with the term carbon utilization (GO:0015976) expressed differentially at 12 hai, simultaneously with transcripts encoding proteins of the photosynthetic apparatus. Out of 76 DEGs annotated with photosynthesis terms (GO:0015979 and GO:0009765), sugarcane repressed 77.63% at 12 hai. Later (24 hai), 30 of these genes were upregulated. Additionally, 48 hai showed an increased expression of 64 genes. Downregulation of genes associated with photosynthesis is a general phenomenon observed in defense response shortly after infection  (Berger et al. 2007; Bilgin et al. 2010; Gervasi et al. 2018; Nejat et al. 2015; Teixeira et al. 2014). We identified the same response in our data, where repression occurred for the majority of genes 12 hai, then later (48 hai) restored the basal expression levels (Fig. 2c).

Because circadian rhythm could potentially influence the observed repression, we investigated the contrast between 24 hai and the basal expression level in the control. This comparison revealed 9978 downregulated genes, of which 88.34% were also downregulated at 12 hai. Only 30 DEGs were upregulated at 12 hai and downregulated at 24 hai against control. We identified 84 DEGs annotated with the photosynthesis terms at 24 hai against control, 89.29% of which with the same profile as at 12 hai.

Terms associated with chlorophyll, chloroplast and photosystems were highly represented in our observations from 12, 24 and 48 hai (Supplementary File 1). Various studies demonstrated the negative impact of pathogen infection in the expression of photosynthesis-related genes (Bilgin et al. 2010). The physiological impact of orange rust measured by chlorophyll content and net photosynthetic rate of sugarcane leaves revealed a reduction according to susceptibility to P. kuehnii (Zhao et al. 2011). Other studies, such as gene expression analysis using microarray in Arabidopsis infected by P. syringae, revealed the repression of nuclear genes with chloroplast functions in the early stages of infection (de Torres Zabala et al. 2015). From 6 to 8 hai there was a decrease in photosynthetic \(\hbox {CO}_2\) accumulation, reducing net photosynthesis in comparison to mock and the infection with a non-virulent P. syringae strain (de Torres Zabala et al. 2015). Similar responses identified in susceptible cucumber and rice plants have also shown repression of genes associated with chlorophyll, chloroplast and photosynthesis functions (Ahmed et al. 2017; Xia et al. 2017).

One of the most consistent defense responses activated in plants resistant to pathogens is the oxidative burst. ROS can either directly suppress pathogen growth or act as signal molecules to induce defense mechanisms such as the expression of plant pathogenesis-related proteins, activate localized cell death and cell wall reinforcement (Demidchik 2015; Inzé and Montagu 1995; Lehmann et al. 2015; Mittler et al. 2011; Santa Brigida et al. 2016). In all the time points evaluated in our experiment, P. kuehnii altered the overall balance of redox reactions in sugarcane. We identified enrichment of the terms GO:0055114; GO:0016491; GO:0016705; GO:0045454 in genes up and downregulated in inoculated plants (Supplementary File 1). Oxidoreductase activity (GO:0016491 and GO:0016705) included 743 downregulated genes at 12 hai, 588 of them upregulated at 48 hai. Redox homeostasis (GO:0045454) represented by redoxins were downregulated at 12 and upregulated at 24 hai.

We identified other possible connections between the oxidative pathways and the results discussed previously. Examples of these are the involvement of ROS with Hsp70 accumulation  (Piterková et al. 2013) and for the activity of calmodulins, amine oxidases and other types of peroxidases  (Demidchik 2015; Liebthal et al. 2018; Gervasi et al. 2018; He et al. 2018). Amine oxidases were detected as repressed 12 hai in our analysis showed an LFC less than − 2.0 until 24 hai (Supplementary File 3-O). The expression profiles of genes associated with cellular response to the oxidative burst followed a pattern similar to that for genes involved in photosynthesis described before (Fig. 2d). Approximately 70% of all the genes in this category differed in expression 48 hai (Supplementary File 1). The term response to oxidative stress (GO:0006979) composed mainly by peroxidases was assigned to genes up and downregulated in about equal amounts (Supplementary File 3-E). To protect their cells, in response to ROS, plants induce enzymes such as catalases (CAT), superoxide dismutases (SOD), and glutathione transferases (GST) (Demidchik 2015). We did not detect a common expression pattern of genes encoding these enzymes in our study (Supplementary File 3 L-N). Our data suggests that ROS may be involved in signaling rather than directly affecting the pathogen, and that this response was late probably because of the susceptibility of the sugarcane genotype. Similar results were described for sugarcane infected with smut (Peters et al. 2017).

qRT-PCR validation of differentially expressed genes

We selected 11 transcripts for validation based on their differential expression analysis and their biological relevance to the infection process. They were selected among the genes related to lignin biosynthesis, hormones, transcription factors, response to stress, response to oxidative stress and photosynthesis. The first group of genes showed significant differential expression at all time points and were functionally assigned as orthologs of probable cinnamyl alcohol dehydrogenase 8D, calmodulin binding protein2, putative ethylene response sensor 2 and glutathione transferase. The second group presented differential expression in at least two time points and were functionally assigned as orthologs of catalase, probable WRKY transcription factor 70, Superoxide dismutase, small heat shock protein, Glucan endo-1,3-beta-glucosidase, bZIP transcription factor and chloroplast ribulose-1,5-bisphosphate carboxylase/oxygenase small subunit. The only non-DEG selected was a Proteasome subunit beta type, used as internal control along with GAPDH. Each primer pair was designed to ensure effective amplification of a unique transcript (Supplementary Table 2). Overall, we performed 55 experiments. Of these, only seven did not confirm the RNA-Seq differential expression results. The correlation of both evaluation strategies was high (\(r=0.85\)), showing good correspondence between RNA-Seq and qRT-PCR results (Supplementary Figs. 2 and 3).

We only detected two situations in disagreement between both techniques. The first one was the detection of significance in only one technique, and second, we detected significant changes in opposite directions for each technique. For the first situation, RNA-Seq results detected significant small fold changes for BZIP 12 hai, and SOD and GST 24 hai, not detectable in our qPCR analysis. Otherwise, qPCR detected a significant variation of GLU 48 hai. We believe the latter case was due to the variation of CPMs among biological replicates, which resulted in higher residual variance. The second situation was seen for HEAT at 48 hai, which showed repression using RNA-Seq data and upregulation using RT-qPCR. Both the forward and reverse primers aligned to four isoforms of this gene, but only the forward primer aligned against the longest isoform. There is a possibility that these primers amplified other heat shock transcripts.

Successful infection of susceptible sugarcane genotype

Biotrophic fungi typically germinate on the leaf surface, form appressoria to invade the host, feed from the host cells using haustoria and propagate producing numerous uredospores  (Lo Presti et al. 2015). Oppositely, plants avoid tissue colonization using the first line of defense mechanism (PTI) responsible for the recognition of pathogen molecular patterns (PAMPs), resulting in the activation of defense responses. A second and more specific defense mechanism, ETI, occurs when resistance proteins recognize pathogen effectors resulting in unsuccessful tissue colonization (Jones and Dangl 2006). One of the most used strategies to capture the dynamics of these molecular events is to detect changes in gene expression following the progress of plant infection (Loarce et al. 2016). In the sugarcane-P. kuehnii pathosystem, appressoria and pustules were observed (Gomez 2013) at 24 hai and 12 dai, respectively, and symptoms detected 9 dai, following sporulation two days later. Our analyses revealed specific genes induced at the beginning of the infection process. Among these, we identified genes involved in cell wall reinforcement (lignification), and encoding chitinases and other PR proteins. These results suggested that sugarcane can recognize the presence of the microorganism. However, we also identified repressed genes indicating an inefficient and/or late plant counteraction. Some of these are genes encoding proteins related to the hypersensitive and defense responses, CaM-binding, photosynthesis-related, PPIases, \(\beta\)-1,3-glucanase, thiamine-thiazole biosynthesis, PALs and CADs.

The description of waves of expression in the various molecular events following the interaction was a hallmark in our work. The overall impact of P. kuehnii on the cultivar SP89-1115 involved recognition of the pathogen shortly after inoculation (12 hai), repressing the expression of genes involved in this biotic stress and a gradual return to basal levels (24 hai) (Fig. 1). Next, 24 hai we identified repression of SA potentially favoring pathogen growth, and again a return to normal levels 48 hai. Biotrophic pathogens, in turn, change the signals to avoid host cell death (Kazan and Lyons 2014). For instance, the biotrophic fungus Ustilago maydis secretes the chorismate mutase into the plant cell, modulating the SA pathway during the infection in maize  (Djamei et al. 2011; Patkar and Naqvi 2017). Genes related to ABA metabolism were downregulated in our analysis at 12 hai and remained repressed at 24 hai (Fig. 1, Supplementary File 2). ABA induces resistance favoring signal transduction associated with defense response (Fang et al. 2016; Loarce et al. 2016). These modulations of plant hormones are usually known as crucial to establishing defense signaling against many biotrophs (Glazebrook 2005; O’Connell and Panstruga 2006).

Early responses are essential to initiate a successful defense strategy in sugarcane resulting in resistance (McNeil et al. 2018; Peters et al. 2017; Que et al. 2014). Our study unveiled two possible causes associated with SP89-1115 susceptibility, the inefficiency of the defense response, and the effectivity of P. kuehnii to disable the existing defense mechanisms after initial contact.