{'fraction_non_alphanumeric': 0.0425531914893617, 'fraction_numerical': 0.01741641337386018, 'mean_word_length': 5.39599533437014, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Differential gene expression associated with a floral scent polymorphism in the evening primrose Oenothera harringtonii (Onagraceae)

Background Plant volatiles play an important role in both plant-pollinator and plant-herbivore interactions. Intraspecific polymorphisms in volatile production are ubiquitous, but studies that explore underlying differential gene expression are rare. Oenothera harringtonii populations are polymorphic in floral emission of the monoterpene (R)-(−)-linalool; some plants emit (R)-(−)-linalool (linalool+ plants) while others do not (linalool- plants). However, the genes associated with differential production of this floral volatile in Oenothera are unknown. We used RNA-Seq to broadly characterize differential gene expression involved in (R)-(−)-linalool biosynthesis. To identify genes that may be associated with the polymorphism for this trait, we used RNA-Seq to compare gene expression in six different Oenothera harringtonii tissues from each of three linalool+ and linalool- plants. Results Three clusters of differentially expressed genes were enriched for terpene synthase activity: two were characterized by tissue-specific upregulation and one by upregulation only in plants with flowers that produce (R)-(−)-linalool. A molecular phylogeny of all terpene synthases identified two putative (R)-(−)-linalool synthase transcripts in Oenothera harringtonii, a single allele of which is found exclusively in linalool+ plants. Conclusions By using a naturally occurring polymorphism and comparing different tissues, we were able to identify candidate genes putatively involved in the biosynthesis of (R)-(−)-linalool. Expression of these genes in linalool- plants, while low, suggests a regulatory polymorphism, rather than a population-specific loss-of-function allele. Additional terpene biosynthesis-related genes that are up-regulated in plants that emit (R)-(−)-linalool may be associated with herbivore defense, suggesting a potential economy of scale between plant reproduction and defense. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08370-6.

Keywords: linalool, RNA-Seq, terpene synthase, volatile, transcriptome, floral scent, floral fragrance Background Floral fragrance is a complex trait that is typically composed of tens to hundreds of volatile organic compounds (VOCs) whose ecological roles include the attraction of pollinators, manipulation of pollinator behavior, and defense against herbivores [1][2][3][4][5][6][7]. Considerable qualitative and quantitative variation in floral fragrance has been documented, with over 1700 floral volatiles having been described from more than 900 angiosperm species [8]. In addition to substantial interspecific variation, floral fragrance has also been shown to differ intra-specifically [2,9]. The evolutionary and ecological mechanisms that maintain intraspecific variation in floral fragrance constitute an active area of research and have been attributed to differences in community composition, pollinator preferences, phenotypic plasticity, and genetic drift [10][11][12][13].
The biosynthetic pathways that produce floral volatiles are increasingly wellcharacterized [14], and recent studies have used transcriptomics-based approaches to identify homologous genes from these pathways in non-model plant species [15][16][17][18]. Despite these advances [19,20], relatively few studies have identified the genetic mechanisms that impact how these pathways modify floral scent or how selection acts on the resulting intraspecific variation.
Volatile terpenes comprise a major source of biosynthetic diversity in scents, primarily 10carbon monoterpenes and 15-carbon sesquiterpenes [21]. Terpene synthase (TPS) enzymes, which are involved in the biosynthesis of these secondary metabolites, are encoded by a moderately large gene family in angiosperms [22]. The TPS genes exist in seven to eight subfamilies that show either lineage-specific (e.g. exclusively gymnosperm TPSs) or molecule/lineage-specific (e.g. monocot sesquiterpene TPSs) affinities. Despite the widespread characterization of the gene family and their products across plant lineages [22,23], the mechanisms maintaining intraspecific variation in the production of volatile terpenes remain poorly understood. The functional implications of intraspecific terpenoid variation are epitomized by the culinary herb Thymus vulgaris (Lamiaceae), whose distinct "chemotypes" differ in herbivore defense, frost sensitivity, and allelopathic interactions with neighboring plants [24,25]. Genetic crosses helped to document an epistatic network of mendelian loci responsible for volatile chemotypes in T. vulgaris [26,27]. However, beyond a few such model systems, the genetic controls of chemotype variation remain largely unknown, and this knowledge is needed to understand how conflicting selective forces shape chemical polymorphism. To date, there have been few opportunities to study the genetic underpinnings of volatile terpenoid polymorphism in floral scent.
Oenothera; [30]), O. acutissima (sect. Lavauxia; [31] Linalool in floral tissues can function as a pollinator attractant or as a defense compound [32][33][34][35], whereas linalool emitted by vegetative tissues has been implicated in both direct and indirect plant defenses [34]. Flowers of O. harringtonii are pollinated by the widespread hawkmoth species Hyles lineata and Manduca quinquemaculata [36]. Parallel studies suggest that these hawkmoths are likely attracted to the strong floral scent and high visual contrast of O. harringtonii flowers [37,38]. In addition to its role in plant reproduction, the primary pollinator H. lineata is also an herbivore of the plant, with its caterpillars feeding on leaves and floral tissues. Furthermore, the developing seeds of O. harringtonii are consumed by caterpillars of the moth genus Mompha [39], whose adults oviposit on flowers at dusk [40]. Considering these selective forces, it is likely that both floral scent and volatiles emitted by vegetative tissues play important roles in the complex fitness landscape of this species. However, understanding how these selective forces mold volatile terpene variation requires the identification of genes, including (R)-(-)-linalool synthase and related terpene synthases, that are associated with scent polymorphism in O. harringtonii. In this way, identification and subsequent manipulation of the (S)-(+) linalool synthase gene in a wild tobacco, Nicotiana attenuata, led to novel insights about geographic variation in selective pressures and the importance of terpenoid volatile background in a tri-trophic/plant defense context [41].
Within the family Onagraceae, the gene encoding (S)-(+)-linalool synthase was identified from Clarkia breweri and Clarkia concinna in the first study of biosynthetic genes responsible for floral scent [42]. However, (R)-(-)-linalool synthase has not been identified in any evening primrose species to date, and only a single, partial TPS sequence has been identified thus far in Oenothera [43], one of the most species-rich lineages within the family [44]. In this study, we use RNA-Seq to identify differentially expressed genes (DEGs) correlated with linalool production in O. harringtonii by sequencing transcriptomes of six different tissues from three biological replicates of each of two chemotypes: (1) flowers that emit (R)-(-)-linalool (linalool+) and (2) flowers that do not emit (R)-(-)-linalool (linalool-). We use this information, in combination with a phylogenetic analysis of terpene synathases, to identify a gene that likely encodes (R)-(-)-linalool synthase in O. harringtonii and identify an allele that is specific to linalool+ individuals. This study includes the first transcriptomes assembled for O. harringtonii, providing critical genomic information to inform ongoing studies of intraspecific floral scent variation within this species. In addition, the preliminary identification of (R)-(-) linalool synthase and other differentially expressed transcripts is a first step towards a broader, comparative analysis of the molecular evolution of floral scent among species of Oenothera that differ in mating system, pollination mode, life history strategy and, thus, the direction and intensity of natural selection (see [45]). As such, results of the present study provide a key foundational step in understanding the ubiquity and relevance of (R)-(-)-linalool throughout flowering plants.

Floral Fragrance Analysis
Floral scent profiles were collected from replicated linalool+ and linalool-plants from which expressed transcripts were subsequently harvested; volatile chemical data are summarized in Table 1. Across all individuals, the floral scent was dominated by the monoterpene olefin (Ε)β-ocimene, with small amounts of other monoterpenes ((Ζ)-β-ocimene and β-myrcene) and sesquiterpenes (β-caryophyllene and (Ε,Ε)-farnesol). The flowers of three individuals emitted large amounts of (R)-(-)-linalool (linalool+), whereas those from another three individuals did not emit any (Fig. 1A). This pattern is consistent with phenotypic data from the source populations for these plants (Additional File 4), confirming that linalool chemotypes breed true in greenhouse-grown plants (see Methods). Other volatiles (ocimene epoxide, methyl geranate, isophytol, α-humulene, and methyl benzoate) were emitted in small amounts by some individuals but not others, but the presence of these compounds varied independently of linalool chemotype (as they did in field-collected samples; Additional File 4) and are therefore likely due to

Transcriptome Sequencing, Assembly, and Annotation
RNA was successfully sequenced from six tissues (stamens, stigma/style, petals, hypanthium, bud, and leaf) for all individuals, except for individuals D and F, for which bud and leaf samples were not obtained (see Methods). Across all 32 samples, 1.57 billion paired reads were sequenced (Additional file 1). On average, 95.77% of reads passed quality trimming, and 99.04% of these high-quality reads were retained after filtering for plastid genes and rRNA. In total, 94.86%, or 1.49 billion reads, were retained after filtering. After normalization, a total of 33 million reads were assembled into a reference assembly consisting of 489,895 contigs and 244,284 components. Through alignment and abundance estimation using RSEM, 72% of cleaned reads aligned to the reference assembly. Median contig length was 393 base pairs, average contig length was 676.8 base pairs, and transcript N50 was 1047 base pairs. As determined using the Trinotate pipeline [46], 27.6% of contigs had either a blastx or blastp hit to Swiss-Prot, and 25.5% of contigs had a gene ontology term assigned from either blast or pfam.
In a blastx search of the assembly against an Arabidopsis database [47], 174,068 contigs had a hit in the Arabidopsis thaliana protein database; 20,480 of these were unique Arabidopsis protein hits. In a tBLASTn search of Arabidopsis proteins against the assembly, 31064 Arabidopsis peptides had a hit and 14,614 of these were unique. In total, there were 11,754 reciprocal best hits between the transcriptome assembly and the Arabidopsis proteome, with an ortholog hit ratio of 0.34 and an average collapse factor of 2.13.

Differential Expression and GO Enrichment Analysis
Principal components analysis was conducted to test whether it is appropriate to group biological replicates together during subsequent differential expression analysis. In general, biological replicates of the same tissue and chemotype clustered together (Fig. 2). Additionally, samples of the same linalool chemotype grouped together, suggesting significant gene expression differences between the chemotypes. In total, edgeR identified 6,577 differentially expressed genes (DEGs) among the six tissues and between linalool+ and linalool-individuals, with an overall pattern suggesting that tissue identity explains most of the observed differential expression (Fig. 3). All DEGs were grouped into 12 clusters by cutting the DEG dendrogram at 65% of its height (Fig. 4). Gene ontology (GO) enrichment analysis was performed for each cluster to identify categories of genes that are overrepresented (see Additional file 2 for a complete list of enriched GO terms for each cluster).
Of the twelve clusters, two (clusters 10 and 11) were characterized by differential expression of genes between linalool+ and linalool-individuals, suggesting that genes in these clusters may be involved in linalool biosynthesis. Genes in cluster 10 (n=43) were upregulated in all tissues from linalool+ individuals, whereas genes within cluster 11 (n=46) were downregulated in all tissues from linalool+ individuals. In general, both of these clusters were characterized by several genes that lack annotations in the Trinotate report (26/43 and 38/46 genes for each cluster, respectively), either because they are unique or significantly differ from a homolog in the reference database.
To obtain an overview of the functions of genes within each cluster, we conducted GO enrichment. Cluster 10 is enriched for terpene synthase activity (Table 2), and transcripts in this cluster with this higher-level GO term (TPS activity) include those annotated as terpene synthase 10 (DN74292_c2_g1), γ-terpinene synthase/(+)-α-terpineol synthase (TRINITY_DN62538_c3_g2) and santalene synthase (DN76019_c5_g5). Terpene synthase 10 is known to produce (R)-(-)-linalool as a primary product in Arabidopsis thaliana, with β-myrcene and (E)-β-ocimene as minor products [48]. γ-terpinene and (+)-α-terpineol are monoterpenes, and although santalene is a sesquiterpene, terpene synthases responsible for its biosynthesis are classified with the TPS-b monoterpene synthases [49]. γ-terpinene and santalene are not present in the floral scent of O. harringtonii, whereas α-terpineol is only present in small amounts, inconsistently, in some populations (e.g. FLO; Additional File 4). Thus, either the appropriate biosynthetic enzymes are not translated from these transcripts or, if they are, their eventual products are either metabolized into other secondary compounds, or they simply are not volatilized. These results suggest that linalool+ individuals have the capacity to synthesize terpenes in all six tissues included here; therefore, differences between the linalool+ and linaloolchemotypes include terpene synthase expression in both floral and vegetative tissues. In addition to being common components of floral fragrance, terpenes are known to be involved in constitutive and induced defenses in both floral and non-floral tissues [4,50]. Terpenes are also involved in indirect defense, as they can attract members of the third trophic level [51][52][53][54][55].

Overexpression of TPS genes across all tissues including leaves suggests that O. harringtonii
populations emitting R-(-) linalool from flowers may also be better-defended.
Cluster 11 did not have any significant GO term enrichments, possibly because only a small number of the transcripts in the cluster were able to be annotated. The few genes in cluster 11 that do have annotations are not associated with terpene biosynthesis (Additional File 2). Two additional clusters contained significantly enriched GO terms related to terpene biosynthesis: clusters 1 (n=232) and 6 (n=542). Cluster 1 comprised genes upregulated in buds and petals but downregulated in leaves. Significant (FDR > 0.05) GO term enrichments for this cluster include isoprenoid biosynthetic process, terpenoid biosynthetic process, and isopentenyl diphosphate biosynthetic process, among others (Table 2 and Additional file 2). Cluster 6 comprised genes upregulated in petals and downregulated in leaves. Significant GO term enrichments for this cluster include terpene synthase activity, sesquiterpene synthase activity, isoprenoid biosynthetic process, and terpenoid biosynthetic process (Table 2 and Additional file 2). Neither cluster showed substantial differentiation between the two linalool chemotypes, which, in combination with the enriched GO terms for each cluster, is an indication that the genes contained within those clusters may be involved in the biosynthesis of floral scent volatiles that are shared between the two chemotypes. Thus, the expression patterns of these two clusters show that certain terpenoid scent compounds are likely to be produced primarily in the petals of O.
harringtonii, and that buds may be synthesizing floral scent or floral scent precursors 24-48 hours prior to anthesis.

Phylogenomic Analysis of Terpene Synthases and (R)-(-)-linalool Polymorphisms
To further characterize and identify terpene synthases in the transcriptome, regardless of their expression pattern, a profile Hidden Markov Model (pHMM) for terpene synthase (TPS) proteins was created based on a multiple sequence alignment of TPS sequences from van Schie et al. [23]. The pHMM successfully identified 40 TPS inferred protein sequences in the O.
harringtonii reference assembly (Fig. 5). While it is likely that this is an overestimate due to the co-assembly of multiple individuals for the reference transcriptome, the 40 sequences were successfully placed into five of the seven TPS subfamilies with high confidence. According to the TPS subfamily circumscription in van Schie et al. [23], we found: (1) two distinct sequences of diterpene synthases (TPSc); (2) a single sequence of (S)-(+)-linalool synthase (TPSf); (3) 11 sequences of a second group of diterpene synthases (TPSe) with seven distinct groups of genes (components arising from the Trinity assembly); (4) eight sequences of sesquiterpene synthases (TPSa) with six distinct components, and (5) 18 sequences of monoterpene synthases (TPSb). No sequences were reconstructed within the TPSd or TPSg subfamilies (Fig. 5); the former (TPSd) comprises exclusively gymnosperm terpene synthases and the latter (TPSg) comprises a group of acyclic monoterpene synthases related to ocimene synthases in Antirrhinum [23].
Genes belonging to TPSc, TPSe, and TPSf largely showed no differential expression between either tissue or chemotype (Fig. 5). Two genes within the TPSe clade belong to DEG cluster 3 ( Fig. 4 and 5) and were upregulated in leaves, suggesting that they may play a role in defense against herbivory rather than floral fragrance. Some sesquiterpene synthase genes (the TPSa clade) were upregulated in floral tissues of linalool+ individuals (see cluster 6 in Fig. 4; To-date, the most intensive genetic and physiological study of floral fragrance emission in Onagraceae has focused on a pair of sibling species in the genus Clarkia [42,43,56,57]. Most Clarkia species are bee-pollinated or autogamous [58], and range from unscented to mildly scented, with the exception of C. breweri, which produces a strong floral fragrance rich in (S)-(+)-linalool as a derived trait associated with a shift to hawkmoth pollination [57]. Pichersky et al. [56] found that (S)-(+)-linalool is primarily emitted from the petals and, to a lesser extent, the pistil of C. breweri flowers, whereas comparative analyses of the closest relative, C. concinna, revealed linalool production only in stigmatic tissues. The authors inferred that increased levels of linalool synthase (LIS) gene expression, expansion of tissues in which it is expressed, and enzymatic upregulation all contributed to the evolutionary gain-of-function of strong (S)-(+)linalool emission from the flowers of C. breweri [29]. Our phylogenomic analysis placed a single sequence of (S)-(+)-linalool synthase (TPSf) from O. harringtonii as sister to (S)-(+)linalool synthase from Clarkia (Fig. 5). While we cannot confirm that this putative (S)-(+)linalool synthase in O. harringtonii is non-functional, it is clear that no detectable (S)-(+)linalool is emitted from the flowers ( Fig. 1B; Table 1), the transcript is not differentially expressed, and it is expressed at a low level.
While both clades are strongly supported, their relationships to other TPSb genes from other species are not well resolved. Genes from O. harringtonii within the TPSb 1 clade show some evidence of upregulation in floral tissues and bud (clusters 1 and 6) with no evidence of chemotype-related differential expression. As the two chemotypes differ primarily in the presence of a single monoterpene, (R)-(-)-linalool, we would expect that only one of the TPSb subclades would include a gene that is upregulated in linalool+ floral tissues. Two contigs from clade TPSb 2 belong to cluster 10, which showed strong differential expression related to chemotype (see rows annotated as cluster 10 in Figure 5). Although cluster 10 includes genes that are differentially expressed (upregulated in linalool+ individuals) for all tissues including leaves, the signal for upregulation of these two transcripts in leaves is relatively weak. In fact, we detect no (R)-(-)-linalool emitted from the leaves of individuals with linalool+ flowers (Additional file 3). Furthermore, the TPSb 2 subclade is sister (with moderate bootstrap support of 68) to a clade that includes other linalool synthases, including the (R)-(-)-linalool synthase gene from Artemesia annua [59]. When the phylogeny, differential expression, and floral fragrance data are taken together, they strongly suggest that the two transcripts belonging to cluster 10 in citriodora, or is significantly diverged in sequence (sequence similarity to the B. citriodora CDS is 64%). Further sequencing and functional analysis would be needed to determine which is the case. As RNA-Seq has been shown to reliably predict genes that are significantly associated with a treatment [60][61][62][63] and given the overwhelming difference between the abundance of reads mapped to the candidate (R)-(-)-linalool synthase transcripts in linalool+ and linaloolindividuals, we suggest that RT-qPCR validation is unnecessary in this case. Future studies that focus exclusively on the function of this particular gene in Oenothera in more detail may benefit from the sequences described here, as a template for primer design or as a candidate locus in association studies.
To investigate if there are differences in the coding region of (R)-(-)-linalool synthase associated with either chemotype, we compared SNPs between individuals. Calling SNPs for each individual revealed that all linalool+ individuals are homozygous and share the same haplotype which matched the reference (merged consensus) sequence, with the exception of two SNPs at the 5' and 3' ends of the coding region (Fig. 6), which may be due to sequencing or assembly error (regions where read depth decreases). Conversely, all linalool-individuals are heterozygous, and do not share the same haplotype as each other or the linalool+ individuals. Terpenes are known to be regulated both transcriptionally and post-transcriptionally. In Arabidopsis, for example, differences in the amount of (E)-β-ocimene and (E,E)-α-farnesene emitted between two ecotypes is controlled by allelic variation of two terpene synthase genes, whereby one of the two genes in each ecotype produces a non-functional transcript [64].
Transcription factors regulating a variety of terpenes have also been identified; in cotton WRKY1 positively regulates (+)-ẟ-cadinene synthase (CAD1) [65], and in kiwifruit, NAC transcription factors were found to regulate the expression of TPS1 in ripe fruit [66]. Posttranscriptional regulation includes the action of miRNAs, which have been identified in plants such as Ocimum basilicum [67] and Ferula gummosa [68]. The existence of SNPs between transcripts from each chemotype may suggest differences in functionality of each linalool synthase; however, few of the SNPs are shared among linalool-plants, such that the regulation of (R)-(-)-linalool synthase in O. harringtonii may be associated only with upregulation and the allele found in linalool+ individuals. Further sequence analysis and more extensive sampling is required to determine whether this apparent allelic diversity shows strong evidence of population structure. Additionally, a high-quality draft genome sequence for O. harringtonii would allow us to better detect variation that may be associated with regulatory regions or copy number variation.

Conclusions
This study presents the first transcriptome assembly and analysis of differential expression in Oenothera harringtonii, a species with intraspecific variation for the emission of (R)-(-)-linalool in its floral fragrance. In this study, RNA-Seq was used to compare transcriptomes of six different Oenothera harringtonii tissues from six individuals (three linalool+, three linalool-). 6577 differentially expressed genes were identified, and these DEGs were sorted into 12 clusters. Three clusters were significantly enriched for GO terms related to floral fragrance biosynthesis. Two of these clusters showed upregulation in petals, suggesting that floral scent biosynthesis is mostly localized to petals. The third cluster showed upregulation in all linalool+ tissues and contained a significantly enriched GO term for terpene synthase.

Plant Material
Parental individuals were grown from seed collected in natural populations of

Floral Scent Analysis by GC-MS
Prior to analysis, hexane-eluted floral headspace samples were concentrated to a uniform volume of 50μl using gaseous N2, and 5μl of 0.03% toluene in hexane (23ng) was added as an internal standard. One µl aliquots of each sample were injected into a Shimadzu GC-17A gas chromatograph equipped with a Shimadzu QP5000 quadrupole electron ionization (EI) mass spectrometer (Shimadzu Scientific Instruments, Columbia, MD) as a detector. All analyses were made using splitless injections on a polar GC column (0.25mm diameter, 30m length, 0. harringtonii (Florence) were chromatographed using the same protocol, comparing retention times for the leading edge of the linalool peak with those of the authentic standards.

RNA Extraction
RNA was extracted from the same flowers used in fragrance analysis. While still attached to the plant, each flower was dissected into four parts: (1) stamens, (2) petals, (3) stigma and style, and (4) hypanthium and sepals. Additionally, RNA was extracted from flower buds  hours before anthesis and from young leaf tissues from the center of the rosette (excluding individuals D and F, two linalool-individuals for which bud and leaf samples could not be obtained due to insufficient growth post initial sampling). After tissues were excised from the plant, they were immediately placed into either 15ml conical tubes (petals, hypanthium, bud) or 2ml microcentrifuge tubes (stigma and style, stamens, leaf) and submerged in liquid nitrogen.
Assembly quality was evaluated using several different metrics. Percent of reads mapped to the assembly using bowtie2 v2.2.6 [75] was compared across all samples to ensure that reads mapped back evenly across all individuals and tissue types. Mean contig length, median contig length, and contig N50 were calculated using the "TrinityStats.pl" script provided with Trinity.
Annotation statistics were calculated according to O'Neil and Emrich [79], including the number of hits and unique hits in a BLASTx and tBLASTn (v2.2.30) search to the Arabidopsis thaliana proteome, number of reciprocal best hits, average ortholog hit ratio, and average collapse factor.

Differential Expression Analysis
Differential expression analysis followed the protocol recommended by Trinity [77,78] and scripts provided with Trinity; these scripts use R v3.3.0 [89] and several R packages (edgeR in each sample were mapped to the assembled contigs to determine gene expression levels using RSEM v1.2.15 [90]. In this analysis, there were six tissue types (petals; stigma and style; stamens; hypanthium and sepals; leaf; bud), two conditions (linalool+ and linalool-), and three biological replicates (excluding leaf and bud samples for linalool-, which only had one biological replicate).
Expression levels were calculated with the R package edgeR v3.12.0 [91,92] to determine differentially expressed genes (DEGs). Pairwise comparisons between linalool+ and linalool-floral organs, among linalool-floral organs, and among linalool+ tissue types were conducted; a gene was considered differentially expressed if it had at least a four-fold change in expression and a p-value less than 0.001. This approach captures differences between floral organ expression and chemotype expression while minimizing the number of DEGs extracted that are not involved in floral fragrance biosynthesis. Clusters of DEGs with a similar expression pattern were created by cutting the DEG dendrogram at 65% of its height.

GO Enrichment Analysis and Terpene Biosynthesis Gene Identification
Gene ontology (GO) enrichment analysis, performed with scripts provided in the Trinity software package, was used to identify gene categories that are overrepresented in each cluster relative to the assembly. Phylogenetic inference was used to identify and classify terpene synthases (TPS) in the reference transcriptome assembly. To classify putative terpene synthases into the seven described subfamilies a-f [23,42,93,94], all of the sequences included in van Schie et al. [23] were downloaded and verified by reproducing a tree that circumscribed the seven TPS subfamilies. The verified alignment was then used to create a profile Hidden Markov Model using HMMER v3.1b (hmmer.org). The O. harringtonii reference assembly was searched for sequences explained by the pHMM using hmmscan (e-value 1e-10). All putative terpene synthesis protein sequences were aligned using MAFFT v7.310 [95] (--maxiterate 1000 -localpair), and a phylogeny was inferred using RAxML v8.2.11 [96] with the PROTGAMMAWAG model of sequence evolution and 100 fast bootstrap replicates.

SNP Analysis
To determine the distribution of SNPs in the putative (R)-(-)-linalool synthase coding sequence, a putative full length transcript was created by merging partial transcripts identified by the phylogenetic and DEG analyses. Three transcripts were merged that formed a clade (DN62538_c3_g1_i1, DN62538_c3_g2_i2, and DN74292_c2_g1_i11) and could be aligned with only a single mismatch. A web BLAST search against the GenBank database of non-redundant proteins resulted in alignment over approximately 52% of the length of (R)-(-)-linalool synthase from the lemon myrtle Backhousia citriodora (Myrtaceae), suggesting that the consensus sequence produced by merging the transcripts may not represent a full length transcript of (R)-(-)-linalool synthase. However, the merged O. harringtonii (R)-(-)-linalool synthase transcript and the lemon myrtle transcript (GenBank accession AB438045) were only 64% identical and it is possible that the full-length transcript in Oenothera harringtonii may be shorter. SNPs were then identified by mapping reads to the merged transcript from each RNA-Seq library. All reads identified in the DEG analysis as mapping to the three unmerged transcripts were extracted and variants were discovered using the GATK best practices workflow for germline SNPs (https://software.broadinstitute.org/gatk/best-practices/workflow?id=11145).
Briefly, reads were re-mapped to the merged transcript and duplicate mappings were marked and removed before variants were called using HaplotypeCaller in GATK. Phasing of alleles within individuals was done with WhatsHap [97], after which two alternate phased haplotypes were extracted using bcftools (https://samtools.github.io/bcftools/). Each of the phased allele sequences, along with several linalool synthase accessions from the web BLAST search, were added to the previously generated alignment of terpene synthase sequences using the --addalignment feature in MAFFT. A new gene family phylogeny was generated using RAxML as previously described.

Consent for Publication
Not applicable.

Availability of Data and Material
The Illumina RNA-Seq reads supporting the conclusions of this article are available from the NCBI Sequence Read Archive (BioProject PRJNA432694). The reference transcriptome assembly, mapping results, terpene synthase HMM and the alignment of terpene synthase sequences for phylogenetic analysis, and files used for the SNP analysis are currently available from the following link: https://doi.org/10.5061/dryad.np5hqbzrz.

Competing Interests
The authors declare that they have no competing interests.   Table S1) and retention time (min). The identity of all compounds except isophytol (*) and ocimene epoxide (*) was confirmed through cochromatography and mass spectral comparison with synthetic standards.    harringtonii transcript is shown on the right, including its cluster number (CL) if the gene was differentially expressed.