Constitutive and insect‐induced transcriptomes of weevil‐resistant and susceptible Sitka spruce

Abstract Spruce weevil (Pissodes strobi) is a significant pest of regenerating spruce (Picea) and pine (Pinus) forests in North America. Weevil larvae feed in the bark, phloem, cambium, and outer xylem of apical shoots, causing stunted growth or mortality of young trees. We identified and characterized constitutive and weevil‐induced patterns of Sitka spruce (Picea sitchensis) transcriptomes in weevil‐resistant (R) and susceptible (S) trees using RNA sequencing (RNA‐seq) and differential expression (DE) analyses. We developed a statistical model for the analysis of RNA‐seq data from treatment experiments with a 2 × 3 factorial design to differentiate insect‐induced responses from the effects of mechanical damage. Across the different comparisons, we identified two major transcriptome contrasts: A large set of genes that was constitutively DE between R and S trees, and another set of genes that was DE in weevil‐induced S‐trees. The constitutive transcriptome unique to R trees appeared to be attuned to defense, while the constitutive transcriptome unique to S trees was enriched for growth‐related transcripts. Notably, a set of transcripts annotated as “fungal” was detected consistently in the transcriptomes. Fungal transcripts were identified as DE in the comparison of R and S trees and in the weevil‐affected DE transcriptome of S trees, suggesting a potential microbiome role in this conifer‐insect interaction.

