Transcriptomic Profiling Reveals Shared Signalling Networks Between Flower Development and Herbivory-Induced Responses in Tomato

Most flowering plants must defend themselves against herbivores for survival and attract pollinators for reproduction. Although traits involved in plant defence and pollinator attraction are often localised in leaves and flowers, respectively, they will show a diffuse evolution if they share the same molecular machinery and regulatory networks. We performed RNA-sequencing to characterise and compare transcriptomic changes involved in herbivory-induced defences and flower development, in tomato leaves and flowers, respectively. We found that both the herbivory-induced responses and flower development involved alterations in jasmonic acid signalling, suppression of primary metabolism and reprogramming of secondary metabolism. We identified 411 genes that were involved in both processes, a number significantly higher than expected by chance. Genetic manipulation of key regulators of induced defences also led to the expression changes in the same genes in both leaves and flowers. Targeted metabolomic analysis showed that among closely related tomato species, jasmonic acid and α-tomatine are correlated in flower buds and herbivory-induced leaves. These findings suggest that herbivory-induced responses and flower development share a common molecular machinery and likely have coevolved in nature.


INTRODUCTION
Most plants have to protect themselves from insect herbivory for survival and attract pollinators for reproduction. While the defence against insect herbivores is often mediated by toxic chemical compounds and physical barriers in leaves (Wu and Baldwin, 2010), pollinator attraction is mostly achieved by various floral signals, such as colours and scents produced by flowers (Fenster et al., 2004;Sheehan et al., 2012;Fattorini and Glover, 2020). The spatial separation of tissues that mediate plant defence and pollinator attraction resulted in often isolated studies on the molecular mechanisms and evolution of traits involved in defence and floral signals, but see examples in Strauss and Irwin (2004); Pashalidou et al. (2020). However, if plant defence and floral signals are regulated by the same molecular machinery, they will show diffuse coevolution due to pleiotropy (Iwao and Rausher, 1997;Strauss and Irwin, 2004;Agrawal and Stinchcombe, 2009;Wise and Rausher, 2013). For example, in Cardamine hirsuta, allelic variations at the major floral repressor Flowering Locus C affect both flowering time and plant defences, such as total glucosinolates and stress-related phytohormones in leaves (Rasmann et al., 2018). Therefore, in a hypothetic population, the flowering time of C. hirsuta will evolve under joint selection imposed by both pollinators and herbivores. The same is also true for defences. Indeed, recent experimental evolution studies in Brassica rapa found that floral signals and plant defence evolved via diffuse coevolution, likely due to pleiotropic effects (Zu et al., 2016;Ramos and Schiestl, 2019). For example, plants under selection by bee pollinators evolved increased floral attractiveness, but the process was compromised by the presence of herbivores. However, the underlying molecular mechanisms remain unclear.
Interestingly, similar molecular and metabolic processes are involved in herbivory-induced defence responses, as well as during flower development and opening. JA signalling is required for flower development Yuan and Zhang, 2015;Huang et al., 2017). For example, in Arabidopsis thaliana, the JA-deficient dad1 mutant shows slower flower development, delayed anther dehiscence and reduced viability of pollen grains (Ishiguro et al., 2001). In wild tobacco, the floral limbs cannot fully open when JA biosynthesis is interrupted (Stitz et al., 2014). A similar pattern is also found in Chinese cabbage (Brassica campestris; Peng et al., 2019) and rice (Xiao et al., 2014). Recently, a study in tomato suggests that the JA-mediated flower development and opening are mediated by SlMYB21 (Niwa et al., 2018). Flower development and opening also involve serious changes in the primary and secondary metabolisms. During flower development, the chromoplasts gradually increase (Giuliano et al., 1993) resulting in a progressively elevated biosynthesis of carotenoids and other coloured substances to produce colourful flowers. In addition, flower opening coincides with petal expansion, which involves hydrolysis of stored carbohydrates (Bieleski, 1993;Bieleski et al., 2000;Vergauwen et al., 2000;Yap et al., 2008), and release of volatile organic compounds from petals and nectar (González-teuber and Heil, 2009;Heil, 2011;Dudareva et al., 2013).
Recent studies have shown that herbivory-induced defence responses and flower development not only involve a similar phytohormonal regulation and metabolic rearrangement, but also that they can even share the same molecular machinery. For instance, CYP82G1, which is expressed in both herbivoryinduced leaves and flowers, is required for the biosynthesis of a homoterpene emitted from both respective tissues in Arabidopsis (Lee et al., 2010). Similarly, in wild tobacco, tissue-specific expression of NaTPS38 is required for the emission of (E)-αbergamotene from leaves after herbivory and from flowers Xu et al., 2020). Likewise, changes in diterpene glycoside biosynthesis alter both the anti-herbivore defence and flower development in the wild tobacco (Li et al., 2021). In tomato and potato, altering the biosynthesis of steroidal glycoalkaloids (SGA), one of the major chemical defences in Solanum suppresses flower development (Itkin et al., 2011;Umemoto et al., 2016). However, it remains unclear to what extent herbivory-induced responses and flower development share the same molecular machinery at the genome-wide level.
Here, we aim to address this challenge using comparative transcriptomic profiling. To this end, we sequenced the transcriptomes of control and tobacco hornworm treated leaves, as well as flower tissues at two different developmental stages in three genotypes of domesticated tomato (WT, JA biosynthesis mutant spr8 and the transgenic Prosystemin (PS) over-expression line 35S::PS). The results revealed that a large number of genes, higher than expected by chance, was shared between herbivoryinduced responses and flower development. In addition, targeted metabolomic analysis showed correlated changes of JA or α-tomatine in flower buds and herbivory-induced leaves among six closely related species, suggesting diffuse coevolution due to gene pleiotropy in tomatoes.

