Allelic differences of clustered terpene synthases contribute to correlated intraspecific variation of floral and herbivory‐induced volatiles in a wild tobacco

: Plant volatile emissions can recruit predators of herbivores for indirect defense and attract pollinators to aid in pollination. Although volatiles involved in defense and pollinator attraction are primarily emitted from leaves and flowers, respectively, they will co‐evolve if their underlying genetic basis is intrinsically linked, due either to pleiotropy or to genetic linkage. However, direct evidence of co‐evolving defense and floral traits is scarce. We characterized intraspecific variation of herbivory‐induced plant volatiles (HIPVs), the key components of indirect defense against herbivores, and floral volatiles in wild tobacco Nicotiana attenuata. We found that variation of (E)‐฀‐ocimene and (E)‐฀‐bergamotene contributed to the correlated changes in HIPVs and floral volatiles among N. attenuata natural accessions. Intraspecific variations of (E)‐฀‐ocimene and (E)‐฀‐bergamotene emissions resulted from allelic variation of two genetically co‐localized terpene synthase genes, NaTPS25 and NaTPS38, respectively. Analyzing haplotypes of NaTPS25 and NaTPS38 revealed that allelic variations of NaTPS25 and NaTPS38 resulted in correlated changes of (E)‐฀‐ocimene and (E)‐฀‐bergamotene emission in HIPVs and floral volatiles in N. attenuata. Together, these results provide evidence that pleiotropy and genetic linkage result in correlated changes in defenses and floral signals in natural populations, and the evolution of plant volatiles is probably under diffuse selection. the key components of indirect defense against herbivores, and ﬂoral volatiles in wild tobacco Nicotiana attenuata . (cid:1) We found that variation of ( E )- b -ocimene and ( E )- a -bergamotene contributed to the correlated changes in HIPVs and ﬂoral volatiles among N. attenuata natural accessions. Intraspeciﬁc variations of ( E )- b -ocimene and ( E )- a -bergamotene emissions resulted from allelic variation of two genetically co-localized terpene synthase genes, NaTPS25 and NaTPS38 , respectively. Analyzing haplotypes of NaTPS25 and NaTPS38 revealed that allelic variations of NaTPS25 and NaTPS38 resulted in correlated changes of ( E )- b -ocimene and ( E )- a -bergamotene emission in HIPVs and ﬂoral volatiles in N. attenuata . (cid:1) Together, these results provide evidence that pleiotropy and genetic linkage result in correlated changes in defenses and ﬂoral signals in natural populations, and the evolution of plant volatiles is probably under diffuse selection.


