Many conifers have distributions that span wide ranges in both biotic and abiotic conditions, but the basis of response to biotic stress has received much less attention than response to abiotic stress. In this study, we investigated the gene expression response of lodgepole pine (Pinus contorta) to attack by the fungal pathogen Dothistroma septosporum, which causes Dothistroma needle blight (DNB), a disease that has caused severe climate-related outbreaks in northwestern British Columbia. We inoculated tolerant and susceptible pines with two D. septosporum isolates and analyzed the differentially expressed genes, differential exon usage, and co-expressed gene modules using RNA-seq data. We found a rapid and strong transcriptomic response in tolerant lodgepole pine samples inoculated with one D. septosporum isolate, and a late and weak response in susceptible samples inoculated with another isolate. We mapped 43 of the DEG- or gene-module-identified genes to the reference plant-pathogen interaction pathway deposited in KEGG database. These genes are present in PAMP-triggered and effector-triggered immunity pathways, including genes encoding mitogen-activated protein kinase and disease resistance protein. Genes comprising pathways and gene modules had signatures of strong selective constraint, while the highly expressed genes in tolerant samples appear to have been favored by selection to counterattack the pathogen. We identified candidate resistance genes that may respond to D. septosporum effectors. Taken together, our results show that gene expression response to D. septosporum infection in lodgepole pine varies both among tree genotypes and pathogen strains, and involves both known candidate genes and a number of genes with previously unknown functions.

The tolerant and susceptible lodgepole pine (Pinus contorta) samples used in this study were created by crossing parent trees that were ranked for foliage retention in the presence of Dothistroma septosporum. The progeny test site was planted in the year 2000 near Kispiox, British Columbia, Canada. The test included 121 families (117 wind-pollinated families and 4 full-sib families) and eight wild-stand seedlots. The test was established using a randomized complete block design with eight blocks, and families were planted in four-tree row plots within each block. Infection by D. septosporum was confirmed by a field pathologist prior to assessment. Each tree was assessed at age 12 for the percent of foliage retention on a scale of 0 (complete defoliation) to 100% (complete foliage retention). At this site, trees at age 12 should have foliage at ground level (100% foliage retention), because crown closure hasn’t yet occurred. D. septosporum infection in the lower crown of trees will cause foliage loss, thus assessing the percent of foliage retention is a good indication of the ability of a tree to tolerate D. septosporum. The mean foliage retention of the 121 families on site was 48% (ranged from 34% to 66%) and the mean of the eight wild seedlots was 49%. A model was run to test the difference between family means using the R package “lmer” (Bates et al., 2015). This model incorporated families as fixed effects, blocks and block × family as random effects, and foliage retention at age 12 of all trees in the 125 open pollinated families and seedlots (117 wind-pollinated families and eight wild-stand seedlots) as the dependent variable. A multiple comparison test was performed using the Tukey’s adjustment implemented by the R package “emmeans” (Lenth, 2020). We found a significant difference (P-value < 0.00) between two families with high foliage retention and two families with low foliage retention. Two parents of families (exposed to D. septosporum in a field trial) with high foliage retention were crossed and two with low foliage retention were crossed to create the tolerant and susceptible families used for this study. The offspring of the two tolerant parent trees had a mean foliage retention of 60%, and the offspring of the two susceptible parent trees had a mean foliage retention of 41%. The controlled crosses were made in Orchard 230 at the Kalamalka Seed Orchard in Vernon, British Columbia. Seeds from crossed families were harvested and used to produce seedlings that were inoculated in our experiment. Seedlings from tolerant parents were named as T, and those from susceptible parents as S.

T and S seedlings, between one- and two-year-old, were inoculated with two D. septosporum isolates (D1 & D2) separately with 4 – 6 replicates. D1 was isolated from infected hybrid P. contorta × P. banksiana in Alberta, Canada; while D2 was isolated from infected P. contorta in Smithers, British Columbia, Canada. D. septosporum conidia were harvested from colonies grown for 10 – 16 days on Dothistroma sporulation medium plates (Bradshaw et al., 2000) by suspension in sterile distilled water. Each seedling was inoculated with a standard inoculum of approximately 3 mL of 1.6 × 10⁶ conidia mL^–1, using a household trigger action atomizer. The plants were left on the inoculation bench until the needles were dry (~ 45 minutes), then wrapped in transparent plastic bags and placed in a growth chamber (16 hours daylight at 20°C; 8 hours night at 12°C) for 48 h. A mister spraying tap water was activated every hour for 3 minutes for the duration of the trial to maintain needle wetness. Control seedlings that were inoculated with water were also prepared. Seedlings were rated according to the level of disease observed on a linear scale from 0 to 100% of symptomatic needles, i.e. necrotic needles with red bands and/or fruiting bodies, at the end of the experiment (12 weeks post inoculation). The numbers of samples rated for tolerant seedlings are: 5 — Control, 7 — inoculated with D1, 7 — inoculated with D2; the numbers of samples rated for susceptible seedlings are: 7 — Control, 6 — inoculated with D1, 6 — inoculated with D2. Samples for RNA sequencing (RNA‐seq) were taken at 4 days, 8, and 12 weeks post inoculation when fungal surface growth, new immature lesions and sporulating lesions were expected and/or observed, respectively (Table S7). 8 – 12 necrotic needles with disease lesions were preferentially sampled. Needles were flash frozen in liquid nitrogen, and stored at -80°C.

