Science, Faculty of
DSpace
UBCV
Lo, Theodora
Master of Science - MSc
University of British Columbia
Abiotic and biotic stresses associated with climate change have been identified as a dominant cause of forest tree mortality in boreal forests. Some tree populations may have the capacity for rapid adaptation or migration to keep pace with changing environmental conditions. One species of interest is black spruce (Picea mariana [Mill.] B.S.P.) as adaptive variation in relation to climate change has previously been reported for this transcontinental North American conifer. As exhibited in studies of other economically important forest trees, genomic resources play a critical role in advancing our understanding of the genomic basis of adaptive variation; however, such resources are lacking for black spruce, with the few available being predominantly transcriptome-related. This thesis describes the first genome assembly of P. mariana with a reconstructed genome size of 18.3 Gbp and NG50 scaffold length of 36.0 kbp. A total of 66,332 protein-coding sequences were predicted in silico and annotated based on sequence homology. To showcase the value of this new genomic resource, phylogenetic and comparative genomics analyses were performed. Phylogenetic trees were estimated from the nuclear and organelle genome sequences of P. mariana and five other spruces. The phylogenetic tree estimated from mitochondrial genome sequences agrees with biogeography; specifically, P. mariana was strongly supported as a sister lineage to P. glauca and three other taxa found in western North America, followed by the European P. abies. In contrast, mixed topologies with weaker statistical support were obtained in phylogenetic trees estimated from nuclear and chloroplast genome sequences, indicative of ancient reticulate evolution affecting these two genomes. Clustering of protein-coding sequences from the six Picea taxa and two Pinus species resulted in 34,776 orthogroups, 560 of which appeared to be specific to P. mariana. Analysis of these specific orthogroups and dN/dS analysis of positive selection signatures for 497 single-copy orthogroups indicate gene functions mostly related to plant development and stress response. The P. mariana genome assembly and annotation provides a valuable resource for forest genetics research and applications in this broadly distributed species, especially in relation to climate adaptation.
THE GENOME OF BLACK SPRUCE: GENOME ANNOTATION & ANALYSES by Theodora Lo B.Sc., The University of British Columbia, 2020 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Bioinformatics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2023 © Theodora Lo, 2023 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled: The Genome of Black Spruce: Genome Annotation & Analyses submitted by Theodora Lo in partial fulfilment of the requirements for the degree of Master of Science in Bioinformatics Examining Committee: Dr. Inanç Birol, Professor, Department of Medical Genetics, UBC Supervisor Dr. Joerg Bohlmann, Professor, Department of Botany, UBC Supervisory Committee Member Dr. Steven Hallam, Professor, Department of Microbiology & Immunology, UBC Supervisory Committee Member Dr. Marco Todesco, Assistant Professor, Department of Botany, UBC Additional Examiner iii Abstract Abiotic and biotic stresses associated with climate change have been identified as a dominant cause of forest tree mortality in boreal forests. Some tree populations may have the capacity for rapid adaptation or migration to keep pace with changing environmental conditions. One species of interest is black spruce (Picea mariana [Mill.] B.S.P.) as adaptive variation in relation to climate change has previously been reported for this transcontinental North American conifer. As exhibited in studies of other economically important forest trees, genomic resources play a critical role in advancing our understanding of the genomic basis of adaptive variation; however, such resources are lacking for black spruce, with the few available being predominantly transcriptome-related. This thesis describes the first genome assembly of P. mariana with a reconstructed genome size of 18.3 Gbp and NG50 scaffold length of 36.0 kbp. A total of 66,332 protein-coding sequences were predicted in silico and annotated based on sequence homology. To showcase the value of this new genomic resource, phylogenetic and comparative genomics analyses were performed. Phylogenetic trees were estimated from the nuclear and organelle genome sequences of P. mariana and five other spruces. The phylogenetic tree estimated from mitochondrial genome sequences agrees with biogeography; specifically, P. mariana was strongly supported as a sister lineage to P. glauca and three other taxa found in western North America, followed by the European P. abies. In contrast, mixed topologies with weaker statistical support were obtained in phylogenetic trees estimated from nuclear and chloroplast genome sequences, indicative of ancient reticulate evolution affecting these two genomes. Clustering of protein-coding sequences from the six Picea taxa and two Pinus species resulted in 34,776 orthogroups, 560 of which appeared to be specific to P. mariana. Analysis of these specific orthogroups and dN/dS analysis of positive selection signatures for 497 single-copy iv orthogroups indicate gene functions mostly related to plant development and stress response. The P. mariana genome assembly and annotation provides a valuable resource for forest genetics research and applications in this broadly distributed species, especially in relation to climate adaptation. v Lay Summary Increased tree mortality due to climate change-induced drought and temperature stress have occurred in many tree species and will likely escalate as conditions worsen. Black spruce (Picea mariana [Mill.] B.S.P.) is of particular interest as it is a widely distributed species in the North American boreal forest and local populations can exhibit genetic variation for adaptive traits. This thesis describes the first genome assembly and annotation of black spruce. I utilized this new genomic resource to investigate the evolutionary relationships between P. mariana and five other spruces and identified potential genes related to stress adaptation. The discordant phylogenies obtained from nuclear and chloroplast genomes compared to the mitochondrial phylogeny strongly suggests interspecific gene flow via hybridization is affecting the former two genomes. Furthermore, comparative genomics analyses highlighted the importance of genes related to plant development and stress response. vi Preface This work was completed under the supervision of Dr. Inanç Birol from the Bioinformatics Technology Lab (BTL) at Canada’s Michael Smith Genome Sciences Centre (GSC). The task of assembling and annotating the black spruce genome was initially conceived by Dr. Ashley Thomson at Lakehead University and Dr. Inanc Birol, while I proposed the downstream analyses in conjunction with Dr. Inanc Birol to further expand this project. A more detailed breakdown of contributions is described below. A version of this material has been published: Lo T, Coombe L, Gagalova KK, Marr A, Warren RL, Kirk H, Pandoh P, Zhao Y, Moore RA, Mungall AJ, Ritland C, Pavy N, Jones SJM, Bohlmann J, Bousquet J, Birol I, Thomson A. Assembly and annotation of the black spruce genome provides insights on spruce phylogeny and evolution of stress response. G3 (Bethesda). 2023 Oct 24:jkad247. doi: 10.1093/g3journal/jkad247. TL annotated the black spruce nuclear genome, performed all analyses, and drafted the manuscript. LC assembled both the black spruce nuclear and mitochondrial genomes and performed repeat-masking on the nuclear genome prior to annotation. KKG trained the AUGUSTUS model and assisted in compiling protein evidence used in the annotation process. JBou and NP shared valuable biological insights on the analyses results and contributed to writing related sections. IB supervised the project. vii Table of Contents Abstract ......................................................................................................................................... iii Lay Summary ................................................................................................................................ v Preface ........................................................................................................................................... vi Table of Contents ........................................................................................................................ vii List of Tables ................................................................................................................................ ix List of Figures ................................................................................................................................ x List of Symbols ............................................................................................................................. xi List of Abbreviations .................................................................................................................. xii Acknowledgements .................................................................................................................... xiii Dedication ................................................................................................................................... xiv Chapter 1: Introduction ............................................................................................................... 1 1.1 Climate change and projected impacts on the boreal forest ........................................... 1 1.2 Genomic approach to understanding climate change tolerance ...................................... 3 1.3 Black spruce and existing genomic/transcriptomic resources ........................................ 4 1.4 Thesis objectives ............................................................................................................. 6 Chapter 2: Methods ...................................................................................................................... 7 2.1 DNA extraction and sequencing ..................................................................................... 7 2.2 Mitochondrial genome assembly .................................................................................... 8 viii 2.3 Nuclear genome assembly .............................................................................................. 9 2.4 Identification and annotation of protein-coding sequences .......................................... 10 2.5 Phylogenetic analyses ................................................................................................... 12 2.6 Comparative genomics analyses ................................................................................... 13 2.6.1 Analysis of P. mariana specific orthogroups ........................................................ 13 2.6.2 Identifying positively selected protein-coding sequences in P. mariana.............. 13 Chapter 3: Results and discussion ............................................................................................. 16 3.1 Genome sequencing, assembly and annotation ............................................................ 16 3.2 Phylogenetic analyses ................................................................................................... 22 3.3 Comparative genomics.................................................................................................. 26 3.3.1 P. mariana specific orthogroup analysis ............................................................... 26 3.3.2 Protein-coding sequences under positive selection ............................................... 32 Chapter 4: Conclusion ................................................................................................................ 37 Bibliography ................................................................................................................................ 42 Appendices ................................................................................................................................... 57 Appendix A Genome sequencing, assembly and annotation supplementary material .... 57 Appendix B Phylogenetic and comparative genomics analyses supplementary material 60 B.1 Picea mariana specific orthogroup analysis supplementary material .................. 61 B.2 Positive selection analysis supplementary material .............................................. 63 ix List of Tables Table 3.1 P. mariana genome assembly and annotation metrics. ................................................ 18 Table 3.2 Assembly and annotation statistics for all spruce and pine nuclear genomes used in downstream analyses. ................................................................................................................... 21 x List of Figures Figure 3.1 Assembly and annotation workflow. ........................................................................... 17 Figure 3.2 Assessment of P. mariana annotation quality (in bold face) relative to other Picea and Pinus taxa. ..................................................................................................................................... 22 Figure 3.3 MashTree phylogenies from the (A) nuclear, (B) chloroplast, and (C) mitochondrial genomes of P. mariana (in bold face) as well as five other spruce and two pine taxa. ................ 24 Figure 3.4 Distribution of (A) EggNOG functional categories and (B) best matching PFAM domains (n > 5). ............................................................................................................................ 28 Figure 3.5 Functional classes of the 45 positive selected protein-coding sequences. .................. 33 xi List of Symbols df degrees of freedom ω dN/dS xii List of Abbreviations CTAB cetyl trimethylammonium bromide EST expressed sequence tag HMW high molecular weight LRT likelihood ratio test LTR long terminal repeat MMP matrix metalloproteinases PET potential evapotranspiration RBH reciprocal best hit SNP single nucleotide polymorphism SPEI standardized precipitation evapotranspiration index TE transposable element TPM transcripts per million xiii Acknowledgements I would like to thank my supervisor Dr. Inanc Birol and my supervisory committee members, Dr. Joerg Bohlmann and Dr. Steven Hallam for the feedback and opportunities provided throughout the course of my degree. Special thanks to Dr. Jean Bousquet and Dr. Nathalie Pavy for their guidance and invaluable discussions as I, one who had little experience with plant biology in general, navigated the results of my analyses. Additionally, I would like to thank my colleagues at the BTL lab, particularly Lauren Coombe, Rene L. Warren, and Dr. Kristina K. Gagalova for all their support and suggestions. This project was supported through the CanSeq150 Initiative lead by S.J.M. Jones, and the Spruce-Up project co-lead by J. Bohlmann and J. Bousquet with funding from Genome Canada, Genome BC and Genome Quebec (grant number 243FOR). xiv Dedication To my family, thank you for your attempts to help decipher my thoughts whenever I was thinking out loud. I still very much appreciated the blank stares and occasional nods! All jokes aside, this would not have been possible without your endless support.1 Chapter 1: Introduction 1.1 Climate change and projected impacts on the boreal forest The boreal forest is one of the largest biomes on Earth, encompassing approximately 30% of the global forest area with Canada’s boreal forests comprising 28% (Brandt et al. 2013). These forests are characterized by frequent natural disturbances including fires, insect outbreaks, and storms, as well as severe winters where freezing temperatures often occur for 6-8 months of the year (Gauthier et al. 2015). Due to these harsh conditions, boreal ecosystems have low plant diversity, typically dominated by cold-tolerant, fire-adapted coniferous trees such as those in the Picea (spruce) and Pinus (pine) genera (Gauthier et al. 2015; Kayes and Mallik 2020). These characteristics have enabled the boreal biome to be a globally important economic asset, accounting for 33% and 26% of all sawnwood and paper exports, respectively (Burton et al. 2010). Furthermore, boreal forest soils are a significant contributor to global carbon storage and consequently have been considered a major carbon sink (Pan et al. 2011), with most being stored in peatlands and permafrost soils attributed to the cold, wet conditions slowing organic matter decay (Tarnocai et al. 2009; Deluca and Boisvenue 2012). Unfortunately, boreal forests are one of the ecosystems most affected by climate change. In Canada, this is already evident—the annual average temperature increased by 1.7 °C from 1948-2016, which is just over twice the global average of 0.8°C, and precipitation increased by 20% between 1948 and 2012 (Zhang et al. 2019). Conditions will only worsen as climate models generated by Zhang et al. (2019) have predicted continued increase in temperatures across all provinces and seasons, with a change in mean annual temperature ranging from 1.5-2.3°C as we approach mid-21st century. Within a similar timeframe, rainfall is generally expected to increase 2 as well by less than 10%, though there is more temporal and spatial variation in comparison to temperature. These climate extremes will further exacerbate the severity of various abiotic and biotic stresses including drought, forest fires, and pest outbreaks. Precipitation and soil water content are not necessarily positively correlated as elevated temperature will enhance losses via evapotranspiration. Jeong, Sushama, and Khaliq (2014) assessed drought risk in regions all over North America based on the calculated standardized precipitation evapotranspiration index (SPEI), which assesses drought risk based on the difference between precipitation and potential evapotranspiration (PET). Despite the projected increases in precipitation, their study indicated increased drought risk in southern Canada, particularly in regions below 50°N latitude due to small increases in precipitation combined with elevated temperature and thus, PET. This, in addition to the climate extremes mentioned previously, will naturally lead to recurring forest fires. Although serotinous species such as Picea mariana (black spruce) and Pinus banksiana (Jack pine) have adapted to thrive after fires by heat-triggered cone opening, these species are still at risk as the increased frequency of fires prevent trees from reaching reproductive maturity and consequently, dispersal of seeds, prior to the onset of subsequent fires (Pinno et al. 2013; Baltzer et al. 2021). Similarly, increased mortality has been observed in species that were previously successful in avoiding the worst of biotic outbreaks. For instance, P. mariana is relatively resistant to Choristoneura fumiferana (spruce budworm) compared to other species such as Balsam fir, the main host tree of the spruce budworm, partly due to asynchrony between budworm emergence and black spruce bud burst; however, warmer temperatures have not only resulted in a northern shift of spruce budworm outbreaks, where black spruce resides, but have also accelerated black spruce bud development favouring synchrony with larvae (Bellemin-Noël et al. 2021). 3 1.2 Genomic approach to understanding climate change tolerance There are three possible outcomes when forest tree populations are faced with climate change: 1) migration by shifting ranges, typically along altitudinal or latitudinal gradients, to track their ecological niches, 2) persisting in their current location by adapting to the new conditions, or 3) extirpation, the local extinction of a species in a particular geographic area (Aitken et al. 2008). Studies have shown that the rate of environmental change is exceeding the migration capacity of many species (Dyderski et al. 2018; Aubin et al. 2018); thus, recent research is more focused on assessing adaptive capacity as well as identifying traits associated with adaptation and the genetic architecture underlying these traits. There are several mechanisms contributing to the adaptive capacity of a species (reviewed in Royer-Tardif et al. 2021). Boreal tree species exhibit large population sizes and high population-level genetic diversity (Aitken 2008), increasing the chance of genotypes within the population that are adapted to future conditions (Royer-Tardif et al. 2021). The high fecundity of these species also increases the probability of yielding adapted individuals. Given these favourable traits, trees inhabiting this region are expected to have high adaptive capacity. The classic approach for demonstrating local adaptation is the use of provenance trails–seeds are sourced from different regions and planted in a common area to reveal genetic variance. However, this method is often not feasible due to limited resources, time, and space needed for the trees to grow. More recently, as sequencing costs have decreased exponentially, methods have scaled up to surveying genetic variability in whole genomes across populations, most of which aim to identify patterns of adaptive DNA polymorphisms. The increasing accessibility of 4 sequencing technologies has also driven major advances in assembly and annotation tools, leading to the development of genomic resources for numerous economically important trees including those residing in the boreal ecosystem such as Norway spruce (P. abies), white spruce (P. glauca), Engelmann spruce (P. engelmannii), Sitka spruce (P. sitchensis), and hybrid interior spruce (P. glauca x engelmannii x sitchensis) (Nystedt et al. 2013; Gagalova et al. 2022). These resources can contribute greatly to studies attempting to infer the genetic basis of adaptation. An annotated reference genome can not only provide more accurate single nucleotide polymorphism (SNP) calls, as reads from various individuals can be mapped to the reference followed by removal of any erroneous reads, but can also locate the candidate adaptive variants relative to gene structures, facilitating interpretation of results. 1.3 Black spruce and existing genomic/transcriptomic resources Black spruce is a dominant species in Canada’s boreal forests and is considered one of the most important commercial species for its high-quality pulpwood (Mullin et al. 2011). In addition to traits contributing to the high adaptive capacity of boreal species mentioned previously, black spruce is a transcontinental species–growing in a broad band spanning from the cold, dry conditions of Alaska to the warm, wet maritime provinces. Its ability to thrive in these diverse environments suggests populations have developed growth traits genetically adapted to the local climate (Rossi and Bousquet 2014; Sniderhan et al. 2018; Napier et al. 2020). Studies aimed at identifying the genotypes could prove useful in selecting for and/or breeding individuals with desirable traits. 5 P. mariana, P. glauca, P. engelmannii, P. sitchensis and P. rubens are all important economic species in Canada. Besides P. mariana and P. rubens, annotated genome assemblies for the remaining three species and interior spruce, a hybrid complex of P. glauca and P. engelmannii (Gagalova et al. 2022), as well as extensive gene catalogs for both P. glauca (Rigault et al. 2011) and P. sitchensis (Ralph et al. 2008) are available. Of particular importance is the P. glauca gene catalog consisting of 23,589 unique full-length insert cDNAs, which have been utilized in various studies investigating genomic variation in P. glauca itself, but also other spruces. This includes the identification of 583 SNPs potentially playing adaptive roles in temperature and precipitation variation (Prunier et al. 2011) as well as five SNPs related to drought stress in P. mariana (Napier et al. 2020). Despite the advances made in these studies using resources from closely related species, there is still a need for black spruce resources as morphological and ecological differences between P. mariana and other spruce all point to species-specific traits that may not be entirely accounted for in other resources (Weng and Jackson 2000; Nicklen et al. 2021). To date, several P. mariana resources are available, though most are transcriptome-related including a high-quality expressed sequence tag (EST) resource consisting of 4,594 sequences (Mann et al. 2013), a 97,095 SNP catalog from exome capture with probes designed from a P. glauca transcriptome (Pavy et al. 2016), and a set of 24 polymorphic EST-based microsatellite DNA markers (Rajora and Mann 2021). However, these resources are generated from reduced-representation sequencing methods; thus, besides the inherent inability to include non-coding sequences, these resources are limited by missing genomic content. For instance, ESTs often account for approximately 60% of an organism’s genes due to the under-representation of rare transcripts (Rudd 2003), whereas the SNP catalog only represents transcripts in particular tissue types and developmental stages as a result of probe design (Parla 6 et al. 2011; Lowe et al. 2017). Therefore, future studies would benefit from utilizing these resources in combination with whole genome sequencing and assembly. 1.4 Thesis objectives This thesis aims to provide an annotated P. mariana genome assembly that will facilitate forest genomics research, enabling improved conservation and breeding strategies under climate change. More specifically, the objectives for my thesis work is as follows: 1. Provide an annotated P. mariana genome assembly 2. Utilize this novel genomic resource to: a. Explore the evolutionary relationships between P. mariana and five other spruce species b. Identify potential genes involved in stress response 7 Chapter 2: Methods 2.1 DNA extraction and sequencing High molecular weight (HMW) DNA was extracted from a black spruce needle sample (genotype 40-10-1) collected in Thunder Bay, Ontario (50° 57' 39.96'' N, 90° 27' 20.16'' E; elevation, 741 m) by Bio S&T (http://www.biost.com/, Montreal, QC, Canada) using the cetyl trimethylammonium bromide method and HMW genomic DNA extraction protocol as detailed in the Chromium Genome Reagent Kits Version 2 User Guide (PN-120229). Assessment of DNA integrity by pulsed field gel electrophoresis indicated DNA sizes were concentrated at the 20 KB to 250 KB range. A total amount of 60 g HMW DNA was dissolved in 10 mM Tris Cl (pH 8), then sent to Canada’s Michael Smith Genome Sciences Centre (https://bcgsc.ca/, Vancouver, BC, Canada) for library preparation and sequencing. Four Illumina-compatible libraries, two with estimated fragment sizes of 400 bp and the other two with estimated fragment sizes of 800 bp, were constructed. Following analysis of DNA on 8% PAGE gels, DNA fragments were excised and eluted from gel slices overnight at 4°C in 300 µl of elution buffer, then purified using a Spin-X Filter Tube (Fisher Scientific) and ethanol precipitation. Genome libraries were prepared using a modified paired-end tag protocol supplied by Illumina Inc and assessed using an Agilent DNA 1000 series II assay and Nanodrop 7500 spectrophotometer. Additionally, one microfluidic partitioned library was constructed using the 10x Genomics Chromium system (10x Genomics, Pleasanton, CA, USA). GEMs (Gel beads-in-EMulsion) were produced by combining DNA, Master Mix, and partitioning oil in the 10x Genomics Chromium 8 Controller instrument with the micro-fluidic Genome Chip (PN-120216) (10x Genomics). The DNA in each GEM underwent isothermic amplification as a barcode was added to each fragment. Barcoded fragments then underwent Illumina library construction (as per the Chromium Genome Reagent Kits Version 2 User Guide (PN-120229)). The resulting library was assessed for quality using the Agilent 2100 Bioanalyzer (Santa Clara, CA, USA) and a DNA 1000 assay. The resulting 10x Genomics and Illumina-compatible libraries were sequenced on an Illumina HiSeqX instrument yielding five lanes of paired-end 150 bp reads and four lanes of paired-end 250 bp reads, respectively (Table A.1). 2.2 Mitochondrial genome assembly Given that mitochondrial genome assemblies are available for a number of spruce taxa (Jackman et al. 2015, 2020; Sullivan et al. 2020; Gagalova et al. 2022), the P. mariana mitochondrial genome was also assembled here to be used for downstream analyses. Adapter-trimming was performed on all reads using Trimadap vr11 (Li 2014), then assembled with ABySS v2.1.0 (Jackman et al. 2017) at various k-mer sizes (k= 64, 72, 80, 88, 96, 104, 116, 122, 128) and k-mer multiplicity thresholds (kc=3, 4). The assembly with the highest NG50 (k=116; kc=3) was determined by QUAST v5.0.2 (Gurevich et al. 2013) and mitochondrial DNA sequences were extracted based on BWA-MEM alignments of the scaffolds to a reference interior spruce mitochondrial genome (GenBank accessions MK697696–MK697708) (Jackman et al. 2015). These sequences were then error-corrected with Tigmint v1.1.2 (Jackman et al. 2018), and 9 passed to ARCS v1.0.6 (Yeo et al. 2018) and LINKS v1.8.7 (Warren, Yang, et al. 2015) (m=4-20000; k=20; l=10; a=0.1) for scaffolding. 2.3 Nuclear genome assembly Longer read lengths are highly advantageous in genome assembly as they are able to span low complexity regions. Given sufficient coverage, it is possible to generate longer pseudo-reads from the available paired-end reads through read-merging, which attempts to fill the gap between each read pair and extend the reads outwards based on overlapping k-mers. Two rounds of read-merging were performed on the adapter-trimmed reads. The first round consisted of cascading Konnector runs (Vandervalk et al. 2015), where the 10x Genomics and Illumina HiSeq reads were merged using k values ranging from 115-75 and 235-75, respectively, both with a step size of -10. Reads that were unable to be merged in the first round were subjected to a second round of merging with abyss-mergepairs (Jackman et al. 2017). The longer pseudo-reads and any remaining unmerged reads were then assembled using ABySS v2.2.5 (k=96, 112, 128, 144, 160; kc=3, 4). To determine which assemblies to move forward with, the quality of each assembly was assessed by abyss-fac (Jackman et al. 2017). Of the assemblies with reconstruction values on par with the expected genome size (approximately 20 Gbp), those with the highest NG50 values were passed to ntJoin v1.0.3 (Coombe et al. 2020) to perform iterative assembly-guided scaffolding runs, each with the following parameters: no_cut=True; k=32; w=250; reference_weights=’2’. Following scaffolding, a round of misassembly correction was performed on the resulting assembly using Tigmint v1.1.2 with span=2 and 10x Genomics Chromium reads provided as evidence. This was followed by another round of scaffolding using ARCS v1.1.1 and LINKS v1.8.6 (c=3; l=3; a=0.9; z=3000; s=90). Gaps introduced 10 during scaffolding were filled with Sealer v2.2.3 (Paulino et al. 2015) (L=150; P=10l; k=75,85,95,105,115), yielding the final genome assembly. 2.4 Identification and annotation of protein-coding sequences Prior to annotating the genome assembly, repeat regions were masked using RepeatMasker v4.1.1 (Chen 2004). A custom P. mariana repeat library was constructed using both LTR_retriever v2.9.0 (Ou and Jiang 2018) and RepeatModeler v2.0.1 (Flynn et al. 2020), then supplied as input to RepeatMasker to supplement the RepBase v22.08 repeat library (Bao et al. 2015). Subsequently, gene models were identified in the repeat-masked genome assembly using BRAKER v2.1.6 (Brůna et al. 2021) (--softmasking –etpmode), providing both protein sequences and RNA-seq alignments as evidence as well as an AUGUSTUS model that was pre-trained with BUSCO v4.1.4 (Stanke et al. 2006; Simão et al. 2015) (--long). Specifically, a total of 3,463,432 unique proteins from OrthoDB v10 Viridiplantae database (Kriventseva et al. 2019), UniProtKB/Swiss-Prot plant entries (The Uniprot Consortium 2019), P. glauca manual annotations (Warren, Keeling, et al. 2015), and high-quality proteins from annotations found in at least three of four other spruce taxa (Gagalova et al. 2022) selected using a reciprocal best hit (RBH) approach were used (Table A.2). Reciprocal BLAST searches were performed between P. engelmannii, the North American spruce with the most annotations, and each of the other three North American spruces (Gagalova et al. 2022). Proteins found in at least two of the three reciprocal BLAST search results were considered common amongst the North American spruces and included as evidence. To compile the RNA evidence, publicly available P. mariana Illumina HiSeq 2100 bp reads derived from seedling tissue were obtained (SRA accessions: SRR9595774, SRR9595777) (Shao et al. 2019) (Table A.3) and subjected to adapter trimming as 11 well as quality filtering with fastp v0.23.1 (Chen et al. 2018). By default, reads with at least 40% of bases having a Phred quality score < 15 were filtered out. Contaminant filtering was performed on the remaining reads with BioBloom tools v2.3.3 (Chu et al. 2014), where Bloom filters were created from Aphid, Archaea, Bacteria, Fungi, Protozoa, and Viral reference sequences obtained from RefSeq (O'Leary et al. 2016). Reads that did not have any hits to any of the Bloom filters, and thus not identified as contaminants, were aligned to the repeat-masked genome assembly using HISAT2 v2.2.0 (Kim et al. 2015) with –max-intronlen 1000000 (Table A.4). Following BRAKER annotation, complete protein-coding sequences, defined as those with start and stop codons, containing introns  10 bp as identified by GAG v2.0.1 (Geib et al. 2018), were functionally annotated using EnTAP v0.10.8 (Hart at al. 2020) – an annotation pipeline that assigns functions based on similarity search hits to user-selected databases as well as EggNOG and InterProScan hits. EnTAP was run in protein mode (--runP) with P. mariana provided as the target species and Aphid, Bacteria, and Fungi as contaminant taxa. Similarity searches were performed against the Gymno 1.0 PLAZA database (Proost et al. 2009), OrthoDB v10 Viridiplantae database, UniprotKB/Swissprot plant entries (The Uniprot Consortium 2019), and Uniref90 database (Suzek et al. 2015). In addition to these, the UniProtKB/TrEMBL database (The Uniprot Consortium 2019) was also included for contamination-screening purposes. Annotations with similarity search and/or ontology hits were assessed for the presence of Pfam domains found in gag and pol genes. Corresponding annotations were removed on the basis that the internal region of all long-terminal repeat (LTR) retrotransposons consists of these two genes 12 and thus, were likely missed by RepeatMasker. Furthermore, as fragmentation of a gene over multiple scaffolds could be common in draft assemblies, fragmented genes were identified and flagged if the following criteria were met: 1) TPM < 1, as determined by SALMON v1.3.0 (Patro et al. 2017) and 2) in at least two sets of read pairs, the paired-end reads mapped to genes on different scaffolds (Table A.5). All 883 flagged fragments were grouped into genes based on shared read pair mappings and the longest fragment of each group was then selected as the representative gene fragment for downstream comparative genomics analyses. The quality of the final assembly and annotation were assessed using abyss-fac and BUSCO v5.2.2 with the OrthoDB v10 embryophyta dataset in protein mode. 2.5 Phylogenetic analyses Three phylogenetic trees were constructed using P. mariana genome sequences and publicly available nuclear, chloroplast, and mitochondrial genome sequences of five other spruce taxa [Norway spruce (P. abies), white spruce (P. glauca), Engelmann spruce (P. engelmannii), Sitka spruce (P. sitchensis) and the tri-hybrid interior spruce (P. glauca x engelmannii x sitchensis)], and two pines [loblolly pine (P. taeda) and sugar pine (P. lambertiana)], using MashTree v1.2.0 (Katz et al. 2019) in bootstrap mode with the following parameters: mashtree_bootstrap.pl –reps 100 –numcpus 24 *.fa -- -–genomeSize [nuclear: 20000000000, chloroplast: 120000, mitochondrion: 5000000] > *.dnd (Table B.1). The resulting trees were visualized with MEGA11 (Tamura et al. 2021) and rooted using the pine species as the outgroup. 13 2.6 Comparative genomics analyses The use of different workflows can yield different annotations and thus, varying results in downstream analyses (Venkatraman et al. 2021). To account for this, the annotation files for P. mariana and the aforementioned taxa (P. abies, P. glauca, P. engelmannii, P. sitchensis, interior spruce, P. taeda, and P. lambertiana) were subjected to the same filtering steps prior to comparative genomics analyses – only complete genes with lengths  1 kbp and introns  10bp were kept. 2.6.1 Analysis of P. mariana specific orthogroups The longest transcript per gene was supplied to OrthoFinder v2.5.4 (Emms and Kelly 2019) and run with default settings to identify phylogenetic hierarchical orthogroups, referred to as orthogroups hereinafter. To verify the orthogroups that appeared specific to P. mariana, MMseqs2 v14.7e284 (Steinegger and Söding 2017) was used to perform RBH analysis on transcripts within these orthogroups against the high-quality transcripts annotated in four North American Picea taxa (Gagalova et al. 2022). Those with RBHs in at least one of the four taxa was removed from analyses pertaining to P. mariana specific orthogroups. 2.6.2 Identifying positively selected protein-coding sequences in P. mariana We used the ratio of nonsynonymous to synonymous substitutions (ω = dN/dS) to detect signatures of positive selection (Kimura 1977; Kryazhimskiy and Plotkin 2008). Given that synonymous mutations tend to be neutral or nearly neutral, when dN > dS, natural selection is favouring changes in protein-coding sequences more than neutral expectations and thus, the fixed mutations likely provide a fitness advantage (Kimura 1977). For instance, such changes 14 would be beneficial in stressful environmental conditions, such as those imposed by climate fluctuations (Hoffmann and Hercus 2000). Detecting signals of positive selection along a lineage based on pairwise estimations of ω is reputed to be a stringent procedure given that positive selection is only detected if the ω averaged over all amino acid sites of a protein is > 1 (Yang and dos Reis 2011). An alternative and more powerful approach is the use of branch-site tests as it allows variation in ω among branches as well as sites and thus, permits detection of positive selection acting on a few amino acids, which often is the case (Zhang et al. 2005). Given a multiple sequence alignment of genes and a species tree with a prespecified set of branches to test for the presence of positive selection on, referred to as foreground branches hereinafter, the parameters of null and alternative models are estimated. The difference between the two models pertains to the predefined foreground branches, where the null model does not allow sites with dN/dS > 1 on any branches while the alternative model does. A likelihood ratio test (LRT) is then performed to determine if the alternative model fits significantly better than the null model. Amongst the six spruce and two pine taxa, a total of 497 single-copy orthogroups were determined by OrthoFinder and used for this analysis. Primary sequence errors can contribute to misalignments and inaccurate branch length estimates, potentially resulting in overestimated ratios of nonsynonymous to synonymous substitution rates (ω = dN/dS) and consequently, false detection of positively selected genes (Schneider et al. 2009; Di Franco et al. 2019). To minimize the effects of such errors, non-homologous sequences and repeats were trimmed and/or masked in each single-copy orthogroup by PREQUAL v1.02 (Whelan et al. 2018) – a pre-15 alignment filtering tool. Due to the presence of annotation errors in four orthogroups, these could not be processed by PREQUAL and were excluded from the analysis. Multiple sequence alignments were generated for the remaining 493 orthogroups using the codon-aware aligner MACSE v2.06 (Ranwez et al. 2018), then provided as input with the previously generated MashTree species topology to CODEML, a program implemented in PAML v4.9j (Yang 2007). To test for positive selection in protein-coding sequences along the P. mariana lineage, branch-site models A1 (null model) and A (alternative model) were used with the branch leading to P. mariana set as the foreground. LRTs were conducted to compare the null and alternative model (at α = 0.05) by first calculating the LRT test statistic, then obtaining p-values from the cumulative chi-square distribution function with degrees of freedom, df = 1. 16 Chapter 3: Results and discussion 3.1 Genome sequencing, assembly and annotation Sequencing one library constructed using the 10x Genomics Chromium system, two Illumina-compatible libraries with an estimated fragment size of 400 bp, and another two Illumina-compatible libraries with an estimated fragment size of 800 bp yielded 2.2 billion 2150 bp reads, 252 million 2250 bp reads, and 341 million 2250 bp reads, respectively, totaling to a fold-coverage of approximately 46 (Table A.1). All reads were used in the assembly and initial scaffolding stages. Various AbySS k and kc parameter combinations were explored, optimizing for the NG50 length metric. The resulting assemblies were merged with subsequent, iterative ntJoin runs, leveraging the contiguity of an assembly from one iteration for the next. Contiguity evaluation of each assembly indicated that the AbySS assembly with k=176 and kc=3 had one of the highest NG50 values (7.2 kbp) and thus, was used as the target for ntJoin runs. Many of the remaining kc=3 assemblies had relatively high NG50s in comparison to those assembled with kc=4 and thus, were provided as references. The long-range information in linked reads was used by Tigmint and ARCS for misassembly correction and further scaffolding (Figure 3.1), ultimately yielding a final assembly consisting of 1,560,767 contigs with a NG50 length of 36.0 kbp and reconstruction of 18.3 Gbp (Table 3.1). 17 Figure 3.1 Assembly and annotation workflow. Using both Illumina HiSeq and 10x Genomics Chromium reads, the P. mariana genome was assembled using the pipeline shown on the left. Repeats were masked in the resulting assembly before being passed through the annotation pipeline shown on the right, yielding protein-coding gene models and corresponding functional annotations. Illumina HiSeq reads10x Genomics Chromium readsKonnectorMerged?abyss-mergepairs ABySSntJoinTigmintARCSLINKSScaffold-level assemblyRepeatMaskerBRAKER2EnTAPPost-processingGene annotationsProtein evidenceRNA-seq evidenceAUGUSTUS modelProtein databasesYesNo18 Table 3.1 P. mariana genome assembly and annotation metrics. Statistic Value Assembly size (Gbp) 18.3 Number of scaffolds (> 500 bp) 1,549,050 Scaffold N50 length (bp) 45,468 Scaffold NG50 length (bp) 35,958 Repeat content (%) 79.1 LTR elements (%) 49.3 Number of protein-coding sequences 66,332 Number of transcripts 68,738 Mean protein-coding sequence length (bp) 4,165 Mean CDS length (bp) 827 Mean exon length (bp) 256 Mean intron length (bp) 1,612 Conifer genomes have a high proportion of repetitive elements, particularly transposable elements (TEs), with the majority of these being LTRs (De La Torre et al. 2014). As the presence of these repetitive sequences can lead to false evidence for annotations, it is standard practice to repeat-mask the genome prior to annotation (Yandell and Ence 2012). By supplying a P. mariana custom repeat library to RepeatMasker, 79.1% of the genome was identified as repetitive sequences whereas 49.3% consisted of LTR elements. This is consistent with findings in other spruce genomes, with reported proportions greater than 70% of highly repetitive DNA (Gagalova et al. 2022). The repeat-masked genome was then passed to BRAKER followed by EnTAP for identification of gene models and functional annotation, respectively (Figure 3.1). A 19 total of 66,332 protein-coding sequences and 68,738 transcripts were annotated. Evaluation of annotation quality with BUSCO indicated a total of 416 (25.7%) complete BUSCOs in the annotation. Compared to the published nuclear genome assemblies and annotations of five other spruce and two pine species (Table B.1), the P. mariana assembly had a lower contiguity and number of complete BUSCOs than most of them (Table 3.2). However, this was expected due to the different sequencing technologies and coverage of data that was used for the various assemblies. Whereas the more contiguous genomes were assembled using reads with 80-110 fold coverage and, in some cases, benefitted from the use of long reads, P. mariana was assembled with short and linked reads sequenced at 46x fold coverage (Gonzalez-Ibeas et al. 2016; Zimin et al. 2017; Gagalova et al. 2022). Despite this limited sequencing coverage, the reconstruction of the P. mariana genome is consistent with the genome sizes of other spruce species, which were reported to be approximately 20 Gbp (Gagalova et al. 2022). The numbers of annotated protein-coding sequences and transcripts must be interpreted with caution when compared with the number of genes identified in previous spruce genome sequence studies given the different annotation pipelines and downstream filtering methods used to obtain the final annotations. For instance, BRAKER was used for predicting gene models in P. mariana, whereas a combination of AUGUSTUS and EuGene were used for P. abies and the MAKER-P pipeline was used for the other North American spruce species (Nystedt et al. 2013; Warren, Keeling, et al. 2015; Gagalova et al. 2022). Furthermore, the presence of frequent long introns interspacing exons in spruce genes (Nystedt et al. 2013; Stival Sena et al. 2014) coupled with the lower NG50 in the present study compared to Gagalova et al. 2022, can lead to the annotation of partial protein-20 coding sequences. Nonetheless, the higher number of annotated protein-coding sequences, despite applying filters such as completeness (presence of start and stop codons) and minimum intron length of 10 bp (Table A.5), was matched by a higher number of complete BUSCOs (Table 3.2, Figure 3.2). 21 Table 3.2 Assembly and annotation statistics for all spruce and pine nuclear genomes used in downstream analyses. NG50 was calculated using abyss-fac with -G 20e9. The fourth and seventh columns indicate the quality of the assemblies and annotations, respectively, as assessed by BUSCO v5. The second last column refers to the number of complete protein-coding sequences with lengths  1 kbp and intron lengths  10 bp. Species and genome assembly version NG50 length (bp) Reconstruction (Gbp) Number of complete BUSCOs (-m genome) Number of protein-coding sequences Number of transcripts Number of complete BUSCOs (-m protein) Number of protein-coding sequences used for analysis Reference Picea mariana v1 35,958 18.27 444 (27.4%) 66,332 68,738 416 (25.7%) 35,719 This study Interior spruce v5 121,714 20.14 796 (49.4%) 28,943 58,801 291 (18.0%) 28,878 Gagalova et al. 2022 Picea abies v1 500* 9.93 570 (35.3%) 70,736 178,049 176 (10.9%) 23,233 Nystedt et al. 2013 Picea engelmannii v1 355,449 20.75 736 (45.6%) 35,226 60,224 283 (17.6%) 35,226 Gagalova et al. 2022 Picea glauca v2 131,339 21.58 784 (48.6%) 30,140 55,309 290 (18.0%) 30,140 Gagalova et al. 2022 Picea sitchensis v1 38,458 18.22 603 (37.4%) 30,324 57,709 294 (18.2%) 30,324 Gagalova et al. 2022 Pinus lambertiana v1.5 2,702,679 19.82 1012 (62.7%) 38,518 38,518 1,183 (73.3%) 31,547 Gonzalez-Ibeas et al. 2016 Pinus taeda v2.01 103,592 20.43 819 (50.7%) 51,751 51,751 673 (41.7%) 37,499 Zimin et al. 2017 * Reported value is the not the true NG50 value as the reconstruction of the genome is significantly lower than the expected genome size of 20e9.22 Figure 3.2 Assessment of P. mariana annotation quality (in bold face) relative to other Picea and Pinus taxa. The reads used for each genome assembly is in parentheses, where “hybrid” refers to the use of both short and long reads. Annotation quality was determined using BUSCO v5 in protein mode. 3.2 Phylogenetic analyses To gain insights on the evolutionary relationships of P. mariana in relation to other spruce species, phylogenetic analysis was performed with six spruce and two pine species for which complete genome sequences were available (Table B.1), using pines as outgroups for these analyses. Phylogenetic trees were also constructed using chloroplast and mitochondrial genome sequences. Given the uniparental mode of inheritance of organellar genomes in Pinaceae (paternal for the chloroplast genome and maternal for the mitochondrial genome), their P. taeda v2.01(hybrid)P. lambertiana v1.5(hybrid)P. sitchensis v1(hybrid)P. glauca v2(short reads)P. engelmannii v1(hybrid)P. abies v1(short reads)Interior spruce v5(short reads)P. mariana v1(short reads)0 20 40 60 80 100%BUSCOs Complete (C) and single−copy (S) Complete (C) and duplicated (D) Fragmented (F) Missing (M)BUSCO Assessment Results23 phylogenetic tree topologies are observed to be different, in concordance with previous studies based on specific chloroplast and mitochondrial regions (Bouillé et al. 2011; Ran et al. 2015; Gagalova et al. 2022). The use of chloroplast, mitochondrial, and biparentally-inherited nuclear genomes should thus provide a more complete picture of evolutionary relationships among these six spruce species. In the nuclear phylogeny, P. mariana together with the European P. abies, are sister groups to the other North American spruces included in this analysis, namely, P. glauca, P. engelmannii, P. sitchensis and the hybrid interior spruce (Figure 3.3A). Whereas P. glauca, P. engelmannii, and P. sitchensis frequently hybridize with each other, giving rise to interior spruce (Gagalova et al. 2022), attempted crosses between these taxa and P. mariana have either failed or had low success rate, indicative of higher evolutionary divergence of P. mariana relative to the others (Ran et al. 2006). However, the position of P. mariana relative to P. abies remained ambiguous given the weak support of P. abies as a sister group to P. glauca, P. engelmannii, P. sitchensis and interior spruce. 24 Figure 3.3 MashTree phylogenies from the (A) nuclear, (B) chloroplast, and (C) mitochondrial genomes of P. mariana (in bold face) as well as five other spruce and two pine taxa. MashTree estimates distances based on k-mer sketches of the genomes, then constructs phylogenetic trees using the neighbour-joining method. For all three phylogenies, the trees were rooted using the pine species as outgroups. Branch lengths are proportional to the number of substitutions per nucleotide position. Values at the nodes indicate bootstrap values. A similar phylogenetic tree topology was obtained with the chloroplast and mitochondrial genomes (Figure 3.3), but the positions of P. sitchensis and P. abies conflicted among phylogenies. In the chloroplast phylogeny (Figure 3.3B), the position of European P. abies as a sister group to North American P. glauca, P. engelmannii, P. sitchensis and interior spruce was BACPicea marianaPicea marianaPicea mariana0.02 nucleotide substitutions per site0.01 nucleotide substitutions per site0.025 nucleotide substitutions per site25 weakly supported. However, in the mitochondrial phylogeny (Figure 3.3C), there was strong support for P. mariana as a sister group to these other North American species. Furthermore, in the mitochondrial phylogeny and contrary to the chloroplast phylogeny, there was strong support for P. sitchensis as a sister group to P, glauca, P engelmannii and interior spruce, with the latter confirmed as a natural tri-hybrid of the former three species (Gagalova et al. 2022). This is indicative of incomplete reproductive isolation and reduced interspecific genetic divergence among them. The topology of the chloroplast phylogeny placed P. sitchensis as an outgroup to all other spruce taxa, conflicting with the two other phylogenies and biogeography. This odd positioning had previously been observed and proposed to result from ancient reticulate evolution affecting the P. sitchensis chloroplast genome (Bouillé et al. 2011; Sullivan et al. 2017; Gagalova et al. 2022). Though weakly supported, the odd positioning of the European P. abies as a sister group to P. glauca and P. engelmannii on the chloroplast phylogeny, and to P. glauca, P. engelmannii and P. sitchensis on the nuclear phylogenetic tree, might also be indicative of reticulation between P. abies and P. glauca through ancient gene flow (Chen et al. 2010). Such contact could date back to an interglacial period approximately 400,000 years ago when spruce species were dominant in Greenland, in particular P. abies (de Vernal and Hillaire-Marcel 2008). Such a tree topology was not observed on the mitochondrial phylogeny where a strong geographic structure placing P. mariana as the sister group to all other North American spruces was observed and expected from ancient geographical speciation (Bouillé et al. 2011). Indeed, spruce mitochondrial genomes are dispersed by spruce seeds only, which disseminate across much smaller distances than pollen and thus, chloroplast genomes. As such, the mitochondrial tree may be more indicative of vertical 26 descent associated with phylogeographic speciation, while the chloroplast tree would integrate horizontal transfers by pollen though more or less ancient reticulate evolution between already genetically distinct lineages (Gérardi et al. 2010; Bouillé et al. 2011; Sullivan et al. 2017; Gagalova et al. 2022). As for the nuclear phylogeny, given the biparental inheritance of the nuclear genome, it well reflects a blend of both types of evolution. 3.3 Comparative genomics 3.3.1 P. mariana specific orthogroup analysis In a recent study, Gagalova et al. (2022) showed that the allopatric and ecologically divergent P. glauca and P. sitchensis had mostly distinct sets of rapidly evolving genes under positive selection. Despite being distinct, the genes were largely related to stress and stimuli response. Although the natural ranges of P. mariana and P. glauca are mostly sympatric, they are reproductively isolated and have adapted to distinct climatic and ecological niches. Black spruce is often found on cold, wet, nutrient-poor soils, whereas white spruce tends to inhabit sites with warmer, well-drained soils (Nicklen et al. 2021). Further delineating and investigating protein-coding sequences found specific to P. mariana may provide insights on how it has adapted to its rather specific ecological niche among the spruce species analyzed here. To compare the protein-coding sequence content among the six spruce and two pine species, and identify those specific to P. mariana, an orthogroup analysis was performed using OrthoFinder. In the context of the present taxa sampling, this analysis resulted in the identification of 34,776 orthogroups in total, with 2,869 being shared amongst all species and 560 found specific to P. mariana (Table B.2). Of the 560 P. mariana specific orthogroups, 109 of the 1,900 protein-27 coding sequences had RBHs with genes in at least one of the four North American spruces and thus, were removed from subsequent analyses. On further analysis of the remaining P. mariana specific protein-coding sequences, 792 were assigned functional categories based on orthologous hits in the EggNOG database. Majority of the predicted functions were unknown (Figure 3.4A), with many protein-coding sequences containing a domain of unknown function (DUF4283); however, the EamA and F-box domains were also highly prevalent (Figure 3.4B, Figure B.1) and these gene families have been partially characterized in Viridiplantae. The EamA gene family is involved in transport (Västermark et al. 2011) and sequences belonging to the F-box superfamily are responsible for controlling diverse biological processes including growth, development, and abiotic stress tolerance (Abd-Hamid et al. 2020; Rao et al. 2021). 28 Figure 3.4 Distribution of (A) EggNOG functional categories and (B) best matching PFAM domains (n > 5). Functional category annotations were assigned based on hits to entries in the EggNOG database, while PFAM annotations were assigned via hits in EggNOG and/or InterPro databases as part of the EnTAP annotation pipeline. If a protein-coding sequence was assigned to multiple functional categories, each category was counted separately. 29 Apart from those classified as function unknown, there was a relatively large presence of sequences in the following categories: signal transduction mechanisms; post-translational modification, protein turnover, chaperones; transcription; carbohydrate transport and metabolism; secondary metabolites biosynthesis, transport and catabolism (Figure 3.4A). Among those classified under signal transduction mechanisms, Pkinase domains were particularly prevalent (Figure 3.4B, Figure B.1). NB-ARC and LRR 1 domains, both found in NLR proteins which are widely known to play a central role in resilience to biotic and abiotic stresses (Van Ghelder et al. 2019; Ence et al. 2022), were present in high numbers as well. Interestingly, a study aimed at developing a repertoire of conifer NLR genes identified more in P. mariana compared to P. glauca (Van Ghelder et al. 2019), thereby supporting our finding of numerous NB-ARC and LRR 1 domain-containing sequences unique to P. mariana. Additionally, Lectin LegB domains were often found in sequences associated with signal transduction mechanisms. Though these domains are also found in L-type LecPs and LecRPs, they have been predominantly identified and characterized in L-type LecRLKs (Bellande et al. 2017), with various studies implicating their critical role in pathogen response (Singh and Zimmerli 2013; Lannoo and Van Damme 2014) and a few identifying links to abiotic stress response and developmental processes (Wan et al. 2008; Li et al. 2014; Zhao et al. 2019). Domains found in sequences annotated with the post-translational modification, protein turnover, chaperones EggNOG functional category included PG binding 1, AAA, and UCH (Figure 3.4B, Figure B.1). Many sequences containing the PG binding 1 domain were annotated with an InterPro domain in matrix metalloproteinases (MMPs) – Pept M10 metallopeptidase. Few plant MMPs have been characterized in detail, but they are known to remodel the extracellular matrix 30 during growth, development, and stress response (Marino and Funk 2011). This has also been shown in conifers. In particular, a recent study assessing genetic variation among drought resilient Norway spruce trees identified a SNP in a MMP family gene that had significant association with wood density (R2 = 0.25-0.26) (Trujillo-Moya et al. 2018). In contrast, the AAA superfamily has much broader functionality, ranging from proteolytic degradation, molecular remodelling, and homologous recombination (Snider et al. 2008). Myb Dna-binding domain was the most frequent domain in the transcription EggNOG functional category (Figure 3.4B, Figure B.1), accounting for 11 of the 58 sequences assigned to that category. MYB transcription factors are divided into clades based on the number of MYB repeats present. Due to the high prevalence of R2R3-type MYB domain proteins in plants (Bedon et al. 2007), research has been focused on this group. Previously, eight genes were identified in P. mariana and divided into three broad subfamilies based on homology to previously characterized genes – anthocyanin or proanthocyanidin biosynthesis, phlobaphene biosynthesis, and regulation of epidermal cell development (Xue et al. 2003). Furthermore, tissue specific expression patterns of 18 R2R3-type MYBs identified in P. glauca suggested roles in lignin biosynthesis as well as terpene and flavonoid synthesis (Bedon et al. 2007). Over 120 R2R3-MYB genes have been reported in Angiosperms (Riechmann et al. 2000). Given the large size of the R2R3-MYB family, it is likely that those previously identified make up only a fraction of the R3R3-MYB clade in spruce, let alone those that have yet to be identified in other MYB clades. Following Myb Dna-binding, the Med26 domain, which plays a possible role in transcription elongation, (Mathur et al. 2011) and RNA pol Rpb2 2, the domain characteristic of 31 genes encoding the second-largest subunit of DNA-dependent RNA polymerase II, also appeared frequently in the transcription EggNOG functional category. Protein-coding sequences related to carbohydrate transport and metabolism were often categorized under energy production and conversion (Figure B.1). Within these two EggNOG functional categories, UDPGT, Alde H, RuBisCO large, and MFS 1 domains were the most represented (Figure 3.4B, Figure B.1). UDPGT is a domain often found in UDP-glycosyltransferases, which have been shown to play a role in cell wall lignification (Lin et al. 2016). The p450 domain was most prevalent in the secondary metabolite biosynthesis, transport and catabolism EggNOG functional category (Figure 3.4B, Figure B.1). Cytochrome P450 oxygenases (CYPs) have been well studied in plants due to their important role in expanding the structural diversity of terpenes, plant specialized secondary metabolites that provide protection against both biotic and abiotic threats (Bathe and Tissier 2019). This includes important roles of CYP genes in oleoresin defenses of conifers and in interactions of conifers with insect pests and pathogens (Celedon and Bohlmann 2019; Chiu and Bohlmann 2022). The CYP superfamily is one of the largest family of enzymes in plants with eleven CYP clans identified in vascular plants (Bak et al. 2011), each of which likely include expansions within subfamilies as observed in the CYP71 and CYP85 clans (Warren, Keeling, et al. 2015). Therefore, although CYP genes have been annotated in spruces previously (Warren, Keeling, et al. 2015), the p450 domain-containing P. mariana specific protein-coding sequences may be CYP enzymes that have not been annotated in the species included in this analysis. 32 In all, 999 P. mariana specific protein-coding sequences were not assigned functional categories as there were no orthologous hits in the EggNOG database, but 122 of these sequences were annotated with PFAM domains based on InterPro hits. The most prevalent PFAM domain amongst these was the TCP domain, which has been found in transcription factors regulating plant growth and development (Cubas et al. 1999). 3.3.2 Protein-coding sequences under positive selection To identify positively selected protein-coding sequences in P. mariana, branch-site tests were conducted on the 497 single-copy orthogroups identified between all eight conifer taxa considered in this study. Of these, 45 (9.1%) were identified as positively selected (Table B.3), 38 of which could be annotated with one or more functional categories based on matches with known proteins in the Araport11 database (Cheng et al. 2017) and/or literature search highlighting a diverse set of molecular functions (Table B.4, Table B.5). The most represented functional category was development (12 sequences), followed by response to stress (eight sequences) and growth (seven sequences) (Figure 3.5). Among them, several protein-coding sequences appear essential either for plant survival or plant adaptation to environmental conditions as illustrated by the following examples. 33 Figure 3.5 Functional classes of the 45 positive selected protein-coding sequences. In total, 38 sequences had a match with an Arabidopsis protein from the Araport11 database (blastp evalue < e-20) and could be associated to a functional class based on homology searches and/or literature searches. If a protein-coding sequence was assigned to multiple functional classes, each class was counted separately. Whereas “not assigned” refers to sequences without hits in the Araport11 database, “uncharacterized” refers to those that had matches but the functions of those Arabidopsis proteins have not been fully characterized yet. Twelve protein-coding sequences are involved in development, seven of which are only assigned to that category. These include two sequences whose homologs regulate flowering time in Arabidopsis: the developmental protein FRIGIDA (Michaels et al. 2004), also involved in temperature adaptation (Tabas-Madrid et al. 2018), and the serine/threonine-protein kinase WNK1, which regulates flowering time by modulating the photoperiod pathway (Wang et al. 024681012DevelopmentStress responseGrowthNot assignedMetabolismRegulationUncharacterizedMitotic cell cycleRNA processingRedox regulationFunctional categoryNumber of protein−coding sequences34 2008). Furthermore, there is an aspartate/glutamate/uridylate kinase family protein that has been identified to be a trihelix/aa-kinase chimera associated with leaf development (Kuromori et al. 2006; Kaplan-Levy et al. 2012). The other four sequences include an ARID domain-containing protein, which is a seed-specific transcription factor in Arabidopsis (Zheng et al. 2014), a NAC transcription factor involved in xylem formation (Endo et al. 2015), ARF19, which is responsible for regulating various auxin-mediated development processes including lateral root formation (Okushima et al. 2005, 2007), and TMKL1, whose function has yet to be fully characterized but has been suggested to mediate vascular tissue development (Wu et al. 2016). The remaining five of 12 sequences involved in development were also associated with other functional classes. One protein-coding sequence plays a role in RNA processing as well as development – RRD1 participates in mitochondrial mRNA deadenylation, which is fundamental in controlling early lateral root organogenesis (Otsuka et al. 2021). There are also two protein-coding sequences involved in both growth and development. One is an ankyrin repeat-containing protein essential for chloroplast biogenesis (Shen et al. 2010), while the other is a member of the IQ67-domain (IQD) protein family known to regulate plant growth, lateral organ polarity, and basal defense response against environmental cues (Bürstenbinder et al. 2017; Barda and Levy 2022). The last two sequences are involved in stress response and development: a DCD (Development and Cell Death) domain-containing protein (Tenhaken et al. 2005) and AAC2, which encodes a mitochondrial ADP/ATP carrier that has been suggested to play a role in the mechanisms of ABA-mediated stress response (Kharenko et al. 2011). 35 A total of eight sequences encode proteins related to stress response. Three are involved in plant immunity or plant immunity regulation: a member of the pleiotropic drug resistance family, the PBL27 serine-threonine kinase involved in a cascade that leads to chitin-induced immunity (Kawasaki et al. 2017) and a homolog of the AtWRKY3 transcription factor regulating the expression of the pathogen induced gene npr1 (Yu et al. 2001). PLIP2 is a glycerolipid A1 lipase that provides tolerance to various abiotic stresses, including cold, by ABA-mediated synthesis of jasmonate and oxolipins (Wang et al. 2018). Four of the eight stress-related protein-coding sequences are assigned to another functional category along with stress response, two of which have been mentioned previously – a DCD domain containing protein and AAC2, both contributing to plant development as well. Besides those, a sequence encoding DGD1 is additionally annotated with metabolism as it is responsible for the synthesis of DGDG, a glycolipid critical for the stabilization of chloroplast membranes, thereby conferring thermotolerance (Kobayashi et al. 2014). APUM8 is also annotated with RNA-processing and has been implicated in salt stress tolerance (Huang et al. 2018). Five protein-coding sequences have direct or indirect roles in growth, cell elongation and plant height in Arabidopsis: a filament-like protein (Chen et al. 2016), an actin-binding FH2 (formin 2) family protein (Vidali et al. 2009), DRACULA2 (Gallemí et al. 2016), CNGC8, which has been shown to be essential for pollen tube tip growth (Tunc-Ozdemir et al. 2013), and KORRIGAN2 (Mølhøj et al. 2001), whose homolog is also associated to early growth in pine (Cabezas et al. 2015). 36 The functional categories of the high dN/dS protein-coding sequences found in the present study overlap well with those reported by Buschiazzo et al. (2012) and De La Torre et al. (2015), both of which use estimations of dN/dS to detect positive selection. Eight gene families implicating high dN/dS values were found in common between the present study and that of De La Torre et al. (2015), including auxin responsive protein, eukaryotic aspartyl protease, alkaline and neutral invertase, ABC transporter, cyclic nucleotide gated channel (CNGC), FRIGIDA-like proteins, ankyrin repeat-containing protein, and the O-Glycosyl hydrolases family 17. However, the specific protein-coding sequences with high dN/dS values were different between studies, including the fact that here, we considered only single-copy genes in the dN/dS analysis. These studies should be viewed as complementary, two being based on transcriptome data and our study on genomic data. Furthermore, there was a large overlap between the abundant functional categories found in this study and those presented in Gagalova et al. (2022), where the most rapidly evolving genes shared by P. sitchensis and P. glauca were identified by estimating the ratios of nonsynonymous to synonymous SNPs (SNP A/S) for each species. Similar to the branch-site tests used in this thesis, the estimation of SNP A/S values also aims to identify signals of selection by comparing synonymous and non-synonymous substitutions. The main difference between the two methods is that the former detects selection by examining changes between species whereas the latter relies on variation within the same species; thus, the two methods are complementary with branch-site tests best suited for detecting selection occurring over longer time periods and the SNP A/S approach for detecting recent selection (Aguileta et al. 2008). 37 Chapter 4: Conclusion The recent surge in genome assemblies for numerous economically important forest trees such as P. engelelmannii, P. glauca, P. sitchensis (Gagalova et al. 2022) , and now, P. mariana, showcases the rapid improvement in this field; however, assembly of conifer genomes continues to be a challenge for current tools with few chromosome-scale assemblies achieved due to the large genome sizes, typically around 20 Gb, and high repetitive content characteristic of this taxonomic class (De La Torre et al. 2014). Indeed, improvements in read lengths and accuracy coupled with decreased costs of sequencing in the past decade have resulted in more contiguous assemblies. This is evident in the lower NG50 of the P. mariana assembly, which was obtained with short reads alone, compared to that of other Picea genome assemblies presented in Gagalova et al. (2022), which used a combination of short and long read technologies and had NG50s ranging from 38 kbp to 355 kbp. The true potential of longer reads is further showcased in the recent assemblies of giant sequoia, coast redwood, and Chinese pine, where the utilization of long-read sequencing and Hi-C scaffolding yielded chromosome-scale scaffolds (Scott et al. 2020; Neale et al. 2022; Niu et al. 2022). While significant advancements have been made in the accuracy and contiguity of conifer genome assemblies, relatively less progress has been made in annotating these genomes. Besides the propagation of assembly errors, which could lead to frame shifts and truncated or fused gene models, the high repeat content as well as large introns make the task of structural annotation challenging. Approximately 75% of conifer genomes are composed of repeat sequences (Ahuja and Neale 2005), with majority of these being TEs (De La Torre et al. 2014). TEs are particularly problematic for annotation tools as insertion of these elements in close proximity to genes may result in parts of the gene model being masked, thus affecting gene annotation. On the other hand, annotation tools may split a gene into several 38 truncated partial gene models if interspersed by long introns. MAKER and the more recently developed, BRAKER, are state of the art tools used for eukaryotic genome annotation. BRAKER has been shown to yield annotations with higher BUSCO completeness scores, sensitivity, and precision than MAKER when evaluated on published Arabidopsis and Populus annotations (Vuruputoor et al. 2023). Though not a direct comparison of pipelines as BRAKER was used for annotating P. mariana and the other North American spruce species were annotated with MAKER (Gagalova et al. 2022), similar results were observed—a larger number of genes were annotated in the P. mariana genome accompanied by a higher BUSCO completeness score. The difference in number of annotations is likely due to fundamental differences between the two tools. MAKER utilizes EST and/or assembled transcripts as evidence, potentially resulting in under-representation of rare transcripts and propagation of assembly errors, respectively, whereas BRAKER utilizes RNA-seq evidence. Besides, there are no current standard practices for post-processing, which is a necessary step to select for high confidence gene annotations, further contributing to discrepancies in annotations. Therefore, as sequencing technologies, assembly tools, and annotation pipelines continue to improve, related genomic resources will need to be updated to ensure high-quality data for various downstream applications. Nevertheless, the annotated genome assembly proves to be a valuable resource that can be leveraged for analyses. To date, there are two hypotheses on the origin of the Picea genus. Based on morphological trait analyses, Wright (1955) proposed that Picea originated in northeast Asia and the genus migrated to North America via the Beringia land bridge, which connected Siberia to Alaska in geological periods with lower sea levels. This is also supported by several phylogenies constructed from mitochondrial DNA (Lockwood et al. 2013). On the other hand, 39 most phylogenies based on chloroplast and nuclear molecular data support a North American origin for the genus (Ran et al. 2006, 2015; Bouillé et al. 2011). Unlike the geographic origin, phylogenetic relationships between lineages are still unresolved with conflicts arising depending on the combinations of nuclear and organelle DNA used (Lockwood et al. 2013; Ran et al. 2015). With the P. mariana nuclear genome now available, a unique opportunity to investigate the phylogenetic relationships alongside phylogenies based on chloroplast and mitochondrial genomes arose. Substantial genetic differentiation was detected between the black spruce nuclear and organellar genomes and those of other North American spruce species, reflecting the reproductive isolation of black spruce from the others and the more ancient divergence of its lineage. The unusual positioning of P. sitchensis as an outgroup to all other spruce taxa in the chloroplast phylogeny alone is characteristic of ancient reticulate evolution. Though phylogenies including more spruce taxa would provide further insight, similar results in other phylogenetic studies suggest that the chloroplast genome of a likely Asian spruce taxon introgressed P. sitchensis from long distance dispersal through airborne pollen (Bouillé et al. 2011; Sullivan et al. 2017; Gagalova et al. 2022) or by direct contact at a time where P. sitchensis could have been parapatric to an Asian taxon through the Beringia land bridge during or before the Pleistocene. Ancient reticulate evolution was also detected between P. abies and P. glauca, though weakly supported by the nuclear and chloroplast phylogenies. The occurrences of at least two signatures of ancient reticulate evolution of the spruce chloroplast and nuclear genomes supports evidence of incomplete lineage sorting (Bouillé and Bousquet 2005) and slow speciation characterized by delayed establishment of reproductive isolation influenced by the large effective population sizes and airborne dispersion of pollen over potentially long distances in these archaic tree species (Bousquet et al. 2021). 40 Analysis of specific biological functions and positive selection revealed an abundance of genes implicated in development and stress responses, whether for biotic or abiotic factors. This pattern is largely shared with observations from previous studies of molecular adaptation in spruces (e.g. Hornoy et al. 2015; Yeaman et al. 2016; Depardieu et al. 2021; Gagalova et al. 2022) and is congruent with the often-harsh environments or climates characterizing the ecological niches occupied by spruce species, particularly black spruce as it is found mostly in the boreal forests of North America. While the P. mariana specific protein-coding sequences identified could represent potential novel genes related to its unique ecological traits, it is also possible that some appear specific due to sequence divergence leading to failed inference of orthologs and/or unannotated sequences in the other spruce and pine species. In particular, within the MYB gene family, the amino acid sequences outside of the MYB DNA-binding domains are highly divergent (Bedon et al. 2010), which may be challenging for multiple sequence alignment tools and thus, complicate the process of gene tree and orthogroup inference. Furthermore, several highly abundant PFAM domains in P. mariana specific protein-coding sequences were also found to be abundant in unannotated P. lambertiana transcripts, namely DUF4283, PKinase, NB-ARC, WD40, and LRR 1 (Gonzalez-Ibeas et al. 2016). Therefore, as tools and conifer genome annotations continue to improve, further analyses will be needed to corroborate the initial findings presented as part of the P. mariana specific orthogroup and positive selection analyses. Ultimately, validation of whether these candidate genes genuinely play fundamental roles in development and stress response is necessary before establishing their biological function. 41 With black spruce predicted to decline due to increasing temperature across its natural range (Thomson et al. 2009), and the increasing frequency of extreme weather events such as late cold spells and drought episodes under mid-northern latitudes (Benomar et al. 2022; Laverdière et al. 2022), much remains to be understood about the specific mechanisms and molecular basis of genetic adaptation in this ecologically and economically important species. The resources and research presented in this thesis would aid the forest research community in gaining a better understanding of the different mechanisms and capacities for local adaptation to climate warming and increasing instability, thus contributing to the development of molecular-assisted breeding strategies and guidelines to orientate assisted migration. 42 Bibliography Abd-Hamid N-A, Ahmad-Fauzi M-I, Zainal Z, Ismail I. 2020. Diverse and dynamic roles of F-box proteins in plant biology. Planta. 251(3):68. doi:10.1007/s00425-020-03356-8. Aguileta G, Refrégier G, Yockteng R, Fournier E, Giraud T. 2009. Rapidly evolving genes in pathogens: methods for detecting positive selection and examples among fungi, bacteria, viruses and protists. Infect Genet Evol. 9(4):656–670. doi:10.1016/j.meegid.2009.03.010. Ahuja MR, Neale DB. 2005. Evolution of genome size in conifers. Silvae Genet. 54(1–6):126–137. doi:10.1515/sg-2005-0020. Aitken SN, Yeaman S, Holliday JA, Wang T, Curtis-McLane S. 2008. Adaptation, migration or extirpation: climate change outcomes for tree populations. Evol Appl. 1(1):95–111. doi:10.1111/j.1752-4571.2007.00013.x. Asaf S, Khan AL, Khan MA, Shahzad R, Lubna, Kang SM, Al-Harrasi A, Al-Rawahi A, Lee I-J. 2018. Complete chloroplast genome sequence and comparative analysis of loblolly pine (Pinus taeda L.) with related species. PLoS One. 13(3):e0192966. doi:10.1371/journal.pone.0192966. Aubin I, Boisvert-Marsh L, Kebli H, McKenney D, Pedlar J, Lawrence K, Hogg EH, Boulanger Y, Gauthier S, Ste-Marie C. 2018. Tree vulnerability to climate change: improving exposure-based assessments using traits as indicators of sensitivity. Ecosphere. 9(2):e02108. doi:10.1002/ecs2.2108. Bak S, Beisson F, Bishop G, Hamberger B, Höfer R, Paquette S, Werck-Reichhart D. 2011. Cytochromes P450. Arab Book. 9:e0144. doi:10.1199/tab.0144. Baltzer JL, Day NJ, Walker XJ, Greene D, Mack MC, Alexander HD, Arseneault D, Barnes J, Bergeron Y, Boucher Y, et al. 2021. Increasing fire and the decline of fire adapted black spruce in the boreal forest. Proc Natl Acad Sci USA. 118(45):e2024872118. doi:10.1073/pnas.2024872118. Bao W, Kojima KK, Kohany O. 2015. Repbase update, a database of repetitive elements in eukaryotic genomes. Mob DNA. 6(1):11. doi:10.1186/s13100-015-0041-9. Barda O, Levy M. 2022. IQD1 involvement in hormonal signaling and general defense responses against Botrytis cinerea. Front Plant Sci. 13:845140. doi:10.3389/fpls.2022.845140. Bathe U, Tissier A. 2019. Cytochrome P450 enzymes: a driving force of plant diterpene diversity. Phytochemistry. 161:149–162. doi:10.1016/j.phytochem.2018.12.003. Bedon F, Bomal C, Caron S, Levasseur C, Boyle B, Mansfield SD, Schmidt A, Gershenzon J, Grima-Pettenati J, Séguin A, et al. 2010. Subgroup 4 R2R3-MYBs in conifer trees: gene family expansion and contribution to the isoprenoid- and flavonoid-oriented responses. J Exp Bot. 61(14):3847–3864. doi:10.1093/jxb/erq196. 43 Bedon F, Grima-Pettenati J, Mackay J. 2007. Conifer R2R3-MYB transcription factors: sequence analyses and gene expression in wood-forming tissues of white spruce (Picea glauca). BMC Plant Biol. 7(1):17. doi:10.1186/1471-2229-7-17. Bellande K, Bono J-J, Savelli B, Jamet E, Canut H. 2017. Plant lectins and lectin receptor-like kinases: how do they sense the outside? Int J Mol Sci. 18(6):1164. doi:10.3390/ijms18061164. Bellemin-Noël B, Bourassa S, Despland E, De Grandpré L, Pureswaran DS. 2021. Improved performance of the eastern spruce budworm on black spruce as warming temperatures disrupt phenological defences. Glob Chang Biol. 27(14):3358–3366. doi:10.1111/gcb.15643. Benomar L, Bousquet J, Perron M, Beaulieu J, Lamara M. 2022. Tree Maladaptation Under Mid-Latitude Early Spring Warming and Late Cold Spell: Implications for Assisted Migration. Front Plant Sci. 13:920852. doi:10.3389/fpls.2022.920852. Bouillé M, Bousquet J. 2005. Trans-species shared polymorphisms at orthologous nuclear gene loci among distant species in the conifer Picea (Pinaceae): implications for the long-term maintenance of genetic diversity in trees. Am J Bot. 92(1):63–73. doi:10.3732/ajb.92.1.63. Bouillé M, Senneville S, Bousquet J. 2011. Discordant mtDNA and cpDNA phylogenies indicate geographic speciation and reticulation as driving factors for the diversification of the genus Picea. Tree Genet Genomes. 7(3):469–484. doi:10.1007/s11295-010-0349-z. Bousquet J, Gérardi S, de Lafontaine G, Jaramillo-Correa JP, Pavy N, Prunier J, Lenz P, Beaulieu J. 2021. Spruce population genomics. In: Raja OP, editor. Population genomics. Cham: Springer International Publishing. 64p. Brandt JP, Flannigan MD, Maynard DG, Thompson ID, Volney WJA. 2013. An introduction to Canada’s boreal zone: ecosystem processes, health, sustainability, and environmental issues. Environ Rev. 21(4):207–226. doi:10.1139/er-2013-0040. Brůna T, Hoff KJ, Lomsadze A, Stanke M, Borodovsky M. 2021. BRAKER2: automatic eukaryotic genome annotation with GeneMark-EP+ and AUGUSTUS supported by a protein database. NAR Genom Bioinform. 3(1):lqaa108. doi:10.1093/nargab/lqaa108. Bürstenbinder K, Mitra D, Quegwer J. 2017. Functions of IQD proteins as hubs in cellular calcium and auxin signaling: a toolbox for shape formation and tissue-specification in plants? Plant Signal Behav. 12(6):e1331198. doi:10.1080/15592324.2017.1331198. Burton P, Bergeron Y, Bogdanski B, Juday G, Kuuluvainen T, Mcafee B, Ogden A, Teplyakov V, Alfaro R, Francis DA, et al. 2010. Sustainability of boreal forests and forestry in a changing environment. In: Mery G, Katila P, Galloway G, Alfaro RI, Kanninen M, Lobovikov M, Varjo J, editors. Forests and society - responding to global drivers of change. Vienna (Austria): International Union of Forest Research Organizations. p. 249–282. 44 Buschiazzo E, Ritland C, Bohlmann J, Ritland K. 2012. Slow but not low: genomic comparisons reveal slower evolutionary rate and higher dN/dS in conifers compared to angiosperms. BMC Evol Biol. 12(1):8. doi:10.1186/1471-2148-12-8. Cabezas JA, González-Martínez SC, Collada C, Guevara MA, Boury C, de María N, Eveno E, Aranda I, Garnier-Géré PH, Brach J, et al. 2015. Nucleotide polymorphisms in a pine ortholog of the Arabidopsis degrading enzyme cellulase KORRIGAN are associated with early growth performance in Pinus pinaster. Tree Physiol. 35(9):1000–1006. doi:10.1093/treephys/tpv050. Celedon JM, Bohlmann J. 2019. Oleoresin defenses in conifers: chemical diversity, terpene synthases and limitations of oleoresin defense under climate change. New Phytol. 224(4):1444–1463. doi:10.1111/nph.15984. Chen J, Källman T, Gyllenstrand N, Lascoux M. 2010. New insights on the speciation history and nucleotide diversity of three boreal spruce species and a tertiary relict. Heredity (Edinb). 104(1):3–14. doi:10.1038/hdy.2009.88. Chen L, Peng Y, Tian J, Wang X, Kong Z, Mao T, Yuan M, Li Y. 2016. TCS1, a microtubule-binding protein, interacts with KCBP/ZWICHEL to regulate trichome cell shape in Arabidopsis thaliana. PLoS Genet. 12(10):e1006266. doi:10.1371/journal.pgen.1006266. Chen N. 2004. Using RepeatMasker to identify repetitive elements in genomic sequences. Curr Protoc Bioinformatics. Chapter 4:Unit 4.10. doi:10.1002/0471250953.bi0410s05. Chen S, Zhou Y, Chen Y, Gu J. 2018. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinform. 34(17):i884–i890. doi:10.1093/bioinformatics/bty560. Cheng C-Y, Krishnakumar V, Chan AP, Thibaud-Nissen F, Schobel S, Town CD. 2017. Araport11: a complete reannotation of the Arabidopsis thaliana reference genome. Plant J. 89(4):789–804. doi:10.1111/tpj.13415. Chiu CC, Bohlmann J. 2022. Mountain pine beetle epidemic: an interplay of terpenoids in host defense and insect pheromones. Annu Rev Plant Biol. 73:475–494. doi:10.1146/annurev-arplant-070921-103617. Chu J, Sadeghi S, Raymond A, Jackman SD, Nip KM, Mar R, Mohamadi H, Butterfield YS, Robertson AG, Birol I. 2014. BioBloom tools: fast, accurate and memory-efficient host species sequence screening using bloom filters. Bioinform. 30(23):3402–3404. doi:10.1093/bioinformatics/btu558. Coombe L, Nikolić V, Chu J, Birol I, Warren RL. 2020. ntJoin: Fast and lightweight assembly-guided scaffolding using minimizer graphs. Bioinform. 36(12):3885–3887. doi:10.1093/bioinformatics/btaa253. Coombe L, Warren RL, Jackman SD, Yang C, Vandervalk BP, Moore RA, Pleasance S, Coope RJ, Bohlmann J, Holt RA, et al. 2016. Assembly of the Complete Sitka spruce chloroplast 45 genome using 10X Genomics’ GemCode sequencing data. PLoS One. 11(9):e0163059. doi:10.1371/journal.pone.0163059. Cubas P, Lauter N, Doebley J, Coen E. 1999. The TCP domain: a motif found in proteins regulating plant growth and development. Plant J. 18(2):215–222. doi:10.1046/j.1365-313x.1999.00444.x. De La Torre AR, Birol I, Bousquet J, Ingvarsson PK, Jansson S, Jones SJM, Keeling CI, MacKay J, Nilsson O, Ritland K, et al. 2014. Insights into Conifer Giga-Genomes. Plant Physiol. 166(4):1724–1732. doi:10.1104/pp.114.248708. De La Torre AR, Lin Y-C, Van de Peer Y, Ingvarsson PK. 2015. Genome-wide analysis reveals diverged patterns of codon bias, gene expression, and rates of sequence evolution in Picea gene families. Genome Biol Evol. 7(4):1002–1015. doi:10.1093/gbe/evv044. Deluca TH, Boisvenue C. 2012. Boreal forest soil carbon: distribution, function and modelling. Forestry. 85(2):161–184. doi:10.1093/forestry/cps003. Depardieu C, Gérardi S, Nadeau S, Parent GJ, Mackay J, Lenz P, Lamothe M, Girardin MP, Bousquet J, Isabel N. 2021. Connecting tree-ring phenotypes, genetic associations and transcriptomics to decipher the genomic architecture of drought adaptation in a widespread conifer. Mol Ecol. 30(16):3898–3917. doi:10.1111/mec.15846. de Vernal A, Hillaire-Marcel C. 2008. Natural variability of Greenland climate, vegetation, and ice volume during the past million years. Science. 320(5883):1622–1625. doi:10.1126/science.1153929. Di Franco A, Poujol R, Baurain D, Philippe H. 2019. Evaluating the usefulness of alignment filtering methods to reduce the impact of errors on evolutionary inferences. BMC Evol Biol. 19(1):21. doi:10.1186/s12862-019-1350-2. Dyderski MK, Paź S, Frelich LE, Jagodziński AM. 2018. How much does climate change threaten European forest tree species distributions? Glob Chang Biol. 24(3):1150–1163. doi:10.1111/gcb.13925. Emms DM, Kelly S. 2019. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genom Biol. 20(1):238. doi:10.1186/s13059-019-1832-y. Ence D, Smith KE, Fan S, Gomide Neves L, Paul R, Wegrzyn J, Peter GF, Kirst M, Brawner J, Nelson CD, et al. 2022. NLR diversity and candidate fusiform rust resistance genes in loblolly pine. G3 (Bethesda). 12(2):jkab421. doi:10.1093/g3journal/jkab421. Flynn JM, Hubley R, Goubert C, Rosen J, Clark AG, Feschotte C, Smit AF. 2020. RepeatModeler2 for automated genomic discovery of transposable element families. Proc Natl Acad Sci USA. 117(17):9451–9457. doi:10.1073/pnas.1921046117. 46 Gagalova KK, Warren RL, Coombe L, Wong J, Nip KM, Yuen MMS, Whitehill JGA, Celedon JM, Ritland C, Taylor GA, et al. 2022. Spruce giga-genomes: structurally similar yet distinctive with differentially expanding gene families and rapidly evolving genes. Plant J. 111(5):1469–1485. doi:10.1111/tpj.15889. Gallemí M, Galstyan A, Paulišić S, Then C, Ferrández-Ayela A, Lorenzo-Orts L, Roig-Villanova I, Wang X, Micol JL, Ponce MR, et al. 2016. DRACULA2 is a dynamic nucleoporin with a role in regulating the shade avoidance syndrome in Arabidopsis. Development. 143(9):1623–1631. doi:10.1242/dev.130211. Gauthier S, Bernier P, Kuuluvainen T, Shvidenko AZ, Schepaschenko DG. 2015. Boreal forest health and global change. Science. 349(6250):819–822. doi:10.1126/science.aaa9092. Geib SM, Hall B, Derego T, Bremer FT, Cannoles K, Sim SB. 2018. Genome Annotation Generator: a simple tool for generating and correcting WGS annotation tables for NCBI submission. GigaScience. 7(4):giy018. doi:10.1093/gigascience/giy018. Gérardi S, Jaramillo-Correa JP, Beaulieu J, Bousquet J. 2010. From glacial refugia to modern populations: new assemblages of organelle genomes generated by differential cytoplasmic gene flow in transcontinental black spruce. Mol Ecol. 19(23):5265–5280. doi:10.1111/j.1365-294X.2010.04881.x. Gonzalez-Ibeas D, Martinez-Garcia PJ, Famula RA, Delfino-Mix A, Stevens KA, Loopstra CA, Langley CH, Neale DB, Wegrzyn JL. 2016. Assessing the gene Content of the megagenome: sugar pine (Pinus lambertiana). G3 (Bethesda). 6(12):3787–3802. doi:10.1534/g3.116.032805. Gurevich A, Saveliev V, Vyahhi N, Tesler G. 2013. QUAST: quality assessment tool for genome assemblies. Bioinform. 29(8):1072–1075. doi:10.1093/bioinformatics/btt086. Hart AJ, Ginzburg S, Xu M, Fisher CR, Rahmatpour N, Mitton JB, Paul R, Wegrzyn JL. 2020. EnTAP: bringing faster and smarter functional annotation to non-model eukaryotic transcriptomes. Mol Ecol Resour. 20(2):591–604. doi:10.1111/1755-0998.13106. Hoffmann AA, Hercus MJ. 2000. Environmental stress as an evolutionary force. BioScience. 50(3):217–226. doi:10.1641/0006-3568(2000)050[0217:ESAAEF]2.3.CO;2. Hornoy B, Pavy N, Gérardi S, Beaulieu J, Bousquet J. 2015. Genetic adaptation to climate in white spruce involves small to moderate allele frequency shifts in functionally diverse genes. Genome Biol Evol. 7(12):3269–3285. doi:10.1093/gbe/evv218. Huang K-C, Lin W-C, Cheng W-H. 2018. Salt hypersensitive mutant 9, a nucleolar APUM23 protein, is essential for salt sensitivity in association with the ABA signaling pathway in Arabidopsis. BMC Plant Biol. 18:40. doi:10.1186/s12870-018-1255-z. Jackman SD, Coombe L, Chu J, Warren RL, Vandervalk BP, Yeo S, Xue Z, Mohamadi H, Bohlmann J, Jones SJM, et al. 2018. Tigmint: correcting assembly errors using linked reads from large molecules. BMC Bioinform. 19(1):393. doi:10.1186/s12859-018-2425-6. 47 Jackman SD, Coombe L, Warren RL, Kirk H, Trinh E, MacLeod T, Pleasance S, Pandoh P, Zhao Y, Coope RJ, et al. 2020. Complete mitochondrial genome of a gymnosperm, Sitka spruce (Picea sitchensis), indicates a complex physical structure. Genome Biol Evol. 12(7):1174–1179. doi:10.1093/gbe/evaa108. Jackman SD, Vandervalk BP, Mohamadi H, Chu J, Yeo S, Hammond SA, Jahesh G, Khan H, Coombe L, Warren RL, et al. 2017. ABySS 2.0: resource-efficient assembly of large genomes using a Bloom filter. Genome Res. 27(5):768–777. doi:10.1101/gr.214346.116. Jackman SD, Warren RL, Gibb EA, Vandervalk BP, Mohamadi H, Chu J, Raymond A, Pleasance S, Coope R, Wildung MR, et al. 2015. Organellar genomes of white spruce (Picea glauca): assembly and annotation. Genome Biol Evol. 8(1):29–41. doi:10.1093/gbe/evv244. Jeong DI, Sushama L, Naveed Khaliq M. 2014. The role of temperature in drought projections over North America. Clim Change. 127(2):289–303. doi:10.1007/s10584-014-1248-3. Kaplan-Levy RN, Brewer PB, Quon T, Smyth DR. 2012. The trihelix family of transcription factors – light, stress and development. Trends Plant Sci. 17(3):163–171. doi:10.1016/j.tplants.2011.12.002. Katz LS, Griswold T, Morrison SS, Caravas JA, Zhang S, den Bakker HC, Deng X, Carleton HA. 2019. Mashtree: a rapid comparison of whole genome sequence files. J Open Source Softw. 4(44):1762. doi:10.21105/joss.01762. Kawasaki T, Yamada K, Yoshimura S, Yamaguchi K. 2017. Chitin receptor-mediated activation of MAP kinases and ROS production in rice and Arabidopsis. Plant Signal Behav. 12(9):e1361076. doi:10.1080/15592324.2017.1361076. Kayes I, Mallik A. 2020. Boreal forests: distributions, biodiversity, and management. In: Leal Filho W, Azul AM, Brandli L, Lange Salvia A, Wall T, editors. Life on land. Cham: Springer International Publishing. p. 1–12. https://doi.org/10.1007/978-3-319-71065-5_17-1. Kharenko OA, Boyd J, Nelson KM, Abrams SR, Loewen MC. 2011. Identification and characterization of interactions between abscisic acid and mitochondrial adenine nucleotide translocators. Biochem J. 437(1):117–123. doi:10.1042/BJ20101898. Kim D, Langmead B, Salzberg SL. 2015. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 12(4):357–360. doi:10.1038/nmeth.3317. Kimura M. 1977. Preponderance of synonymous changes as evidence for the neutral theory of molecular evolution. Nature. 267(5608):275–276. doi:10.1038/267275a0. Kobayashi K, Fujii S, Sasaki D, Baba S, Ohta H, Masuda T, Wada H. 2014. Transcriptional regulation of thylakoid galactolipid biosynthesis coordinated with chlorophyll biosynthesis during the development of chloroplasts in Arabidopsis. Front Plant Sci. 5:272. doi:10.3389/fpls.2014.00272. 48 Kriventseva EV, Kuznetsov D, Tegenfeldt F, Manni M, Dias R, Simão FA, Zdobnov EM. 2019. OrthoDB v10: sampling the diversity of animal, plant, fungal, protist, bacterial and viral genomes for evolutionary and functional annotations of orthologs. Nucleic Acids Res. 47(D1):D807–D811. doi:10.1093/nar/gky1053. Kryazhimskiy S, Plotkin JB. 2008. The population genetics of dN/dS. PLoS Genet. 4(12):e1000304. doi:10.1371/journal.pgen.1000304. Kuromori T, Wada T, Kamiya A, Yuguchi M, Yokouchi T, Imura Y, Takabe H, Sakurai T, Akiyama K, Hirayama T, et al. 2006. A trial of phenome analysis using 4000 Ds-insertional mutants in gene-coding regions of Arabidopsis. Plant J. 47(4):640–651. doi:10.1111/j.1365-313X.2006.02808.x. Lannoo N, Van Damme EJM. 2014. Lectin domains at the frontiers of plant defense. Front Plant Sci. 5:397. doi:10.3389/fpls.2014.00397. Laverdière J-P, Lenz P, Nadeau S, Depardieu C, Isabel N, Perron M, Beaulieu J, Bousquet J. 2022. Breeding for adaptation to climate change: genomic selection for drought response in a white spruce multi-site polycross test. Evol Appl. 15(3):383–402. doi:10.1111/eva.13348. Li C-H, Wang G, Zhao J-L, Zhang L-Q, Ai L-F, Han Y-F, Sun D-Y, Zhang S-W, Sun Y. 2014. The receptor-like kinase SIT1 ,ediates salt sensitivity by activating MAPK3/6 and regulating ethylene homeostasis in rice. Plant Cell. 26(6):2538–2553. doi:10.1105/tpc.114.125187. Li H. 2014. lh3/trimadap. [accessed 2022 Aug 10]. https://github.com/lh3/trimadap. Lin D, Coombe L, Jackman SD, Gagalova KK, Warren RL, Hammond SA, Kirk H, Pandoh P, Zhao Y, Moore RA, et al. 2019a. Complete chloroplast genome sequence of a white spruce (Picea glauca, genotype WS77111) from eastern Canada. Microbiol Resour Announc. 8(23):e00381-19. doi:10.1128/MRA.00381-19. Lin D, Coombe L, Jackman SD, Gagalova KK, Warren RL, Hammond SA, McDonald H, Kirk H, Pandoh P, Zhao Y, et al. 2019b. Complete chloroplast genome sequence of an Engelmann spruce (Picea engelmannii, genotype Se404-851) from western Canada. Microbiol Resour Announc. 8(24):e00382-19. doi:10.1128/MRA.00382-19. Lin J-S, Huang X-X, Li Q, Cao Y, Bao Y, Meng X-F, Li Y-J, Fu C, Hou B-K. 2016. UDP-glycosyltransferase 72B1 catalyzes the glucose conjugation of monolignols and is essential for the normal cell wall lignification in Arabidopsis thaliana. Plant J. 88(1):26–42. doi:10.1111/tpj.13229. Lo T, Coombe L, Lin D, Warren RL, Kirk H, Pandoh P, Zhao Y, Moore RA, Mungall AJ, Ritland C, et al. 2020. Complete chloroplast genome sequence of a black spruce (Picea mariana) from eastern Canada. Microbiol Resour Announc. 9(39):e00877-20. doi:10.1128/MRA.00877-20. 49 Lockwood JD, Aleksić JM, Zou J, Wang J, Liu J, Renner SS. 2013. A new phylogeny for the genus Picea from plastid, mitochondrial, and nuclear sequences. Mol Phylogenet Evol. 69(3):717–727. doi:10.1016/j.ympev.2013.07.004. Lowe R, Shirley N, Bleackley M, Dolan S, Shafee T. 2017. Transcriptomics technologies. PLoS Comput Biol. 13(5):e1005457. doi:10.1371/journal.pcbi.1005457. Mann IK, Wegrzyn JL, Rajora OP. 2013. Generation, functional annotation and comparative analysis of black spruce (Picea mariana) ESTs: an important conifer genomic resource. BMC Genom. 14:702. doi:10.1186/1471-2164-14-702. Marino G, Funk C. 2012. Matrix metalloproteinases in plants: a brief overview. Physiol Plant. 145(1):196–202. doi:10.1111/j.1399-3054.2011.01544.x. Mathur S, Vyas S, Kapoor S, Tyagi AK. 2011. The mediator complex in plants: structure, phylogeny, and expression profiling of representative genes in a dicot (Arabidopsis) and a monocot (rice) during reproduction and abiotic Stress. Plant Physiol. 157(4):1609–1627. doi:10.1104/pp.111.188300. Michaels SD, Bezerra IC, Amasino RM. 2004. FRIGIDA-related genes are required for the winter-annual habit in Arabidopsis. Proc Natl Acad Sci. 101(9):3281–3285. doi:10.1073/pnas.0306778101. Mølhøj M, Jørgensen B, Ulvskov P, Borkhardt B. 2001. Two Arabidopsis thaliana genes, KOR2 and KOR3, which encode membrane-anchored endo-1,4-beta-D-glucanases, are differentially expressed in developing leaf trichomes and their support cells. Plant Mol Biol. 46(3):263–275. doi:10.1023/a:1010688726755. Mullin TJ, Andersson B, Bastien J-C, Beaulieu J, Burdon RD, Dvorak WS, King JN, Kondo T, Krakowski J, Lee SJ, et al. 2011. Economic importance, breeding objectives and achievements. In: Plomion C, Bousquet J, Kole C, editors. Genetics, genomics, and breeding of conifers. New York (NY): CRC Press. 456 p. Napier JD, de Lafontaine G, Hu FS. 2020. Exploring genomic variation associated with drought stress in Picea mariana populations. Ecol Evol. 10(17):9271–9282. doi:10.1002/ece3.6614. Neale DB, Zimin AV, Zaman S, Scott AD, Shrestha B, Workman RE, Puiu D, Allen BJ, Moore ZJ, Sekhwal MK, et al. 2022. Assembled and annotated 26.5 Gbp coast redwood genome: a resource for estimating evolutionary adaptive potential and investigating hexaploid origin. G3 (Bethesda). 12(1):jkab380. doi:10.1093/g3journal/jkab380. Nicklen EF, Roland CA, Ruess RW, Scharnweber T, Wilmking M. 2021. Divergent responses to permafrost and precipitation reveal mechanisms for the spatial variation of two sympatric spruce. Ecosphere. 12(7):e03622. doi:10.1002/ecs2.3622. 50 Niu S, Li J, Bo W, Yang W, Zuccolo A, Giacomello S, Chen X, Han F, Yang J, Song Y, et al. 2022. The Chinese pine genome and methylome unveil key features of conifer evolution. Cell. 185(1):204-217.e14. doi:10.1016/j.cell.2021.12.006. Nystedt B, Street NR, Wetterbom A, Zuccolo A, Lin Y-C, Scofield DG, Vezzi F, Delhomme N, Giacomello S, Alexeyenko A, et al. 2013. The Norway spruce genome sequence and conifer genome evolution. Nature. 497(7451):579–584. doi:10.1038/nature12211. Okushima Y, Fukaki H, Onoda M, Theologis A, Tasaka M. 2007. ARF7 and ARF19 regulate lateral root formation via direct activation of LBD/ASL genes in Arabidopsis. Plant Cell. 19(1):118–130. doi:10.1105/tpc.106.047761. Okushima Y, Overvoorde PJ, Arima K, Alonso JM, Chan A, Chang C, Ecker JR, Hughes B, Lui A, Nguyen D, et al. 2005. Functional genomic analysis of the auxin response factor gene family members in Arabidopsis thaliana: unique and overlapping functions of ARF7 and ARF19. Plant Cell. 17(2):444–463. doi:10.1105/tpc.104.028316. O’Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, Rajput B, Robbertse B, Smith-White B, Ako-Adjei D, et al. 2016. Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 44(D1):D733–D745. doi:10.1093/nar/gkv1189. Otsuka K, Mamiya A, Konishi M, Nozaki M, Kinoshita A, Tamaki H, Arita M, Saito M, Yamamoto K, Hachiya T, et al. 2021. Temperature-dependent fasciation mutants provide a link between mitochondrial RNA processing and lateral root morphogenesis. Elife. 10:e61611. doi:10.7554/eLife.61611. Ou S, Jiang N. 2018. LTR_retriever: a highly accurate and sensitive program for identification of long terminal repeat retrotransposons. Plant Physiol. 176(2):1410–1422. doi:10.1104/pp.17.01310. Pan Y, Birdsey RA, Fang J, Houghton R, Kauppi PE, Kurz WA, Phillips OL, Shvidenko A, Lewis SL, Canadell JG, et al. 2011. A large and persistent carbon sink in the world’s forests. Science. 333(6045):988–993. Parla JS, Iossifov I, Grabill I, Spector MS, Kramer M, McCombie WR. 2011. A comparative analysis of exome capture. Genome Biol. 12(9):R97. doi:10.1186/gb-2011-12-9-r97. Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. 2017. Salmon: fast and bias-aware quantification of transcript expression using dual-phase inference. Nat Methods. 14(4):417–419. doi:10.1038/nmeth.4197. Paulino D, Warren RL, Vandervalk BP, Raymond A, Jackman SD, Birol I. 2015. Sealer: a scalable gap-closing application for finishing draft genomes. BMC Bioinform. 16(1):230. doi:10.1186/s12859-015-0663-4. 51 Pavy N, Gagnon F, Deschênes A, Boyle B, Beaulieu J, Bousquet J. 2016. Development of highly reliable in silico SNP resource and genotyping assay from exome capture and sequencing: an example from black spruce (Picea mariana). Mol Ecol Resour. 16(2):588–598. doi:10.1111/1755-0998.12468. Pinno BD, Errington RC, Thompson DK. 2013. Young jack pine and high severity fire combine to create potentially expansive areas of understocked forest. For Ecol Manag. 310:517–522. doi:10.1016/j.foreco.2013.08.055. Proost S, Van Bel M, Sterck L, Billiau K, Van Parys T, Van de Peer Y, Vandepoele K. 2009. PLAZA: a comparative genomics resource to study gene and genome evolution in plants. Plant Cell. 21(12):3718–3731. doi:10.1105/tpc.109.071506. Prunier J, Laroche J, Beaulieu J, Bousquet J. 2011. Scanning the genome for gene SNPs related to climate adaptation and estimating selection at the molecular level in boreal black spruce. Mol Ecol. 20(8):1702–1716. doi:10.1111/j.1365-294X.2011.05045.x. Rajora OP, Mann IK. 2021. Development and characterization of novel EST-based single-copy genic microsatellite DNA markers in white spruce and black spruce. Mol Biol Rep. 48(3):2963–2971. doi:10.1007/s11033-021-06231-1. Ralph SG, Chun HJE, Kolosova N, Cooper D, Oddy C, Ritland CE, Kirkpatrick R, Moore R, Barber S, Holt RA, et al. 2008. A conifer genomics resource of 200,000 spruce (Picea spp.) ESTs and 6,464 high-quality, sequence-finished full-length cDNAs for Sitka spruce (Picea sitchensis). BMC Genom. 9:484. doi:10.1186/1471-2164-9-484. Ran J-H, Shen T-T, Liu W-J, Wang P-P, Wang X-Q. 2015. Mitochondrial introgression and complex biogeographic history of the genus Picea. Mol Phylogenet Evol. 93:63–76. doi:10.1016/j.ympev.2015.07.020. Ran J-H, Wei X-X, Wang X-Q. 2006. Molecular phylogeny and biogeography of Picea (Pinaceae): Implications for phylogeographical studies using cytoplasmic haplotypes. Mol Phylogenet Evol. 41(2):405–419. doi:10.1016/j.ympev.2006.05.039. Ranwez V, Douzery EJP, Cambon C, Chantret N, Delsuc F. 2018. MACSE v2: toolkit for the alignment of coding sequences accounting for frameshifts and stop codons. Mol Biol Evol. 35(10):2582–2584. doi:10.1093/molbev/msy159. Rao V, Virupapuram V. 2021. Arabidopsis F-box protein At1g08710 interacts with transcriptional protein ADA2b and imparts drought stress tolerance by negatively regulating seedling growth. Biochem Biophys Res Commun. 536:45–51. doi:10.1016/j.bbrc.2020.12.054. Riechmann JL, Heard J, Martin G, Reuber L, Jiang C, Keddie J, Adam L, Pineda O, Ratcliffe OJ, Samaha RR, et al. 2000. Arabidopsis transcription factors: genome-wide comparative analysis among eukaryotes. Science. 290(5499):2105–2110. doi:10.1126/science.290.5499.2105. 52 Rigault P, Boyle B, Lepage P, Cooke JEK, Bousquet J, MacKay JJ. 2011. A white spruce gene catalog for conifer genome analyses. Plant Physiol. 157(1):14–28. doi:10.1104/pp.111.179663. Rossi S, Bousquet J. 2014. The bud break process and its variation among local populations of boreal black spruce. Front Plant Sci. 5:574. doi:10.3389/fpls.2014.00574. Royer-Tardif S, Boisvert-Marsh L, Godbout J, Isabel N, Aubin I. 2021. Finding common ground: toward comparable indicators of adaptive capacity of tree species to a changing climate. Ecol and Evol. 11(19):13081–13100. doi:10.1002/ece3.8024. Rudd S. 2003. Expressed sequence tags: alternative or complement to whole genome sequences? Trends Plant Sci. 8(7):321–329. doi:10.1016/S1360-1385(03)00131-6. Schneider A, Souvorov A, Sabath N, Landan G, Gonnet GH, Graur D. 2009. Estimates of positive Darwinian selection are inflated by errors in sequencing, annotation, and alignment. Genome Biol Evol. 1:114–118. doi:10.1093/gbe/evp012. Scott AD, Zimin AV, Puiu D, Workman R, Britton M, Zaman S, Caballero M, Read AC, Bogdanove AJ, Burns E, et al. 2020. A reference genome sequence for giant sequoia. G3 (Bethesda). 10(11):3907–3919. doi:10.1534/g3.120.401612. Shao C-C, Shen T-T, Jin W-T, Mao H-J, Ran J-H, Wang X-Q. 2019. Phylotranscriptomics resolves interspecific relationships and indicates multiple historical out-of-North America dispersals through the Bering Land Bridge for the genus Picea (Pinaceae). Mol Phylogenet Evol. 141:106610. doi:10.1016/j.ympev.2019.106610. Shen G, Kuppu S, Venkataramani S, Wang J, Yan J, Qiu X, Zhang H. 2010. Ankyrin repeat-containing protein 2A is an essential molecular chaperone for peroxisomal membrane-bound ascorbate peroxidase3 in Arabidopsis. Plant Cell. 22(3):811–831. doi:10.1105/tpc.109.065979. Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. 2015. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 31(19):3210–3212. doi:10.1093/bioinformatics/btv351. Singh P, Zimmerli L. 2013. Lectin receptor kinases in plant innate immunity. Front Plant Sci. 4:124. doi:10.3389/fpls.2013.00124. Snider J, Thibault G, Houry WA. 2008. The AAA+ superfamily of functionally diverse proteins. Genome Biol. 9(4):216. doi:10.1186/gb-2008-9-4-216. Sniderhan AE, McNickle GG, Baltzer JL. 2018. Assessing local adaptation vs. plasticity under different resource conditions in seedlings of a dominant boreal tree species. AoB Plants. 10(1):ply004. doi:10.1093/aobpla/ply004. Stanke M, Keller O, Gunduz I, Hayes A, Waack S, Morgenstern B. 2006. AUGUSTUS: ab initio prediction of alternative transcripts. Nucleic Acids Res. 34(suppl_2):W435–W439. doi:10.1093/nar/gkl200. 53 Steinegger M, Söding J. 2017. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat Biotechnol. 35(11):1026–1028. doi:10.1038/nbt.3988. Stival Sena J, Giguère I, Boyle B, Rigault P, Birol I, Zuccolo A, Ritland K, Ritland C, Bohlmann J, Jones S, et al. 2014. Evolution of gene structure in the conifer Picea glauca: a comparative analysis of the impact of intron size. BMC Plant Biol. 14(1):95. doi:10.1186/1471-2229-14-95. Sullivan AR, Eldfjell Y, Schiffthaler B, Delhomme N, Asp T, Hebelstrup KH, Keech O, Öberg L, Møller IM, Arvestad L, et al. 2020. The mitogenome of Norway spruce and a reappraisal of mitochondrial recombination in plants. Genome Biol Evol. 12(1):3586–3598. doi:10.1093/gbe/evz263. Sullivan AR, Schiffthaler B, Thompson SL, Street NR, Wang X-R. 2017. Interspecific plastome recombination reflects ancient reticulate evolution in Picea (Pinaceae). Mol Biol Evol. 34(7):1689–1701. doi:10.1093/molbev/msx111. Suzek BE, Wang Y, Huang H, McGarvey PB, Wu CH, UniProt Consortium. 2015. UniRef clusters: a comprehensive and scalable alternative for improving sequence similarity searches. Bioinformatics. 31(6):926–932. doi:10.1093/bioinformatics/btu739. Tabas-Madrid D, Méndez-Vigo B, Arteaga N, Marcer A, Pascual-Montano A, Weigel D, Xavier Picó F, Alonso-Blanco C. 2018. Genome-wide signatures of flowering adaptation to climate temperature: regional analyses in a highly diverse native range of Arabidopsis thaliana. Plant Cell Environ. 41(8):1806–1820. doi:10.1111/pce.13189. Tamura K, Stecher G, Kumar S. 2021. MEGA11: molecular evolutionary genetics analysis version 11. Mol Biol Evol. 38(7):3022–3027. doi:10.1093/molbev/msab120. Tarnocai C, Canadell JG, Schuur EAG, Kuhry P, Mazhitova G, Zimov S. 2009. Soil organic carbon pools in the northern circumpolar permafrost region. Global Biogeochem Cycles. 23(2). doi:10.1029/2008GB003327. Tenhaken R, Doerks T, Bork P. 2005. DCD – a novel plant specific domain in proteins involved in development and programmed cell death. BMC Bioinform. 6(1):169. doi:10.1186/1471-2105-6-169. The UniProt Consortium. 2019. UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 47(D1):D506–D515. doi:10.1093/nar/gky1049. Thomson AM, Riddell CL, Parker WH. 2009. Boreal forest provenance tests used to predict optimal growth and response to climate change: 2. Black spruce. Can J For Res. 39(1):143–153. doi:10.1139/X08-167. Trujillo-Moya C, George J-P, Fluch S, Geburek T, Grabner M, Karanitsch-Ackerl S, Konrad H, Mayer K, Sehr EM, Wischnitzki E, et al. 2018. Drought sensitivity of Norway spruce at the species’ warmest fringe: quantitative and molecular analysis reveals high genetic variation among and within provenances. G3 (Bethesda). 8(4):1225–1245. doi:10.1534/g3.117.300524. 54 Tunc-Ozdemir M, Rato C, Brown E, Rogers S, Mooneyham A, Frietsch S, Myers CT, Poulsen LR, Malhó R, Harper JF. 2013. Cyclic nucleotide gated channels 7 and 8 are essential for male reproductive fertility. PLoS One. 8(2):e55277. doi:10.1371/journal.pone.0055277. Van Ghelder C, Parent GJ, Rigault P, Prunier J, Giguère I, Caron S, Stival Sena J, Deslauriers A, Bousquet J, Esmenjaud D, et al. 2019. The large repertoire of conifer NLR resistance genes includes drought responsive and highly diversified RNLs. Sci Rep. 9(1):11614. doi:10.1038/s41598-019-47950-7. Vandervalk BP, Yang C, Xue Z, Raghavan K, Chu J, Mohamadi H, Jackman SD, Chiu R, Warren RL, Birol I. 2015. Konnector v2.0: pseudo-long reads from paired-end sequencing data. BMC Med Genom. 8(3):S1. doi:10.1186/1755-8794-8-S3-S1. Västermark Å, Almén MS, Simmen MW, Fredriksson R, Schiöth HB. 2011. Functional specialization in nucleotide sugar transporters occurred through differentiation of the gene cluster EamA (DUF6) before the radiation of Viridiplantae. BMC Evol Biol. 11:123. doi:10.1186/1471-2148-11-123. Venkatraman M, Fleischer RC, Tsuchiya MTN. 2021. Comparative analysis of annotation pipelines using the first Japanese white-eye (Zosterops japonicus) genome. Genome Biol Evol. 13(5):evab063. doi:10.1093/gbe/evab063. Vidali L, van Gisbergen PAC, Guérin C, Franco P, Li M, Burkart GM, Augustine RC, Blanchoin L, Bezanilla M. 2009. Rapid formin-mediated actin-filament elongation is essential for polarized plant cell growth. Proc Natl Acad Sci. 106(32):13341–13346. doi:10.1073/pnas.0901170106. Vuruputoor VS, Monyak D, Fetter KC, Webster C, Bhattarai A, Shrestha B, Zaman S, Bennett J, McEvoy SL, Caballero M, et al. 2023. Welcome to the big leaves: best practices for improving genome annotation in non‐model plant genomes. Appl Plant Sci. 11(4):e11533. doi:10.1002/aps3.11533. Wan J, Patel A, Mathieu M, Kim S-Y, Xu D, Stacey G. 2008. A lectin receptor-like kinase is required for pollen development in Arabidopsis. Plant Mol Biol. 67(5):469–482. doi:10.1007/s11103-008-9332-6. Wang K, Guo Q, Froehlich JE, Hersh HL, Zienkiewicz A, Howe GA, Benning C. 2018. Two abscisic acid-responsive plastid lipase genes involved in jasmonic acid biosynthesis in Arabidopsis thaliana. Plant Cell. 30(5):1006–1022. doi:10.1105/tpc.18.00250. Wang Y, Liu K, Liao H, Zhuang C, Ma H, Yan X. 2008. The plant WNK gene family and regulation of flowering time in Arabidopsis. Plant Biol (Stuttg). 10(5):548–562. doi:10.1111/j.1438-8677.2008.00072.x. Warren RL, Keeling CI, Yuen MMS, Raymond A, Taylor GA, Vandervalk BP, Mohamadi H, Paulino D, Chiu R, Jackman SD, et al. 2015. Improved white spruce (Picea glauca) genome 55 assemblies and annotation of large gene families of conifer terpenoid and phenolic defense metabolism. Plant J. 83(2):189–212. doi:10.1111/tpj.12886. Warren RL, Yang C, Vandervalk BP, Behsaz B, Lagman A, Jones SJM, Birol I. 2015. LINKS: scalable, alignment-free scaffolding of draft genomes with long reads. GigaScience. 4(1):35. doi:10.1186/s13742-015-0076-3. Weng C, Jackson ST. 2000. Species differentiation of North American spruce (Picea) based on morphological and anatomical characteristics of needles. Can J Bot. 78(11):1367–1383. doi:10.1139/b00-111. Whelan S, Irisarri I, Burki F. 2018. PREQUAL: detecting non-homologous characters in sets of unaligned homologous sequences. Bioinformatics. 34(22):3929–3930. doi:10.1093/bioinformatics/bty448. Wright JW. 1955. Species crossability in spruce in relation to distribution and taxonomy. For Sci. 1(4):319–49. Wu Y, Xun Q, Guo Y, Zhang J, Cheng K, Shi T, He K, Hou S, Gou X, Li J. 2016. Genome-wide expression pattern analyses of the Arabidopsis leucine-rich repeat receptor-like kinases. Molecular Plant. 9(2):289–300. doi:10.1016/j.molp.2015.12.011. Xue B, Charest PJ, Devantier Y, Rutledge RG. 2003. Characterization of a MYBR2R3 gene from black spruce (Picea mariana) that shares functional conservation with maize C1. Mol Gen Genomics. 270(1):78–86. doi:10.1007/s00438-003-0898-z. Yandell M, Ence D. 2012. A beginner’s guide to eukaryotic genome annotation. Nat Rev Genet. 13(5):329–342. doi:10.1038/nrg3174. Yang Z. 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 24(8):1586–1591. doi:10.1093/molbev/msm088. Yang Z, dos Reis M. 2011. Statistical properties of the branch-site test of positive selection. Mol Biol Evol. 28(3):1217–1228. doi:10.1093/molbev/msq303. Yeaman S, Hodgins KA, Lotterhos KE, Suren H, Nadeau S, Degner JC, Nurkowski KA, Smets P, Wang T, Gray LK, et al. 2016. Convergent local adaptation to climate in distantly related conifers. Science. 353(6306):1431–1433. doi:10.1126/science.aaf7812. Yeo S, Coombe L, Warren RL, Chu J, Birol I. 2018. ARCS: scaffolding genome drafts with linked reads. Bioinformatics. 34(5):725–731. doi:10.1093/bioinformatics/btx675. Yu D, Chen C, Chen Z. 2001. Evidence for an important role of WRKY DNA binding proteins in the regulation of NPR1 gene expression. Plant Cell. 13(7):1527–1540. doi:10.1105/tpc.010115. 56 Zhang J, Nielsen R, Yang Z. 2005. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol Biol Evol. 22(12):2472–2479. doi:10.1093/molbev/msi237. Zhang X, Flato G, Kirchmeier-Young M, Vincent L, Wan H, Wang X, Rong R, Fife J, Li G, Karin VV. 2019. Changes in temperature and precipitation across Canada. In: Bush E, Lemmen DS, editors. Canada’s changing climate report. Ottawa (ON): Government of Canada. p. 112-193. Zhao J-L, Zhang L-Q, Liu N, Xu S-L, Yue Z-L, Zhang L-L, Deng Z-P, Burlingame AL, Sun D-Y, Wang Z-Y, et al. 2019. Mutual regulation of receptor-like kinase SIT1 and B’κ-PP2A shapes the early response of rice to salt stress. Plant Cell. 31(9):2131–2151. doi:10.1105/tpc.18.00706. Zheng B, He H, Zheng Y, Wu W, McCormick S. 2014. An ARID domain-containing protein within nuclear bodies is required for sperm cell formation in Arabidopsis thaliana. PLoS Genet. 10(7):e1004421. doi:10.1371/journal.pgen.1004421. Zimin AV, Stevens KA, Crepeau MW, Puiu D, Wegrzyn JL, Yorke JA, Langley CH, Neale DB, Salzberg SL. 2017. An improved assembly of the loblolly pine mega-genome using long-read single-molecule sequencing. Gigascience. 6(1):1–4. doi:10.1093/gigascience/giw016. 57 Appendices Appendix A Genome sequencing, assembly and annotation supplementary material Table A.1 Sequencing data used for assembly of Picea mariana nuclear genome. Platform Read length (bp) Nominal fragment size (bp) Number of libraries Number of sequencing lanes Fold coverage SRA accessions 10x Genomics Chromium 2x150 400 1 5 31.15x SRX890468 Illumina HiSeqX 2x250 400 2 2 6.30x SRX9351292, SRX9351293 Illumina HiSeqX 2x250 800 2 2 8.50x SRX9351294, SRX9351295 58 Table A.2 Summary of protein evidence provided to BRAKER. Data was obtained from OrthoDB v10 Viridiplantae database, UniProtKB/Swissprot plant entries, Picea glauca manual annotations of terpene synthase and cytochrome P450 genes, and common protein found in at least three of four North American spruces included in the analyses (Picea glauca, interior spruce, Picea engelmannii and Picea sitchensis), as identified by reciprocal-best hit analysis. Specifically, reciprocal BLAST searches were performed between Picea engelmannii, the North American spruce with the most annotations, and each of the other three North American spruces. Proteins found in at least two of the three reciprocal BLAST search results were considered common amongst the North American spruces and included as evidence. Database Number of proteins OrthoDB v10 Viridiplantae 3,510,742 UniProtKB/SwissProt plant entries 40,925 P. glauca manual annotations 298 Common proteins in Canadian spruce annotations 29,204 Total number of unique proteins 3,463,432 59 Table A.3 RNA-seq libraries to be used as annotation evidence for BRAKER.SRA accessions Tissue type Collection site Platform Read length (bp) References SRR9595774 seedling Sheffield's Seed Company, USA Illumina HiSeq 2000 2x100 Shao et al. 2019 SRR9595777 seedling Sheffield's Seed Company, USA Illumina HiSeq 2000 2x100 Shao et al. 2019 Table A.4 Statistics of RNA-seq alignments provided to BRAKER as annotation evidence. After obtaining the data from SRA, quality and contaminant filtering was performed prior to aligning the reads to the genome. SRR9595774 SRR9595777 Total number of reads 51,739,556 58,954,918 Number of reads after fastp & contaminant-filtering 43,337,856 53,418,516 Number of aligned reads 38,663,631 47,320,896 Number of multi-mapped reads 5,166,472 5,705,162 Table A.5 Number of protein-coding genes and transcripts after each stage in the Picea mariana annotation pipeline. Stage in annotation pipeline Number of protein-coding genes Number of transcripts BRAKER 146,065 150,036 Remove incomplete (no start and stop codon) and/or intron lengths < 10 bp 121,974 125,377 EnTAP 67,373 69,782 Remove gag and pol Pfam domains 66,332 68,738 Identified as fragmented 883 1,070 60 Appendix B Phylogenetic and comparative genomics analyses supplementary material Table B.1 NCBI accessions of nuclear, chloroplast and mitochondrial genome assemblies used in analyses. Species Nuclear genome Chloroplast genome Mitochondrial genome References Picea mariana JASDQU000000000 MT261462 10.5281/zenodo.7828188 This study, Lo et al. 2020 Interior spruce ALWZ000000000 NC_028594 LKAM01000001.1-LKAM01000036.1 Gagalova et al. 2022, Jackman et al. 2015 Picea abies GCA_900067695 NC_021456 MN642623-MN642626 Nystedt et al. 2013 Picea engelmannii WSFP000000000 NC_041067 10.5281/zenodo.7828263 Gagalova et al. 2022, Lin et al. 2019b Picea glauca JZKD000000000 MK174379 10.5281/zenodo.7828293 Gagalova et al. 2022, Lin et al. 2019a Picea sitchensis SNQJ000000000 KU215903 MK697696-MK607708 Gagalova et al. 2022, Coombe et al. 2016, Jackman et al. 2020 Pinus lambertiana GCA_001447015 NC_011156 https://treegenesdb.org/FTP/Genomes/Pila/mito/plambertiana.mito.scafSeq.gz Gonzalez-Ibeas et al. 2016 Pinus taeda GCA_000404065 KY964286 NC_039746 Zimin et al. 2017, Asaf et al. 2018 61 B.1 Picea mariana specific orthogroup analysis supplementary material Table B.2 Summary of OrthoFinder results. The longest transcript of complete protein-coding sequences with lengths  1 kbp and intron lengths  10 bp were passed as input to OrthoFinder. Species and genome assembly versions Number of protein-coding sequences Number of protein-coding sequences assigned to orthogroups Number of orthogroups Number of species-specific orthogroups Reference Picea mariana v1 35,719 32,936 19,330 560 This study Picea abies v1 23,233 21,494 16,975 139 Nystedt et al. 2013 Interior spruce v5 28,878 16,638 27,309 188 Gagalova et al. 2022 Picea engelmannii v1 35,226 33,104 18,180 253 Gagalova et al. 2022 Picea glauca v2 30,140 28,696 16,852 269 Gagalova et al. 2022 Picea sitchensis v1 30,324 28,381 17,661 197 Gagalova et al. 2022 Pinus lambertiana v1.5 31,547 29,395 18,006 793 Gonzalez-Ibeas et al. 2016 Pinus taeda v2.01 37,499 33,764 17,772 1,164 Zimin et al. 2017 62 Figure B.1 Best matching PFAM domains (n > 5) in Picea mariana specific protein-coding sequences. The annotations within each bar indicates the distribution of EggNOG functional categories for sequences containing the PFAM domain. Functional category annotations were assigned based on hits to entries in the EggNOG database, while PFAM annotations were assigned via hits in EggNOG and/or InterPro databases as part of the EnTAP annotation pipeline. Bars in grey represent sequences that were not assigned functional categories as they were annotated with InterPro hits only. 5V; 1Q7O11C16S12S7S7L; 1K6P3J; 4S6A6Z11T8T3K; 3S7G11K10T; 2S9Q6O,W40T6K2G; 7C18C,G3T,U; 1S; 4U3Beta_HSDK_transKH_1KinesinMed26RNA_pol_Rpb2_2AAAF−boxMFS_1Helicase_CLRR_1PG_binding_1WD40Kelch_1p450RuBisCO_largeTCPAldedhLectin_legBMyb_DNA−bindingEamANB−ARCUDPGTPkinaseDUF42830 20 40CountPFAM domainNo EggNOG functional category assignedEggNOG functional category assigned63 B.2 Positive selection analysis supplementary material Table B.3 Statistics for the 45 positively selected protein-coding sequences in Picea mariana. Log-likelihood values (lnLs) of the null (A1) and alternative (A) branch-site models were determined by CODEML, then used to calculate the log-likelihood ratio test statistic (LRT statistic) and corresonding p-values (P-value) for each sequence. Picea mariana mRNA ID A1 lnL A lnL LRT statistic P-value G9G70_30239.1 -1883.985411 -1881.34589 5.27905 0.021583521 G9G70_36237.1 -8424.182251 -8377.05373 94.25705 2.77E-22 G9G70_62668.1 -3509.896272 -3473.02342 73.745702 8.89E-18 G9G70_4240.1 -4349.999715 -4237.79674 224.405944 9.89E-51 G9G70_60414.1 -2345.748409 -2323.14698 45.202852 1.78E-11 G9G70_45605.1 -13995.41974 -13983.9371 22.965346 1.65E-06 G9G70_2442.1 -6389.133201 -6330.89047 116.485466 3.72E-27 G9G70_26879.1 -9569.712431 -9553.04382 33.33722 7.75E-09 G9G70_2048.1 -2205.524985 -2157.45088 96.148216 1.07E-22 G9G70_725.1 -7382.129136 -7377.89189 8.474494 0.003601606 G9G70_22547.2 -1867.689918 -1865.58135 4.21713 0.040017794 G9G70_1487.1 -1692.799197 -1688.20379 9.190812 0.00243233 G9G70_28286.1 -2655.389319 -2652.97462 4.829404 0.027978301 G9G70_13089.1 -2902.246647 -2873.55782 57.377652 3.60E-14 G9G70_7440.1 -2495.682065 -2492.27264 6.818858 0.009020023 G9G70_53953.1 -2627.791758 -2574.93945 105.704624 8.56E-25 G9G70_2518.1 -3203.725414 -3199.77429 7.902244 0.004937351 G9G70_11476.1 -3466.607601 -3459.78464 13.645928 0.000220719 G9G70_4241.1 -4970.573594 -4943.9602 53.226788 2.97E-13 G9G70_47211.1 -3954.809791 -3951.12711 7.365364 0.006649214 G9G70_24781.1 -4150.063763 -4145.01018 10.107172 0.001476934 G9G70_44353.1 -3103.671687 -3094.5708 18.201774 1.99E-05 G9G70_38193.1 -3123.498902 -3117.31136 12.37509 0.0004351 64 G9G70_27262.1 -12009.81685 -11967.942 83.74976 5.62E-20 G9G70_22542.2 -4835.502147 -4833.11403 4.776244 0.028854992 G9G70_35728.1 -2274.69198 -2271.73682 5.910322 0.01505242 G9G70_17075.1 -6441.992153 -6436.5487 10.886912 0.000968462 G9G70_43159.1 -3459.779524 -3454.81322 9.932614 0.001623758 G9G70_49722.3 -1973.613614 -1965.91579 15.395654 8.72E-05 G9G70_52308.1 -5891.597539 -5856.56249 70.070108 5.72E-17 G9G70_13510.1 -3998.658056 -3995.60513 6.105846 0.013473544 G9G70_23673.1 -3259.995687 -3221.75205 76.487282 2.22E-18 G9G70_12776.1 -3090.044261 -3075.43279 29.222934 6.45E-08 G9G70_37904.1 -11660.02206 -11643.2491 33.545884 6.96E-09 G9G70_31361.1 -1983.290722 -1958.28596 50.009528 1.53E-12 G9G70_35560.1 -3726.459873 -3712.10892 28.701912 8.44E-08 G9G70_15482.1 -2451.96402 -2447.31509 9.297866 0.002294209 G9G70_20090.1 -7500.354631 -7498.2114 4.286454 0.03841721 G9G70_13346.1 -3815.288524 -3808.88122 12.81461 0.000343923 G9G70_3741.1 -6062.68774 -6058.74623 7.883022 0.004990103 G9G70_9963.1 -3016.585576 -3014.10315 4.964852 0.025867529 G9G70_8291.1 -7134.504935 -7132.31045 4.388974 0.036172075 G9G70_3643.1 -3613.271458 -3550.36161 125.819698 3.37E-29 G9G70_58800.1 -7045.65254 -7043.04466 5.215762 0.022383036 G9G70_54074.1 -3799.519435 -3796.75883 5.521214 0.018787217 65 Table B.4 BLAST best hits of positive selected protein-coding sequences in Picea mariana with Arabidopsis genes (Araport 11 database). Picea mariana mRNA ID Arabidopsis thaliana mRNA ID pident qstart qend length qcovs evalue G9G70_725.1 AT1G02660.1 45.7 184 810 644 75 8.97E-157 G9G70_1487.1 AT1G68000.1 59.9 4 230 227 98 1.38E-96 G9G70_2048.1 AT1G28470.1 72.3 36 211 177 62 6.16E-91 G9G70_2442.1 AT1G19220.1 73.6 1 474 477 94 0 G9G70_2518.1 AT5G22850.1 55.1 76 527 459 86 1.41E-178 G9G70_3643.1 AT3G11670.1 61.2 2 518 520 97 0 G9G70_3741.1 AT1G10390.2 51.3 52 949 947 92 0.00E+00 G9G70_4240.1 AT1G71070.1 58.3 1 276 276 58 8.21E-121 G9G70_4241.1 AT1G22240.1 57.2 379 712 334 47 2.70E-132 G9G70_7440.1 AT5G13490.1 71.5 7 372 368 98 6.48E-178 G9G70_8291.1 AT1G19835.3 42.1 78 614 554 53 3.13E-103 G9G70_9963.1 AT5G18610.3 65.8 3 505 509 98 0 G9G70_11476.1 AT5G22510.1 82.2 1 411 411 100 0 G9G70_12776.1 AT2G32910.1 42.2 109 281 180 55 2.10E-39 G9G70_13089.1 AT2G17390.1 45.0 15 273 271 76 2.31E-70 G9G70_13346.1 AT3G24660.1 39.7 22 558 556 95 1.38E-99 G9G70_13510.1 AT2G03340.1 39.6 63 499 447 67 1.52E-89 G9G70_15482.1 AT3G13780.1 35.6 1 117 118 30 1.54E-12 G9G70_17075.1 AT1G05910.1 62.2 51 876 857 94 0 G9G70_20090.1 AT3G48050.2 29.9 817 1249 465 50 7.52E-28 G9G70_22542.2 AT3G04910.1 55.5 17 512 528 56 0 G9G70_22547.2 AT2G22780.1 78.6 11 356 350 97 0 G9G70_23673.1 AT5G03040.2 42.9 99 317 219 39 2.71E-43 G9G70_24781.1 AT5G20850.2 45.3 20 114 95 25 2.76E-14 G9G70_26879.1 AT2G26910.1 64.0 49 998 1010 92 0 66 G9G70_27262.1 AT3G44050.1 48.2 171 1070 964 65 0 G9G70_28286.1 AT3G57070.1 33.6 29 418 405 93 8.37E-56 G9G70_30239.1 AT1G19780.2 67.9 1 299 315 99 8.03E-151 G9G70_31361.1 AT3G15610.1 53.3 6 224 259 95 3.31E-86 G9G70_35560.1 AT4G11400.1 34.3 11 333 329 85 5.33E-42 G9G70_35728.1 AT1G09850.1 52.3 39 361 325 86 3.14E-110 G9G70_36237.1 AT2G25050.1 64.8 1 400 401 35 0 G9G70_37904.1 AT1G02990.5 44.3 1315 1411 97 22 4.73E-11 G9G70_38193.1 AT3G22440.1 27.8 29 545 565 92 2.98E-50 G9G70_43159.1 AT1G65610.1 54.3 7 601 606 96 0 G9G70_44353.1 AT3G10030.2 43.7 79 394 339 80 2.48E-71 G9G70_45605.1 AT1G63020.4 27.5 247 787 596 36 2.87E-48 G9G70_47211.1 AT1G75310.1 23.3 184 298 116 20 1.7 G9G70_49722.3 AT1G03250.4 36.8 91 245 155 59 1.05E-20 G9G70_52308.1 AT1G15240.1 40.5 4 761 765 98 3.57E-159 G9G70_53953.1 AT2G42040.1 73.8 239 280 42 13 5.24E-15 G9G70_54074.1 AT3G25430.1 37.9 88 259 177 60 7.30E-28 G9G70_58800.1 AT1G28520.4 26.3 556 647 99 9 0.42 G9G70_60414.1 AT3G21870.1 51.7 25 228 209 61 1.10E-61 G9G70_62668.1 AT2G01630.3 68.6 84 500 417 76 0 67 Table B.5 Descriptions of BLAST best hit of positive selected protein-coding sequences in Picea mariana with Arabidopsis genes and best hit found against PFAM families and domains. Picea mariana mRNA ID Arabidopsis thaliana mRNA ID Arabidopsis thaliana gene name Description of Arabidopsis thaliana sequence in Araport 11 PFAM ID G9G70_725.1 AT1G02660.1 PLASTID LIPASE2, PLIP2 PLIP2 is a glycerolipid A1 lipase with substrate preference for monogalactosyldiacylglycerol. Expression is induced by ABA. Lipase_3 G9G70_1487.1 AT1G68000.1 ATPIS, ATPIS1, PHOSPHATIDYLINOSITOL SYNTHASE, PHOSPHATIDYLINOSITOL SYNTHASE 1, PIS1 phosphatidylinositol synthase 1 CDP-OH_P_transf G9G70_2048.1 AT1G28470.1 ANAC010, NAC DOMAIN CONTAINING PROTEIN 10, NAC010, SECONDARY WALL-ASSOCIATED NAC DOMAIN PROTEIN 3, SND3 NAC domain containing protein 10 NAM G9G70_2442.1 AT1G19220.1 ARF11, ARF19, AUXIN RESPONSE FACTOR 19, AUXIN RESPONSE FACTOR11, IAA22, INDOLE-3-ACETIC ACID INDUCIBLE 22 Encodes an auxin response factor that contains the conserved VP1-B3 DNA-binding domain at its N-terminus and the Aux/IAA-like domains III and IV present in most ARFs at its C-terminus. The protein interacts with IAA1 (yeast two hybrid) and other auxin response elements such as ER7 and ER9 (yeast one hybrid). ARF19 protein can complement many aspects of the arf7 mutant phenotype and , together with ARF7, is involved in the response to ethylene. In the arf7 arf19 double mutant, several auxin-responsive genes (e.g. IAA5, LBD16, LBD29 and LBD33) are no longer upregulated by auxin. Auxin_resp 68 G9G70_2518.1 AT5G22850.1 not available Eukaryotic aspartyl protease family protein Asp G9G70_3643.1 AT3G11670.1 ATTS02, DGD1, DIGALACTOSYL DIACYLGLYCEROL DEFICIENT 1 Responsible for the final assembly of galactolipids in photosynthetic membranes. Provides stability to the PS I core complex (e.g. subunits PsaD, PsaE). Glycos_transf_1 G9G70_3741.1 AT1G10390.2 DRA2, DRACULA2, NUP98A DRA2 is a homolog of mammalian nucleoporin 98 and a likely component of the nuclear pore complex in Arabidopsis. It positively participates in the control of the hypocotyl elongation response to plant proximity and control of shade induced gene expression. Nucleoportin which redundantly inhibits flowering together with Nup98b through multiple pathways including clock, photoperiod, and age pathways. Gates flowering in a CONSTANS (CO)-independent mode and bypasses the CO checkpoint in photoperiodic signaling and integrated signals from multiple pathways to directly target FLOWERING LOCUS T (FT) for flowering control. Nucleoporin2 G9G70_4240.1 AT1G71070.1 not available Core-2/I-branching beta-1,6-N-acetylglucosaminyltransferase family protein Branch G9G70_4241.1 AT1G22240.1 APUM8, PUM8, PUMILIO 8 Encodes a member of the Arabidopsis Pumilio (APUM) proteins containing PUF domain (eight repeats of approximately 36 amino acids each). PUF proteins regulate both mRNA stability and translation through sequence-specific binding to the 3' UTR of target mRNA transcripts. PUF 69 G9G70_7440.1 AT5G13490.1 AAC2, ADP/ATP CARRIER 2 AAC2, ADP/ATP CARRIER 2 Mito_carr G9G70_8291.1 AT1G19835.3 TCS1, TRICHOME CELL SHAPE 1 TCS1 encodes a coiled-coil domain protein that binds to microtubules and co-localizes with the cortical microtubules. Mutants have defects in trichome branching and hypocotyl elongation. TCS1 interacts with ZWI and appears to be involved in microtubule assembly. FPP G9G70_9963.1 AT5G18610.3 PBL27 Encodes a receptor-like cytoplasmic kinase that is an immediate downstream component of the chitin receptor CERK1 and contributes to the regulation of chitin-induced immunity. Pkinase G9G70_11476.1 AT5G22510.1 A/N-INVE, ALKALINE/NEUTRAL INVERTASE E, ARABIDOPSIS ALKALINE/NEUTRAL INVERTASE E, AT-A/N-INVE, INV-E Encodes a chloroplast-targeted alkaline/neutral invertase that is implicated in the development of the photosynthetic apparatus and nitrogen assimilation in seedlings to control the sucrose to hexose ratio. GDE_C G9G70_12776.1 AT2G32910.1 not available DCD (Development and Cell Death) domain protein Dev_Cell_Death G9G70_13089.1 AT2G17390.1 AKR2B, ANKYRIN REPEAT-CONTAINING 2B Highly homologous to AKR2A. Involved in chloroplast biogenesis. Double mutants of AKR2A and AKR2B have yellow leaves, significantly reduced chloroplast proteins, and no thylakoid membranes. The mRNA is cell-to-cell mobile. Ank G9G70_13346.1 AT3G24660.1 TMKL1, TRANSMEMBRANE KINASE-LIKE 1 member of Receptor kinase-like protein family LRR_1 70 G9G70_13510.1 AT2G03340.1 WRKY DNA-BINDING PROTEIN 3, WRKY3 Encodes WRKY DNA-binding protein 3 (WRKY3). WRKY G9G70_15482.1 AT3G13780.1 not available SMAD/FHA domain-containing protein zf-CCHH G9G70_17075.1 AT1G05910.1 BRAT1, BROMODOMAIN AND ATPASE DOMAIN-CONTAINING PROTEIN 1 Encodes an acetylated histone-binding protein BRAT1. BRAT1 forms a complex with BRP1 to prevent transcriptional silencing. AAA G9G70_20090.1 AT3G48050.2 'SHUTTLE' IN CHINESE, SUO Encodes a large protein with N-terminal bromo-adjacent homology (BAH) and transcription elongation factor S-II (TFS2N) domains and two C-terminal GW (glycine and tryptophan) repeats. It is nuclear and colocalizes with the processing-body component DCP1 in the cytoplasm. SOU is a component of the miRNA pathway and is involved in translational repression. Med26 G9G70_22542.2 AT3G04910.1 ATWNK1, WITH NO LYSINE (K) KINASE 1, WNK1, ZIK4 Serine/threonine protein kinase, whose transcription is regulated by circadian rhythm. Pkinase G9G70_22547.2 AT2G22780.1 PEROXISOMAL NAD-MALATE DEHYDROGENASE 1, PMDH1 encodes an peroxisomal NAD-malate dehydrogenase that is involved in fatty acid beta-oxidation through providing NAD to the process of converting fatty acyl CoA to acetyl CoA. Ldh_1_N G9G70_23673.1 AT5G03040.2 IQ-DOMAIN 2, IQD2 Member of IQ67 (CaM binding) domain containing family. IQ G9G70_24781.1 AT5G20850.2 ATRAD51, RAD51 Encodes a homolog of yeast RAD51. Its mRNA is most abundant in early flower buds and is expressed at high levels in exponentially growing cells in suspension cultures and is induced in response to gamma radiation. Also involved in defense gene transcription during plant immune responses. Rad51 71 G9G70_26879.1 AT2G26910.1 ABCG32, ATABCG32, ATP-BINDING CASSETTE G32, ATPDR4, PDR4, PEC1, PERMEABLE CUTICLE 1, PLEIOTROPIC DRUG RESISTANCE 4 Encodes a member of the PLEIOTROPIC DRUG RESISTANCE family of ATP binding cassette transporters. Required for the formation of a functional cuticle. ABC2_membrane G9G70_27262.1 AT3G44050.1 KINESIN-12E Kinesin involved in spindle assembly. Kinesin G9G70_28286.1 AT3G57070.1 not available Glutaredoxin family protein Glutaredoxin G9G70_30239.1 AT1G19780.2 ATCNGC8, CNGC8, CYCLIC NUCLEOTIDE GATED CHANNEL 8 Encodes a member of the cyclic nucleotide gated channel (CNGC) family that is essential for male reproductive fertility. cNMP_binding G9G70_31361.1 AT3G15610.1 not available Transducin/WD40 repeat-like superfamily protein WD40 G9G70_35560.1 AT4G11400.1 not available ARID/BRIGHT DNA-binding , ELM2 domain and myb-like DNA-binding domain-containing protein ARID G9G70_35728.1 AT1G09850.1 XBCP3, XYLEM BARK CYSTEINE PEPTIDASE 3 Arabidopsis thaliana papain-like cysteine peptidase Inhibitor_I29 G9G70_36237.1 AT2G25050.1 not available actin-binding FH2 (formin 2) family protein FH2 G9G70_37904.1 AT1G02990.5 not available hypothetical protein zf-CW G9G70_38193.1 AT3G22440.1 not available FRIGIDA-like protein Frigida G9G70_43159.1 AT1G65610.1 ARABIDOPSIS THALIANA GLYCOSYL HYDROLASE 9A2, ATGH9A2, ATKOR2, KOR2, KORRIGAN 2 Six-hairpin glycosidases superfamily protein Glyco_hydro_9 G9G70_44353.1 AT3G10030.2 not available aspartate/glutamate/uridylate kinase family protein AA_kinase 72 G9G70_45605.1 AT1G63020.4 NRPD1, NRPD1A, NUCLEAR RNA POLYMERASE D 1A, NUCLEAR RNA POLYMERASE D1A, POL IVA, SDE4, SILENCING MOVEMENT DEFICIENT 2, SMD2 Encodes one of two alternative largest subunits of a putative plant-specific RNA polymerase IV (aka RNA polymerase D). Required for posttranscriptional gene silencing. RNA_pol_Rpb1_R G9G70_47211.1 AT1G75310.1 AUL1, AUXILIN-LIKE 1 auxin-like 1 protein R3H G9G70_49722.3 AT1G03250.4 not available R3H domain protein involved in ribosome biogenesis. R3H-assoc G9G70_52308.1 AT1G15240.1 SNX3, SORTING NEXIN 3 Encodes a member of the Arabidopsis sorting nexin family. PX G9G70_53953.1 AT2G42040.1 not available WRC protein WRC G9G70_54074.1 AT3G25430.1 RRD1 Polynucleotidyl transferase, ribonuclease H-like superfamily protein CAF1 G9G70_58800.1 AT1G28520.4 ATVOZ1, VASCULAR PLANT ONE ZINC FINGER PROTEIN, VOZ1 VOZ transcription factor which acts as positive regulator of several salt-responsive genes. Functionally redundant in salt stress with VOZ2. Rif1_N G9G70_60414.1 AT3G21870.1 CYCLIN P2;1, CYCP2;1 cyclin p2 Cyclin G9G70_62668.1 AT2G01630.3 not available O-Glycosyl hydrolases family 17 protein X8
Thesis/Dissertation
eng
Bioinformatics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Birol, Inanc
The genome of black spruce : genome annotation & analyses
Text