Transcriptome Analysis of Poplar during Leaf Spot Infection with Sphaerulina spp.

Diseases of poplar caused by the native fungal pathogen Sphaerulina musiva and related species are of growing concern, particularly with the increasing interest in intensive poplar plantations to meet growing energy demands. Sphaerulina musiva is able to cause infection on leaves, resulting in defoliation and canker formation on stems. To gain a greater understanding of the different responses of poplar species to infection caused by the naturally co-evolved Sphaerulina species, RNA-seq was conducted on leaves of Populus deltoides, P. balsamifera and P. tremuloides infected with S. musiva, S. populicola and a new undescribed species, Ston1, respectively. The experiment was designed to contain the pathogen in a laboratory environment, while replicating disease development in commercial plantations. Following inoculation, trees were monitored for disease symptoms, pathogen growth and host responses. Genes involved in phenylpropanoid, terpenoid and flavonoid biosynthesis were generally upregulated in P. balsamifera and P. tremuloides, while cell wall modification appears to play an important role in the defense of P. deltoides. Poplar defensive genes were expressed early in P. balsamifera and P. tremuloides, but their expression was delayed in P. deltoides, which correlated with the rate of disease symptoms development. Also, severe infection in P. balsamifera led to leaf abscission. This data gives an insight into the large differences in timing and expression of genes between poplar species being attacked by their associated Sphaerulina pathogen.


Introduction
There are numerous important diseases caused by members of the fungal family Mycosphaerellaceae that affect plants, and forest trees are no exception. Many Populus species and hybrids are susceptible to leaf and stem diseases caused by Sphaerulina species including leaf spot of balsam poplar (P. balsamifera), eastern cottonwood (P. deltoides) and black poplar (P. nigra) caused by S. populicola, S. musiva and S. populi, respectively, and stem cankers of hybrid poplars caused by S. musiva. Diseases caused by S. musiva were reported to reduce tree growth through premature defoliation during severe leaf spot infections [1][2][3]. Canker infections, caused by S. musiva are considered to be the most severe diseases in hybrid poplar plantations of North America [1][2][3][4][5]. Additionally, new and yet to be described species of 24 h at 4°C prior to planting in high density trays filled with Pro-Mix BX. After 6 weeks, seedlings were transferred to 3 L pots and grown for a further 3 months in the greenhouse. The 5-month-old trees were transferred to Conviron growth chambers set at 22°C and 16 h light with a photosynthetic flux of 90-100 μmol m -2 s -1 and incubated for 4 weeks to reach a final age of 6 months.
One monosporal Sphaerulina isolate from each species was utilized. Isolates were prepared from spore suspensions stored in 25% glycerol at -80°C. All isolates were originally obtained from pycnidia emerging from leaf spots, but were collected at different locations and at different times: the S. musiva isolate (ND05.02) was originally obtained from a P. deltoides leaf in Notre-Dame du Nord, Quebec, Canada in 2002 [21]. S. populicola isolate 02.72A was isolated from P. trichocarpa in Marion County, Oregon, USA in 2002 [21] and the uncharacterized Ston1 isolate was obtained from leaf spots on P. tremuloides in Stoneham, Québec, Canada in 2003 [6]. All fungal isolates were confirmed to be members of their species by PCR-based detection [6,21]. To obtain conidia for foliar application, storage cultures were seeded to potato dextrose agar (PDA; Difco 213400) plates and grown at 22°C in the dark for 8 days. Conidia were harvested by flooding the plates (20 per isolate) with sterile distilled water, then filtered through 100-μl nylon mesh sterile Falcon 1 cell strainers (Fisher Scientific, Waltham, MA, USA), centrifuged at 4000 x g for 10 min, washed with 40 ml of sterile distilled water and resuspended in sterile distilled water to 1 x 10 5 conidia L -1 prior to application.