Plant Materials and RNA Isolation
Tomato (Solanum lycopersicum L.) cv. Castlemart (CM) was used as the wild type (WT). The mutant plants 35S::PS and spr8 used in this study were previously described (Yan et al., 2013). Six wild tomato species used in this study included S. arcanum, S. chilense, S. chmielewskii, S. corneliomulleri, S. habrochaites and S. neorickii. Plants were grown from seeds in an insect-free glasshouse at ETHZ (Lindau-Eschikon, canton Zurich, Switzerland). They were potted in 1-L pots using fresh soil (Ricoter Substrate 214, Ricoter Erdaufbereitung AG, Aarberg, Switzerland) 1 and fertiliser granules (Gartensegen, Hauert HBG Dünger AG, Grossaffoltern, Switzerland). 2 Plants were watered two-four times per week. All plants were grown in a greenhouse with a temperature of 20°C at daytime and 16°C at night, with 12 h light at 18 kLux and 50% relative humidity.
The leaf samples were collected 6 weeks after germination. Herbivory treatments were performed using Manduca sexta, which is a specialist herbivore of Solanaceae. Manduca sexta eggs were obtained from an in-house colony of the Max Planck Institute for Chemical Ecology, Jena, Germany. One M. sexta neonate was placed on each of four selected leaves per plant. After ~24 h, three well-damaged leaves were selected per plant, the damage level on each of those leaves was recorded and they were separately collected and immediately frozen in liquid nitrogen. On average, ~24 h herbivory resulted in ~10% leaf damage. Comparable leaves from untreated plants were used as control. Three replicates (plants) were analysed for each genotype and treatment.
Flower sampling was conducted 11-12 weeks after the herbivory sampling (when most of the plants started flowering). To synchronise flowering time, we transferred all plants to 3-L pots and shifted to a fully controlled climate chamber with a temperature of 22°C at daytime and 18°C at night, with 10 h light 4-5 weeks before sampling. To avoid the confounding effects of herbivory on flower development, we only collected flower samples from the control plants (plants that were not treated with herbivory). Different flower stages of the same plants were collected on the same day (all between 13:00-16:00) and frozen in liquid nitrogen. The two developmental stages were identified according to a previous study (Dobritzsch et al., 2015). The fully opened flowers were considered as stage 5, whereas the flower bud located next the fully opened flowers on the same inflorescence was considered as stage 4.
Three replicates were collected for each genotype per developmental stage. In total, 18 flower samples (3 genotypes * 3 replicates * 2 developmental stages) were collected. Both leaf and flower samples were stored in a −80°C freezer until RNA extraction.
The mapping rate ranged from 90.5 to 95.2% of each library (on average, 93.2%). Following mapping, StringTie (v1.0.4; Pertea et al., 2015) was used to reconstruct and assemble the transcripts. The expression abundance of each gene was calculated using RSEM (v1.2.12; Li and Dewey, 2011) and presented as transcripts per million reads (TPM). The raw sequencing data of 36 samples were submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) 3 under accession number PRJNA600385.

Differential Gene Expression Analysis
The R package edgeR (v3.26.8; Robinson et al., 2010) was used to identify differentially expressed genes (DEGs). Raw count data for each pair of comparison groups were used as input. Genes with counts greater than 10 in at least three samples were considered in the downstream analysis, while all others were defined as lowly expressed genes and removed. The criterion for being considered a DEG was a gene with an absolute fold change of more than two (|log 2 FoldChange| ≥ 1) and a false discovery rate (FDR) ≤ 0.05.
Herbivory-responsive genes were defined as genes that show differential expression between control samples and herbivory samples (control vs. herbivory) in WT plants. Floweringassociated genes were defined as genes that show differential expression between the two developmental stages in WT plants.

Permutation Test
Permutation tests were used to determine whether the number of overlapping genes is significantly higher than the expected number under the null model (i.e. by chance). The analysis was carried out using R. In brief, we first simulated the expected number of overlapping genes under the assumption that the induced responses and flower development are independent processes. To this end, we randomly selected a certain number of genes (same as the number of DEGs) from each of the background gene lists and calculated how many genes were found in both processes. We repeated this process 10,000 times (denoted as N). We then compared how many times (n) the simulated results were higher than the observed value and estimated the value of p with the following formula: p = (n + 1) / (N + 1). To determine the overlap between herbivory-induced genes and flowering-associated genes, we used all genes expressed in leaves and flowers of WT plants, respectively, as the background gene list. To determine the overlap between LOXD/PS-regulated genes in leaves and flowers, we used all genes expressed in the respective tissues of WT, 35S::PS and spr8 plants as the background dataset. The expression data and R-scripts are deposited on figshare. 4