Introduction
Plant volatile emissions mediate interactions between plants and their friends and foes at a distance, which significantly affect plant fitness in nature (Dicke, 1994;Baldwin, 2010). When attacked by herbivores, many plants emit herbivory-induced plant volatiles (HIPVs) from their leaves, which attract predators or parasitoids of herbivores, as well as deter herbivore oviposition, thus reducing herbivore loads and increasing the chance of plant survival (Turlings & W€ ackers, 2004;Halitschke et al., 2008;Degenhardt et al., 2009;Allmann & Baldwin, 2010;Schuman et al., 2012). For many flowering plants, floral scent is a key signal that attracts pollinators to transfer pollen, and hence is critical for their reproductive success (Raguso, 2008). Interestingly, many of the same volatile organic compounds are found in both floral volatiles and HIPVs, including many terpenoids (Tholl, 2006;Raguso, 2016).
Terpenoids represent the largest class of plant-derived natural products (Boutanaev et al., 2015), and are produced from terpene precursors: isopentenyl diphosphate, geranyl diphosphate or farnesyl diphosphate, by their 'signature enzymes' -terpene synthases (TPSs; Tholl, 2006;Degenhardt et al., 2009). Because terpene precursors are commonly found in different plant tissues and species, changes in terpenoid abundance, in both floral volatiles and HIPVs, are largely determined by the expression and function of TPSs (van Schie et al., 2006;Tholl, 2006;Nagegowda, 2010;Irmisch et al., 2012). Most previous studies have investigated only floral and foliar expression of TPSs and terpenoid emissions in isolation. However, Lee et al. (2010) found that a single homoterpene synthase, expressed in both Arabidopsis flowers and herbivore-induced leaves, resulted in the biosynthesis of (E,E)-4,8,12-trimethyltrideca-1,3,7,11-tetraene in both HIPVs and floral volatiles. Furthermore, recent genomic studies have shown that many TPSs are co-localized in the genome and form different clusters (defined as two or more genes located near to each other in the genome that were evolved from the same common ancestral gene or involved in the same biosynthetic pathway) (Falara et al., 2011;Chen et al., 2019Chen et al., , 2020, which can result in strong genetic linkages among TPSs that are expressed in different tissues. These studies suggest that pleiotropy or genetic linkage can result in the co-evolution of floral and foliar terpenoid emissions, which may lead to correlated changes in plant-pollinator and plant-herbivore interactions. However, direct evidence of this inference remains scarce. In Nicotiana attenuata, a wild tobacco native to the North American Great Basin Desert, many of the same terpenoids are found in both HIPVs and floral volatiles (Kessler & Baldwin, 2007;Steppuhn et al., 2008;Wu et al., 2008;Gaquerel et al., 2009;Schuman et al., 2009). Analysing the HIPVs in five natural accessions revealed that herbivory-induced (E)-b-ocimene and (E)-a-bergamotene emissions are variable in N. attenuata (Steppuhn et al., 2008;Wu et al., 2008;Schuman et al., 2009). Through manipulation of individual HIPVs in plants, studies showed that (E)-a-bergamotene can significantly increase the attraction of predatory bugs to N. attenuata (i.e. Geocoris spp.) and improve plant fitness under herbivore attack (Kessler & Baldwin, 2001;Halitschke et al., 2008;Schuman et al., 2012). Moreover, increased (E)-a-bergamotene emission in leaves by ectopic expression of a maize TPS gene in N. attenuata also suggested that (E)-a-bergamotene may interact with green leaf volatiles to improve plant fitness in populations of emitting or nonemitting plants (Schuman et al., 2015).
Although N. attenuata is a predominantly self-fertilizing species, it opportunistically attracts pollinators for outcrossing (Sime & Baldwin, 2003;Kessler et al., 2008;Guo et al., 2019). Analysis of the floral volatiles of N. attenuata revealed that both (E)-b-ocimene and (E)-a-bergamotene are also present (Kessler & Baldwin, 2007). While the emission of (E)-a-bergamotene in corolla tubes increases Manduca sexta moth-mediated pollination success (Zhou et al., 2017), (E)-b-ocimene can elicit distinct electrophysiological responses in the M. sexta antennae and may act as a repellent to the moths (Fraser et al., 2003;Kessler & Baldwin, 2007). However, whether the emission of (E)-b-ocimene and (E)-a-bergamotene in HIPVs and floral volatiles show correlated changes among different genotypes remains unknown.
Here, we first screened the floral volatiles and HIPVs of natural accessions collected from different habits within the range of N. attenuata's geographical distribution. We found that (E)-bocimene and (E)-a-bergamotene in HIPVs and floral volatiles showed correlated changes among N. attenuata genotypes. Using both forward and reverse genetic tools, we characterized the genetic mechanisms underlying the variation of (E)-b-ocimene and (E)-a-bergamotene in HIPVs and floral volatiles, respectively. The results showed that the within-species variation of (E)-b-ocimene and (E)-a-bergamotene results from expression changes in two clustered TPS genes in the TPS-b clade located on chromosome 2: NaTPS25 and NaTPS38, respectively. To further examine whether these two TPS genes contributed to the correlated changes of HIPVs and floral volatiles among N. attenuata genotypes, we characterized the haplotypes of NaTPS25 and NaTPS38 in N. attenuata and their associations with HIPV and floral volatiles. The results demonstrated that allelic differences in the NaTPS25 and NaTPS38 gene cluster result in correlated changes in both floral volatiles and HIPVs in natural genotypes of N. attenuata.

Plant materials and growth conditions
Seeds of N. attenuata Torrey ex. Watson (Solanaceae) natural accessions (genotypes) were collected from a selection of habitats within the range of the species in the southwestern United States (Supporting Information Table S1) and inbred for one generation in the glasshouse in Jena (Germany) except for two inbred lines: Utah and Arizona, which have been self-fertilized for 30 and 22 generations in the glasshouse, respectively. To develop the advanced inter-cross recombinant inbred line (AI-RIL) population, Utah genotype (self-fertilized for 30 generations) and Arizona genotype (self-fertilized for 22 generations) were crossed to generated F 1 plants, which were then self-fertilized to generate 150 F 2 plants (Zhou et al., 2017). Briefly, in each generation from F 2 to F 6 , we intercrossed c. 150 progeny using a random mating and equal contribution crossing design. For generation F 7 , two seeds from each of the crosses at F 6 were germinated and used for the single-seed descendant inbreeding process. In total, five generations of inbreeding were conducted (from F 6 to F 11 ).
Seeds were germinated following the protocol described by Kr€ ugel et al. (2002): seeds were sterilized, placed on Gamborg's B5 medium and kept under LD (16 h : 8 h, light : dark) in a growth chamber (Percival, Perry, IA, USA) for 10 d, and then transferred to small TEKU pots for another 10 d in the glasshouse. Subsequently, plants were transferred to 1 liter pots and kept in the glasshouse (26 AE 1°C; 16 h : 8 h, light : dark). For the virus-induced gene silencing (VIGS) experiments, plants in 1 liter pots were transferred to a climate chamber under a constant temperature of 26°C, a 16 h : 8 h (light : dark) light regime and 65% relative humidity . Watering, fertilization and light regimes were as previously described Schuman et al., 2014).

Measuring leaf volatiles among 25 natural genotypes and AI-RILs
To measure the emission of leaf volatiles from the 25 natural genotypes, plants in the elongation stage (50 d after germination) were used. Simulated herbivory treatments and leaf volatile collection were conducted as described in detail in a previous study (Zhou et al., 2017). Briefly, the simulated herbivory treatment was carried out at 08:00 h, the leaves were wounded with a pattern wheel, and 20 µl diluted (1 : 5 dilution with water) M. sexta oral secretion was supplied to the wounds. To obtain temporally resolved comparative data for leaf volatiles, small pieces of polydimethylsiloxane tubes (see Kallenbach et al., 2014) were incubated with the damaged leaves from 12:00-16:00 h, 16:00-20:00 h and 20:00-06:00 h (the following morning). To measure the effects of jasmonic acid signaling on the induction of leaf volatiles in the 25 genotypes, plants in the elongation stage (55 d after germination) were treated with methyl jasmonate (MeJA).

New Phytologist
At 08:00 h, 150 µg MeJA dissolved in 40 µl lanolin paste was deposited at the base of one stem leaf of each plant. After 24 h of the MeJA treatment, the leaf was enclosed in a plastic cup (polyethylene terephthalate) together with two polydimethylsiloxane tubes, and the polydimethylsiloxane tubes were collected 24 h later and kept at À20°C until analysis. Three biological plant replicates were sampled for each group.
The measurement of leaf volatiles from the AI-RIL plants (55 d after germination) has been previously described (Zhou et al., 2017). Briefly, simulated herbivory treatment was carried out at 08:00 h, and the leaf was immediately enclosed in a plastic cup where two polydimethylsiloxane tubes were incubated as described above. The polydimethylsiloxane tubes were collected after 24 h (08:00 h the next day), and immediately kept at À20°C until analysis.

QTL mapping
To identify the genetic basis underlying HIPV variations in N. attenuata, we used a forward genetic mapping approach. The quantitative trait locus (QTL) mapping was described in detail in our previous study (Zhou et al., 2017) and was conducted with the Rp a c k a g eQTLREL (Cheng et al., 2011) following the tutorial. In brief, the genotype information of AI-RIL plants and the linkage map were obtained from a dataset reported earlier (Zhou et al., 2017;Guo et al., 2020). In total, we obtained 25 469 single nucleotide polymorphism (SNP) markers among 248 individuals. These markers were grouped into 12 linkage groups. The peak area of herbivore-induced (E)-b-ocimene was log-transformed. Samples with missing genotype or phenotype information were removed. The variance of (E)-b-ocimene within the population was estimated via 'estVC', and then used for the genome-wide scan together with the genotypes. The empirical threshold was estimated based on 500 permutations.
Sequencing and assembling the leaf transcriptome of the Arizona genotype To identify candidate genes that are associated with (E)-bocimene biosynthesis, we sequenced the transcriptomes of both control and M. sexta oral secretion-induced rosette leaves using the Illumina HiSeq 2000 system (San Diego, CA, USA). The growth conditions and simulated herbivory treatments were as described above. Samples from untreated plants were used as control. All samples were collected at 12:00 h. Total RNA was extracted from c. 100 mg of plant tissue using TRIzol (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's protocol. All RNA samples were subsequently treated with DNase-I (Fermentas, Burlington, ON, Canada) to remove all genomic DNA contamination. The mRNA was enriched using an mRNA-seq sample preparation kit (Illumina), and c. 200 bp insertion size libraries were constructed using the Illumina whole transcriptome analysis kit following the manufacturer's standard protocol (HiSeq system; Illumina). All libraries were then sequenced on the Illumina HiSeq 2000 at the sequencing core facility at Max Planck Institute for Molecular Genetics (Berlin, Germany). After removing the adaptor sequences and low-quality reads, in total we obtained 25 297 814 and 19 124 690 paired-end 100 bp reads for control and induced samples, respectively. The raw reads have been deposited in the NCBI SRA repository (SRX1804553 and SRX1804540).
The RNA-seq reads from control and induced samples were pooled for de novo assembly of the leaf transcriptome from the Arizona genotype using the TRINITY tool kit (Grabherr et al., 2011). The minimum contig length was set at 200 bp. In total, we assembled 62 132 transcripts with N50 length equal to 1364 bp.

Measuring abundance of NaTPS25 transcripts in natural accessions
To measure the abundance of NaTPS25 transcripts in the 25 natural accessions, plants in the elongation stage (50 d after germination) were used. Simulated herbivory treatment was carried out on one stem leaf of each plant at 08:00 h. Eight hours after treatment, the leaf was harvested and immediately kept at À80°C until use. Total RNA was isolated from the leaves using TRIzol reagent (Thermo Fisher Scientific) and then treated with DNase-I (Fermentas) following the manufacturers' instructions, to remove contamination from genomic DNA. One microgram of total RNA from each sample was reverse transcribed into cDNA using SuperScript II reverse transcriptase following the manufacturer's instructions (Thermo Fisher Scientific). The relative transcript accumulation of NaTPS25 was measured using reverse transcriptase PCR (RT-PCR) on a Stratagene MX3005P PCR cycler (Stratagene, La Jolla, CA, USA). The elongation factor-1A gene, NaEF 1 a (accession number D63396), was used as the internal standard for normalization (Oh et al., 2013). PCR amplification was performed with Phusion Green High-Fidelity DNA polymerase (Thermo Fisher Scientific) using 40 PCR cycles to enable the detection of rare transcripts of NaTPS25. The primers used are given in Table S2.
Silencing NaTPS25 using VIGS, and measurement of HIPVs from VIGS plants To investigate the function of NaTPS25 in vivo, VIGS based on the tobacco rattle virus was used to transiently knock down expression of the NaTPS25 gene in N. attenuata as previously described . A 163 bp fragment of the open reading frame region of NaTPS25 was amplified by PCR with the primers listed in Table S2. PCR fragments were separated by agarose gel electrophoresis, excised and purified by a gel band purification kit (Amersham Biosciences, Little Chalfont, UK) according to the manufacturer's instructions. Subsequently, they were digested with BamH1 and Sal1 and inserted into plasmid pTV00. After sequencing, the VIGS constructs and pTV00 (empty vector, EV) were separately transformed into Agrobacterium. Additionally, the pTVPDS vector, harboring a part of the sequence of a phytoene desaturase, was used to monitor progress of the gene silencing.
The VIGS experiments were carried out on the Arizona genotype. Inoculation was performed 25 d after germination, and volatile trapping another 20 d later, after establishment of the viral vector. Before any treatment, three leaves (one per plant) were sampled to test silencing efficiency. Afterwards another three leaves (one per plant) were treated with 75 µg MeJA dissolved in 20 µl lanolin paste on the base of leaves at 12:00 h. After 1 h, the outer portion of these leaves was harvested and frozen to test the silencing efficiency, after which the whole plant was again trapped for volatile collection. Lanolin paste (20 µl) was used for the control group.

Cloning and heterologous expression of NaTPS25
To examine the enzymatic function of NaTPS25 in vitro, a truncated version of NaTPS25_AZ (missing the first 13 codons) was cloned from cDNA of Arizona genotype plants using the primer pair TPS25f/r (Table S2). The resulting PCR fragment was cloned into the pET200/D-TOPO vector (Invitrogen) and the Escherichia coli strain BL21 Codon Plus (Invitrogen) was used for heterologous expression. The protein expression and bacterial extract preparation were done as described in a previous study (Zhou et al., 2017). Enzyme assays were performed in a Teflonsealed, screw-capped 1 ml GC glass vial containing 50 ll of the bacterial extract and 50 ll assay buffer with 10 µM geranyl diphosphate, 10 mM MgCl 2 , 0.2 mM Na 2 WO 4 and 0.1 mM NaF. A solid phase micro-extraction (SPME) fiber coated with a 100 lm thick layer of polydimethylsiloxane (Supelco, Belafonte, PA, USA) was placed into the headspace of the vial where it was incubated for 60 min at 30°C, before it was inserted into the injector of the gas chromatograph to analyze the absorbed reaction products. GC-MS analysis was conducted using an Agilent 6890 Series gas chromatograph coupled to an Agilent 5973 quadrupole mass-selective detector. The gas chromatograph was operated with a DB-5MS column (30 m 9 0.25 mm 9 0.25 µm; Agilent, Santa Clara, CA, USA).
Targeted resequencing of the cDNA and genomic regions of NaTPS38 and confirmation of the inversion of NaTPS25 in Utah genotype To sequence the NaTPS38 locus from all the accessions and characterize its allelic variations, we amplified the full-length open reading frame using high-fidelity Phusion DNA Polymerase (Thermo Fisher Scientific) and the primer combination TPS38f and TPS38r (Table S1). The following cycle parameters were used: 98°Cf o r 30 s; 35 cycles of 98°C for 10 s, 52°Cf o r2 0sa n d7 2°Cf o r 1 . 5m i n ;a n dafi n a ls t e po f7 2°C for 5 min. PCR products were separated by agarose gel electrophoresis, excised and purified with a PCR clean-up gel extraction kit (Macherey-Nagel, D€ uren, Germany), then cloned into the pJET 1.2/blunt vector (Thermo Fisher Scientific) and transformed into E. coli (TOP10 Competent Cells; Invitrogen) according to the manufacturers' instructions. Plasmid extraction was performed using the NucleoSpin Plasmid Kit (Macherey-Nagel). Sequencing reactions were performed using the BigDye Terminator mix (Thermo Fisher Scientific). The obtained sequences were manually corrected, and data analysis was executed with GENEIOUS 6.0.5. The cDNA sequences have been deposited in the GenBank Nucleotide Sequence Database (accession nos. MN400683-MN400708).
To examine the genetic architecture of NaTPS25, we isolated genomic DNA from 10-d-old seedlings using a cetrimonium bromide (CTAB) extraction buffer (Bubner et al., 2004). Genomic DNA purity and integrity were evaluated by amplifying the elongation factor-1A with the primer set shown in Table S2. Sequencing approaches for the comparison between Arizona genotype and traced fragments in Utah genotype were conducted as described above. PCR fragment amplifications were carried out using Phusion Green DNA Polymerase (Thermo Fisher Scientific) in 20 µl PCRs (91 Phusion Green HF Buffer, 10 mM dNTPs, 3% DMSO and 5 mM primers). The PCR program (30 s at 98°C, 35 cycles of 10 s at 98°C, 30 s at 60°C and 2 min at 72°C, with a final elongation of 5 min at 72°C) varied in annealing temperature (72°C) according to the primer set (Table S2). We validated the presence of the intact NaTPS25 allele with primer set I (PrS I), and the presence of exon 6-7 with PrS II. Moreover, we utilized two primer combinations to provide genetic information for the inversion event of exons 6 and 7. PrS II was used to detect exons 6 and 7, and primer PrS III to amplify a fragment specific to this inversion.

Measuring the abundance of NaTPS38 transcripts in hybrid plants
To characterize whether the variations of transcript abundance of NaTPS38 were due to changes in cis-o rtrans-regulatory elements, we measured the abundance of NaTPS38 transcripts in the two parental inbred lines and their hybrid lines in F 1 . Plants in their elongation stage (50 d after germination) were used. Simulated herbivory treatment was carried out at 08:00 h on one stem leaf of each plant. Eight hours after treatment, leaves were harvested and flash-frozen in liquid nitrogen, then kept at À80°C until processing. Total RNA isolation, cDNA synthesis and quantitative PCR (qPCR) amplification were carried out as described above. The primers used to specifically amplify the alleles are shown in the Supporting Information (Table S2).

Statistical analysis
All statistical analyses were performed in R 3.5.1 (R Core Team, 2016). The functions 'prcomp', 't.test' and 'anova' from the STATS package were used for principal components analysis, Student's t-tests and ANOVA analysis, respectively. The raw data and data analysis R scripts are deposited in the figshare repository: https://doi.org/10.6084/m9.figshare.11848137.v1.

Results
(E)-b-ocimene and (E)-a-bergamotene contributed to associated changes in HIPVs and floral volatiles among N. attenuata natural accessions We first investigated HIPV variations in N. attenuata by sampling the herbivory-induced leaf headspace from 25 genotypes that were New Phytologist (2020) 228: 1083-1096 Ó2020 The Authors New Phytologist Ó2020 New Phytologist Trust www.newphytologist.com

Research
New Phytologist collected from different habitats within the species range. We found two monoterpenes, (E)-b-ocimene and linalool, and one sesquiterpene, (E)-a-bergamotene, to be highly variably released among the studied genotypes ( Fig. 1a-c). In addition, previous studies have shown that the emission of all three compounds in leaves was significantly induced by herbivory (Halitschke et al., 2000), supporting their classification as HIPVs. Jasmonate signaling is involved in regulating the emission of many HIPVs in N. attenuata, and its accumulation is variable among wild genotypes (Schuman et al., 2009). To examine whether the observed HIPV variations were attributable to differences in herbivory-induced jasmonate signaling, we supplied a similar level of MeJA to the leaves of these genotypes. Under similarly elevated levels of jasmonate signaling in MeJA-treated plants, the emission of the three terpenoids remained variable and were the major contributors to the variation in HIPVs among the investigated genotypes (Fig. 1d). This suggests that intraspecific variations of HIPVs are largely independent of jasmonate signaling.
In addition to a flower-specific volatile, benzoyl acetone (Guo et al., 2020), the three terpenoids in HIPVs, (E)-bocimene, linalool and (E)-a-bergamotene, were also found in N. attenuata floral volatiles (Kessler & Baldwin, 2007). We reasoned that if the same compound in HIPVs and floral volatiles was regulated by the same genetic mechanism, its abundance would show correlated changes in amounts of HIPVs and floral volatiles among different genotypes. To examine this, we further measured the floral volatiles of the 25 genotypes and compared the emissions of (E)-b-ocimene, linalool and (E)-a-bergamotene emission in the HIPVs and floral volatiles. We found that linalool, which was found in both floral volatiles and HIPVs, showed no correlation among the different genotypes, consistent with a recent study (He et al., 2019). However, the variations of (E)-a-bergamotene and (E)-b-ocimene emission in HIPVs and floral volatiles were significantly correlated ( Fig. 1;   Leaves (peak area × 10 6 ) Flowers (peak area * 10 6 ) (E)-α-bergamotene 57  Changes in NaTPS25 and NaTPS38 expression are associated with intraspecific variations of (E)-b-ocimene and (E)a-bergamotene emissions, respectively, in both HIPVs and floral volatiles To characterize the genetic mechanisms underlying the correlated changes of (E)-a-bergamotene and (E)-b-ocimene emission in HIPVs and floral volatiles, we used both forward and reverse genetic tools. We have previously established an AI-RIL population that was derived from crossing the Arizona and Utah genotypes whose (E)-b-ocimene and (E)-a-bergamotene emissions differed in HIPVs and floral volatiles. Using this population, we identified a terpene synthase, NaTPS38, on linkage group 2 that controls the emissions of (E)-a-bergamotene in both HIPVs and floral volatiles (Zhou et al., 2017). Among different N. attenuata genotypes, the changes in NaTPS38 expression were correlated with (E)-a-bergamotene emissions in both HIPVs and floral volatiles (Zhou et al., 2017).
Here, using a similar approach, we found that variations in (E)-b-ocimene emission in HIPVs were mapped to the same locus as NaTPS38 on linkage group 2 of N. attenuata ( Fig. 2a; Zhou et al., 2017). This location overlaps with a terpene synthase cluster, which contains a subclade of TPS-b monoterpene synthase genes, including TPS25, TPS27 and TPS38 (Falara et al., 2011). Among these genes, TPS25 and TPS27 in tomato were highly similar to the grape VvTPS47 that produces (E)-b-ocimene. However, after searching the draft genome of N. attenuata (created using the Utah genotype), we did not find any genes that showed high similarity to either TPS25 or TPS27. Sequencing the transcriptome of the Arizona genotype allowed us to identify NaTPS25, which showed high similarity to TPS25 and TPS27 in tomato.
To examine the enzymatic function of NaTPS25, we measured the TPS activity of NaTPS25 in vitro by heterologous expression of the gene in E. coli. The recombinant NaTPS25 protein converts the substrate geranyl diphosphate into (E)-b-ocimene as its only product (Fig. 2b). To examine the function of NaTPS25 in vivo, we silenced the expression of NaTPS25 in N. attenuata (Arizona genotype) using VIGS. When NaTPS25 expression was silenced, levels of both constitutive and MeJA-induced (E)-bocimene emission in HIPVs were strongly reduced (Fig. 2c,d). Together, these results suggest that herbivore-induced (E)-bocimene is associated with the expression of NaTPS25 in HIPVs.
We further investigated whether the expression changes in NaTPS25 are associated with (E)-b-ocimene in HIPVs and floral volatiles. We first measured the presence of NaTPS25 transcripts among the 25 genotypes by RT-PCR. Interestingly, we observed presence and absence variations of NaTPS25 at the level of transcripts, which correlate with the emission of (E)-b-ocimene from induced leaves among genotypes (Fig. 3a). While NaTPS25 transcripts can be found in all 15 genotypes that emitted (E)-bocimene after simulated herbivory, they were not found in any of the other 10 genotypes that did not emit detectable levels of (E)-  (Fig. 3a). We then used qPCR to measure the transcript abundance of NaTPS25 in flowers in order to compare them to emissions of (E)-b-ocimene in the respective floral volatile profiles of each genotype. We found that among different genotypes, emissions of (E)-b-ocimene in floral volatiles were significantly correlated with NaTPS25 transcript abundance (P = 0.013, Fig. 3b). Interestingly, while most of the genotypes that had a high abundance of NaTPS25 transcripts also had a high floral (E)-b-ocimene emission, a few genotypes (e.g. genotype Utah, 43 and 351) showed very low levels of NaTPS25 transcript abundance, yet still emitted high levels of (E)-b-ocimene in their flowers (Fig. 3b). These results suggest that in addition to NaTPS25, other genes are involved in floral (E)-b-ocimene emission. Consistently, in the N. attenuata genome (Utah genotype), there is another member of the TPS-b subfamily, of which several homologs from Vitis vinifera and Arabidopsis thaliana can synthesize (E)-b-ocimene in vitro using geranyl diphosphate as the substrate (Huang et al., 2010;Martin et al., 2010). Furthermore, one homolog of TPS19 (NIATv7_g29416), which is a member of TPS-e/f subfamily that can produce (E)-b-ocimene using neryl pyrophosphate substrate (Matsuba et al., 2013), was found highly to be expressed in the flowers of N. attenuata (http://nadh.ice. mpg.de/NaDH-eFP/cgi-bin/efpWeb.cgi; Brockm€ oller et al., 2017).
Together with our previous study (Zhou et al., 2017), these results suggest that expression changes in NaTPS25 and NaTPS38 are associated with intraspecific variations of (E)-bocimene and (E)-a-bergamotene emissions, respectively, in HIPVs and floral volatiles.

Two distinct genetic mechanisms resulted in changes in allelic expression of NaTPS25 and NaTPS38
To investigate the genetic mechanisms underlying the allelic expression differences in NaTPS25 and NaTPS38, we analyzed sequence polymorphisms of NaTPS25 and NaTPS38 among different genotypes.
Consistent with its presence and absence of detectable transcripts levels, we found a structural variation in the NaTPS25 gene among different genotypes. For genotypes that emit (E)-bocimene and express NaTPS25 after simulated herbivory or MeJA treatment, the complete NaTPS25 gene can be amplified from genomic DNA with primers that bind to regions containing the start codon and stop codon (Fig. 4a,b). There was no obvious gene size variation among these genotypes (Fig. 4b). However, for the accessions that failed to emit (E)-b-ocimene, NaTPS25 cannot be amplified using the same primers. Designing exonspecific primers and comparing the genomic sequences of NaTPS25 revealed that the Utah genotype contains two fragments of NaTPS25, disrupted by a c. 9 kb long terminal repeat (LTR)-gypsy transposon. The first fragment contains only exons 6 and 7, whereas the second fragment contains partial sequences of exons 2 and 7 (Fig. 4a). Using two primer pairs that can specifically amplify these two fragments among different genotypes, we found that most of the genotypes that failed to emit (E)-bocimene after simulated herbivory shared the same disrupted gene structure (Fig. 4c). Furthermore, additional genomic changes to the NaTPS25 fragments were found in some genotypes, such as genotype 421, in which the joined fragments    New Phytologist (2020) 228: 1083-1096 Ó2020 The Authors New Phytologist Ó2020 New Phytologist Trust www.newphytologist.com

Research
New Phytologist between exon 2 and exon 7 could not be amplified, and genotypes 176, 341 and 384, in which the exon 6 and 7 fragments could not be amplified. These results showed that allelic differences of NaTPS25, caused by LTR-gypsy transposon insertions, contributed to the variations of (E)-b-ocimene emission in the HIPVs among different genotypes.
NaTPS38 transcripts were found in all genotypes but differed in their abundances (Zhou et al., 2017). To test whether polymorphic sites at a cis-regulatory region and/or changes in transregulatory elements contributed to the observed variations in the expression of NaTPS38, we measured the abundance of NaTPS38 transcripts in the plants from Arizona and Utah genotypes and their F 1 progeny. Because all the trans factors are shared amongst the F 1 progeny, differences in allelic transcript abundance between parents and F 1 indicate changes in relative transregulatory activity, whereas differences in relative transcript abundance between two alleles in F 1 indicate changes in relative cisregulatory activity (Wittkopp, 2013). After simulated herbivore elicitation, while transcript levels of the NaTPS38 Arizona allele did not differ between parental and F 1 plants (Fig. 5, P > 0.1), levels of the NaTPS38 Utah allele were significantly lower in F 1 plants than in the Utah genotype (P = 0.0079, Student's t-test), suggesting that the Arizona genotype carries a trans-suppressor of NaTPS38. In addition, the relative transcript abundances of Arizona and Utah alleles in F 1 plants were significantly different (Fig. 5,P = 0.012), indicating that the cisregulatory activity of NaTPS38 is different between Arizona and Utah genotypes. We further sequenced the upstream regulatory region of NaTPS38 in both Arizona and Utah genotypes. In comparison to the Utah genotype, the Arizona genotype has a 129 bp insertion and one SNP (change from A to G) within the 1 kb upstream region of NaTPS38 (Fig. S1). In summary, polymorphisms at both trans-and cis-regulatory regions contribute to allelic differences of herbivory-induced NaTPS38 expression and therefore (E)-a-bergamotene emission in N. attenuata.
Allelic variations of NaTPS25 and NaTPS38 resulted in correlated changes of (E)-b-ocimene and (E)-a-bergamotene emission in HIPVs and floral volatiles in N. attenuata We further explored whether the correlated changes in (E)-bocimene and (E)-a-bergamotene emission in HIPVs and floral volatiles are due to allelic variations of NaTPS25 and NaTPS38.
To this end, we first sequenced the full-length open reading frame of NaTPS38 from the 25 natural genotypes. None of these contained either a premature stop codon or an insertion/deletion (Fig. 6a). The NaTPS38 transcripts could be classified into three different haplotypes (YGI, HGV and YEI) that resulted from three common polymorphic nonsynonymous substitution sites: Y32H, G280E and I450V. Mapping these amino acid changes to the protein model of NaTPS38 suggested that Y32H and G280E are not involved in the formation of the active site cavity, although I450V is part of a loop structure near the entrance of the active site cavity (Fig. S2). However, the substitution I450V was found both in genotypes that emit high (e.g. genotype 304) and low levels (e.g. genotype 85) of herbivory-induced (E)-abergamotene, indicating that this amino acid change does not affect the enzymatic function of the NaTPS38 protein (Fig. 6a). We then compared (E)-a-bergamotene emission among natural genotypes that carry different haplotypes. The genotypes that carry the YEI haplotype showed lower (E)-a-bergamotene emissions from both flowers and herbivory-induced leaves than those that harbor the YGI haplotype (Fig. 6b). Consistently, abundance of NaTPS38 transcripts from the YEI haplotype carriers is also lower than those with the YGI haplotype (Fig. S3). Interestingly, the HGV haplotype showed intermediate (E)-a-bergamotene emission in both flowers and herbivory-induced leaves (Fig. 6b). While the floral transcript abundance of NaTPS38 in HGV haplotype carriers is similar to that in YEI haplotype carriers, in terms of herbivory-induced leaf abundances, they are similar to YGI haplotype carriers (Fig. S3).
Because NaTPS38 and NaTPS25 are co-localized in the genome, we hypothesized that genetic variation of NaTPS25 is also linked to the three different haplotypes of NaTPS38. Mapping NaTPS25 presence/absence information in herbivory-induced leaves to the three haplotypes of NaTPS38 showed that 80% (eight out of 10) of individuals that carried the NaTPS38-YGI haplotype also had no expression of NaTPS25, whereas individuals that carried the NaTPS38-HGV and NaTPS38-YEI haplotypes all had detectable NaTPS25 transcript levels (Fig. 6a). Consistently, in both flowers and herbivory-induced leaves, Green arrows refer to primers that can specifically amplify the Utah allele and blue arrows refer to primers that can specifically amplify the Arizona allele. Primers (gray arrows) that do not distinguish the two alleles were also used to measure total transcript abundance of NaTPS38. (b) Expression differences of NaTPS38 alleles between F 0 and F 1 plants. Plants were treated with wounding plus oral secretion (n = 5). Mean and SE are shown. Student's t-test was used to determine differences between F 0 and F 1 plants. ns, no significant difference was found. (c) Expression differences of alleles in F 1 individuals. Plants were treated with wounding plus oral secretion (n = 5). Each dot represents an individual plant. Difference between two alleles was determined by a paired Student's t-test. Blue and green colors indicate alleles from Arizona and Utah genotypes, respectively.

Ó2020 The Authors
New Phytologist Ó2020 New Phytologist Trust New Phytologist (2020) 228: 1083-1096 www.newphytologist.com individuals that carried the NaTPS38-YGI haplotype also accumulated only very low levels of NaTPS25 transcripts and produced a lower amount of (E)-b-ocimene, whereas individuals that carried the NaTPS38-YEI haplotype accumulated high levels of NaTPS25 transcripts and emitted high levels of (E)-b-ocimene (Figs 3a, 6b, S3c). Again, the NaTPS38-HGV haplotype plants showed intermediate NaTPS25 transcript abundance but high levels of (E)-b-ocimene emission in both flowers and leaves (Figs 3a,6b,S3c). These results reveal that the linked allelic differences of NaTPS38 and NaTPS25 in N. attenuata contributed to correlated changes in both floral volatiles and HIPVs (Fig. 7).

Discussion
Intrinsic genetic linkage and pleiotropy have been proposed to explain correlated variation in floral and defense traits (Johnson  Fig. 6 Allelic variations of the terpene cluster are associated with differences in floral volatiles and herbivory-induced plant volatiles (HIPVs) among Nicotiana attenuata natural accessions. (a) Three different haplotypes found from the open reading frame of NaTPS38 from N. attenuata accessions. The amino acid sequence of NaTPS38 from the Utah genotype is set as the reference, and the different haplotype loci are marked in red. The presence of an intact NaTPS25 gene in the natural accessions is shown on the right (+, present; À, not present). (b) The emission of (E)-a-bergamotene and (E)-b-ocimene in both flowers and herbivory-induced leaves are different among three NaTPS38 haplotypes. The data for (E)-a-bergamotene emission were extracted from our previous study (Zhou et al., 2017). Mean and SE are shown. Statistical differences among haplotypes were determined by ANOVA (P < 0.05).   et al., 2015). However, direct genetic evidence for this is scarce. Here, by exploring natural variation and using both forward and reverse genetics tools, we demonstrate that allelic differences in two clustered terpene synthase genes, NaTPS25 and NaTPS38, result in correlated intraspecific variations of both HIPVs and floral volatiles in N. attenuata.
We found that variations in the emissions of two terpenoids, (E)-b-ocimene and (E)-a-bergamotene, showed correlated changes in HIPVs and floral volatile profiles (Fig. 1). (E)-bocimene is often found among HIPVs and floral headspace in different species (Farre-Armengol et al., 2017). While (E)-bocimene in HIPVs can act as an important chemical cue for attracting natural enemies of phytophagous insects among plant species, it did not increase the attraction of predators of M. sexta larvae in N. attenuata in nature (Kessler & Baldwin, 2001). Floral (E)-b-ocimene emission is also common among flowering plants and can act as an attractant to generalist pollinators, such as bees and beetles (Farre-Armengol et al., 2017). Although N. attenuata often attracts M. sexta as its pollinator, in some populations (e.g. Arizona), bees are also found to visit N. attenuata flowers (D. Kessler, personal communication). Therefore, it is plausible that floral (E)-b-ocimene emission in N. attenuata might increase beemediated pollination success. Different from (E)-b-ocimene, (E)a-bergamotene emission has only been reported in a few species, such as maize and N. attenuata. In both maize and N. attenuata, (E)-a-bergamotene in HIPVs is important for indirect defense function (Kessler & Baldwin, 2001;Degen et al., 2004). Although it remains unclear whether (E)-a-bergamotene is also emitted in floral tissues in maize, the present study showed that the flowers of N. attenuata emit (E)-a-bergamotene, which increases M. sexta moth-mediated pollination success (Zhou et al., 2017). Therefore, the correlated changes of (E)-b-ocimene and (E)-a-bergamotene in HIPVs and floral volatiles can result in concerted adaptation to pollinators and herbivores in N. attenuata.
The changes of (E)-b-ocimene and (E)-a-bergamotene emissions in HIPVs and floral headspace are probably due to allelic variations of two clustered terpene synthase genes, NaTPS25 and NaTPS38. This is supported by four independent lines of evidence. First, forward genetic mapping showed that the same locus is associated with the emissions of (E)-b-ocimene and (E)a-bergamotene in HIPVs. Although it was difficult to find the responsible gene in the QTL directly from the fragmented N. attenuata genome assembly, synteny information between N. attenuata and tomato showed that the QTL probably contains a TPS cluster that includes TPS25 and TPS38. Second, using reverse genetic approaches, we demonstrated that the expression of NaTPS25 and NaTPS38 were required for (E)-b-ocimene and (E)-a-bergamotene emission in HIPVs ( Fig. 2; Zhou et al., 2017), respectively. Among different natural accessions, the expression levels of NaTPS25 and NaTPS38 were furthermore correlated with (E)-b-ocimene and (E)-a-bergamotene emission, respectively, in both floral headspace and HIPVs ( Fig. 2b; Zhou et al., 2017). Third, in vitro enzymatic assays demonstrated the specific terpene synthase functions of both enzymes ( Fig. 2; Zhou et al., 2017). Fourth, allelic variation of NaTPS38 and NaTPS25 were tightly linked and formed three main haplotypes. Individuals carrying different haplotypes showed correlated changes in (E)-b-ocimene and (E)-a-bergamotene emissions in floral headspace and HIPVs (Fig 6).
While the correlated changes of (E)-a-bergamotene in HIPVs and floral volatiles are due to pleiotropy of NaTPS38, the correlated changes of (E)-b-ocimene in HIPVs and floral volatiles are probably due to genetic linkage. This is because some genotypes do not emit (E)-b-ocimene in HIPVs due to loss-of-function mutations in NaTPS25, but they still emit (E)-b-ocimene in flowers, although at relatively low levels. A plausible scenario is that a transcriptional regulator of another TPS that is involved in floral (E)-b-ocimene emissions is in strong genetic linkage with the NaTPS25/38 locus. If this is the case, the correlation between emissions of (E)-b-ocimene in the floral headspace and in HIPVs might disappear when the genetic linkage is interrupted, for example among genotypes that harbor much more frequent recombination than we have analyzed or among closely related species.
We found two different genetic mechanisms are associated with allelic expression changes in NaTPS25 and NaTPS38. While loss-of-function mutations mediated by a transposable element insertion resulted in expression variations of NaTPS25, changes at both cis-andtrans-regulatory elements contributed to the different expression of NaTPS38 among genotypes. The differences in underlying genetic mechanisms might be due to differences in fitness optima of (E)-a-bergamotene (one optimum peak, thus under stabilizing selection) and (E)-b-ocimene (two optimum peaks, thus under disruptive selection). However, it is currently unclear which stresses may favor loss-of-function in NaTPS25. Another explanation is that NaTPS25 and NaTPS38 have different pleiotropic functions. As silencing NaTPS38 resulted in the complete absence of (E)-a-bergamotene in the HIPVs and floral volatile profiles (Zhou et al., 2017), NaTPS38 is probably the only key enzyme that is responsible for the biosynthesis of (E)-a-bergamotene in both HIPVs and floral volatiles in N. attenuata. However, while silencing NaTPS25 resulted in strong decreases in (E)-b-ocimene emission in HIPVs, genotypes (such as Utah) with loss-of-function in NaTPS25 still emitted (E)-b-ocimene from flowers. This suggests that the biosynthesis of (E)-b-ocimene in flowers might be a result of several genes with similar functions in N. attenuata. Therefore, the level of functional pleiotropy for NaTPS25 is lower than for NaTPS38. As a consequence, a loss-of-function mutation in NaTP25 might be associated with less distinct fitness costs than for one in NaTPS38. This is consistent with the pattern observed from the genetic mechanisms underlying changes of floral color and scent (Hoballah et al., 2007;Wessinger & Rausher, 2012;Xu et al., 2012;Sas et al., 2016;Sheehan et al., 2016), in which loss-of-function mutations are more common than regulatory element changes in genes with less pleiotropy in their functions.
The intrinsic link between floral volatiles and HIPVs in natural populations of N. attenuata implies that intraspecific variation of plant volatiles can be maintained due to pleiotropy and/or genetic linkage. For example, when an N. attenuata population is under directional selection for high amounts of floral (E)-a-Ó2020 The Authors New Phytologist Ó2020 New Phytologist Trust New Phytologist (2020) 228: 1083-1096 www.newphytologist.com bergamotene emission, either imposed by M. sexta moths or hummingbirds (both show strong preference for (E)-a-bergamotene (see Kessler & Baldwin, 2007)), individuals that carry the YGI haplotype would have higher fitness and the frequency of the YGI haplotype in the population would increase. As a consequence, the level of herbivory-induced (E)-a-bergamotene and (E)-b-ocimene emissions would increase in the population. In such a scenario, changes in HIPVs in these populations are not due to selection on HIPVs, but rather are due to the pleiotropy or genetic linkage of the YGI haplotype, which is also associated with high levels of herbivory-induced (E)-a-bergamotene and low levels of (E)-b-ocimene emission. The same logic may also apply to how selection on HIPVs could result in changes in floral volatiles.
The presence of genetic linkage and pleiotropy between HIPVs and floral volatiles suggest that plant volatiles might evolve under diffuse selection imposed from herbivores and pollinators in N. attenuata (Iwao & Rausher, 1997;Strauss & Irwin, 2004;Agrawal & Stinchcombe, 2009;Wise & Rausher, 2013). Diffuse coevolution can be common for plant resistance to multiple-herbivore communities (Wise & Rausher, 2013) and recent genetic analyses have discovered correlated changes of chemical defense compounds among tissues (Chan et al., 2011;Keith & Mitchell-Olds, 2019) which were due to genetic linkage or pleiotropy, similar to the two terpenoids we found here. It is plausible that many other metabolic traits also show correlated changes among tissues, although future studies using comparative metabolomics and transcriptomics approaches are needed to further examine this. Diffuse interactions with multiple-herbivore communities and correlated changes of chemical defenses among tissues often constrain the evolution of plant resistance (Wise & Rausher, 2013;Keith & Mitchell-Olds, 2019). Because plant-herbivore and plant-pollinator interactions are often studied in isolation (Lucas-Barbosa, 2016), the extent to which diffuse interactions exist within plant-herbivore and plant-pollinator interactions, and their evolutionary consequences, remain largely unclear. One recent study using an experimental evolution approach has demonstrated a diffuse coevolution of floral signals and plant direct defenses in Brassica rapa (Ramos & Schiestl, 2019), although this study did not investigate genetic correlation of floral signals and plant defenses. The studies and results reported here highlight the importance of investigating trait evolution in a community context (Strauss & Irwin, 2004). Prospective studies deploying experimental evolution approaches, in which the multigenerational fitness of different genotypes is determined under manipulated selection agents in their native habitat, will be needed to provide complementary understanding on how adaptive traits evolve in nature.

Supporting Information
Additional Supporting Information may be found online in the Supporting Information section at the end of the article.   Table S1 GPS coordinates of the location where the 25 Nicotiana attenuata natural accessions were originally collected.

Table S2
The primer sets used in this study.
Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.
New Phytologist is an electronic (online-only) journal owned by the New Phytologist Trust, a not-for-profit organization dedicated to the promotion of plant science, facilitating projects from symposia to free access for our Tansley reviews and Tansley insights. Regular papers, Letters, Research reviews, Rapid reports and both Modelling/Theory and Methods papers are encouraged. We are committed to rapid processing, from online submission through to publication 'as ready' via Early View -our average time to decision is <26 days. There are no page or colour charges and a PDF version will be provided for each article.
The journal is available online at Wiley Online Library. Visit www.newphytologist.com to search the articles and register for table of contents email alerts.
If you have any questions, do get in touch with Central Office (np-centraloffice@lancaster.ac.uk) or, if it is more convenient, our USA Office (np-usaoffice@lancaster.ac.uk) For submission instructions, subscription and all the latest information visit www.newphytologist.com New Phytologist (2020) 228: 1083-1096 Ó2020 The Authors New Phytologist Ó2020 New Phytologist Trust www.newphytologist.com