Experimental design and measurements
Eighteen clones for each tree species were prepared for experimentation. Trees from each species were separated into different growth chambers to prevent cross contamination by different Sphaerulina species. Using a Crown spray tool (North American Professional Products, Vaughan, Ontario, Canada), each tree was sprayed with 10 mL of atomized spore suspension (1 x 10 5 conidia mL -1 ). Three trees from each species were sprayed with water to serve as controls. To further prevent spore cross contamination, inoculation of each species was done on separate days. The leaves from each tree were harvested following a timeline such that trees were only sampled once. The top three fully expanded leaves were excised, pooled and flash frozen in liquid nitrogen for storage prior to DNA and RNA extraction. Samples were taken 0, 1, 4, 8, 15 and 30 days after inoculation (DAI). All leaves on trees were monitored daily for signs of disease such as the occurrence of leaf spots, leaf abscission and plant height. Disease was scored based on a four-point scoring system for leaf spots on each tree prior to harvest, where 1 represents 1-5% of leaf area, 2 represents 6-30% of leaf area and 3 represents over 30% of leaf area. A score of 0 was assigned when no leaf disease symptoms were observed. Differences in plant height and in the number of leaves were statistically analyzed at 35 DAI by analysis of variance (ANOVA) and Tukey's pair-wise HSD test with JMP 8.0 (SAS Institute Inc., Cary, North Carolina, USA).

DNA and RNA extraction
Both DNA and RNA were extracted from the leaves following inoculation. DNA was used to track the growth of the fungi through quantitative PCR (qPCR), while RNA was used for transcriptome analysis and reverse transcript-qPCR (RT-qPCR). Leaves were homogenized in a liquid nitrogen chilled mortar and pestle, and were then stored at -80°C to be used for both RNA and DNA extraction. DNA was isolated from 100 mg of tissue for each sample using the DNeasy Plant Mini Kit (Qiagen, Toronto, Ontario, Canada). RNA was also isolated from 100 mg of tissue using the RNeasy Plant Mini Kit (Qiagen) with in-column RNase-Free DNase Set (Qiagen) treatment.

Transcriptome sequencing
RNA was extracted from each tree sampled and assayed for quality with a 2100 Bioanalyzer (Agilent, Santa Clara, California, USA). All RNA samples had a ratio of absorbance at 260 and 280nm (A260/280) above 2.0 and no peaks indicating the presence of genomic DNA were detected. When further analyzed using a 2100 Bioanalyzer, all samples had a RNA integrity number (RIN) greater than 7. Libraries for sequencing were prepared by Brian Boyle at the Institut de biologie intégrative et des systèmes (IBIS) of Université de Laval, Quebec, Quebec, Canada. Each cDNA library was barcoded and six samples were pooled together. In total, six lanes of pooled samples were single-end sequenced using an Illumina HiSeq 2500 platform at the McGill University and Génome Québec Innovation Centre (Montreal, Quebec, Canada).

RNA-seq pipeline
The sequencing libraries obtained from Génome Québec were analyzed using RobiNA software version 1.2.4 (build656) pipeline [22]. First, the data was quality checked for base call quality, homopolymers, kmer frequency, base call frequency and overrepresented sequences. Data was then trimmed by removing the adapter sequences, low quality reads as well as reads shorter than 20 base pairs. The data were then aligned separately to the P. trichocarpa V3.0 210 genome build using BOWTIE [23] set to a seed length of 28 base pairs allowing up to two mismatches. Unaligned reads were saved for further analysis. The mapped libraries were then analyzed for digital gene expression (DGE) using edgeR [24] with a Benjamini-Hochberg correction and false discovery rate (FDR) threshold of 0.05.
Gene annotations for significantly expressed poplar genes with a log2 fold change > 1/-1 were obtained from the P. trichocarpa V3.0 210 peptide data set. Significantly expressed poplar genes were functionally annotated using BLAST2GO version 2.8 [25]. Fisher's exact test was performed on a list of significantly differentially expressed genes at each time point with a log2 fold changes > 2.0. GO terms were screened to generate the most specific terms with a FDR < 0.05. Benjamini-Hochberg multiple testing correction of FDR was also performed. MapMan version 3.5.1 was used to display differentially expressed poplar genes (log2 fold change > 1/-1) according to their metabolic pathways. [26]. The Wilcoxon rank sum test was used to identify groups of genes in specific functional bins whose gene expression had greater variation then the entire dataset [26].
Raw RNA-seq files as well as read count and digital gene expression analysis matrixes for P. trichocarpa primary genes are publicly available on NCBI Gene Expression Omnibus (GEO) (accession number: GSE67697).