Annotation of Putative Gene Functions
Functional annotation of tomato genes was performed by eggNOG-mapper v2 5 (Huerta-Cepas et al., 2017, 2019) using the ITAG 3.2 protein sequences and predicted ORF protein sequences of the novel genes (resulting from transcript assembly) as input sequences, followed by the use of the R package AnnotationForge (v1.26.0; Carlson and Pagès, 2019). GO enrichment analysis was performed by R package clusterProfiler (v3.12.0; Yu et al., 2012) with all expressed genes as the background gene list. Enriched gene ontology (GO) terms with both value of p and value of q ≤ 0.05 were selected. The top 10 enriched GO terms (at level four) were visualised as bar plots.
To identify genes involved in phytohormone biosynthesis and signalling, photosynthesis and secondary metabolism, we conducted a functional annotation with the reference proteins (ITAG 3.2) and predicted ORFs of the novel genes using Mercator4 v2.0 6 (Schwacke et al., 2019). The genes involved in the different pathways were searched in MapMan (v3.6.0RC1; Thimm et al., 2004) based on the Mercator4 annotation and tomato metabolic pathway database. 7

Quantifying Phytohormones
Extraction was done similar as described by Bonaventure et al. (2011). In brief, approximately 10-15 mg fine-ground freezedried plant material was extracted twice with ethyl acetate (Fisher Chemical). The supernatants were collected, combined and evaporated to dryness. After reconstitution in 70% methanol (Fisher Chemical) in water, the samples were analysed on an AB Sciex QTRAP 6500+ equipped with an Ion Drive™ Turbo V source with TurboIonSpray Probe, which was operated in negative ionisation mode multi-reaction-monitoring modus. Settings were as follows: curtain gas, 45; collision gas, high; IonSpray voltage, −4,500; temperature, 500°C; ion source gas 1, 50; ion source gas 2, 70; and entrance potential, −10. For quantification, we added 10 ng D 6 -ABA (OlChemIm), 10 ng D 5 -JA (Sigma-Aldrich) per sample as internal standard during the first extraction step. JA-Ile was quantified based on the D 5 -JA standard and an empirically determined conversion factor (m JA-Ile = 0.274 × m D5-JA ).

Quantifying Specialized Metabolites
Tissue was lyophilized for at least 24 h and homogenised in GenoGrinder at 1,000 strokes per minute for 5 min. Metabolites were extracted from approximately 2-20 mg of lyophilized tissue in 1.5 ml 80% methanol containing 27 μg of umbelliferone (internal standard) for 2 h and constant shaking. Extracts were centrifuged for 40 min at 3,500 rpm, and 100 μl of supernatant was collected and analysed by LC-MS-QTOF. Extracts were analysed in AutoMSMS, negative mode using a Zorbax SB-C18 column (1.8 μm, 2.1 × 100 mm; Agilent). Quantification of α-tomatine was based on the mz1032.5385 (Rogachev and Aharoni, 2012) and the external standard.

Gene Expression Validation
To validate the expression changes observed from RNA-seq, we selected 10 genes that are putatively related to the JA pathway or defences and measured their expression using qPCR. For each sample, ~1 μg total RNA was used to synthesise firststrand cDNA using the GoScript™ Reverse Transcription Mix and Oligo(dT) kit (Promega, United States). The qPCRs reactions (in 10 μl) were prepared according to the manufacturer's instruction of UltraSYBR Mixture (High ROX; CWBIO, China). The qPCRs were performed on a StepOnePlus™ Real-Time PCR System (Applied Biosystems, United States) with following PCR parameters: predenaturation at 95°C for 10 min, 40 cycles of denaturation at 95°C for 15 s and annealing/elongation at 60°C for 1 min and melting curve was carried out in the end of each PCR with the default program to validate the specificity of primers. Primer sequences of the 10 selected genes are shown in Supplementary Table 1. Tomato UBI3, a common internal reference gene for analysing gene expression changes affected by the JA-signalling pathway or biotic stresses (Fowler et al., 2009), was used as a reference gene. The expression at stage 4 from WT was used for calibration. Relative expression values were calculated using the 2 −ΔΔCt method (Livak and Schmittgen, 2001). Three biological replicates (samples from three plants) and two technical replicates were used.

