A long-read and short-read transcriptomics approach provides the first high-quality reference transcriptome and genome annotation for Pseudotsuga menziesii (Douglas-fir)

Abstract

Douglas-fir (Pseudotsuga menziesii) is native to western North America. It grows in a wide range of environmental conditions and is an important timber tree. Although there are several studies on the gene expression responses of Douglas-fir to abiotic cues, the absence of high-quality transcriptome and genome data is a barrier to further investigation. Like for most conifers, the available transcriptome and genome reference dataset for Douglas-fir remains fragmented and requires refinement. We aimed to generate a highly accurate, and complete reference transcriptome and genome annotation. We deep-sequenced the transcriptome of Douglas-fir needles from seedlings that were grown under nonstress control conditions or a combination of heat and drought stress conditions using long-read (LR) and short-read (SR) sequencing platforms. We used 2 computational approaches, namely de novo and genome-guided LR transcriptome assembly. Using the LR de novo assembly, we identified 1.3X more high-quality transcripts, 1.85X more "complete" genes, and 2.7X more functionally annotated genes compared to the genome-guided assembly approach. We predicted 666 long noncoding RNAs and 12,778 unique protein-coding transcripts including 2,016 putative transcription factors. We leveraged the LR de novo assembled transcriptome with paired-end SR and a published single-end SR transcriptome to generate an improved genome annotation. This was conducted with BRAKER2 and refined based on functional annotation, repetitive content, and transcriptome alignment. This high-quality genome annotation has 51,419 unique gene models derived from 322,631 initial predictions. Overall, our informatics approach provides a new reference Douglas-fir transcriptome assembly and genome annotation with considerably improved completeness and functional annotation.

Keywords: Pseudotsuga menziesii var. glauca; Pseudotsuga menziesii var. menziesii; de novo assembly; NovaSeq; PacBio Iso-Seq; coastal Douglas-fir; full-length isoform; functional annotation; genome annotation; interior Douglas-fir; long noncoding RNA; reference transcriptome; transcription factors.

Fig. 1.

Comparison of quality of de novo and reference genome-guided assembly of Douglas-fir LR-generated transcriptome. a) Length versus number of unique transcripts and b) transcriptome completeness score. c) Number and percentage of unique transcripts with functional annotation.

Fig. 2.

Summary of gene family and gene ontology assignments in de novo assembled Douglas-fir LR transcriptome. Top ten a,b) taxonomic groups versus unique transcripts, c) GO biological processes, and d) GO molecular function terms versus number of GO terms assigned. White and gray bars represent data from de novo and genome-guided transcriptome assembly, respectively.

Fig. 3.

Trait values of transcripts predicted as lncRNAs in Douglas-fir. a) Mean transcript length, b) ORF length, c) GC content, d) Fickett testscore, e) Hexamer score, and f) lncRNA score of lncRNAs and all other assembled transcripts including protein-coding transcripts.

Fig. 4.

Genome annotation processes and BUSCO completeness. (a–c) Flow chart summarizes the steps taken to obtain genome annotations (b) Annotation v2 (pre-filter) and Annotation v2 and (c) Annotation v1. (d) BUSCO scores for genome assembly* used in genome annotation steps and de novo assembled transcriptome^§ are also shown.

Fig. 5.

Genome annotation evaluation. a) Gene and intron length distribution across genome annotation approaches, and the transcriptome alignment (de novo assembled prior to alignment). The log-scaled values for gene length and intron length reflect improvements in contiguity with the addition of Iso-Seq data. b) Reciprocal BLAST-style analysis was conducted at 2 coverage values for the total set of genes produced from each method. The 50% coverage of the target/query and 80% coverage of target and query are shown. The numerical value at the end of each bar represents the percentage of total sequences that were functionally annotated at that coverage value.