Total RNA was extracted from frozen needles. Needles were manually ground in liquid nitrogen using a mortar and pestle, and approx. 60mg of powdered tissue was used for the RNA extraction. Total RNA was extracted from each of 90 samples using Spectrum Plant Total RNA Kit (Sigma Aldrich, St. Louis, MO, USA) in a random order with respect to treatment to avoid batch effects. After eliminating traces of genomic DNA by treatment with On-Column DNase I Digest Set (Sigma, Aldrich, St. Louis, MO, USA), RNA quality was evaluated using an automated capillary gel electrophoresis and quantity was assessed using a Bioanalyzer 2100 with RNA 6000 Nano Labchips (Agilent Technologies, Santa Clara, CA, USA). A260/A280 and A260/A230 ratios were measured using a Nanodrop fluorometer (Nanodrop Products, Wilmington, DE, USA). Total RNA concentration was determined by Qubit 2 using Qubit RNA HS Assay Kit (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). For all the samples,  the A260/A280 ratios were 1.94 or higher, A260/A230 ratios were 1.2 or higher, with the RNA integrity number (RIN) ranging from 6.8 to 8.0. Standard paired end 100bp Illumina RNA-seq libraries were prepared using NEB mRNA stranded library preparation kits (New England Biolabs, Ipswich, MA, USA) at Centre d'expertise et de services Génome Québec and sequenced on 5 lanes of an Illumina HiSeq4000 instrument.

Raw reads were processed with FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/), including clipping adaptors, filtering artifacts, and keeping reads with a minimum quality score of 20. De novo transcriptome assembly was performed using Trinity v2.5.1 (Grabherr et al., 2011) with an input file gathering reads from the library with the highest number of reads for each of the 18 treatments.  The de novo assembled transcripts were then aligned to a previous lodgepole pine transcriptome (Yeaman et al., 2014) using blastn with an E-value threshold of 1e-10. Newly assembled transcripts that could be aligned to the previous transcriptome were removed, and the remaining transcripts were added to the previously assembled transcriptome. This combined transcriptome was further examined to filter out low-expression and contaminant genes. Filtered RNA-seq reads (from each of the 90 libraries) were mapped to the combined transcriptome using the software RSEM with default parameters (Li and Dewey, 2011). Transcripts with high expression levels were identified from RSEM mapping output as those having a transcript per million (TPM) value ≥ 1.00 in all replicates of any sample type (T or S or control) at any harvest time point per D. septosporum isolate.

To remove any eventual contamination, the combined transcriptome was aligned against the non-redundant (nr) protein database (downloaded on March 24, 2018) using the blastx function and “–more-sensitive” mode implemented by the software DIAMOND (Buchfink et al., 2015). Query transcripts with high similarity (E-value threshold of 1e-10) to the green plant subjects were then selected. Additionally, the combined transcriptome was aligned against the loblolly pine reference genome (Neale et al., 2014) using the software GMAP-2017-06-20 (Wu and Watanabe, 2005) to identify transcripts with high similarity (identity ≥ 75% and coverage ≥ 75%) to the loblolly pine genes. Transcripts with high similarity to green plants or loblolly pine were combined into one single file. To remove redundant genes, the longest isoform with a length ≥ 300bp for each gene was selected, and transcripts with high expression levels (TPM value ≥ 1.00) were retained. Putative D. septosporum genes were removed from the combined transcriptome by searching for transcripts assigned to “Dothistroma” or “Dothideomycetidae” in a blastx search against the nr database, and by using GMAP-2017-06-20 to identify transcripts with a high similarity (identity ≥ 50% and coverage ≥ 50%) to D. septosporum genes in the D. septosporum reference genome (de Wit et al., 2012). At the end, we obtained a lodgepole pine reference transcriptome that includes 27,371 transcripts.

To annotate the lodgepole pine reference transcriptome with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms, the transcripts were translated to protein sequences using TransDecoder (https://github.com/TransDecoder/TransDecoder/wiki), and then were annotated using eggNOG4.5 (Huerta-Cepas et al., 2016), InterProScan 5 (Jones et al., 2014), and Blast2Go (Götz et al., 2008). All annotations were combined. To make a uniform KEGG annotation format, EC numbers were converted to KOs according to the hierarchical text file, which was downloaded from https://www.genome.jp/kegg-bin/get_htext?br08902. 

The file Pinus_contorta_reference_transcriptome includes 27,371 lodgepole pine (Pinus contorta) genes (in FASTA format) that were analyzed in the manuscript.

The files Pcontorta_annotation_GO, Pcontorta_annotation_KO, and Pcontorta_annotation_description contain the annotated GO (Gene Ontology) terms, KO (KEGG Orthology ) terms, and potential gene function separately. The first column contains the gene identifiers which can be found in the file Pinus_contorta_reference_transcriptome. The second column contains the annotation terms.

The files AM_list, DEG_list, DEU_list, pathogenpathway_gene, includes genes identfied from different analyses. The AM_list file includes the genes of the gene modules that were related to D. septosporum infection phenotypes; DEG_list includes gene identified as differentially expressed genes; DEU_list includes the genes containing differential exon usage; pathogenpathway_gene includes the genes that can be mapped to plant-pathogen interaction pathway. 

T — tolerant lodgepole pine samples; S — susceptible lodgepole pine samples; D1 — D. septosporum isolate 1; D2 — D. septosporum isolate 2; 4d — 4 days post inoculation; 8w — 8 weeks post inoculation; 12w — 12 weeks post inoculation.