Herbivory-Induced Transcriptomic Responses in Tomato Leaves
To investigate the herbivory-induced transcriptomic responses in tomato, we sequenced and compared transcriptomes of leaves from both control and herbivory-induced samples. In comparison with control samples, herbivory increased the transcript levels of 806 genes and decreased the levels of 1,028 genes (FDR ≤ 0.05, fold change ≥ 2; Figure 1A; Supplementary Table 2). Gene ontology (GO) enrichment analysis suggests that M. sexta feeding increased the expression of genes enriched in the functional terms: response to wounding, jasmonic acid metabolic process, regulation of hormone levels, small molecule catabolic process and hormone-mediated signalling pathway ( Figure 1B). Genes that were suppressed by M. sexta feeding were enriched in the GO terms: photosynthesis, generation of precursor metabolites and energy, electron transport chain, cofactor metabolic process and pigment biosynthetic process ( Figure 1B).
Annotating gene functions using MapMan revealed that M. sexta feeding increased the transcript levels of many phytohormone biosyntheses and signalling genes, in particular of the JA and ET pathway (Supplementary Table 3). In the JA biosynthesis pathway, four LOXs, three AOSs, AOC, OPR3, OPCL1 and ACX1A were all significantly upregulated by herbivory (Supplementary Table 3). The same pattern was also found for several JA signalling genes: eight JASMONATE ZIM-domain (JAZ) genes, MYC2, and the targets of MYC2: MTB1, MTB2 and MTB3 (Supplementary Table 3). Consistently, herbivory also increased the levels of JA and JA-Ile in the leaves (Supplementary Figure 1) and upregulated JA-regulated defence genes, such as CDI, TD2, LAP2, ARG2 and PIs (Supplementary Table 3). In the ethylene (ET) biosynthesis pathway, M. sexta feeding increased the transcript levels of SAMS2, ACO1 and ACO2 (Supplementary Table 3). A similar pattern was also found for the ET signalling genes, such as ETR4, EIL1, EBF1/2 and several ERFs (Supplementary Table 3).
Transcription factors (TFs) are key elements for coordinated transcriptomic responses. In total, we found that 111 TFs from 29 families showed differential expression between control and M. sexta-attacked leaves ( Figure 1E; Supplementary Table 5). Among them, bHLH, AP2-EREBP, MYB, WRKY and TIFY were the most frequent TF families. Among the differentially expressed bHLH TFs, 13 out of 19 were upregulated by herbivory, most of which are involved in regulating JA signalling. The TIFY family contains mostly JAZs, which are known targets but also suppressors of JA signalling. All DEGs from the TIFY family were upregulated by M. sexta feeding. The six DEGs from the Sigma70-like family were all suppressed by herbivory.
In summary, these results show that while feeding by M. sexta elicited JA and ET biosynthesis and signalling, as well as the biosynthesis of terpenoids and phenylpropanoids, it suppressed the photosynthesis and primary metabolisms in tomato leaves. The process is likely mediated by expression changes of a few key TF families.

Transcriptomic Reprogramming During Flower Development in Tomato
To investigate transcriptomic changes during flower development, we collected and sequenced the transcriptomes of flower samples at two stages: stage 4 (flower bud shortly before opening) and stage 5 (fully opened flower). In total, 4,151 genes showed differential gene expression between the two stages, including 2,000 upregulated genes and 2,151 downregulated genes (Figure 2A; Supplementary Table 6). Among them, we found several previously reported genes involved in flower development and opening processes, such as the regulator of floral organ formation SEPALLATA3 (SEP3), the gene involved in flower development WUSCHEL (WUS) and the stamen differentiation regulator Tomato MADS-box gene 6 (TM6).
Gene ontology (GO) enrichment analysis showed that the upregulated genes are mainly enriched in the biological processes: regulation of pollen tube growth, cell surface receptor signalling pathway, pollen tube development, developmental cell growth and response to wounding ( Figure 2B). The downregulated genes are most enriched in the GO terms: plant-type cell wall organisation or biogenesis, carbohydrate metabolic process, cell wall macromolecule metabolic process, cell wall organisation and secondary metabolite biosynthetic process ( Figure 2B).
Functional annotation further revealed that during flower development, many DEGs are involved in phytohormone biosynthesis and signalling (Supplementary Table 7 Table 7). We further performed permutation tests to identify the key metabolic pathways that were altered during flower development. The results show that not only the primary metabolism, in particular the biosynthesis of sugars, fatty acids and amino acids, but also several secondary metabolite pathways were significantly altered during flower development ( Table 2). This includes the biosynthesis of phenylpropanoids, plant cell wall structures, polysaccharides and nitrogen-containing secondary compounds.
Detailed analysis of these metabolic pathways revealed that genes involved in carbohydrate metabolism were downregulated during flower development, especially the genes for the biosynthesis of sucrose, starch and nucleotide sugar (Supplementary Table 4). Genes involved in steroidal glycoalkaloid (SGA) biosynthesis, including SSR2, SMO3, HYD1, DWF5-2, GAMEs and the key regulator JRE4, showed higher expression in open flowers than in flower buds ( Figure 2C). In the phenylpropanoid biosynthesis, several genes involved in flavonoid biosynthesis, including PAL4, 4CL5, 4CL, Chalcone synthase (CHS2), Chalcone isomerase (CHI1 and CHIL), flavanone 3-hydroxylase (F3H) and flavonol synthase (FLS), showed higher expression at stage 4 than at stage 5. The same expression pattern was observed for genes involved in lignin biosynthesis ( Figure 2C). In contrast, flavonoid-3′-hydroxylase (F3′H_3), flavonoid-3′5′-hydroxylase (F3′5′H), dihydroflavonol 4-reductase (DFR) and anthocyanidin synthase (ANS_2) showed the opposite pattern ( Figure 2C). FIGURE 1 | expression abundance was log 10 -transformed and z-score normalised. (B) GO enrichment analysis for herbivory-induced (red column) and herbivoryrepressed genes (blue column). Only the top 10 most enriched biological process terms are shown in the plots. (C) Heat map showing the expression patterns of genes involved in photosynthesis. The internal three circles represent the control leaves, while the external three circles represent herbivore-infested leaves. Each treatment had three replicates. Gene expression levels were log 10 -transformed and z-score normalised. (D) Herbivory-responsive genes involved in the terpenoid, SGA and phenylpropanoid pathways. Only genes showing significantly differential expression between the two treatments are shown. Gene names in blue represent genes downregulated by feeding, while gene names in red represent upregulated genes. IPP, isopentenyl diphosphate; DMAPP, dimethylallyl diphosphate; GA-3P, glyceraldehyde-3-phosphate; GPP, geranyl diphosphate; FPP, farnesyl diphosphate; and GGPP, geranylgeranyl diphosphate; for the full names of the indicated genes please refer to Supplementary Table 13. (E) Number of TFs which were up-or downregulated by herbivory sorted by TF family. Only TF families with at least three differentially expressed TFs were shown.  We further identified 259 TFs (from 43 families) that were differentially expressed between the two flowering stages (Figure 2D; Supplementary Table 8). Among them, more than 56% (146/259) were upregulated, showing higher expression levels at stage 5 (open flowers) than stage 4. The most frequent TF family was MYB, which contained 34 DEGs (18 upregulated and 16 downregulated genes). This included the previously characterised MYB21, which is required for flower opening and the development of reproductive organs (Niwa et al., 2018). Among other frequent TF families, AP2-EREBP, NAC, bHLH and WRKY contained more upregulated than downregulated genes, whereas C2H2, MADS, TCP and TIFY contained more downregulated genes. Interestingly, herbivory upregulated 9 TIFY genes, whereas 8 TIFY genes were downregulated during flower development.
Together, these results reveal that flower development is mostly associated with a downregulation of JA biosynthesis and upregulation of ET signalling, as well as a suppression of carbohydrate and phenylpropanoid metabolism and increasing SGA biosynthesis.

