Genomic and metabolic profiling of two tomato contrasting cultivars for tolerance to Tuta absoluta

The raw reads from Illumina RNA-Seq produced in our previous work (D’Esposito et al. 2021) were cleaned and filtered for quality. The high-quality reads were aligned to tomato genome reference SL3.0 by A.I.R platform (https://transcriptomics.sequentiabiotech.com↗). The BAM files generated for tolerant and susceptible genotypes were used for SNP calling using BCFtools (Li 2011; https://samtools.github.io/bcftools/bcftools.html↗). SNPs were filtered for quality using the following parameters: allele frequency (AF) higher than 0.75, minimum quality (QUAL) and minimum genotype quality (GQ) equal to 30, depth between 5 and 100 and minimum mapping quality (MQ) equal to 20. The program SnpEff (Cingolani et al. 2012; http://SnpEff.sourceforge.net/↗) was used to annotate private variants for each genotype based on their genomic locations and to predict variant effects.

The functional annotation of genes affected by variants with high (modification of gene products), moderate (missense variations due to changes amino acid codons) and modifier (non-coding variants or variants affecting non-coding genes) effect was performed using MapMan (Thimm et al. 2004). Assignment to MapMan classes was performed using the Mercator pipeline. Tomato pathway annotations for genes with variants were retrieved by the Plant Metabolic Network Database (Hawkins et al. 2021), TomatoCyc version 5.0.1 based on ITAG3.2 annotation.

The initial strain of T. absoluta was collected in 2017 in tomato greenhouses located in Battipaglia (Salerno, Italy). It was continuously reared at the Istituto per la Protezione Sostenibile delle Piante (IPSP) inside bug dorms^® isolators on tomato plants (cultivar ‘San Marzano nano’) at the following conditions: temperature of 24 ± 2 °C, relative humidity (RH) of 65 ± 5%, and photoperiod of 16L:8D.

The tomato genotypes BR221 (T = tolerant) and PS650 (S = susceptible) were grown in a glasshouse under the following conditions: temperature of 24 ± 2 °C, relative humidity (RH) of 65 ± 5%, and photoperiod of 16L:8D. Five-week-old plants, with 4–6 completely expanded leaves and a height of 18 cm, were placed into a single mesh cage (60 × 60x180 cm; Vermandel, Hulst, The Netherlands). Three to five days old mated females of T. absoluta were released into the cage keeping a 1:1 ratio between them and the plants. After two days (oviposition period), the females were removed from the cage by an insect aspirator. Plants were kept for 12–14 days in the same conditions reported above until the larvae hatched. Eleven replicates for each uninfested plants and plants infested by larvae were individually placed inside a 20-L glass jar for a 3-h volatile extraction from the headspace of those plants.

The volatile compounds were collected in Tenax traps and stored at − 20 °C for later volatile profiling by means of gas chromatography. After trapping on Tenax (30 mg) and carboxen (30 mg) packed tubes, the samples were analyzed by CIS4–TDU–GC/MS. Gerstel TDU (Gerstel, Mülheim, Germany) was heated at 300 °C for 7 min under a helium stripping flow of 30 mL min⁻¹. The TDU unit was directly assembled over the PTV injector (CIS4, Gerstel) with a liner-in-liner coupling, which eliminates the carryover effect and analyte loss. During this stage, the CIS4 was cooled to − 20 °C by computer-controlled liquid CO₂ pulsed flow. After cryo-trapping on a Tenax packing liner, the PTV was quickly ramped to 260 °C for desorption and the analyte was transferred to CIS4. An Agilent 7890 GC equipped with a 5975 MSD was used for the analysis, all from Agilent Tech (Palo Alto, CA, USA). Helium was used as the carrier gas, and the flow was kept constant at 1.2 mL/min. The chromatographic settings were as follows: injector in splitless mode set at 260 °C, J&W Innowax column (50 m, 0.20 mm i.d., 0.4 µm df); oven temperature program: initial temperature 40 °C for 1 min, then 10 °C min⁻¹ increase until 130 °C, then 5 °C min⁻¹ increase until 210 °C, then 20 °C min⁻¹ increase until 260 °C, hold time 3 min. The mass spectrometer was operating with an electron ionization of 70 eV, in scan mode in the m/z range 29–330, at three scans s⁻¹. The deconvoluted peak spectra obtained by Agilent Mass-Hunter software were matched against the NIST 11 spectral library for tentative identification. Kovats’ retention indices were calculated for further compound confirmation and compared with those reported in the literature for the chromatographic column used. Authentic standards were also injected to confirm compound identity.

The vast majority of variants was classified in the “modifier” category, including downstream gene variant, intron variant and upstream gene variant, whilst the less abundant corresponded to the “high” variation effect. A good number of missense variations showed a moderate effect (Fig. 3a, b). The most abundant variants with a high impact were the “splice donor,” (796 and 780, in T and S genotypes, respectively) followed by “splice acceptor variants” (320, T and 329, S) (Fig. 3b).

To explore the potential variations with effect at protein and transcript level, we focused on SNPs and InDels, with high, moderate and modifier effect in each tomato genotype. In T genotype 6294 genes were affected by SNPs and or InDels classified in one, two or all the three above mentioned classes, while 12,668 were found in the S. It is worth noting that, although the variants are different between the two tomato genotypes, they can affect common genes.

The prosystemin (Solyc05g051750) had variants with high effect in T genotype while S showed variants with modifier effect. Within JA biosynthesis, common and genotype-specific variants were identified (Supplementary Fig. S2). T-specific genes included two genes located in the HS of the chromosome 12 and HS-III on chromosome 2 (i.e. Solyc12g011040 and Solyc02g093140). In addition, T-specific variants also affected the lipoxygenase D on chromosome 3 (Solyc03g122340, TomLoxD). On the other hand, S genotype showed variants in LOXC.

T genotype also showed a higher number of fatty acid beta-oxidation multifunctional proteins with variants, including Solyc12g007170 in HS at the beginning of the chromosome 12. In this region, T-specific variants were also localized in subtilisin-like proteases (Solyc12g011140, Solyc12g011213). On the chromosome 2, T showed specific variants in Solyc02g094040, ortholog to Arabidopsis MPL1, MYZUS PERSICAE-INDUCED LIPASE 1 (At5g14180) and in genes involved in cell wall metabolism, including a glucan endo-1,3-beta-glucosidase (Solyc02g069700), four CASP-like proteins, three cobra proteins (Solyc02g089115, Solyc02g089120, Solyc02g089130). Hydroxyproline-rich glycoproteins (Solyc04g076410, Solyc04g074165) and rhamnogalacturonate lyases (Solyc04g076630, Solyc04g076640, Solyc04g076650, Solyc04g076660) (Deepak et al. 2007). Among genes with S-specific variants there was the histone-lysine N-methyltransferase, H3 lysine-9 specific SUVH1, Solyc10g077070. Common genes affected by variants in both genotypes include an another subtilisin-like protease (Solyc01g087850) affected by missense variation.

Genes involved in acylsugars biosynthetic pathway and acylsugars transport, were also affected by variants. Solyc03g121540, encoding for ACYLSUCROSE FRUCTOFURANOSIDASE 1 (ASFF1) and Solyc02g093180, encoding for an acyltransferase showed variants in T genotype (Supplementary Fig. S4). In addition, two genes with specific variants in T genotype (Solyc02g061740 and Solyc02g061750) and one with variants in both genotypes (Solyc04g025450), involved in the biosynthesis of acylsugar isoC5-CoA acyl chain precursor were found. In S genotype variants with impact on the ABC transporter Solyc03g005860 were identified.

Both genotypes showed genes with variants in camphor and linalool biosynthesis (Fig. 6). On the chromosome 2 of the T genotype, unique variants were identified in the terpene synthase Solyc02g079840 involved in the (−)-α-terpineol biosynthesis. Variants affecting five 7-deoxyloganetin glucosyltransferases (Solyc04g074340, Solyc04g074350, Solyc04g074370, Solyc04g074380, Solyc05g055840), located in HSs on chromosome 4 and 5, were also identified in T genotype. Differences were observed in the pentacyclic terpene biosynthesis leading to the production of lupeol and amyrin, due to two S-specific genes with variants (Solyc12g006520, Solyc12g006530) and one common to both genotypes (Solyc07g042630). A sesquiterpene synthase with T-specific variants, involved in the production of caryophyllene, was identified (Fig. 6).

Finally, two genes (Solyc02g069920, Solyc02g069925) involved in the production of 1,4-dihydroxy-2-naphthoate I, were polymorphic in T genotype. By contrast, the S genotype showed variants in a gene (Solyc01g108720) involved in methylsalycylate degradation.

The T and S genotypes showed a huge difference in the number of variants, suggesting that they have a divergent genomic background. The distribution of variants along the chromosomes, in both genotypes, showed that SNPs were more frequent in the distal part of chromosomes, which corresponds to regions with higher recombination frequency (Sim et al. 2012) and gene density (The Tomato Genome Consortium 2012). Nevertheless, a dissimilar pattern of variants distribution along the chromosomes, tracing a footprint of different history for the two genotypes, was highlighted. In addition, the identification of contrasting variant pattern between the two genotypes in specific regions of the chromosomes 2, 5, 8 and 12 suggested the presence of important loci involved in response to T. absoluta. The identification of private variants led us to identify important differences between the two genotypes. SNPs can have a great impact on the variation of genomes and the phenotypic traits (Shastry 2009) and the presence of private SNPs in the two genotypes could drive a different response to the insect (Mhoswa et al. 2020; Zhou and Jander 2022).

Several genes with private variants involved in processes associated to plant defense, including stimulus perception, hormone metabolism, cell wall, signaling and secondary metabolism, were identified in the two genotypes. Interestingly, the dissimilar HS regions on chromosome 12 and 2 showed many genes with variants belonging to different classes of receptor kinases that can affect the plant perception of the pathogen. Differences between the two genotypes also regarded the signaling. For example, on the chromosome 2, variants in genes belonging to kinase/phosphatase signaling were found in T genotype. In Solyc02g070260, a phosphatase having a high number of systemin-responsive partners (Ahmad et al. 2019) was affected by several variants. Systemin is a crucial small peptide with important functions in plant wound signaling and defense (Ryan and Pearce 2003). Solyc04g076460, a leucine-rich repeat receptor-like protein involved in plant innate immunity (Zhang et al. 2020b), showed interesting variants in T as well as the kinase SlSERK3B (Solyc01g104970). It is worth mentioning that the silencing of SlSERK3B resulted in enhanced susceptibility to root knot-nematode and non-pathogenic Pseudomonas syringae pv. tomato (Peng and Kaloshian2014). In addition, two genes coding for cyclic nucleotide-gated channels (Solyc02g086990 and Solyc02g088560), which have an important role in plant disease resistance and innate immunity (Ma et al. 2009), were affected by variants in T genotype.

Although both genotypes showed variants in the prosystemin gene (Solyc05g051750), the variants were different and had a different effect on the protein and transcript. Solyc05g051750 is constitutively expressed in plant and prompts, also in the absence of wounding, the synthesis of the wound-inducible defensive proteins, whereas its over-expression in antisense orientation blocks the systemic wound signaling (McGurl et al. 1994). Interestingly, transgenic tomato plants overexpressing a truncated prosystemin modulated the expression of defense-related genes against the lepidopteran pest Spodoptera littoralis (Molisso et al. 2022).

JA is a lipid-derived phytohormone that plays a crucial role in plants defense against herbivorous insects, and its biosynthesis regulation is very complex (Devoto and Turner 2003). The two genotypes differed for key genes involved in JA synthesis and affected by variants. For example, T-specific variants included the lipoxygenase D TomLoxDon chromosome 3 (Solyc03g122340), which is mainly involved in JA biosynthesis and, when overexpressed, leads to enhanced resistance to insect and necrotrophic pathogen (Yan et al. 2013). On the other hand, the S genotype showed variants in LOXC, which is chloroplast-targeted and generates volatile C6 flavor compounds from both linoleic and linolenic acids (Chen et al. 2021).

T-specific variants were localized in subtilisin-like proteases on chromosome 12, but interestingly a subtilisin-like protease located on chromosome 1 (Solyc01g087850), known to be involved in resistance to the tobacco hornworm Manduca sexta, is affected by variants in both genotypes (Meyer et al. 2016). Finally, variants were also found on Solyc02g094040, the ortholog of Arabidopsis MPL1, MYZUS PERSICAE-INDUCED LIPASE 1 (At5g14180), involved in defense against the green peach aphid (Louis et al. 2010).

The cell wall is a complex structure subject to dynamic remodeling that often determines the outcome of the interactions between plants and pathogens (Bellincampi et al. 2014). The presence of variants in cell wall related genes might suggest differences between the two genotypes in cellulose deposition or cell wall remodeling in response to T. absoluta. For example, on chromosome 2, the T genotype showed variants affecting a glucan endo-1,3-beta-glucosidase (Solyc02g069700), well known to be implicated in plant defense (Balasubramanian et al. 2012); four CASP-like proteins able to form a mechanical barrier to trap pathogens at the infection site (Lee et al. 2019); three cobra proteins (Solyc02g089115, Solyc02g089120, Solyc02g089130) involved in cell wall composition alteration (Roudier et al. 2005). Variants in the T genotype were also displayed by genes encoding for plant cell wall components involved in defense against pathogens, such as hydroxyproline-rich glycoproteins (Solyc04g076410 and Solyc04g074165) and rhamnogalacturonatelyases (Solyc04g076630, Solyc04g076640, Solyc04g076650, Solyc04g076660) (Deepak et al. 2007).

Several TF genes with a role in biotic stress response, showing variants in the T genotype, were located on chromosome 2. Solyc02g037530, an auxin response factor 8B that promotes JA production (Nagpal et al. 2005; Liu et al. 2014); two zinc finger proteins LSD1, involved in the plant immune system (Solyc02g078270 and Solyc02g069720) (Alves et al. 2014), three WRKY including Solyc02g094270 (Zhang et al. 2020a), Solyc02g088340 SlWRKY3, which act as a positive regulator of induced resistance in response to nematode invasion and infection (Chinnapandi et al. 2019), Solyc02g093050 that is ortholog to the Arabidopsis ATWRKY15, induced by herbivores (van Aken et al. 2016; Rushton); and the MYB-SlMIXTA-like gene (Solyc02g088190) involved in trichome formation (Galdon-Armero et al. 2020; Ying et al. 2020) but also in the regulation and the production and storage of specialized toxic metabolites (Ying et al. 2020). In addition, HD-Zip TF (Solyc08g066500), a negative regulatory of the lignin biosynthesis (Liu et al. 2021), and ULTRAPETALA (Solyc12g010755), a regulator of the biotic and abiotic stress response (Tyler et al. 2019), also showed variants in T genotype.

By contrast, S-specific variants affected a Woolly gene (Wo), encoding a HD-Zip protein essential for trichome formation (Yang et al. 2011) and Solyc12g005830, SlHDZIV8/similar to the HDG2, HOMEODOMAIN GLABROUS 2 gene, participating in trichome development in Arabidopsis (Marks et al. 2009). S-specific variants were also identified in genes involved in the trichome formation: Solyc10g077070, a gene characterized by Fonseca et al. (2022) as HAIRPLUS (HAP), which controls glandular trichome density in tomato plants, and the CUTIN DEFICIENT2 (Solyc01g091630), well known gene that regulates cuticle deposition and the formation of glandular type-VI trichomes. CUTIN DEFICIENT2 variant alleles promoted reduced trichome density and lower volatile terpene production (Nadakuduti et al. 2012).

Many genes with variants were identified in secondary metabolism and were involved in the biosynthesis of VOCs, including phenylpropanoids, flavonoids and terpenoids, which play a key role in direct or indirect defense (War et al. 2012; Cheynier et al. 2013; Kessler 2017; Ameye et al. 2018; Erb and Kliebenstein 2020). In the cultivated tomato, acyl sugars, flavonoids, and terpenes are major secondary compounds produced in type I, IV and VI gland bearing trichomes in response to herbivore attack (Schilmiller et al. 2009; McDowell et al. 2011; Kang et al. 2014; Bergau et al. 2015; Balcke et al. 2017).

Caffeoyl-CoA O-methyltransferase (CCoAOMT) is essential in lignin biosynthesis, an important barrier that protects against pests and pathogens (Liu et al. 2018). Notably, a CCoAOMT (Solyc02g093230) affected by variants in T genotypes was located in a QTL region implicated in resistance against pathogens and herbivores (Vosman et al. 2018, 2019).

In both genotypes, the identification of genes involved in acylsugar biosynthesis and export, affected by variants, could suggest a diversity in the acylsucrose types and in acylglucose abundance. Genes involved in acylsugars biosynthetic pathway are well known to be activated in the trichomes of solanaceous plants to combat herbivores and pathogens (Fobes et al. 1985; Kroumova et al. 2016; Moghe et al. 2017). Interestingly, T variants were found in Solyc03g121540, encoding for ACYLSUCROSE FRUCTOFURANOSIDASE 1 (ASFF1), a trichome gland cell–expressed invertase that cleaves the glycosidic bond of P-type acylsucroses to generate acylglucoses (Leong et al. 2019). It is worth to note that ASFF1 gene is located in the QTL AG3.2 region at the bottom of chromosome 3 (Leckie et al. 2013; Leong et al. 2019). In addition, the presence of variants in genes involved in the biosynthesis of acylsugar isoC5-CoA acyl chain precursor could lead to a different use of isoC5-CoA as a donor, leading to an accumulation of diverse sets of acylsucrose structures (Fan et al. 2019). The S genotype could be affected in acylsugar export due to variants in the ABC transporter Solyc03g005860, previously associated with acylsugar exudation (Mandal et al. 2020). Acylsugar transport could be critical in determining how much acylsugar is produced and secreted, with significant consequences on plant defense. Instead, the T variants discovered in Solyc02g093180 could suggest an impact on the trichome production of acylsugars. This gene, encoding for an acyltransferase involved in acyl sugar biosynthesis (Fan et al. 2016), is located in a QTL region at the bottom of the chromosome 2 (Vosman et al. 2019), related to the presence of trichome type IV (Wf-1).

Among the plant secondary metabolites involved in biotic stress tolerance, terpenoids are the most diverse (Gershenzon and Dudareva 2007; Coppola et al. 2018; Boncan et al. 2020). The presence of variants in genes involved in the biosynthesis of different classes of terpenes in both genotypes could drive a metabolic diversity that can affect the response to T.absoluta. For example, T unique variants in the terpene synthase Solyc02g079840, located on the chromosome 2, were found. Such gene is involved in the (−)-α-terpineol biosynthesis, a compound with insecticidal properties (Khaleel et al. 2018) as well in 7-deoxyloganetin glucosyltransferases, enzymes involved in the synthesis of the monoterpene-derived compounds, called iridoids glycosides, that have a potential role in defense against herbivores (Puttick and Bowers 1988; Biere et al. 2004; War et al. 2018). In addition, both genotypes had variants in genes involved in the production of lupeol and amyrin, which are pentacyclic terpene with a major role as precursors for specialized triterpenoid metabolites, involved in plant defense and development (Cárdenas et al. 2019). T showed also specific variants in a gene producing caryophyllene that functions as a signal in the plant defense against herbivores (Köllner et al. 2008). Finally, a BHLH transcription factor (Solyc09g083360) involved in terpene biosynthesis and resistance against cotton bollworm and B. cinerea (Cao et al. 2022) and Solyc12g09990, scarecrow-like 3 (SlSCL3) involved in the production of monoterpenes and sesquiterpenes (Yang et al. 2021) showed variants in S genotype.

VOCs play an important role in the interactions of tomato genotypes with insect pests and also with their antagonists (Gontijo et al. 2019) and some of the compounds identified could be considered markers to separate genotypes resistant to T. absoluta infestation.

All compounds identified in the present work are in accordance with literature data concerning the VOC released by tomato plants upon herbivore feeding (i.e. T. absoluta larvae) (Proffit et al. 2011; Anastasaki et al. 2015, 2018; Catola et al. 2018; Milonas et al. 2019; Subramani et al. 2021; Ayelo et al. 2022; Deletre et al. 2022; Miano et al. 2022). In particular, we recorded a blend of volatiles emitted by tomato plants infested by Tuta larvae quantitatively different with respect to undamaged control plants as already reported (Anastasaki et al. 2015, 2018; Milonas et al. 2019; Ayelo et al. 2022). We highlighted an increased level of δ3-carene, ɣ-terpinene, ocimene and methyl-salicylate in both S and T genotypes plants and a higher level of ɑ-pinene, β-pinene, linalool, cis-3-hexenol and nonanal in T genotype plants infested by T. absoluta larvae compared to uninfested ones (Anastasaki et al. 2015, 2018; Milonas et al. 2019; Ayelo et al. 2022). Almost all compounds identified in this study are recognized by T. absoluta antennae (Anastasaki et al. 2018; Miano et al. 2022) and some of them have been recently evaluated for their possible use as repellents towards adult pests (Essoung et al. 2020; Miano et al. 2022). For example, α-pinene, β-pinene are oviposition deterrent for several pests of stored food products (Regnault-Roger and Hamraoui 1995; Ferrarini et al. 2008; Chaubey 2012) and the potato tuber moth, Phthorimaea operculella (Sharaby et al. 2009). In line with these evidences, we found a significant reduction of α-pinene (− 63.5%), β-pinene (− 92.6%) released by S genotype plants following larval feeding. Conversely, we recorded a significant increase of these two compounds (+ 155.9% and 158.0%, respectively) in the blend released by T genotype plants upon larval feeding. Similarly, in T we found a significantly higher release of compounds with repellence proprieties toward T. absoluta adults (Essoung et al. 2020; Miano et al. 2022) such as camphene, eucalyptol and camphor these (+ 275.0%, + 1287.1% and + 310.1%, respectively). Changes in the level of these marker VOCs could be associated with the tolerance/susceptibility to the tomato leaf miner.

DD was centrally involved in gene functional annotation, data interpretation and in manuscript writing. AG and CGA were involved in variant identification and annotation. PC and EG performed infestation experiments, VOC collection, analysis and interpretation. GC and MM performed VOC analysis and interpretation. MRE conceived the study and was mainly involved in data interpretation and in manuscript writing. All authors read and approved the final manuscript.

Below is the link to the electronic supplementary material. Supplementary file1 (PDF 215 KB) Supplementary file2 (XLSX 4465 KB)