on the bark, phloem, cambium, and outer xylem of the PYAS, which disrupts the flow of water and nutrients and destroys apical growth.
The endophase is the most destructive life stage of the weevil until it reaches pupation and eventually emergence of young adults. Young trees between 2 and 20 years of age are particularly vulnerable to the weevil, but this range can substantially vary depending on factors such as the level of weevil pressure and rate of tree growth.
The introgression of weevil resistant (R) spruce genotypes into forest landscapes is currently the primary means for tree breeding programs to mitigate weevil damage (King & Alfaro, 2009;King et al., 2011;Kiss & Yanchuk, 1991;. In the Pacific Northwest, Sitka spruce (P. sitchensis) is particularly susceptible to this insect (King & Alfaro, 2009). However, weevil-resistant Sitka spruce genotypes have been discovered. Durable resistance involves a synergistic defense syndrome of at least two defense traits: abundance of cortical stone cells and terpene-rich oleoresin (Krokene, 2016;Whitehill, Henderson, Schuetz, et al., 2016;Whitehill, Henderson, Strong, et al., 2016;. While these two traits appear to provide much of the effective defense against weevil attack, additional factors such as early induced responses may contribute to the outcome of spruce-weevil interactions . As part of unveiling the genomic portrait of spruce-weevil interactions, we have used RNA-seq to capture the host transcriptome component. Analysis of RNA-seq data in plant-insect interactions often focuses on differentially expressed (DE) genes and may apply complex factorial experimental designs that target unique aspects of the ecological interaction. For spruce-insect systems, various treatments have been used to differentiate host responses to artificial wounding from molecular signatures affected by real insect attack (Byun-McKay et al., 2003Franceschi et al., 2002;Mageroy, Christiansen, et al., 2020;Mageroy, Wilkinson, et al., 2020;Martin et al., 2002;Ralph et al., 2006Ralph et al., , 2008Schiebe et al., 2012;. Genes that respond to mechanical wounding can be removed from DE analyses with the goal to identify insect-specific responses. Statistical models used for transcriptome analyses often rely on open-source software packages in R. For example, R-based approaches were used for the DE analysis in RNA-seq datasets to explore the molecular underpinnings of methyl jasmonate-induced resistance against bark beetles in Norway spruce (P. abies; Mageroy, Wilkinson, et al., 2020).
Here, we expanded the analyses of a previously described RNAseq resource for Sitka spruce-weevil interactions . Specifically, in this study, we provide a detailed analysis of insect-specific transcriptome responses and include microbial signatures in our overall analysis. The results from this analysis highlighted two biologically relevant contrasts: genes that are DE between resistant (R) and susceptible (S) trees and genes that are DE between S trees challenged by weevil larvae relative to mechanical wounding or untreated controls. The DE genes identified in the first contrast were defined as the "constitutive difference" (CD); the DE genes from the second contrast were defined as "weevil-induced" (WI). In previous work, signatures of fungi in spruce RNA-seq datasets were removed from DE analyses (Delhomme et al., 2014(Delhomme et al., , 2015. This is despite fungal transcripts being consistently detected in replicated RNA-seq libraries . In the transcriptome analysis described here, we purposely kept genes annotated as non-plant, as they may provide clues towards a role of a microbiome associated with the host tree or with the insect in the spruce-weevil interaction.

| Plant materials
Origins and growing conditions of clonally propagated R and S Sitka spruce genotypes were described previously Whitehill, Henderson, Schuetz, et al., 2016;Whitehill, Henderson, Strong, et al., 2016;. In brief, the R genotype H898 derives from a highly weevil resistant population of Sitka spruce located in the low elevation mainland of British

| Experimental design
The experimental design was previously described in . Here, we present the experimental procedures for downstream analyses used in the present study ( Figure S1). Twelve R and twelve S trees were used for RNA-seq experiments with four trees (biological replicates) used for each genotype × treatment combination. Treatments started on May 20th, 2013 (465 GDD) and concluded July 2nd, 2013 (923 GDD). Trees were arranged in a completely randomized design. One of three treatments were applied within a section of 0-2 cm from the tip of the PYAS to mimic the site of natural oviposition by weevils in the field (Whitehill, Henderson, Schuetz, et al., 2016). The average PYAS length where treatments were applied was 26.6 ± 2.2 cm (R) and 34.4 ± 2.0 cm (S). Needles within this section were removed with a razor blade to facilitate treatment applications. The three treatments were Columbia, Canada (50°16′00″N; 119°16′18″W). Egg isolations (Whitehill, Henderson, Strong, et al., 2016) and collection conditions  were described previously. The top 2 cm of the PYASs were wrapped with parafilm to protect treatment sites. Trees were kept under natural conditions outside at the University of British Columbia Vancouver Point Grey campus for the duration of the experiment. PYAS samples were collected as illustrated in Figure S1 at the time when weevil larvae entered the final stage of development in S trees. The experiment was performed over a 50-day period .

| Collection of bark samples
PYASs were harvested and bark samples, which consisted of outer bark, cortex, phloem, cambium, and traces of outer xylem, isolated on July 2nd at 923 growing degree-days (GDD). Details regarding developmental stage of weevil larvae were previously described . Briefly, R tree larvae were in the 3rd and S tree larvae were in the 4th instar of development at the time of sample collection. Samples for control treatments (1) and AOC treatments (2) were harvested at the top 0-2 cm of the PYAS; samples for gallery treatments (3) were collected by excising the bark surrounding the final weevil-feeding site as described in Whitehill et al. (2019; Figure S1). Weevil gallery tissues used for RNA-seq included the lower 2 cm of a larvae gallery with a 1 cm border of adjacent bark tissue including the larval feeding margin after removal of the insect and fecal materials. Gallery tissues from a given individual tree were pooled and treated as a single biological replicate. All samples were flash frozen in liquid N 2 and stored at −80℃ until RNA extraction.

| Transcriptome resource and DE analysis
The complete methods for RNA extraction, Illumina HiSeq library preparation, sequencing on the HiSeq2000 platform, raw data processing, quality assessments, and de novo transcriptome assembly were previously described ; except, for the present analyses, contigs with annotations to fungal, insect, and other (bacterial, virus, etc.) origins were not removed from the final assembly. Raw reads are available at NCBI SRA (BioProject PRJNA398042). DE analyses were performed using the voom/ limma package (Law et al., 2014) in R on coding sequences (CDS), predicted by TransDecoder (Haas et al., 2013), with expression of at least 1 count per million (cpm) in at least two libraries and quantified using Sailfish (version 0.10.0) with default parameters (Patro et al., 2014). Statistically significant (adjusted p ≤ .05) CDS with a | log 2 fold change | ≥2 were annotated with BLASTP against the UniProt database with an e-value cutoff at 1e −20 (Tables S1 and   S2). Gene ontology annotations were extracted from corresponding best hit UniProt ID. Taxonomic lineage information of the best hit was extracted from the NCBI Taxonomy Database. Complete GO term annotations and taxonomic lineage are available in Table   S2.

| RNA-seq gene expression model for a 2 × 3 factorial design
To remove the effect of wounding on contigs identified as DE in gallery treatments (3), analyses were modeled as a 2 × 3 factorial design for two genotypes and three treatment combinations. A full analysis using the 2 × 3 factorial model revealed only two comparisons had large differences in gene expression, one being the comparison between control treatments (1) of resistant (R) and susceptible (S) trees defined as constitutive differences (CD) and the other being DE genes between AOC (2) and gallery (3) treatments defined as weevil-induced (WI) in S trees. A complete description of the model and analysis is available at https://github.com/myuen/ White_Pine_ Weevil_DE. DE transcripts were identified as WI using the interactive model by removing the AOC wounding effect (2) from the compound interactive effect of the gallery treatment (3) and are assumed to represent genes affected specifically by weevil activity.

| Comparisons of DE transcriptome contrasts
Comparisons of statistically significant (adjusted p ≤ .05 and |log2 fold change| ≥2) DE contigs from the two main contrasting comparisons were performed. Genes DE in the CD contrast were compared against the WI contrast. Genes were grouped by expression patterns into eight different sets (Table S2). Genes with higher constitutive expression exclusively in CD R trees were labeled Set 1, while genes with lower expression exclusively in CD S trees were labeled Set 2. WI genes exclusively up-regulated or down-regulated in S trees were labeled Sets 3 and 4, respectively. Four additional contrasts were found at the intersection of the two-major CD versus WI contrasts. These were defined as higher CD expression in R trees and WI up-regulated in S trees (Set 5), higher CD expression in S trees and WI up-regulated in S trees (Set 6), higher CD expression in R trees and WI down-regulated in S trees (Set 7), and higher CD expression in S trees and lower WI expression in S trees (Set 8). Heatmaps were generated for DE contigs in Sets 5-8 using the D3HEATMAP package (Cheng, 2016)

| Overall patterns of DE contigs
RNA-seq data were generated from 24 cDNA libraries comprising PYAS bark tissue from R and S trees under three different treatment conditions, control (1), AOC (2), and gallery (3). The overall patterns of gene expression between genotype and treatment groups are shown in Figure S2. In the control transcriptomes, 4468 genes were DE between R and S trees. Of these genes, 2254 had higher transcript levels in the R genotype and 2214 had higher transcript levels in the S genotype (Table 1). The genes of this contrast were defined as "constitutive difference" (CD) between R and S trees and annotated using UniProt terms (3731 genes- Table S2) or GO terms (942 genes- Table S2). In the "weevil-induced" (WI) transcriptomes, It is important to note that the experiment was designed to capture transcriptome signatures in prolonged treatments of 50 days that resemble the mature stage of a natural weevil infestation. This design may not have captured the early and transiently induced response in R and S trees. Early induced responses have been the topic of previous work (Lippert et al., 2007;Ralph et al., 2006). Effects of AOC alone were removed using a 2 × 3 factorial model to identify DE responses likely to be specific to weevil activity. WI genes were annotated using UniProt terms (3710 genes) or GO terms (1407 genes; Table S2).

| Comparison of constitutive different (CD) and weevil induced (WI) transcriptomes
We compared the two contrasts to identify unique and shared genes between the CD (Figure 1; x-axis) and WI (y-axis) transcriptomes.
We found eight sets of genes that differed by expression patterns, including 3579 genes unique to the CD contrast (Sets 1 and 2), 4830 genes unique to the WI contrast (Sets 3 and 4), and 889 genes shared between the CD and WI contrasts (Sets 5-8). To further explore potential biological functions of DE genes in Sets 1 and 2, Gene Ontology (GO) biological process terms were evaluated. A majority of genes fell into six GO term categories and accounted for 38.2% and 42.7% of all annotated genes with assigned GO terms for Set 1 (R) and Set 2 (S), respectively. The six GO term categories include signal transduction (21.9% Set 1, 26.7%
For CD genes with higher expression in S trees (Set 2), nine of the "defense" terms and thirteen "growth" terms were identified representing 2.9% and 6.2%, respectively. Fisher's exact test was used to evaluate statistically significant enriched GO terms within Set 1 and Set 2. Two terms were enriched in Set 1 and another two terms in Set 2. The growth term, photosynthesis light reaction, was significantly enriched in the S tree dataset ( The superscript letters indicate which of the numbers are higher expressed or upregulated and correspond to the text in the same row. suggested that the constitutive transcriptome signature of R trees is more attuned to "defense", while the constitutive transcriptome signature of S trees is more enriched towards 'growth'.  (Table S2). The most abundant GO term annotation in Set 3 was "carbohydrate metabolic process" (14.6%) while "biosynthetic process" (5.1%) was the most abundant term in Set 4. WI genes up-regulated in Set 3 consisted of eight "defense" and two "growth" GO terms representing 10% and 0.6%, respectively, of all annotated genes. WI genes down-regulated in Set 4 contained 15 "defense" and nine "growth" terms representing 4% and 2%, respectively. A Fisher's exact test revealed 16 (Set 3) and 21 (Set 4) GO term annotations were significantly enriched.

| Genes unique to the WI contrast
Notably, the term "chitin catabolic process" was enriched in Set 3, consistent with observations of fungal transcripts in this set. The GO terms "aromatic compound biosynthetic process" and "oxylipin biosynthetic process", which are relevant to plant defense, were significantly enriched in Set 3. Conversely, GO terms enriched and down-regulated in Set 4 include processes related to cellular damage (i.e., DNA repair, negative regulation of cellular process) and growth (i.e., sexual reproduction, auxin-activated signaling pathway; Table 2).

| Genes shared between CD and WI contrasts
Genes of Sets 5-8 were shared between the CD and WI transcriptomes, with four different expression patterns (Figure 1): Set 5: CD higher expressed in R trees and WI up-regulated in S trees.
Set 6: CD higher expressed in S trees and WI up-regulated in S trees.
Set 7: CD higher expressed in R trees and WI down-regulated in S trees.
Set 8: CD higher expressed in S trees and WI down-regulated in S trees.
Genes of Sets 5, 6 and 7 may contribute to weevil defense or resistance. Volcano plots show expression patterns of DE transcripts Sets 5-8 (Figure 4c-f). The majority of genes were annotated as viridiplantae ( Figure 4a). However, 39% of annotated genes in Set 5 were identified as fungal, while fungal genes were not detected in Sets 6, 7 and 8. This is noteworthy because Set 5 comprises genes that were constitutively more highly expressed in R trees and were up-regulated in S trees in response to long-term weevil activity and may be of particular interest for defense.
A heatmap of the normalized GO biological process terms for Sets 5-8 revealed distinct differences (Figure 4b). In Set 5, the most abundant GO term was "carbohydrate metabolic process" (15% of annotated transcripts). Only 3 GO annotations were available for Set 6 and include "signal transduction", "response to oxidative stress", and "hydrogen peroxide catabolic process". In Set 7 and Set 8, genes annotated as "signal transduction" (25% and 23%, respectively, of annotated transcripts) were the most abundant GO term. ( Figure 5d).

| D ISCUSS I ON
Tree breeding for insect resistance stands to benefit from the application of genomics-based tools, such as genetic markers identified from genome and transcriptome analyses or genomic selection to overcome challenges associated with screening for pest resistance (Beaulieu et al., 2020;Isik, 2014;Lenz et al., 2019;Neale & Kremer, 2011;Plomion et al., 2011). Previous comparisons of R and S trees elucidated defense mechanisms that contribute to weevil resistance, such as oleoresin terpenes and stone cells (Hall et al., 2011;Roach et al., 2014;Whitehill, Henderson, Schuetz, et al., 2016;Whitehill, Henderson, Strong, et al., 2016;. The present transcriptome analysis was performed using a statistical model to identify weevil specific responses and assess differences in gene expression between R  (Miller et al., 2005;Ralph et al., 2006Ralph et al., , 2008. Features in the rapid early induced response in R trees, such as induced terpenoids and traumatic resin ducts, in addition to differences in constitutive defenses, may contribute to the resistance of R trees.
We used a 2 × 3 factorial design to differentiate insect specific responses. Wounding treatment mimicked the mechanical component of oviposition, although it would not mimic effects of biological elicitors associated with eggs or oviposition fluids (Hilker et al., 2005). This treatment served as a control for mechanical effects of oviposition (Gara & Wood, 1989;. The statistical model was parameterized to the experimental design to remove the effect of mechanical wounding from weevil-specific responses. We identified DE genes as CD between R and S trees or as WI in S trees. One set of genes of particular interest (Set 1) are those that were constitutively expressed at higher levels in R trees. These genes may contribute to effective constitutive defense barriers in R trees against feeding or ovipositing adult weevil or against eggs and developing larvae. Conversely, the lack of or lower expression of such genes may contribute to absence of effective constitutive defense in S trees, which in their geographic location of origin are not co-evolved with weevils. Set 1 genes were enriched for GO F I G U R E 4 Intersection of CD and WI genes. (a) Taxonomic annotation of Set 5-8 genes. The distribution of taxonomic groups was normalized to the total number of DE genes within a subset and grouped into three categories that include viridaeplantae, fungi, and other. (b) Hierarchical clustering analysis (heatmap) of Set 5-8 genes with assigned GO term annotations. GO term categories are normalized to the total number of genes that received a GO annotation. Dark blue indicates a majority proportion of GO term representatives are found within a particular set while red indicates no representatives were present. Intermediate lighter blue and red colors indicate at least 1 representative is present. term annotations associated with 'defense' compared to 'growth'.
In contrast, genes constitutively higher expressed in S trees (Set 2) were enriched for 'growth' GO terms. According to the growthdifferentiation balance hypothesis, well-defended plants, such as the Sitka spruce R genotype, invest more in constitutive defenses (Herms & Mattson, 1992).
Set 5 comprised another interesting group of genes, which was constitutively expressed at higher levels in R trees and up-regulated in response to weevil activity in S trees. Almost 40% of the genes in this set are annotated as fungal origin. Reproducible fungal transcriptome signatures may hint towards a role of tree-or insect-associated microbiomes in the spruce-weevil interaction. Fungal endophytes can contribute to a defensive mutualism against herbivores (Saikkonen et al., 2010). Toxic metabolites produced by endophytic fungi have been implicated in providing protection from spruce budworm in red spruce (P. rubens; Clark et al., 1989). While fungal transcripts are commonly observed in conifer transcriptomes, associations with spruce weevil interactions have not been previously documented. Weevils are closely related to bark beetles (Shin et al., 2018), which commonly harbor symbiotic fungi (Bracewell & Six, 2015). It is possible that weevil resistance of spruce not only relies on the induced transcriptome response of the plant, but also on gene expression patterns of hostassociated fungi.

ACK N OWLED G EM ENTS
We thank Dr. Carol Ritland for project management and Ms.
Angela Chiang for help with materials and laboratory management.
The work was supported with funds to JB from Genome Canada, Genome British Columbia and Genome Quebec in support of the Spruce-Up LSARP project, and with funds to JB from the Natural

Sciences and Engineering Research Council (NSERC) of Canada. JB
is a UBC Distinguished University Scholar.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.