The Shared Transcriptomic Responses in Flower Development and Leaf Herbivory
The similar metabolic responses found during flower development and leaf herbivory could be due to either a contribution of different genes from the same family or the same gene in these two seemingly distinct biological processes. Indeed, in the phenylpropanoid pathway, we observed different members of the same gene family that were involved in herbivory-induced responses and flower development. For example, in the anthocyanin biosynthesis ANS_1 (Solyc10g076660) was upregulated by herbivory, while ANS_2 (Solyc08g080040) was upregulated during flower development. Similar patterns were also found in the PAL, 4CL and C3H families (Figures 1D,  2C). This suggests that some of the metabolic responses that are elicited both by herbivory and during flower development are achieved via independent regulation of different duplicated gene copies.
A broader comparison revealed that among the 1,834 herbivory-induced genes, 411 were also found differentially    Table 9). Permutation tests reveal that the observed overlap is significantly higher than by chance (p < 0.001; 10,000 permutations). GO enrichment analysis showed that the overlapping DEGs are enriched in the GO terms: response to wounding and jasmonic acid metabolic process. Among the 411 overlapping genes, 276 showed increased transcript levels after M. sexta feeding, of which around 2/3 were downregulated during flower development ( Figure 3B). For example, genes in the JA biosynthesis (LOXs, AOS, AOC and OPCL1), JA signalling (JAZ3, JAZ4, JAZ6, JAZ9, JAZ11and MTB2) and JA-dependent defences (ARG2, TD2, CDI, LAP2, PPO-F and PIs) were upregulated by herbivory in leaves and downregulated during flower development ( Figure 3C). The same pattern was also found for various genes of the MVA pathway (ACAT3), flavonoid biosynthesis (4CL) and lignin biosynthesis (CCOAOMT; Figures 1D, 2C). In contrast, several genes of the ET pathway (ACO2, EBF1/2, ERF D2, ERF F4 and ERF F5) were upregulated both by M. sexta feeding and during flower development (Supplementary Tables 3 and 7). The other 135 overlapping genes were suppressed by M. sexta feeding in leaves, and around half of them were also downregulated during flower development ( Figure 3B). For example, several genes involved in carbohydrate metabolism and carotenoid biosynthesis were repressed by both herbivory and flower development (Figures 1D, 2C; Supplementary  Table 4). Interestingly, JAR4, a key regulator of the SGA biosynthesis, was downregulated by herbivore feeding but upregulated during flower development.
More than 7% (30 out of 411) of the overlapping genes are transcription factors (Supplementary Table 10). Among them, the five most frequent families are: AP2-EREBP, bHLH, TIFY, MYB and NAC ( Figure 3D). While most of the shared TFs (26 out of 30) were upregulated by M. sexta feeding in leaves, their expression patterns varied during flower development. More than half of the shared AP2-EREBP TFs (six out of seven) and bHLH TFs (three out of five) and all of the NAC TFs (three) were also upregulated during flower development. In contrast, all shared TIFY TFs (five) and MYB TFs (three) were suppressed in open flowers (stage 5).