Quantitative PCR
To determine the relative abundance of fungal DNA within the leaf samples, DNA was extracted from the leaves of inoculated plants and qPCR was used to amplify and determine the abundance of the actin-1 fungal gene (S. musiva transcript 147624, S. populicola transcript 133327 and Ston1 NCBI genome: 33152, Contig00301:46340-46570) relative to poplar α-tubulin-1 gene (Potri.002G111900.1). The Ston1 actin-1 gene was identified as the top BLASTN match to S. musiva transcript 147624. Differences in relative fungal abundance between time points was statistically analyzed by ANOVA and Tukey's pair-wise HSD test.
RT-qPCR was used to validate the expression profile of 13 poplar genes identified by the RNA-seq analyses. Ten genes with induced expression in one or more of the host-pathogen systems following inoculation were chosen, while three genes with stable expression were chosen to represent non-induced reference genes. Prior to RT-qPCR, cDNA libraries were generated from each RNA sample using the QuantiTect Reverse Transcription Kit (Qiagen).
Primers were designed using the Oligo Analyzer/Oligo Explorer program (http://www. genelink.com/tools/gl-oe.asp). Annealing temperature (Tm) values were set at 65°C using the following parameters: 50 mM NaCl concentration and 250 pM DNA. No template control reactions were run on the primer pairs to detect dimer formation. A list of primers used in this study is shown in S1 Table. Gene expression was analyzed for each of the sample using the Light Cycler 480 (Roche). All reactions were performed in a final volume of 10 μl and contained 1X QuantiTect SYBR 1 Green PCR Master Mix (Qiagen), 0.5 μM of each primer, and 1 μl of template DNA, which is equivalent to 10 ng of total RNA. PCR thermocycling conditions were set at 95°C for 15 min, 50 cycles at 95°C for 15 sec, 58°C for 30 sec, and 65°C for 90 sec. Fluorescent readings were taken at the end of each cycle, and the specificity of amplification as well as the absence of primer dimers were confirmed with a melting curve analysis at the end of each reaction. To correct for technical variation in RNA extraction, total RNA quantification, reverse transcription, and RT-qPCR reactions as well as biological variation expression data were normalized against the geometric mean of three reference genes, UBi 10 (Potri.014G115100.1), Eif 4 (Potri.006G225700.1), and MKK 2 (Potri.018G050800.1), by geNORM VBA applet for Microsoft Excel (http://medgen.ugent.be/jvdesomp/genorm) [27]. For all experiments, fold change is expressed as treatment relative to the control calculated using the 2 -ΔΔCt method [28]. Standard deviation related to the biological variation within one line was calculated in accordance with error propagation rules. Statistical analysis of gene expression between 1 and 15 DAI was conducted using Student's paired t-test at p < 0.05.

Results
RNA-seq was successfully used to associate disease symptoms with molecular plant responses. Large differences were observed between the different plant pathogen interactions regarding the rate of disease development, observed phenotypes during infection and molecular pathways significantly regulated during infection.

Disease development
Leaf spots developed in all poplar-Sphaerulina pathosystems, but with distinct temporal and spatial patterns within individual trees and between tree host species. Leaf spots occurred earlier and more severely on lower, more mature leaves after inoculation. During inoculation, stems were also coated in atomized conidia; however, no stem disease was observed in any pathosystem. The level of disease severity between the pathosystems was different. Fig 1 shows both healthy non-inoculated leaf and an infected leaf 30 DAI. The occurence of first disease symptoms varied between pathosystems. In P. deltoides, leaf spots were first observed 14 DAI as minute localized discolouration of the epidermis on the adaxial leaf surface. In P. balsamifera and P. tremuloides, the visual symptoms appeared earlier, at 8 and 9 DAI, respectively. The disease severity scores were relatively low for P. deltoides with 1 and P. tremuloides with 2. Leaf spot symptoms were much greater in P. balsamifera and disease severity scores reached the maximum on the four-point scale (Fig 2). Pathogen growth was quantified by qPCR using genomic DNA extracted from infected leaves (Fig 2). In P. deltoides and P. tremuloides, the amount of fungal DNA significantly decreased from the initial spore application at later time points (Fig 2). In P. balsamifera, pathogen growth did not significantly increase from 1 to 15 DAI, but significantly increased by 30 DAI (Fig 2).
The lowest leaves on all inoculated P. balsamifera clones began to abscise at 20 DAI, which was 12 days after appearance of the first leaf spots. By 35 DAI, the majority of leaves on the inoculated plants had abscised (Fig 3A). At 35 DAI, there was a significant difference in the number of leaves in the control vs inoculated plants ( Fig 3B). No significant leaf abscission was observed in the other two pathosystems. No statistical differences in plant height at 35 DAI were detected between inoculated and controls plants in all tree species (Fig 3C). Leaf abscission appears to have had no effect on plant growth, but this may be due to the short nature of this 35 day study. The timing of leaf harvest had no significant effect on leaf number or plant height at 35 DAI.

Mapping statistics of transcriptome analysis
In total, over 100 billion base pair (bp) of sequenced information in 1.06 billion 100 bp reads was obtained from all combined libraries (Table 1). Each individual library ranged from 9.5-46 million reads. These reads were mapped to the P. trichocarpa genome version 3. Summary mapping statistics for each time point in each pathosystem are summarized in Table 1, and each individual library is shown in S2 Table. A significant difference in read mapping was observed between pairwise analyses of read mapping from each poplar species. Reads generated from P. balsamifera had the highest percentage of mapping to the poplar genome at 55.5%, followed by 53.5% for P. deltoides and 42.6% for P. tremuloides. As expected, read mapping percentage was highest with the closest related poplar species, P. balsamifera and lowest with the more distantly related species to P. trichocarpa [29]. Mapping variation between biological replicates in each pathosystem is shown in S1 Fig, where the lowest variation between biological replicates at each time point was observed in P. balsamifera. An average of 24947 ± 720 transcripts were mapped with a minimum sum of 10 reads from replicates for all species and time points with no significant variation in the number of genes. The poplar genome is known to be highly heterozygous [30], and these variations could also have an effect on read mapping.
Two individual libraries from P. tremuloides had poor read mapping statistics, with only 33.1% and 7.5% of reads successfully mapping at 4 DAI and 15 DAI, respectively. These two libraries both had consistent base quality scores above 30 for the majority of reads; however, a kmer enrichment test revealed different conserved 5 bp sequences present at multiple locations reaching an enrichment factor up to 10%. Due to this, potential library construction or sequencing error resulted in 40% and 80% of the reads in one sample at 4 DAI and 15 DAI to be excluded from mapping. Despite the reduction of raw reads, 12 million and 1.8 million reads were still successfully and accurately mapped. These datasets were included in the downstream statistical analysis, though this resulted in a decreased confidence in gene detection due to the substantially variation in read depth (S2 Table).

Differential gene expression analysis
Species and time specific sets of genes are differentially expressed in each pathosystem. In total, over 5000 poplar genes were identified with differential expression between inoculated and healthy tissues. Contrasts for poplar gene mapping were conducted between non-inoculated healthy plants and infected plants at 1, 4 and 15 DAI. The Venn diagrams shown in S2 and S3 Figs demonstrate the unique and shared differentially expressed genes at each time point and between each pathosystem. In all poplar species, the number of observed differentially expressed genes was higher later during infection (S3 Fig). The number of differentially expressed genes identified was higher at 4 DAI than at 1 DAI in P. tremuloides; however, in the other two plant species, there were fewer differentially expressed genes observed at this midpoint. At 15 DAI, over 1000 upregulated and downregulated genes were identified at each time point and significantly regulated genes were shared between different time points. A number of genes were commonly differentially expressed between host species at various time points (S3 Fig).

Functional annotation of differentially expressed genes
Different functionally annotated genes and gene families were expressed at different time points between the three pathosystems. Pathway analysis focused on the examination of secondary metabolic pathways to assess the major host defensive responses. MapMan graphical secondary metabolic pathway was used to assess the plant response to fungal inoculation ( Table 2). Wilcoxon rank sum test was used to identify significantly regulated pathways at each time point. Gene expression data for each significantly regulated pathway are shown in the heat maps of Fig 4 for each species, where each box represents a gene annotating to the specific bin ID to the given pathway. Only genes that annotated to MapMan secondary metabolite bins are shown and only a representative fraction of lignin biosynthesis genes are shown. Overall, genes involved with anthocyanin and simple phenol biosynthesis are generally downregulated, while genes of the terpenoid, lignin or flavonoid biosynthesis are generally upregulated. In P. deltoides, which showed the lowest level of leaf spot infection, the significant regulation in simple phenol biosynthesis pathways was the only observed change occurring before leaf spot emergence. After leaf spot formation at 15 DAI, genes involved in lignin biosynthesis were observed ( Table 2).
In the more moderately infected P. tremuloides trees, terpenoid biosynthesis genes, including those in the mevalonate pathway that produces the isoprenoid precursors needed for  terpenoid biosynthesis, were upregulated at 1 DAI (Table 2). However, these pathways were no longer observed to be differentially regulated midway through infection at 4 DAI, when fungal biomass was decreasing. After leaf spot formation, the mevalonate pathway was again activated along with the biosynthesis of sulfur containing secondary metabolites, while genes involved in anthocyanin biosynthesis were generally downregulated. Terpenoid biosynthesis gene upregulation was also observed in the P. balsamifera clones, but with a different timing of expression than in P. tremuloides. The terpenoid pathway was differentially regulated at 4 DAI until leaf spot formation, but was not observed early during infection when only the general downregulation of anthocyanin biosynthesis was observed (Table 2). Also, flavonoid chalcone biosynthesis was significantly regulated both before and after leaf spot formation in correlation with terpenoid biosynthesis. P. balsamifera clones were highly susceptible to leaf spot despite the expression of these genes; therefore, the differential expression of genes involved in terpenoid and anthocyanin biosynthesis may indicate an insufficient defense response to S. populicola.
Fisher's exact test was used successfully to identify specific GO terms at different time points during leaf spot infection. Similar results were observed between MapMan and BLAST2GO analysis. No specific GO terms were overrepresented in P. deltoides at 1 DAI or 4 DAI.  (S4-S6 Figs). The respective pathogen appeared to be detected soon after inoculation in P. tremuloides and P. balsamifera, as a number of GO terms associated with defense responses to fungi, innate immune responses, defense involving deposition of callose, oxidation-reduction processes and transcription factors are overrepresented at 1 DAI. Few or no GO terms were observed to be overrepresented at 4 DAI in all pathosystems. However, at 15 DAI, all poplar species were expressing a number of fungal defensive response reactions such as PR4-like endochitinase, lignin biosynthesis protein, and CC-NBS-LRR and TMV resistance protein encoding genes.
GO terms associated with hormone signaling were also overrepresented in all species at 15 DAI. The ethylene-mediated signaling pathway was overrepresented in P. deltoides

Gene expression related to visible phenotypes and defense
Two major phenotypes were observed in poplar following infection. A more precise examination of gene expression reveals that multiple lignin biosynthesis genes of P. deltoides are significantly upregulated after leaf spot formation (Fig 5). Among these genes were a single caffeoyl-CoA 3-O-methyltransferase (CCoAMT) that was more than 6 log2-fold upregulated, three cinnamoyl CoA reductases (CCR) were more than 4 log2-fold upregulated and three cinnamylalcohol dehydrogenases (CAD) that were more than 4 log2-fold upregulated. The highest overexpressed CAD gene, potri.009g062800.1, shares close similarity with the defense-related A. thaliana ELI3-2 gene [31].
Another observable phenotype was the senescence and abscission of leaves in P. balsamifera following infection with S. populicola (Fig 3). To determine an underlying molecular cause for this phenotype, gene expression data was screened for the presence of annotations associated with leaf abscission and senescence. A small subset of genes putatively involved in leaf senescence were overexpressed in P. balsamifera at 4 DAI (Table 3) prior to visible disease symptoms. By 15 DAI, two additional senescence-related genes were significantly overexpressed: Potri.001G112600 and Potri.010G166200. These were significantly overexpressed in all three species at 15 DAI compared to the control. In this study, only leaf tissue was examined and petioles were excluded; thus, abscission zone tissue-specific genes were not screened.
Upregulation of genes with annotation relating to plant defense were detected at 1 DAI. Among these genes upregulated were chitinases, β-glucosidases and NBS-LRR resistance proteins (Table 4). Chitinase activity was the most common defensive function observed at 1 DAI, representing 8 of 37 shared upregulated defensive genes, though no individual chitinase genes were commonly regulated between pathosystems. Two β-glucosidases were commonly  upregulated at 1 DAI between P. deltoides and P. balsamifera. After 1 DAI, those genes were no longer significantly upregulated in P. deltoides, but were upregulated in P. balsamifera throughout the entire time course. In P. tremuloides, one β-glucosidase (Potri.004G110200) was significantly upregulated at 4 DAI, and it remained upregulated at 15 DAI. The expression of a group of plant proteases belonging to the Kunitz-type were significantly regulated in all pathosystems (Table 5). In P. deltoides, one upregulated and two downregulated Kunitz trypsin inhibitors were detected at 15 DAI. In P. balsamifera, three Kunitz trypsin inhibitors were observed upregulated early during infection, but none were upregulated at 4 DAI. At 15 DAI, 15 genes belonging to this protease family were overexpressed. In P. tremuloides, four or more Kunitz trypsin inhibitors were upregulated at all of the sampled time points. A number of different transcription factors showed differential expression throughout the infection. Of note were WRKY domain containing transcription factors that were differentially expressed in all pathosystems at 15 DAI (Table 4). Overall, more than one gene with significant expression mapped to each of the WRKY75, WRKY72, and WRKY70 families. However, overexpression of WRKY72 genes was not detected in P. deltoides.

Comparison of RNA-seq analysis to RT-QPCR
RT-qPCR was conducted using a small subset of genes for comparison between two time points: 0 and 15 DAI. Of the ten genes examined, seven genes followed the expected pattern of upregulation, downregulation or no change; two genes (Potri.001G449800 and Potri.010G106900) showed expression patterns different than expected in one poplar species, and one gene (Potri.013G041900) showed large inconsistency in observed gene expression in two pathosystems (Fig 6). All three genes that showed inconsistent expression between RNAseq and RT-QPCR have close sequence similarity to other members in their respective gene family in P. trichocarpa. As such, it is possible that transcripts from other gene family members were co-amplified. All primers were designed using the P. trichocarpa reference sequences, but the use of more accurate species specific template libraries for primer design may allow for more precise RT-QPCR results. The lack of correspondence between RT-qPCR and RNAseq analysis may also be due to errors in read mapping related to the use of a non-species specific genome template or other issues generated with this data analysis pipeline.

Discussion
The two unique host responses to leaf spot infection that observed. The first was the differential expression of many lignin biosynthesis in genes in P. deltoides. Lignification of cell walls surrounding sites inoculated with the Dothideomycetes pathogens was reported following Zymoseptoria tritici and Phaeosphaeria nodorum infection in wheat [32] and S. musiva stem canker in hybrid poplar [33,34]. During S. musiva canker formation on hybrid poplar, several layers of ligno-suberin necrophylactic periderm can form between the hyphal colonized tissues [32,34].
The presence or absence of these tissue layers may represent an accurate marker for disease resistance in hybrid poplar [34]. Susceptible hybrid clones were not observed to produce a necrophylactic periderm, while these barriers were observed to form rapidly in resistant trees inoculated without pre-wounding [34]. Other transcriptome studies in poplar have detected genes involved in lignin biosynthesis in response to wounding [35] and rust infection by both M. larici-populina and M. medusa [15]. However, no significant expression of genes directly involved in lignin biosynthesis was reported in either resistant or susceptible hybrid poplars in response to S. musiva leaf infection [19], though the upregulation of laccases encoding genes (diphenol oxidases) was observed. Laccases were also observed to be upregulated in all pathosystems at 15 DAI. Senescence and abscission of diseased leaves in P. balsamifera was a unique phenotype. Genes putatively involved in senescence were overexpressed early in infection prior to the occurrence of any visible leaf symptoms. Leaf senescence and abscission in response to disease caused by Drepanopeziza populi and M. medusae was reported in another poplar species, P. angustifolia, particularly if plants were pre-treated with Penicillium spp. endophytes [8,36]. In blueberry, leaf senescence and abscission were reported following severe leaf infection by the Sphaerulina relative Septoria albopunctata [37]. In these studies, the pathogens caused leaves to senescence before abscission. Leaf abscission could be an indicator of severe infection or a triggered difference response. Within a natural forest setting, senescence and abscission of infected leaves early in the growing season could prove to be beneficial by removing a source of inoculum. Though leaf senescence and abscission was not observed in P. deltoides and P. tremuloides, some senescence-related genes were significantly overexpressed at later stages of infection.
The timing of necrotic spot appearance was fairly consistent between the different pathosystems, ranging from 8 to 14 DAI. This is within the range of necrotic leaf spot formation reported for various plant species infected with other Mycosphaerellaceae family pathogens. Previous studies have indicated the appearance of significant symptoms at approximately 20-32 days in wheat and blueberry following inoculation with Z. tritici and S. albopunctata, respectively [37,38]. Also, leaf spot symptoms on hybrid poplar leaf discs inoculated with S. musiva became visible in as early as 6 days [39]. Different inoculation conditions or differences in disease epidemiology may account for these observed differences.
Major differences in gene expression were observed in host gene expression of three poplar species during leaf spot infection by co-evolved Sphaerulina pathogens. However, RNA-seq data from infected leaves did show some common trends among the large differences in observed phenotypes. One common observation between all species was the expression of genes involved in defensive reactions early during infection, particularly chitinases and β-glucosidases (Table 3). Chitin has a primary role in pathogen detection in plants [40] and the upregulation of chitinase genes was previously reported after inoculation of hybrid poplar with the rust pathogen M. medusae [17]. Chitinases were also reported to be upregulated in hybrid poplar leaves 4 DAI with S. musiva [19]. In both P. deltoides and P. balsamifera, two β-glucosidase were identified as being moderately upregulated, indicating the possible start of cell wall modifications to defend against the invading pathogens. β-glucosidases are known to activate a group of compounds known as 'phytoanticipins', which are inert chemical defenses [41].
Kunitz trypsin inhibitor encoding genes are another group of defense related genes that were upregulated following challenge by Sphaerulina leaf pathogens in their poplar hosts (Table 4). Kunitz protease inhibitors can protect plant cells by blocking the action of pathogen proteases [42]. These genes were shown to be highly upregulated in hybrid poplar 24 h following wounding or insect herbivory [35,43,44]. Leaf infection of hybrid poplar by M. medusae induced a small upregulation of two Kunitz trypsin inhibitors at 1 DAI; however, these genes were downregulated at later stages of infection as the pathogen began to heavily colonize the leaves [17]. The early expression of trypsin inhibitors may indicate that cell damage occurred early in infection, though this was not quantified in this study.
Differential expression of many transcription factors was detected. WRKY family transcription factors were expressed by all three poplar species. A gene (Potri.015g099200) with similarity to Arabidopsis WRKY75 was significantly upregulated at 15 DAI in all three host species (Table 4). WRKY75 is involved in defense against both bacterial and fungal pathogens [45], and in Arabidopsis, loss of function for this gene showed delayed leaf senescence [46]. Overexpression of a single WRKY70 gene was also detected in all pathosystems (Potri.013G090300). In Arabidopsis, WKRY70 was identified as a key gene used in modulation of plant defense responses through suppression of jasmonic acid signalling [47]. These genes may play an important role in resistance to Sphaerulina pathogens and as the inducers of leaf abscission in P. balsamifera.
The functional annotations of differentially expressed genes revealed different hormone signalling pathways being activated in the different pathosystems. Ethylene signalling was observed in P. deltoides, and jasmonic acid and salicylic acid signalling were observed in P. balsamifera. Multiple hormone pathways were activated in P. tremuloides. Changes in hormone signalling gene expression were typically observed at 15 DAI with only ethylene signalling being detected in P. tremuloides during early infection. Some possible links can be hypothesized between hormone signalling and the observed plant phenotypes. Ethylene is known to play an important role in the inducible production of lignin [48], while jasmonic acid signalling is known to induce terpenoid production in plants [49]. The activation of these signalling pathways correlates with the observed upregulation of lignin and terpenoid biosynthesis genes in the different pathosystems. Different hormones have been implicated in defense against different pathogen classes. Ethylene and jasmonic acid signalling are typically thought to be involved in the defense against necrotrophic pathogens, while salicylic acid is an essential signalling hormone for defense against biotrophic and hemi-biotrophic pathogens [50,51]. The activation of different hormone signalling pathways may be an indicator of the hemi-biotrophic Sphaerulina pathogens at different stages of their infection cycle.
The data presented here shows an in depth examination of gene expression in different poplars with leaf spots induced by the different Sphaerulina species. Intraspecific variation in disease resistance was not examined; instead, what is presented here is a snapshot of the possible gene expression during Sphaerulina leaf spot infection in poplar. A number of commonly expressed genes were identified between the species examined among some large differences observed in whole pathway analysis and the associated plant phenotypes. The libraries generated here will be valuable for future studies. One limitation of this study was the use of the P. trichocarpa genome, as the poplar species we studied have different degrees of evolutionary distance from the model tree. A re-examination of this data may be conducted in the future when the precise genomes of these poplar species are completed and publically released and may also be utilized to annotate genes in future genome studies. An in-depth analysis of pathogen gene expression is also currently underway and will be published separately.  Table. List of primers used in this study for validation of RNAseq data and measurement of fungal growth by relative quantification of plant and fungal genes. Table. Individual statistics for read mapping to P. trichocarpa genome.