DOI: 10.5061/dryad.1vhhmgqps
title: Comparative gene expression analysis reveals mechanism of Pinus contorta response to the fungal pathogen Dothistroma septosporum

description: 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.
description: The tolerant and susceptible lodgepole pine (<i>Pinus
 contorta</i>) samples used in this study were created by crossing
 parent trees that were ranked for foliage retention in the presence of
 <i>Dothistroma septosporum</i>. 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 <i>D.
 septosporum</i> 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.
 <i>D. septosporum</i> 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 <i>D.
 septosporum</i>. 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
 (<i>P</i>-value &lt; 0.00) between two families with high
 foliage retention and two families with low foliage retention. Two parents
 of families (exposed to <i>D. septosporum </i>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 <i>D. septosporum</i> isolates (D1
 & D2) separately with 4 – 6 replicates. D1 was isolated from
 infected hybrid <i>P. contorta</i> × <i>P.
 banksiana</i> in Alberta, Canada; while D2 was isolated from
 infected <i>P. contorta</i> in Smithers, British Columbia,
 Canada. <i>D. septosporum</i> conidia were harvested from
 colonies grown for 10 – 16 days on <i>Dothistroma</i>
 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<sup>6</sup> conidia
 mL<sup>–1</sup>, 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. <i>De novo</i> 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 <i>de novo</i> 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 <i>D. septosporum</i> 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 <i>D. septosporum</i> genes were removed from the combined transcriptome by searching for transcripts assigned to “<i>Dothistroma</i>” 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 <i>D. septosporum</i> genes in the <i>D. septosporum</i> 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. 
description: The file Pinus_contorta_reference_transcriptome includes
 27,371 lodgepole pine (<em>Pinus contorta</em>) 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 <i>D. septosporum </i>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
 — <i>D. septosporum </i>isolate 1; D2 — <i>D.
 septosporum </i>isolate 2; 4d — 4 days post inoculation; 8w — 8
 weeks post inoculation; 12w — 12 weeks post inoculation.