Signalling Networks Regulate Both Induced Defences and Flower Development
The correlated transcriptomic changes during flower development and herbivory-induced responses in leaves suggest that the same genes and signalling cascades regulate both processes. To test this hypothesis, we further compared the transcriptomes of WT plants and mutant plants in which herbivory-induced responses were genetically manipulated. We reasoned that if a signalling network is involved in regulating both processes, the genetic manipulation of its key elements will result in expression changes of similar genes in both processes. Here, we focus on the comparisons between WT and two mutant/ transgenic plants: spr8 and 35S::PS. In spr8 plants, the expression of LOXD (a JA biosynthesis gene) is altered due to a mutation, which results in significantly lower levels of herbivory-induced JA and JA-Ile than in WT plants (Yan et al., 2013;Supplementary Figure 2). In 35S::PS plants, PS is overexpressed by a 35S promoter from the Cauliflower mosaic virus, resulting in high levels of endogenous JA and JA-regulated defence responses, even without herbivore attacks (Chen et al., 2006).
First, we compared both the leaf and flower transcriptomes of WT plants and spr8 mutants. Between leaves from WT and spr8 plants, only 105 genes were differentially expressed in control samples (without herbivory), while 1,253 genes were differentially expressed in herbivory-induced samples (Figure 4A), which included most of the well-studied herbivoryinduced defence responses, such as TD2, ARG2 and PIs. This is consistent with a previous study showing that the herbivoryinduced LOXD expression and JA signalling are impaired in the spr8 mutant (Yan et al., 2013). Functional enrichment analysis indicated that these 1,253 genes are enriched in the GO processes of: photosynthesis, generation of precursor metabolites and energy and response to wounding (Supplementary Figure 3A). Between flowers of WT and spr8 plants, the transcriptomes differed most at stage 4, where transcript levels of LOXD were significantly different ( Figure 4A). In total, we found 1,741 DEGs at this stage that are involved in pollen tube growth, response to wounding and pollination (Figure 4A, Supplementary Figure 3B). Among them, 166 were also found differentially expressed in herbivory-induced leaves (Figure 4B; Supplementary  Table 11). Permutation tests showed that the observed number of genes (166) found to be differentially expressed between WT and spr8 plants both in leaves and flowers is significantly higher than expected by chance (p < 0.001; 10,000 permutations). These 166 genes are primarily enriched in the biological process response to wounding and contain mostly wellcharacterised herbivory-induced JA-signalling and antiherbivore defence genes, such as LOXs, AOS, JAZs, TD2, ARG2 and LAP2. In addition, the gene expression changes between WT and spr8 plants in leaves and flowers were highly similar. For example, most genes downregulated in herbivoryinduced leaves of spr8 plants were also downregulated in the flowers of spr8 plants ( Figure 4C).
Second, we compared the leaf and flower transcriptomes of WT and 35S::PS plants. Comparing WT and 35S::PS plants revealed 1,368 and 1,116 DEGs in control and herbivory-induced leaves, respectively ( Figure 4D). This is consistent with the fact that the differences in transcript levels of PS between WT and 35S::PS plants were larger under control conditions than under herbivory (Figure 4D). GO analysis showed that the 1,368 genes are highly related to photosynthesis, generation of precursor metabolites and energy and pigment biosynthetic process (Supplementary Figure 3C). In flowers, the largest transcriptomic difference between WT and 35S::PS plants was found at developmental stage 5 (open flowers), where PS showed significantly different transcript levels ( Figure 4D). In total, 1,260 genes were differentially expressed, which are enriched in GO terms of: secondary metabolite biosynthetic process, regulation of hormone levels and response to nitrogen compound (Supplementary Figure 3D). Among them, 190 were differentially expressed in leaves (under control conditions; Figure 4E; Supplementary Table 12). Permutation tests showed that the observed number of genes (190) that were found both differentially expressed in leaves and flowers of WT and 35S::PS plants is significantly higher than expected by chance (p < 0.001; 10,000 permutations). These 190 genes are mostly involved in photosynthesis and herbivory-induced defence responses, such as PIs and LAP. While half of the genes activated by the over-expression of PS in leaves also showed increased expression in flowers, most genes that were suppressed by over-expression of PS in leaves also showed decreased transcript levels in flowers ( Figure 4F).
To validate the results from RNA-seq, we further measured the expression of 10 selected genes using qPCR, in both WT and spr8 flowers in the stage 4 and 5. The expression patterns of the selected genes observed from RNA-seq are overall highly consistent with qPCR (Supplementary Figure 4).
Together, these results reveal that herbivory-induced responses in leaves and flower development are partially regulated by the same signalling cascades.
Frontiers in Plant Science | www.frontiersin.org

Diffuse Coevolution of Defensive Metabolites in Tomato Leaves and Flowers
The shared signalling network predicts that defensive metabolites can show diffuse coevolution in flower development and herbivory-induced responses. To test this prediction, we analysed the concentration of several specialised metabolites and phytohormones (α-tomatine, kaempferol, rutin, JA, JA-Ile and abscisic acid) in six closely related wild tomato species. Among these wild tomatoes, only JA and α-tomatine, two key defence-related metabolites in tomato, showed correlated variations between flowers and herbivory-induced leaves (p < 0.05; after correction for phylogeny; developmental stages, JA only co-varied between the induced leaves and flower buds. Interestingly, while previous studies have shown that α-tomatine is regulated by JA (Cárdenas et al., 2016;Abdelkareem et al., 2017;Montero-Vargas et al., 2018;Nakayasu et al., 2018), no correlation was found between the levels of herbivory-induced JA and α-tomatine (p > 0.1), suggesting additional genetic mechanisms are also involved in herbivory-induced α-tomatine changes.

DISCUSSION
Both herbivory-induced responses in leaves and flower development are highly relevant processes for plant reproduction and fitness. However, they were often studied in isolation. Using a comparative transcriptomics approach, we found that the two seemly distinct biological processes require a similar transcriptional regulation: modulation of JA signalling, suppression of primary metabolisms and reprogramming of secondary metabolisms. While in some pathways different duplicated gene copies specifically responded to one or the other process, further genetic manipulations revealed that herbivory-induced responses in leaves and flower development are partially regulated by shared signalling networks. Together with correlated changes in key defence metabolites among six closely related tomato species, these results suggest that floral traits and herbivory-induced defences likely have evolved under diffuse selection due to pleiotropy.
In tomato leaves, M. sexta herbivory increased the transcript levels of genes involved in the JA and ET pathways (Supplementary Table 3), as well as the secondary metabolism, but decreased the transcript levels of genes involved in primary metabolism. This pattern is consistent with herbivory-induced responses in wild tobacco (Nicotiana attenuata; Diezel et al., 2009), Arabidopsis thaliana (De Vos et al., 2005;Ehlting et al., 2008), cotton (Gossypium hirsutum; Huang et al., 2015), maize and rice (Qi et al., 2018).
The elicitation of JA and ET biosynthesis and signalling can lead to reprogramming of the specialised metabolism, which is a widespread phenomenon in plants (Bennett and Wallsgrove, 1994;War et al., 2012;Fürstenberg-Hägg et al., 2013;Aljbory and Chen, 2018). Some of these induced responses result in increased plant defences. For example, in tomato, the increased expression of the anti-digestive proteins (such as PIs and TDs) can reduce the digestive capacity of herbivores (Chung and Felton, 2011), and the increased expression of terpene synthase genes (TPS5, TPS25 and TPS46) is involved in the release of volatiles that can attract their natural enemies ( Figure 1D; Paudel et al., 2019). Interestingly, although SGAs are known as effective defence compounds against biotic attackers because of their cytotoxic and anti-nutritional properties (Roddick, 1996;Milner et al., 2011), most of the SGA biosynthesis genes were not altered by herbivory in our experiments (except DWF7-2 and JRE4, which were downregulated by feeding; Figure 1D). In the tomato jre4 mutant, SGA accumulation is decreased and susceptibility to Spodoptera littoralis larvae is increased (Nakayasu et al., 2018). Because JRE4 is JA-inducible, it is reasonable to assume that JRE4 is part of the herbivoryinduced plant defence responses (Nakayasu et al., 2018). However, our data show that the expression of JRE4 was suppressed by M. sexta feeding, which could be due to the induction of the ET pathway that can suppress the JA-induced JRE4 expression in tomato leaves (Nakayasu et al., 2018). In addition to chemical defences, herbivore attack can also alter the cell wall metabolism. For example, in maize, the attack by Sesamia nonagrioides upregulates cell wall reorganisation and biogenesis, which is accompanied by increased lignin content (Rodríguez et al., 2012). Consistently, M. sexta-infested tomato leaves also show increased expression levels of genes involved in the phenylpropanoid pathway and lignin biosynthesis (PAL5, 4CL, 4CL2, C3H1, CCOAOMT, CCOAOMT2 and HCT; Figure 1D).
We found a majority of genes involved in carbohydrate metabolism and photosynthesis that showed a significantly decreased expression in herbivore-attacked leaves ( Figure 1C; Supplementary Table 4). Herbivory-induced suppression of photosynthesis and carbohydrate metabolism is common among different plant taxa after attack by herbivores from different feeding guilds (Halitschke et al., 2003;Zhu-salzman et al., 2004;Ralph et al., 2006;Bilgin et al., 2010), with the exception of mirid bugs feeding on wild tobacco (Halitschke et al., 2011). Similarly, herbivory and artificially increased endogenous jasmonate levels can lead to decreased carbohydrate concentrations (Babst et al., 2005;Machado et al., 2013Machado et al., , 2015Machado et al., , 2017Ferrieri et al., 2015). Therefore, our results fall in line with previous observations regarding changes in carbohydrate metabolism following herbivory.
During flower development, we also observed transcriptional changes of genes involved in the JA and ET pathways, the primary metabolism and the secondary metabolism. The observed changes in JA signalling during flower development are consistent with previous studies that showed the importance of JA signalling for flower development (Yuan and Zhang, 2015;Huang et al., 2017). Previous studies additionally showed that JA and JA-Ile levels peak at the bud stage and gradually decline to basal levels in fully opened flowers (Dobritzsch et al., 2015;Niwa et al., 2018). Consistently, we observed a higher expression of JA biosynthesis genes at stage 4 than at stage 5 ( Figure 3C;  Supplementary Table 7). In addition to JA, ET also plays a fundamental role in flower development in Arabidopsis, rice, China rose (Hibiscus rosa sinensis) and tomato (Iqbal et al., 2017). We observed that the transcript levels of ET biosynthesis and signalling genes were higher in fully opened flowers than in flower buds (Supplementary Table 7). This is consistent with the fact that ET regulates flower senescence and fruit ripening (Lanahan et al., 1994;Tieman et al., 2000;Alexander and Grierson, 2002;Shinozaki et al., 2015).
Flower development involves several metabolic shifts. First, flower development requires secondary cell wall biogenesis (Shalit-Kaneh et al., 2019). Our analysis showed that the expression of genes involved in lignin biosynthesis, including PAL4, 4CL, 4CL5, C3H2, CCOAOMT, CCOAOMT3, COMT and CAD, was suppressed in fully opened flowers in comparison with flower buds (Figure 2C). Similarly, in snapdragon flowers, downregulated genes during petal development were notably enriched for enzymes of secondary cell wall biogenesis (Muhlemann et al., 2012). Second, flower development involves petal expansion, which involves hydrolysis of stored carbohydrates (Bieleski, 1993;Bieleski et al., 2000;Vergauwen et al., 2000;Yap et al., 2008). This might explain the observed higher expression of genes involved in starch and sucrose degradation in fully opened flowers (stage 5) compared to flower buds (stage 4; Supplementary Table 4). Third, flower development also involves the increased biosynthesis of pigments and volatiles that are critical for pollinator attraction. We found that the expression of genes involved in late steps of carotenoids biosynthesis (e.g. CYCB, CYP97B3, NXD1, NSY2 and NSY3) was mostly suppressed in fully opened flowers (Figure 2C), while the expression of genes involved in anthocyanin biosynthesis, including F3′5′H, DFR and ANS_2, was increased. In addition, a key gene in flavonol biosynthesis, FLS, was suppressed in fully opened flowers ( Figure 2C). This suggests that the metabolism of flavonoids in tomato flowers shifts towards the production of anthocyanins during the opening process, which is similar to the pattern observed in Eustoma grandifloru (Kawabata et al., 2009). Interestingly, the expression of TPS7, a monoterpene synthase that mainly catalyses the formation of β-myrcene and limonene in tomato (Falara et al., 2011;Zhou and Pichersky, 2020), was reduced in fully opened flowers in comparison with flower buds. Surprisingly, most of the SGAs biosynthesis genes were upregulated during flower development, which might explain why the concentration of SGAs is higher in flowers than in other tissues (Friedman and Levin, 1995). The higher abundance of SGAs in flowers might also be related to the phytotoxicity of the intermediate metabolites in SGA biosynthesis (Itkin et al., 2011).
Although herbivory-induced responses and flower development are two distinct biological processes, they partially share the same regulatory mechanisms. We identified 411 genes that were affected by both herbivory and flower development ( Figure 3A). Among them, many genes are associated with phytohormone biosynthesis and signalling (e.g. JA and ET), carbohydrate metabolism, the MVA and MEP pathways, as well as the phenylpropanoid pathway. Further, the manipulation of key regulators (LOXD and PS) of the herbivory-induced response confirmed that the same signalling machinery is involved in both processes. The shared signalling machinery of both the herbivory-induced responses and flower development suggests diffuse coevolution of plant defence and floral signals (Jacobsen and Raguso, 2018;Rusman et al., 2019;Ramos and Schiestl, 2020). Interestingly, in the phenylpropanoid pathway, we observed that different duplicated gene copies specifically responded to one or the other process (Figures 1D, 2C), suggesting that plants can also escape from pleiotropic effects via gene duplications. However, a significantly higher number of genes than expected by chance is shared by both induced defence responses and flower development, indicating that some of the pleiotropic effects are likely adaptive. However, it is worth noticing that the observed pleiotropy is likely an underestimate. Due to the use of complete flower buds and flowers, we likely have missed the gene expression changes that are specific in different floral organs, which were marked due to pooling effects. Future studies that analyse the transcriptome of each floral organ will provide further insights into shared pleiotropy between induced responses and flower development.
Consistent with the transcriptomic data, we found the content of JA and α-tomatine was correlated between flower buds and herbivory-induced leaves among closely related wild tomato species (Figures 5A,B), suggesting that herbivory-induced defences and floral traits coevolved due to the same molecular regulatory networks. However, the sample size in this study is relatively small, which might limit the power to detect other co-evolving traits. Future studies that include large sample size and systematically compare metabolomes and transcriptomes between flowers and induced leaves will provide a better understanding on the extent to which diffuse coevolution shapes the evolution of defences and floral signals in nature.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm. nih.gov/, PRJNA600385.

AUTHOR CONTRIBUTIONS
SX conceived the study. SX, RZ and YS planned the research. LK, YW, MS, JF, TS and HP performed the experiment and analysed the data. SX, RZ, CD and YS contributed the resources. SX and LK wrote the manuscript. All authors contributed to the article and approved the submitted version.